import npy2cube
import numpy as np
import scipy
from scipy import special
from scipy import misc
from scipy.special import factorial as fact
import ipyvolume as ipv
from IPython.display import display,Image
import os
формула волновой функции:
\begin{equation*} \psi_{nlm}(r,\theta,\phi)=\sqrt{\Bigl({\frac{2}{na_0}}\Big)^3\frac{(n-l-1)!}{2n(n+l)!}}e^{\frac{-p}{2}}p^lL_{n-l-1}^{2l+1}(p)Y_l^m(\theta,\phi) \end{equation*}
где р это
\begin{equation*} p=\frac{2r}{na_0} \end{equation*}
Стоит заметить, что в функции w в функции R отсутствует кое-какие множители. Исправим это и также прокомментируем остальные функции:
Спойлер : после того, как я это исправил и построил куб-файлы , были проблемы с масштабом , которые так и не получилось решить. Решил оставить такой вариант, потому что без домножения на корень получается более менее похожие орбитали на orca.
Также пробовал через ipyvolume, но почему-то он ipv.show() почему-то не отображала построений...
def w(n,l,m,d,points):
#n,l,m - quantum numbers, d - parameter for a grid
#l - orbital number (0,1,2,3,..,n-1)
#m - magnetic number (-l,..+l)
#n - basic number(1,2,3)
x,y,z = numpy.mgrid[-d:d:points*1j,-d:d:points*1j,-d:d:points*1j]
#go to spherical coordinates and calculate r,theta,phi parameters:
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)
#Bohr radius
a0 = 0.53 #* 10**(-10)#1.
coef_fac=numpy.sqrt(scipy.special.factorial(n-l-1) / (2*n*scipy.special.factorial(n+l))* (2/n/a0)**3)
#scaling for good representations
p = lambda x,y,z,n: 2*r(x,y,z)/n/a0*2**n
#radial dependence
R = lambda p,n,l: p**l * numpy.exp(-p/2) * scipy.special.genlaguerre(n-l-1,2*l+1)(p)
#wave function: R func multiplied by spherical harmonics
WF = lambda p,theta,phi,n,l,m: R(p,n,l) * scipy.special.sph_harm(m,l,phi,theta)
#absolute value of wave function
absWF = lambda p,theta,phi,n,l,m: numpy.absolute(WF(p,theta,phi,n,l,m))**2
return absWF(p(x,y,z,n),theta(x,y,z),phi(x,y,z),n,l,m)
#we create 30 (x,y,z) points in 3D coordinates. Each coordinate has a range from -d to d .
#Therefore step is equal to 2*d/(num_of_points-1)
points=30
d=3
step=2.*d/(points-1)
path="C:/Users/Mister_Ruvim/Desktop/hw4/cube/"
# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
for l in range(0,n):
for m in range(0,l+1,1):
grid= w(n,l,m,d,points)
name="fold/" + '%s-%s-%s' % (n,l,m) +'.cube'
print(name)
# для сохранения нужно задать координаты старта grid и шаг по каждому направлению
npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),cube_f=name)
f = open ('script.txt', 'w')
f.write('cmd.volume_ramp_new(\'ramp007\', [\
-0.03, 0.00, 0.00, 1.00, 0.00, \
-0.035, 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.045, 0.00, 0.00, 1.00, 0.00, \
])\n')
for n in range(1,4):
for l in range(0,n):
for m in range(0,l+1):
name='%s-%s-%s' % (n,l,m)
full_name = name+'.cube'
vol_name = name+'_volume'
f.write('load C:/Python27/Scripts/fold/%s, %s\n' % (full_name, name))
f.write('volume %s, %s\n' % (vol_name, name))
f.write('volume_color %s, ramp007\n' % (vol_name))
f.close()
Отобразим полученные картинки: ...почему-то некоторые из объемов не отобразились с помощью скрипта, но вручную при помощи action:volume(default) они отобразились.
path="C:/Python27/Scripts/pic/"
for f in os.listdir("pic/"):
print f
display(Image(path+f))
Создадим файл для orca:
%%writefile h_2.in
! UHF SVP XYZFile
%plots Format Cube
MO("H-0.cube",0,0);
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
*
Дальше с помощью WinSCP я скопировал на сервер h.inp, и запустил с помощью orca, предварительно скопировав экзешник в свою папку (mrPuk). Затем с помощью того же WinSCP выкачал куб-файлы.
Код в консоле:¶
cd mkdir ~/hse/mrPuk
cp -r /home/shad/progs/bin/orca ~/hse/mrPuk
cd ~/hse/mrPuk
orca h.in
Построим орков:
f = open ('script_orc.txt', 'w')
f.write('cmd.volume_ramp_new(\'ramp008\', [\
-0.045, 0.00, 0.00, 1.00, 0.00, \
-0.035, 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.045, 0.00, 0.00, 1.00, 0.00, \
])\n')
for i in range(0,5):
name='%s' % (i)
full_name = "H-"+name+'.cube'
vol_name = "H-"+name+'_volume'
f.write('load C:/Python27/Scripts/H/%s, %s\n' % (full_name, name))
f.write('volume %s, %s\n' % (vol_name, name))
f.write('volume_color %s, ramp008\n' % (vol_name))
f.close()
Получилось следующее: (картинки сохранял вручную)
path="C:/Python27/Scripts/pictures/"
for f in os.listdir("pictures/"):
print f
display(Image(path+f))
Схожесть, конечно, присутствует, однако заметно, что orca куда лучше строит