Add tests

[ohno]

The normalization and the orthogonality of eigenfunctions have already been tested. We need additional tests for the parameter dependence, the energy, the matching of Hamiltonian and eigenfunctions, etc..

• [x] other unit systems (:HydrogenAtom)
• [x] other unit systems (:CoulombTwoBody)
• [x] other unit systems (:PoschlTeller)
• [x] calculating energy by numerical integration (:DeltaPotential)
• [x] calculating energy by numerical integration (:PoschlTeller)
• [x] calculating energy by numerical integration (RigidRotor)
• [x] calculating energy by numerical integration (InfinitePotentialWell3D)
• [x] virial theorem (:SphericalOscillator)
• [x] virial theorem (:CoulombTwoBody)
• [x] recurrence relation of energies (the number of bound states v.s. maximum quantum number) (:MorsePotential)

[lhapp27]

It is not clear to me what is meant with "other unit systems".

[ohno]

Thank you for your efforts.

I meant "other parameters". I want to make sure that it is correct not only in atomic unit system but also in other systems of units (e.g. IS unit system).

atomic unit system:

HA = HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)

IS unit system:

a₀ = 5.29177210903e-11   # m  # https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
Eₕ = 4.3597447222071e-18 # J  # https://physics.nist.gov/cgi-bin/cuu/Value?hr
ℏ  = 1.054571817e-34     # J s # https://physics.nist.gov/cgi-bin/cuu/Value?hbar
mₑ = 9.1093837015e-31  # kg # https://physics.nist.gov/cgi-bin/cuu/Value?me
HA = HydrogenAtom(Z=1, Eₕ=Eₕ, a₀=a₀, mₑ=mₑ, ℏ=ℏ)

[lhapp27]

Thank you, but it is still not 100% clear what you intend to add for it. Let's take your example of the Hydrogen Atom:

HA = HydrogenAtom(Eₕ = eh, Z = z)

Then as of the code

function E(model::HydrogenAtom; n=1)
  Z = model.Z
  Eₕ = model.Eₕ
  return -Z^2/(2*n^2) * Eₕ
end

E(HA, n=n0) will provide the result of -z^2/(2*n0^2) * eh for any unit system of the user. E.g. by using eh=7 you will get -7/(2*n0^2).

Therefore the users can use any unit system they like.

[ohno]

Please see this commit https://github.com/ohno/Antique.jl/commit/4c96f2453eea62e050b8b579285ba9709a335c48 . The parameter dependency of the potential energy was wrong. I means that we want to find such mistakes.

julia> using Antique
julia> H = HydrogenAtom(Eₕ=1.0)
julia> V(H,1.0)
-1.0

julia> H = HydrogenAtom(Eₕ=2.0)
julia> V(H,1.0)
-1.0 # wrong!
-2.0 # correct!

In general, It is difficult to test energy using only energy. We should calculate the expected value of Hamiltonian to compare with energy. Based on the Virial theorem, I use the expected value of potential energy instead of the Hamiltonian:

using Test
using Printf
using QuadGK
using Antique

# https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
# https://physics.nist.gov/cgi-bin/cuu/Value?hr
# https://physics.nist.gov/cgi-bin/cuu/Value?hbar
# https://physics.nist.gov/cgi-bin/cuu/Value?me
# https://physics.nist.gov/cgi-bin/cuu/Value?hrev

