Appendix 2: Visualizing Atomic Orbitals

Note

This appendix assumes a future version of SymPy for the Z_lm() function. This function has been temporarily defined in a code cell below to provide this feature until the next SymPy release.

The visualization of atomic orbitals and orbital information is important enough topics in chemistry to warrant specific attention. This appendix focuses on different methods of visualizing various aspects of atomic orbitals and tools to assist in this task. This content is not included in the chapter on plotting with matplotlib because this appendix heavily utilizes various libraries such as SymPy, interact, and NumPy not yet introduced before Chapter 03. While this appendix is written to be standalone as much as possible, knowledge of matplotlib, including surface plots, will be helpful along with NumPy and SymPy basics.

Atomic orbitals are described by a wavefunction, \(\Psi (n, l, m)\), which is the product of the radial wavefunction, \(R(n, l)\), and the angular wavefunction, \(Y(l, m)\). Each atomic orbital has a different wavefunction \(\Psi\), but they sometimes share common radial wavefunctions.

\[ \Psi(n, l, m) = R(n, l)Y(l, m) \]

The radial wavefunction depends upon the principal (n) and angular (\(l\)) quantum numbers and provides information about the wavefunction or electron probability at various distances from the nucleus. The radial wavefunction is independent of the direction. The angular wavefunction describes the direction of the orbital with respect to the spherical coordinate angles and depends upon the angular (\(l\)) and magnetic (\(m\) or \(m_l\)) quantum numbers. We will first visualize the radial and angular components individually before combining them into a more complete picture of atomic orbitals.

We will use NumPy and matplotlib heavily in this chapter, and we will make heavy use of the SymPy library for convenient functions in its hydrogen module. These are all imported below.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import sympy
from sympy.physics.hydrogen import R_nl, Psi_nlm #, Z_lm

# delete this cell and replace with actual Z_lm after next SymPy release
from sympy.functions.special.spherical_harmonics import Znm
def Z_lm(l, m, phi, theta):
    return Znm(l, m, theta, phi).expand(func=True)

Note

The SymPy library is introduced in Chapter 8 and provides mainly tools for symbolic mathematics along with other tools for wavefunctions, harmonic oscillators, biomechanics, etc.

Radial Wavefunctions

Because the radial wavefunctions are independent of direction, they can be represented effectively on a simple 2D plot. The toughest part is coding the equations for every combination of n and \(l\). The good news is that the SymPy library includes a function, R_nl(), in the Hydrogen Wavefunction (sympy.physics.hydrogen) module that provides this functionality. This function takes the principle quantum number (n), angular quantum number (\(l\)), radius in Bohrs (r), and atomic number (Z). A Bohr equals about 52.9 pm.

R_nl(n, l, m, r, Z=1)

We can evaluate the function for any hydrogen-like atomic orbitals such as the 3p orbitals (n = 3 and \(l\) = 1) at 4.0 Bohrs.

R_nl(3, 1, 4.0, Z=1)

\[\displaystyle 0.0173561901639985 \sqrt{6}\]

SymPy prefers to return results in exact form, so it includes \(\sqrt6\) in this particular result. To get a float answer, use the evalf() method.

Note

The evalf() method can take an optional argument for the precision number such as evalf(5) for 5 digits of precision.

R_nl(3, 1, 4.0, Z=1).evalf()

\[\displaystyle 0.0425138097805085\]

It might now be interesting to evaluate this radial function at a range of distances and plot them. This function does not support taking multiple radii, so you have two options below.

1. Iterate through a list or array of radii and evaluate this function one radius at a time.

2. Convert the R_nl() function to a function that can accept an array using the lambdify() method.

Both approaches are demonstrated below.

# first approach - iterate through iterable
radii = np.linspace(0, 30, 200)
R_eval = [R_nl(3, 1, r, Z=1) for r in radii]

plt.plot(radii, R_eval)
plt.hlines(0, 0, 30, colors='r', linestyles='dashed')
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Evaluated Radiable Wave');

# second approach - lambdify
r = sympy.symbols('r') # create SymPy symbol

# create a numpy compatible function using lambdify
R_3p = sympy.lambdify(r, R_nl(3, 1, r, Z=1), modules='numpy')
radii = np.linspace(0, 30, 200)

plt.plot(radii, R_3p(radii))
plt.hlines(0, 0, 30, colors='r', linestyles='dashed')
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Evaluated Radiable Wave');

