Theory

FewBodyECG.jl solves non-relativistic few-body Hamiltonians by expanding the wave function in explicitly correlated Gaussians and evaluating all matrix elements analytically. The derivation below follows the shifted-Gaussian generating-function route used for rank-0, rank-1, and rank-2 matrix elements.

Variational reduction

We start from the Schrödinger equation

\[\hat{H}|\psi\rangle = \epsilon|\psi\rangle.\]

Approximate the wave function in a finite, non-orthogonal basis:

\[|\psi\rangle = \sum_{i=1}^{n} c_i |G_i\rangle.\]

Multiplying from the left by $\langle G_k|$ gives

\[\sum_i \langle G_k|\hat{H}|G_i\rangle c_i = \epsilon \sum_i \langle G_k|G_i\rangle c_i.\]

In matrix form,

\[Hc = \epsilon N c, \qquad H_{ki} = \langle G_k|\hat{H}|G_i\rangle, \qquad N_{ki} = \langle G_k|G_i\rangle.\]

The lowest eigenvalue is a variational upper bound for the represented Hamiltonian. Improving the basis can only lower or leave unchanged the best energy.

Coordinates and Hamiltonian

The package removes center-of-mass motion by transforming particle coordinates to Jacobi coordinates $\mathbf{x}$. In those coordinates the Hamiltonian has the form

\[\hat{H} = -\frac{\partial}{\partial \mathbf{x}} \Lambda \frac{\partial}{\partial \mathbf{x}^{T}} + \sum_{i<j}\frac{Z_i Z_j}{|\omega_{ij}^{T}\mathbf{x}|} + \sum_\alpha V_\alpha \exp\left[-\gamma_\alpha(\omega_\alpha^T\mathbf{x})^2\right].\]

Here $\Lambda$ is the reduced mass matrix in Jacobi coordinates and each $\omega$ vector selects one pair distance from $\mathbf{x}$.

Shifted Gaussian generator

The shifted correlated Gaussian is

\[\langle \mathbf{x}|A,\mathbf{a}\rangle = \exp\left(-\mathbf{x}^{T}A\mathbf{x}+\mathbf{a}^{T}\mathbf{x}\right),\]

where $A$ is symmetric positive definite and $\mathbf{a}$ is a column of shift vectors. The positive definiteness makes the basis function square integrable; the off-diagonal entries of $A$ encode pair correlations.

For a ket $(A,\mathbf{a})$ and bra $(B,\mathbf{b})$ define

\[C = A + B, \qquad R = C^{-1}, \qquad \mathbf{v} = \mathbf{a} + \mathbf{b}.\]

Their product is another Gaussian, so completing the square gives the overlap

\[M(C,\mathbf{v}) = \langle B,\mathbf{b}|A,\mathbf{a}\rangle = \left(\frac{\pi^n}{\det C}\right)^{3/2} \exp\left(\frac{1}{4}\mathbf{v}^{T}R\mathbf{v}\right).\]

This shifted overlap is the generating function. Polynomial prefactors and position moments are obtained by differentiating it with respect to the shift vectors and then setting shifts to zero.

Shifted operator matrix elements

The shifted kinetic-energy matrix element can be written compactly as

\[\langle B,\mathbf{b}|\hat{K}|A,\mathbf{a}\rangle = M(C,\mathbf{v}) \left[ 6\operatorname{Tr}(B\Lambda A R) + (\mathbf{b}-BR\mathbf{v})^T \Lambda (\mathbf{a}-AR\mathbf{v}) \right].\]

For a Coulomb term selected by $\omega$, set

\[\beta = \frac{1}{\omega^T R\omega}, \qquad \rho = \frac{1}{2}\omega^T R\mathbf{v}.\]

Then

\[\left\langle B,\mathbf{b}\left| \frac{1}{|\omega^T\mathbf{x}|} \right|A,\mathbf{a}\right\rangle = M(C,\mathbf{v})\,\frac{\operatorname{erf}(\sqrt{\beta}\rho)}{\rho},\]

with the finite zero-shift limit

\[2\sqrt{\frac{\beta}{\pi}}\,M(C,0).\]

