Hydrogen atomic orbitals

Symbolic generator and 3D viewer of hydrogen-like atomic orbitals for any valid set of quantum numbers n, l and m, with a configurable atomic number Z. The tool does three things:

• Derives the wavefunction Psi_{n,l,m} symbolically with sympy, in the complex basis or in the real one (the familiar px, py, dxy combinations).
• Prints the energy, the radial part, the angular part and the full wavefunction to the terminal, along with their LaTeX source, ready to paste into a paper or a lecture note.
• Renders the orbital in 3D with pyvista using four different techniques, including a solid surface that carries the phase of a complex orbital as color.

_{--n 5 --l 4 --m 4 --plot-method isosurface-phase --limit 42 --n-points 150. The surface is a level set of |Psi|^2 and the color is arg(Psi), which completes four full turns around the z axis, one per unit of m.}

Under every figure are the options that produced it. The complete commands, render size included, live in scripts/make_figures.sh, which regenerates the whole set in one go.

Getting started

The virtual environment is not part of the repository, so the first step is to build it:

git clone https://github.com/jofavalle/atomic-orbitals-of-hydrogen-simulation.git
cd atomic-orbitals-of-hydrogen-simulation
python3 -m venv .venv
.venv/bin/pip install -r requirements.txt

Every command below is run from the repository root, so that python -m hydrogen_orbitals can find the package. Either call the interpreter inside the environment directly, as shown here, or activate the environment first with source .venv/bin/activate (.venv\Scripts\activate on Windows) and then use a plain python.

To confirm the installation, run the test suite:

.venv/bin/pytest

Quick start

# Real 2px orbital, selected by its short alias
.venv/bin/python -m hydrogen_orbitals --name 2px

# The same state, giving n, l and m explicitly
.venv/bin/python -m hydrogen_orbitals --n 2 --l 1 --m 1 --real

# Complex orbital (n=3, l=2, m=+2) as a point cloud with an animated phase
.venv/bin/python -m hydrogen_orbitals --n 3 --l 2 --m 2

# The same complex orbital as a solid surface with the phase painted on it
.venv/bin/python -m hydrogen_orbitals --n 3 --l 2 --m 2 --plot-method isosurface-phase

# Print the wavefunction and its LaTeX source without opening any window
.venv/bin/python -m hydrogen_orbitals --name 3dz2 --print-wavefunction --latex --no-plot

# Write the scene to a file instead of opening a window
.venv/bin/python -m hydrogen_orbitals --n 4 --l 3 --m 3 --plot-method isosurface-phase \
    --screenshot orbital.png --window-size 1000 750

Quantum numbers

• n (principal): integer >= 1. Sets the energy E_n = -Z^2/(2n^2) in Hartree atomic units and the typical size of the orbital, which grows roughly as n^2 Bohr radii.
• l (azimuthal, angular momentum): integer with 0 <= l < n. Sets the angular shape: l=0 is s, l=1 is p, l=2 is d, l=3 is f, l=4 is g, and so on.
• m (magnetic): integer with -l <= m <= l. In the complex basis it indexes the exp(i m phi) factor; in the real basis its sign picks the cosine-like combination (m > 0) or the sine-like one (m < 0).
• Z (atomic number, optional, --Z): 1 for hydrogen, 2 for He+, and so on. It rescales the radial part, contracting the orbital as Z grows.

Selecting an orbital

There are two equivalent ways to pick a state.

1. Short alias (--name), available for the well known real orbitals:

l	alias (m>0, cosine)	alias (m=0)	alias (m<0, sine)
0	-	s	-
1	px (m=+1)	pz (m=0)	py (m=-1)
2	dxz (m=+1), dx2-y2 (m=+2)	dz2 (m=0)	dyz (m=-1), dxy (m=-2)
3	fxz2 (m=+1), fzx2-y2 (m=+2), fx3-3y2 (m=+3)	fz3 (m=0)	fyz2 (m=-1), fxyz (m=-2), fy3-3x2 (m=-3)

Prefix the alias with n, as in --name 4dxy or --name 5fz3. There are no aliases for l >= 4 (g, h, ...); use --n/--l/--m for those.

