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

Цель занятия: опираясь на уравнение построить электронную плотность в одно-электронном атоме и сравнить с известными программами.

Импортируем необходимые модули, функцию Андрея Демкива npy2cube прогоним через 2to3. Работает!

python D:\pymol\anaconda\Tools\scripts\2to3.py --output-dir=C:\Users\днс\npy3cube -W -n C:\Users\днс\npy2cube

In [1]:
from pymol import cmd, stored
import numpy
import scipy.special
import scipy.misc
from math import factorial
import matplotlib.pyplot as plt
from IPython.display import display,Image
In [7]:
from npy3cube import npy2cube

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

In [8]:
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 = 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 - Боровский радиус, радиус ближайшей к ядру орбиты электрона атома водорода в модели атома, предложенной Нильсом Бором
    a0 = 1
    
    #R - Радиальная часть волновой функции
    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)
    
    #Волновая функция - произведение coef*R*угловая функция (сферическая гармоника)
    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)

Cube-файлы со значениями волновых функций для первыx трёх уровней:

In [10]:
step = 1
#перебор квантовых чисел 0,1,2,3
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(grid,(-d,-d,-d),(step,step,step), name+'.cube')

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

Визуализируем орбитали в pymol:

In [11]:
orbit = open ('orbit.pml', 'w')

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

#записываем в файл инструкции по открыванию cube'ов и окрашиванию

1-0-0
2-0-0
2-1-0
2-1-1
3-0-0
3-1-0
3-1-1
3-2-0
3-2-1
3-2-2
In [15]:
import __main__

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

import pymol 
pymol.finish_launching()
from pymol import cmd, stored
cmd.load("orbit.pml")

PyMOL>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,     ])
PyMOL>load 1-0-0.cube, 1-0-0
PyMOL>volume 1-0-0_vol, 1-0-0
PyMOL>volume_color 1-0-0_vol, ramp007
PyMOL>load 2-0-0.cube, 2-0-0
PyMOL>volume 2-0-0_vol, 2-0-0
PyMOL>volume_color 2-0-0_vol, ramp007
PyMOL>load 2-1-0.cube, 2-1-0
PyMOL>volume 2-1-0_vol, 2-1-0
PyMOL>volume_color 2-1-0_vol, ramp007
PyMOL>load 2-1-1.cube, 2-1-1
PyMOL>volume 2-1-1_vol, 2-1-1
PyMOL>volume_color 2-1-1_vol, ramp007
PyMOL>load 3-0-0.cube, 3-0-0
PyMOL>volume 3-0-0_vol, 3-0-0
PyMOL>volume_color 3-0-0_vol, ramp007
PyMOL>load 3-1-0.cube, 3-1-0
PyMOL>volume 3-1-0_vol, 3-1-0
PyMOL>volume_color 3-1-0_vol, ramp007
PyMOL>load 3-1-1.cube, 3-1-1
PyMOL>volume 3-1-1_vol, 3-1-1
PyMOL>volume_color 3-1-1_vol, ramp007
PyMOL>load 3-2-0.cube, 3-2-0
PyMOL>volume 3-2-0_vol, 3-2-0
PyMOL>volume_color 3-2-0_vol, ramp007
PyMOL>load 3-2-1.cube, 3-2-1
PyMOL>volume 3-2-1_vol, 3-2-1
PyMOL>volume_color 3-2-1_vol, ramp007
PyMOL>load 3-2-2.cube, 3-2-2
PyMOL>volume 3-2-2_vol, 3-2-2
PyMOL>volume_color 3-2-2_vol, ramp007
In [17]:
#делаю скриншоты

orbitales = {0: 's', 1: 'p', 2: 'd'}

_, axarr = plt.subplots(5, 2, figsize=(10, 21.33))

row, col = 0, 0

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)
            f = plt.imread(name + '.png')
            axarr[row, col].axis('off')
            axarr[row, col].imshow(f)
            axarr[row, col].set_title(str(n)+orbitales[l], {'fontsize': 35}, pad=10)
            col += 1
            if col == 2:
                row += 1
                col = 0

plt.tight_layout()

Вроде похоже, сравним с расчётами, сделанными с помощью psi4.

In [ ]:
import psi4
In [7]:
g = '''
H
'''

m = psi4.geometry(g)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True, 'reference':'uhf'})
#https://en.wikipedia.org/wiki/Unrestricted_Hartree%E2%80%93Fock

energy, wave=psi4.energy('scf/cc-pvtz', return_wfn = True, molecule = m )
In [6]:
psi4.cubeprop(wave)
In [8]:
energy
Out[8]:
-0.49980981130184243
In [1]:
import __main__

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

import pymol 
pymol.finish_launching()
from pymol import cmd, stored
cmd.load("volume.pml")

У меня возникли проблемы с volume, то от орбиталей остаются куски, то окрашивается куб.. Выглядит эпично:

In [5]:
from IPython.display import display, Image
display(Image('epic.jpg'))

В итоге, всё выданное программой в порядке выдачи (названия у файлов вышли весьма замысловатые):

In [7]:
display(Image('all_together.jpg'))

Весьма похоже!