Вычисление орбиталей водорода¶
0. Импорты и макрос¶
In [52]:
import numpy as np
import scipy.special
import scipy.misc
import math
import matplotlib.pyplot as plt
from textwrap import dedent
import time
Функция от Андрея Демкива:
In [60]:
def npy2cube(grid, start, step, cube_f):
"""Write .cube file from a grid.
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
"""
bohr_to_angs = 0.529177210859 # const
start = [x / bohr_to_angs for x in start]
step = [x / bohr_to_angs for x in step]
with open(cube_f, 'w') as cube_file:
# HEADER
cube_file.write('CPMD CUBE FILE.\n'
'OUTER 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 grid_item in np.nditer(grid, order='C'):
if i < 5:
cube_file.write('%f ' % (float(grid_item)))
i += 1
elif i == 5:
cube_file.write('%f\n' % (float(grid_item)))
i = 0
1. Орбитальные волновые функции своими силами¶
1) Зададим волновую функцию¶
По-идее, нам надо бы визуализировать плотность вероятности, однако она получается какой-то странной (для 2-1-1 она имеет форму бублика). Я честно пытался это исправить, но так и не смог. Однако, простое суммирование реальной и комплексой части дает замечательную картинку.
In [58]:
def w(n, l, m, d):
x, y, z = np.mgrid[-d:d:30j, -d:d:30j, -d:d:30j]
# Convert to Spherical coordinates.
r = np.sqrt(x**2 + y**2 + z**2)
theta = -np.arcsin(z / r) + math.pi/2
phi = np.arctan2(y, x) + math.pi
# Bohr radius.
a0 = 1.
# Normalization coefficient.
norm = math.sqrt((2/n/a0)**3 * math.factorial(n-l-1)
/ (2*n*math.factorial(n+l)))
# Radial function.
R = (2*r/n/a0)**l * np.exp(-r/n/a0) * \
scipy.special.genlaguerre(n-l-1, 2*l+1)(2*r/n/a0)
# Radial function * spherical harmonics --> Wave function.
WF = norm * R * scipy.special.sph_harm(m, l, phi, theta)
# Probability density.
absWF = np.absolute(WF)**2
return np.real(WF) + np.imag(WF)
2) Запишем в .cube файлы функции для первых трех уровней¶
In [61]:
for n in range(1, 4):
d = 7 * 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')
3) Визуализируем в пимоле (некоторые орбитали в разрезе)¶
In [2]:
import pymol
pymol.finish_launching(['pymol'])
from pymol import cmd
In [62]:
cmd.reinitialize()
# Color orbitals (x, r, g, b, a).
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])
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')
In [63]:
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}, width={1920//3}, height={1080//3}')
In [64]:
orb_names = {0: 's', 1: 'p', 2: 'd'}
_, axarr = plt.subplots(7, 2, figsize=(15, 32))
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].imshow(f)
axarr[row, col].axis('off')
axarr[row, col].set_title(f'{name} -- {n}{orb_names[l]}',
{'fontsize': 30}, pad=10)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()
2) Используем орку для построения орбиталей¶
На вход орке был предоставлен файл со следующим наполнением:
! 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);
MO("H-10.cube",10,0);
MO("H-11.cube",11,0);
MO("H-12.cube",12,0);
MO("H-13.cube",13,0);
end
* xyz 0 2
H 0 0 0
*
Посмотрим на ее выдачу в пимоле:
In [74]:
cmd.reinitialize()
# Color orbitals (x, r, g, b, a).
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])
for n in range(14):
name = f'H-{n}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp007')
In [75]:
slab = ['H-0', 'H-1', 'H-2', 'H-5', 'H-13', 'H-14']
for n in range(14):
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}, width={1920//3}, height={1080//3}')
In [76]:
names = ['1s', '2s', '2p', '2p', '2p', '3s', '3d', '3d', '3d', '3d', '3d', '3p', '3p', '3p']
names = (n for n in names)
_, axarr = plt.subplots(7, 2, figsize=(15, 32))
row, col = 0, 0
for n in range(14):
name = f'H-{n}'
f = plt.imread(f'{name}.png')
axarr[row, col].imshow(f)
axarr[row, col].axis('off')
axarr[row, col].set_title(f'{name} -- {next(names)}', {'fontsize': 30}, pad=10)
col += 1
if col == 2:
row += 1
col = 0
plt.tight_layout()
В целом, орка выдает похожие орбитали. Бросается в глаза то, что размеры орбиталей в ней странные, но это все из-за одинаковых координат в .cube файле (видимо orca никак их не подгоняет под реальные размеры орбиталей). Также d орбитали имеют странную форму, но совсем немного.