2. Explicit quantum numbers (--n --l --m), valid for any l, with an optional --real to ask for the real combination instead of the complex one.

The choice of basis is not a matter of taste, and it changes the picture. Both figures below are n = 4, l = 3, |m| = 3 at the same iso level:

Real basis: --name 4fx3-3y2	Complex basis: --n 4 --l 3 --m 3

The real combination is proportional to cos(3 phi), so its density vanishes on six half planes and the orbital breaks into six lobes, colored by the sign of Psi. The complex state carries exp(3 i phi), whose modulus is 1: the density does not depend on the azimuth at all, which is why the same state is a smooth torus, and the information that the real basis spent on shape is here carried by the color.

Visualization methods

All four methods take their shape from the probability density |Psi|^2. What changes is how that density is turned into geometry and what the color means.

point-cloud (default): draws a Monte Carlo point cloud built by rejection sampling, so the density of points is proportional to |Psi|^2. Complex orbitals are colored by their phase arg(Psi) with a cyclic HSV map, and that phase is animated. It is the method that shows the interior of the orbital best, since it is transparent by construction.

_{--n 6 --l 3 --m 2 --limit 62 --n-candidates 20000000. This state has two radial nodes, and the cloud is the method that lets the inner shells and the hollow between them be seen through the outer one.}

isosurface-phase: draws the solid surface |Psi|^2 = const and paints the phase of Psi on it, so a complex orbital can be read as a solid body without giving up its phase. Around the z axis the color completes exactly |m| full turns, which is the phase factor exp(i m phi) made visible. The phase is animated just like in point-cloud. For a real wavefunction the phase collapses to a sign and the surface is colored by lobe, so the method is useful in both bases.

_{--n 6 --l 4 --m 2 --limit 60 --n-points 150. The l - |m| = 2 polar nodes cut the surface into the three rings stacked along z, and the color turns twice around the axis, once per unit of m.}

isosurface-density: the same surface in a single flat color. Useful for the textbook shape of an orbital, with no phase information at all.

_{--name 4fxyz --plot-method isosurface-density --limit 30 --n-points 140.}

isosurface-signed: two surfaces, at +c and -c of the real Psi, colored cyan and magenta by lobe sign. It only makes sense for real orbitals; asked for a complex one it warns on the console and falls back to isosurface-density.

_{--name 4fxyz --plot-method isosurface-signed --limit 30 --n-points 140. The same orbital as the figure above, so that the method is the only difference between the two: the eight lobes sit on the corners of a cube and alternate in sign, which is what the flat gray surface cannot tell.}

Notes on the phase surface

Two details of isosurface-phase are worth knowing, because they explain its options and its cost.

The phase is re-evaluated, never interpolated. arg(Psi) is a cyclic quantity: two neighbouring points can sit at +3.13 and -3.13 radians, which is a gap of 0.02 in phase but of almost 2 pi in value. Interpolating the angle across that branch cut tears the color map with a visible seam. The tool instead evaluates Psi exactly at the surface vertices, which is both exact and cheap, since a surface holds of the order of 1e4 vertices against the n_points^3 evaluations already spent on the grid.

The lighting is deliberately flat. With standard shading, the faces turned away from the light go dark and their hue becomes unreadable, which defeats the purpose of the method. --surface-ambient (default 0.45) sets the ambient term, which does not depend on the surface normal and therefore puts a floor under the brightness of every face. Raise it if the phase is still hard to read, lower it for a stronger sense of relief.

Use --surface-opacity below 1.0 to see through the outer shell, which matters for orbitals with radial nodes (n - l - 1 > 0), where the surface has several nested components and the outermost one hides the rest.

Opaque (default)	--surface-opacity 0.35

_{--n 6 --l 2 --m 2 --plot-method isosurface-phase --iso-frac 0.02 --limit 70 --n-points 150, with and without the opacity. This state has n - l - 1 = 3 radial nodes, so its surface is four nested shells and the outermost one is the only thing an opaque render shows.}

Finally, the color is interpolated between vertices, so a coarse grid facets the color gradient when |m| is large. If the bands look blocky, raise --n-points.

CLI reference

Run .venv/bin/python -m hydrogen_orbitals --help for the full list. The most relevant options:

Option	Description
--name ALIAS	Select by short alias (implies a real orbital).
--n --l --m	Select explicitly by quantum numbers.
--real	Force the real basis when using --n/--l/--m.
--Z	Atomic number (default 1, hydrogen).
--print-wavefunction	Print energy, radial part, angular part and full Psi to the terminal.
--latex	Print the LaTeX source of Psi, the radial part and the angular part.
--plot-method	point-cloud (default), isosurface-phase, isosurface-density or isosurface-signed.
--no-plot	Do not open any window, useful together with --print-wavefunction or --latex.
--animate / --no-animate	Force or disable the phase animation (default: on for complex orbitals with m != 0).
--spin-speed	Base angular speed of the animation (rad/step; the actual speed is m * spin-speed).
--limit	Half-size of the visualization box, in Bohr radii (defaults to a value scaling with n^2).
--n-points	Resolution of the 3D grid for the isosurface-* methods.
--n-candidates	Number of Monte Carlo candidates for point-cloud (defaults to a value scaling with n^2).
--iso-frac	Fraction of the maximum density used as the isosurface level.
--surface-opacity	Opacity of the isosurface-phase surface (default 1.0).
--surface-ambient	Ambient light of the isosurface-phase surface (default 0.45).
--grid-step	Spacing of the reference grid drawn on the planes through the origin (default: limit/4).
--parallel-projection	Drop the perspective, so equal distances along an axis are equal on screen.
--axis-labels	Write the tick values on the axes, where the reference planes cross them. The box itself carries no numbers.
--screenshot PATH	Render off screen and write the scene to PATH, without opening a window.
--gif PATH	Write one full turn of the phase animation to PATH as a GIF.
--gif-frames, --gif-fps	Frames in that turn (default 40) and playback rate (default 20).
--window-size W H	Size of the render in pixels, for the window and for a capture alike.

The scene draws its reference grid on the planes x = 0, y = 0 and z = 0, not only on the outer faces of the box. That is what makes an angular node readable: the node of 3d (m = 1), for instance, lies on z = 0, and the grid plane runs through the gap between its two lobes. Under the default perspective camera a plane still shortens with depth, so add --parallel-projection when you want to read a coordinate off the axes without that distortion, at the cost of the depth cue.

Perspective camera (default)	--parallel-projection

_{--n 3 --l 2 --m 1 --plot-method isosurface-phase --no-animate --axis-labels --limit 20 --n-points 130, with and without the flag. The angular node of this state lies on z = 0, exactly where a reference plane runs, and it is the gap between the two lobes.}

The box states the name of each axis and carries no numbers. A number written on an edge names its plane from a far corner, and since a plane projects to a band on screen rather than to a line, the two only meet in a view looking straight down an axis: the label reads as displaced everywhere else. --axis-labels puts the values where that cannot happen, on the points where the planes cross the axes, at the price of a busier picture, with a number sitting in front of the orbital at every crossing.

What the phase animation does and does not mean

For a stationary eigenstate, |Psi(r,theta,phi,t)|^2 = |Psi(r,theta,phi)|^2 does not depend on time: a pure orbital does not spin in the sense of its probability density changing. What the animation shows is the direction and relative magnitude of the phase gradient of exp(i m phi), which is directly related to the azimuthal probability current and to the sign of the angular momentum along z:

• The angular speed is m * spin-speed. The direction of rotation follows the sign of m (counterclockwise seen from +z when m > 0, consistent with L_z = m*hbar > 0), and the relative speed between different values of m is proportional to |m|.
• The absolute time scale remains illustrative. There is no single real speed to show, because there is no observable time evolution in |Psi|^2. What is preserved, and is physically meaningful, is the direction and the proportionality to m.

Printing the wavefunction and exporting LaTeX

.venv/bin/python -m hydrogen_orbitals --name 3dxy --print-wavefunction --no-plot
.venv/bin/python -m hydrogen_orbitals --n 4 --l 2 --m 1 --latex --no-plot

