Вычисление атомных орбиталей водорода¶
Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами
Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol.
- Запустите Jupiter Notebook на kodomo и если надо настройте туннель (см занятие Хемоинформатика) .
- В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции
1
2import numpy
3import scipy.special
4import scipy.misc
- Также вам понадобиться функция от Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py , скачайте его в рабочую директорию и подключите.
1
2import npy2cube
- Давайте зададим волновую функцию, попробуйте предоставить формулу в отчёте
1
2def w(n,l,m,d):
3
4 x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
5
6 r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
7 theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
8 phi = lambda x,y,z: numpy.arctan(y/x)
9
10 a0 = 1.
11
12 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)
13 WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
14 absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
15
16 return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
- Вставьте в код комментарии про каждую внутреннюю функцию (lambda)
- Давайте рассчитаем значения для первых трех уровней. Функция w выдает трехмерный массив из 30*30*30 элементов с неким шагом (или grid).
- определите шаг grid при заданном диапозоне от -d до d
1
2d=.....
3step= .....
4
5# Зададим цикл по перебору квантовых чисел
6
7for n in range(0,4):
8 for l in range(0,n):
9 for m in range(0,l+1,1):
10 grid= ....
11 name='%s-%s-%s' % (n,l,m)
12 # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
13 npy2cube.npy2cube(grid,(-d,....),(step,.....),name+'.cube')
- В результате работы скрипта появятся cube файлы, попробуйте их открыть в PyMol и построить volume для них или воспользуйтесь Ipyvolume (https://github.com/maartenbreddels/ipyvolume).
- Попробуйте изменить окраску с помощью panel в colors
- Давайте создадим скрипт для PyMol для визуализации всех файлов
1
2### Откуда эти цифры?
3pml='''
4### cut below here and paste into script ###
5cmd.volume_ramp_new('ramp007', [\
6 0.005, 0.00, 0.00, 1.00, 0.00, \
7 0.01, 0.00, 1.00, 1.00, 0.20, \
8 0.015, 0.00, 0.00, 1.00, 0.00, \
9 ])
10### cut above here and paste into script ###
11'''
12
13for n in range(0,4):
14 for l in range(0,n):
15 for m in range(0,l+1,1):
16 name='%s-%s-%s' % (n,l,m)
17 # Загрузить cube файл
18 # Отобразить электронную плотность (hint: volume)
19 # Покрасить ее разумным образом
20
21# запишите переменную в файл
22v=open(.....
- Модифицируйте скрипт как Вам нравится для наилучшего отображения
- Загрузите изображения в Ipython Notebook
1
2from IPython.display import display,Image
3display(Image=(file))
- Рассчитаем орбитали в программе Orca. Создадим текстовый файл h.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
*
запустим Orca
1
2export PATH=${PATH}:/tmp
3orca h.inp
Альтернативный способ визуализации орбиталей¶
- Сохраним плотность в файл hdf5
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)
- Создадим файл py, который будет запускаться из PyMol
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()