API Reference

solve_ECG(operators, n=50; kwargs...) -> SolverResults

Build an ECG basis of n Rank0Gaussian functions using stochastic greedy search and return the ground-state energy.

Candidate Gaussians are generated from a quasi-random sequence (Halton by default). Each candidate is accepted if it is linearly independent from the existing basis (normalised overlap < threshold) and does not make the overlap matrix ill-conditioned. The ground-state energy after each accepted function is stored in SolverResults.energies.

Arguments

• operators : Vector{<:Operator} — kinetic + Coulomb operators (see KineticOperator, CoulombOperator).
• n : target number of basis functions (default 50).

Keyword arguments

Example

using FewBodyECG
masses = [1.0e15, 1.0]   # hydrogen atom (fixed nucleus)
Λmat = Λ(masses)
_, U = _jacobi_transform(masses)
w = U' * [1.0, -1.0]
ops = Operator[KineticOperator(Λmat); CoulombOperator(-1.0, w)]
sr = solve_ECG(ops, 30; scale=1.0, verbose=false)
println(sr.ground_state)   # ≈ -0.5 Ha

solve_ECG_variational(operators, n; kwargs...) -> SolverResults

Minimise a variational loss over the parameters of a Rank0Gaussian ECG basis using any OptimKit.jl optimisation algorithm.

Two loss types are available via the loss_type keyword:

• :energy (default) — minimise the lowest generalised eigenvalue λ_min of Hc = λSc. By the variational principle λ_min ≥ E₀ for any basis, so its global minimum is the exact ground-state energy. This is the recommended choice and is equivalent to the standard Rayleigh-Ritz variational principle applied to the Gaussian parameters. Gradients are computed via ForwardDiff and the Hellmann-Feynman theorem (differentiating through H and S only, not the eigensolver).

• :trace — minimise Tr(S⁻¹H), the sum of all generalised eigenvalues. For a basis already close to the physical optimum this can accelerate convergence, but for a randomly initialised basis whose upper eigenvalues are large and positive the optimizer may find a degenerate near-zero minimum (all Gaussians spread to infinity). Prefer :trace only when warm-starting from a stochastic result.

The A matrices of the basis functions are parameterised through their Cholesky factors (log-diagonal) which keeps them positive-definite for any parameter vector. Shift vectors s are also included in the optimised parameter vector (unconstrained), so the full variational freedom of each Gaussian is exploited.

After optimisation the full generalised eigenvalue problem Hc = λSc is solved once to produce the ground-state energy and eigenvectors, so all SolverResults-based utilities (ψ₀, correlation_function, etc.) work unchanged.

Arguments

• operators : Vector{<:Operator} — kinetic + Coulomb operators
• n : number of basis functions (default 50)

Keyword arguments

Example

