Atomic Orbitals

AO.py module

The AO module provides the exact, analytical hydrogenic atomic orbitals (radial parts and real spherical harmonics, as implemented in SymPy’s sympy.physics.hydrogen and sympy.Znm) together with a set of plotting functions to visualize them: radial curves, 2D angular polar diagrams, density-weighted electron clouds, and interactive 3D nodal surfaces.

Physical background

For the hydrogen atom (or any hydrogenic ion of nuclear charge \(Z\)), the time-independent Schrödinger equation

\[\hat{H}(r,\theta,\varphi)\,\psi(r,\theta,\varphi) = E\,\psi(r,\theta,\varphi)\]

is separable in spherical coordinates. Each stationary state is a product wavefunction

\[\psi_{n\ell m}(r,\theta,\varphi) = R_{n\ell}(r)\, Y_\ell^m(\theta,\varphi)\]

where:

  • \(R_{n\ell}(r)\) is the radial part, labeled by the principal quantum number \(n\) and the azimuthal quantum number \(\ell\) (\(0 \le \ell \le n-1\)). It has \(n-\ell-1\) radial nodes (spherical surfaces where \(R_{n\ell}=0\)).

  • \(Y_\ell^m(\theta,\varphi)\) is the angular part, a spherical harmonic labeled by \(\ell\) and the magnetic quantum number \(m\) (\(-\ell \le m \le \ell\)). The module uses the real combinations \(Z_\ell^m\) (SymPy’s Znm), which are the familiar real-valued \(p_x\), \(p_y\), \(p_z\), \(d_{z^2}\), \(d_{x^2-y^2}\), … orbitals used in chemistry. \(Y_\ell^m\) has angular nodes: cones (\(\theta\) = const) and/or planes (\(\varphi\) = const) where the function changes sign.

The electron probability density is \(|\psi_{n\ell m}|^2\), and the radial probability density — the probability of finding the electron in a thin spherical shell at distance \(r\), integrated over all angles — is

\[P_{n\ell}(r) = r^2 R_{n\ell}(r)^2\]

This is the quantity plotted by plot_radial() when density=True.

All lengths are expressed in units of the Bohr radius \(a_0\) (scaled by \(1/Z\) for hydrogenic ions with nuclear charge \(Z \ne 1\)), consistent with SymPy’s convention for R_nl and Psi_nlm.

Quick start

from pyphyschemtools import AO

# Radial parts of the 1s, 2s and 2p orbitals
AO.plot_radial([(1, 0), (2, 0), (2, 1)])

# Radial probability density r^2 R_nl(r)^2 instead of R_nl(r)
AO.plot_radial([(1, 0), (2, 0), (2, 1)], density=True)

# 2D angular polar diagram of a 3d_(x^2-y^2)-like orbital (l=2, m=2), xy plane
AO.plot_polar(l=2, m=2, plane="xy")

# Density-weighted electron cloud of the 2p_z orbital, xz projection
AO.plot_electron_cloud(n=2, l=1, m=0, plane="xz")

# Interactive 3D nodal surfaces of the 3d_z2 orbital, with the electron
# cloud overlaid
AO.plot_orbital_3D(n=3, l=2, m=0, show_cloud=True)

Every plotting function accepts show_formula=True to display (via IPython) the symbolic SymPy expression of the radial and/or angular part being plotted — useful for checking a result against a textbook formula, or for teaching.

Set save_as=... on any function to save the figure to disk (a static image for plot_radial, plot_polar and plot_electron_cloud; a static image or an interactive .html for plot_orbital_3D, selected by the file extension).

Orbital naming convention

Orbitals are identified throughout the module by (n, l) or (n, l, m) tuples, and labeled on plots with the usual spectroscopic notation (1s, 2s, 2p, 3d, …) via the l = 0, 1, 2, 3, 4, 5 -> s, p, d, f, g, h correspondence.

API reference

Atomic orbitals of the hydrogen atom and hydrogenic ions: radial parts, real spherical harmonics, and their visualization (radial plots, 2D polar cross-sections, electron cloud projections, and interactive 3D nodal surfaces).

Public functions: - plot_radial : radial part R_nl(r) or radial probability density - plot_polar : 2D polar diagram of the angular part, in a chosen plane - plot_electron_cloud : electron cloud of an orbital (density-weighted point sampling), projected onto a chosen plane - plot_orbital_3D : interactive 3D nodal surfaces of an orbital (radial nodes as spheres, angular nodes as cones and planes), with the electron cloud as an optional overlay

pyphyschemtools.AO.plot_electron_cloud(n, l, m, plane='xz', Z=1, n_points=20000, point_size=3, alpha=0.25, color_pos='tab:blue', color_neg='tab:red', size=(6, 6), show_formula=False, save_as=None)[source]

Draw the electron cloud of a hydrogenic orbital psi_nlm, projected onto a plane, as a density-weighted point cloud: points are sampled with probability proportional to |psi_nlm|^2 (rejection sampling in 3D, see _sample_orbital_cloud), then projected onto the chosen plane. Both radial nodes (R_nl(r) = 0) and angular nodes (Z_l^m = 0) appear as empty regions, since no points are sampled there – no separate geometric reconstruction is needed.

