API Reference

Operators

Accumulates kinetic and potential terms to build a Hamiltonian, including NumericalPotential terms whose radial quadrature is handled internally. 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

Adding terms

push!(ops, term) appends term to ops in place and returns ops. ops + term returns a new Operators and leaves ops unchanged, so ops += term rebinds ops to the extended copy. Both accept every term form shown below.

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)
ops += (r -> -exp(-r^2), numerical, 1, 2)

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 or build_hamiltonian_matrix. Use coulomb_weights to retrieve the Jacobi-frame weight vectors for manual basis construction.

Operator

Alias for FewBodyHamiltonians.Operator, exported so raw operator vectors can be typed as Operator[...] alongside the Operators builder.

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.

GaussianOperator(coefficient, γ, w)

Two-body Gaussian potential $V(r_{ij}) = \text{coefficient} \cdot e^{-\gamma r_{ij}^2}$ operator, where $r_{ij} = |w^T \mathbf{r}|$ is the inter-particle distance in Jacobi coordinates selected by the weight vector w.

The matrix element reduces to an overlap with a shifted exponent matrix: $S' = A + B + \gamma\, w w^T$, making evaluation exact and free of special functions.

Fields

• coefficient : coupling constant (negative for attractive well).
• γ : inverse-square range parameter ($\gamma > 0$).
• w : weight vector in Jacobi coordinates selecting the pair.

OscillatorOperator(coefficient, w)

Harmonic (oscillator) two-body potential $V = \text{coefficient}\cdot|w^T\mathbf{r}|^2$, where $w^T\mathbf{r}$ is the inter-particle coordinate selected by the Jacobi weight vector w.

Fields

• coefficient : coupling constant.
• w : weight vector in Jacobi coordinates selecting the pair.

ManyBodyGaussianOperator(coefficient, W)

Many-body Gaussian interaction $V = \text{coefficient}\cdot\exp(-\mathbf{r}^T W\,\mathbf{r})$ with a symmetric positive-definite exponent matrix W acting on all Jacobi coordinates at once (e.g. a repulsive regulator).

Fields

• coefficient : coupling constant.
• W : symmetric positive-definite exponent matrix.

NumericalPotential(f, w; rtol = 1e-8, atol = 0.0, maxevals = 10_000)

Numerical radial pair potential $f(|w^T r|)$ in Jacobi coordinates.

The callable f receives a nonnegative scalar distance. Matrix elements are evaluated by analytically reducing the ECG product to the radial coordinate selected by w, followed by adaptive numerical quadrature.

Fields

• f : user-supplied callable evaluated as f(r).
• w : Jacobi-coordinate weight vector selecting the pair coordinate.
• rtol : relative quadrature tolerance.
• atol : absolute quadrature tolerance.
• maxevals : maximum number of quadrature function evaluations.

numerical

Marker used by the system-aware numerical-potential shorthand: ops += (f, numerical, i, j).

GaussianTensorOperator(coefficient, γ, w, i, j; traceless = true)

Gaussian-form tensor interaction coupling the coordinate wᵀr (range γ) to the spins on sites i and j. With traceless = true the rank-2 spatial tensor rₐr_b − ⅓r²δₐ_b is used.

GaussianSpinOrbitOperator(coefficient, γ, w, i, j)

Gaussian-form spin-orbit interaction coupling the orbital motion in the coordinate wᵀr (range γ) to the total spin Sᵢ + Sⱼ. Produces complex Hermitian matrix elements for shifted Gaussians.

SpinProjection

Spin-½ projection eigenstates: up (+½) and down (−½).

SpinState(projections)

Direct-product spin-½ state: one SpinProjection per particle site.

SpinGaussian(orbital, spin)

