import numpy
import math
import scipy.special
import scipy.misc
from IPython.display import display,Image
Функция Андрея Демкива:
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
Волновая функция:
def w(n,l,m,d):
#n,l,m – квантовые числа
#координаты x, y, z
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.arcsin(z/r(x,y,z)) + math.pi/2
#угол фи
phi = lambda x,y,z: numpy.arctan2(y, x) + math.pi
#радиус Бора
a0 = 1
#R – радиус функция
#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)
#модуль волновой функции
absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
final = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
return numpy.real(final) + numpy.imag(final)
Расчитаем значения для первых четырех энергетических уровней:
#шаг grid при заданном диапазоне от -d до d
for n in range(0,5):
d = 10*n
step = d/10
for l in range(0,n):
for m in range(0,l+1,1):
grid= w(n, l, m, d)
name='%s-%s-%s' % (n,l,m)
npy2cube(grid,(-d, -d, -d),
(step, step, step),
name+'.cube')
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd,stored
from IPython.display import Image
cmd.do('''
reini
bg_color black
''')
Цветовая гамма для отображения орбиталей:
cmd.volume_ramp_new('colors01', [-0.015, 1.00, 1.00, 0.00, 0.50,
-0.01, 1.00, 0.00, 0.00, 0.20,
-0.005, 0.00, 0.00, 1.00, 0.00,
0.005, 0.00, 0.00, 1.00, 0.00,
0.01, 1.00, 0.00, 0.00, 0.20,
0.015, 1.00, 1.00, 0.00, 0.50])
Загружаем все орбитали, присваивая им соответсвующие имена и применяя цветовую палитру:
for n in range(0,5):
for l in range(0,n):
for m in range(0,l+1,1):
name = f'{n}-{l}-{m}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='colors01')
Некоторые из полученных орбиталей:
Image('1-0-0.png')
Image('2-1-1.png')
Image('3-1-0.png')
Image('3-2-0.png')
Image('3-2-1.png')
Image('4-1-0.png')
Image('4-3-0.png')
with open('h.inp', 'w') as file:
file.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, в результате получаем файлы формата .cube, которые можно открыть в Pymol:
cmd.do('''
reini
bg_color black
''')
Цветовая гамма:
cmd.volume_ramp_new('colors02', [-0.015, 0.50, 0.00, 1.00, 0.50,
-0.01, 0.00, 1.00, 0.0, 0.30,
-0.005, 0.00, 0.00, 1.00, 0.00,
0.005, 0.00, 0.00, 1.00, 0.00,
0.01, 0.50, 1.00, 0.00, 0.30,
0.015, 0.00, 0.00, 1.00, 0.50])
Открываем файлы:
for i in range(1,5):
cmd.load(f'H-{i}.cube')
cmd.volume(f'H-{i}-vol', f'H-{i}', ramp='colors02')
Примеры полученных орбиталей:
Image('H-1.png')
Image('h-3.png')
Из графиков видно, что форма первых орбиталей, полученных с помощью Orca, примерно совпадает с формой орбиталей, полученных с помощью кода. Есть разница в соотношении размеров орбиталей:
Результат работы кода:
Image('pr.png')
Результат работы Orca:
Image('orca.png')