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

Задание: https://kodomo.fbb.msu.ru/wiki/2016/8/MolSim/task4

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

Установим зависимости

 conda install -c anaconda numpy

Скачаем функцию npy2cube Андрея Демкива, работает в Python 2.x

 wget http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py
 mv npy2cube.py term8/pr4
In [13]:
from pymol import cmd, stored
from xmlrpc.client import ServerProxy
import numpy
import scipy.special
import scipy.misc
from math import factorial
import matplotlib.pyplot as plt

# Работает только под Python 2.x
import npy2cube

Мы рассчитываем ордитали abinitio, задаем волновую функцию.

In [1]:
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)      ##расстояние (переход от декартовых координат в сферические)
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z)) ##тета и фи - 2 угла, вместе с r описывающие положение точки в сферических координатах 
    phi = lambda x,y,z: numpy.arctan(y/x)
    
    a0 = 1
        ##корень с факториалами перед остальной частью волновой функции
    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)
    ##радиальная функция, в уравнении волновой функции это e^(-ro/2)*ro^l*обобщенный полином Лагерра
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    ##Волновая функция - произведение coef*R*угловая функция (сферическая гармоника)
    ##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-файлы с орбиталями.

In [9]:
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.

In [10]:
import pymol 
from pymol import cmd, stored
In [11]:
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
In [12]:
cmd.reinitialize()
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()

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 [9]:
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>load 2-0-0.cube, 2-0-0
PyMOL>volume 2-0-0_vol, 2-0-0
PyMOL>load 2-1-0.cube, 2-1-0
PyMOL>volume 2-1-0_vol, 2-1-0
PyMOL>load 2-1-1.cube, 2-1-1
PyMOL>volume 2-1-1_vol, 2-1-1
PyMOL>load 3-0-0.cube, 3-0-0
PyMOL>volume 3-0-0_vol, 3-0-0
PyMOL>load 3-1-0.cube, 3-1-0
PyMOL>volume 3-1-0_vol, 3-1-0
PyMOL>load 3-1-1.cube, 3-1-1
PyMOL>volume 3-1-1_vol, 3-1-1
PyMOL>load 3-2-0.cube, 3-2-0
PyMOL>volume 3-2-0_vol, 3-2-0
PyMOL>load 3-2-1.cube, 3-2-1
PyMOL>volume 3-2-1_vol, 3-2-1
PyMOL>load 3-2-2.cube, 3-2-2
PyMOL>volume 3-2-2_vol, 3-2-2
In [18]:
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()

Получили вот такое орбитали. Похоже, школьный учебник по химии не врет.

Теперь посмотрим, как выглядят орбитали в Orca.

In [14]:
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
*''')

Запускаем в терминале orca h.inp

Посмотрим результаты Orca в PyMol...... К сожалению, я вижу только черный экран.

Судя по картинкам, которые построили другие ребята, формы орбиталей в Orca очень похожи на построенные нами, но выглядят чуть более четко оерценными, возможно, они рассчитаны более точно.

In [ ]: