iTranslated by AI

The content below is an AI-generated translation. This is an experimental feature, and may contain errors. View original article

💯

Introducing Antique.jl

2023/12/11に公開

Julia Advent Calendar 2023 Day 11🎄

Introduction

I have developed and would like to introduce a package called "Antique.jl" which is useful for students, teachers, and researchers involved in quantum mechanics. It can be utilized for presentation materials for students, teaching materials for quantum mechanics for teachers, and benchmarks for calculation methods for researchers. Please read on for specific examples of how to use it.

Antique.jl

https://ohno.github.io/Antique.jl/

This package currently supports four models: Infinite Potential Well, Harmonic Oscillator, Morse Potential, and Hydrogen Atom. It provides modules that summarize analytical solutions (eigenvalues, eigenfunctions) and potentials for the Schrödinger equation in these models. For all models, tests are conducted using numerical integration (QuadGK.jl) or a symbolic computation system (Symbolics.jl). We plan to sequentially add models from this list. If you have any requests, please post them in the Issues, and I will prioritize them.

Why is it necessary?

In short, "because implementing analytical solutions is just a pain." It would be putting the cart before the horse if you spent all your time on benchmark implementation and testing when you are trying to create new methods. My intention is to save that trouble by wrapping them into a package.

Road to Release

The story goes back to January last year (2022). I published three articles titled Visualization with Julia Series.

As mentioned earlier, the purpose of writing these articles was to provide benchmarks for calculation methods, but I also used the figures to explain quantum mechanical concepts (such as zero-point vibration, the Born-Oppenheimer approximation, and vibrational excited states) in conference and master's thesis presentations. Antique.jl provides improved versions of the analytical solutions (eigenvalues and eigenfunctions) and potentials defined in those articles, organized into modules.

If you read the articles, you'll notice that I particularly struggled with handling the Associated Laguerre polynomials and Associated Legendre polynomials included in the eigenfunctions of the hydrogen atom. These polynomials are most often defined using Rodrigues' formula, but since that formula involves differential operations and is difficult to implement directly, implementations typically use closed forms or recurrence relations. However, since most quantum mechanics textbooks do not list the closed forms, implementation is challenging. While there are already formula collections and packages for orthogonal polynomials, the verification process (testing) took an enormous amount of time because it was often unclear which definition (Rodrigues' formula) the implementation was based on, or the definition differed slightly from the textbook (since they are solutions to differential equations, some variation in normalization or sign is allowed), making it impossible to use the closed forms as-is.

To solve this problem, in last year's (2022) Advent Calendar, I wrote an article on how to automatically test whether the expansion of Rodrigues' formula matches the closed form using Symbolics.jl. I am actually utilizing the testing techniques described in that article.

Additionally, in the "Visualization with Julia" series, I verified the orthogonality of eigenfunctions, eigenvalues, and the virial theorem using numerical integration. Antique.jl employs a testing strategy that combines both numerical integration and symbolic computation. If you are interested, please see the Introduction to Testing Techniques. Next, I will introduce two use cases.

Installation

Simply run the following in your REPL or Jupyter Notebook. The usage will be explained through the use cases.

using Pkg; Pkg.add("Antique")

Use Case: Quantum Mechanics Teaching Materials

While the original purpose was to provide benchmarks for method development, it is also extremely useful as a teaching material for quantum mechanics. First, let's obtain the module for the hydrogen atom.

using Antique
H = HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0)

Now, let's calculate the wavelengths of the Lyman series.

# Lyman series n -> 1
wavelength(ΔE) = 1 / (ΔE * 2.1947463136320e-2) # nm, https://physics.nist.gov/cgi-bin/cuu/CCValue?hrminv
for i in 2:7
    println("$i -> 1\t", wavelength( E(H,n=i) - E(H,n=1) ), " nm")
end

