iTranslated by AI

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

☀️

Visualization with Julia: Hydrogen Atom

2022/01/09に公開

In this third installment of the "Visualization with Julia" series, we will visualize the solutions to the Schrödinger equation for the hydrogen atom.

Spectrum
Radial distribution

$l=0$	$l=1$	$l=2$

Packages

The following six packages need to be installed beforehand using the package mode. In the notebook, we declare using Plots, using Printf, using LaTeXStrings, using Polynomials, and using SpecialPolynomials. Note that since Plots.jl and GR.jl often interfere with each other, we will declare them later.

Packages

# using Pkg
# Pkg.add("GR")
# Pkg.add("Plots")
# Pkg.add("Printf")
# Pkg.add("LaTeXStrings")
# Pkg.add("Polynomials")
# Pkg.add("SpecialPolynomials")
# using GR
using Plots
using Printf
using LaTeXStrings
using Polynomials
using SpecialPolynomials

Schrödinger Equation

The time-independent Schrödinger equation for a hydrogen-like atom is

\left[ -\frac{\hbar^2}{2m_\mathrm{e}}\nabla_\mathrm{e}^2 -\frac{\hbar^2}{2M}\nabla_\mathrm{n}^2 - \frac{Ze^2}{4\pi\varepsilon_0 r} \right] \psi_\mathrm{tot}(\pmb{r}_\mathrm{e},\pmb{r}_\mathrm{n}) = E_\mathrm{tot} \psi_\mathrm{tot}(\pmb{r}_\mathrm{e},\pmb{r}_\mathrm{n})

where $m_\mathrm{e}$ and $M$ are the masses of the electron and nucleus, respectively; $\pmb{r}_\mathrm{e}$ and $\pmb{r}_\mathrm{n}$ are their positions; $\nabla_\mathrm{e}^2$ and $\nabla_\mathrm{n}^2$ are their respective Laplacians; and their charges are $-e$ and $Ze$ . By using the reduced mass $\frac{1}{\mu}:=\frac{1}{m_\mathrm{e}}+\frac{1}{M}$ and the relative motion $\pmb{r}:=\pmb{r}_\mathrm{e}-\pmb{r}_\mathrm{n}$ to separate the translational degrees of freedom, we obtain

\left[ -\frac{\hbar^2}{2\mu}\nabla^2 - \frac{Ze^2}{4\pi\varepsilon_0 r} \right] \psi(\pmb{r}) = E \psi(\pmb{r})

For details on this separation, please refer to Harada (2007), Takayanagi (2000), Kawai (2002), Jensen (2007), etc. Furthermore, approximating $M\rightarrow\infty$ leads to $\mu\rightarrow m_\mathrm{e}$ , resulting in

\left[ -\frac{\hbar^2}{2m_\mathrm{e}}\nabla^2 - \frac{Ze^2}{4\pi\varepsilon_0 r} \right] \psi(\pmb{r}) = E \psi(\pmb{r})

As is well known, the eigenvalues for this Schrödinger equation are

