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

import numpy as np
import math
import scipy.special
import scipy.misc
from npy2cube import npy2cube
import time

Зададим волновую функцию

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.arctan(y/x)
    
    # Боровский радиус
    a0 = 1.
    
    # Радиальная функция, 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)
    
    # Нормировочный коэффициент
    base = lambda n,l: np.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
    # Волновая функция, m - магнитное квантовое число
    WF = lambda r,theta,phi,n,l,m: base(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

    # Из опыта коллег стало известно, что орибтали получаются более похожими на истинные,
    # если возвращать сумму истинной и мнимой частей.
    res = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    return np.real(res) + np.imag(res)

Зададим цикл по перебору квантовых чисел

for n in range(1, 4):
    d = 10 * n
    step = d / 10
    for l in range(0,n):
        for m in range(-l,l+1,1):
            grid= w(n, l, m, d)
            name='%s_%s_%s' % (n,l,m)
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')

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

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

Теперь визуализиурем полученные файлы в Pymol.

cmd.reinitialize()
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.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, 0.00, 0.20, 
      0.015, 0.00, 0.00, 0.00, 0.00])
names = []
for n in range(1, 4):
    for l in range(0, n):
        for m in range(-l, l+1):
            def_name = f'{n}_{l}_{m}'
            names.append(def_name)
            cmd.load(f'{def_name}.cube')
            cmd.volume(f'{def_name}-volume', def_name, ramp='colors')
            cmd.hide()
            cmd.show('volume', f'{def_name}-volume')
            cmd.reset()
            cmd.turn('x', 55)
            cmd.turn('y', 35)
            cmd.turn('z', 45)
            time.sleep(0.3)
            cmd.do(f'''
            png {def_name}.png, dpi=300
            ''')

Соберем все картинки на одно изображение.

import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

numbers = ['1_0_0', '2_0_0', '2_1_0', '2_1_1', '2_1_-1', '3_0_0', '3_1_0', '3_1_1', '3_1_-1', '3_2_0', '3_2_1', '3_2_-1', '3_2_2', '3_2_-2']
img_paths = [f'{i}.png' for i in numbers]
zaborchik = []
for img_path in img_paths:
    zaborchik.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(30,20))
columns = 3
for i, image in enumerate(zaborchik):
    f.add_subplot(int(len(zaborchik) / columns + 1), columns, i + 1)
    plt.imshow(image)
    plt.xticks([])
    plt.yticks([])
    plt.title(f'{orbitals[i]}')

plt.show()

Теперь сделаем все то же, но теперь воспольземся psi4.

import psi4

psi4.core.set_output_file('output.dat')
g = '''
H
'''
m = psi4.geometry(g)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True, 'reference':'uhf'})
ener, wave=psi4.energy('scf/cc-pvtz', return_wfn = True, molecule = m )
psi4.set_options(
{
  'scf_type': 'df',
  'g_convergence': 'gau_tight',
  'freeze_core': 'true',
  'CUBEPROP_TASKS': ['ORBITALS'],
  'CUBEPROP_FILEPATH': './',
  'CUBEPROP_ORBITALS': list(range(1, 17)),
}
)
psi4.cubeprop(wave)

Создадим список с именами файлов, полученных при работе psi4.

import os
fil=[i for i in os.listdir() if 'Psi_a' in i]

Снова визуализируем все.

cmd.reinitialize()
cmd.volume_ramp_new('colors', [-0.015, 0.00, 0.00, 0.00, 0.00, 
      -0.01,  1.00, 0.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, 0.00, 0.20, 
      0.015, 0.00, 0.00, 0.00, 0.00])
names = []
for n in fil:
    def_name=n.split('.')[0]
    cmd.load(f'./{def_name}.cube')
    cmd.volume(f'{def_name}-volume', def_name, ramp='colors')
    cmd.hide()
    cmd.show('volume', f'{def_name}-volume')
    cmd.reset()
    cmd.turn('x', 55)
    cmd.turn('y', 35)
    cmd.turn('z', 45)
    time.sleep(0.3)
    cmd.do(f'''
    png {def_name}.png, dpi=300
    ''')

import glob
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

zaborchik = []
fils=[i.split('.')[0]+'.png' for i in fil]
for img_path in fils:
    zaborchik.append(mpimg.imread(img_path))

f = plt.figure(figsize=(30,20))
columns = 3
for i, image in enumerate(zaborchik):
    f.add_subplot(int(len(zaborchik) / columns + 1), columns, i + 1)
    plt.imshow(image)
    plt.xticks([])
    plt.yticks([])

plt.show()

Во-первых, из-за того, что psi4 создает очень странные названия файлов, не получается автоматически расположить картинки в нужном порядке и с нужными названиями.

Во-вторых, мы видим, что psi4 смог нам построить такие же орбитали, как и наш рассчет функцией. Однако еще появились какие-то странные "кубические" орбитали.