Задание 4¶
Предварительные ласки.
In [1]:
import numpy as np
import scipy.special
import scipy.misc
In [2]:
from tqdm import tqdm, trange
In [3]:
import time
In [4]:
from IPython.display import display,Image
Запишем функцию для расчёта электронной плотности на сетке.
In [5]:
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
# Номировочный коэффициет
a = np.sqrt((2 / n / a0) ** 3 * scipy.special.factorial(n - l - 1) / (2 * n * scipy.special.factorial(n + l)))
# Радиальная функция
R = lambda r, n, l: a * (2 * r / n / a0) ** l * np.exp(-r / n / a0) * scipy.special.genlaguerre(n - l - 1, 2 * l + 1)(2 * r / n / a0)
# Волновая функция
WF = lambda r, theta, phi, n, l, m: 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)
Возьмём функцию Андрея Демкива и дополним вывод чисел до точности 16 знаков после запятой, так как у нас маленькие величины электронной плотности.
In [6]:
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 %.16f %.16f %.16f\n" %(start[0], start[1], start[2]))
cube_file.write(" %i %.16f 0.000000 0.000000\n" %(grid.shape[0], step[0]))
cube_file.write(" %i 0.000000 %.16f 0.000000\n" %(grid.shape[1], step[1]))
cube_file.write(" %i 0.000000 0.000000 %.16f\n" %(grid.shape[2], step[2]))
cube_file.write(" 1 0.000000 %.16f %.16f %.16f\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("%.16f " %(float(grid[x, y, z])))
i += 1
elif i == 5:
cube_file.write("%.16f\n" %(float(grid[x, y, z])))
i = 0
return 0
Рассчёт плотностей¶
Зададим размер ячейки и шаг сетки, затем переберём все квантовые числа для n от 1 до 3 включительно.
Ввиду того, что рассчитанные значения электронной плотности получаются очень маленькими, PyMol не может их красиво отобразить. Поэтому мы будем брать корень электронной плотности (просто модуль волновой функции).
In [7]:
for n in trange(1, 4):
d = 7 * n
step = d / 15
for l in trange(n):
for m in trange(-l, l + 1):
grid = np.sqrt(w(n, l, m, d))
name = f'{n}-{l}-{m}'
npy2cube(grid,
(- d, - d, - d),
(step, step, step),
name + '.cube')
Визуализация орбиталей в PyMol¶
In [8]:
import __main__
__main__.pymol_argv = ['pymol', '-x']
import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [9]:
cmd.bg_color('black')
cmd.set('antialias', 2)
cmd.set('ray_trace_mode', 3)
Сделаем красивую цветовую схему.
In [10]:
cmd.volume_ramp_new('ramp009', [-0.015, 1.00, 1.00, 1.00, 0.50,
-0.01, 1.00, 0.50, 0.00, 0.20,
-0.005, 0.00, 0.00, 0.50, 0.00,
0.005, 0.00, 0.00, 0.50, 0.00,
0.01, 1.00, 0.50, 0.00, 0.20,
0.015, 1.00, 1.00, 1.00, 0.50])
Загрузим орбитали.
In [11]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', name, ramp='ramp009')
Отрисуем орбитали в срезе.
In [57]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 90)
cmd.turn('y', 67.5)
name = f'{n}-{l}-{m}'
cmd.clip('slab', 5)
time.sleep(0.5)
cmd.do(f'png volumes\\{name}.png, 1080, 1080')
Посмотрим на орбитали.
In [58]:
quantum_numbers = [(n, l, m) for n in range(1, 4) for l in range(n) for m in range(-l, l + 1)]
In [59]:
plot_grid = int(np.ceil(np.sqrt(len(quantum_numbers))))
In [60]:
import matplotlib.pyplot as plt
In [61]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Сравним с реальными данными.
In [67]:
display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/e/e7/Hydrogen_Density_Plots.png/1024px-Hydrogen_Density_Plots.png'))
Вроде ок. Проблема с дырками в центре плотностей вызвана непонятно чем.
Теперь посмотрим на формы орбиталей.
In [68]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
name = f'{n}-{l}-{m}'
time.sleep(0.5)
cmd.do(f'png volumes\\{name}.png, 1080, 1080')
In [67]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Всё время какие-то бублики, а мы помним, что так не должно быть.
Попробуем рисовать сумму реальной и мнимой частей.
In [214]:
def w_sum(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
# Номировочный коэффициет
a = np.sqrt((2 / n / a0) ** 3 * scipy.special.factorial(n - l - 1) / (2 * n * scipy.special.factorial(n + l)))
# Радиальная функция
R = lambda r, n, l: a * (2 * r / n / a0) ** l * np.exp(-r / n / a0) * scipy.special.genlaguerre(n - l - 1, 2 * l + 1)(2 * r / n / a0)
# Волновая функция
WF = lambda r, theta, phi, n, l, m: 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
calc_WF = WF(r(x, y, z),
theta(x, y, z),
phi(x, y, z),
n, l, m)
# возвращаем сумму реальной и комплексной частей.
return np.real(calc_WF) + np.imag(calc_WF)
In [215]:
for n in trange(1, 4):
d = 7 * n
step = d / 15
for l in trange(n):
for m in trange(-l, l + 1):
grid = w_sum(n, l, m, d)
name = f'{n}-{l}-{m}'
npy2cube(grid,
(- d, - d, - d),
(step, step, step),
name + '_sum.cube')
In [216]:
cmd.reinitialize()
cmd.bg_color('black')
cmd.set('antialias', 2)
cmd.set('ray_trace_mode', 3)
Сделаем красивую цветовую схему.
In [217]:
cmd.volume_ramp_new('ramp009', [-0.015, 1.00, 1.00, 1.00, 0.50,
-0.01, 1.00, 0.50, 0.00, 0.20,
-0.005, 0.00, 0.00, 0.50, 0.00,
0.005, 0.00, 0.00, 0.50, 0.00,
0.01, 1.00, 0.50, 0.00, 0.20,
0.015, 1.00, 1.00, 1.00, 0.50])
Загрузим орбитали.
In [218]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.load(f'{name}_sum.cube')
cmd.volume(f'{name}-vol', f'{name}_sum', ramp='ramp009')
Отрисуем орбитали в срезе.
In [219]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 90)
cmd.turn('y', 67.5)
name = f'{n}-{l}-{m}'
cmd.clip('slab', 5)
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_sum.png, 1080, 1080')
In [220]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}_sum.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Кажется, что получилось лучше. По крайней мере, меньше дырок.
Теперь посмотрим на формы орбиталей.
In [221]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
name = f'{n}-{l}-{m}'
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_sum.png, 1080, 1080')
In [222]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}_sum.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Здесь уже не бублики, а более привычные нам бобы. Что же верно? Посмотрим, что выдаст ORCA.
ORCA¶
Закинули такой файл:
! 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 4
H 0 0 0
*
Теперь визуализируем орбитали.
In [228]:
cmd.reinitialize()
cmd.bg_color('black')
cmd.set('antialias', 2)
cmd.set('ray_trace_mode', 3)
Сделаем красивую цветовую схему.
In [229]:
cmd.volume_ramp_new('ramp009', [-0.015, 1.00, 1.00, 1.00, 0.50,
-0.01, 1.00, 0.50, 0.00, 0.20,
-0.005, 0.00, 0.00, 0.50, 0.00,
0.005, 0.00, 0.00, 0.50, 0.00,
0.01, 1.00, 0.50, 0.00, 0.20,
0.015, 1.00, 1.00, 1.00, 0.50])
Загрузим орбитали.
In [230]:
for n in trange(14):
name = f'H-{n}'
cmd.load(f'{name}.cube')
cmd.volume(f'{name}-vol', f'{name}', ramp='ramp009')
Отрисуем срезы орбиталей.
In [232]:
for n in trange(14):
name = f'H-{n}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.clip('slab', 0)
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_orca.png, 1080, 1080')
In [224]:
len(quantum_numbers)
Out[224]:
In [225]:
orca_quantum_numbers = quantum_numbers[:6] + quantum_numbers[9:] + quantum_numbers[6: 9]
In [227]:
len(orca_quantum_numbers)
Out[227]:
In [233]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for n in range(14):
plt.subplot(plot_grid, plot_grid, n + 1)
name = f'H-{n}'
numbers = orca_quantum_numbers[n]
img = plt.imread(f'volumes\{name}_orca.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'{numbers}', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Посмотрим теперь на внешний вид орбиталей.
In [234]:
for n in trange(14):
name = f'H-{n}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_orca.png, 1080, 1080')
In [235]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for n in range(14):
plt.subplot(plot_grid, plot_grid, n + 1)
name = f'H-{n}'
numbers = orca_quantum_numbers[n]
img = plt.imread(f'volumes\{name}_orca.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'{numbers}', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Видимо, "дырки" в орбиталях — это окей. Правильно получается складывать реальную и мнимую части волновой функции, чтобы получить хорошие формы орбиталей. Орбитали у ORCA немного другие по размерам.
Посмотрев на содержимое выходных файлов ORCA, можно увидеть, что для отображения чисел используется научная нотация, степени чисел около -16. Видимо, точность в 16 знаков для float недостаточна, поэтому мы тоже будем выводить в научной нотации. Будем рисовать квадрат модуля волновой функции.
In [121]:
def npy2cube_sci(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 %.6E %.6E %.6E\n" %(start[0], start[1], start[2]))
cube_file.write(" %i %.6E 0.000000 0.000000\n" %(grid.shape[0], step[0]))
cube_file.write(" %i 0.000000 %.6E 0.000000\n" %(grid.shape[1], step[1]))
cube_file.write(" %i 0.000000 0.000000 %.6E\n" %(grid.shape[2], step[2]))
cube_file.write(" 1 0.000000 %.6E %.6E %.6E\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("%.6E " %(grid[x, y, z]))
i += 1
elif i == 5:
cube_file.write("%.6E\n" %(grid[x, y, z]))
i = 0
return 0
Получим электронную плотность и запишем её.
In [122]:
for n in trange(1, 4):
d = 7 * n
step = d / 15
for l in trange(n):
for m in trange(-l, l + 1):
grid = w(n, l, m, d)
name = f'{n}-{l}-{m}'
npy2cube_sci(grid,
(- d, - d, - d),
(step, step, step),
name + '_sci.cube')
Покажем в PyMol.
In [195]:
cmd.reinitialize()
cmd.bg_color('black')
cmd.set('antialias', 2)
cmd.set('ray_trace_mode', 3)
Сделаем красивую цветовую схему, чтобы она была чувствительна к низким значениям.
In [196]:
cmd.volume_ramp_new('ramp009', [-0.015e-2, 1.00, 1.00, 1.00, 0.50,
-0.01e-2, 1.00, 0.50, 0.00, 0.20,
-0.005e-2, 0.00, 0.00, 0.50, 0.00,
0.005e-2, 0.00, 0.00, 0.50, 0.00,
0.01e-2, 1.00, 0.50, 0.00, 0.20,
0.015e-2, 1.00, 1.00, 1.00, 0.50])
Загрузим орбитали.
In [197]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.load(f'{name}_sci.cube')
cmd.volume(f'{name}-vol', f'{name}_sci', ramp='ramp009')
Отрисуем орбитали в срезе.
In [204]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 90)
cmd.turn('y', 67.5)
name = f'{n}-{l}-{m}'
cmd.clip('slab', 5)
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_sci.png, 1080, 1080')
In [205]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}_sci.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
Ура, оно получилось!
Теперь посмотрим на формы орбиталей.
In [200]:
for n in trange(1, 4):
for l in trange(n):
for m in trange(-l, l + 1):
name = f'{n}-{l}-{m}'
cmd.hide()
cmd.show('volume', f'{name}-vol')
cmd.reset()
cmd.turn('x', 45)
cmd.turn('y', 45)
cmd.turn('z', 45)
name = f'{n}-{l}-{m}'
time.sleep(0.5)
cmd.do(f'png volumes\\{name}_sci.png, 1080, 1080')
In [201]:
plt.figure(figsize=(5 * plot_grid, 5 * plot_grid))
for i, numbers in enumerate(quantum_numbers):
plt.subplot(plot_grid, plot_grid, i + 1)
n, l, m = numbers
name = f'{n}-{l}-{m}'
img = plt.imread(f'volumes\{name}_sci.png')
plt.imshow(img)
text_params = {'ha': 'center', 'va': 'center', 'family': 'sans-serif',
'fontweight': 'bold', 'fontsize': '18'}
plt.text(900, 1000, f'({n}, {l}, {m})', color='white', **text_params)
plt.axis('off')
plt.tight_layout()
А получаются те же бублики, а не бобы.
Выводы¶
- Почему-то квадрат модуля волновой функции даёт неправильные формы орбиталей, что заметно на n > 2.
- Обычная сумма реальной и мнимой частей волновой функции даёт похожие на правильные формы орбиталей.