import numpy
import scipy.special
import scipy.misc
from IPython.display import display,Image
from npy2cube import npy2cube
import math
def w(n,l,m,d):
# n,l,m – квантовые числа
# d – шаг
x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
# Переход к сферическим координатам, где r – радиус дистанции, thata, phi – углы относительно осей
r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
theta = lambda x,y,z: -numpy.arcsin(z/r(x,y,z)) + math.pi/2
phi = lambda x,y,z: numpy.arctan2(y, x) + math.pi
# Боровский радиус, наиболее вероятная дистанция между ядром и электроном
a0 = 1
# R – радиальная часть волновой функции, n, l – квантовые числа
# scipy.special.genlaguerre - обобщённый полином Лагерра
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)
# Добавочный коэффициент
k = lambda n,l: numpy.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
# Волновая функция, scipy.special.sph_harm – сферическая гармоника
WF = lambda r,theta,phi,n,l,m: k(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
# Плотность вероятности WF
absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
return numpy.real(result) + numpy.imag(result)
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
Вычисляем первые три уровня...
for n in range(1, 4):
d = 10 * n
step = d / 10
for l in range(0,n):
for m in range(-l,l+1,1):
grid= w(n, l, m, d)
name='%s_%s_%s' % (n,l,m)
npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')
Рассчитаем орбитали в программе Orca на kodomo. Создадим текстовый файл h.inp:
with open('orca.inp', 'w') as f:
f.write('''! UHF QZVPP XYZFile
%plots Format Cube
''')
for i in range(0,14):
f.write(f'MO("H-{i}.cube",{i},0);\n')
f.write('''end
* xyz 0 4
H 0 0 0
* ''')
Далее запустим следующие команды через Putty: export PATH=${PATH}:/srv/databases/orca
orca h.inp
C визуализацей при помощи скриптов возникли какие-то проблемы, поэтому сначала вручную были получены изображения отдельных орбиталей, а потом собраны в одну картинку при помощи Adobe Illustrator. Кроме того, возникали проблемы с отображением Volume. В результате были получены картинки ниже.
Image(filename='pr4_orbitals_npy2cube.png')
Image(filename='pr4_orbitals_orca.png')
display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/e/e7/Hydrogen_Density_Plots.png/1024px-Hydrogen_Density_Plots.png'))
Мы получили что-то странное. К сожалению, мне не удалось разобраться, в чём проблема. Предположу, что это какие-то ошибки в визуализации посредством Pymol. Volume отображается как-то неадекватно. Кажется, что обрезана половина изображения.
В целом оказалось так, что некоторые орбитали совпадают, и есть некоторые, не похожие на то, что должно быть. Многие орбитали, построенные orca не совпадают с теми, что построил наш скрипт.