The electron probability density can be found by calculating \(R^2\) where \(R\) is the radial wavefunction, and the radial probability is \(R^2r^2\) where \(r\) is the distance from the nucleus.

plt.plot(radii, R_3p(radii)**2 * radii**2)
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Radial Probability, $R^2r^2$');

The reason we multiply the probability density by the square of the radial wavefunction, \(r^2\), is to account for the greater surface area of a sphere (\(A_{sphere} = 4\pi r^2\)) the larger the radius. We are effectively carrying out the calculation depicted below. We divide the sphere surface area by 4\(\pi\) to normalize the integration making the probability over all space total to one.

Show code cell source Hide code cell source

fig = plt.figure(figsize=(10, 4))

ax1 = fig.add_subplot(1, 3, 1)
ax1.plot(radii, R_3p(radii)**2, color='C0')
ax1.set_xlabel('Radius (Bohrs)')
ax1.set_title('1s Probability Density')

ax2 = fig.add_subplot(1, 3, 2)
ax2.plot(radii, 4 * sympy.pi * radii**2, color='C1')
ax2.set_xlabel('Radius (Bohrs)')
ax2.set_ylabel('Surface Area / 4$\\pi$ (Bohrs$^2$)')
ax2.set_title('Sphere Surface Area Over 4$\\pi$')

ax3 = fig.add_subplot(1, 3, 3)
ax3.plot(radii, radii**2 * R_3p(radii)**2, color='C2')
ax3.set_xlabel('Radius (Bohrs)')
ax3.set_ylabel('Radial Probability')
ax3.set_title('1s Radial Probability')

fig.tight_layout()
fig.subplots_adjust(wspace=0.5)
fig.text(0.31, 0.5, 'X', fontdict={'size': 'xx-large'})
fig.text(0.66, 0.5, '=', fontdict={'size': 'xx-large'});

One of the uses of these radial plots is to compare the radial probability of multiple different orbitals on the same axes like below for the fourth row of the periodic table. This can be used, for example, to discuss the valence electron configurations of Cr and Cu.

r = sympy.symbols('r') # create SymPy symbol

# create a numpy compatible function using lambdify
R_3s = sympy.lambdify(r, R_nl(4, 0, r, Z=1), modules='numpy')
R_3p = sympy.lambdify(r, R_nl(4, 1, r, Z=1), modules='numpy')
R_3d = sympy.lambdify(r, R_nl(3, 2, r, Z=1), modules='numpy')

radii = np.linspace(0, 45, 200)

plt.plot(radii, R_3s(radii)**2 * radii**2, label='4s')
plt.plot(radii, R_3p(radii)**2 * radii**2, label='4p')
plt.plot(radii, R_3d(radii)**2 * radii**2, label='3d')
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Radial Probability ($R^2r^2$)');
#plt.ylim(0, 0.2)
plt.legend();

The probability \(R^2r^2\) can be integrated using the sympy.integrate() function which accepts the function or mathematical expression to be integrated and a tuple that contains the variable, the min, and the max values.

sympy.integrate(f(x), (x, min, max))

For example, we can integrate the R_nl() for the 2s orbital from 0 to 3.0 Bohrs like below.

\[\displaystyle 0.473330547984585\]

Let’s test that the radial probability is normalized by integrating from zero to infinity.

Note

The sympy.oo is the SymPy variable for infinity.

sympy.integrate(R_nl(2, 0, r, Z=1)**2 * r**2, (r, 0, sympy.oo)).evalf()

\[\displaystyle 1.0\]

Angular Wavefunctions

The other component of \(\Psi\) is the angular wavefunctions which provides directional information about an orbital. The angular equations can be coded by hand or we can also use the Y_lm() or Z_lm() spherical harmonics wavefunctions from sympy.physics.hydrogen to assist us. The difference between these two functions is that Y_lm() may return a complex expression whereas Z_lm() will return the real-valued angular wavefunction. Because our goal is to visualize the wavefunctions, we will restrict ourselves to the latter here. The angular wavefunction provides information in all directions, so we will plot this information in 3D.

Warning

The plot of angular wavefunction does not include the radial information so does not fully describe the shape of atomic orbitals. Do not interpret the angular plots below as the actual shape of atomic orbitals even though they resemble them.

