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

Загрузим все необходимые модули:

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

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

def w(n,l,m,d):
    #n,l,m - квантовые числа
    #d - шаг
    
    #координаты в Декартовой системе
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    
    #переходим в сферическую систему координат
    #(удобнее, тк наша система имеет сферическую симметрию)
    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2)
    theta = lambda x,y,z: -numpy.arcsin(z/r(x,y,z)) + math.pi/2
    phi = lambda x,y,z: numpy.arctan2(y, x) + math.pi
    
    a0 = 1
    
    #радиальная часть волновой функции 
    #scipy.special.genlaguerre - обобщённый полином Лагерра степени n-l-1 ; ρ = 2*r/n/a0
    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)
    
    #дополнительный коэффициент
    b = lambda n,l: numpy.sqrt((math.factorial(n-l-1)/(2*n*math.factorial(n+l)))*(2/n/a0)**3) 
    
    #волновая функция 
    #scipy.special.sph_harm - сферическая гармоника
    WF = lambda r,theta,phi,n,l,m: b(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
    #плотность вероятности волновой функции
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
    
    out = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    return numpy.real(out) + numpy.imag(out)

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

# Зададим цикл по перебору квантовых чисел
for n in range(0,4):
    d = 10 * n
    step = d / 15
    for l in range(0,n):
        for m in range(0,l+1,1):
            grid= w(n,l,m,d)
            name='%s-%s-%s' % (n,l,m)
            npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')

Создадим скрипт для pymol для визуализации .cube файлов:

v = open ('v.pml', 'w')
v.write('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, \
    ])\n')

for n in range(0,4):
    for l in range(0,n):
        for m in range(0,l+1,1):
            name='%s-%s-%s' % (n,l,m)
            lname = name+'.cube'
            volname = name+'_vol'
            v.write('load %s, %s\n' % (lname, name))
            v.write('volume %s, %s\n' % (volname, name))
            v.write('volume_color %s, ramp007\n' % (volname))
v.close()

from multiprocessing import Process

def launch_pymol():
    !pymol  -R &> /dev/null
Process(target=launch_pymol).start()

import xmlrpc.client as xmlrpclib
cmd = xmlrpclib.ServerProxy('http://localhost:9123')

cmd.reinitialize()
cmd.load("v.pml")

from IPython.display import Image
Image('v.png', width = 700)

 #Для s-орбиталей изображены также срезы

Рассчитаем орбитали в 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 на kodomo:

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

Визуализация выдачи в PyMOL:

orca = open ('orca.pml', 'w')
orca.write('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, \
    ])\n')
for n in range(1,5):
        name='H-%s' % n
        volname = name+'_vol'
        orca.write('load %s, %s\n' % (name+'.cube', name))
        orca.write('volume %s, %s\n' % ( volname, name))
        orca.write('volume_color %s, ramp007\n' % ( volname))
orca.close()

cmd.reinitialize()
cmd.load("orca.pml")

Image('orca.png', width = 700)

 В целом результаты выдачи ORCA схожи с полученными выше. Форма орбителей одинакова, однако соотношение их размеров в ORCA другое. По форме орбители H1 затруднительно определить ее тип, но если посмотреть на срез, то становится понятным, что это 2-0-0 (2s). H1, H3, H4 соответствуют 2p орбиталям.