# output
2 -> 1	121.50227341098497 nm
3 -> 1	102.51754319051858 nm
4 -> 1	97.20181872878798 nm
5 -> 1	94.92365110233202 nm
6 -> 1	93.73032520275984 nm
7 -> 1	93.02517808028537 nm

All wavelengths match the values on Wikipedia. Next, let's plot the radial distribution functions.

using Plots
plot(xlabel="\$r~/~a_0\$", ylabel="\$r^2|R_{nl}(r)|^2~/~a_0^{-1}
$", ylims=(-0.01,0.55), xticks=0:1:20, size=(480,400), dpi=300)
for n in 1:3
  for l in 0:n-1
    plot!(0:0.01:20, r->r^2*Antique.R(H,r,n=n,l=l)^2, lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1], label="\$n = $n, l=$l\$")
  end
end
plot!()

This matches well with this graph. Note that these are not real-valued orbitals like $2\mathrm{p_x}$ or $2\mathrm{p_y}$ . You would need to take linear combinations of the above eigenfunctions to eliminate the imaginary parts. Finally, let's calculate the wavelength of the hyperfine splitting (21 cm line).

# constants: https://doi.org/10.1103/RevModPhys.93.025010
e  = 1.602176634e-19    # C      https://physics.nist.gov/cgi-bin/cuu/Value?e
h  = 6.62607015e-34     # J Hz-1 https://physics.nist.gov/cgi-bin/cuu/Value?h
c  = 299792458          # m s-1  https://physics.nist.gov/cgi-bin/cuu/Value?c
a0 = 5.29177210903e-11  # m      https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
µ0 = 1.25663706212e-6   # N A-2  https://physics.nist.gov/cgi-bin/cuu/Value?mu0
µB = 9.2740100783e-24   # J T-1  https://physics.nist.gov/cgi-bin/cuu/Value?mub
µN = 5.0507837461e-27   # J T-1  https://physics.nist.gov/cgi-bin/cuu/Value?mun
ge = 2.00231930436256   #        https://physics.nist.gov/cgi-bin/cuu/Value?gem
gp = 5.5856946893       #        https://physics.nist.gov/cgi-bin/cuu/Value?gp

# calculation: https://doi.org/10.1119/1.12733
δ = abs(ψ(H,0,0,0))^2
ΔE = 2 / 3 * µ0 * µN * µB * gp * ge * δ * a0^(-3)
println("1/π    = ", 1/π)
println("<δ(r)> = ", δ, " a₀⁻³")
println("<δ(r)> = ", δ * a0^(-3), " m⁻³")
println("ΔE = ", ΔE, " J")
println("ν = ΔE/h = ", ΔE / h * 1e-6, " MHz")
println("λ = hc/ΔE = ", h*c/ΔE*100, " cm")

# output
1/π    = 0.3183098861837907
<δ(r)> = 0.3183098861837908 a₀⁻³
<δ(r)> = 2.1480615849063944e30 m⁻³
ΔE = 9.427622831641132e-25 J
ν = ΔE/h = 1422.8075794882932 MHz
λ = hc/ΔE = 21.070485027063118 cm

Here, following Griffiths (1982), we calculated the energy difference, frequency, and wavelength of the hyperfine splitting perturbatively using the value of the ground-state wave function (squared absolute value) at $r=0$ . As known by the name 21 cm line, the calculated wavelength is approximately 21 cm. The 1/π shown in the output corresponds to the value of $\langle \delta^3(\pmb{r}) \rangle_\mathrm{1S} = |\psi_\mathrm{1S}(\pmb{0})|^2 = \frac{1}{\pi} a_0^{-3}$ excluding the units. I like how a clean number like $1/\pi$ appears for the probability density at the origin of the hydrogen atom (the probability that the proton and electron exist at the same position).

In addition to these examples, tests for the orthogonality of eigenfunctions, the virial theorem, and calculations of expectation values for various physical quantities are implemented and running. All of these would serve as good topics for exercises or assignments.

Use Case: Benchmarks for Calculation Methods

