### 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

a0 = 5.291772109217e-11;

$$\rho=\frac{2r}{n a_0}$$
In [2]:
# The radial coordinate
ρ(r,n) = 2r/(n*a0);

$$R^{n,l} ( \rho ) = L^{2 l+1}_{n-l-1} ( \rho ) \cdot e^{- \rho /2} \cdot \rho ^l$$

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;

$$Y^m_l(θ,ϕ) = (-1)^m \cdot P^m_l (\cos(θ)) \cdot e^{i m \phi}$$

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;

$$N(n,l) = \sqrt{\left(\frac{2}{n}\right)^3 \frac{(n-l-1)!}{2n(n+l)!}}$$
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: $$\Psi(\rho,\theta,\phi) = N(n, l) \cdot R^{n,l} (\rho) \cdot Y^m_l(\theta,\phi)$$

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]: