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

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

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

• Запустите Jupiter Notebook на kodomo и если надо настройте туннель (см занятие Хемоинформатика) .
• В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции

import numpy
import scipy.special
import scipy.misc

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

import npy2cube

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

 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)

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

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

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

 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')

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

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

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

 
 ### Откуда эти цифры? 
 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(.....

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

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

from IPython.display import display,Image
display(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, это исполняемый файл /home/shad/progs/bin/orca .

export PATH=${PATH}:/tmp
orca h.inp

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

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

• Сохраним плотность в файл 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()