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

In [ ]:
import numpy
import scipy.special
import scipy.misc
In [8]:
import npy2cube
In [14]:
from IPython.display import Image

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

In [9]:
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))
    phi = lambda x,y,z: numpy.arctan(y/x)
    
    # Радиус Бора
    a0 = 1.
    
    # Нормировочный коэффициент
    k = numpy.sqrt((2/n/a0) ** 3 * scipy.special.factorial(n-l-1) / (2*n*scipy.special.factorial(n+l)))
    # Радиальная функция
    R = lambda r,n,l: k * (2*r/n/a0)**l * numpy.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0)
    # Волновая функция
    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 absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
In [10]:
# Зададим цикл по перебору квантовых чисел

for n in range(0,4):
    d = n * 10
    step = d / 10
    for l in range(0,n):
        for m in range(0,l+1,1):
            grid= numpy.sqrt(w(n,l,m,d))
            name='%s-%s-%s' % (n,l,m)
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            npy2cube.npy2cube(grid,(-d,-d,-d),(step,step,step),name+'.cube')

В результате получились файлы в формате CUBE. Посмотрим на них в 3DMol. Вот, например, один из них:

In [23]:
import py3Dmol
view = py3Dmol.view()
alpha = open('1-0-0.cube','r').read()
view.addVolumetricData(alpha, "cube", {'isoval': 0.04, 'color': "red", 'opacity': 0.5})
view.zoomTo()
view.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Но лучше все-таки откроем их в PyMol

In [1]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd, stored
In [2]:
cmd.volume_ramp_new('ramp007', [\
      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, \
    ])

# С помощью panel можно делать другие, интересные окраски, но по-моему этой схемы для данного задания вполне достаточно
In [51]:
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)
            cmd.load(f'{name}.cube')
            cmd.volume(f'{name}-vol', name, 'ramp007')
            cmd.turn('x', 90)
            cmd.clip('slab', 5)
            cmd.do(f'png class4_{name}.png')
            cmd.reset()
            cmd.hide()

Посмотрим теперь, что получилось:

In [68]:
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)
            display(name)
            display(Image(filename=f'./class4_pngs/class4_{name}.png', width=200))

'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'

Кажется, все выглядит более-менее нормально

Теперь попробуем рассчитать орбитали в программе Psi4

In [2]:
import psi4
import numpy as np
psi4.core.set_output_file('output.dat')
In [9]:
# Геометрия молекулы:
g = '''
H
'''
In [10]:
# Рассчитаем энергию b-орбитали
m = psi4.geometry(g)
psi4.set_options({"maxiter": 200, "fail_on_maxiter" :  True, 'reference':'uhf'})
# Выбрали неограниченный метод Хартри-Фока

ener, wave=psi4.energy('scf/cc-pvtz', return_wfn = True, molecule = m )
In [11]:
ener
Out[11]:
-0.49980981130184254
In [12]:
psi4.cubeprop(wave)

Получился набор файлов cube. Они имеют очень странные замысловатые названия. Поэтому довольно сложно определить, что есть что. Так или иначе, вручную в PyMol отредактируем результаты и посмотрим, что вышло (и думаю, что обойдемся теперь без разрезов)

In [15]:
print('1s-орбиталь')
Image(filename=f'./class4_pngs/psi_a1.png', width=300)

1s-орбиталь
Out[15]:
In [18]:
print('2s-орбиталь')
Image(filename=f'./class4_pngs/psi_a2_new.png', width=300)

2s-орбиталь
Out[18]:
In [17]:
print('2p-орбиталь')
Image(filename=f'./class4_pngs/psi_a5.png', width=300)

2p-орбиталь
Out[17]:
In [19]:
print('3s-орбиталь')
Image(filename=f'./class4_pngs/psi_a6.png', width=300)

3s-орбиталь
Out[19]:
In [20]:
print('3p-орбиталь')
Image(filename=f'./class4_pngs/psi_a12.png', width=300)

3p-орбиталь
Out[20]:
In [21]:
print('3d-орбиталь')
Image(filename=f'./class4_pngs/psi_a9.png', width=300)

3d-орбиталь
Out[21]:
In [22]:
print('3dz^2-орбиталь')
Image(filename=f'./class4_pngs/psi_a11.png', width=300)

3dz^2-орбиталь
Out[22]:

Изображения, полученные нами ранее, по большому счету вроде бы совпадают с теми, что мы получили сейчас с помощью Psi4