--print-wavefunction shows the energy E_n, the radial part R_{n,l}(r), the angular part and the full wavefunction, formatted with sympy's unicode pretty printer. --latex prints the LaTeX source of those same expressions.

Saving images and animations

# A still image, rendered off screen: no window opens. This is the figure at the top.
.venv/bin/python -m hydrogen_orbitals --n 5 --l 4 --m 4 --plot-method isosurface-phase \
    --limit 42 --n-points 150 \
    --screenshot docs/images/hero-5g-phase.png --window-size 1100 800

# One full turn of the phase, as a GIF: the animation below
.venv/bin/python -m hydrogen_orbitals --n 4 --l 3 --m 3 --plot-method isosurface-phase \
    --limit 30 --n-points 130 \
    --gif docs/images/phase-animation.gif --window-size 600 450

Both options render off screen, so they need no window and work over ssh on a machine with no display, provided the OpenGL libraries are there. Directories missing from the path are created. --screenshot also switches the animation off, since a still image has nothing to animate.

--gif writes exactly one turn of the color map, with the closing frame left out, so the file loops without a stutter when it restarts. It needs a complex orbital with m != 0 and a method that carries phase, point-cloud or isosurface-phase; ask for it outside those cases and the tool says which condition failed instead of writing a file of identical frames.

Every figure in this README is generated by bash scripts/make_figures.sh, which is also the place to look for the exact parameters behind each one.

Performance and tuning

• The typical size of an orbital grows roughly as n^2. The defaults for --limit, --n-candidates and --n-points already scale with n, but for n >= 5, or simply for more detail, it pays to raise --n-candidates (for point-cloud) or --n-points (for the isosurface-* methods), at the cost of time and memory.
• The isosurface-* methods build a full 3D grid, that is n_points^3 evaluations, so they are heavier on memory than point-cloud at high resolution. Raising --iso-frac gives a smaller, simpler surface.

Project layout

hydrogen_orbitals/       Main package
  wavefunctions.py         Symbolic physics (sympy): Psi_nlm, energy, aliases
  viz_utils.py             Shared scene, reference grid and phase-to-color helpers
  plotting.py              The four visualization methods (pyvista)
  cli.py                   Command line interface
  __main__.py              Enables `python -m hydrogen_orbitals ...`
tests/                   Test suite, runs without a display
scripts/
  make_figures.sh          Regenerates every figure in this README
docs/images/             Those figures
requirements.txt         Dependencies (numpy, sympy, pyvista, imageio, matplotlib, pytest)

Troubleshooting

• No window opens, or an error mentions OpenGL, GLX or the display. The machine needs working OpenGL libraries and a graphical environment (X11 or Wayland). On a bare Debian or WSL install these are often missing: libgl1 and libglx-mesa0 provide them. On a headless server there is nothing to open a window on, so use --screenshot or --gif, which render off screen; that path additionally needs OSMesa (libosmesa6).
• ModuleNotFoundError for numpy, sympy or pyvista. The virtual environment is not active, or the dependencies are not installed. Re-run .venv/bin/pip install -r requirements.txt.
• The isosurface comes out empty. No point on the grid reaches --iso-frac of the maximum density. Lower --iso-frac or raise --n-points.
• The surface has elliptical holes on the faces of the box, showing its inside. The box is too small: the density still exceeds the iso level when it reaches the wall, and marching cubes cannot close a surface where it has no voxel on the far side, so it leaves the mesh open there. The hole belongs to the box, not to the orbital, which is why it is worth catching: it looks like a node. The tool prints a warning naming how many vertices sit on a face; raise --limit until it stops.
• The computation takes too long for large n. Lower --n-candidates or --n-points, or raise --iso-frac for a simpler surface.

Scope and limitations

• Covers stationary states of the hydrogen-like atom with configurable Z, in the real or complex basis, for any valid n, l, m.
• Does not cover superpositions of states (time dependent wave packets), relativistic corrections, or genuinely multi-electron atoms.

License

Released under the MIT License. See LICENSE for the full text.