Домашнее задание 4¶
Выполнил: Ким Адамейко, группа мАДБМ16
In [4]:
! wget -nv -nc http://kodomo.fbb.msu.ru/~golovin/ipynb/npy2cube.py
from npy2cube import *
In [186]:
import time
import __main__
import os
import numpy as np
import scipy.special
import scipy.misc
from IPython.display import Image, display, HTML
from __future__ import print_function, division
from xmlrpclib import ServerProxy
cmd = ServerProxy(uri="http://localhost:9123/RPC2")
cmd.do("cd " + os.getcwd())
def img_row(img_data): # tuples of (filename, title)
img_files, titles = zip(*img_data)
display(HTML("<table>" +
('<tr>' + ''.join(['<th style="text-align: center;">' + title + '</th>'
for title in titles]) + '</tr>' if titles else '') +
'<tr>' + ''.join(['<td><img src="' + img + '.png?v=' + str(int(time.time())) + '"/></td>'
for img in img_files])+'</tr></table>'))
Часть 1. Волновая функция для атома водорода¶
В полярных координатах $(r,\vartheta, \varphi)$:
$${\large \psi_{n\,l\,m}(r,\vartheta, \varphi)} = \sqrt{\left(\dfrac{2}{n a_{0}}\right)^{\!\!3} \dfrac{(n-l-1)!}{2n\cdot(n+l)!}}\ {\large e^{-\frac{p}{2}}\ p^{l}\ L^{2l+1}_{n-l-1}(p)\ Y^{m}_{l}(\vartheta, \varphi),\quad p = \dfrac{2r}{na_{0}}}$$ где:
- $L^{2l+1}_{n-l-1}(p) -$ обобщенный полином Лагерра степени $n-l-1;$
- $ Y^{m}_{l}(\vartheta, \varphi) -$ сферическая гармоника степени $l$ и порядка $m$;
- $a_0 = \dfrac{4 \pi \varepsilon_0 \hbar^2}{\mu e^2} -$ приведенный радиус Бора;
- $\varepsilon_0 - $ электрическая постоянная;
- $e - $ элементарный заряд;
- $n, l, m\ - $ основные квантовые числа
- $n$ – основное число $(1,2,3, \dots)$
- $l$ – орбитальное число $(0,1,2,\dots n-1)$
- $m$ – магнитное число $(-l, \dots, +l)$
Часть 2. Ручной расчёт орбиталей¶
Определим необходимые функции для расчётов и воспользуемся функцией() npy2cube из соответствующего файла, предоставленного нам заранее:
In [229]:
def w(n, l, m, d):
# Неразреженный 3D массив точек, величина компл.ч. (30j) задаёт число точек между границами (такой синтаксис)
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 = 1.0
# Радиальная часть волновой функции
R = lambda r,n,l: (2.0*r/n/a0)**l * np.exp(-r/n/a0) * scipy.special.genlaguerre(n-l-1,2*l+1)(2.0*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: np.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)
In [236]:
d = 30
step= float(2.*d/29)
cdir = './cube/'
! mkdir -p {cdir}
qns = []
for n in range(1, 4):
for l in range(0, n):
for m in range(0, l+1):
data = w(n, l, m, d)
qns.append('%d-%d-%d' % (n, l, m))
filename = '%s%d-%d-%d.cube' % (cdir, n, l, m)
npy2cube(data, (-d,-d, -d), (step,step,step), filename) # парам: коорд. "самого отриц.точки" куба и шаг
!ls {cdir}
In [237]:
!head {cdir}2-1-1.cube
Сделаем необходимые настройки и осуществим рендеринг объемов с помощью Pymol (на локальном компьютере)
In [269]:
# настройка цветового градиента, для правильного отображения он должен быть симметричен относительно нуля
# числа: x-координата, RGB-код цвета, прозрачность
cmd.do(
"""
cmd.volume_ramp_new('ramp', [\
-0.05, 0.36, 0.70, 1.00, 0.0, \
-0.03, 0.50, 0.00, 1.00, 0.05, \
-0.02, 0.92, 0.00, 0.80, 0.00, \
0.02, 0.92, 0.00, 0.80, 0.00, \
0.03, 0.50, 0.00, 1.00, 0.05, \
0.05, 0.36, 0.70, 1.00, 0.0, \
])
""")
In [270]:
imgdir = './images/manualcomp/'
!mkdir -p {imgdir}
for qn in qns:
cmd.reinitialize()
cmd.do('load ./cube/{0}.cube, {0}'.format(qn))
cmd.do('volume {0}_vol, {0}'.format(qn))
cmd.do('volume_color {0}_vol, ramp'.format(qn))
time.sleep(0.5)
if qn in ['2-1-0', '3-1-0', '3-2-0']: # некоторые картинки нужно повернуть для лучшего отражения ситуации
cmd.turn('y', 90)
cmd.turn('x', -45)
cmd.do('png {0}{1}.png, width={2}, height={3}'.format(imgdir, qn, 400, 300))
time.sleep(0.5)
cmd.delete('all')
In [271]:
for i in range(0, len(qns), 2):
img_row(zip(map(lambda x: imgdir + x, qns[i:i+2]), qns[i:i+2]))
| 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 |
|---|---|
 |  |
Часть 3. Расчёт орбиталей в ORCA¶
Следующие две ячейки запускались на сервере. Затем результаты расчётов в виде cube-файлов архивировались и загружались на локальный компьютер, где снова отрисовывались в Pymol. К сожалению, удалось получить корректные изображения только для первых четырёх конфигураций квантовых чисел.
In [ ]:
%%writefile h.in
! UHF SVP XYZFile
%plots Format Cube
MO("./cube-orca/1-0-0.cube",0,0);
MO("./cube-orca/2-0-0.cube",1,0);
MO("./cube-orca/2-1-0.cube",2,0);
MO("./cube-orca/2-1-1.cube",3,0);
end
* xyz 0 4
H 0 0 0
*
In [ ]:
import os
curpath = os.getcwd()
! mkdir -p './cube-orca'
! export PATH=${PATH}:/{curpath}
! orca h.in > h.out
! tail h.out -n 2
! tar -zcf cube-orca.tar.gz cube-orca
Следующие ячейки запускались на локальном компьютере
In [212]:
cdir2 = './cube-orca/'
! mkdir -p {cdir2}
! tar -zxf cube-orca.tar.gz -C {cdir2} --strip-components=1
In [282]:
imgdir2 = './images/orca/'
!mkdir -p {imgdir2}
for qn in qns[:4]:
cmd.reinitialize()
cmd.do('load ./cube-orca/{0}.cube, {0}'.format(qn))
cmd.do('volume {0}_vol, {0}'.format(qn))
cmd.do('volume_color {0}_vol, ramp'.format(qn))
time.sleep(0.5)
if qn in ['2-1-0', '2-1-1']:
cmd.turn('y', 90)
cmd.turn('x', 45)
cmd.do('png {0}{1}.png, width={2}, height={3}'.format(imgdir2, qn, 400, 300))
time.sleep(0.5)
cmd.delete('all')
In [283]:
for i in range(0, 4, 2):
img_row(zip(map(lambda x: imgdir2 + x, qns[i:i+2]), qns[i:i+2]))
| 1-0-0 | 2-0-0 |
|---|---|
 |  |
| 2-1-0 | 2-1-1 |
|---|---|
 |  |
Наблюдается некоторое сходство с картинками, полученными с помощью "ручных" расчётов. Также заметна разница в масштабах. Возможно, это можно настроить с помощью каких-либо параметров ORCA, но я не успел разобраться в этом.