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>

In [1]:
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)

Реализация волновой функции

In [2]:
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-файлов.

In [3]:
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')
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/npy2cube.py:39: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f " %(float(grid[x, y, z])))
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/npy2cube.py:42: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f\n" %(float(grid[x, y, z])))
In [2]:
#запустим 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)
In [5]:
#Прорисуем один из файлов. 
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)
Out[5]:
In [26]:
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)
In [27]:
for name in names:
    display(Image('/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/SB4/'+name+'.png'))

Построение орбиталей в ORCA и сравнение
Эта ячейка была запущена в питона на shadbox, чтобы не мучиться с переносами файла на сервер

In [8]:
%%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
*
Overwriting h2.inp

Код на 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

In [9]:
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')

Можно сказать, что определенная схожесть есть

In [10]:
for adress in files:
    print(adress)
    display(Image(adress))
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/H-1.png
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/H-2.png
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/H-3.png
/Users/ElenTevanyan/Documents/Study/Структурная биоинформатика/H-4.png