API reference

BVM(; kwargs...)

Alias for BayesianVariationalMethod.

Basis is an abstract type.

BasisSet(basis1, basis2, ...)

\[\{ \phi_1, \phi_2, \phi_3, \cdots \}\]

The basis set is the input for Rayleigh-Ritz method. You can define the basis set like this:

The concrete element type is preserved when all basis functions have a common type, avoiding abstract Basis dispatch in matrix-construction loops.

\[\begin{aligned} \phi_1(r) &= \exp(-13.00773 ~r^2), \\ \phi_2(r) &= \exp(-1.962079 ~r^2), \\ \phi_3(r) &= \exp(-0.444529 ~r^2), \\ \phi_4(r) &= \exp(-0.1219492 ~r^2). \end{aligned}\]

BS = BasisSet(
  SimpleGaussianBasis(13.00773),
  SimpleGaussianBasis(1.962079),
  SimpleGaussianBasis(0.444529),
  SimpleGaussianBasis(0.1219492),
)

BayesianVariationalMethod(; max_basis, pool_size=100, tuple_size=3,
                          initial_samples=100, batch_size=50, rounds=4,
                          search_size=10_000, beta=2.0, abstol=1e-8,
                          patience=3, overlap_tol=1e-10)

Configure Bayesian selection of fixed-size groups from a finite candidate basis. max_basis limits the accepted basis dimension; pool_size is the temporary candidate pool; tuple_size is the number of functions adopted per outer step; initial_samples, batch_size, rounds, and search_size control the Gaussian-process search budget; and beta weighs posterior uncertainty in the lower-confidence-bound acquisition. abstol is an absolute energy improvement in the Hamiltonian's units, patience counts consecutive negligible steps, and overlap_tol rejects nearly linearly dependent subsets.

ComplexGaussianBasis(a=1, ω=1, component=:cos, l=0, m=0)

One normalized real component of a complex-range Gaussian pair:

\[\phi^{\cos}_{nlm}=N^{\cos}_{nl}r^l e^{-a_n r^2} \cos(\omega a_n r^2)Y_l^m(\hat{\boldsymbol r}),\]

or the corresponding sine function when component=:sin. Equivalently, the pair is formed from exponents (1±iω)aₙ. Use ComplexGaussianBasisSet(r₁, rₙ, n; ω, l, m) to construct both members at every geometrically spaced range.

The analytic overlap, kinetic, constant, power-law, Coulomb, linear, and Gaussian-potential matrix elements are supported by the ordinary Rayleigh–Ritz solve interface.

References: E. Hiyama, Y. Kino, and M. Kamimura, Prog. Part. Nucl. Phys. 51, 223–307 (2003), and E. Hiyama and M. Kamimura, Front. Phys. 13, 132106 (2018).

Constant(constant=1)

\[+ c\]

ContractedBasis([c1, c2, ...], [primitive1, primitive2, ...])

\[\phi' = \sum_i c_i \phi_i\]

A contracted basis is a nonempty linear combination of PrimitiveBasis objects. The numbers of coefficients and primitive bases must match. Numeric coefficient types are promoted to a common type.

The components are stored in tuples, making their number part of the type and preserving the concrete type of every primitive basis. This lets Julia specialize and inline evaluation without dispatching through an abstract PrimitiveBasis container. Tuple inputs are recommended when constructing a basis in performance-sensitive code; vectors are also accepted and converted once during construction.

The fields are named coefficients and primitives. The aliases c and φ are retained for compatibility.

Coulomb(coefficient=1)

\[+ \frac{a}{r}\]

Custom(f)

\[+ f(r)\]

DatabaseEntry(hamiltonian, energy)

A benchmark problem stored in the database. hamiltonian is ready to be passed to a solver and energy is its reference energy.

Delta(coefficient=1)

\[+ a δ(r)\]

Exponential(coefficient=1, exponent=1)

\[+ a \exp(- b r)\]

FiniteDifferenceMethod(Δr=0.1, rₘₐₓ=50.0, R=Δr:Δr:rₘₐₓ, l=0, direction=:c, solver=:LinearAlgebra)

Gaussian(coefficient=1, exponent=1)

\[+ a \exp(- b r^2)\]

GaussianBasis(a=1, l=0, m=0)

\[\phi_{ilm}(r, θ, φ) = N _{il} r^l \exp(-a_i r^2) Y_l^m(θ, φ)\]

This normalized primitive is used by the Gaussian expansion method (GEM). N_{il} normalizes the primitive, while l and m specify its spherical harmonic. For a central Hamiltonian, combine primitives with the same l and m and pass the basis set to the existing Rayleigh–Ritz solver.

For example, a $p$-wave basis can be constructed from geometrically spaced exponents as follows:

exponents = TwoBody.geometric(0.1, 10.0, 20)
BS = BasisSet((GaussianBasis(a; l=1, m=0) for a in exponents)...)

GeometricBasisSet(basistype, r₁, rₙ, n; nₘᵢₙ=1, nₘₐₓ=n)

This is a basis set with exponents generated by geometric(r₁, rₙ, n; nₘₐₓ=n, nₘᵢₙ=1). You can define the 20 real-range Gaussian basis functions used for the hydrogen atom in Appendix A.2 and Table VII of E. Hiyama, M. Kamimura, Front. Phys. 13, 132106 (2018) like this:

\[ r_1 = 0.1, r_{n_\mathrm{max}} = 80.0, n_\mathrm{max} = 20.\]

BS = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20)

Hamiltonian(operator1, operator2, ...)

\[\hat{H} = \sum_i \hat{o}_i\]

The Hamiltonian is the input for each solver. This is an example for the non-relativistic Hamiltonian of hydrogen atom in atomic units:

\[\hat{H} = - \frac{1}{2} \nabla^2 - \frac{1}{r}\]

H = Hamiltonian(
  Kinetic(hbar = 1, m = 1),
  Coulomb(coefficient = -1),
)

Kinetic(hbar=1, m=1)

\[-\frac{\hbar^2}{2m} \nabla^2\]

KineticTerm <: Operator is an abstract type.

Laplacian(coefficient=1)

\[+ a\nabla^2\]

Linear(coefficient=1)

\[+ ar\]

Operator is an abstract type.

PotentialTerm <: Operator is an abstract type.

PowerLaw(coefficient=1, exponent=1)

\[+ ar^n\]

PowerSlaterBasis(n=0, a=1)

An unnormalized power-Slater basis function for an s wave:

\[\phi_{n,a}(r) = r^n \exp(-ar).\]

n is an integer power and a is the exponential parameter. For a square-integrable basis function, choose parameters for which the required matrix elements are finite (normally $n \ge 0$ and $a > 0$).

PrimitiveBasis <: Basis is an abstract type.

QuanticsTensorTrainMethod(; quantics=10, r₀=0.0, rₘₐₓ=50.0, l=0,
                           tolerance=1e-10, maxbonddim=64,
                           maxoperatorbonddim=128, maxiter=100, sweeps=4,
                           deflation_shift=10.0)

Configure a radial QTT grid with $2^{\mathtt{quantics}}$ uniformly spaced interior points. The excluded endpoints r₀ and rₘₐₓ impose zero-value (Dirichlet) boundaries. maxbonddim and maxoperatorbonddim bound state and MPO ranks; deflation_shift is the penalty for each previously computed state.

RelativisticCorrection(c=1, m=1, n=2) The p^{2n} term of the Taylor expansion:

\[\begin{aligned} \sqrt{p^2 c^2 + m^2 c^4} =& m \times c^2 \\ &+ 1 / 2 / m \times p^2 (n=1) \\ &- 1 / 8 / m^3 / c^2 \times p^4 (n=2) \\ &+ 1 / 16 / m^5 / c^4 \times p^6 (n=3) \\ &- 5 / 128 / m^7 / c^6 \times p^8 (n=4) \\ &+ \cdots \end{aligned}\]

Use c = 137.035999177 (from 2022 CODATA) in the atomic units.

RelativisticKinetic(c=1, m=1)

\[\sqrt{p^2 c^2 + m^2 c^4} - m c^2\]

Use c = 137.035999177 (from 2022 CODATA) in the atomic units.

RestEnergy(c=1, m=1)

\[m c^2\]

Use c = 137.035999177 (from 2022 CODATA) in the atomic units.

SimpleGaussianBasis(a=1)

Position-Space

\[\phi_i(\pmb{r}) = \exp(-a_i r^2)\]

Momentum-Space

\[\phi_{i}(\pmb{k}) = \frac{1}{(2a_i)^{\frac{3}{2}}} \exp(-k^2/4a_i)\]

Proof (Fourier Transform)

