-
Notifications
You must be signed in to change notification settings - Fork 0
/
3D-Orbitals.py
103 lines (50 loc) · 1.93 KB
/
3D-Orbitals.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy.integrate as integrate
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina')
import ipywidgets as widgets
import ipyvolume as ipv
import scipy.special as spe
def psi_R(r,n=15,l=8):
coeff = np.sqrt((2.0/n)**3 * spe.factorial(n-l-1) /(2.0*n*spe.factorial(n+l)))
laguerre = spe.assoc_laguerre(2.0*r/n,n-l-1,2*l+1)
return coeff * np.exp(-r/n) * (2.0*r/n)**l * laguerre
r = np.linspace(0,100,1000)
R = psi_R(r,n=15,l=8)
plt.plot(r, R**2, lw=3)
plt.xlabel('$r [a_0]$',fontsize=20)
plt.ylabel('$R_{nl}(r)$', fontsize=20)
plt.grid('True')
nmax=70
@widgets.interact(n = np.arange(1,nmax,1), l = np.arange(0,nmax-1,1))
def plot_radial(n=15,l=8):
r = np.linspace(0,250,10000)
psi2 = psi_R(r,n,l)**2 * (r**2)
plt.plot(r, psi2, lw=2, color='red')
plt.xlabel('$r [a_0]$')
plt.ylabel('$R_{nl}(r)$')
rmax = n**2*(1+0.5*(1-l*(l+1)/n**2))
plt.xlim([0, 2*rmax])
def psi_ang(phi,theta,l=8,m=0):
sphHarm = spe.sph_harm(m,l,phi,theta)
return sphHarm.real
phi, theta = np.linspace(0, np.pi, 100), np.linspace(0, 2*np.pi, 100)
phi, theta = np.meshgrid(phi, theta)
Ylm = psi_ang(theta,phi,l=2,m=0)
x = np.sin(phi) * np.cos(theta) * abs(Ylm)
y = np.sin(phi) * np.sin(theta) * abs(Ylm)
z = np.cos(phi) * abs(Ylm)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
fcolors = (Ylm - Ylm.min())/(Ylm.max() - Ylm.min())
ax.plot_surface(x, y, z, facecolors=cm.seismic(fcolors), alpha=0.3)
cset = ax.contour(x, y, z,20, zdir='z',offset = -1, cmap='summer')
cset = ax.contour(x, y, z,20, zdir='y',offset = 1, cmap='winter' )
cset = ax.contour(x, y, z,20, zdir='x',offset = -1, cmap='autumn')
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)