input:

  • n, l, m = quantum numbers of the orbital (0 <= l <= n-1, -l <= m <= l)

  • plane = projection plane: “xy”, “xz” or “yz” (default: “xz”)

  • Z = nuclear charge (default: 1 for hydrogen)

  • n_points = number of sampled points

  • point_size, alpha = marker size and transparency of each point

  • color_pos, color_neg = colors for points where psi_nlm is positive / negative

  • size = (width, height) in inches for the figure

  • show_formula = if True, display (via IPython) the symbolic expressions of R_nl(r) and Z_l^m(theta, phi)

  • save_as = file path to save the figure to (passed to pyphyschemtools.core’s save_fig); not saved if None (default)

output: None (the figure is drawn and shown, and saved if save_as is given)

pyphyschemtools.AO.plot_orbital_3D(n, l, m, Z=1, prob_fraction=0.7, show_nodes=True, show_cloud=False, n_points=20000, point_size=2, color_pos='blue', color_neg='red', grid_points=60, width=900, height=900, show_formula=False, save_as=None)[source]

Draw the 3D shape of a hydrogenic orbital psi_nlm as isosurfaces: the surface(s) where |psi_nlm| = iso_fraction * max(|psi_nlm|), colored by the sign of psi_nlm. For orbitals with radial nodes (n-l-1 > 0), this naturally shows nested lobes – a smaller inner surface inside a larger outer one (e.g. 2s, 3p, 4d, …), with no extra work needed: the isosurface algorithm finds every disconnected region on its own.

Optionally, the exact nodal surfaces (show_nodes=True) and/or the density-weighted electron cloud (show_cloud=True) can be overlaid.

input:

  • n, l, m = quantum numbers of the orbital (0 <= l <= n-1, -l <= m <= l)

  • Z = nuclear charge (default: 1 for hydrogen)

  • prob_fraction = fraction of the total electron probability enclosed by the lobe surfaces (default: 0.7 – higher shows fatter lobes)

  • show_nodes = if True, also draw the exact nodal surfaces (radial spheres, angular cones/planes); default: True

  • show_cloud = if True, also draw the sampled electron cloud; default: False

  • n_points = number of sampled points for the cloud (only used if show_cloud=True)

  • point_size = marker size for the cloud points

  • color_pos, color_neg = colors for the lobes/cloud points where psi_nlm is positive / negative

  • grid_points = resolution (per axis) of the 3D grid used for the lobe isosurfaces

  • width, height = size of the figure, in pixels

  • show_formula = if True, display (via IPython) the symbolic expressions of R_nl(r) and Z_l^m(theta, phi)

  • save_as = file path to save the figure to, as a static image (requires kaleido) or an interactive .html; format taken from the extension. Not saved if None (default)

output: None (the figure is shown interactively, and saved if save_as is given)

pyphyschemtools.AO.plot_polar(l, m, plane='xy', density=False, show_formula=False, n_points=400, size=(6, 6), save_as=None)[source]

Draw the 2D polar diagram of the angular part of an atomic orbital (the real spherical harmonic Z_lm, independent of n) in a plane containing the origin.

A sweep angle gamma (0 to 2*pi) runs once around the chosen plane; at each gamma the corresponding direction (theta, phi) is computed, and the radius plotted is |Z_lm(theta, phi)| (or its square if density=True). The two signs of Z_lm are filled in different colors, so the lobes of p, d, … orbitals and their phase change are visible; the radius goes to zero at a node (sign change).

input:

  • l = orbital angular momentum quantum number (l >= 0)

  • m = magnetic quantum number (-l <= m <= l)

  • plane = projection plane: “xy”, “xz” or “yz” (default: “xy”)

  • density = if False (default), plot |Z_lm|; if True, plot the angular probability density Z_lm^2

  • show_formula = if True, display (via IPython) the symbolic expression of the quantity plotted (and its squared/density form if density=True); default: False

  • n_points = number of points used to sample the sweep angle

  • size = (width, height) in inches for the figure

  • save_as = file path to save the figure to (passed to pyphyschemtools.core’s save_fig); not saved if None (default)

output: None (the figure is drawn and shown, and saved if save_as is given)

pyphyschemtools.AO.plot_radial(orbitals, Z=1, density=False, show_formula=False, r_max=None, n_points=500, size=(9, 5), save_as=None)[source]

Plot the radial part R_nl(r) (default) or the radial probability density r^2 R_nl(r)^2 of one or several hydrogenic atomic orbitals, as a function of r (in units of the Bohr radius a0, scaled by 1/Z).

input:

  • orbitals = list of (n, l) tuples, one curve per orbital, e.g. [(1, 0), (2, 0), (2, 1)] for 1s, 2s, 2p

  • Z = nuclear charge (default: 1 for hydrogen)

  • density = if False (default), plot R_nl(r); if True, plot the radial probability density r^2 R_nl(r)^2 (probability of finding the electron in a thin spherical shell at r)

  • show_formula = if True, display (via IPython) the symbolic expression of the quantity plotted (and its squared/density form if density=True); default: False

  • r_max = largest r shown, in units of a0/Z; if None (default), set automatically from the largest n among orbitals

  • n_points = number of points used to sample r

  • size = (width, height) in inches for the figure

  • save_as = file path to save the figure to (passed to pyphyschemtools.core’s save_fig); not saved if None (default)

output: None (the figure is drawn and shown, and saved if save_as is given)