The physical charge factor $Z_iZ_j$ is multiplied in by CoulombOperator. A Gaussian pair potential is even simpler:

\[V_0 e^{-\gamma(\omega^T\mathbf{x})^2} \quad\Longrightarrow\quad C \mapsto C + \gamma\,\omega\omega^T.\]

Its matrix element is $V_0$ times the overlap evaluated with that modified quadratic form.

Rank-0 Gaussians

Rank-0, or s-wave, Gaussians are the zero-shift limit:

\[\langle \mathbf{x}|A\rangle = \lim_{\mathbf{a}\to 0}\langle \mathbf{x}|A,\mathbf{a}\rangle = \exp(-\mathbf{x}^{T}A\mathbf{x}).\]

The rank-0 overlap is

\[M_0 = \langle B|A\rangle = \left(\frac{\pi^n}{\det(A+B)}\right)^{3/2}.\]

The corresponding kinetic and Coulomb elements are

\[\langle B|\hat{K}|A\rangle = 6M_0\,\operatorname{Tr}(B\Lambda A(A+B)^{-1}),\]

and

\[\left\langle B\left|\frac{1}{|\omega^T\mathbf{x}|}\right|A\right\rangle = 2\sqrt{\frac{\beta}{\pi}}\,M_0, \qquad \beta = \frac{1}{\omega^T(A+B)^{-1}\omega}.\]

These are the formulas used for the high-level stochastic and variational solvers, which currently sample rank-0 bases.

Tensor prefactors

Rank-1 and rank-2 Gaussians are generated by taking Taylor coefficients of the shifted Gaussian:

\[\langle \mathbf{x}|(\mathbf{u})A\rangle = (\mathbf{u}^T\mathbf{x})e^{-\mathbf{x}^{T}A\mathbf{x}},\]

\[\langle \mathbf{x}|(\mathbf{u}\mathbf{v})A\rangle = (\mathbf{u}^T\mathbf{x})(\mathbf{v}^T\mathbf{x}) e^{-\mathbf{x}^{T}A\mathbf{x}}.\]

For example, the rank-1 overlap comes from the $O(\mathbf{u}\mathbf{v})$ term in the shifted overlap:

\[\langle (\mathbf{v})B|(\mathbf{u})A\rangle = \frac{1}{2}\mathbf{v}^{T}(A+B)^{-1}\mathbf{u}\,M_0.\]

Kinetic and Coulomb rank-1/rank-2 formulas are obtained the same way: expand the shifted matrix element, keep the coefficient with the required shift order, and set the remaining shifts to zero. This is why one shifted formula can generate the s-, p-, and d-wave matrix elements used by the power-user matrix layer.

Basis construction

Once the matrix elements are analytic, the numerical problem is choosing a useful basis:

SVM draws quasi-random candidates and keeps the one that lowers the target eigenvalue.
Refine revisits existing basis slots and tries replacements.
Variational jointly optimizes all rank-0 parameters with LBFGS.
GrowVariational alternates growth and continuous optimization.

Stochastic methods are cheap and robust but can saturate under a fixed sampling scale. Gradient methods cost more but move the Gaussian parameters continuously after the sampled basis has found the right region.

Hydrogen check

For hydrogen in Hartree units,

\[\hat{H} = -\frac{1}{2}\nabla^2 - \frac{1}{r}, \qquad E_n = -\frac{1}{2n^2}.\]

Rank-0, rank-1, and rank-2 Gaussian bases target the lowest s-, p-, and d-wave states. The corresponding exact energies are

\[E_{1s}=-\frac{1}{2}, \qquad E_{2p}=-\frac{1}{8}, \qquad E_{3d}=-\frac{1}{18}.\]

The hydrogen example compares these values against the analytical energies reported by Antique.jl.

References

D. V. Fedorov, A. F. Teilmann, M. C. Østerlund, and T. L. Norrbohm, "Explicitly Correlated Gaussians with Tensor Pre-factors: Analytic Matrix Elements," Few-Body Systems 65, 75 (2024). doi:10.1007/s00601-024-01945-x.
Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems, Springer, 1998.
Antique.jl, analytical solutions for solvable quantum-mechanical models: github.com/ohno/Antique.jl.