Rayleigh-Ritz Method

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, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4)

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

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\]

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

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)

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

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}}\]

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.

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

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::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::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(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}\]

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::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::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::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::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::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::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.\]