In [37]:
import numpy
import scipy.special
import scipy.misc
import npy2cube
from IPython.display import display,Image
from xmlrpclib import ServerProxy
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
cmd.reinitialize()
Волновая функция:
$${\psi}_{n,l,m}(\rho,\theta,\phi) = \sqrt{\left({\frac{2}{na_0}}\right)^3\frac{(n-l-1)!}{2n(n+l)!}}e^{-\frac{\rho}{2}}{{\rho}^l}L_{n-l-1}^{2l+1}(\rho)Y_{l}^{m}(\theta,\phi)$$
- $\rho = \frac{2r}{na_0}$,
- $n$, $l$ и $m$ - квантовые числа: основное число $n$, орбитальное число $l$ и магнитное число $m$
In [19]:
def w(n,l,m,d):
x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2) #сферические координаты
theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
phi = lambda x,y,z: numpy.arctan(y/x)
a0 = 1.0
R = lambda r,n,l: (2.0*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2.0*r/n/a0) #радиальная часть волновой функции
WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta) #умножение радиальной части на угловую
absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2 #абсолютное значение волновой функции
return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
Расчет значений для первых трех уровней и создание Cube-файлов.¶
In [20]:
d = 30
step = float(2.*d/27)
In [21]:
# Зададим цикл по перебору квантовых чисел
for n in range(1,4):
for l in range(0,n):
for m in range(0,l+1):
grid = w(n,l,m,d)
name='%s-%s-%s' % (n,l,m)
npy2cube.npy2cube(grid,(-d, -d, -d),(step, step, step),name+'.cube')
In [11]:
#import os
#os.getcwd()
Out[11]:
Тест для одного файла, сохранение и вывод на экран¶
In [28]:
cmd.volume_ramp_new('ramp007', [\
0.002, 0.00, 0.00, 1.00, 0.00, \
0.01, 0.00, 1.00, 1.00, 0.20, \
0.015, 0.00, 0.00, 1.00, 0.00, \
])
cmd.set("ray_volume")
cmd.load('/home/kvlml/Desktop/hwgol/hw4/1-0-0.cube')
cmd.volume('vol', '1-0-0')
cmd.volume_color ('vol','ramp007')
cmd.png('/home/kvlml/Desktop/hwgol/hw4/png/'+'100test.png')
In [29]:
from IPython.display import display,Image
display(Image('/home/kvlml/Desktop/hwgol/hw4/png/'+'100test.png'))
Проделаем это же для всех¶
In [34]:
cmd.set("ray_volume")
cmd.volume_ramp_new('ramp007', [\
-0.035, 0.00, 0.00, 1.00, 0.00, \
-0.03, 0.00, 1.00, 1.00, 0.20, \
-0.025, 0.00, 1.00, 0.00, 0.00, \
0.025, 0.00, 1.00, 0.00, 0.00, \
0.03, 0.00, 1.00, 1.00, 0.20, \
0.035, 0.00, 0.00, 1.00, 0.00, \
])
names = []
for n in range(0,4):
for l in range(0,n):
for m in range(0,l+1,1):
name='%s-%s-%s' % (n,l,m)
file = '/home/kvlml/Desktop/hwgol/hw4/'+name+'.cube'
cmd.load(file)
cmd.volume('vol', name)
cmd.volume_color ('vol','ramp007')
cmd.png('/home/kvlml/Desktop/hwgol/hw4/png/'+name+'.png')
names.append(name)
In [35]:
for name in names:
display(Image('/home/kvlml/Desktop/hwgol/hw4/png/'+name+'.png'))
Сравнение с ORCA¶
In [4]:
%%writefile h3.inp
! UHF SVP XYZFile
%plots Format Cube
MO("H-1.cube",1,0);
MO("H-2.cube",2,0);
MO("H-3.cube",3,0);
MO("H-4.cube",4,0);
end
* xyz 0 4
H 0 0 0
*
In [5]:
#orca h3.inp
In [ ]:
cmd.reinitialize()
cmd.set("ray_volume")
cmd.volume_ramp_new('ramp007', [\
-0.05, 0.00, 0.00, 1.00, 0.00, \
-0.01, 0.00, 1.00, 1.00, 0.20, \
-0.015, 0.00, 0.00, 1.00, 0.00, \
0.005, 0.00, 0.00, 1.00, 0.00, \
0.01, 0.00, 1.00, 1.00, 0.20, \
0.015, 0.00, 0.00, 1.00, 0.00, \
])
names = []
for i in range(1,5):
name='H-'+'%s'%i
file = '/home/kvlml/Desktop/hwgol/hw4/'+name+'.cube'
cmd.load(file)
cmd.volume('vol', name)
cmd.volume_color ('vol','ramp007')
cmd.png('/home/kvlml/Desktop/hwgol/hw4/png/'+name+'.png')
names.append(name)
In [44]:
for name in names:
print(name)
display(Image('/home/kvlml/Desktop/hwgol/hw4/png/'+name+'.png'))
