Virtual Structural Biology Group » Учебные курсы » Структурная биоинформатика в ВШЭ

h1. Вычисление атомных орбиталей водорода 

Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами

Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol. 

* Запустите Jupiter Notebook или "Colab":https://colab.research.google.com/drive/1Jp9IzcD2YyHgrg7BScdZqF19LBLBygNy?usp=sharing .

* В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции

<pre><code class="python">import numpy

import scipy.special

import scipy.misc

</code></pre>

* Также вам понадобиться функция от Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py , скачайте его в рабочую директорию и подключите.

<pre><code class="python">import npy2cube

</code></pre>

* Давайте зададим волновую функцию, попробуйте предоставить формулу в отчёте

<pre><code class="python">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)

</code></pre>

* Вставьте в код комментарии про каждую внутреннюю функцию (lambda)

* Давайте рассчитаем значения для первых трех уровней. Функция w выдает трехмерный массив из 30*30*30 элементов с неким шагом (или grid). 

* определите шаг grid при заданном диапозоне от -d до d

<pre><code class="python">d=.....

step= .....

# Зададим цикл по перебору квантовых чисел

for n in range(0,4):

    for l in range(0,n):

        for m in range(0,l+1,1):

            grid= .... 

            name='%s-%s-%s' % (n,l,m)

            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению

            npy2cube.npy2cube(grid,(-d,....),(step,.....),name+'.cube')

</code></pre>  

* В результате работы скрипта появятся cube файлы, попробуйте их открыть в PyMol и построить volume для них или воспользуйтесь  Ipyvolume (https://github.com/maartenbreddels/ipyvolume).

* Попробуйте изменить окраску с помощью panel в colors

* Давайте создадим скрипт для PyMol для визуализации всех файлов

<pre><code class="python"> 

### Откуда эти цифры? 

pml='''

### cut below here and paste into script ###

cmd.volume_ramp_new('ramp007', [\

      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, \

    ])

### cut above here and paste into script ###

'''

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)

            # Загрузить cube файл

            # Отобразить электронную плотность (hint: volume)

            # Покрасить ее разумным образом 

# запишите переменную в файл

v=open(.....

</code></pre>

* Модифицируйте скрипт как Вам нравится для наилучшего отображения

* Загрузите изображения в Ipython Notebook

<pre><code class="python">from IPython.display import display,Image

display(Image=(file))

</code></pre>

* Рассчитаем орбитали в модуле Psi4

<pre>

import psi4

import numpy as np

psi4.core.set_output_file('output.dat')

</pre>

. Создадим input:

<pre>

g = '''

         0 1

         C     0.000000     0.000000     0.000000

'''

</pre>

* запустим

<pre>

m = psi4.geometry(g)

psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True})

ener=psi4.energy('scf/cc-pvtz', molecule = m )

Расчёт орбиталей в Psi4:  информация  об api https://psicode.org/psi4manual/1.4.0/api/psi4.core.cubeproperties

</pre>

 * Сравните визуально Ваши орбитали и рассчитанные программой Orca 

h2. Альтернативный способ визуализации орбиталей

* Сохраним плотность в файл hdf5

<pre>

import h5py

step = [cell[0]/100,cell[1]/100,cell[2]/100]

origin = [0,0,0]

with h5py.File('cube.hdf5', 'w') as f:

    oset = f.create_dataset("origin", data = origin)

    sset = f.create_dataset("step", data = step)

    dset = f.create_dataset("grid", data = e)

</pre>

* Создадим файл py, который будет запускаться из PyMol

<pre>

from chempy.brick import Brick

from pymol import cmd

import numpy as np

import h5py

F = h5py.File('cube.hdf5', 'r')

step = F['step'][()]

origin = F['origin'][()]

data = F['grid'][()] 

b = Brick.from_numpy(data, step, origin)

cmd.load_brick(b, 'cone')

volname =  'cone_volume'

cmd.volume(volname, 'cone')

F.close()

</pre>