Atomic Orbitals

This notebook is based on Christina C. Lee's blog post (see https://albi3ro.github.io/M4/prerequisites/Atomic-Orbitals.html). A similar solution for Python can be found here: https://github.com/jheinen/gr/blob/master/examples/orbitals.py

What's below is a bunch of function definitions that make our calculations easier later on. Here we utilize the GNU scientific library (GSL) to calculate the special functions.

In [1]:
using GSL
using GR

# Bohr radius
a0 = 5.291772109217e-11;
\begin{equation} \rho=\frac{2r}{n a_0} \end{equation}
In [2]:
# The radial coordinate
ρ(r,n) = 2r/(n*a0);
\begin{equation} R^{n,l} ( \rho ) = L^{2 l+1}_{n-l-1} ( \rho ) \cdot e^{- \rho /2} \cdot \rho ^l \end{equation}

where $L^{2 l+1}_{n-l-1}(\rho)$ is the generalized Laguerre polynomial.

In [3]:
# The Radial dependence
function R(n::Int,l::Int,ρ::Real)
    sf_laguerre_n(n-l-1,2*l+1,ρ) * e^(-ρ/2)*ρ^l
end;
\begin{equation} Y^m_l(θ,ϕ) = (-1)^m \cdot P^m_l (\cos(θ)) \cdot e^{i m \phi} \end{equation}

where $P^m_l (\cos (\theta))$ is the associated Legendre polynomial.

In [4]:
# The θ and ϕ dependence
function Yml(m::Int,l::Int,θ::Real,ϕ::Real)
    (-1.0)^m * sf_legendre_Plm(l,abs(m),cos(θ)) * e^(im*m*ϕ)
end;
\begin{equation} N(n,l) = \sqrt{\left(\frac{2}{n}\right)^3 \frac{(n-l-1)!}{2n(n+l)!}} \end{equation}
In [5]:
# A normalization: This is dependent on the choice of polynomial representation
function norm(n::Int,l::Int)
    sqrt((2/n)^3 * factorial(n-l-1)/(2n*factorial(n+l)))
end;

Guess a solution with separated radial and angular variables: \begin{equation} \Psi(\rho,\theta,\phi) = N(n, l) \cdot R^{n,l} (\rho) \cdot Y^m_l(\theta,\phi) \end{equation}

In [6]:
# Generates an Orbital Funtion of (r,θ,ϕ) for a specificied n,l,m.
function Orbital(n::Int,l::Int,m::Int)
    # we make sure l and m are within proper bounds
    if l > n || abs(m) > l
        throw(DomainError())
    end
    Ψ(ρ,θ,ϕ) = norm(n,l) * R(n,l,ρ) * abs(Yml(m,l,θ,ϕ))
    Ψ
end;

In calculate_electronic_density(), we create a square cube, convert those positions over to spherical coordinates, and calculate the electronic density.

In [7]:
function calculate_electronic_density(n,l,m)
    N = 50
    r = 1e-9
    s = linspace(-r,r,N)
    x,y,z = meshgrid(s,s,s)
    
    ϕ,θ,r = cart2sph(x,y,z)
    # for reasons of compatibility with MATLAB, θ is the counterclockwise
    # angle in the x-y plane measured in radians from the positive x-axis
    θ -= π/2
    Ψ = Orbital(n,l,m)

    Ψv = zeros(Float32,N,N,N)
    for i in 1:N
        for j in 1:N
            for k in 1:N
                Ψv[i,j,k] = abs(Ψ(ρ(r[i,j,k],n),θ[i,j,k],ϕ[i,j,k]))
            end
        end
    end
    Ψv
end;

Use GR inline() graphics to create an animation on the fly. Visualize and rotate a 3d Orbital (corresponds to $n=3, l=2, m=0$) using isosurface().

In [8]:
inline("mov")

Ψv = calculate_electronic_density(3, 2, 0)

for alpha in 30:210
    isosurface(Ψv, isovalue=0.25, rotation=alpha)
end

Show the resulting movie in the browser.

In [9]:
GR.show()
Out[9]: