API Reference

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

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.

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 (−½).

up :: SpinProjection

Spin-½ projection eigenstate with eigenvalue +½.

down :: SpinProjection

Spin-½ projection eigenstate with eigenvalue −½.

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 for compatibility and mapped to the z component.

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.

solve(ops, alg::Method = 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, Variational, GrowVariational, 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.

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

Variational(basis; scale = :auto, maxiter = 500, gtol = 1e-6, gradient = AutoDiff())

Joint LBFGS optimisation of all Gaussian parameters (widths via log-Cholesky encoding, plus shifts). Cold-starts from a quasi-random basis unless solve(...; init = sol) provides one.

GrowVariational(basis; candidates = 10, scale = :auto, maxiter_step = 100, gtol = 1e-6)

Per-step selection followed by joint LBFGS of the whole current basis (SVM-style sequential growth).

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.

AutoDiff()

ForwardDiff-based gradient backend (Hellmann–Feynman theorem).

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, energies, report)

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

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

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

Per-step target-state energy history — concatenated across stages, or of stage i. Ready for plotting (see also plot(sol)).

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

Return the cumulative solver-step indices 1:length(energies(sol)) together with the per-step target-state energy history = energies(sol), ready for plotting a convergence curve. See also energies 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 jacobi_transform). Obtained from wavefunction; plot with plot(ψ; coord = i) or sample with radial_profile.

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

Sample the radial density r²|ψ(r)|² along Jacobi coordinate coord on the physical half-line r ≥ 0 (the other coordinates held at zero), returning (r, density). When normalize = true the density is scaled so that its trapezoidal integral over [0, rmax] equals 1.

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

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.

build_overlap_matrix(basis)

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

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.

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

default_scale(masses)

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