for HA in [
  HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0),
  HydrogenAtom(Z=1, Eₕ=1.0, a₀=2.0, mₑ=1.0, ℏ=1.0),
  HydrogenAtom(Z=2, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0),
  HydrogenAtom(Z=2, Eₕ=27.211386245988, a₀=1.0, mₑ=1.0, ℏ=1.0),
  HydrogenAtom(Z=2, Eₕ=4.3597447222071e-18, a₀=5.29177210903e-11, mₑ=9.1093837015e-31, ℏ=1.054571817e-34),
]
  println(HA)
  @testset "<ψₙ|V|ψₙ> / 2 = Eₙ" begin
    println(" n |        analytical |         numerical ")
    println("-- | ----------------- | ----------------- ")
    for n in 1:10
      analytical = E(HA, n=n) * 2
      numerical  = quadgk(r -> 4*π*r^2 * conj(ψ(HA,r,0,0, n=n)) * V(HA,r) * ψ(HA,r,0,0, n=n), 0, HA.a₀*50*n, maxevals=10^3)[1]
      acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5)
      @test acceptance
      @printf("%2d | %17.12e | %17.12e %s\n", n, analytical, numerical, acceptance ? "✔" : "✗")
    end
  end
end

[ohno]

@lhapp27 Please check the dependence not only on λ but also on m, ℏ, and x₀ in :PoschlTeller potential.

[ajarifi]

I will try to check the energy for Deltapotential and Rigidrotor.

[ajarifi]

For the delta potential, when I compute the eigenenergy numerically, there is an error or inaccuracy in computing the derivative. This is because there is a cusp at the origin. I can compare the eigenenergy and the expectation value of the Hamiltonian analytically, not numerically.

[ajarifi]

Perhaps it is okay like this.

@testset "<ψₙ|H|ψₙ> = ∫ψₙ*Hψₙdx = Eₙ" begin
function ψHψ(DP)
    κ = DP.m * DP.α/ DP.ℏ^2
    T = -DP.ℏ^2 / (2 * DP.m) *( κ^2 - 2κ * ψ(DP,0)^2 )
    V = -DP.α * ψ(DP, 0)^2
    return T + V
end

println("  α |   m |   ℏ |        analytical |         numerical ")
println("--- | --- | --- | ----------------- | ----------------- ")
for α in [0.1, 0.5, 2.0]
for m in [0.1, 0.5, 2.0]
for ℏ in [0.1, 0.5, 2.0]
    DP = DeltaPotential(α=α, m=m, ℏ=ℏ)
    analytical = E(DP)
    numerical = ψHψ(DP)
    acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-4) : isapprox(analytical, numerical, rtol=1e-4)
    @printf("%.1f | %.1f | %.1f | %17.12f | %17.12f %s\n", α, m, ℏ,analytical, numerical, acceptance ? "✔" : "✗") 
end
end
end
end

[ajarifi]

Also for Rigidrotor, I computed the expectation value analytically. I hope this is okay.

@testset "<ψₙ|H|ψₙ> = ∫ψₙ*Hψₙdx = Eₙ" begin
    function ψHψ(RR;l)
        μ = (1/RR.m₁ + 1/RR.m₂)^(-1)
        I = μ * RR.R^2
        YY = real(quadgk(φ -> quadgk(θ -> conj(Y(RR,θ,φ,l=l,m=0)) * Y(RR,θ,φ,l=l,m=0) * sin(θ), 0, π, maxevals=10)[1], 0, 2π, maxevals=10)[1])
        T = RR.ℏ^2/2/I * l*(l+1) * YY
        return T 
    end

    ℏ = 1
    println(" m₁ |  m₂ |  R  |  l  |       analytical  |         numerical ")
    println("--- | --- | --- | --- | ----------------- | ----------------- ")
    for m₁ in [0.5, 2.0]
    for m₂ in [0.5, 2.0]
    for R in [0.5, 2.0]
    for l in [1, 2, 3]
        RR = RigidRotor(m₁=m₁, m₂=m₂, R=R, ℏ=ℏ)
        analytical = E(RR,l=l)
        numerical = ψHψ(RR,l=l)
        acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-4) : isapprox(analytical, numerical, rtol=1e-4)
        @printf("%.1f | %.1f | %.1f | %.1f | %17.12f | %17.12f %s\n", m₁, m₂, R, l,analytical, numerical, acceptance ? "✔" : "✗") 
    end
    end
    end
    end
end