Python Forum

Plotting of wave functions for an electrons quantum numbers (n,l ,ml) using spherical harmonics.
The first plot shows the probability of different distances which integral always is 1.

from scipy import *
from pylab import *
from numpy import *
import sys
import matplotlib.pyplot as plt
from scipy import integrate

R10 = lambda r: 4*exp(-2*r)*r**2                 
R20 = lambda r: 0.5*(1-r+r**2/4)*exp(-r)*r**2
R30 = lambda r: 1/(9**2*3)*((6-4*r+4*r**2/9)*exp(-r/3))**2*r**2

r = linspace(0,35,1000)
r100 = R10(r)
r200 = R20(r)
r300 = R30(r)

plt.plot(r, r100, label='|R10|²r²')
plt.plot(r, r200, label='|R20|²r²')
plt.plot(r, r300, label='|R30|²r²')
plt.title('Probability density')
plt.xlabel('aₒr')
plt.ylabel('P(aₒr)')
plt.legend()
plt.show()

I1 = integrate.quad(R10, 0, np.inf)[0]
I2 = integrate.quad(R20, 0, np.inf)[0]
I3 = integrate.quad(R30, 0, np.inf)[0]

print('Normalization integral from 0 to infinity')
print('<100|100> =', round(I1,8))
print('<200|200> =', round(I2,8))
print('<300|300> =', round(I3,8))

Calculation of the expectation value for the electron-nucleous distance in unit of a Bohr radius.

from scipy import *
from pylab import *
from numpy import *
import sys
import matplotlib.pyplot as plt
from scipy import integrate

R10 = lambda r: 4*exp(-2*r)*r**3                  
R20 = lambda r: 0.5*(1-r+r**2/4)*exp(-r)*r**3
R30 = lambda r: 1/(9**2*3)*((6-4*r+4*r**2/9)*exp(-r/3))**2*r**3

r = linspace(0,35,1000)
r100 = R10(r)
r200 = R20(r)
r300 = R30(r)

I1 = integrate.quad(R10, 0, np.inf)[0]
I2 = integrate.quad(R20, 0, np.inf)[0]
I3 = integrate.quad(R30, 0, np.inf)[0]

print('Integral from 0 to infinity')
print('<100|r|100> =', round(I1,8),'aₒ')
print('<200|r|200> =', round(I2,8),'aₒ')
print('<300|r|300> =', round(I3,8),'aₒ')

Expectation value for the potential energy

from scipy import *
from pylab import *
from numpy import *
import sys
import matplotlib.pyplot as plt
from scipy import integrate

R10 = lambda r: 4*exp(-2*r)*r                  
R20 = lambda r: 0.5*(1-r+r**2/4)*exp(-r)*r
R30 = lambda r: 1/(9**2*3)*((6-4*r+4*r**2/9)*exp(-r/3))**2*r

r = linspace(0,20,1000)
r100 = R10(r)
r200 = R20(r)
r300 = R30(r)

plt.plot(r, r100, label='|R10|²r')
plt.plot(r, r200, label='|R20|²r')
plt.plot(r, r300, label='|R30|²r')
plt.title('Expectation value for the potential energy')
plt.xlabel('r')
plt.ylabel('e²/4πεₒaₒ')
plt.legend()
plt.show()

I1 = integrate.quad(R10, 0, np.inf)[0]
I2 = integrate.quad(R20, 0, np.inf)[0]
I3 = integrate.quad(R30, 0, np.inf)[0]

print('Integral from 0 to infinity')
print('<r100|Epot|r100> =', round(I1,8),'e²/4πεₒaₒ ≈', round(I1*27.2, 2), 'eV')
print('<r200|Epot|r200> =', round(I2,8),'e²/4πεₒaₒ ≈',  round(I2*27.2, 2), 'eV')
print('<r300|Epot|r300> =', round(I3,8),'e²/4πεₒaₒ ≈',  round(I3*27.2, 2), 'eV')