HW 4

Решение уравнения Шрёдингера для одноэлектронного атома:

$$\psi_{nlm}(r, \vartheta, \varphi) = \sqrt{{\left(\frac{2}{n a_0}\right)}^3\frac{(n-l-1)!}{2n(n+l)!}} \: e^{- \rho / 2} \: {\rho}^l \: L_{n-l-1}^{2l+1}(\rho) \: Y_l^m (\vartheta, \varphi)$$

$\rho = \frac{2r}{n a_0};$

При подстановке $\rho$ получаем:

$$\psi_{nlm}(r, \vartheta, \varphi) = \sqrt{{\left(\frac{2}{n a_0}\right)}^3\frac{(n-l-1)!}{2n(n+l)!}} \: e^{- r / na_0} \: {\left(\frac{2r}{na_0}\right)}^l \: L_{n-l-1}^{2l+1}{\left(\frac{2r}{na_0}\right)} \: Y_l^m (\vartheta, \varphi)$$

$L_{n-l-1}^{2l+1}(\rho)$ - обобщённый полином Лагерра степени $n-l-1;$

$Y_l^m (\vartheta, \varphi)$ - сферическая гармоника;

$n, \: l, \: m$ - квантовые числа:

  • $n$ - основное число, отражающее номер уровня $(1,2,3, \dots)$
  • $l$ - орбитальное число, отражающее тип орбитали $(0,1,2, \dots ,n-1)$
  • $m$ - магнитное число $(-l, \dots, +l)$
In [1]:
import os
import numpy as np
import scipy.special
import scipy.misc
import npy2cube as n2c
import pymol
from pymol import cmd
from IPython.display import display, Image

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

In [2]:
def w(n,l,m,d):

    # The meshgrid function is useful for creating coordinate arrays to vectorize function evaluations over a grid. 
    # nj - step length, integer part of its magnitude is interpreted as specifying 
    # the number of points to create between the -d and d values
    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.arctan(y/x)
    # боровский радиус
    a0 = 0.529
    
    # множитель с корнем
    sqrt_part = lambda n,l: np.sqrt((2/n/a0)**3 * scipy.special.factorial(n-l-1)/scipy.special.factorial(n+l)/2/n)
    # произведение всех множителей волнового уравнения, за исключением сферической гармоники
    R = lambda r,n,l: sqrt_part(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)
    # волновая функция: произведение определенной выше функции R и сферической гармоники
    WF = lambda r,theta,phi,n,l,m: 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

    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Сгенерируем .cube файлы для отображения:

In [3]:
d = 15
step = 1


for n in range(0, 3):
    for l in range(0, n):
        # посмотрим также, что будет при отрицательных значениях m
        for m in range(-l, l+1, 1):
            name='%s-%s-%s' % (n,l,m)
            n2c.npy2cube(w(n,l,m,d), (-d,-d,-d), (step,step,step), name + '.cube')

Скрипт для отображения в pymol:

In [4]:
with open('draw_mgrids.pml', 'w') as f:    
    for n in range(0, 3):
        for l in range(0, n):
            for m in range(-l, l+1, 1):
                name = '%s-%s-%s' % (n, l ,m)
                f.write('load /Users/marinavol/Desktop/hse/drug_design/hw4/%s, %s\n' % (name + '.cube', name)) 
                f.write('volume %s, %s\n' % (name + '_v', name))
                f.write('disable all\n')
                f.write('enable %s\n' % (name + '_v'))
                f.write('zoom\n')
                f.write('png %s.png, dpi=500\n' % (name))

Визуализация результатов для первых трех уровней:

In [5]:
for f in [f for f in os.listdir("/Users/marinavol/Desktop/hse/drug_design/hw4/") if f.endswith(".png")]:
    print f
    display(Image("/Users/marinavol/Desktop/hse/drug_design/hw4/" + f))

1-0-0.png

2-0-0.png

2-1--1.png

2-1-0.png

2-1-1.png

Создадим входной файл для пакета ORCA и сгенерим изображения.

In [6]:
# ssh -p 22025 shad@vsb.fbb.msu.ru
# mkdir ./hse/volosnikova
# cd ./hse/volosnikova
# nano h.inp
# orca h.inp
# scp -P 22025 shad@vsb.fbb.msu.ru:/home/shad/hse/volosnikova/*.cube ./Desktop/hse/drug_design/hw4/cubes/

Входной файл h.inp

In [7]:
# # ! UHF SVP XYZFile
# %plots Format Cube
#  MO("1-0-0.cube",1,0);
#  MO("2-0-0.cube",2,0);
#  MO("2-1-0.cube",3,0);
#  MO("2-1-1.cube",4,0);
# end
# * xyz 0 4
#  H 0 0 0
# *

Получили .cube файлы.

Визуализируем с помощью pymol.

In [8]:
with open('draw_mgrids_orca.pml', 'w') as f:    
    for n in range(0, 3):
        for l in range(0, n):
            for m in range(0, l+1, 1):
                name = '%s-%s-%s' % (n, l ,m)
                f.write('load /Users/marinavol/Desktop/hse/drug_design/hw4/cubes/%s, %s\n' % (name + '.cube', name)) 
                f.write('volume %s, %s\n' % (name + '_v', name))
                f.write('disable all\n')
                f.write('enable %s\n' % (name + '_v'))
                f.write('zoom\n')
                f.write('png %s_orca.png, dpi=500\n' % (name))
In [10]:
for f in [f for f in os.listdir("/Users/marinavol/Desktop/hse/drug_design/hw4/cubes/") if f.endswith(".png")]:
    print f
    display(Image("/Users/marinavol/Desktop/hse/drug_design/hw4/cubes/" + f))

1-0-0_orca.png

2-0-0_orca.png

2-1-0_orca.png

2-1-1_orca.png