Практикум 4. Вычисление атомных орбиталей водорода

In [1]:
import numpy as np
import scipy.special
import scipy.misc
In [2]:
from IPython.display import display, Image
In [3]:
def w(n,l,m,d):
    
    x,y,z = np.mgrid[-d:d:30j,-d:d:30j,-d:d:30j]
    
    # Ну тут понятно - переходим в радиальную систему координат
    # Радиус...
    r = lambda x,y,z: np.sqrt(x**2+y**2+z**2)
    # Угол тета...
    theta = lambda x,y,z: np.arccos(z/r(x,y,z))
    # Угол фи...
    phi = lambda x,y,z: np.arctan2(y,x)
    
    """
    В сто сорок солнц закат пылал,
    В июнь катилось лето, --
    Два синус тета эр квадрат
    Дэ эр дэ фи дэ тэта.
    """
    
    a0 = 1
    
    # Радиус функция
    R = lambda r,n,l: (2*r/n/a0)**l * np.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
    # Коэффициент нормировки
    b = lambda n,l: np.sqrt((scipy.special.factorial(n-l-1)/(2*n*scipy.special.factorial(n+l)))*(2/n/a0)**3)
    # Волновая функция
    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: np.absolute(WF(r,theta,phi,n,l,m))**2

    ### может тут чего то не хватает? - да, выводим не WF, а absWF
    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Копипаста функции Андрея Демкива

In [4]:
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 [5]:
for n in range(1, 4):
    d = 10 * n
    step = d / 10
    for l in range(n):
        for m in range(0,l+1,1):
            grid = np.sqrt(w(n, l, m, d))
            name = name='%s-%s-%s' % (n,l,m)
            npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')

Теперь сделаем скрипт для PyMol:

In [6]:
with open('open_cube.pml', 'w') as scr:
    scr.write('cmd.volume_ramp_new(\'ramp007\', \
                               [-0.015, 0.00, 1.00, 1.00, 0.50,  \
                                -0.01,  0.20, 0.80, 1.00, 0.20,  \
                                -0.005, 0.40, 0.60, 1.00, 0.00,  \
                                 0.005, 0.60, 0.40, 1.00, 0.00,  \
                                 0.01,  0.80, 0.20, 1.00, 0.20,  \
                                 0.015, 1.00, 0.00, 1.00, 0.50]) \n')
                    
    for n in range(1,4):
        for l in range(n):
            for m in range(0,l+1):
                name = '%s-%s-%s' % (n,l,m)
                cube_name = name + '.cube'
                volume_name = name + '_vol'
                scr.write(f'load {cube_name}, {name}\n')
                scr.write(f'volume {volume_name}, {name}\n')
                scr.write(f'volume_color {volume_name}, ramp007\n')
In [7]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [8]:
cmd.do("bg_color black")
cmd.load("open_cube.pml")
In [14]:
cmd.do('png 3-2-0.png')
In [15]:
Image('3-2-0.png')
Out[15]:

Выглядит неплохо... Но очень замыленно

Посмотрим в срезе:

In [19]:
cmd.do('png 3-2-0-cut.png')
In [21]:
Image('3-2-0-cut.png')
Out[21]:

Запустим ORCA

In [30]:
cmd.load('3-2-1.cube', 'H-13')
In [31]:
cmd.do('volume H-13_vol, H-13')
cmd.do('volume_color H-13_vol, ramp007')
In [34]:
cmd.do('png H-13.png')
In [35]:
Image('H-13.png')
Out[35]:
In [36]:
cmd.do('png H-13_cut.png')
In [37]:
Image('H-13_cut.png')
Out[37]:

Орбиталь сильно отличается.

А если мы посмотрим на первые 4 - отличий почти нет =(