2024
2023
2022
2021
2020
2018
2017
2016Установка
- Запустите терминал в юпитере
conda activate psi4 python -m ipykernel install --user --name psi4 --display-name "Python (Psi4)"
Вычисление атомных орбиталей водорода
Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами
Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol.
- В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции
Также вам понадобиться функция от Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py , скачайте его в рабочую директорию, адаптируйте и подключите.
1 import npy2cube
- Давайте зададим волновую функцию
1 def w(n,l,m,d):
2
3 x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
4
5 r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
6 theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
7 phi = lambda x,y,z: numpy.arctan(y/x)
8
9 a0 = 1.
10
11 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)
12 WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
13 absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
14
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).
1 # определите шаг grid при заданном диапозоне от -d до d
2 d=.....
3 step= .....
4
5 # Зададим цикл по перебору квантовых чисел
6
7 for 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 файлы, оценим их в 3Dmol
import py3Dmol
view = py3Dmol.view()
alpha = open('1-0-0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.04, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()Попробуйте их открыть в PyMol и построить volume для них.
- Попробуйте изменить окраску с помощью panel в colors
Давайте создадим скрипт для PyMol для визуализации всех файлов
1 ### Откуда эти цифры?
2 pml='''
3 ### cut below here and paste into script ###
4 cmd.volume_ramp_new('ramp007', [\
5 0.005, 0.00, 0.00, 1.00, 0.00, \
6 0.01, 0.00, 1.00, 1.00, 0.20, \
7 0.015, 0.00, 0.00, 1.00, 0.00, \
8 ])
9 ### cut above here and paste into script ###
10 '''
11
12 for n in range(0,4):
13 for l in range(0,n):
14 for m in range(0,l+1,1):
15 name='%s-%s-%s' % (n,l,m)
16 # Загрузить cube файл
17 # Отобразить электронную плотность (hint: volume)
18 # Покрасить ее разумным образом
19
20 # запишите переменную в файл
21 v=open(.....
- Модифицируйте скрипт как Вам нравится для наилучшего отображения
Можно посмотреть орбитали ввиде isosurface в nglviwer. Посмотрите презентиацию Евы Смородины https://docs.google.com/presentation/d/1S3E_YleBMWxl-ottIMsQAjPhMrQtmzSzFQvOOsPymY0/edit?usp=sharing
- Загрузите изображения в Jupyter Notebook
- Рассчитаем орбитали в программе Psi4 :
import psi4 import numpy as np psi4.core.set_output_file('output.dat') - Зададим геометрию
g = '''
0 1
C 0.000000 0.000000 0.000000
'''- Расчитаем энергию b орбитали
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- Сравните визуально Ваши орбитали и рассчитанные программой Psi4