DEMO: H wavefunctions

Requirements | Importing libraries

import numpy as np
import  matplotlib.pyplot as plt
import plotly.graph_objects as go
from scipy.special import sph_harm, genlaguerre, factorial

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

try:
    from google.colab import output
    output.enable_custom_widget_manager()
    print('All good to go')
except:
    print('Okay we are not in Colab just proceed as if nothing happened')

Okay we are not in Colab just proceed as if nothing happened

1. Describing a radial function \(R_{nl}(r)\)

\[R_{nl} = \sqrt{ \Big( \frac{2}{n a_0} \Big)^3 \frac{(n-l-1)!}{2n(n+l)!}} \cdot e^{-\frac{r}{n a_0}} \cdot \Big( \frac{2 r}{n a_0} \Big)^l \cdot L^{2l+1}_{n-l-1}\Big( \frac{2 r}{n a_0} \Big)\]

def radial_function(r, 
                    n=1, 
                    l=0):
    '''Rnl(r) normalized radial function
    r: radial distance Float (0, inf)
    n: principal quantum number Int (1,2,3... inf)
    l: angular quantum number Int (0,1,2,... n-1)
    '''
    
    a0 = 1 # set Bohr radius to 1 (Hartree units)

    prefactor = np.sqrt( ((2 / n * a0) ** 3 * (factorial(n - l - 1))) / (2 * n * (factorial(n + l))) )

    laguerre = genlaguerre(n - l - 1, 2 * l + 1)

    p = 2 * r / (n * a0)
    

    return  prefactor * np.exp(-p / 2) * (p ** l) * laguerre(p)

r = np.linspace(0, 20, 1000)

Rnl = radial_function(r, n=2, l=0)

#plt.plot(r, Rnl)
#plt.plot(r, Rnl**2)
plt.plot(r, r**2 * Rnl**2)
plt.xlabel(r'$r [a_0]$')
plt.ylabel(r'$R_nl(r) r^2$')

Text(0, 0.5, '$R_nl(r) r^2$')

Check properties of radial distributon functions

r = np.linspace(0, 100, 10000)
dr = r[1]-r[0]

Rnl = radial_function(r, n=5, l=3)

pr = r**2 * Rnl**2 

norm = np.trapz(pr * dr) # Check normalizaton
#cdf  = np.cumsum(pr * dr) 

print(norm)

0.9999999052314471

np.trapz(pr*r*dr)

31.4999902089603

2. Describing an angular function | Spherical harmonic \(Y_{l}^{m}(\theta, \varphi)\)

• We will make use of sph_harm function from Scipy

• We can also build up spherical harmonics using Associated Legendre function using lpvm

\[ Y_{lm}(\phi,\theta) = \sqrt{\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!} } P_{lm}(cos \phi) \cdot e^{im\theta} \]

sph_harm(0, 0, np.pi, np.pi) # test a few spherical harmonics

(0.28209479177387814+0j)

def plot_Yml(l, m):
    '''Visualizing spherical harmonics using sph_harm funcion from Scipy special function library
     Note that the name of angles is different from the notation adopted in QM textbooks!
     l: angular quantum number (0,1,2,...)
     m: magnetic quantum number (-l,...l)
     '''

    #Creade grid of phi and theta angles for ploting surface mesh
    phi, theta = np.linspace(0, np.pi, 100), np.linspace(0, 2*np.pi, 100)
    phi, theta = np.meshgrid(phi, theta)

    #Calcualte spherical harmonic with given m and l
    Ylm = sph_harm(m, l, theta, phi) 
    R=abs(Ylm)
    
    # Let's normalize color scale
    fcolors    = Ylm.real
    fmax, fmin = fcolors.max(), fcolors.min()
    fcolors    = (fcolors - fmin)/(fmax - fmin)

    # Since we want to plot on cartesian reference frame we will use cartesian coordiniates x, y, z using R as the absolute value of Yml
    # Try a plot without R part. 
    fig = go.Figure(data=[go.Surface(x=R*np.sin(phi) * np.cos(theta), 
                                     y=R*np.sin(phi) * np.sin(theta), 
                                     z=R*np.cos(phi), 
                                     surfacecolor=fcolors,
                                     colorscale='balance')])

    # Show the plot
    fig.update_layout(title=fr'$Y_{l, m}$', autosize=False,
                      width=700, height=700,
                      margin=dict(l=65, r=50, b=65, t=90)
    )
    fig.show()

plot_Yml(l=1, m=0)

3. Describing the normalized probability as wavefunction squared \(|\psi _{nlm}(r,\theta ,\varphi)|^2\)

\[\psi_{nlm} = R_{nl}(r) \cdot Y_{l}^{m}(\theta, \varphi)\]

def normalized_wavefunction(n, 
                            l, 
                            m):

    '''Ψnlm(r,θ,φ) normalized wavefunction
    by definition of quantum numbers n, l, m and a bohr radius augmentation coefficient'''

    # Compute where radial function has nearly decayed to zero and determine cutoff distance
    r     = np.linspace(0, 500, 10000)
    pr = radial_function(r, n, l)**2 * r**2 * (r[1]-r[0])
    max_r = r[ np.where(np.cumsum(pr) >0.95)[0][0] ] 

    # Set coordinates grid to assign a certain probability to each point (x, y) in the plane
    x = y = np.linspace(-max_r, max_r, 501)
    x, y = np.meshgrid(x, y)
    r = np.sqrt((x ** 2 + y ** 2))

    # Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ)
    psi = radial_function(r, n, l) * sph_harm(m, l, 0, np.arctan(x / (y + 1e-7)))
    psi_sq = np.abs(psi)**2

    # Plot orbitals
    fig, ax = plt.subplots()
    ax.contour(psi_sq, 40, cmap='RdBu', extent=[-max_r, max_r,-max_r, max_r])
    ax.set_title(r'$|\psi_{{({0}, {1}, {2})}}|^2$'.format(n, l, m), fontsize=15)
    ax.set_aspect('equal')
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    ax.set_xlim(-max_r, max_r)
    ax.set_ylim(-max_r, max_r)

normalized_wavefunction(n=3, 
                            l=2, 
                            m=1)

Full Atomic orbitals in 3D

xyz = np.linspace(-10, 10, 51)
x,y,z = np.meshgrid(xyz, xyz, xyz, sparse=False)

r = np.sqrt(x**2 + y**2 + z**2)
phi = np.arctan2(y+1e-10, x)
theta = np.where( np.isclose(r, 0.0), np.zeros_like(r), np.arccos(z/r) )

n,l,m=4,1,0
psi = radial_function(r, n, l) * sph_harm(m, l, phi, theta) 

/tmp/ipykernel_2166/925093407.py:6: RuntimeWarning:

invalid value encountered in divide

fig= go.Figure(data=go.Isosurface(
    x=x.flatten(),
    y=y.flatten(),
    z=z.flatten(),
    value=abs(psi).flatten(),
    colorscale='RdBu',
    isomin=-.75*abs(psi).max(),
    isomax=.75*abs(psi).max(),
    surface_count=6,
    opacity=0.5,
    caps=dict(x_show=False,y_show=False,z_show=False)
))

fig.show()