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

In [1]:
import numpy
import scipy.special
import scipy.misc
from IPython.display import Image
In [3]:
import npy2cube
In [42]:
def npy2cube(grid, start, step, cube_f):
    '''
    PARAMETERS:
        grid : numpy array
            3-dimentional array, containing grid data
        start : tuple
            format: (x, y, z), coordinates of cube start point
        step: tuple
            format: (x, y, z), step size on 3 axes 
        cube_f: string
             name of output .cube file
    
    RETURNS:
        void    
    '''
    
    cube_string = ""
    bohr_to_angs = 0.529177210859 #const
    start = list(map(lambda x: x/bohr_to_angs, start))
    step = list(map(lambda x: x/bohr_to_angs, step))

    with open(cube_f, "w") as cube_file:
        ###HEADER###
        cube_file.write(" CPMD CUBE FILE.\nOUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n")
        cube_file.write("    1 %f %f %f\n" %(start[0], start[1], start[2]))
        cube_file.write("  %i    %f   0.000000     0.000000\n" %(grid.shape[0], step[0]))
        cube_file.write("  %i    0.000000    %f    0.000000\n" %(grid.shape[1], step[1]))
        cube_file.write("  %i    0.000000    0.000000    %f\n" %(grid.shape[2], step[2]))
        cube_file.write("    1    0.000000    %f %f %f\n" %(start[0], start[1], start[2]))

        ###DATA###
        i = 0
        for x in range(grid.shape[0]):
            for y in range(grid.shape[1]):
                for z in range(grid.shape[2]):
                    if i < 5:
                        cube_file.write("%f " %(float(grid[x, y, z])))
                        i += 1
                    elif i == 5:
                        cube_file.write("%f\n" %(float(grid[x, y, z])))
                        i = 0
    return 0
In [16]:
def w(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) # длина вектора из (x, y, z) = радиальное расстояние
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z)) # зенитный угол
    phi = lambda x,y,z: numpy.arctan(y/x) # азимутальный угол
    
    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)
     # wave function
    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)

    out = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    return numpy.real(out) + numpy.imag(out)
In [43]:
# цикл по перебору квантовых чисел
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)  # трехмерный массив из 30*30*30 элементов
            name='%s-%s-%s' % (n,l,m)
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            npy2cube(grid,(-d, -d, -d),(step,step,step),name+'.cube')

<ipython-input-42-7eb3ff519de2>:37: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f " %(float(grid[x, y, z])))
<ipython-input-42-7eb3ff519de2>:40: ComplexWarning: Casting complex values to real discards the imaginary part
  cube_file.write("%f\n" %(float(grid[x, y, z])))

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

In [44]:
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()
In [45]:
from multiprocessing import Process

def launch_pymol():
    !pymol  -R &> /dev/null
Process(target=launch_pymol).start()
In [ ]:
import xmlrpc.client as xmlrpclib
cmd = xmlrpclib.ServerProxy('http://localhost:9123')
cmd.reinitialize()
cmd.load("v.pml")
In [49]:
Image('4_0.png', width = 700)
Out[49]:

Рассчитаем орбитали в ORCA. Создадим файл h.inp, запустим ORCA на kodomo.

In [52]:
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")
In [53]:
Image('4_1.png', width = 700)
Out[53]:

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