An explicitly correlated Gaussian with an attached direct-product spin state. orbital is a Rank0Gaussian; spin is a SpinState. Only introduced to support the tensor and spin-orbit interactions; central operators factor through the spin overlap.

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} + \operatorname{tr}(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 supervector of size $n_{\text{dim}} \times 3$; row i is the Cartesian shift of Jacobi coordinate i. A length-N vector is accepted and mapped to the z component.

Rank1Gaussian(A, a, s)

Rank-1 (p-wave-like) ECG basis function $g(\mathbf{r}) = (a \cdot \mathbf{r}) \exp(-\mathbf{r}^T A\,\mathbf{r} + \operatorname{tr}(s^T \mathbf{r}))$ with $a \cdot \mathbf{r} = \operatorname{tr}(a^T \mathbf{r})$.

The polarization a and the shift s are N × 3 supervectors (row = Jacobi coordinate, column = Cartesian component), where N = size(A, 1). Each accepts a length-N vector, which is placed in the z component. A, a, and s are converted to a common element type.

Rank2Gaussian(A, a, b, s)

Rank-2 (d-wave-like) ECG basis function $g(\mathbf{r}) = (a \cdot \mathbf{r})(b \cdot \mathbf{r}) \exp(-\mathbf{r}^T A\,\mathbf{r} + \operatorname{tr}(s^T \mathbf{r}))$.

The polarizations a, b and the shift s are N × 3 supervectors, as for Rank1Gaussian; a length-N vector is placed in the z component. Orthogonal Cartesian polarizations such as a = [1 0 0], b = [0 1 0] give pure d-wave (xy) character.

BasisSet(functions)

A collection of GaussianBase functions that form the variational basis.

Supports the standard container interface: length, iteration, and indexing (including begin/end), with eltype reporting the concrete basis-function type G.

Fields

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

solve(ops, alg::SolverMethod = SVM();
      state = 1, tol = 1e-4, window = 20, init = nothing, verbose = false)

Solve the few-body eigenproblem defined by ops (an Operators builder or a raw Vector{<:Operator}) with algorithm alg — one of SVM, Refine, GVM, DynamicGVM, or a Pipeline composed with →.

Problem-level keywords: state targets the state-th eigenvalue, tol (absolute, Hartree) and window define the stochastic saturation criterion, init warm-starts from a previous Solution.

Returns a Solution carrying energies, the basis, S-orthonormal coefficients, and an honest ConvergenceReport.

SolverMethod

Abstract supertype of all solver algorithms. A method is a small struct of algorithm-level options; problem-level options (state, tol, window, init, verbose) live on solve. Adding a new method = defining a new subtype plus solve/step! methods — pure multiple dispatch.

StochasticMethod <: SolverMethod

Abstract supertype for stochastic basis-selection methods such as SVM and Refine.

GradientMethod <: SolverMethod

Abstract supertype for gradient-optimization methods such as GVM and DynamicGVM.

SVM(basis; candidates = 25, scale = :auto, sampler = HaltonSample(), indep_tol = 1e-4)

Suzuki–Varga stochastic selection (Sect. 4.2.5). At each of basis steps, candidates quasi-random Gaussians are drawn and scored in O(k²) by the incremental whitened eigensolver; the best admissible one is committed. candidates = 1 is the accept-first strategy. scale = :auto resolves via default_scale from the system's masses.

Refine(sweeps; candidates = 25, scale = :auto, sampler = HaltonSample(), indep_tol = 1e-4)

Suzuki–Varga cyclic refinement (Sect. 4.2.6, steps r1–r4): revisit each basis function in turn, draw candidates replacements, keep the best of {current, candidates}. Requires an existing basis (init = or a pipeline).

GVM([basis]; scale = nothing, optimizer = LBFGS(maxiter = 500, gradtol = 1e-6))

Joint gradient optimisation of all Gaussian parameters (widths via log-Cholesky encoding, plus shifts) using ForwardDiff/Hellmann–Feynman gradients. optimizer accepts an OptimKit LBFGS, ConjugateGradient, or GradientDescent instance and owns settings such as maxiter, gradtol, and verbosity. A cold start requires basis; a warm start infers it from init when omitted. scale controls only cold-start sampling and must be omitted for warm starts.

DynamicGVM(basis; candidates = 10, scale = :auto,
           optimizer = LBFGS(maxiter = 100, gradtol = 1e-6))

Per-step selection followed by joint gradient optimisation of the whole current basis (SVM-style sequential growth). basis is the final basis size, including any functions supplied through init. optimizer accepts the same OptimKit algorithm instances as GVM and is reused at every growth step.

Pipeline(stages)
alg₁ → alg₂ → alg₃

Composition of methods run left to right; each stage warm-starts from the previous stage's result. Built with the → operator (\to<tab>).

alg₁ → alg₂

Compose two solver methods into a left-to-right Pipeline.

Solution

Result of solve. Fields: E (eigenvalues of the final basis, ascending), basis::BasisSet, coefficients (generalized eigenvectors, cᵀSc = I), operators, state (target eigenstate), stages (length 1 unless a Pipeline ran), convergence (final report). sol.E₀ is the target-state energy E[state].

ConvergenceReport

What a solver run can honestly certify.

• converged::Bool
• criterion::Symbol — :saturation (stochastic: ΔE over the last window additions below tol), :stationarity (gradient tolerance met), :max_steps, or :early_stop
• ΔE::Float64 — tail energy change (Ha)
• tol::Float64, window::Int (0 for gradient methods)
• gradnorm — final gradient norm (nothing for stochastic methods)
• cond_S::Float64 — final overlap condition number
• notes::Vector{String} — caveats and early-stop explanations

StageResult(method, history, report)

One pipeline stage: the method that ran, its per-step target-state energy history, and its convergence report.

converged(sol::Solution) -> Bool
converged(report::ConvergenceReport) -> Bool

energy(sol::Solution; state = sol.state) -> Float64

Eigenvalue state of the final basis, sol.E[state]. With the default state this equals sol.E₀, the energy of the state the solver targeted.

energy_history(sol::Solution)            -> Vector{Float64}
energy_history(sol::Solution, i::Integer)

Per-step target-state energy history — concatenated across stages, or of stage i. Ready for plotting (see also plot(sol)). For a single stage, returns its stored energy history without copying. The eigenvalues of the final basis are sol.E; see energy for a single one.

convergence(sol::Solution) -> (steps, history)

Return the cumulative solver-step indices 1:length(energy_history(sol)) together with the per-step target-state energy history = energy_history(sol), ready for plotting a convergence curve. See also energy_history and plot(sol).

wavefunction(sol::Solution; state = sol.state) -> Wavefunction

Build the callable Wavefunction for the given state from a Solution's basis and generalized-eigenvector coefficients.

Wavefunction

Callable variational wavefunction ψ(r) = Σᵢ cᵢ gᵢ(r) in Jacobi coordinates (mass-weighted: the package's Jacobi transform normalises each relative coordinate by √μ — see FewBodyECG.jacobi_transform). Obtained from wavefunction; plot with plot(ψ; coord = i) or sample with radial_profile.

ψ(r) accepts either an N × 3 matrix of Cartesian positions (row = Jacobi coordinate, column = x, y, z) or a length-N vector, which places every Jacobi coordinate on the z axis: ψ(v) == ψ([0 0 v[1]; …]).

radial_profile(ψ::Wavefunction; coord = 1, direction = (0, 0, 1),
               rmax = 10.0, npoints = 400, normalize = true)

Sample the radial density r²|ψ|² along Jacobi coordinate coord on the physical half-line r ≥ 0, placing that coordinate at r * direction (the Cartesian direction is normalized to unit length) and holding the other coordinates at zero. Returns (r, density). When normalize = true the density is scaled so that its trapezoidal integral over [0, rmax] equals 1.

Polarized basis functions are not spherically symmetric: for example an xy-polarized Rank2Gaussian vanishes along the z axis, so choose a direction such as (1, 1, 0) for it.

Because r²|ψ|² is defined only for non-negative radial distance, no mirrored negative-r branch is produced.

FewBodyECG.build_hamiltonian_matrix(basis, operators)

Return the Hamiltonian matrix assembled from all operator matrix elements over basis. operators may be an Operators builder or a vector of operator terms.

FewBodyECG.build_overlap_matrix(basis)

Return the ECG overlap matrix S with entries <g_i|g_j> for a BasisSet.

FewBodyECG.solve_generalized_eigenproblem(H, S; max_condition=1e12, regularization=0)

Solve the symmetric generalized eigenproblem H*c = E*S*c, returning eigenvalues and S-orthonormal eigenvectors.

FewBodyECG.Λ(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.

FewBodyECG.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.

FewBodyECG.default_scale(masses)

Return the default Gaussian length scale inferred from the lightest finite particle mass in atomic units.

FewBodyECG.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)

FewBodyECG.up :: SpinProjection

Spin-½ projection eigenstate with eigenvalue +½.

FewBodyECG.down :: SpinProjection

Spin-½ projection eigenstate with eigenvalue −½.