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

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

In [2]:
import numpy
import scipy.special
import scipy.misc
from IPython.display import Image
import math
import npy2cube

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

In [2]:
def w(n,l,m,d):
    # n - principal quantum number
    # l - angular momentum quantum number
    # m - magnetic quantum number
    # d - step
    
    x,y,z = numpy.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    # создаем куб размером 30х30х30 в декартовой системе координат
    
    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)
    # переход к сферическим координатам (r, theta, phi)
    
    a0 = 1
    # боровский радиус - радиус ближайшей к ядру орбиты электрона атома водорода
    
    Normalization_factor = lambda n,l: numpy.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
    # Нормировочный коэффициент
    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: 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
    # Плотность вероятности волновой функции

    return WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
In [3]:
d = 100
step = 1.0

for n in range(0,4):
    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')

npy2cube.py:39: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f " %(float(grid[x, y, z])))
npy2cube.py:42: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f\n" %(float(grid[x, y, z])))

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

Загружаем PyMol

In [ ]:
import __main__

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

import pymol 
pymol.finish_launching()

from pymol import cmd

Скрипт для визуализации в PyMol

In [5]:
pml='''
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, 1.00, 0.00, 0.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, \
    ])
'''
with open('ramp007.pml', 'w') as f:
    f.write(pml)


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)
            with open('run_%s.pml' %name, 'w') as f:
                f.write('load ramp007.pml\n')
                f.write('load ' + name + '.cube\n')
                f.write('volume 1, ' + name + ', ramp007\n')
                f.write('draw\n')
                f.write('png ' + name + '.png')
                try:
                    Image(filename = name + '.png', width=500)
                except IOError:
                    continue
In [6]:
with open('color.pml', 'w') as color:

    color.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)
                color.write('load %s, %s\n' % (name+'.cube', name))
                color.write('volume %s, %s\n' % (name+'_vol', name))
                color.write('volume_color %s, ramp007\n' % (name+'_vol'))
In [7]:
cmd.load("color.pml")
In [8]:
cmd.zoom()
In [9]:
Image(filename='3-1-0.png')
Out[9]:
In [10]:
Image(filename='3-2-1.png')
Out[10]:
In [11]:
Image(filename='3-2-2.png')
Out[11]:
In [12]:
Image(filename='3-1-0.png')
Out[12]:
In [13]:
Image(filename='all.png')
Out[13]:

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

Теперь рассчитаем орбитали в программе orca. Для этого создаем текстовый файл h.inр

In [4]:
with open('h.inp', 'w') as f:
    f.write('''! 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
*''')

На kodomo запускаем orca с помощью нижеследующих команд

In [15]:
export PATH=${PATH}:/srv/databases/orca
orca h.inp

Рассчитываем орбитали в программе Orca и визуализируем все в PyMol

In [3]:
with open('orca.pml', 'w') as color:

    color.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,6):
        name='H-%s' % n
        color.write('load %s, %s\n' % (name+'.cube', name))
        color.write('volume %s, %s\n' % (name+'_vol', name))
        color.write('volume_color %s, ramp007\n' % (name+'_vol'))
In [4]:
from IPython.display import display,Image
cmd.load("orca.pml")
In [5]:
Image(filename='H-1(2s).png')
Out[5]:
In [6]:
Image(filename='H-2(2p).png')
Out[6]:
In [7]:
Image(filename='H-3(2p).png')
Out[7]:
In [8]:
Image(filename='H-4(2p).png')
Out[8]:

Результат работы orca - это четыре орбитали, одна аналогична орбитали с квантовыми числами 1-0-0, три другие - орбитали 2-1-0