Вычисление атомных орбиталей водорода

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

Используем функцию из скрипта npy2cube.py от Андрея Демкива

def npy2cube(grid, start, step, cube_f):

    bohr_to_angs = 0.529177210859  # const
    start = [x / bohr_to_angs for x in start]
    step = [x / bohr_to_angs for x in step]

    with open(cube_f, 'w') as cube_file:
        cube_file.write('CPMD CUBE FILE.\n'
                        'OUTER 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]))


        i = 0
        for grid_item in np.nditer(grid, order='C'):
            if i < 5:
                cube_file.write('%f ' % (float(grid_item)))
                i += 1
            elif i == 5:
                cube_file.write('%f\n' % (float(grid_item)))
                i = 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)

Считаем значения для первых 3х уровней c помощью перебора квантовых чисел и запись в файл .cube.

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):
            grid = w(n, l, m, d)
            name = f'{n}_{l}_{m}'
            npy2cube(grid, (-d,-d,-d), (step,step,step), f'{name}.cube')

Визуализируем в Pymol

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

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

cmd.reinitialize()
# #задаём цветовую схему
cmd.volume_ramp_new('ramp007', [-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(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')
            cmd.hide()
            cmd.show('volume', f'{name}-vol')
            cmd.reset()
            cmd.turn('x', 90)
            cmd.clip('slab', 0)
            time.sleep(0.3)
            cmd.do(f'png clip_{name}, width={480}, height={270}')
            
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')
            cmd.hide()
            cmd.show('volume', f'{name}-vol')
            cmd.reset()
            cmd.turn('x', 45)
            cmd.turn('y', 45)
            cmd.turn('z', 45)
            time.sleep(0.3)
            cmd.do(f'png {name}, width={480}, height={270}')

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'clip_{i}.png' for i in numbers]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

# поверхности

img_paths = [f'{i}.png' for i in numbers]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

insides = '''! 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
* '''

with open('hydrogen.inp', 'w') as f:
    f.write(insides)

С выдачей Orca повторяем то же, что и с выдачей npy2cube.

cmd.reinitialize()
cmd.volume_ramp_new('ramp007', [-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')
    cmd.volume(f'H-{n}-vol', f'H-{n}', ramp='ramp007')
    cmd.hide()
    cmd.show('volume', f'H-{n}-vol')
    cmd.reset()
    cmd.turn('x', 90)
    cmd.clip('slab', 0)
    time.sleep(0.5)
    cmd.do(f'png clip_orca_{n}, width={480}, height={270}')
    
for n in range(14):
    cmd.load(f'H-{n}.cube')
    cmd.volume(f'H-{n}-vol', f'H-{n}', ramp='ramp007')
    cmd.hide()
    cmd.show('volume', f'H-{n}-vol')
    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={480}, height={270}')

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

# срезы Orca
img_paths = [f'clip_orca_{i}.png' for i in range(14)]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

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

# поверхности Orca
img_paths = [f'orca_{i}.png' for i in range(14)]
images = []
for img_path in img_paths:
    images.append(mpimg.imread(img_path))

orbitals = ['s1','s2','p2','p2','p2','s3','p3','p3','p3','d3','d3','d3','d3','d3']
f = plt.figure(figsize=(20,10))
columns = 4
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'{orbitals[i]}')

plt.show()

Выдачb получились похожими, только форма орбиталей более сложная у Orca. Скорее всего, Orca поточнее.