QuantumOptics.jl can calculate the eigenvalues and eigenfunctions of a harmonic oscillator. Let's test whether the calculated eigenvalues match the values from Antique.jl.

# Parameters
k = 1.0
m = 1.0
ℏ = 1.0

# https://ohno.github.io/Antique.jl/dev/HarmonicOscillator/#Usage-and-Examples
# using Pkg; Pkg.add("Antique")
using Antique
HO = HarmonicOscillator(k=k, m=m, ℏ=ℏ)

# https://qojulia.org/
# using Pkg; Pkg.add("QuantumOptics")
using QuantumOptics
basis = PositionBasis(-5, 5, 100)
x = position(basis)
p = momentum(basis)
H = p^2*ℏ^2/2/m + k/2*DenseOperator(x^2)
energies, states = eigenstates((H+dagger(H))/2, 10)

# Antique.jl vs QuantumOptics.jl
using Printf
println(" n\tAntique.jl(exact)  \tQuantumOptics.jl")
for i in keys(energies)
    @printf("%2d\t%.15f\t%.15f\n", i-1, E(HO,n=i-1), energies[i])
end

 n	Antique.jl(exact)  	QuantumOptics.jl
 0	0.500000000000000	0.499999999920220
 1	1.500000000000000	1.500000003671307
 2	2.500000000000000	2.499999911435342
 3	3.500000000000000	3.500001221372221
 4	4.500000000000000	4.499986371258386
 5	5.500000000000000	5.500098711913332
 6	6.500000000000000	6.499320233969457
 7	7.500000000000000	7.502927850989741
 8	8.500000000000000	8.485825750231648
 9	9.500000000000000	9.536571503980930

We can see that QuantumOptics.jl is able to calculate the harmonic oscillator with fairly high precision. Here, the spatial basis range is set from $-5$ to $5$ , but if you narrow this, the accuracy decreases significantly. Also, not just the basis range, but since higher excited states have more nodes, the results will suffer without sufficient grid resolution. Of course, the accuracy also changes depending on the parameters. Please give it a try.

Let's also plot the probability density function just in case. It seems that states[n].data is normalized by $\sum |\psi(x)|^2=1$ rather than $\sum |\psi(x)|^2 \Delta x=1$ , so we need to divide it by $\Delta x$ .

using Plots
default(lw=2)
X = samplepoints(basis)
Δx = (maximum(X) - minimum(X)) / length(X)
for i in 0:9
    plot(xlabel="x", ylabel="ψ(x)", title="v = $i")
    plot!(-5:0.01:5, x -> abs(ψ(HO,x,n=i))^2, label="Antique.jl(exact)")
    scatter!(X, abs2.(states[i+1].data) / Δx, label="QuantumOptics.jl", ms=3, msw=0)
    plot!() |> display
end

The results match very well. We can also see how the description worsens in higher excited states.

Conclusion

If you are developing a new method, please use this as a benchmark. For requests for new models, etc., please head here. I look forward to your stars 🌟

https://github.com/ohno/Antique.jl

https://ohno.github.io/Antique.jl/dev/

For Developers

Introduction to Documentation

 Introduction to DocumentationThe documentation for each model summarizes the following three points:
A self-contained explanation of the model
Usage and concrete examples
Tests using symbolic computation and numerical integration
For example, the following is an example of the harmonic oscillator.
There are explanations for the Hamiltonian, eigenvalues, and eigenfunctions, and it also introduces the definition and closed form of the Hermite polynomials included in the eigenfunctions, along with concrete examples.
Usage and concrete examples are also summarized in the documentation. antique(model_name, parameters) returns a module, and since each module contains the potential V(x), eigenvalues E(), and eigenfunctions ψ(x), you can call and use them.

It might seem a bit odd that test execution results are also included in the documentation, but I think the fact that you can check what kind of tests are being performed, including the mathematical formulas, is an initiative worth appreciating. I will introduce it in the next section.

