User:Geek3/hydrogen-cut

hydrogen-cut is a Python script that creates images of 2D cuts through hydrogen atomic wavefunctions. For 3D orbitals see my other script hydrogen.

About hydrogen-cut

hydrogen-cut uses the SciPy library to compute the analytic functions of atomic orbitals. Images are written with misc.imsave, which requires the Python Imaging Library.

There are different ways in which orbitals are plotted, namely

• $\psi$ using colors to depict the phase
• $|\psi |$ white on black background
• $r|\psi |^{2}$ white on black background, which gives a brightness proportional to the electron density in each antinode if the orbital is spherically symmetric. Here, basically each pixel collects the probability density along a ring around the center through the pixel and perpendicular to the image plane. The circumference of that ring is 2πr, therefore the collected probability density is proportional to $r|\psi |^{2}$.
• $|x|\cdot |\psi |^{2}$ white on black background, which gives a brightness proportional to the electron density in each antinode if the orbital is symmetric around the z-axis. Here, basically each pixel collects the probability density along a ring around the z-axis through the pixel and perpendicular to the image plane. The circumference of that ring is 2π|x|, therefore the collected probability density is proportional to $|x|\cdot |\psi |^{2}$.

However:

• The probability density $|\psi |^{2}$ is not plotted here, because it drops too steeply towards the outer antinodes that these would be nearly invisible.
• $r^{2}|\psi |^{2}$ is also not shown. It gives the total probability density at radius r, but having additionally the larger area at larger r would over-represent the outer antinodes. $r^{2}|\psi |^{2}$ is more suited for 1D plots.

Code

#!/usr/bin/python3
# -*- coding: utf8 -*-

import numpy as np
import imageio
from scipy import special
from math import *

def Rnl(n, l, r):
    '''
    radial part of the wavefunction. r may be an array
    '''
    rho = np.abs(r) * 2 / n
    L = special.eval_genlaguerre(n - l - 1, 2 * l + 1, rho)
    p = np.product(np.arange(n - l, n + l + 1.0))
    return sqrt((2/n)**3 / (2*n*p)) * np.exp(-rho/2) * rho**l * L

def psinlm(n, l, m, x, y, z):
    '''
    hydrogen atom wavefunction. x, y, z may be arrays of equal length
    '''
    r = np.sqrt(x**2 + y**2 + z**2)
    phi = np.arctan2(y, x)
    theta = np.arctan2(np.sqrt(x**2 + y**2), z)
    R = Rnl(n, l, r)
    Ylm = special.sph_harm(m, l, phi, theta)
    return R * Ylm

def image_psi(n, l, m, imgsize):
    x = 3 * n**2 * ((np.arange(imgsize) + 0.5) / imgsize * 2 - 1)
    xv, zv = np.meshgrid(x, -x)
    psi = np.real(psinlm(n, l, m, xv, np.zeros_like(xv), zv))
    psi /= np.max(np.abs(psi))
    img_r = np.interp(psi, [-1, -.5, 0, .5, 1], [.5, .92, 1, 0, 0])
    img_g = np.interp(psi, [-1, -.5, 0, .5, 1], [0, 0, 1, 0, 0])
    img_b = np.interp(psi, [-1, -.5, 0, .5, 1], [0, 0, 1, .97, .5])
    img = np.stack((img_r, img_g, img_b), axis=2)
    img = np.clip(img * 256, 0, 255).astype('uint8')
    return img

def image_abspsi(n, l, m, imgsize):
    x = 3 * n**2 * ((np.arange(imgsize) + 0.5) / imgsize * 2 - 1)
    xv, zv = np.meshgrid(x, -x)
    psi = psinlm(n, l, m, xv, np.zeros_like(xv), zv)
    img = np.hypot(psi.real, psi.imag)
    img = np.clip(img * 256 / np.max(img), 0, 255).astype('uint8')
    return img

def image_rpsi2(n, l, m, imgsize):
    x = (2.2+1.1/n) * n**2 * ((np.arange(imgsize) + 0.5) / imgsize * 2 - 1)
    xv, zv = np.meshgrid(x, -x)
    psi = psinlm(n, l, m, xv, np.zeros_like(xv), zv)
    img = psi.real**2 + psi.imag**2
    img *= np.hypot(xv, zv)
    img = np.clip(img * 256 / np.max(img), 0, 255).astype('uint8')
    return img

def image_xpsi2(n, l, m, imgsize):
    x = (2.2+1.1/n) * n**2 * ((np.arange(imgsize) + 0.5) / imgsize * 2 - 1)
    xv, zv = np.meshgrid(x, -x)
    psi = psinlm(n, l, m, xv, np.zeros_like(xv), zv)
    img = psi.real**2 + psi.imag**2
    img *= np.abs(xv)
    img = np.clip(img * 256 / np.max(img), 0, 255).astype('uint8')
    return img

for n in range(1, 5+1):
    for l in range(n):
        for m in range(l + 1):
            fname = 'atomic-orbital-cut_psi_n{}l{}m{}.png'.format(n, l, m)
            print('drawing', fname)
            image = image_psi(n, l, m, imgsize=1600)
            imageio.imsave(fname, image)
            
            fname = 'atomic-orbital-cut_abs(psi)_n{}l{}m{}.png'.format(n, l, m)
            print('drawing', fname)
            image = image_abspsi(n, l, m, imgsize=1600)
            imageio.imsave(fname, image)
            
            if m == 0:
                fname = 'atomic-orbital-cut_r*abs(psi)**2_n{}l{}m{}.png'.format(n, l, m)
                print('drawing', fname)
                image = image_rpsi2(n, l, m, imgsize=1600)
                imageio.imsave(fname, image)
            else:
                fname = 'atomic-orbital-cut_x*abs(psi)**2_n{}l{}m{}.png'.format(n, l, m)
                print('drawing', fname)
                image = image_xpsi2(n, l, m, imgsize=1600)
                imageio.imsave(fname, image)

Images

See Wikimedia Commons search for a list of created images.