using FewBodyECG, OptimKit
masses = [1.0e15, 1.0, 1.0]
Λmat   = Λ(masses)
J, U   = _jacobi_transform(masses)
w_list = [[1,-1,0],[1,0,-1],[0,1,-1]]
w_raw  = [U'*w for w in w_list]
ops    = Operator[KineticOperator(Λmat);
                  [CoulombOperator(c,w) for (c,w) in zip([-1.,-1.,1.], w_raw)]...]

# Default L-BFGS (convenience kwargs control it)
sr = solve_ECG_variational(ops, 30; scale=1.0, max_iterations=300, verbose=false)
println(sr.ground_state)

# Pass any OptimKit algorithm directly
sr2 = solve_ECG_variational(ops, 30; scale=1.0,
                            optimizer=LBFGS(; maxiter=500, gradtol=1e-8, verbosity=2))
println(sr2.ground_state)

# Warm-start refinement with conjugate gradient
sr0   = solve_ECG(ops, 30; scale=1.0, verbose=false)
basis0 = BasisSet(Rank0Gaussian[sr0.basis_functions...])
sr_cg = solve_ECG_variational(ops, 30; initial_basis=basis0, loss_type=:trace,
                              optimizer=ConjugateGradient(; maxiter=200, verbosity=0))
println(sr_cg.ground_state)

solve_ECG_sequential(operators, n=50; kwargs...) -> SolverResults

Build an ECG basis of n Rank0Gaussian functions using sequential variational optimisation (the Stochastic Variational Method, SVM).

At each step k = 1, …, n:

1. n_candidates quasi-random Rank0Gaussian functions are generated.
2. The candidate giving the lowest ground-state energy (evaluated without further optimisation) is appended to the current basis.
3. All k × n_per parameters of the combined basis are jointly optimised by an L-BFGS minimisation of the chosen loss_type.

Compared to solve_ECG_variational, which optimises all n functions simultaneously from a cold start, the sequential approach:

• Avoids the high-dimensional landscape of a full cold start.
• Produces a monotone non-increasing convergence curve (energies[k] after each step).
• Closely mirrors the SVM algorithm of the ECG literature (see e.g. Suzuki & Varga 1998).

Arguments

• operators : Vector{<:Operator} — kinetic + Coulomb operators
• n : number of basis functions (default 50)

Keyword arguments

Example

using FewBodyECG
masses = [1.0e15, 1.0, 1.0]   # H⁻ (fixed nucleus + two electrons)
Λmat   = Λ(masses)
_, U   = _jacobi_transform(masses)
w_list = [[1,-1,0],[1,0,-1],[0,1,-1]]
w_raw  = [U'*Float64.(w) for w in w_list]
ops    = Operator[KineticOperator(Λmat);
                  [CoulombOperator(c,w) for (c,w) in zip([-1.,-1.,1.], w_raw)]...]

sr = solve_ECG_sequential(ops, 30; scale=1.0, verbose=false)
println(sr.ground_state)   # converges toward -0.52775... Ha

Operators

Accumulates KineticOperator and CoulombOperator terms to build a Hamiltonian. Handles Jacobi coordinate transforms internally so that callers can work with physical particle indices rather than Jacobi-frame weight vectors.

Constructors

Operators()                    # system-unaware; add pre-built operators with `+=`
Operators(masses)              # system-aware; enables string/index shorthand
Operators(masses, charges)     # fully automatic; enables `ops += "Coulomb"` shorthand

System-aware interface

Particle indices follow the original ordering of masses. All Jacobi transforms are computed internally.

ops = Operators([m₁, m₂, m₃])
ops += "Kinetic"
ops += ("Coulomb", 1, 2, +1.0)   # pair (1,2) with coupling coefficient +1.0
ops += ("Coulomb", 1, 3, -1.0)

When charges are also supplied, the fully-automatic shorthand ops += "Coulomb" adds all $N(N-1)/2$ pairwise terms with coefficients $q_i q_j$:

# Helium atom: nucleus (Z=2), two electrons
ops = Operators([1e15, 1.0, 1.0], [+2, -1, -1])
ops += "Kinetic"
ops += "Coulomb"   # adds (1,2)→-2, (1,3)→-2, (2,3)→+1 automatically

System-unaware interface (pre-built operators)

ops = Operators()
ops += KineticOperator(Λmat)
ops += CoulombOperator(-1.0, w)

Both interfaces can be mixed freely. Pass ops directly to solve_ECG, solve_ECG_variational, solve_ECG_sequential, or build_hamiltonian_matrix. Use coulomb_weights to retrieve the Jacobi-frame weight vectors for manual basis construction.

coulomb_weights(ops::Operators) -> Vector{Vector{Float64}}

Return the Jacobi-frame weight vectors for every CoulombOperator in ops, in the order they were added. Useful for manual basis construction:

w_jac = coulomb_weights(ops)
A = _generate_A_matrix(bij, w_jac)

KineticOperator(K)
KineticOperator(masses)

Kinetic-energy operator in Jacobi coordinates.

When constructed from a mass vector the Jacobi-transformed kinetic-energy matrix $\Lambda = J M^{-1} J^T / 2$ is computed automatically via Λ.

Fields

• K : symmetric $n_{\text{dim}} \times n_{\text{dim}}$ kinetic-energy matrix ($\Lambda$).

CoulombOperator(coefficient, w)

Two-body Coulomb ($1/r_{ij}$) interaction operator.

The inter-particle distance is $|w^T \mathbf{r}|$ where w is a weight vector in Jacobi coordinates selecting the pair $(i,j)$. Construct w by transforming the charge-difference vector with the inverse Jacobi matrix: w = U' * charge_vector.

Fields

• coefficient : coupling constant (e.g. $q_i q_j$; negative for attraction).
• w : weight vector in Jacobi coordinates.

SolverResults

Container returned by all ECG solvers (solve_ECG, solve_ECG_variational, solve_ECG_sequential).

Fields

• basis_functions : accepted basis functions.
• n_basis : number of accepted basis functions.
• operators : the operator list used during the solve.
• method : solver symbol (:quasirandom, :random, :variational, :sequential).
• sampler : quasi-/pseudo-random sampler used for basis generation.
• length_scale : Gaussian width scale passed at construction.
• ground_state : energy of the target eigenstate (ground state by default).
• state : which eigenstate was targeted (1 = ground state, 2 = first excited, …).
• energies : energy after each accepted basis function (stochastic) or after each step (sequential).
• eigenvectors : list of eigenvector matrices; eigenvectors[end][:, state] is the target-state coefficient vector.
• fg_history : monotone-decreasing objective value after each gradient evaluation (variational solvers).

convergence(sr::SolverResults) -> (indices, energies)

Return the greedy build-up convergence curve from a stochastic solve_ECG run.

Returns (1:n_basis, sr.energies): the energy after each basis function was added. For variational results energies contains only one entry [ground_state]; use convergence_history instead.

convergence_history(sr::SolverResults) -> (indices, energies)

Return the per-objective-call convergence history.

For solve_ECG_variational results this is the cumulative-minimum energy at every primal fg evaluation, giving a monotone non-increasing curve suitable for plotting optimisation progress. The x-axis is the fg-call index.

For solve_ECG results this mirrors convergence.

correlation_function(sr; rmin=0.01, rmax=10.0, npoints=400,
                     coord_index=1, normalize=true)

Compute the one-body radial density $\rho(r) = r^2 |\psi_0(r)|^2$ along a single Jacobi coordinate.

Arguments

• sr : SolverResults from either solver.
• rmin, rmax : radial grid range (a.u.).
• npoints : number of grid points.
• coord_index : which Jacobi coordinate to scan (default 1).
• normalize : if true, normalise so that $\int \rho(r)\,dr = 1$.

Returns (r_grid, ρ) as plain Vector{Float64}.

ψ₀(r, c, basis_fns)
ψ₀(r, sr; state=1)

Evaluate the ground-state wavefunction at Jacobi-coordinate point r.

The wavefunction is the linear combination $\psi_0(\mathbf{r}) = \sum_i c_i \exp(-\mathbf{r}^T A_i \mathbf{r} + \mathbf{s}_i^T \mathbf{r})$.

When called with a SolverResults, the eigenvector for the requested state (default 1, i.e. the ground state) is used automatically.

Λ(masses) -> Symmetric matrix

Compute the kinetic-energy matrix in Jacobi coordinates for a system with the given particle masses (in atomic units).

Returns the symmetric matrix $\Lambda = J M^{-1} J^T / 2$, where $J$ is the Jacobi transformation matrix and $M = \operatorname{diag}(m_i)$. Pass the result directly to KineticOperator.

_jacobi_transform(masses) -> (J, U)

Compute the Jacobi coordinate transformation matrix J and its pseudo-inverse U for a system with the given particle masses.

Returns (J, U) where:

• J is the $(N-1) \times N$ matrix mapping particle coordinates to Jacobi relative coordinates (centre-of-mass motion is factored out).
• U = \operatorname{pinv}(J) is the $N \times (N-1)$ back-transformation.

The weight vectors for CoulombOperator are constructed as U' * charge_vector.

GaussianBase

Abstract supertype for all explicitly correlated Gaussian basis functions. Concrete subtypes differ by the rank of the polynomial prefactor: Rank0Gaussian (plain Gaussian), Rank1Gaussian (linear prefactor), Rank2Gaussian (quadratic prefactor).

Rank0Gaussian(A, s)

Basis function $g(\mathbf{r}) = \exp(-\mathbf{r}^T A\,\mathbf{r} + \mathbf{s}^T\mathbf{r})$.

Fields

• A : symmetric positive-definite $n_{\text{dim}} \times n_{\text{dim}}$ matrix controlling the Gaussian width and correlations.
• s : shift vector $\mathbf{s} \in \mathbb{R}^{n_{\text{dim}}}$; controls the location of the Gaussian maximum.

Rank1Gaussian(A, a, s)

Rank-1 (p-wave-like) ECG basis function with linear prefactor.

a can be either:

• a vector of length size(A,1) (single polarization component), or
• a matrix of size size(A,1) × ncomp (multi-component polarization).

Rank2Gaussian(A, a, b, s)

Rank-2 (d-wave-like) ECG basis function with quadratic prefactor.

a and b can each be either vectors or matrices. Their first dimension must match size(A,1). For matrix polarizations, a and b must have the same number of columns (ncomp), enabling multi-component pure d-wave channels.

BasisSet(functions)

A collection of GaussianBase functions that form the variational basis.

Fields

• functions : Vector{G} of basis functions, all of the same concrete GaussianBase subtype G.