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

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

Введем функцию npy2cube:

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

def w(n,l,m,d):
    #n,l,m – квантовые числа
    #d – шаг
    
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j] 
    #координаты в декартовой системе (30 на 30 на 30)
    
    #переход к сферическим координатам, r – радиус дистанции, thata, phi – углы относительно осей
    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
    
    #R – радиальная часть волновой функции, n, l – квантовые числа
    #scipy.special.genlaguerre - обобщённый полином Лагерра
    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)
    
    #добавочный коэффициент
    k = lambda n,l: numpy.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: numpy.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 numpy.real(result) + numpy.imag(result)

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

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

for n in range(1,4):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            d = 10 * n
            step=2*d/(30-1)
            grid=w(n, l, m, d)
            name='%s-%s-%s' % (n,l,m)
            numpy2cube.npy2cube(grid, (-d,)*3, (step,)*3, name+'.cube')

Откроем cube-файлы в PyMol и построим для них volume:

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

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

out = open ('volume_color.pml', 'w')

out.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'
            vname = name+'_vol'
            out.write('load %s, %s\n' % (lname, name))
            out.write('volume %s, %s\n' % (vname, name))
            out.write('volume_color %s, ramp007\n' % (vname))
out.close()

cmd.load("volume_color.pml")

cmd.reinitialize()

Покажу некоторые из полученных орбиталей:

from IPython.display import Image, display

2р-орбиталь

Image('C:/Users/Nastya/Desktop/mm/12.png')

3р-орбиталь

Image('C:/Users/Nastya/Desktop/mm/13.png')

3d-орбиталь

Image('C:/Users/Nastya/Desktop/mm/14.png')

Рассчитаем орбитали в программе Orca на kodomo:

with open('hydrogen.inp', 'w') as f:
    f.write('''! 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 2
H 0 0 0
* ''')

out = open ('volume_color_orca.pml', 'w')

out.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 file in range(1,14):
    name='%s-%s' % ('H', file)
    lname = name+'.cube'
    vname = name+'_vol'
    out.write('load %s, %s\n' % (lname, name))
    out.write('volume %s, %s\n' % (vname, name))
    out.write('volume_color %s, ramp007\n' % (vname))
out.close()

cmd.load("volume_color_orca.pml")

Получили изображение:

Image('C:/Users/Nastya/Desktop/mm/all.png')

Покажем некоторые отдельные орбитали:

Image('C:/Users/Nastya/Desktop/mm/15.png')

2р-орбиталь

Image('C:/Users/Nastya/Desktop/mm/16.png')

3-d орбиталь

Таким образом, результаты наших подсчетов и выдачи Orca в некоторых случаях совпадают. Но есть различия в форме орбитали, например, для 3d.