In [ ]:
# Kernel: Python (psi4)
In [ ]:
#Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами
#Работа разделена на две части: расчёт плотности в Ipython Notebook с сохранением в формате CUBE и визуализация в Pymol.
#В нашем notebook cначала загрузим модули scipy и numpy для эффективной работы с массивами и содержащими нужные функции
In [57]:
import numpy
import scipy.special
import scipy.misc
import matplotlib.pyplot as plt
from IPython.display import Image
from xmlrpc.client import ServerProxy
import os, sys
In [58]:
# Также вам понадобиться функция от Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py
In [59]:
import npy2cube
In [42]:
#Зададим волновую функцию
In [52]:
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), absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
In [53]:
#Вставьте в код комментарии про каждую внутреннюю функцию (lambda)
#Давайте рассчитаем значения для первых трех уровней. Функция w выдает трехмерный массив из 30*30*30 элементов с неким шагом (или grid).
In [54]:
# определите шаг grid при заданном диапозоне от -d до d
step= 0.5
d = 15
# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
for l in range(0,n):
for m in range(0,l+1,1):
grid, prob= w(n,l,m,d)
name='%s-%s-%s' % (n,l,m)
npy2cube.npy2cube(prob,(-d,-d,-d),(step,step,step),name+'.cube')
In [55]:
# В результате работы скрипта появятся cube файлы, оценим их в 3Dmol
In [56]:
import py3Dmol
view = py3Dmol.view()
alpha = open('3-0-0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.04, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
In [50]:
#Сделаем картинки volume в PyMol из cube-файлов и посмотрим, что получилось
In [61]:
Image(filename='1-0-0.png') # 1s
Out[61]:
In [62]:
Image(filename='2-0-0.png') # 2s
Out[62]:
In [63]:
Image(filename='2-1-0.png') # 2p(z)
Out[63]:
In [64]:
Image(filename='2-1-1.png') # 2p (x or y)
Out[64]:
In [65]:
Image(filename='3-0-0.png') # 3s
Out[65]:
In [66]:
Image(filename='3-1-0.png') # 3p (z)
Out[66]:
In [67]:
Image(filename='3-1-1.png') # 3p (x or y)
Out[67]:
In [68]:
Image(filename='3-2-0.png') # 3d (z^2)
Out[68]:
In [69]:
Image(filename='3-2-1.png') # 3d (xz or yz)
Out[69]:
In [70]:
Image(filename='3-2-2.png') # 3d (x^2-y^2 or xy)
Out[70]:
In [26]:
#Рассчитаем орбитали в программе Psi4 :
In [71]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [72]:
# Зададим геометрию
g = '''
0 1
C 0.000000 0.000000 0.000000
'''
In [73]:
#Расчёт орбиталей в Psi4: информация об api https://psicode.org/psi4manual/1.4.0/api/psi4.core.cubeproperties
In [78]:
# Расчитаем энергию b орбитали
m = psi4.geometry(g)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" : True})
ener, wave=psi4.energy('scf/cc-pvtz', molecule = m, return_wfn=True)
In [79]:
#Сравните визуально Ваши орбитали и рассчитанные программой Psi4
In [82]:
cubeprop_tasks = ['ORBITALS']
psi4.set_options(
{
'scf_type': 'df',
'g_convergence': 'gau_tight',
'freeze_core': 'true',
'cubeprop_tasks': cubeprop_tasks,
'cubeprop_orbitals': list(range(1, 11)),
}
)
ener, wave = psi4.energy('scf/cc-pvtz', molecule = m, return_wfn=True)
psi4.cubeprop(wave)
In [ ]:
# Посмотрим полученные файлы в PyMol
In [ ]:
# Не получается правильно отрисовать volume, так как cubeprop не принимает numpy.absolute(wave)**2, поэтому отобразим в dots
In [83]:
Image(filename='1.png') #1s
Out[83]:
In [84]:
Image(filename='2.png') #2s
Out[84]:
In [85]:
Image(filename='3.png') # 2p
Out[85]:
In [86]:
Image(filename='4.png') # 2p
Out[86]:
In [87]:
Image(filename='5.png') # 2p
Out[87]:
In [88]:
Image(filename='6.png') # 3p
Out[88]:
In [89]:
Image(filename='7.png') # 3p
Out[89]:
In [90]:
Image(filename='8.png') # 3p
Out[90]:
In [91]:
Image(filename='9.png') # 3s
Out[91]:
In [95]:
Image(filename='10.png') # 3d
Out[95]:
In [ ]: