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

Расчёт плотности

In [1]:
import numpy
import scipy.special
import scipy.misc
from IPython.display import display,Image
import math

Подгружаю функцию из файла npy2cube.py:

In [2]:
import npy2cube

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

In [3]:
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 - радиус, theta и phi - углы)  
    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  #боровский радиус
    #Допкоэффициент
    k = lambda n,l: numpy.sqrt((math.factorial(n-l-1)*(2/n/a0)**3)/(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)
    #Волновая функция, где scipy.special.sph_harm - сферическая гармоника
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta) * k(n, l)
    #Плотность вероятности волновой функции WF
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2
    
    return numpy.real(WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)) + numpy.imag( WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m))

Рассчитываю значения для первых трех уровней (на выход подаются cube файлы):

In [4]:
step = 1 #шаг
#Задаю цикл по перебору квантовых чисел
for n in range(0,4):
    d = 10*n
    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

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

In [5]:
#Записываю в файл
v = open ('task1.pml', 'w')
#Цифрами задаю цвет слоя в формате: значение электронной плотности, rgb, прозрачность
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)
            v.write('load %s, %s\n' % (name+'.cube', name))       # Загрузить cube файл
            v.write('volume %s, %s\n' % (name+'_vol', name))      # Отобразить электронную плотность
            v.write('volume_color %s, ramp007\n' % (name+'_vol')) #Покрасить ее разумным образом 
v.close()
In [6]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [7]:
from IPython.display import display,Image
cmd.load("task1.pml")

В PyMol я вручную сохранила изображения, некоторые из которых представлены здесь:

In [8]:
from IPython.display import display,Image
Image('orb.jpg', width=500)
Out[8]:

Расчет в программе Orca

В данном задании необходимо рассчитать орбитали в программе Orca. Для этого создаю текстовый файл h.inp:

In [9]:
h = open ('h.inp', 'w')
h.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
*''')
h.close()

На kodomo я запустила orca с помощью данных команд:

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

Визуализация выданных данных в PyMol:

In [10]:
#Записываю в файл
v = open ('orc.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(1,5):
            name='H-%s' % (n)
            v.write('load %s, %s\n' % (name+'.cube', name))       # Загрузить cube файл
            v.write('volume %s, %s\n' % (name+'_vol', name))      # Отобразить электронную плотность
            v.write('volume_color %s, ramp007\n' % (name+'_vol')) #Покрасить ее разумным образом 
v.close()
In [11]:
cmd.reinitialize()
cmd.load("orc.pml")

В PyMol я вручную сохранила изображения, которые представлены здесь:

In [12]:
from IPython.display import display,Image
Image('orc.jpg', width=500)
Out[12]: