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:

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

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, ...], [basis1, basis2, ...])

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

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