E_n = -\frac{m_\mathrm{e} e^4 Z^2}{2n^2(4\pi\varepsilon_0)^2\hbar^2} = -\frac{Z^2}{2n^2} E_\mathrm{h} \tag{2.29'}

(Modified slightly based on Takayanagi (2000)), where $E_\mathrm{h} = \frac{m_\mathrm{e} e^4}{(4\pi\varepsilon_0)^2 \hbar^2}$ is the Hartree energy, an atomic unit. For an explanation of atomic units, I recommend reading IUPAC GreenBook 3.9.2, as it is exceptionally clear and accurate. The eigenfunctions are

\psi_{nlm}(\pmb{r}) = R_{nl}(r) Y_{lm}(\theta,\varphi) \tag{2.10}

R_{nl}(r) = -\sqrt{\frac{(n-l-1)!}{2n[(n+l)!]^3} \left(\frac{2Z}{n a_0}\right)^3} \left(\frac{2Zr}{n a_0}\right)^l \exp \left(-\frac{Zr}{n a_0}\right) L_{n+l}^{2l+1} \left(\frac{2Zr}{n a_0}\right) \tag{2.25'}

Y_{lm}(\theta,\varphi) = i^{|m|+m} \sqrt{\frac{(2l+1)(l-|m|)!}{2(l+|m|)!}} P_l^{|m|} (\cos\theta) \frac{e^{im\varphi}}{\sqrt{2\pi}} \tag{2.15'}

(Modified slightly based on Takayanagi (2000)). Note that the above expressions are based on the definitions of Laguerre polynomials as $L_n(x)=\mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)$ and associated Laguerre polynomials as $L_n^{k}(x) = \frac{\mathrm{d}^k}{\mathrm{d}x^k} L_n(x)$ . Furthermore, while the definition of associated Legendre polynomials $P_k^m(t)$ is generally consistent, the phase of the spherical harmonics $Y_{lm}(\theta,\varphi)$ varies more significantly between different references than that of the radial functions, so extra care is needed.

Potential

V(r;Z=1) = -Z/abs(r)

Eigenvalues

E(n;Z=1) = -Z^2/(2*n^2)

Radial_function

function R(n,l,r;Z=1)
    a0 = 1
    ρ = 2*Z*abs(r)/(n*a0)
    N = -sqrt( factorial(n-l-1)/(2*n*factorial(n+l)^3) * (2*Z/(n*a0))^3 )
    return N*ρ^l * exp(-ρ/2) * L(n+l,2*l+1,ρ)
end

Spherical_harmonics

function Y(l,m,θ,φ)
    N = (im)^(m+abs(m)) * sqrt( (2*l+1)*factorial(l-Int(abs(m))) / (2*factorial(l+Int(abs(m)))) )
    return N * P(l,Int(abs(m)),cos(θ)) * exp(im*m*φ) / sqrt(2*π)
end

Total_wave_function

# ψ(n,l,m,r,θ,φ) = R(n,l,r)*Y(m,l,θ,φ)
ψ(n,l,m,x,y,z) = R(n,l,r(x,y,z))*Y(l,m,θ(x,y,z),φ(x,y,z))

Supplementary Information

A few additional functions are needed, so I will provide them here.

Notes on Radial Functions

The definitions used vary depending on whether the source is Japanese or English, by a mathematician or a physicist, or whether it's an old or new reference. Since this affects library implementation, it is important to understand this. First, as shown in the table below, there are at least four patterns for the radial function $R_{nl}(r)$ .

Pattern	Definition of Laguerre Polynomial $L_n(x)$	$L_n^{k}(x) = L_{n-k}^{(k)}(x)$ or adoption of $L_n^{(\alpha)}(x)$	Denominator of normalization factor	References
1	$\frac{1}{n!} \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)$	$L_{n+l}^{2l+1}(\rho)$	$2n(n+l)!$
2	$\frac{1}{n!} \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)$	$L_{n-l-1}^{(2l+1)}(\rho)$	$2n(n+l)!$	1, 2
3	$\mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)$	$L_{n+l}^{2l+1}(\rho)$	$2n[(n+l)!]^3$	Many
4	$\mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)$	$L_{n-l-1}^{(2l+1)}(\rho)$	$2n[(n+l)!]^3$

Broadly speaking, they can be categorized into four patterns based on: first, whether the definition of the Laguerre polynomial $L_n(x)$ includes $1/n!$ ; and second, whether associated Laguerre polynomials $L_n^{k}(x)$ or generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ are used.

Pattern 3 is common and can be found in "McQuarrie & Simon Physical Chemistry" p.224, "Atkins' Physical Chemistry (10th Ed)" p.381, Yoshiya Harada's "Quantum Chemistry Vol. 1" (2nd Ed, 2019, Shokabo) p.123, Kazuo Takayanagi's "Atomic and Molecular Physics" (2000, Asakura Shoten) p.16, Teiji Kunihiro's "Quantum Mechanics" (Tokyo Tosho, 2018) p.125, Wikipedia (Japanese), etc.

However, Pattern 2 leads to the most natural implementation in SpecialPolynomials.jl. In that library, the definition of the Laguerre polynomial $L_n(x)$ includes $1/n!$ , associated Laguerre polynomials $L_n^{k}(x)$ are not supported, and generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ are supported. Since there are few Japanese references for generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ , I will explain them here.

Polynomial	Equation to satisfy	References
Laguerre polynomials $L_n(x)$	$xy'' + (1-x)y' + ny = 0$
Associated Laguerre polynomials $L_n^k(x)$	$xy'' + (k+1-x)y' + (n-k)y = 0$	Japanese Wikipedia
Generalized Laguerre polynomials $L_n^{(\alpha)}(x)$	$xy'' + (\alpha+1-x)y' + ny = 0$	English Wikipedia

As shown in the table above, Laguerre polynomials $L_n(x)$ , associated Laguerre polynomials $L_n^{k}(x)$ , and generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ all satisfy different equations. Regarding the first issue—whether to include $1/n!$ in the definition of $L_n(x)$ —since they are all solutions to linear differential equations, the difference of a constant multiple is ultimately a matter of preference. However, most libraries seem to adopt the definition that includes $1/n!$ . As for the second issue—whether to use associated Laguerre polynomials $L_n^{k}(x) = L_{n-k}^{(k)}(x)$ or generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ —within the range of $k \in \mathbb{N}$ , they are linked by the relation $L_n^{k}(x) = L_{n-k}^{(k)}(x)$ . Therefore, either choice is fine and doesn't affect the normalization factor. (For the Morse potential, where $k \in \mathbb{R}$ must be considered, generalized Laguerre polynomials are typically used). SpecialPolynomials.jl supports generalized Laguerre polynomials $L_n^{(\alpha)}(x)$ .

