GitHub

TwoBody.jl

Stable Dev Build Status

TwoBody.jl: a Julia package for two-body problems

Install

Run the following code on the REPL or Jupyter Notebook to install this package.

import Pkg; Pkg.add(url="https://github.com/ohno/TwoBody.jl.git")

Usage

Run the following code before each use.

using TwoBody

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

H = Hamiltonian(
  NonRelativisticKinetic(ℏ = 1 , m = 1),
  CoulombPotential(coefficient = -1),
)

The usage depends on the method. Define the basis set for the Rayleigh–Ritz Method:

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

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(H, BS)
# method

Rayleigh–Ritz method with SimpleGaussianBasis
J. Thijssen, Computational Physics 2nd Edition (2013)
https://doi.org/10.1017/CBO9781139171397

# basis function

φ₁(r) = TwoBody.φ(SimpleGaussianBasis(a=13.00773), r)
φ₂(r) = TwoBody.φ(SimpleGaussianBasis(a=1.962079), r)
φ₃(r) = TwoBody.φ(SimpleGaussianBasis(a=0.444529), r)
φ₄(r) = TwoBody.φ(SimpleGaussianBasis(a=0.1219492), r)

# eigenfunction

ψ₁(r) = + 0.096102φ₁(r) + 0.163017φ₂(r) + 0.185587φ₃(r) + 0.073701φ₄(r) 
ψ₂(r) = + 0.119454φ₁(r) + 0.081329φ₂(r) + 0.496216φ₃(r) - 0.205916φ₄(r) 
ψ₃(r) = - 0.010362φ₁(r) + 1.744891φ₂(r) - 0.629196φ₃(r) + 0.097774φ₄(r) 
ψ₄(r) = - 6.155100φ₁(r) + 1.240202φ₂(r) - 0.226412φ₃(r) + 0.030780φ₄(r) 

# eigenvalue

E₁ = -0.4992784056674876
E₂ = 0.11321392045798988
E₃ = 2.592299571959808
E₄ = 21.144365190122507

# others

n 	norm, <ψₙ|ψₙ> = cₙ' * S * cₙ
1	1.0
2	1.0000000000000004
3	1.0
4	0.9999999999999988

n 	error check, |<ψₙ|H|ψₙ> - E| = |cₙ' * H * cₙ - E| = 0
1	1.8318679906315083e-15
2	3.4833247397614286e-15
3	3.552713678800501e-15
4	1.7763568394002505e-14

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)), perturbation = Hamiltonian(), basisset = BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), nₘₐₓ = 4, H = [0.5772684658780091 0.072002466903411 -0.32154049871511875 -0.43612626272313476; 0.072002466903411 0.5070499286923358 -0.9891848263978418 -2.377419924763058; -0.32154049871511875 -0.9891848263978418 -2.6380824346106566 -7.342216931051755; -0.43612626272313476 -2.377419924763058 -7.342216931051755 -17.30516271277891], 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], 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], φ = TwoBody.var"#42#49"{Int64, BasisSet}[TwoBody.var"#42#49"{Int64, BasisSet}(1, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492))), TwoBody.var"#42#49"{Int64, BasisSet}(2, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492))), TwoBody.var"#42#49"{Int64, BasisSet}(3, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492))), TwoBody.var"#42#49"{Int64, BasisSet}(4, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)))], ψ = TwoBody.var"#44#51"{Int64, BasisSet, Matrix{Float64}, Int64}[TwoBody.var"#44#51"{Int64, BasisSet, Matrix{Float64}, Int64}(1, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), [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], 4), TwoBody.var"#44#51"{Int64, BasisSet, Matrix{Float64}, Int64}(2, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), [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], 4), TwoBody.var"#44#51"{Int64, BasisSet, Matrix{Float64}, Int64}(3, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), [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], 4), TwoBody.var"#44#51"{Int64, BasisSet, Matrix{Float64}, Int64}(4, BasisSet(SimpleGaussianBasis(a=13.00773), SimpleGaussianBasis(a=1.962079), SimpleGaussianBasis(a=0.444529), SimpleGaussianBasis(a=0.1219492)), [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], 4)])

The wave function is also good. However, the Gaussian basis does not satisfy the Kato’s cusp condition.

# solve
res = solve(H, BS, info=0)

# benchmark
import Antique
HA = Antique.HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)

# plot
using CairoMakie
fig = Figure(size=(420,300), fontsize=11.5, backgroundcolor=:transparent)
axis = Axis(fig[1,1], xlabel=L"$r / a_0$", ylabel=L"$\psi(r) / a_0^{-3/2}$", ylabelsize=16.5, xlabelsize=16.5, limits=(0,4,0,1.1/sqrt(π)))
lines!(axis, 0..5, r -> abs(res.ψ[1](r)), label="TwoBody.jl")
lines!(axis, 0..5, r -> abs(Antique.ψ(HA,r,0,0)), linestyle=:dash, color=:black, label="Antique.jl")
axislegend(axis, position=:rt, framevisible=false)
fig
Example block output

API reference