Rayleigh–Ritz method
The solver finds the best $c_i$ for the trial wavefunction
\[\varPsi(r) = \sum_i c_i \phi_i(r),\]
to minmize the expectation value of the energy
\[E = \frac{\langle\psi|\hat{H}|\psi\rangle}{\langle\psi|\psi\rangle}.\]
Here, the energy by trial wavefunction is the upper bound for the exact energy.
Exaples
Define the Hamiltoninan. 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}\]
hamiltonian = Hamiltonian(
NonRelativisticKinetic(ℏ = 1 , m = 1),
CoulombPotential(coefficient = -1),
)
Define the basis set:
\[\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}\]
basisset = BasisSet(
SimpleGaussianBasis(13.00773),
SimpleGaussianBasis(1.962079),
SimpleGaussianBasis(0.444529),
SimpleGaussianBasis(0.1219492),
)
You should find
\[E_{n=1} = -0.499278~E_\mathrm{h},\]
which is amazingly good for only four basis functions according to Thijssen(2007). The exact ground-state energy is $-0.5~E_\mathrm{h}$.
julia> solve(hamiltonian, basisset)
n basis function φₙ 1 SimpleGaussianBasis(a=13.00773) 2 SimpleGaussianBasis(a=1.962079) 3 SimpleGaussianBasis(a=0.444529) 4 SimpleGaussianBasis(a=0.1219492) n wavefuntion ψₙ 1 + 0.096102 φ₁ + 0.163017 φ₂ + 0.185587 φ₃ + 0.073701 φ₄ 2 + 0.119454 φ₁ + 0.081329 φ₂ + 0.496216 φ₃ - 0.205916 φ₄ 3 - 0.010362 φ₁ + 1.744891 φ₂ - 0.629196 φ₃ + 0.097774 φ₄ 4 - 6.155100 φ₁ + 1.240202 φ₂ - 0.226412 φ₃ + 0.030780 φ₄ n norm, <ψ|ψ> = c' * S * c 1 0.9999999999999998 2 1.0000000000000004 3 1.0000000000000007 4 0.9999999999999988 n eigenvalue, E 1 -0.4992784056674876 2 0.11321392045798988 3 2.592299571959808 4 21.144365190122507 n expectation value of the Hamiltonian, <ψ|H|ψ> = c' * H * c 1 -0.49927840566748577 2 0.11321392045798652 3 2.592299571959811 4 21.14436519012249 n expectation value of NonRelativisticKinetic(ħ=1, m=1) 1 0.4992783686700055 2 0.8428088332141157 3 4.432656608731447 4 26.465623640332108 n expectation value of CoulombPotential(coefficient=-1) 1 -0.9985567743374912 2 -0.7295949127561296 3 -1.8403570367716342 4 -5.321258450209621 (hamiltonian = Hamiltonian(NonRelativisticKinetic(ħ=1, m=1), CoulombPotential(coefficient=-1)), basisset = BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), E = [-0.4992784056674876, 0.11321392045798988, 2.592299571959808, 21.144365190122507], C = [0.09610151618612488 0.1194538057449333 -0.010362061881687392 -6.155100006789123; 0.16301716963905885 0.08132945379475047 1.7448913023470436 1.2402020851506472; 0.18558698714513683 0.49621626366832666 -0.6291955141735303 -0.22641160819529882; 0.07370076069275631 -0.20591550816511817 0.09777447415099819 0.030779842546714373], S = [0.041964064408426524 0.0961391814715395 0.11285790355607861 0.11704251263182287; 0.0961391814715395 0.7163167080668228 1.4914777365294443 1.8508423236885296; 0.11285790355607861 1.4914777365294443 6.642471010530628 13.060205391889545; 0.11704251263182287 1.8508423236885296 13.060205391889545 46.22866820431064], H = [0.5772684658780091 0.072002466903411 -0.32154049871511875 -0.43612626272313476; 0.072002466903411 0.5070499286923358 -0.9891848263978418 -2.377419924763058; -0.3215404987151187 -0.9891848263978418 -2.6380824346106566 -7.342216931051755; -0.43612626272313476 -2.3774199247630583 -7.342216931051756 -17.30516271277891])
Solver
TwoBody.solve
— Functionsolve(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, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4)
This function is a wrapper for solve(hamiltonian::Hamiltonian, basisset::BasisSet, ...)
.
TwoBody.optimize
— Functionoptimize(hamiltonian::Hamiltonian, basisset::BasisSet; 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\]
Basis Set
TwoBody.BasisSet
— TypeBasisSet(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}\]
basisset = BasisSet(
SimpleGaussianBasis(13.00773),
SimpleGaussianBasis(1.962079),
SimpleGaussianBasis(0.444529),
SimpleGaussianBasis(0.1219492),
)
TwoBody.GeometricBasisSet
— TypeGeometricBasisSet(basistype, r₁, rₙ, n; nₘₐₓ=n, nₘᵢₙ=1)
This is a basis set with exponentials generated by geometric()
.
TwoBody.geometric
— FunctionExponents 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
Basis Functions
TwoBody.SimpleGaussianBasis
— TypeSimpleGaussianBasis(a=1)
\[\phi_i(r) = \exp(-a_i r^2)\]
Note: This basis is not normalized. This is only for s-wave.
TwoBody.GaussianBasis
— TypeGaussianBasis(a=1, l=0, m=0)
\[\phi_i(r, θ, φ) = N _{il} r^l \exp(-a_i r^2) Y_l^m(θ, φ)\]
TwoBody.ContractedBasis
— TypeContractedBasis([c1, c2, ...], [basis1, basis2, ...])
\[\phi' = \sum_i c_i \phi_i\]
Matrix Elements
TwoBody.element
— Functionelement(SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
\[\begin{aligned} S_{ij} = \langle \phi_{i} | \phi_{j} \rangle &= \iiint \phi_{i}^*(r) \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}^{-(\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::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::NonRelativisticKinetic, SGB1::SimpleGaussianBasis, SGB2::SimpleGaussianBasis)
\[\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}\]
or
\[\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}\]
element(o::ConstantPotential, 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::LinearPotential, 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:
\[\begin{aligned} \int_0^{\infty} r^{2n+1} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{n!}{2 a^{n+1}} \end{aligned}\]
element(o::CoulombPotential, 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:
\[\begin{aligned} \int_0^{\infty} r^{2n+1} \exp \left(-a r^2\right) ~\mathrm{d}r = \frac{n!}{2 a^{n+1}} \end{aligned}\]
element(o::PowerLawPotential, 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::GaussianPotential, 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}\]