Structural bionformatics. Homework №4.¶
Домашняя работа делается на своей машине и используется третий питон, поэтому файл npy2cube был скорректирован (map объекты переделаны в list):
start = list(map(lambda x: x/bohr_to_angs, start))]
step = list(map(lambda x: x/bohr_to_angs, step))
Elen Tevanyan
Inspiration resources: </br>
import numpy
import scipy.special
import scipy.misc
import npy2cube
Формулы, задающие волновую функцию <br. /> $$ r (x,y,z)=\sqrt{x^2+y^2+z^2}$$ $$ \theta(x,y,z)=arccos\Big(\frac{z}{r(x,y,z)}\Big)$$ $$ \phi(x,y,z)=arctg\Big(\frac{y}{x}\Big)$$
Волновая функция задается как произведение радиальной и угловой частей: $$ WF = R\cdot Y^m_n(\theta,\phi), $$ где $$R = \Big(\frac{2r}{na_0}\Big)^le^{-\frac{r}{na_0}}L^{2l+1}_{n-l-1}\frac{2r}{na_0},$$ т.е. $$ \psi_{nlm}(r,\theta,\phi) = \Big(\frac{2r}{na_0}\Big)^le^{-\frac{r}{na_0}}L^{2l+1}_{n-l-1}\frac{2r}{na_0}\cdot Y^m_n(\theta,\phi), $$ где
- n – основное число (1,2,3..)
- l – орбитальноечисло (0,1,2..n-1)
- m – магнитноечисло (-l..+l)
Реализация волновой функции
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.
#Радиальная часть функции
R = lambda r,n,l: (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*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-файлов.
d = 30
step = float(2.*d/29)
# Зададим цикл по перебору квантовых чисел
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)
name ='%s-%s-%s' % (n,l,m)
npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
#запустим PyMol и украдем у себя же часть кода из ДЗ2
import time
from IPython.display import display,Image
from xmlrpc.client import ServerProxy
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
defaultImage = '/tmp/pymolimg.png'
def prepareImage(width=300, height=300, sleep=2, filename=defaultImage):
## To save the rendered image
cmd.ray(width, height)
cmd.png('/tmp/pymolimg.png')
time.sleep(sleep)
#Прорисуем один из файлов.
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('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/1-0-0.cube')
cmd.volume('vol', '1-0-0')
cmd.volume_color ('vol','ramp007')
prepareImage()
Image(defaultImage)
cmd.set("ray_volume")
#я задаю два угла рассмотрения, т.к. при одной версии один из cube-файлов не отображается
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.01, 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 = '/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/'+name+'.cube'
cmd.load(file)
cmd.volume('vol', name)
cmd.volume_color ('vol','ramp007')
cmd.png('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/SB4/'+name+'.png')
names.append(name)
for name in names:
display(Image('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/SB4/'+name+'.png'))
Построение орбиталей в ORCA и сравнение
Эта ячейка была запущена в питона на shadbox, чтобы не мучиться с переносами файла на сервер
%%writefile h2.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
*
Код на bash
ssh -p 22025 -L 8888:localhost:8888 shad@vsb.fbb.msu.ru
cp -r ./progs/bin/orca ./hse/tevanyan
cd ./hse/tevanyan
orca h.inp
Нагенерились cube.файлы, которые я утаскиваю к себе на машину, чтобы нарисовать в PyMOL
scp -r -P 22025 shad@vsb.fbb.msu.ru:./hse/tevanyan/H-*.cube ./Study/tmp
cmd.set("ray_volume")
#я задаю два угла рассмотрения, т.к. при одной версии один из cube-файлов не отображается
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, \
])
files = []
for i in range(1,5):
name='H-'+'%s'%i
file = '/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/'+name+'.cube'
cmd.load(file)
cmd.volume('vol', name)
cmd.png('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/'+name+'.png')
files.append('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/'+name+'.png')
Можно сказать, что определенная схожесть есть
for adress in files:
print(adress)
display(Image(adress))
