API reference

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:

\[\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),
)

Constant(constant=1)

\[+ c\]

ContractedBasis([c1, c2, ...], [basis1, basis2, ...])

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

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_i(r, θ, φ) = N _{il} r^l \exp(-a_i r^2) Y_l^m(θ, φ)\]

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

This is a basis set with exponentials generated by geometric(r₁, rₙ, n; nₘₐₓ=n, nₘᵢₙ=1). You can define the same basis set as Table A2 in 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(SimpleGaussianBasis, 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\]

PrimitiveBasis <: Basis is an abstract type.

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.

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.

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::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}}{\mathrm{d}r^{2}}(r) = \frac{\psi(r+\Delta r) - 2f(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}{\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\]

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

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