Ab initio вычисления орбиталей водорода

import numpy as np
import scipy.special
import scipy.misc
from math import factorial

from IPython.display import Image

import npy2cube

import matplotlib.pyplot as plt

def w(n, # главное квантовое число
      l, # орбитальное квантовое число
      m, # магнитное квантовое число
      d, # границы решётки, на которой считаем функцию,
      step
      ):
    
    
    # решётка точек, в которых будем считать
    # x,y,z = np.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]  
    # start:stop:step, комплексный шаг 30j = расставить равномерно 30 точек   
    x,y,z = np.mgrid[-d:d:step,-d:d:step,-d:d:step]  

    
    r = lambda x,y,z: np.sqrt(x**2+y**2+z**2)    # radius to origin
    theta = lambda x,y,z: np.arccos(z/r(x,y,z))  # zenith angle
    phi = lambda x,y,z: np.arctan(y/x)           # azimuthal angle

    a0 = 1.  # боровский радиус
    
    # нормировка
    norm = lambda n,l:  (( factorial(n-l-1) / (2*n*factorial(n+l)) )**0.5)  * ((2/n/a0)**(3/2)) 
    
    # радиальная функция    
    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)
    
    # волновая функция: нормировка * радиальная функция * сферические гармоники
    WF = lambda r,theta,phi,n,l,m: norm(n,l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
    
    # # волновая функция одним выражением
    # WF = lambda r,theta,phi,n,l,m:  ( factorial(n-l-1) / (2*n*factorial(n+l)) )**0.5  * \
    #                                 (2/n/a0)**(3/2)  * \
    #                                 numpy.exp(-r/n/a0) * (2*r/n/a0)**l  * \
    #                                 scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0) * \
    #                                 scipy.special.sph_harm(m,l,phi,theta)        
    
    return WF( r(x,y,z), theta(x,y,z), phi(x,y,z), n, l, m )

Должно получиться что-то такое:

Image(url='https://upload.wikimedia.org/wikipedia/commons/thumb/8/8f/Atomic_orbitals_spdf_m-eigenstates.png/1920px-Atomic_orbitals_spdf_m-eigenstates.png', width=600)

Проекции

Почему-то мой pymol испытывает проблемы с отображением volume орбиталей третьего периода. Возможно, это также проблема с подрезкой по изоповерхностям, но, к сожалению, мой pymol подвисает при попытках сделать что-либо кроме volume>default (в частности, при построении объёма с заданной кастомной ramp или попытке открыть панель управления цветом объёма/ramp). Поэтому ниже для них я привожу dots-репрезентацию вместо volume.

Возможно, это что-то не так с open-source pymol, установленным на WSL. Но и open-source и лицензионный pymol, установленные на Windows, вообще выдавают ошибку при попытке открыть сгенерированные кодом выше cube-файлы.

Для проверки, что считается всё верно, я ниже визуализирую проекции получившихся орбиталей как 2d тепловые карты. В данном случае они рисуются с одинаковым масштабам по пространственным координатам, но с собственной нормализацией цвета для каждого графика.

def plot_orbital(grid, n, l, m):
    fig, axes = plt.subplots(1, 3, figsize=(6, 2))

    axes[0].imshow(np.transpose(grid.sum(axis=2)))
    axes[0].set_xlabel('x')
    axes[0].set_ylabel('y')
    axes[0].set_xticks([])
    axes[0].set_yticks([])

    axes[1].imshow(np.transpose(grid.sum(axis=0)))
    axes[1].set_xlabel('y')
    axes[1].set_ylabel('z')
    axes[1].set_xticks([])
    axes[1].set_yticks([])

    axes[2].imshow(np.transpose(grid.sum(axis=1)))
    axes[2].set_xlabel('x')
    axes[2].set_ylabel('z')
    axes[2].set_xticks([])
    axes[2].set_yticks([])

    fig.suptitle(f'{n}-{l}-{m}', size=10)
    plt.tight_layout()
    fig.show()

d = 10
step = 0.3

# Зададим цикл по перебору квантовых чисел
for n in range(0,4):  # период
    for l in range(0,n):  # форма орбитали: 0=s, 1=p, 2=d, 3=f, 4=g
        for m in range(0, l+1): # ориентация орбитали
            #  вообще-то от -l до l, но если рисовать только плотность, она одинакова для m, симметричных относительно нуля
                        
            wf_grid =  w(n, l, m, d, step)
            
            # квадрат модуля волновой функции — электронная плотность
            prob_grid = np.absolute( wf_grid ) ** 2    
            
            plot_orbital(prob_grid, n, l, m)
            
            npy2cube.npy2cube(prob_grid, (-d, -d, -d),(step, step, step), f'{n}-{l}-{m}.cube')

Volume для периодов 1-2:

# получены скриптом вида
for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            cmd.load(f'{n}-{l}-{m}.cube')
            cmd.volume(f'{n}-{l}-{m}_volume', f'{n}-{l}-{m}')
            cmd.hide('everything', f'{n}-{l}-{m}')
            cmd.png(f'{n}-{l}-{m}_vol2')

for n in range(1,3):
    fig = plt.figure(figsize = (2*n, 2*n))
    for l in range(0,n):
        for m in range(0, l+1): 
            ax = fig.add_subplot(n, n, 1 + l * n +  m )
            ax.imshow(plt.imread(f'{n}-{l}-{m}_vol.png'))
            ax.set_title(f'{n}-{l}-{m}')
            ax.set_xticks([])
            ax.set_yticks([])
plt.tight_layout()

Dots для периодов 1-3

# получены скриптом вида
for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            cmd.load(f'{n}-{l}-{m}.cube')
            cmd.hide('extent', f'{n}-{l}-{m}')
            cmd.show('dots', f'{n}-{l}-{m}')
            cmd.png(f'{n}-{l}-{m}_vol2')

for n in range(1,4):
    fig = plt.figure(figsize = (2*n, 2*n))
    for l in range(0,n):
        for m in range(0, l+1): 
            ax = fig.add_subplot(n, n, 1 + l * n +  m )
            ax.imshow(plt.imread(f'{n}-{l}-{m}_dots_new.png'))
            ax.set_title(f'{n}-{l}-{m}')
            ax.set_xticks([])
            ax.set_yticks([])
    fig.tight_layout()
plt.show()

$P_x$, $P_y$

Можно заметить, что на всех изображениях выше p орбитали выглядят не так, как привычно с курсов химии — вместо трёх гантелек $p_x, p_y, p_z$ — одна гантелька $p_0$(=$p_z$) и два тора $p_{-1}$, $p_1$. Привычные $p_x$ и $p_y$ можно получить как линейные комбинации $p_{-1}$ и $p_1$:

d = 5
step = 0.3
n = 2
l = 1

wf_grid_1 =  w(n, l, -1, d, step)
wf_grid_2 =  w(n, l, 1, d, step)

prob_grid1 = np.absolute( wf_grid_1 + wf_grid_2 ) ** 2
prob_grid2 = np.absolute( wf_grid_1 - wf_grid_2 ) ** 2

plot_orbital(prob_grid2, n, l, 'x')
plot_orbital(prob_grid1, n, l, 'y')

npy2cube.npy2cube(prob_grid1, (-d, -d, -d),(step, step, step), f'2p_xy_1.cube')
npy2cube.npy2cube(prob_grid2, (-d, -d, -d),(step, step, step), f'2p_xy_2.cube')

0

fig = plt.figure(figsize = (5, 5*3))

ax = fig.add_subplot(1, 3, 1)
ax.imshow(plt.imread(f'2-1-0_vol.png'))
ax.set_xticks([])
ax.set_yticks([])

ax = fig.add_subplot(1, 3, 2)
ax.imshow(plt.imread(f'3d_xy_2_vol.png'))
ax.set_xticks([])
ax.set_yticks([])

ax = fig.add_subplot(1, 3, 3)
ax.imshow(plt.imread(f'3d_xy_1_vol.png'))
ax.set_xticks([])
ax.set_yticks([])

plt.tight_layout()

Проверка в Psi

import psi4
import numpy as np
psi4.core.set_output_file('output.dat')

Основываясь на https://sharifsuliman.medium.com/visualizing-frontier-orbitals-with-psi4-and-python-c01f335afaad и https://forum.psicode.org/t/calculating-atomization-energies/900/8

hatom = psi4.geometry('C')
psi4.set_options({'reference': 'uhf'})
eng, wave_function = psi4.energy('b3lyp/6-311G(2DF,P)', molecule=hatom, return_wfn=True)

cubeprop_tasks = []
cubeprop_tasks.append('ORBITALS')
# cubeprop_tasks.append('DENSITY')

psi4.set_options(
{
  'scf_type': 'df',
  'g_convergence': 'gau_tight',
  'freeze_core': 'true',
  'CUBEPROP_TASKS': cubeprop_tasks,
  'CUBEPROP_FILEPATH': './',
  'CUBEPROP_ORBITALS': list(range(1, 15)),
}
)

psi4.cubeprop(wave_function)

# скрипт для получения картинок выглядел примерно так
cmd.show('dots')
cmd.hide('extent')
for i in range(1, 15):
    cmd.disable('*')
    cmd.enable(f'Psi_a_{i}')
    cmd.png(f'Psi_a_{i}_dots')
    
    cmd.disable('*')
    cmd.volume('vol', f'Psi_a_{i}')
    cmd.enable('vol')
    cmd.png(f'Psi_a_{i}_vol')
    cmd.delete('vol')

fig = plt.figure(figsize = (2*3, 2*5))
for N in range(1,15):
        ax = fig.add_subplot(5, 3, N)
        ax.imshow(plt.imread(f'Psi_a_{N}_vol.png'))
        ax.set_title(f'Psi_a_{N}')
        ax.set_xticks([])
        ax.set_yticks([])
fig.tight_layout()
plt.show()

Volume показывает только половину реальной орбитали — вероятно, это как-то связано с распределением их по спинам, отрицательными числами в файле и обозначениями альфа-бета в psi4, и можно было бы покрасить вторую половину через ramp на отрицательные сигма, но это сложно сделать без визуальной панели и я недоразобралась.

Можно также было бы подумать, что DENSITY вместо ORBITALS в параметрах psi4 будет рисовать то что надо (электронную плотность), но это не так, она рисует что-то совсем иное.

Зато привычные симметричные картинки можно увидеть через dots:

fig = plt.figure(figsize = (2*3, 2*5))
for N in range(1,15):
        ax = fig.add_subplot(5, 3, N)
        ax.imshow(plt.imread(f'Psi_a_{N}_dots.png'))
        ax.set_title(f'Psi_a_{N}')
        ax.set_xticks([])
        ax.set_yticks([])
fig.tight_layout()
plt.show()

Видно, что посчитанные молекулярные орбитали (в присутствии многих электронов, т.к. атом углерода) имеют вид как в левой части картинки ниже, в то время как рассчитанные нами индивидуальные атомные (в присутсвии только одного электрона) выглядят как на правой части картинки ниже.

Как было показано выше на примере p-орбиталей, они получаются линейной комбинацией соответствующих атомных.

Image(url='https://upload.wikimedia.org/wikipedia/commons/thumb/2/2a/Atomic_orbitals_spdf_m-eigenstates_and_superpositions.png/1920px-Atomic_orbitals_spdf_m-eigenstates_and_superpositions.png', width=600)