History

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

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

Работа разделена на две части: расчёт плотности в 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 для них.
  • Попробуйте изменить окраску с помощью 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
{{{#!highlight python
from IPython.display import display,Image
display(Image=(file))
}}}
  • xyz 0 4
    H 0 0 0 *
    }}}
запустим Orca

{{{#!highlight bash
export PATH=${PATH}:/tmp
orca h.inp
}}}

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