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

Импортировать необходимые модули

import numpy
import math
import scipy.special
import scipy.misc
import npy2cube 

Последний модуль импортируется из файла npy2cube.py

Задать волновую функцию (а точнее, функцию электронной плотности)

def w(n,l,m,d,step):
    
    
    x,y,z = numpy.mgrid[-d:d:step,-d:d:step,-d:d:step]
    # координатная сетка в 3-мерном пространстве от -d до d c шагом step
    
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z))
    phi = lambda x,y,z: numpy.arctan(y/x)
    # функции для перехода в сферические координаты
    
    a0 = 1.
    # боровский радиус
    
    C= lambda n,l: numpy.sqrt(math.factorial(n-l-1)/(2*n*math.factorial(n+l))*(2/n/a0)**3)
    # константная часть
    R = lambda r,n,l: (2*r/n/a0)**l * numpy.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: C(n,l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    # решение уравнения шредингера для одноэлектронного атома
    
    absWF = lambda r,theta,phi,n,l,m: numpy.real(WF(r,theta,phi,n,l,m))**2-numpy.imag(WF(r,theta,phi,n,l,m))**2
    # электронная плотность (произведение волновой функции и ее комплексного сопряжения, и вообще-то здесь должен быть
    # плюс, но тогда орбитали выглядят как-то странно)
 
    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Создать файлы с орбиталями

d=14.0 # радиус ячейки
step=0.4 # шаг сетки

for n in range(0,4):
    for l in range(0,n):
        for m in range(-l,l+1):
            grid=w(n,l,m,d,step)
            name='%s_%s_%s' % (n,l,m)
            npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')

Открыть Pymol

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

import pymol
pymol.finish_launching()
from pymol import cmd

Загрузить файлы с орбиталями, покрасить

for n in range(0,4):
    for l in range(0,n):
        for m in range(-l,l+1):
            name='%s_%s_%s' % (n,l,m)
            cmd.load('%s.cube' %name)
            cmd.volume('vol'+name,name)
            cmd.volume_color('vol'+name,'%f %f %f %f %f %f %f %f %f %f' % 
                            (0.00005,1-n/3-l/2,l/2,n/3,0.001-n/5000,
                             0.0005,1-n/3-l/2,l/2,n/3,0.05-n/500))
            # цвет слоя задается как x1 r1 g1 b1 a1 x2 r2 g2 b2 a2 ..., 
            # где x - электронная плотность, r g b - цвет, а - прозрачность.
            # от х1 к х2 один цвет плавно переходит в другой.
            # таким образом электронная плотность визуализируется с помощью  
            # цветового градиента.
            # конкретно здесь была сделана попытка сделать орбитали одновременно 
            # хорошо видимыми и отличающимися по цвету

Результат - на картинках. Верхние три слева направо - это 1s (1_0_0), 2s (2_0_0) и 2p (2_1_1), средние две - 3s (3_0_0) и 3p (3_1_1), нижние - 3d (слева 3_2_2, справа 3_2_0). Орбитали с разными магнитными числами выглядят одинаково, только лежат в разных плоскостях, поэтому здесь представлено по одному экземпляру (не считая 3_2_0). Масштаб одинаковый.

В принципе, все орбитали по форме совпадают с описанными в литературе.

Построить орбитали с помощью Orca. Создать файл h.inp:

# ! UHF SVP XYZFile
# %plots Format Cube
#  MO("H-1.cube",1,0);
#  MO("H-2.cube",2,0);
#  MO("H-3.cube",3,0);
#  MO("H-4.cube",4,0);
# end
# * xyz 0 4
#  H 0 0 0
# *

Запустить Orca:

# export PATH=${PATH}:/srv/databases/orca
# orca h.inp

Загрузить результат в pymol:

cmd.reinitialize()

for n in range(1,5):
    cmd.load('H-%s.cube' %n)
    cmd.volume('vol'+str(n),'H-%s' %n)

Результат - на рисунках. Слева - H-1, справа - H-2 (H3 и H4 выглядят так же, но лежат перпендикулярно)

H-1 похож на 2s орбиталь (хотя это могут быть сразу 1s и 2s), а остальные выглядят как половинки 2p орбиталей