Вычисление атомных орбиталей водорода¶
В данном практикуме нам нужно, опираясь на уравнение, построить электронную плотность в одно-электронном атоме и сравнить с известными программами
Импортируем пакеты
In [1]:
import numpy as np
import scipy.special
import scipy.misc
import math
from os import listdir
from IPython.display import display,Image
import pandas as pd
import matplotlib.pyplot as plt
import time
In [2]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
import pymol
pymol.finish_launching()
from pymol import cmd, stored
Модифицируем функцию
In [3]:
#Функция Андрея Демкива http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py
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]:
def w(n,l,m,d):
# n,l,m - квантовые числа: главное, орбитальное, магнитное. d - шаг
# координаты внутри параллелепипеда в декартовой системе координат (30*30*30)
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.arcsin(z/r(x,y,z)) + math.pi/2
phi = lambda x,y,z: np.arctan2(y, x) + math.pi
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 расписан на слайде 14
b = lambda n,l: np.sqrt((math.factorial(n-l-1)/(2*n*math.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
out = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
#возвращаем реальную и мнимую части (иначе скрипт в следующей ячейке будет ругаться)
return np.real(out) + np.imag(out)
Нарисуем орбитали
In [5]:
# проходимся по всем квантовым числам
for n in range(1, 4):
d = 10 * n
# задаем шаг
step = 2 * d / 29
for l in range(0, n):
for m in range(-l, l+1):
# получаем волновую функцию
grid = w(n, l, m, d)
name = f'{n}_{l}_{m}'
# и сохраняем функцию в файл
npy2cube(grid, (-d,)*3, (step,)*3, f'{name}.cube')
In [53]:
cmd.reinitialize()
# задаем цвет орбиталей
cmd.volume_ramp_new('ramp007', [-0.015, 1.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 1.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, 0.00, 1.00, 1.00, 0.20,
0.015, 0.00, 0.00, 1.00, 0.00])
# отрисовываем орбитали в PyMol
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}_{l}_{m}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
cmd.load("file.pml") не сработал, поэтому сделаем ручками
In [54]:
#Сохраняем картинки орбиталей
slab = ['1_0_0', '2_0_0', '2_1_0', '3_0_0', '3_1_0', '3_2_0']
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}_{l}_{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 90)
if name in slab:
cmd.clip('slab', 0)
time.sleep(0.1)
cmd.do(f'png {name}')
orbitales = {0: 's', 1: 'p', 2: 'd'}
_, axarr = plt.subplots(7, 2, figsize=(15, 25))
row, col = 0, 0
for n in range(1, 4):
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}_{l}_{m}'
f = plt.imread(f'{name}.png')
axarr[row, col].axis('off')
axarr[row, col].imshow(f)
axarr[row, col].set_title(f'{n}{orbitales[l]}',
{'fontsize': 25}, pad=8)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()
In [31]:
# ради интереса посчитаем для n = 6
n = 6
d = 5 * n
step = 2 * d / 29
for l in range(0, n):
for m in range(-l, l+1):
grid = w(n, l, m, d)
name = f'{n}_{l}_{m}'
npy2cube(grid, (-d,)*3, (step,)*3, f'{name}.cube')
cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [-0.015, 1.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 1.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, 0.00, 1.00, 1.00, 0.20,
0.015, 0.00, 0.00, 1.00, 0.00])
for l in range(0, n):
for m in range(-l, l+1):
name = f'{n}_{l}_{m}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
In [33]:
Image(filename='pr4.png')
Out[33]:
Расчет орбиталей в программе Orca¶
Создадим файл с интересующимим нас орбиталями для программы Orca
In [34]:
with open('h.inp', 'w') as f:
f.write('''! UHF QZVPP XYZFile
%plots Format Cube
MO("H-0.cube",0,0);
MO("H-1.cube",1,0);
MO("H-2.cube",2,0);
MO("H-3.cube",3,0);
MO("H-4.cube",4,0);
MO("H-5.cube",5,0);
MO("H-6.cube",6,0);
MO("H-7.cube",7,0);
MO("H-8.cube",8,0);
MO("H-9.cube",9,0);
end
* xyz 0 2
H 0 0 0
* ''')
Запустили программу на кодомо с помощью команд export PATH=${PATH}:/srv/databases/orca orca h.inp
В результате работы программы Orca были получены файлы .cube для орбиталей
In [57]:
cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [-0.015, 1.00, 0.00, 0.00, 0.00,
-0.01, 1.00, 1.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, 0.00, 1.00, 1.00, 0.20,
0.015, 0.00, 0.00, 1.00, 0.00])
for n in range(10):
name = f'H-{n}'
cmd.load(f'./orca/{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
In [58]:
slab = ['H-{x}'.format(x=x) for x in range(10)]
for n in range(10):
name = f'H-{n}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
if name in slab:
cmd.clip('slab', 0)
time.sleep(0.1)
cmd.do(f'png {name}')
orbitales = {0: 's', 1: 'p', 2: 'd'}
_, axarr = plt.subplots(5, 2, figsize=(15, 25))
row, col = 0, 0
for n in range(10):
name = f'H-{n}'
f = plt.imread(f'{name}.png')
axarr[row, col].axis('off')
axarr[row, col].imshow(f)
axarr[row, col].set_title(f'{n}{orbitales[l]}',
{'fontsize': 25}, pad=8)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()
Выводы¶
Орбитали 1 и 2 энергенического уровня, в целом, сходятся между orca и нашей функцией. Орбитали высших порядков существенно отличаются. Видно, что "наша" шестая орбиталь совсем непохожа на полученную в orca. Я больше доверяю своим расчетам, потому что они выглядят более детализированными. Также заметна разница в размерах орбиталей. На рисунке ниже сравнение 2p орбиталей, полученных двумя подходами. "Наша" больше)
In [60]:
Image(filename='pr4_1.png')
Out[60]:
In [ ]: