Практикум 4

import numpy as np 
import scipy.special
import scipy.misc
from numpy2cube import npy2cube
import time
import math
import matplotlib.pyplot as plt
import pymol
from pymol.cgo import *
from pymol import cmd,stored

Определим функцию w для расчета квадрата модуля $\psi$

def w(n,l,m,d):
    #n,l,m – квантовые числа
    #d – шаг
    
    #строим трехмерную сетку от -d до d 100 шагов в каждом измерении (потому что так орбитали будут детальнее)
    x,y,z = np.mgrid[-d:d:100j,-d:d:100j,-d:d:100j]
    
    #для каждой точки сетки считаем ее сферические координаты - r theta phi
    r = lambda x,y,z: np.sqrt(x**2+y**2+z**2)
    theta = lambda x,y,z: -np.arcsin(z/r(x,y,z)) + math.pi/2
    phi = lambda x,y,z: np.arctan2(y, x) + math.pi
    
    #Боровский радиус - наиболее вероятная дистанция между ядром и электроном
    a0 = 1
    
    #R – радиальная часть волновой функции, n, l – квантовые числа
    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)
    
    #нормировочный коэффициент
    k = lambda n,l: np.sqrt((2/n/a0)**3 * math.factorial(n-l-1)/(2*n*math.factorial(n+l)))
    
    #волновая функция, scipy.special.sph_harm – сферическая гармоника
    WF = lambda r,theta,phi,n,l,m: k(n, l) * R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    
    #плотность вероятности WF
    absWF = lambda r,theta,phi,n,l,m: np.absolute(WF(r,theta,phi,n,l,m))**2
    
    result = WF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)
    
    return np.real(result) + np.imag(result)

Рассмотрим первые 5 электронных уровней и создадим для них все возможные орбитали и запишем данные в cube файлы

# Зададим цикл по перебору квантовых чисел
for n in range(1,5):
    d = 15*n 
    step = 2*d/(100-1)
    for l in range(0,n):
        for m in range(-l,l+1,1):
            # для сохранения нужно задать координаты старта grid и шаг по каждому направлению
            grid = w(n, l, m, d)
            name = '{n}-{l}-{m}'.format(n=n,l=l,m=m)
            npy2cube(grid,
                     (- d, - d, - d),
                     (step, step, step),
                     name + '.cube')

Посмотрм на то как выглядят орбитали в pymol

import __main__
from IPython.display import Image
__main__.pymol_argv = [ 'pymol', '-x' ] # открыть pymol  без использования внешнего GUI
pymol.finish_launching()

cmd.reinitialize()
# первое число это значение функции ЭП при котором (и меньше ) плотность будет окрашена
# цветом по rgb (следующие три числа)  и прозрачность - alpha последнее число


cmd.volume_ramp_new('my_ramp', [\
     -0.3, 0.0, 0.0, 1.00, 0.6, \
     -0.01, 0.0, 0.0, 1.00, 0.6, \
     -0.005, 0.63, 0.31, 1.00, 0.20, \
     -0.001, 0.73, 0.67, 1.00, 0.0, \
      0.00, 1.00, 1.00, 1.00, 0.00, \
      0.001, 1.00, 0.79, 0.65, 0.0, \
      0.005, 1.00, 0.64, 0.22, 0.20, \
      0.01, 1.00, 0.00, 0.00, 0.6, \
      0.3, 1.00, 0.00, 0.00, 0.6                                
    ])


for n in range(1,5):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            name='%s-%s-%s' % (n,l,m)
            cmd.load(f"{name}.cube")
            cmd.volume(f"{name}_volume", name)
            cmd.volume_color(f"{name}_volume", "my_ramp")
cmd.hide("volume")

Покажем несколько орбиталей в срезе

cmd.hide("volume")
for n in range(1,5):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            cmd.reset()
            time.sleep(1)
            name=f"{n}-{l}-{m}_volume"
            d=-n*18
            cmd.show("volume",f"{name}")
            time.sleep(1)
            cmd.center(f"{name}")
            time.sleep(1)
            cmd.do(f"zoom {name}")
            if l==1:
                if m==0:
                    cmd.turn("y","90")
            if l>1:
                if abs(m)<=1:
                    cmd.turn("y","90")
                    ang=45*m
                    cmd.turn("x",f"{ang}")
            time.sleep(1)
            cmd.do(f"clip near, {d}")
            time.sleep(1)
            cmd.draw()
            cmd.png(f"{name}")
            time.sleep(3)
            cmd.hide("volume",f"{name}")

for n in range(1,5):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            name=f"{n}-{l}-{m}"
            display(name)
            display(Image(f"{name}_volume.png"))

'1-0-0'

'2-0-0'

'2-1--1'

'2-1-0'

'2-1-1'

'3-0-0'

'3-1--1'

'3-1-0'

'3-1-1'

'3-2--2'

'3-2--1'

'3-2-0'

'3-2-1'

'3-2-2'

'4-0-0'

'4-1--1'

'4-1-0'

'4-1-1'

'4-2--2'

'4-2--1'

'4-2-0'

'4-2-1'

'4-2-2'

'4-3--3'

'4-3--2'

'4-3--1'

'4-3-0'

'4-3-1'

'4-3-2'

'4-3-3'

Image("https://upload.wikimedia.org/wikipedia/commons/e/e7/Hydrogen_Density_Plots.png")

Рассчитаные при помощи python орбитали действительно похожи на то, что должно быть. То что между областями высокой электронной плотности не такие видимые провалы, как на картинке из интернета объясняется особенностью нашей цветовой схемы.

Посмотрим на них в объеме

cmd.hide("volume")
for n in range(1,5):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            cmd.reset()
            time.sleep(1)
            name=f"{n}-{l}-{m}_volume"
            d=-n*18
            cmd.show("volume",f"{name}")
            time.sleep(1)
            cmd.center(f"{name}")
            time.sleep(1)
            cmd.do(f"zoom {name}")
            if l==1:
                if m==0:
                    cmd.turn("y","90")
            if l>1:
                if abs(m)<=1:
                    cmd.turn("y","90")
                    ang=45*m
                    cmd.turn("x",f"{ang}")
            time.sleep(1)
            cmd.draw()
            cmd.png(f"{name}_3d")
            time.sleep(5)
            cmd.hide("volume",f"{name}")

---------------------------------------------------------------------------
CmdException                              Traceback (most recent call last)
<ipython-input-65-e596ce278cde> in <module>
      7             name=f"{n}-{l}-{m}_volume"
      8             d=-n*18
----> 9             cmd.show("volume",f"{name}")
     10             time.sleep(1)
     11             cmd.center(f"{name}")

~/anaconda3/envs/pymol3/lib/python3.8/site-packages/pymol/viewing.py in show(representation, selection, _self)
    569 
    570 '''
--> 571         return _showhide(representation, selection, 1, _self)
    572 
    573     def show_as(representation="wire", selection="", _self=cmd):

~/anaconda3/envs/pymol3/lib/python3.8/site-packages/pymol/viewing.py in _showhide(rep, selection, value, _self)
    529 
    530         with _self.lockcm:
--> 531             r = _cmd.showhide(_self._COb, str(selection), int(repn), value)
    532 
    533         if _self._raising(r,_self): raise QuietException

CmdException:  Error: Invalid selection name "1-0-0_volume".
1-0-0_volume<--

for n in range(1,4):
    for l in range(0,n):
        for m in range(-l,l+1,1):
            name=f"{n}-{l}-{m}"
            display(f"{n}_{l}_{m}")
            display(Image(f"{name}_volume_3d.png"))

'1_0_0'

'2_0_0'

'2_1_-1'

'2_1_0'

'2_1_1'

'3_0_0'

'3_1_-1'

'3_1_0'

'3_1_1'

'3_2_-2'

'3_2_-1'

'3_2_0'

'3_2_1'

'3_2_2'

Image("https://i.pinimg.com/originals/4b/fa/2b/4bfa2bee9532b12650b11d3845748aa8.jpg")

Вроде бы похожи

Сравним с выдачей orca

#Запишем файл для orca

with open('h.inp', 'w') as f:
    f.write('''! UHF aug-cc-pV5Z XYZFile