Associated Laguerre polynomials $L_n^k(x)$

While associated Laguerre polynomials $L_n^k(x)$ are not supported in SpecialPolynomials.jl, generalized Laguerre polynomials $L_n^{(\alpha)}$ are supported.

We generate associated Laguerre polynomials using the relationship:

L_n^{k}(x) = L_{n-k}^{(k)}(x)

Also, since the definition of the Laguerre polynomial $L_n(x)$ in SpecialPolynomials.jl is:

L_n(x) = \frac{1}{n!}\mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)

it is necessary to multiply by $n!$ to match the definition used for the radial functions explained earlier:

L_n(x) = \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)

Associated_Laguerre_polynomials_L_n^k(x)

# Provided that the Laguerre polynomial L_n(x) does not include 1/n!
# using SpecialPolynomials

L(n,k,x) = factorial(n) * basis(Laguerre{k}, n-k)(x)

Notes on Spherical Harmonics

The definition of the associated Legendre polynomials $P_k^m(t)$ appears to be standardized, and I have yet to observe differences arising from its definition. However, the definition of the coefficients for the spherical harmonics themselves has many variants, requiring even more caution than the radial functions.

Phase Factor	References
$$i^{	m
$$(-1)^{(m+	m
$(-1)^m~(m>0), 1~(m\leq0)$	This page, Inoki/Kawai. Same as above?
$(-1)^m$	Teiji Kunihiro, "Quantum Mechanics" (Tokyo Tosho, 2018), Wikipedia - Spherical Harmonics (English)

Landau-Lifshitz includes specific expressions for spherical harmonics in the appendix, and since they include complex numbers, they match.
Does the $(-1)^m$ pattern match the above if $P_l^{m}$ is used instead of $P_l^{|m|}$ ? All $|m|$ inside the $\sqrt{\cdot}$ would need to be replaced with $m$ .

There are many other variations. I believe I've put in a good effort, so I'll leave the investigation of other "ecosystems" to the reader. Did you understand that you need to carefully verify the properties whenever you actually use them?

Associated Legendre Polynomials $P_k^m(t)$

Associated Legendre polynomials $P_k^m(t)$ also don't seem to be supported in SpecialPolynomials.jl, so I created a function using polynomial differentiation in collaboration with Polynomials.jl. However, performing polynomial differentiation on every call is very slow, so I decided to pre-calculate and store values up to $k=7$ in memory.

Associated_Legendre_polynomials_P_k^m(t)

# using SpecialPolynomials
# using Polynomials
function P(k,m)
    #if m<0
    #    return (-1)^m * factorial(l-m) / factorial(l+m) * P(k,m,x)
    #end
    p = basis(Legendre, k)
    q = convert(Polynomial, p)
    r = derivative(q, m)
    return x -> (1-x^2)^(m/2) * r(x)
end

LegendreArray = Dict()

for k in 0:7 # Pre-store up to k=7 in memory.
    for m in 0:k
        LegendreArray[k,m] = P(k,m)
    end
end

P(k,m,x) = LegendreArray[k,m](x)

Spherical Coordinates

The transformation between Cartesian and spherical coordinates is as described on Wikipedia:

\begin{cases} x = r\sin\theta\cos\varphi\\ y = r\sin\theta\sin\varphi\\ z = r\cos\theta \end{cases}

\begin{cases} r = \sqrt{x^2+y^2+z^2}\\ \theta = \arccos(z/\sqrt{x^2+y^2+z^2})\\ \varphi = \mathrm{sgn}(y)\arccos(x/\sqrt{x^2+y^2}) \end{cases}

However, you notice issues when implementing this. First, since the domain of $\arccos(x)$ is $-1\leq x \leq 1$ , exceeding this domain results in an error. This problem likely occurs in most programming languages, not just Julia. It is best to extend the domain gracefully by adding periodicity. In addition to extending the domain, some cleverness is required to prevent divergence, such as setting $\theta\sim\arccos(1)=0$ when $x^2+y^2\sim0$ , and $\varphi\sim\arccos(1)=0$ when $y\sim0$ .

x(r,θ,φ) = r*sin(θ)*cos(φ)
y(r,θ,φ) = r*sin(θ)*sin(φ)
z(r,θ,φ) = r*cos(θ)
r(x,y,z)  = sqrt(x^2+y^2+z^2)
θ(x,y,z) = x^2+y^2<1e-9 ? 0.0 : myacos(z/r(x,y,z)) 
φ(x,y,z) = y^2<1e-9 ? 0.0 : sign(y)*myacos(x/sqrt(x^2+y^2))
# θ(x,y,z) = acos(z/r(x,y,z))
# φ(x,y,z) = sign(y)*acos(x/sqrt(x^2+y^2))

The myacos(x) function above is a function that adds periodicity to acos(x).

Extension_of_acos(x)

loop(x) = x<-1 ? loop(x+2) : (1<x ? loop(x-2) : x)
myacos(x) = acos(loop(x))

# Test
println("x\tloop(x)\tacos(x)\tmyacos(x)")
for x in -3:0.2:3
    @printf("%.2f\t%.2f\t%.2f\t%.2f\n", x, loop(x), (try; acos(x); catch; NaN; end), myacos(x) )
end

Output

    x	loop(x)	acos(x)	myacos(x)
    -3.00	-1.00	NaN	3.14
    -2.80	-0.80	NaN	2.50
    -2.60	-0.60	NaN	2.21
    -2.40	-0.40	NaN	1.98
    -2.20	-0.20	NaN	1.77
    -2.00	0.00	NaN	1.57
    -1.80	0.20	NaN	1.37
    -1.60	0.40	NaN	1.16
    -1.40	0.60	NaN	0.93
    -1.20	0.80	NaN	0.64
    -1.00	-1.00	3.14	3.14
    -0.80	-0.80	2.50	2.50
    -0.60	-0.60	2.21	2.21
    -0.40	-0.40	1.98	1.98
    -0.20	-0.20	1.77	1.77
    0.00	0.00	1.57	1.57
    0.20	0.20	1.37	1.37
    0.40	0.40	1.16	1.16
    0.60	0.60	0.93	0.93
    0.80	0.80	0.64	0.64
    1.00	1.00	0.00	0.00
    1.20	-0.80	NaN	2.50
    1.40	-0.60	NaN	2.21
    1.60	-0.40	NaN	1.98
    1.80	-0.20	NaN	1.77
    2.00	0.00	NaN	1.57
    2.20	0.20	NaN	1.37
    2.40	0.40	NaN	1.16
    2.60	0.60	NaN	0.93
    2.80	0.80	NaN	0.64
    3.00	1.00	NaN	0.00

Orbital Names

Finally, I created a function to map quantum numbers to orbital names, which turned out to be surprisingly useful.

# Orbital names
orbital_name(n,l) = "$n$(["s","p","d","f","g","h"][l+1])"

Verification of Expectation Values

To calculate various expectation values, I prepared a simple numerical integration routine and tested it by calculating a Gaussian integral.

Numerical_integration_routine_(1D)

integral(func, x_min, x_max, dx) = sum(func.(x_min:dx:x_max))*dx
# Test (Gaussian integral)
integral(x->exp(-x^2)/sqrt(pi), -10,10,0.1) |> display

Output

1.0

Numerical_integration_routine_(2D)

function integral2D(func, x_min, x_max, dx, y_min, y_max, dy)
    sum = 0.0
    for x in x_min:dx:x_max
        for y in y_min:dy:y_max
            sum += func(x,y)*dx*dy
        end
    end
    return sum
end
# Test (Gaussian integral)
integral2D((x,y)->exp(-x^2-y^2)/sqrt(pi^2), -10,10,0.1, -10,10,0.1) |> display

Output

    0.9999999999999957

Normalization Condition

Since it is often difficult to determine whether eigenfunctions are declared correctly just by looking at them, I will first perform a simple calculation to check if the normalization condition is satisfied. This also serves to verify that the library usage is correct, that there are no bugs, and that there are no errors in the literature. Because this is numerical integration, it will not be exactly 1, but it is fine as long as it is approximately 1.

For the radial part, the following equation should hold:

\int_0^\infty |R_{nl}(r)|^2 r^2 \mathrm{d}r = 1

Normalization_condition_for_radial_functions

println("\t∫r^2R^2(r)dr")
for n in 1:5
    for l in 0:n-1
       @printf("%s\t%.17f\n", orbital_name(n,l), integral(r->r^2*abs(R(n,l,r))^2, 0, 100, 0.01))
    end
end

Output

    	∫r^2R^2(r)dr
    1s	0.99999999933335459
    2s	0.99999999991666833
    2p	0.99999999999999989
    3s	0.99999999997530897
    3p	1.00000000000000000
    3d	1.00000000000000000
    4s	0.99999999998488820
    4p	0.99999999999693834
    4d	0.99999999999879452
    4f	0.99999999999977918
    5s	0.99999930806578163
    5p	0.99999948583752063
    5d	0.99999972631969702
    5f	0.99999990537286676
    5g	0.99999998382258726

For the angular part, the following equation should hold:

\int_0^\pi \int_0^{2\pi} |Y_{lm}(\theta,\varphi)|^2 \sin(\theta) \mathrm{d}\theta \mathrm{d}\varphi = 1

Normalization_condition_for_spherical_harmonics

println("l\tm\t∫|Y(θ,φ)|^2sinθdθdφ")
for l in 1:3
    for m in -l:l
       @printf("%d\t%d\t%.17f\n", l, m, integral2D((θ,φ)->abs(Y(l,m,θ,φ))^2*sin(θ), 0,π,0.01, 0,2*π,0.01))
    end
end

Output

    l	m	∫|Y(θ,φ)|^2sinθdθdφ
    1	-1	1.00108459213941758
    1	0	1.00106961810439521
    1	1	1.00108459213941758
    2	-2	1.00108459204766698
    2	-1	1.00108459250535242
    2	0	1.00105963486609806
    2	1	1.00108459250535242
    2	2	1.00108459204766698
    3	-3	1.00108459204825762
    3	-2	1.00108459204782041
    3	-1	1.00108459332884259
    3	0	1.00104965071179763
    3	1	1.00108459332884259
    3	2	1.00108459204782041
    3	3	1.00108459204825762

I am concerned that the errors are somewhat large, but this is likely due to the poor description of the boundaries in the integration domain. However, as shown next, it seems to satisfy the orthonormalization condition for now, so we will proceed.

Orthogonality

While we're at it, let's also verify the orthogonality.

\int_0^\pi \int_0^{2\pi} Y^\ast_{lm}(\theta,\varphi) Y_{l'm'}(\theta,\varphi) \sin(\theta) \mathrm{d}\theta \mathrm{d}\varphi = \delta_{ll'} \delta_{mm'} \tag{2.38}

Orthonormality_of_spherical_harmonics

println("l\tl'\t∫Y*(θ,φ)Y(θ,φ)sinθdθdφ")
m = 0
for l in 1:4
    for ll in 1:4
       @printf("%d\t%d\t%.17f\n", l, ll, integral2D((θ,φ)->conj(Y(l,m,θ,φ))*Y(ll,m,θ,φ)*sin(θ), 0,π,0.01, 0,2*π,0.01))
    end
end

Output

    l	l'	∫Y*(θ,φ)Y(θ,φ)sinθdθdφ
    1	1	1.00106961810439521
    1	2	-0.00001297905216095
    1	3	-0.00002287377512273
    1	4	-0.00001741363401505
    2	1	-0.00001297905216095
    2	2	1.00105963486609806
    2	3	-0.00001982616166688
    2	4	-0.00003348500577979
    3	1	-0.00002287377512273
    3	2	-0.00001982616166688
    3	3	1.00104965071179763
    3	4	-0.00002660021073039
    4	1	-0.00001741363401505
    4	2	-0.00003348500577979
    4	3	-0.00002660021073039
    4	4	1.00103966527685562

Virial Theorem

In the Coulomb potential, since $2\langle T \rangle = -\langle V \rangle$ holds, it immediately follows that $E = \langle T \rangle + \langle V \rangle = \langle V \rangle/2$ . We will verify the virial theorem by calculating $\langle V \rangle$ using the following equation:

\langle V \rangle = \int_0^\infty V(r) |R_{nl}(r)|^2 r^2 \mathrm{d}r

Verification_of_the_virial_theorem

println("\tE\t\t\t6V7/2")
for n in 1:5
    for l in 0:n-1
        @printf("%s\t%.17f\t%.17f\n", orbital_name(n,l), E(n), integral(r->V(r)*r^2*abs(R(n,l,r))^2, 0.01, 100, 0.01)/2)
    end
end

Output

    	E			<V>/2
    1s	-0.50000000000000000	-0.49998333366666148
    2s	-0.12500000000000000	-0.12499791670312468
    2p	-0.12500000000000000	-0.12500000000173608
    3s	-0.05555555555555555	-0.05555493828212155
    3p	-0.05555555555555555	-0.05555555555616521
    3d	-0.05555555555555555	-0.05555555555555557
    4s	-0.03125000000000000	-0.03124973958770489
    4p	-0.03125000000000000	-0.03125000000025629
    4d	-0.03125000000000000	-0.03124999999999411
    4f	-0.03125000000000000	-0.03124999999999892
    5s	-0.02000000000000000	-0.01999986331948914
    5p	-0.02000000000000000	-0.01999999751095307
    5d	-0.02000000000000000	-0.01999999867474622
    5f	-0.02000000000000000	-0.01999999954163147
    5g	-0.02000000000000000	-0.01999999992160406

Expectation Value $\langle r \rangle$

Various expectation values are introduced on this site and in Takayanagi (2000) p. 20, among other sources.

Expectation_value_<i>

r(n,l;Z=1) = (3*n^2-l*(l+1))/2

println("\texact\t\t\t6r7")
for n in 1:6
    for l in 0:n-1
        @printf("%s\t%2.17f\t%2.17f\n", orbital_name(n,l), r(n,l), integral(r->r^3*abs(R(n,l,r))^2, 0.1, 1000, 0.01))
    end
end

Output

    	exact			<r>
    1s	1.50000000000000000	1.49993035023037469
    2s	6.00000000000000000	5.99999130697251903
    2p	5.00000000000000000	4.99999999535692119
    3s	13.50000000000000000	13.49999742501150735
    3p	12.50000000000000000	12.49999999836981246
    3d	10.50000000000000000	10.49999999999993605
    4s	24.00000000000000000	23.99999891378351080
    4p	23.00000000000000000	22.99999999927470284
    4d	21.00000000000000000	20.99999999999996803
    4f	18.00000000000000000	18.00000000000000355
    5s	37.50000000000000000	37.49999944388246575
    5p	36.50000000000000000	36.49999999961975305
    5d	34.50000000000000000	34.49999999999998579
    5f	31.50000000000000000	31.50000000000001066
    5g	27.50000000000000000	27.50000000000000000
    6s	54.00000000000000000	53.99999967818066438
    6p	53.00000000000000000	52.99999999977713117
    6d	51.00000000000000000	50.99999999999999289
    6f	48.00000000000000000	48.00000000000000711
    6g	44.00000000000000000	16094.43124042773706606
    6h	39.00000000000000000	14506950.32743391394615173
"

Visualization

Energy Eigenvalues

First, let's plot the potential and the eigenvalues.

Potential_and_Eigenvalues

# using Plots

plot(xlabel="$r$", ylabel="$V(r)$", xlims=(-20,20), ylims=(-0.55,0.05), legend=:bottomright) #, xticks=-15:1:15)
plot!(-50:0.1:50, r->V(r), label="V(r)", lc="#000000", lw=2)
for n in 1:10
    plot!([1/E(n), -1/E(n)], fill(E(n),2), lc=n, lw=2, label="n = $n")
end
plot!()

Energy levels

Spectral Series

Spectra can be derived from the differences between energy eigenvalues, allowing for direct comparison with experiments. Please refer to here.

Wavelength_Calculation

# https://physics.nist.gov/cgi-bin/cuu/CCValue?hrminv
Eh2m1 = 2.1947463136320*1e7
Eh2cm1 = 2.1947463136320*1e5
Eh2nm1 = 2.1947463136320*1e-2

frequency(n,m;Z=1) = (E(n,Z=Z) - E(m,Z=Z))*Eh2cm1
wavelength(n,m;Z=1) = 1/((E(n,Z=Z) - E(m,Z=Z))*Eh2nm1)

println("\tLymann\t\tBalmer\t\tPaschen")
println("n -> m\tm=1\t\tm=2\t\tm=3")
@printf("n = %-3d\t%.2f\tnm\n", 2, wavelength(2,1))
@printf("n = %-3d\t%.2f\tnm\t%.2f\tnm\n", 3, wavelength(3,1), wavelength(3,2))
for i in [4,5,6,7,8,9,10,Inf]
    @printf("n = %-3d\t%.2f\tnm\t%.2f\tnm\t%.2f\tnm\n", i, wavelength(i,1), wavelength(i,2), wavelength(i,3))
end

Output

    	Lymann		Balmer		Paschen
    n -> m	m=1		m=2		m=3
    n = 2  	121.50	nm
    n = 3  	102.52	nm	656.11	nm
    n = 4  	97.20	nm	486.01	nm	1874.61	nm
    n = 5  	94.92	nm	433.94	nm	1281.47	nm
    n = 6  	93.73	nm	410.07	nm	1093.52	nm
    n = 7  	93.03	nm	396.91	nm	1004.67	nm
    n = 8  	92.57	nm	388.81	nm	954.35	nm
    n = 9  	92.27	nm	383.44	nm	922.66	nm
    n = 10 	92.05	nm	379.69	nm	901.25	nm
    n = Inf	91.13	nm	364.51	nm	820.14	nm

# using Plots

Plots.plot(size=(800, 150), yshowaxis=false, grid=false, xscale=:log10, ylims=(0,3))
xticks!([100,1000,10000], ["100nm", "1000nm", "10000nm"])

for j in 1:6
    plot!(annotations=(wavelength(j+1,j), 2.9, (["Lyman", "Balmer ", "Paschen ", "Brackett ", "Pfund ", "Humphreys "][j], 8, 0.0, :right)))
    plot!(annotations=(wavelength(j+1,j), 2.5, (@sprintf("~%.2f", wavelength(j+1,j)), 8, 0.0, :right)))
    for i in j+1:15
        if i==j+1
            plot!(fill(wavelength(i,j),2), [0,2.3], label="", lc=j)
        else    
            plot!(fill(wavelength(i,j),2), [0,2], label="", lc=j)
        end
    end
end

plot!()

Spectrum

Radial Functions

The radial density functions are plotted below. Please compare them with this page or Wikipedia to verify the positions of the nodes and antinodes.

Plotting_Radial_Density_Functions

# using Plots

plot(xlabel="$r$", ylabel="$r^2|R_{nl}(r)|^2$", ylims=(-0.01,0.55), xticks=0:1:21)
for n in 1:3
    for l in 0:n-1
        plot!(0:0.1:21, r->r^2*R(n,l,r)^2, label=orbital_name(n,l), lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1])
    end
end
plot!()

Radial distribution

Spherical Harmonics

Since spherical harmonics are complex-valued, some ingenuity is required to create real-valued functions like the commonly seen $2p_{xy}$ or $3d_{xy}$ . This explanation was very helpful regarding specific expressions and rendering methods. The actual plotting was easily achieved using the PyPlot backend of Plots.jl. For phase color-coding, I specified fill_z=C referring to this article.

Plotting_Spherical_Harmonics

pyplot()

function my_surface(N, M, f; xlims=(-0.2,0.2), ylims=(-0.2,0.2), zlims=(-0.2,0.2), clims=(-0.5,0.5), title="")
    
    Θ = range(0, stop=π,   length=N)
    Φ = range(0, stop=2*π, length=M)

    p(θ,φ) = abs(f(θ,φ))^2
    
    X = zeros(N, M)
    Y = zeros(N, M)
    Z = zeros(N, M)
    C = zeros(N, M)

    for j in 1:M
        for i in 1:N
            X[i,j] = p(Θ[i],Φ[j]) * sin(Θ[i])*cos(Φ[j])
            Y[i,j] = p(Θ[i],Φ[j]) * sin(Θ[i])*sin(Φ[j])
            Z[i,j] = p(Θ[i],Φ[j]) * cos(Θ[i])
            C[i,j] = real(f(Θ[i],Φ[j]))
        end
    end
    
    plt = surface(X, Y, Z, fill_z=C, xlims=xlims, ylims=ylims, zlims=zlims, clims=clims, title=title, c=:Spectral)
    plot(plt) |> display
    return plt
end

f(θ,φ) = Y(0,0,θ,φ)
my_surface(100, 200, f, title=L"1s")

f(θ,φ) = (Y(1,-1,θ,φ) - Y(1,1,θ,φ)) / sqrt(2)
my_surface(100, 200, f, title=L"2p_x")

f(θ,φ) = (Y(1,-1,θ,φ) + Y(1,1,θ,φ)) * im / sqrt(2)
my_surface(100, 200, f, title=L"2p_y")

f(θ,φ) = Y(1,0,θ,φ)
my_surface(100, 200, f, title=L"2p_z")

f(θ,φ) = (Y(2,-2,θ,φ) - Y(2,2,θ,φ)) * im / sqrt(2)
my_surface(100, 200, f, title=L"d_{xy}")

f(θ,φ) = (Y(2,-1,θ,φ) + Y(2,1,θ,φ)) * im / sqrt(2)
my_surface(100, 200, f, title=L"d_{yz}")

f(θ,φ) = (Y(2,-1,θ,φ) - Y(2,1,θ,φ)) / sqrt(2)
my_surface(100, 200, f, title=L"d_{zx}")

f(θ,φ) = (Y(2,-2,θ,φ) + Y(2,2,θ,φ)) / sqrt(2)
my_surface(100, 200, f, title=L"d_{x^2-y^2}")

f(θ,φ) = Y(2,0,θ,φ)
my_surface(100, 200, f, title=L"d_{3z-r^2}")

Spherical harmonics

Orbital Visualization

The spherical harmonic plots above only represent the dependence on the angles $\theta$ and $\varphi$ . Since I couldn't find a good way to render isosurfaces that also include the contribution of the radial part, I will use GR.jl to render isosurfaces from a 3D array A[:,:,:]. Since Plots.jl and GR.jl often conflict, it is better to declare using GR within a module.

Convert_functions_to_arrays

module Orbital
    using GR
    GR.__init__()

    function isosurface(x_min, x_max, dx, func; isovalue=NaN, rotation=180 , tilt=150)
        s = x_min:dx:x_max # Sample points at dx intervals for x_min ≦ x ≦ x_max in all 3 axes
        N = length(s)
        array = zeros(N,N,N)
        for k in 1:N
            for j in 1:N
                for i in 1:N
                    array[i,j,k] = func(s[i],s[j],s[k])
                end
            end
        end
        if isnan(isovalue) || isovalue<minimum(array) || maximum(array)<isovalue
            isovalue = sum(array)/length(array)
        end
        GR.isosurface(array, isovalue=isovalue, rotation=rotation , tilt=tilt) |> display
    end
end

# # Test (Plotting a d-type Gaussian orbital)
# f(x,y,z) = (x*y*exp(-(x^2+y^2+z^2)))^2
# Orbital.isosurface(-2, 2, 0.1, f)

L"1s" |> display
f(x,y,z) = abs(ψ(1,0,0,x,y,z))^2
Orbital.isosurface(-10, 10, 0.5, f)

L"2p_x" |> display
f(x,y,z) = abs(ψ(2,1,1,x,y,z) + ψ(2,1,-1,x,y,z))^2 / 2
Orbital.isosurface(-10, 10, 0.5, f)

L"2p_y" |> display
f(x,y,z) = abs(ψ(2,1,1,x,y,z) - ψ(2,1,-1,x,y,z))^2 / 2
Orbital.isosurface(-10, 10, 0.5, f)

L"2p_z" |> display
f(x,y,z) = abs(ψ(2,1,0,x,y,z))^2
Orbital.isosurface(-10, 10, 0.5, f)

L"3d_{xy}" |> display
f(x,y,z) = abs(ψ(3,2,-2,x,y,z) - ψ(3,2,2,x,y,z))^2 / 2
Orbital.isosurface(-20, 20, 1.0, f)

L"3d_{yz}" |> display
f(x,y,z) = abs(ψ(3,2,-1,x,y,z) + ψ(3,2,1,x,y,z))^2 / 2
Orbital.isosurface(-20, 20, 1.0, f)

L"3d_{zx}" |> display
f(x,y,z) = abs(ψ(3,2,-1,x,y,z) - ψ(3,2,1,x,y,z))^2 / 2
Orbital.isosurface(-20, 20, 1.0, f)

L"3d_{x^2-y^2}" |> display
f(x,y,z) = abs(ψ(3,2,-2,x,y,z) + ψ(3,2,2,x,y,z))^2 / 2
Orbital.isosurface(-20, 20, 1.0, f)

L"3d_{3z^2-r^2}" |> display
f(x,y,z) = abs(ψ(3,2,0,x,y,z))^2
Orbital.isosurface(-20, 20, 1.0, f)

$1s$
1s orbital
$2p_x$
2p_x orbital
$2p_y$
2p_y orbital
$2p_z$
2p_z orbital
$3d_{xy}$
$3d_{xy} orbital$
$3d_{yz}$
$3d_{yz} orbital$
$3d_{zx}$
$3d_{zx} orbital$
$3d_{x^2-y^2}$
$3d_{x^2-y^2} orbital$
$3d_{3z^2-r^2}$
$3d_{3z^2-r^2} orbital$

Of course, there are also 2s and 3s, but they are omitted because the differences are difficult to see without a library that can display cross-sections or represent transparency.

Conclusion

I caught a glimpse of the complex circumstances surrounding special functions. It might be a conflict between physicists wanting to stick to established definitions and mathematicians wanting to replace them with more generalized concepts. I will summarize how to use special functions in Julia in another article.

I didn't expect that rendering isosurfaces would be this challenging even for the simplest atom, the hydrogen atom. If a library for easily plotting isosurfaces is developed, it will likely advance the use of Julia in computational physics and chemistry.

Environment

versioninfo()

Output

    Julia Version 1.6.2
    Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
    Platform Info:
      OS: Windows (x86_64-w64-mingw32)
      CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
      WORD_SIZE: 64
      LIBM: libopenlibm
      LLVM: libLLVM-11.0.1 (ORCJIT, haswell)

References

The Jupyter Notebook data that served as the basis for this article is available at the link below.

iTranslated by AI

Visualization with Julia: Hydrogen Atom

Packages

Schrödinger Equation

Supplementary Information

Notes on Radial Functions

Associated Laguerre polynomials $L_n^k(x)$

Notes on Spherical Harmonics

Associated Legendre Polynomials $P_k^m(t)$

Spherical Coordinates

Orbital Names

Verification of Expectation Values

Normalization Condition

Orthogonality

Virial Theorem

Expectation Value $\langle r \rangle$

Visualization

Energy Eigenvalues

Spectral Series

Radial Functions

Spherical Harmonics

Orbital Visualization

Conclusion

Environment

References

Visualization with Julia Series

Discussion

iTranslated by AI

Packages

Schrödinger Equation

Supplementary Information

Notes on Radial Functions

Associated Laguerre polynomials L_n^k(x)

Notes on Spherical Harmonics

Associated Legendre Polynomials P_k^m(t)

Spherical Coordinates

Orbital Names

Verification of Expectation Values

Normalization Condition

Orthogonality

Virial Theorem

Expectation Value \langle r \rangle

Visualization

Energy Eigenvalues

Spectral Series

Radial Functions

Spherical Harmonics

Orbital Visualization

Conclusion

Environment

References

Visualization with Julia Series

Discussion

Associated Laguerre polynomials $L_n^k(x)$

Associated Legendre Polynomials $P_k^m(t)$

Expectation Value $\langle r \rangle$