There are multiple conventions for spherical coordinates. We will use the SciPy/SymPy convention of using theta (\(\theta\)) for the azmuthal (i.e., direction on xy-plane) and phi (\(\phi\)) as the polar angle (i.e., angle from the positive z-axis) for plotting the angular wavefunctions. Below, we plot the \(d_{z^2}\) orbital by coding the angular wavefunction expression by hand.

Tip

See section 3.6.3 for guidance on plotting surfaces in 3D.

# generate mesh grid of theta and phi values
theta, phi = np.meshgrid(np.linspace(0, np.pi, 51),
                         np.linspace(0, 2 * np.pi, 101))

# convert angles to xyz values of a sphere, r = 1
x = np.sin(theta) * np.sin(phi)
y = np.sin(theta) * np.cos(phi)
z = np.cos(theta)

# multiply xyz values by angular wavefunction
dz2 = np.sqrt((5 / 16) * np.pi) * (3 * np.cos(theta)**2 - 1)
X, Y, Z = x * dz2, y * dz2, z * dz2

fig = plt.figure(figsize = (10,6))
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')
ax.set_aspect('equal') # sets aspect ratio to equal

Alternatively, we can use the Z_lm() function from sympy.physics.hydrogen to generate the angular wavefunction based on the angular and magnetic quantum numbers.

Z_lm(l, m, phi, theta)

SymPy functions cannot calculate wavefunctions for an array of angles like NumPy functions can, but fortunately SymPy functions can be converted to NumPy function using the lambdify() method. Just provide the lambdify() method with a collection of argument variables for the wavefunction as SymPy symbols, the wavefunction, and modules='numpy', and it returns a new function.

# from sympy.physics import Z_lm

theta, phi = np.meshgrid(np.linspace(0, np.pi, 51),
                         np.linspace(0, 2*np.pi, 101))

x = np.sin(theta) * np.sin(phi)
y = np.sin(theta) * np.cos(phi)
z = np.cos(theta)

# create a numpy function
p, t = sympy.symbols('p t')
f = sympy.lambdify((p, t), Z_lm(2, 0, p, t), modules='numpy')

# multiply xyz values by wave angular wavefunction
f_pt = f(phi, theta)
X, Y, Z = x * f_pt, y * f_pt, z * f_pt