%plots Format Cube
MO("H0.cube",0,0);
MO("H1.cube",1,0);
MO("H2.cube",2,0);
MO("H3.cube",3,0);
MO("H4.cube",4,0);
MO("H5.cube",5,0);
MO("H6.cube",6,0);
MO("H7.cube",7,0);
MO("H8.cube",8,0);
MO("H9.cube",9,0);
MO("H10.cube",10,0);
MO("H11.cube",11,0);
MO("H12.cube",12,0);
MO("H13.cube",13,0);
MO("H14.cube",14,0);
MO("H15.cube",15,0);
MO("H16.cube",16,0);
MO("H17.cube",17,0);
MO("H18.cube",18,0);
MO("H19.cube",19,0);
MO("H20.cube",20,0);
end
* xyz 0 2
H 0 0 0
* ''')

Загрузим cube файлы в pymol и посмотрим на орбитали 2го энергетического уровня

cmd.reinitialize()
# первое число это значение функции ЭП при котором (и меньше ) плотность будет окрашена
# цветом по rgb (следующие три числа)  и прозрачность - alpha последнее число


cmd.volume_ramp_new('my_ramp', [
     -0.3, 0.0, 0.0, 1.00, 0.9, \
     -0.021, 0.0, 0.0, 1.00, 0.9, \
     -0.019, 0.63, 0.31, 1.00, 0.50, \
     -0.018, 0.73, 0.67, 1.00, 0.0, \
      0.00, 1.00, 1.00, 1.00, 0.00, \
      0.018, 1.00, 0.79, 0.65, 0.0, \
      0.019, 1.00, 0.64, 0.22, 0.50, \
      0.021, 1.00, 0.00, 0.00, 0.9,
      0.3, 1.00, 0.00, 0.00, 0.9                    
    ])


for n in range(0,21):
    cmd.load(f"H{n}.cube")
    cmd.volume(f"H{n}_volume", f"H{n}")
    cmd.volume_color(f"H{n}_volume", "my_ramp")
cmd.hide("volume")

Отметим что ввиду неколько других значений функции которые выдает orca пришлось переделать цветовую схему которая сама по себе исказит форму некоторых орбиталей.

Нарисуем орбитали

cmd.hide("volume")
for i in range(0,21):
            cmd.reset()
            time.sleep(1)
            name=f"H{i}_volume"
            cmd.show("volume",f"{name}")
            time.sleep(1)
            cmd.center(f"{name}")
            time.sleep(1)
            cmd.do(f"zoom {name}, complete=1")
            time.sleep(1)
            cmd.draw()
            cmd.png(f"{name}_orca")
            time.sleep(5)
            cmd.hide("volume",f"{name}")

---------------------------------------------------------------------------
CmdException                              Traceback (most recent call last)
<ipython-input-77-5d025e81e09f> in <module>
      7             time.sleep(1)
      8             name=f"H{i}_volume"
----> 9             cmd.show("volume",f"{name}")
     10             time.sleep(1)
     11             cmd.center(f"{name}")

~/anaconda3/envs/pymol3/lib/python3.8/site-packages/pymol/viewing.py in show(representation, selection, _self)
    569 
    570 '''
--> 571         return _showhide(representation, selection, 1, _self)
    572 
    573     def show_as(representation="wire", selection="", _self=cmd):

~/anaconda3/envs/pymol3/lib/python3.8/site-packages/pymol/viewing.py in _showhide(rep, selection, value, _self)
    529 
    530         with _self.lockcm:
--> 531             r = _cmd.showhide(_self._COb, str(selection), int(repn), value)
    532 
    533         if _self._raising(r,_self): raise QuietException

CmdException:  Error: Invalid selection name "H21_volume".
H21_volume<--

for i in range(0,21):
            name=f"H{i}_volume"
            display(f"{name}")
            display(Image(f"{name}_orca.png"))

'H0_volume'

'H1_volume'

'H2_volume'

'H3_volume'

'H4_volume'

'H5_volume'

'H6_volume'

'H7_volume'

'H8_volume'

'H9_volume'

'H10_volume'

'H11_volume'

'H12_volume'

'H13_volume'

'H14_volume'

'H15_volume'

'H16_volume'

'H17_volume'

'H18_volume'

'H19_volume'

'H20_volume'

Конечно стоит понимать, что цветовая схема несколько исказила сами орбитали (например кольца в 4f орбиталях - не замкнуты), но мы видим что топологически орбитали рассчитанные orca и нами - идентичны. Но конечно они различаются размером, значениями функции и распределением значений функции (иначе не пришлось бы так сильно и непропорционально переделывать цветовую схему чтобы все отобразилось нормально) и немного формой - орбитали orca более округлые. Очевидно что orca считает немного по-другому чем мы.