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

import numpy as np
import pandas as pd
import scipy.special
import scipy.misc
import math
import matplotlib.pyplot as plt
import time

Введем функцию npy2cube (автор: Андрей Демкив)

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 – квантовые числа
    #d – шаг
    
    #координаты в декартовой системе (30*30*30)
    x,y,z = np.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    
    #переход к сферическим координатам, r – радиус дистанции, thata, phi – углы относительно осей
    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 – радиальная часть волновой функции, n, l – квантовые числа
    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)
    
    #добавочный коэффициент
    k = lambda n,l: np.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)
    
    #плотность вероятности WF
    absWF = lambda r,theta,phi,n,l,m: np.absolute(WF(r,theta,phi,n,l,m))**2
    
    result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    
    return np.real(result) + np.imag(result)

Рассчитаем значения для первых трех уровней:

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')

Откроем cube-файлы в PyMol'е и построим для них volume:

import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd,stored

from pymol import cmd,stored
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(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')

#Сохраняем картинки орбиталей
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={480}, height={270}')

orbitales = {0: 's', 1: 'p', 2: 'd'}

_, axarr = plt.subplots(7, 2, figsize=(10, 21.33)) #будет 7 строк по 2 столбца с изображениями

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': 35}, pad=10)
            col += 1
            if col == 2:
                row += 1
                col = 0
plt.tight_layout()

Расчет орбиталей в программе Orca на kodomo

Создадим текстовый файл для работы Orca:

with open('hydrogen.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);
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
* ''')

Запустим Orca на kodomo, "скормив" ей файл "hydrogen.inp". Теперь снова построим изображения в PyMol и сравним результаты:

#Визуализация орбиталей, рассчитанных в Orca, в PyMol
cmd.reinitialize()
#Немного изменим цветовую гамму
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')

#Сохранение изображений

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={480}, height={270}')

orbitales = ['1s', '2s', '2p', '2p', '2p', '3s', '3d', '3d', '3d', '3d', '3d', '3p', '3p', '3p']
orbitales = (n for n in orbitales)

_, axarr = plt.subplots(7, 2, figsize=(10, 21.33))

row, col = 0, 0
for n in range(14):
    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'{next(orbitales)}', {'fontsize': 35}, pad=10)
    col += 1
    if col == 2:
        row += 1
        col = 0
plt.tight_layout()

Выводы:

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