fig = plt.figure(figsize = (10, 6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_aspect('equal') # sets aspect ratio to equal
ax.set_axis_off() # turns off axes and background

We can also visualize the angular component of wavefunctions in 2D using a polar plot, but we can only visualize one angle at a time. Below we will visualize theta and leave phi fixed. Because we are only visualizing in 2D and not sweeping around the phi angles, we need to make theta go from 0 \(\rightarrow\) 2\(\pi\).

l, m = 2, 0
azmuth, polar = sympy.symbols('azmuth polar')
f = sympy.lambdify((polar, azmuth), Z_lm(l, m, polar, azmuth), modules='numpy')

th = np.linspace(0, 2 * np.pi, 200)
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
ax.plot(th, np.abs(f(0, th)))
# orient 0 degrees to up/north
ax.set_theta_zero_location('N');

l, m = 2, 1
azmuth, polar = sympy.symbols('azmuth polar')
f = sympy.lambdify((polar, azmuth), Z_lm(l, m, polar, azmuth), modules='numpy')

th = np.linspace(0, 2 * np.pi, 200)
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
ax.plot(th, np.abs(f(0, th)))
# orient 0 degrees to up/north
ax.set_theta_zero_location('N');

l, m = 2, 2
azmuth, polar = sympy.symbols('azmuth polar')
f = sympy.lambdify((polar, azmuth), Z_lm(l, m, polar, azmuth), modules='numpy')

th = np.linspace(0, 2 * np.pi, 200)
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
ax.plot(th, np.abs(f(0, th)))
# orient 0 degrees to up/north
ax.set_theta_zero_location('N');

The last orbital image is a d-orbital viewed from the side.

Complete Wavefunction

Now we will visualize both angular and radial components together (\(\Psi\)) which is again the product of the radial, \(R(n, l)\) and angular, \(Y(l, m)\) wavefunctions.

\[ \Psi(n, l, m) = R(n, l)Y(l, m) \]

To obtain the entire wavefunction, \(\Psi\), we can either multiply the radial and angular wavefunctions from the previous sections or use the SymPy Psi_nlm() function which makes this task a little more convenient. Orbitals have no edge, so there are multiple ways of representing orbitals including contour plots, isosurfaces, 90% surface plots, scatter plots, and translucent 3D plots. The scatter and contour plot methods are demonstrated below. We will need the probability density, P, of the atomic orbital which is proportional to the product of a wavefunction, \(\Psi\), and its complex conjugate, \(\Psi\)* or the square of the absolute value of a wavefunction.

\[ P = \Psi^*\Psi = |\Psi|^2 \]

First, let’s take a look at the Psi_nlm() function which operates similarly to the other SymPy wavefunctions above. Below, we integrate it over all space returning 1 which tells us that this function is normalized when we include \(r^2sin(\theta)\).

Note

The \(|\Psi|^2\) approach is favored below, but if you want to use \(\Psi^*\Psi\), you can wrap your wavefunction in sympy.conjugate().

azmuth, polar, r = sympy.symbols('azmuth polar r')
wf = Psi_nlm(3, 1, 0, r, azmuth, polar)

# integrate normalized wavefunction over all area
sympy.integrate(wf**2 * r**2 * sympy.sin(polar), 
                (r, 0, sympy.oo), 
                (azmuth, 0, 2 * sympy.pi), 
                (polar, 0, sympy.pi))

\[\displaystyle 1\]

Now let’s visualize an orbital using a scatter plot. We will use a strategy previously reported in J. Chem. Educ., 1990, 67, 42-44 which includes the following steps.

1. Use a random number generator to produce a series of \(r\), \(\theta\), and \(\phi\) values or just \(r\) and \(\theta\) values depending upon dimensions

2. Use the values above to calculate the xyz or yz values

3. Use the above radius and angles to calculate probabilities using the wavefunction

4. Normalize the probabilities by dividing by the maximum probability value across all the data points

5. If each normalized probability is above a random value from 0 \(\rightarrow\) 1, it gets included in the scatter plot

Tip

If plotting a very large number of data points, consider using plt.plot() instead of plt.scatter() because the latter is slower and uses more memory due to its ability to individualize each marker in the plot.

# 2p orbital - 3D simulation

# create wavefunction as python function
r, azmuth, polar = sympy.symbols('r azmuth polar')
wf_sym = Psi_nlm(2, 1, 0, r, azmuth, polar)
wf = sympy.lambdify((r, azmuth, polar), wf_sym, modules='numpy')

# generate random coordinates
rng = np.random.default_rng(seed=21)
n_points = 100000
r = 15 * rng.random(size=(n_points))
polar = np.pi * rng.random(size=(n_points))
azmuth = 2 * np.pi * rng.random(size=(n_points))

x = r * np.sin(polar) * np.sin(azmuth)
y = r * np.sin(polar) * np.cos(azmuth)
z = r * np.cos(polar)

# normalize and create mask
prob_dens = np.abs(wf(r, azmuth, polar))**2
norm_prob = prob_dens / prob_dens.max()
mask = norm_prob > rng.random(n_points)

#plt.plot(y[mask], z[mask], ',');
fig = plt.figure(figsize = (6,6))
ax = fig.add_subplot(1, 1, 1)
ax.scatter(y[mask], z[mask], s=0.1)
ax.set_xlabel('Distance from Nucleus (Bohrs)')
ax.set_ylabel('Distance from Nucleus (Bohrs)');

The plot above makes the shape of the orbital look like the orbital lobes are almost conical, which is not what we typically see in accurate orbital shapes. This is an effect of there being more data points visualized along the vertical axis due to the orbital being thicker there. If we instead reduce the simulation to 2D (i.e., only yz plane), like below, the orbital lobes appear rounder because we are visualizing a slice through the middle of the orbital.

# 2p orbital - 2D simulation

# create wavefunction as python function
r, polar = sympy.symbols('r polar')
wf_sym = Psi_nlm(2, 1, 0, r, 0, polar)
wf = sympy.lambdify((r, polar), wf_sym, modules='numpy')

# generate random coordinates
rng = np.random.default_rng(seed=21)
n_points = 100000
r = 15 * rng.random(size=(n_points))
polar = 2 * np.pi * rng.random(size=(n_points))

x = r * np.sin(polar) * np.sin(0)
y = r * np.sin(polar) * np.cos(0)
z = r * np.cos(polar)

# normalize and create mask
prob_dens = np.abs(wf(r, polar))**2
norm_prob = prob_dens / prob_dens.max()
mask = norm_prob > rng.random(n_points)

fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(1, 1, 1)
ax.scatter(y[mask], z[mask], s=0.1)
ax.set_xlabel('Distance from Nucleus (Bohrs)')
ax.set_ylabel('Distance from Nucleus (Bohrs)');

We can visualize larger orbitals to see more nodes such as in the 3p and 3s orbitals below. We can also color the points based on the sign of the wavefunction before calculating the probability. In the examples below, the color only represents the sign of the wavefunction and not the magnitude of the value.

# 3p orbital

# create wavefunction as python function
r, polar = sympy.symbols('r, polar')
wf_sym = Psi_nlm(3, 1, 0, r, 0, polar)
wf = sympy.lambdify((r, polar), wf_sym, modules='numpy')

# generate random coordinates
rng = np.random.default_rng(seed=21)
n_points = 500000
r = 30 * rng.random(size=(n_points))
polar = 2 * np.pi * rng.random(size=(n_points))

x = r * np.sin(polar) * np.sin(0)
y = r * np.sin(polar) * np.cos(0)
z = r * np.cos(polar)

# normalize and create mask
prob_dens = np.abs(wf(r, polar))**2
norm_prob = prob_dens / prob_dens.max()
mask = norm_prob > rng.random(n_points)

#plt.plot(x[mask], y[mask], 
fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(1, 1, 1)
is_pos = wf(r, polar)[mask] > 0 # test if wavefunc is positive
ax.scatter(y[mask], z[mask], s=0.5, c=is_pos, cmap='coolwarm')
ax.set_xlabel('Distance from Nucleus (Bohrs)')
ax.set_ylabel('Distance from Nucleus (Bohrs)');

# 3s orbital

# create wavefunction as python function
r, polar = sympy.symbols('r, polar')
wf_sym = Psi_nlm(3, 0, 0, r, 0, polar)
wf = sympy.lambdify((r, polar), wf_sym, modules='numpy')

# generate random coordinates
rng = np.random.default_rng(seed=21)
n_points = 1000000
r = 30 * rng.random(size=(n_points))
polar = 2 * np.pi * rng.random(size=(n_points))

x = r * np.sin(polar) * np.sin(0)
y = r * np.sin(polar) * np.cos(0)
z = r * np.cos(polar)

# normalize and create mask
prob_dens = np.abs(wf(r, polar))**2
norm_prob = prob_dens / prob_dens.max()
mask = norm_prob > rng.random(n_points)

fig = plt.figure(figsize = (6, 6))
ax = fig.add_subplot(1, 1, 1)
is_pos = wf(r, polar)[mask] > 0 # test if wavefunc is positive
ax.scatter(y[mask], z[mask], s=0.5, c=is_pos, cmap='coolwarm')
ax.set_xlabel('Distance from Nucleus (Bohrs)')
ax.set_ylabel('Distance from Nucleus (Bohrs)');

A second way to visualize orbitals is through a contour plot. Here we calculate the probability in a mesh of locations and provide the plt.countour() function with the locations and probabilities.

Y, Z = np.meshgrid(np.linspace(-20, 20, 200),
                   np.linspace(-20, 20, 200))

# create wavefunction as python function
r, polar = sympy.symbols('r polar')
wf = Psi_nlm(3, 1, 0, r, 0, polar)
f = sympy.lambdify((r, polar), wf * sympy.conjugate(wf), modules='numpy')
polar = np.arctan(Z / Y)
r = np.sqrt(Y**2 + Z**2)

# calculate probability
prob = np.abs(f(r, polar))**2

plt.contour(Z, Y, prob, levels=[1e-9, 3e-9, 5e-9, 1e-8, 5e-8, 1e-7, 3e-7, 9e-7])
plt.colorbar()
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Distance from Nucleus (Bohrs)');

Y, Z = np.meshgrid(np.linspace(-20, 20, 200),
                   np.linspace(-20, 20, 200))

# create wavefunction as python function
r, polar = sympy.symbols('r polar')
wf = Psi_nlm(3, 2, 0, r, 0, polar)
f = sympy.lambdify((r, polar), wf * sympy.conjugate(wf), modules='numpy')
polar = np.arctan(Z / Y)
r = np.sqrt(Y**2 + Z**2)

# calculate probability
prob = np.abs(f(r, polar))**2

plt.contour(Z, Y, prob, levels=[1e-9, 3e-9, 5e-9, 1e-8, 5e-8, 1e-7, 3e-7, 5e-7])
plt.colorbar()
plt.xlabel('Distance from Nucleus (Bohrs)')
plt.ylabel('Distance from Nucleus (Bohrs)');