python chemistry quantum mechanics sympy 12 June 2017

Radial density and sympy

Plotting quantum mechanical electron densities

A bit of theory

The Schrödinger equation in the center of mass coordinates for the Hydrogen-like ion is given by

\begin{align} \left(-\frac{\hbar}{2m}\nabla^2 + V\right)\psi=E\psi \end{align} the wave function is assumed to be separable into the radial and angular part according to the formula \begin{align} \psi(r, \theta, \phi) = R_{nl}(r)Y(\theta, \phi) \end{align}

the the radial Schrödinger equation is \begin{align} \left(-\frac{\hbar}{2m}\nabla^2_r + \frac{l(l + 1)}{2r^2} - \frac{Z}{r}\right)\psi=E\psi \end{align} where \begin{align} \nabla^2_r = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2 \frac{\partial}{\partial r}\right) = \frac{\partial^2}{\partial r^2} + \frac{2}{r}\frac{\partial}{\partial r} \end{align}

and the solution is $$ R_{nl}(r) = N_{nl}L^{2l+1}_{n+l}\left(\frac{2Zr}{n}\right)\exp\left(\frac{-Zr}{n}\right) $$ where

$$ N_{nl} = -\sqrt{\left(\frac{2Z}{n}\right)^3\frac{\left(n-l-1\right)!}{2n\left[\left(n+l\right)!\right]^3}} $$ and \(L\) are associated Laguerre polynomials

in python

from sympy import S, Symbol, simplify, latex
from sympy.physics.hydrogen import R_nl

def radial_density(n, l, r, Z=1):
    'calculate the radial density'

    r = S(r)
    r_nl = R_nl(n, l, r, Z)
    return r_nl**2*r**2

if you want to try it out in jupyter notebook

from sympy import init_printing
init_printing()