\[\begin{aligned} \phi_{n}(\pmb{k}) &= \frac{1}{\sqrt{2 \pi}^3} \int \phi_{n}(\pmb{r}) \mathrm{e}^{\mathrm{i} \pmb{k} \cdot \pmb{r}} \mathrm{d}\pmb{r} \\ &= \frac{1}{\sqrt{2 \pi}^3} \int \phi_{n}(\pmb{r}) \mathrm{e}^{\mathrm{i} \pmb{k} \cdot \pmb{r}} r^2 \sin (\theta) ~\mathrm{d}r \mathrm{d}\theta \mathrm{d} \varphi \\ &= \frac{1}{\sqrt{2 \pi}^3} \iiint \mathrm{e}^{-\alpha_i r^2} \sqrt{4\pi} Y_{00}(\hat{\pmb{r}}) \left[ 4 \pi \sum_{l'=0}^{\infty} \sum_{m=-l'}^{l'} \mathrm{i}^{l'} j_{l'}(pr) Y_{l'm'}(\hat{\pmb{k}}) Y_{l'm'}^*(\hat{\pmb{r}}) \right] r^2 \sin\theta~ \mathrm{d} r \mathrm{d} \theta \mathrm{d} \varphi \\ &= \frac{1}{\sqrt{2 \pi}^3} 4 \pi \sqrt{4\pi} \sum_{l'=0}^{\infty} \sum_{m=-l'}^{l'} \left[ \mathrm{i}^{l'} Y_{l'm'}(\hat{\pmb{k}}) \int_0^{2 \pi} \int_0^\pi Y_{00}(\hat{\pmb{r}}) Y_{l'm'}^*(\hat{\pmb{r}}) \sin (\theta)~ \mathrm{d} \theta \mathrm{d} \varphi \int_0^{\infty} j_{l'}(pr) \mathrm{e}^{-\alpha_i r^2} r^{2} \mathrm{d}r \right]\\ &= \frac{1}{\sqrt{2 \pi}^3} 4 \pi \sqrt{4\pi} \sum_{l'=0}^{\infty} \sum_{m=-l'}^{l'} \left[ \mathrm{i}^{l'} Y_{l'm'}(\hat{\pmb{k}}) \delta_{0l'} \delta_{0m'} \int_0^{\infty} j_{l'}(kr) \mathrm{e}^{-\alpha_i r^2} r^{2} \mathrm{d}r \right] \\ &= \frac{1}{\sqrt{2 \pi}^3} 4 \pi \sqrt{4\pi} \mathrm{i}^{0} Y_{00}(\hat{\pmb{k}}) \int_0^{\infty} j_{0}(kr) \mathrm{e}^{-\alpha_i r^2} r^{2} \mathrm{d}r \\ &= \frac{1}{2\pi\sqrt{2\pi}} 4 \pi \frac{\sqrt{4\pi}}{\sqrt{4\pi}} \sqrt{\frac{\pi}{2}} \sqrt{\frac{2}{\pi}} \int_0^{\infty} j_{0}(kr) \mathrm{e}^{-\alpha_i r^2} r^{2} ~\mathrm{d}r \\ &= \frac{1}{(2\alpha_i)^{\frac{3}{2}}} \mathrm{e}^{-\frac{k^2}{4 \alpha_i}} \end{aligned}\]

Formula

plane-wave expansion in spherical harmonics:

\[\mathrm{e}^{\mathrm{i} \pmb{k} \cdot \pmb{r}} = 4 \pi \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \mathrm{i}^{l} j_{l}(pr) Y_{lm}(\hat{\pmb{k}}) Y_{lm}^*(\hat{\pmb{r}})\]

special case of spherical harmonics:

\[Y_{00}(\hat{\pmb{r}}) = \frac{1}{\sqrt{4\pi}}\]

orthonormality of spherical harmonics:

\[\int_0^{2\pi} \int_0^\pi Y_{lm}(\hat{\pmb{r}})^* Y_{l'm'}(\hat{\pmb{r}}) \sin(\theta) ~ \mathrm{d} \theta \mathrm{d} \varphi = \delta_{ll'} \delta_{mm'}\]

citation needed:

\[\sqrt{\frac{2}{\pi}} \int r^{l} j_l(kr) \mathrm{e}^{-\alpha r^2} r^{2} \mathrm{d} r = \frac{1}{(2\alpha)^{l+\frac{3}{2}}} k^l e^{-\frac{k^2}{4\alpha}}\]

Tabulated(R, V)

VariationalMonteCarlo(n_steps=10^5, burn_in=10^3, thinning=1, n_walkers=1, δ=0.5, r₀=[1.0, 0.0, 0.0])

Options for variational Monte Carlo with a symmetric, uniform Metropolis proposal. The sampler targets $|\psi(\mathbf{r})|^2$. n_steps is the number of retained samples per walker, burn_in is the number of initial transitions discarded from each walker, thinning is the number of transitions between retained samples, n_walkers is the number of Markov chains, δ is the proposal-box width, and r₀ is the initial position of every walker. Each walker performs burn_in + n_steps * thinning transitions.

VariationalNeuralNetwork(; fdm=nothing, Δr=nothing, rₘₐₓ=nothing,
  R=nothing, l=nothing, direction=nothing, solver=nothing,
  architecture=[2], activation=softplus, init=nothing, optimizer=nothing,
  maxiters=1000, abstol=1e-8, patience=10, every=100)

Options for optimizing a Lux neural network as a radial trial wavefunction. VNN is an abbreviation for VariationalNeuralNetwork. VNN support is activated by loading Lux.jl, Optimisers.jl, and Zygote.jl. If init or optimizer is nothing, the extension uses Lux.glorot_normal or Optimisers.Adam(0.01), respectively.

Yukawa(coefficient=1, exponent=1)

\[+ \frac{a}{r} \exp(- b r)\]

put!(key, hamiltonian, energy)

Add a benchmark problem to the database. key may be a Symbol or string, hamiltonian must be a Hamiltonian, and energy must be real. Registering the same key twice throws an ArgumentError.

ComplexGaussianBasisSet(r₁, rₙ, n; ω=1.0, l=0, m=0)

Construct the 2n real, normalized complex-range Gaussian primitives

\[r^l e^{-\nu_j r^2}\cos(\omega\nu_jr^2),\qquad r^l e^{-\nu_j r^2}\sin(\omega\nu_jr^2),\]

where νⱼ = 1/rⱼ² and the ranges from r₁ through rₙ form a geometric progression. This is the real cos/sin form of the complex-range GEM basis introduced by Hiyama, Kino, and Kamimura (2003) and used for the highly excited hydrogen example of Hiyama and Kamimura (2018).

FC(hamiltonian, basis; g=PowerSlaterBasis(1, 0))

Generate the next Free Complement basis from a PowerSlaterBasis or a BasisSet of power-Slater functions by collecting the basis-function forms in $g(H-E_n)\phi$. Numerical coefficients are ignored. The default scaling function is $g(r)=r$. Duplicate functions and functions singular at the origin are removed.

The Hamiltonian may contain Kinetic, Laplacian, RestEnergy, Constant, Linear, Coulomb, integer-exponent PowerLaw, Exponential, and Yukawa terms. An ArgumentError is thrown when an operator does not map a PowerSlaterBasis to power-Slater functions.

For the hydrogen Hamiltonian, repeated application starting from PowerSlaterBasis(0, 1.5) adds one power of $r$ at each iteration.

db(key::Union{Symbol,AbstractString}) -> DatabaseEntry

Return the benchmark Hamiltonian and reference energy associated with key. The returned Hamiltonian is independent of the stored value and can safely be modified by callers.

Examples

entry = db(:hydrogen)
result = solve(entry.hamiltonian, method)
isapprox(result.values[1], entry.energy)

dbkeys() -> Vector{Symbol}

Return the available database keys in deterministic order.

element(o::Constant, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | c | \phi_{j} \rangle &= c \langle \phi_{i} | \phi_{j} \rangle \\ &= c \iiint \phi_{i}^*(r) \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= c \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^{2} \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= c \times 2\pi \times 2 \times \frac{1!!}{2^{2}} \sqrt{\frac{\pi}{a^{3}}} \\ &= \underline{c \left( \frac{\pi}{\alpha_i + \alpha_j} \right)^{3/2}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{a^{2n+1}}}\]

element(o::Coulomb, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | \frac{1}{r} | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \times \frac{1}{r} \times \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= 2\pi \times 2 \times \frac{0!}{2 (\alpha_i + \alpha_j)} \\ &= \underline{\frac{2\pi}{\alpha_i + \alpha_j}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n+1} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{n!}{2 a^{n+1}}\]

element(o::Custom, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

The matrix element for a custom central potential is evaluated numerically with adaptive Gauss–Kronrod quadrature:

\[\langle \phi_i | V | \phi_j \rangle = 4\pi \int_0^\infty r^2 \phi_i(r) V(r) \phi_j(r)\,\mathrm{d}r.\]

element(o::Gaussian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | \exp(-br^2) | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \times \exp(-br^2) \times \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^2 \mathrm{e}^{-(b+\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= 2\pi \times 2 \times \frac{1!!}{2^{2}} \sqrt{\frac{\pi}{(b + \alpha_i + \alpha_j)^{2\cdot1+1}}} \\ &= \underline{\left( \frac{\pi}{b + \alpha_i + \alpha_j} \right)^{3/2}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{a^{2n+1}}}\]

element(o::Hamiltonian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} H_{ij} &= \langle \phi_{i} | \hat{H} | \phi_{j} \rangle \\ &= \langle \phi_{i} | \sum_k \hat{o}_k | \phi_{j} \rangle \\ &= \sum_k \langle \phi_{i} | \hat{o}_k | \phi_{j} \rangle \\ \end{aligned}\]

element(o::Kinetic, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

Derivation (without Green's identity)

\[\begin{aligned} T_{ij} = \langle \phi_{i} | \hat{T} | \phi_{j} \rangle &= \iiint \mathrm{e}^{-\alpha_i r^2} \left[ -\frac{\hbar^2}{2\mu} \nabla^2 \right] \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= -\frac{\hbar^2}{2\mu} \iiint \mathrm{e}^{-\alpha_i r^2} \left[ \nabla^2 \right] \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= -\frac{\hbar^2}{2\mu} \iiint \mathrm{e}^{-\alpha_i r^2} \left[ -6\alpha_j + 4\alpha_j^2 r^2 \right] \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= -\frac{\hbar^2}{2\mu} \iint \sin\theta ~\mathrm{d}\theta \mathrm{d}\varphi \int \left[ -6\alpha_j + 4\alpha_j^2 r^2 \right] r^2 \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= -\frac{\hbar^2}{2\mu} \cdot 4\pi \left[ -6\alpha_j \mathrm{GGI}(2, \alpha_i + \alpha_j) +4\alpha_j^2 \mathrm{GGI}(4, \alpha_i + \alpha_j) \right] \\ &= -\frac{\hbar^2}{2\mu} \cdot 4\pi \left[ -6\alpha_j \frac{\Gamma\left( \frac{3}{2} \right)}{2 (\alpha_i + \alpha_j)^{\frac{3}{2}}} +4\alpha_j^2 \frac{\Gamma\left( \frac{5}{2} \right)}{2 (\alpha_i + \alpha_j)^{\frac{5}{2}}} \right] \\ &= -\frac{\hbar^2}{2\mu} \cdot 4\pi \left[ -6\alpha_j \frac{ \sqrt{\pi}/2}{2 (\alpha_i + \alpha_j)^{\frac{3}{2}}} +4\alpha_j^2 \frac{3\sqrt{\pi}/4}{2 (\alpha_i + \alpha_j)^{\frac{5}{2}}} \right] \\ &= -\frac{\hbar^2}{2\mu} \cdot 4\pi \left[ \frac{\alpha_j}{\alpha_i + \alpha_j} - 1 \right] \cdot 6 \alpha_j \cdot \frac{\sqrt{\pi}/2}{2 (\alpha_i + \alpha_j)^{\frac{3}{2}}} \\ &= -\frac{\hbar^2}{2\mu} \cdot 4\pi \left[ - \frac{\alpha_i}{\alpha_i + \alpha_j} \right] \cdot 6 \alpha_j \cdot \frac{\sqrt{\pi}/2}{2 (\alpha_i + \alpha_j)^{\frac{3}{2}}} \\ &= \underline{ \frac{\hbar^2}{2\mu} \cdot 6 \cdot \frac{\alpha_i \alpha_j \pi^{\frac{3}{2}}}{(\alpha_i + \alpha_j)^{\frac{5}{2}}} } \end{aligned}\]

Derivation (with Green's identity)

\[\begin{aligned} T_{ij} = \langle \phi_{i} | \hat{T} | \phi_{j} \rangle &= \iiint \mathrm{e}^{-\alpha_i r^2} \left[ -\frac{\hbar^2}{2\mu} \nabla^2 \right] \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= -\frac{\hbar^2}{2\mu} \iiint \mathrm{e}^{-\alpha_i r^2} \nabla^2 \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \frac{\hbar^2}{2\mu} \iiint \left[ \nabla \mathrm{e}^{-\alpha_i r^2} \right] \left[ \nabla \mathrm{e}^{-\alpha_j r^2} \right] ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \frac{\hbar^2}{2\mu} \iiint \left[ -2 \alpha_i r \mathrm{e}^{-\alpha_i r^2} \right] \left[ -2 \alpha_j r \mathrm{e}^{-\alpha_j r^2} \right] ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \frac{\hbar^2}{2\mu} \cdot 4 \alpha_i \alpha_j \iiint \left[ r \mathrm{e}^{-\alpha_i r^2} \right] \left[ r \mathrm{e}^{-\alpha_j r^2} \right] ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \frac{\hbar^2}{2\mu} \cdot 4 \alpha_i \alpha_j \iint \sin\theta ~\mathrm{d}\theta \mathrm{d}\varphi \int r^4 \mathrm{e}^{- (\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= \frac{\hbar^2}{2\mu} \cdot 4 \alpha_i \alpha_j \cdot 4 \pi \cdot \mathrm{GGI}(4, \alpha_i + \alpha_j) \\ &= \frac{\hbar^2}{2\mu} \cdot 4 \alpha_i \alpha_j \cdot 4 \pi \cdot \frac{\Gamma\left( \frac{5}{2} \right)}{2 (\alpha_i + \alpha_j)^{\frac{5}{2}}} \\ &= \frac{\hbar^2}{2\mu} \cdot 4 \alpha_i \alpha_j \cdot 4 \pi \cdot \frac{3\sqrt{\pi}/4}{2 (\alpha_i + \alpha_j)^{\frac{5}{2}}} \\ &= \underline{ \frac{\hbar^2}{2\mu} \cdot 6 \cdot \frac{\alpha_i \alpha_j \pi^{\frac{3}{2}}}{(\alpha_i + \alpha_j)^{\frac{5}{2}}} } \end{aligned}\]

Formula

Green's first identity:

\[\begin{aligned} \iiint_V f \pmb{\nabla}^2 g ~ \mathrm{d}V + \iiint_V \pmb{\nabla} f \cdot \pmb{\nabla} g ~ \mathrm{d}V = \iint_{\partial V} f \pmb{\nabla} g \cdot \pmb{n} ~ \mathrm{d}S \end{aligned}\]

generalized Gaussian integral:

\[\begin{aligned} \mathrm{GGI}(n,b) = \int_0^{\infty} x^{n} \exp \left(-b x^2\right) \mathrm{d}x = \frac{\Gamma\left( \frac{n+1}{2} \right)}{2 b^{\frac{n+1}{2}}} \end{aligned}\]

element(o::Linear, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | r | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \times r \times \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^3 \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= 2\pi \times 2 \times \frac{1!}{2 (\alpha_i + \alpha_j)^{2}} \\ &= \underline{\frac{2\pi}{(\alpha_i + \alpha_j)^2}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n+1} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{n!}{2 a^{n+1}}\]

element(o::PowerLaw, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | r^n | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \times r^n \times \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^{n+2} \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= 2\pi \times 2 \times \frac{\Gamma\left( \frac{n+3}{2} \right)}{2 (\alpha_i + \alpha_j)^{\frac{n+3}{2}}} \\ &= \underline{2\pi\frac{\Gamma\left( \frac{n+3}{2} \right)}{(\alpha_i + \alpha_j)^{\frac{n+3}{2}}}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{n} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{\Gamma\left( \frac{n+1}{2} \right)}{2 a^{\frac{n+1}{2}}}\]

element(o::RestEnergy, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | mc^2 | \phi_{j} \rangle &= mc^2 \langle \phi_{i} | \phi_{j} \rangle \\ &= mc^2 \iiint \phi_{i}^*(r) \phi_{j}(r) ~r^2 \sin\theta ~\mathrm{d}r \mathrm{d}\theta \mathrm{d}\varphi \\ &= mc^2 \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^{2} \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= mc^2 \times 2\pi \times 2 \times \frac{1!!}{2^{2}} \sqrt{\frac{\pi}{a^{3}}} \\ &= \underline{mc^2 \left( \frac{\pi}{\alpha_i + \alpha_j} \right)^{3/2}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{a^{2n+1}}}\]

element(SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} S_{ij} = \langle \phi_{i} | \phi_{j} \rangle &= \int \phi_{i}^*(r) \phi_{j}(r) \mathrm{d} \pmb{r} \\ &= \iiint \mathrm{e}^{-\alpha_i r^2} \mathrm{e}^{-\alpha_j r^2} ~r^2 \sin\theta ~ \mathrm{d} r \mathrm{d} \theta \mathrm{d} \varphi \\ &= \int_0^{2\pi} \mathrm{d}\varphi \int_0^\pi \sin\theta ~\mathrm{d}\theta \int_0^\infty r^{2} \mathrm{e}^{-(\alpha_i + \alpha_j) r^2} ~\mathrm{d}r \\ &= 2\pi \times 2 \times \frac{1!!}{2^{2}} \sqrt{\frac{\pi}{a^{3}}} \\ &= \underline{\left( \frac{\pi}{\alpha_i + \alpha_j} \right)^{3/2}} \end{aligned}\]

Integral Formula:

\[\int_0^{\infty} r^{2n} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{(2n-1)!!}{2^{n+1}} \sqrt{\frac{\pi}{a^{2n+1}}}\]

element(o::Laplacian, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)

\[\begin{aligned} \langle \phi_{i} | \nabla^2 | \phi_{j} \rangle = \underline{ -6 \frac{\alpha_i \alpha_j \pi^{\frac{3}{2}}}{(\alpha_i + \alpha_j)^{\frac{5}{2}}} } \end{aligned}\]

geometric(r₁, rₙ, n::Int; nₘₐₓ::Int=n, nₘᵢₙ::Int=1)

Exponents of Gaussian basis functions are given by geometric progression:

\[\begin{aligned} & v_i = \frac{1}{r_i^2}, \\ & r_i = r_1 a^{i-1}. \end{aligned}\]

This function return array of $\nu_i$:

\[(r_1, r_{n}, n, n_\mathrm{max}) \mapsto (\nu_1, \nu_2, \cdots, \nu_{n-1}, \nu_n, \nu_{n+1}, \cdots, \nu_{n_\mathrm{max}})\]

Usually $n = n_\mathrm{max}$. Set $n<n_\mathrm{max}$ if you want to extend the geometric progression.

Examples:

julia> ν = TwoBody.geometric(0.1, 10.0, 5)
5-element Vector{Float64}:
 100.0
  10.0
   0.9999999999999997
   0.09999999999999996
   0.009999999999999995

julia> ν = TwoBody.geometric(0.1, 10.0, 5, nₘₐₓ = 10)
10-element Vector{Float64}:
 100.0
  10.0
   0.9999999999999997
   0.09999999999999996
   0.009999999999999995
   0.0009999999999999994
   9.999999999999994e-5
   9.999999999999992e-6
   9.999999999999991e-7
   9.999999999999988e-8

local_energy(hamiltonian, wavefunction, position)

Evaluate

\[E_\mathrm{loc}(\mathbf{r}) = \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}.\]

The Laplacian of a real-valued trial wavefunction is evaluated with ForwardDiff.hessian. Non-relativistic kinetic, Laplacian, rest-energy, and potential terms with a defined V method are supported.

matrix(basisset::BasisSet)

This function returns the overlap matrix $\pmb{S}$. The element is written as $S_{ij} = \langle \phi_{i} | \phi_{j} \rangle$.

matrix(hamiltonian::Hamiltonian, basisset::BasisSet)

This function returns the Hamiltonian matrix $\pmb{H}$. The element is written as $H_{ij} = \langle \phi_{i} | \hat{H} | \phi_{j} \rangle$.

matrix(o::Hamiltonian, method::FiniteDifferenceMethod)

The matrix for the Hamiltonian is a sum of matrices for each term,

\[\pmb{H} = \sum_i \pmb{O}_i.\]

matrix(o::Kinetic, method::FiniteDifferenceMethod)

We use the shorthand notation $\psi'(r) = \frac{\mathrm{d}\psi}{\mathrm{d}r}(r)$ and $\psi''(r) = \frac{\mathrm{d}^{2}\psi}{\mathrm{d}r^{2}}(r)$. For the uniform grid spacing ($r_{i+1} = r_{i} + \Delta r$), the finite difference for the first derivative,

\[\frac{\mathrm{d}\psi}{\mathrm{d}r}(r) = \frac{\psi(r+\Delta r) - \psi(r-\Delta r)}{2\Delta r} + O(\Delta r^{2})\]

is written as

\[\left(\begin{array}{ccccc} \psi'(r_1) \\ \psi'(r_2) \\ \psi'(r_3) \\ \psi'(r_4) \\ \vdots \end{array}\right) \simeq \frac{1}{2\Delta r} \left(\begin{array}{ccccc} 0 & 1 & 0 & 0 &\ldots \\ -1 & 0 & 1 & 0 &\ldots \\ 0 & -1 & 0 & 1 &\ldots \\ 0 & 0 & -1 & 0 &\ldots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) \left(\begin{array}{ccccccc} \psi(r_1) \\ \psi(r_2) \\ \psi(r_3) \\ \psi(r_4) \\ \vdots \end{array}\right),\]

and the finite difference for the second derivative,

\[\frac{\mathrm{d}^{2}\psi}{\mathrm{d}r^{2}}(r) = \frac{\psi(r+\Delta r) - 2\psi(r) + \psi(r-\Delta r)}{\Delta r^{2}} + O(\Delta r^{2}).\]

is written as

\[\left(\begin{array}{ccccc} \psi''(r_1) \\ \psi''(r_2) \\ \psi''(r_3) \\ \psi''(r_4) \\ \vdots \end{array}\right) \simeq \frac{1}{\Delta r^2} \left(\begin{array}{ccccccc} -2 & 1 & 0 & 0 & \ldots \\ 1 & -2 & 1 & 0 & \ldots \\ 0 & 1 & -2 & 1 & \ldots \\ 0 & 0 & 1 & -2 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array}\right) \left(\begin{array}{ccccccc} \psi(r_1) \\ \psi(r_2) \\ \psi(r_3) \\ \psi(r_4) \\ \vdots \end{array}\right).\]

Similarly, the matrix for the kinetic energy,

\[\hat{T} = -\frac{\hbar^2}{2\mu} \left[ \frac{\partial^2}{\partial r^2} + \frac{2}{r} \frac{\partial}{\partial r} - \frac{l(l+1)}{r^2} \right]\]

is written as

\[\pmb{T} = - \frac{\hbar^2}{2\mu} \left[ \frac{1}{{\Delta r}^2} \left(\begin{array}{ccccccc} -2 & 1 & 0 & \ldots \\ 1 & -2 & 1 & \ldots \\ 0 & 1 & -2 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) + \left(\begin{array}{ccccccc} 2/r_1 & 0 & 0 & \ldots \\ 0 & 2/r_2 & 0 & \ldots \\ 0 & 0 & 2/r_3 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) \frac{1}{2\Delta r} \left(\begin{array}{ccccccc} 0 & 1 & 0 & \ldots \\ -1 & 0 & 1 & \ldots \\ 0 & -1 & 0 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) - l(l+1) \left(\begin{array}{ccccccc} 1/{r_1}^2 & 0 & 0 & \ldots \\ 0 & 1/{r_2}^2 & 0 & \ldots \\ 0 & 0 & 1/{r_3}^2 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) \right].\]

matrix(operator::Operator, basisset::BasisSet)

This function is used for the expectation values and is not used in computing the Hamiltonian matrix.

This function returns the matrix corresponding to the operator in the given basis set. The element is written as $O_{ij} = \langle \phi_{i} | \hat{o} | \phi_{j} \rangle$.

matrix(o::RestEnergy, method::FiniteDifferenceMethod)

The matrix for the rest energy $mc^2$ is a diagonal matrix,

\[mc^2 \left(\begin{array}{ccccccc} 1 & 0 & 0 & \ldots \\ 0 & 1 & 0 & \ldots \\ 0 & 0 & 1 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right).\]

matrix(o::PotentialTerm, method::FiniteDifferenceMethod)

The matrix for the potential energy $V(r)$ is a diagonal matrix,

\[\pmb{V} = \left(\begin{array}{ccccccc} V(r_1) & 0 & 0 & \ldots \\ 0 & V(r_2) & 0 & \ldots \\ 0 & 0 & V(r_3) & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right).\]

function optimize(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...)

This function minimizes the energy by changing the exponents of the basis functions using Optim.jl.

\[\frac{\partial E}{\partial a_i} = 0\]

optimize(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4, optimizer=Optim.NelderMead())

This a optimizer for 1-basis calculations. This function returns optimize(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info, progress=progress, optimizer=optimizer, options...).

optimize(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4, optimizer=Optim.NelderMead())

This function minimizes the energy by optimizing $r_1$ and $r_n$ using Optim.jl.

\[\frac{\partial E}{\partial r_1} = \frac{\partial E}{\partial r_n} = 0\]

Return the number of binary sites in a QTT.

Contract a QTT at one vector or matrix index without materializing its entries.

Return the left boundary, internal, and right boundary ranks of a QTT.

solve(hamiltonian::Hamiltonian, wavefunction::Function, method::FiniteDifferenceMethod, info=4, nₘₐₓ=4)

solve(hamiltonian, model, method::VariationalNeuralNetwork;
      rng=Random.MersenneTwister(123), parameters=nothing, states=nothing,
      optimizer_state=nothing, trial=(r, value) -> value, info=0)

Minimize the finite-difference Rayleigh quotient with respect to the parameters of a Lux model. The model receives the complete radial grid as a batch and must return one real value per grid point. Pass returned parameters, states, and optionally optimizer_state to continue training.

solve(
    hamiltonian::Hamiltonian,
    candidates::BasisSet,
    method::BayesianVariationalMethod;
    rng=Random.MersenneTwister(123),
    perturbation=Hamiltonian(),
    info=4,
)

Select groups of basis functions from the finite candidate BasisSet with a Tanimoto-kernel Gaussian-process surrogate, then solve the accepted basis with the Rayleigh-Ritz method. Candidate energy evaluations are deterministic; the Bayesian uncertainty describes the discrete search surrogate, not uncertainty in the returned physical energy.

The result is a ResultRayleighRitz with additional method, candidate_basisset, selected_indices, history, n_evaluations, and n_rejected properties. Supply an explicit RNG for reproducibility.

solve(hamiltonian::Hamiltonian, basisset::BasisSet)

This function returns the eigenvalues $E$ and eigenvectors $\pmb{c}$ for

\[\pmb{H} \pmb{c} = E \pmb{S} \pmb{c}.\]

The Hamiltonian matrix is defined as $H_{ij} = \langle \phi_{i} | \hat{H} | \phi_{j} \rangle$. The overlap matrix is defined as $S_{ij} = \langle \phi_{i} | \phi_{j} \rangle$.

solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4)

This a solver for 1-basis calculations. This function returns solve(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info).

solve(hamiltonian::Hamiltonian, method::FiniteDifferenceMethod; perturbation=Hamiltonian(), info=4, nₘₐₓ=4)

This method solve the eigenvalue problem for the Hamiltonian discretized as a sparse matrix with finite difference approximation,

\[\pmb{H} \pmb{\psi} = E \pmb{\psi}.\]

The eigenvalue $E$ is an approximation of the exact energy and the eigenvector $\pmb{\psi}$ is a vector of the approximated values of the exact wavefunction $\psi(r)$ on points of the grid,

\[\pmb{\psi} = \left(\begin{array}{c} \psi(r_1) \\ \psi(r_2) \\ \psi(r_3) \\ \vdots \\ \end{array}\right).\]

solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.MersenneTwister(123), info=0)

Estimate the variational energy by averaging the local energy over Metropolis samples from $|\psi|^2$. A seeded random-number generator can be supplied with rng; by default, a new MersenneTwister(123) is used for a reproducible calculation. Set info to a positive value to print a result summary.

The returned named tuple echoes the input hamiltonian and method and contains E, variance, the naive standard_error, acceptance_rate, the sampling counts n_accepted, n_attempted, n_burn_in_discarded, n_samples, and n_discarded, as well as local_energies and samples. Data from multiple walkers are stored consecutively, with positions in the columns of samples. Non-finite local energies, which can occur at a measure-zero singularity such as the origin of a Coulomb potential, are excluded from the energy statistics and counted in n_discarded.

solve(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4)

This function is a wrapper for solve(hamiltonian::Hamiltonian, basisset::BasisSet, ...).

solve(hamiltonian, method::QuanticsTensorTrainMethod;
      initial=r -> exp(-r), info=4, nₘₐₓ=4)
solve(hamiltonian, initial::Function, method::QuanticsTensorTrainMethod;
      info=4, nₘₐₓ=4)

Compute the lowest nₘₐₓ QTT eigenstates. initial seeds the first DMRG solve; info controls progress output. The result contains energies E, unit-norm tensor trains C, radial functions ψ and u, and the diagnostics overlaps, residuals, histories, and rank_histories.

solve(hamiltonian, method::VariationalNeuralNetwork; kwargs...)

Build the standard Lux model specified by method.architecture and minimize its finite-difference Rayleigh quotient. See the three-argument overload to supply a custom Lux model.

Return the normalized radial part of a real complex-range Gaussian primitive.

Return a normalized Gaussian GEM primitive in position space.

Return the normalized radial part of a Gaussian GEM primitive.

φp(b::GaussianBasis, p) returns the radial momentum-space primitive. With the unitary Fourier convention

\[\widetilde\phi(\boldsymbol p)=(2\pi)^{-3/2}\int e^{i\boldsymbol p\cdot\boldsymbol r}\phi(\boldsymbol r)\,d^3r,\]

the normalized position-space primitive $N_l r^l e^{-a r^2}Y_l^m(\hat{\boldsymbol r})$ becomes

\[\widetilde\phi_{alm}(\boldsymbol p)= i^l N_l\frac{p^l}{(2a)^{l+3/2}}e^{-p^2/(4a)} Y_l^m(\hat{\boldsymbol p}).\]

The two-argument method returns only its radial factor; the four-argument method φp(b, p, θ, ϕ) includes the spherical harmonic. Position and momentum must use reciprocal units (for example, GeV⁻¹ and GeV when ℏ=c=1).

Return the full momentum-space Gaussian primitive.

Return the radial momentum-space Gaussian primitive (unitary Fourier convention).

ψp(result, p; n=1) evaluates a Gaussian-expansion eigenfunction in momentum space by summing result.C[i,n] * φp(result.basisset[i], p). It therefore uses the same unitary Fourier convention as φp. The four-argument form ψp(result, p, θ, ϕ; n=1) also includes the spherical harmonic.

This method currently applies to results built from real-range GaussianBasis primitives. Complex-range cos/sin primitives are intended for the position-space Rayleigh–Ritz calculation.

Evaluate a Rayleigh–Ritz eigenfunction in momentum space.