Pr4: Вычисление атомных орбиталей водорода
In [1]:
import numpy as np
import scipy.special
import scipy.misc
In [2]:
#Описываем волновую функцию
def w(n,l,m,d):
x,y,z = np.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
# Переход на сферическую систему координат
r = lambda x,y,z: np.sqrt(x**2+y**2+z**2) # радиус
theta = lambda x,y,z: np.arccos(z/r(x,y,z)) # угол тета
phi = lambda x,y,z: np.arctan2(y,x) # угол фи
a0 = 1
R = lambda r,n,l: (2*r/n/a0)**l * np.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0) # Радиус функция
b = lambda n,l: np.sqrt((scipy.special.factorial(n-l-1)/(2*n*scipy.special.factorial(n+l)))*(2/n/a0)**3) # Коэффициент нормировки
WF = lambda r,theta,phi,n,l,m: b(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta) # Волновая функция
absWF = lambda r,theta,phi,n,l,m: np.absolute(WF(r,theta,phi,n,l,m))**2 # Модуль волновой функции
return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m) # на выход: модуль
In [3]:
#функция от Андрея Демкива
def npy2cube(grid, start, step, cube_f):
'''
PARAMETERS:
grid : numpy array
3-dimentional array, containing grid data
start : tuple
format: (x, y, z), coordinates of cube start point
step: tuple
format: (x, y, z), step size on 3 axes
cube_f: string
name of output .cube file
RETURNS:
void
'''
cube_string = ""
bohr_to_angs = 0.529177210859 #const
start = list(map(lambda x: x/bohr_to_angs, start))
step = list(map(lambda x: x/bohr_to_angs, step))
with open(cube_f, "w") as cube_file:
###HEADER###
cube_file.write(" CPMD CUBE FILE.\nOUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n")
cube_file.write(" 1 %f %f %f\n" %(start[0], start[1], start[2]))
cube_file.write(" %i %f 0.000000 0.000000\n" %(grid.shape[0], step[0]))
cube_file.write(" %i 0.000000 %f 0.000000\n" %(grid.shape[1], step[1]))
cube_file.write(" %i 0.000000 0.000000 %f\n" %(grid.shape[2], step[2]))
cube_file.write(" 1 0.000000 %f %f %f\n" %(start[0], start[1], start[2]))
###DATA###
i = 0
for x in range(grid.shape[0]):
for y in range(grid.shape[1]):
for z in range(grid.shape[2]):
if i < 5:
cube_file.write("%f " %(float(grid[x, y, z])))
i += 1
elif i == 5:
cube_file.write("%f\n" %(float(grid[x, y, z])))
i = 0
return 0
In [4]:
#Вычисляем первые три уровня
for n in range(1, 4):
d = 10 * n
step = d / 10
for l in range(n):
for m in range(0,l+1,1):
grid = np.sqrt(w(n, l, m, d))
name = name='%s-%s-%s' % (n,l,m)
npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')
In [5]:
#генерируем PyMol скрипт для визуализации полученных cube-файлов
with open('open_cube.pml', 'w') as scr:
scr.write('cmd.volume_ramp_new(\'ramp007\', \
[-0.015, 0.00, 1.00, 1.00, 0.50, \
-0.01, 0.20, 0.80, 1.00, 0.20, \
-0.005, 0.40, 0.60, 1.00, 0.00, \
0.005, 0.60, 0.40, 1.00, 0.00, \
0.01, 0.80, 0.20, 1.00, 0.20, \
0.015, 1.00, 0.00, 1.00, 0.50]) \n')
for n in range(1,4):
for l in range(n):
for m in range(0,l+1):
name = '%s-%s-%s' % (n,l,m)
cube_name = name + '.cube' # название cube-файла
volume_name = name + '_vol' # название volume
scr.write(f'load {cube_name}, {name}\n') # загружаем cube-файл в PyMol
scr.write(f'volume {volume_name}, {name}\n') # Отображаем электронную плотность volume
scr.write(f'volume_color {volume_name}, ramp007\n') # Окрашивание
In [6]:
from IPython.display import display, Image
In [7]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [8]:
cmd.load("open_cube.pml")
In [9]:
# 2s-орбиталь
cmd.do('png 2-0-0.png')
Image('2-0-0.png')
Out[9]:
In [10]:
# 2p-орбиталь
cmd.do('png 2-1-0.png')
Image('2-1-0.png')
Out[10]:
In [11]:
# орбиталь третьего уровня
cmd.do('png 3-2-1.png')
Image('3-2-1.png')
Out[11]:
In [12]:
# генерируем скрипт для запуска Orca
with open('h.inp', 'w') as f:
f.write('''! 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
*''')
# запускаем в kodomo
# на выходе получили 4 cube-файла
In [13]:
# генерируем аналогичный PyMol скрипт
with open('open_orca.pml', 'w') as scr:
scr.write('cmd.volume_ramp_new(\'ramp007\', \
[-0.015, 0.00, 1.00, 1.00, 0.50, \
-0.01, 0.20, 0.80, 1.00, 0.20, \
-0.005, 0.40, 0.60, 1.00, 0.00, \
0.005, 0.60, 0.40, 1.00, 0.00, \
0.01, 0.80, 0.20, 1.00, 0.20, \
0.015, 1.00, 0.00, 1.00, 0.50]) \n')
for i in range(1,5):
name = 'H-%s' % i
cube_name = name + '.cube' # название cube-файла
volume_name = name + '_vol' # название volume
scr.write(f'load orca/{cube_name}, {name}\n') # загружаем cube-файл в PyMol
scr.write(f'volume {volume_name}, {name}\n') # Отображаем электронную плотность volume
scr.write(f'volume_color {volume_name}, ramp007\n') # Окрашивание
In [14]:
cmd.reinitialize()
cmd.load("open_orca.pml")
In [15]:
cmd.do('png H-1.png')
Image('H-1.png')
Out[15]:
In [16]:
cmd.do('png H-2.png')
Image('H-2.png')
Out[16]:
In [17]:
cmd.do('png H-3.png')
Image('H-3.png')
Out[17]:
In [18]:
cmd.do('png H-4.png')
Image('H-4.png')
Out[18]:
In [ ]:
#Видно, что высчитанные программой Orca орбитали лежат во втором уровне
#H-1 является 2s, H-2:4 являются 2pxyz