Task4

Version 4 (Andrey Golovin, 28.09.2020 19:24)

1 1 Andrey Golovin
h1. Вычисление атомных орбиталей водорода 
2 1 Andrey Golovin
3 1 Andrey Golovin
Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами
4 1 Andrey Golovin
5 1 Andrey Golovin
Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol. 
6 1 Andrey Golovin
7 1 Andrey Golovin
* Запустите Jupiter Notebook на kodomo и если надо настройте туннель (см занятие Хемоинформатика) .
8 1 Andrey Golovin
* В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции
9 1 Andrey Golovin
10 1 Andrey Golovin
11 4 Andrey Golovin
<pre><code class="python">import numpy
12 1 Andrey Golovin
import scipy.special
13 1 Andrey Golovin
import scipy.misc
14 1 Andrey Golovin
</code></pre>
15 1 Andrey Golovin
16 1 Andrey Golovin
* Также вам понадобиться функция от Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py , скачайте его в рабочую директорию и подключите.
17 1 Andrey Golovin
18 4 Andrey Golovin
<pre><code class="python">import npy2cube
19 1 Andrey Golovin
</code></pre>
20 1 Andrey Golovin
 
21 1 Andrey Golovin
* Давайте зададим волновую функцию, попробуйте предоставить формулу в отчёте
22 1 Andrey Golovin
23 4 Andrey Golovin
<pre><code class="python">def w(n,l,m,d):
24 1 Andrey Golovin
    
25 1 Andrey Golovin
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
26 1 Andrey Golovin
    
27 1 Andrey Golovin
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
28 1 Andrey Golovin
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
29 1 Andrey Golovin
    phi = lambda x,y,z: numpy.arctan(y/x)
30 1 Andrey Golovin
    
31 1 Andrey Golovin
    a0 = 1.
32 1 Andrey Golovin
    
33 1 Andrey Golovin
    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)
34 1 Andrey Golovin
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
35 1 Andrey Golovin
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
36 1 Andrey Golovin
      
37 1 Andrey Golovin
    return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
38 1 Andrey Golovin
</code></pre>
39 1 Andrey Golovin
40 1 Andrey Golovin
* Вставьте в код комментарии про каждую внутреннюю функцию (lambda)
41 1 Andrey Golovin
42 1 Andrey Golovin
* Давайте рассчитаем значения для первых трех уровней. Функция w выдает трехмерный массив из 30*30*30 элементов с неким шагом (или grid). 
43 1 Andrey Golovin
44 1 Andrey Golovin
45 1 Andrey Golovin
* определите шаг grid при заданном диапозоне от -d до d
46 1 Andrey Golovin
47 4 Andrey Golovin
<pre><code class="python">d=.....
48 1 Andrey Golovin
step= .....
49 1 Andrey Golovin
50 1 Andrey Golovin
# Зададим цикл по перебору квантовых чисел
51 1 Andrey Golovin
52 1 Andrey Golovin
for n in range(0,4):
53 1 Andrey Golovin
    for l in range(0,n):
54 1 Andrey Golovin
        for m in range(0,l+1,1):
55 1 Andrey Golovin
            grid= .... 
56 1 Andrey Golovin
            name='%s-%s-%s' % (n,l,m)
57 1 Andrey Golovin
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
58 1 Andrey Golovin
            npy2cube.npy2cube(grid,(-d,....),(step,.....),name+'.cube')
59 1 Andrey Golovin
</code></pre>  
60 1 Andrey Golovin
      
61 1 Andrey Golovin
* В результате работы скрипта появятся cube файлы, попробуйте их открыть в PyMol и построить volume для них или воспользуйтесь  Ipyvolume (https://github.com/maartenbreddels/ipyvolume).
62 1 Andrey Golovin
63 1 Andrey Golovin
64 1 Andrey Golovin
* Попробуйте изменить окраску с помощью panel в colors
65 1 Andrey Golovin
66 1 Andrey Golovin
* Давайте создадим скрипт для PyMol для визуализации всех файлов
67 1 Andrey Golovin
68 1 Andrey Golovin
<pre><code class="python"> 
69 1 Andrey Golovin
### Откуда эти цифры? 
70 1 Andrey Golovin
pml='''
71 1 Andrey Golovin
### cut below here and paste into script ###
72 1 Andrey Golovin
cmd.volume_ramp_new('ramp007', [\
73 1 Andrey Golovin
      0.005, 0.00, 0.00, 1.00, 0.00, \
74 1 Andrey Golovin
      0.01,  0.00, 1.00, 1.00, 0.20, \
75 1 Andrey Golovin
      0.015, 0.00, 0.00, 1.00, 0.00, \
76 1 Andrey Golovin
    ])
77 1 Andrey Golovin
### cut above here and paste into script ###
78 1 Andrey Golovin
'''
79 1 Andrey Golovin
80 1 Andrey Golovin
for n in range(0,4):
81 1 Andrey Golovin
    for l in range(0,n):
82 1 Andrey Golovin
        for m in range(0,l+1,1):
83 1 Andrey Golovin
            name='%s-%s-%s' % (n,l,m)
84 1 Andrey Golovin
            # Загрузить cube файл
85 1 Andrey Golovin
            # Отобразить электронную плотность (hint: volume)
86 1 Andrey Golovin
            # Покрасить ее разумным образом 
87 1 Andrey Golovin
88 1 Andrey Golovin
# запишите переменную в файл
89 1 Andrey Golovin
v=open(.....
90 1 Andrey Golovin
</code></pre>
91 1 Andrey Golovin
92 1 Andrey Golovin
93 1 Andrey Golovin
* Модифицируйте скрипт как Вам нравится для наилучшего отображения
94 1 Andrey Golovin
95 1 Andrey Golovin
* Загрузите изображения в Ipython Notebook
96 1 Andrey Golovin
97 1 Andrey Golovin
98 1 Andrey Golovin
<pre><code class="python"> 
99 1 Andrey Golovin
from IPython.display import display,Image
100 1 Andrey Golovin
display(Image=(file))
101 1 Andrey Golovin
</code></pre>
102 1 Andrey Golovin
  
103 2 Andrey Golovin
* Рассчитаем орбитали в программе Orca. Создадим текстовый файл h.inp:
104 1 Andrey Golovin
105 1 Andrey Golovin
<pre>
106 1 Andrey Golovin
! UHF SVP XYZFile
107 1 Andrey Golovin
%plots Format Cube
108 1 Andrey Golovin
 MO("H-1.cube",1,0);
109 1 Andrey Golovin
 MO("H-2.cube",2,0);
110 1 Andrey Golovin
 MO("H-3.cube",3,0);
111 1 Andrey Golovin
 MO("H-4.cube",4,0);
112 1 Andrey Golovin
end
113 1 Andrey Golovin
* xyz 0 4
114 1 Andrey Golovin
 H 0 0 0
115 1 Andrey Golovin
*
116 1 Andrey Golovin
</pre>
117 1 Andrey Golovin
118 1 Andrey Golovin
 запустим Orca
119 1 Andrey Golovin
120 1 Andrey Golovin
121 1 Andrey Golovin
<pre><code class="bash"> 
122 1 Andrey Golovin
export PATH=${PATH}:/tmp
123 1 Andrey Golovin
orca h.inp
124 1 Andrey Golovin
</code></pre>
125 1 Andrey Golovin
 * Сравните визуально Ваши орбитали и рассчитанные программой Orca 
126 3 Andrey Golovin
127 3 Andrey Golovin
h2. Альтернативный способ визуализации орбиталей
128 3 Andrey Golovin
 
129 3 Andrey Golovin
* Сохраним плотность в файл hdf5
130 3 Andrey Golovin
<pre>
131 3 Andrey Golovin
import h5py
132 3 Andrey Golovin
step = [cell[0]/100,cell[1]/100,cell[2]/100]
133 3 Andrey Golovin
origin = [0,0,0]
134 3 Andrey Golovin
135 3 Andrey Golovin
with h5py.File('cube.hdf5', 'w') as f:
136 3 Andrey Golovin
    oset = f.create_dataset("origin", data = origin)
137 3 Andrey Golovin
    sset = f.create_dataset("step", data = step)
138 3 Andrey Golovin
    dset = f.create_dataset("grid", data = e)
139 3 Andrey Golovin
140 3 Andrey Golovin
</pre>
141 3 Andrey Golovin
142 3 Andrey Golovin
143 3 Andrey Golovin
* Создадим файл py, который будет запускаться из PyMol
144 3 Andrey Golovin
<pre>
145 3 Andrey Golovin
from chempy.brick import Brick
146 3 Andrey Golovin
from pymol import cmd
147 3 Andrey Golovin
import numpy as np
148 3 Andrey Golovin
import h5py
149 3 Andrey Golovin
150 3 Andrey Golovin
151 3 Andrey Golovin
F = h5py.File('cube.hdf5', 'r')
152 3 Andrey Golovin
153 3 Andrey Golovin
step = F['step'][()]
154 3 Andrey Golovin
origin = F['origin'][()]
155 3 Andrey Golovin
data = F['grid'][()] 
156 3 Andrey Golovin
b = Brick.from_numpy(data, step, origin)
157 3 Andrey Golovin
158 3 Andrey Golovin
159 3 Andrey Golovin
cmd.load_brick(b, 'cone')
160 3 Andrey Golovin
volname =  'cone_volume'
161 3 Andrey Golovin
162 3 Andrey Golovin
cmd.volume(volname, 'cone')
163 3 Andrey Golovin
164 3 Andrey Golovin
F.close()
165 3 Andrey Golovin
</pre>