Introduction to Testing Techniques

 Introduction to Testing TechniquesJulia's flexible type system makes it possible to run tests using both a symbolic computation system (Symbolics.jl) and numerical integration (QuadGK.jl) on the same piece of code. Here, I'll introduce an example of the harmonic oscillator. It verifies whether the Rodrigues' formula and the closed form of Hermite polynomials match using symbolic computation, and the output results are included in the documentation.
The orthogonality of Hermite polynomials is verified using numerical integration:

Design Philosophy (About Old Versions)

In this package, the following two give the same result.

HO = Antique.HarmonicOscillator
@show HO.E(n=0, k=2.0, m=3.0)
@show HO.ψ(0, n=0, k=2.0, m=3.0)

HO = antique(:HarmonicOscillator, k=2.0, m=3.0, ℏ=1.0)
@show HO.E(n=0)
@show HO.ψ(0, n=0)

In the former case, parameters must be specified every time eigenvalues and eigenfunctions are called, whereas in the latter case, they are only specified once at the beginning. I would prefer to use the latter if possible. The latter is designed so that executing antique(model_name, parameters) retrieves the module's source code, replaces the parameters, and returns a module re-evaluated through metaprogramming. I have a feeling something might be wrong with this design philosophy, but it improves both code beauty and usability. I'll explain step-by-step, so first, please take a look at the Harmonic Oscillator module.

module HarmonicOscillator

    # Default
    k = 1.0 # change here!
    m = 1.0 # change here!
    ℏ = 1.0 # change here!

    # Potential
    V(x; k=k, m=m, ω=sqrt(k/m)) = 1/2*k*x^2

    # Energy
    E(; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ) = ℏ*ω*(n+1/2)

    # Wave Function
    function ψ(x; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ)
        A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ)))
        ξ = sqrt(m*ω/ℏ) * x
        return A*H(ξ,n=n)*exp(-ξ^2/2)
    end

    # Hermite polynomials
    H(x; n=0) = factorial(n) * sum(m -> (-1)^m // (factorial(m)  * factorial(n-2*m)) * (2*x)^(n-2*m), 0:Int(floor(n/2)))

end

In this way, the parameters defining the system (e.g., k, m, ℏ for a harmonic oscillator) and its eigenvalues E() and eigenfunctions ψ() are gathered into a single module. In Python, you would typically define methods that refer to parameters within a class and then create instances with different parameters. However, since Julia doesn't have the concept of classes, I debated whether to make the module parameters mutable or have methods within a struct (though I suspect wanting to use classes might go against Julia's design philosophy in the first place). Since making structs hold methods seemed quite difficult, I decided to make the module parameters mutable. While it's not impossible to change parameters, you cannot maintain two systems with different parameters simultaneously. Specifically, as shown in the Examples in the documentation:

julia> H = antique(:HydrogenAtom, Z=1)
julia> H.E(n=1)
-0.5

julia> He⁺ = antique(:HydrogenAtom, Z=2)
julia> He⁺.E(n=1)
-2.0

As used here, it is impossible to have two modules with different parameters, like a hydrogen atom and a helium ion, coexist. Therefore, I took the approach of replacing the change here! parts in the source code and re-evaluating them as a different module using eval(). However, a problem with this approach is that the module is created in Base.Meta instead of within Antique, causing dependency packages (the Morse potential depends on the gamma function in SpecialFunctions.jl) not to be loaded. If SpecialFunctions.jl is not installed,

MO = antique(:MorsePotential, μ=1.0, k=1.0, ℏ=1.0)

will result in an error, so you must use

MO = Antique.MorsePotential
MO.μ = 1.0
MO.k = 1.0
MO.ℏ = 1.0

instead. This dependency loading issue is identified as a future task. If you have any suggestions for a better design, please let me know.

Postscript

~~v0.1.1 has been released!~~

v0.5.1 has been released!

Acknowledgments

This research is a result under the RIKEN Junior Research Associate (JRA) program.