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

Функция npy2cube была загружена отсюда.

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

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

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

Теперь визуализиурем полученные файлы в Pymol. Очень хотелось сделать всё автоматически и сразу, но почему-то не вышло: команды по отдельности в pymol работают, а скриптом нет. Pymol зависает и не отвечает на запросы. Вот такой вот скрипт использовался.

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

import pymol
pymol.finish_launching()
from pymol import cmd,stored
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', 45)
            cmd.turn('y', 45)
            cmd.turn('z', 45)
            cmd.do(f'png {name}, width=600, height=400')

Поскольку автоматически сделать картинки не удалось, а самому их все вручную сохранять не хочется, посмотрим по одному представителю s, p и d орбитали.

Image('s2.png')

Image('p3.png')

Image('d3.png')

Входной файл для программы Orca получили вот так:

with open('orca.inp', 'w') as f:
    f.write('''! UHF QZVPP XYZFile
%plots Format Cube
''')
    for i in range(0,14):
        f.write(f'MO("H-{i}.cube",{i},0);\n')
    f.write('''end
* xyz 0 4
H 0 0 0
* ''')

Затем запустили вот эти команды в терминале и пролучили модели орбиталей

export PATH=${PATH}:/srv/databases/orca

orca h.inp

Скриптом ниже тоже не получилось автоматически визуализировать по той же причине.

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

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

for n in range(14):
    cmd.load(f'H-{n}.cube')
    time.sleep(0.4)
    cmd.volume(f'H-{n}-volume', f'H-{n}', ramp='colors')
    time.sleep(0.4)
    cmd.hide()
    cmd.show('volume', f'H-{n}-volume')
    cmd.reset()
    cmd.turn('x', 45)
    cmd.turn('y', 45)
    cmd.turn('z', 45)
    time.sleep(0.4)
    cmd.do(f'png orca_{n}, width=500, height=500')

Попросил товарища Антона Иззи визуализировать орбитали, полученные оркой, на его компьютере. У него запустилось. Так что можно посмотреть на все орбитали.

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

orcas = [f'orcas/biba_H-{i}.png' for i in range(0,14)]
pics = []
for img_path in img_paths:
    pics.append(mpimg.imread(img_path))

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

plt.show()

К сожалению, из-за каких-то неполадок с pymol не удалось получить все орбитали первым методом. Те орбитали, которые сравнил в pymol вручную, оказались очень схожими с теми, которые получены программой Orca.