В этом таске мы будем работать с волновыми функциями электрона.

Для начала немного теории.

(Если читателю уже известно или же не интересует, откуда берётся точная форма орбиталей электрона, эту часть можно пропустить.)

• Электроны — объекты квантовой природы. Это значит, что их поведение описывается волновыми функциями, квантовыми операторами и прочим математическим аппаратом квантовой теории.
  • Строго говоря, с точки зрения квантовой тероии все объекты имеют квантовую природу.
  • На практике, впрочем, объекты крупного масштаба обычно достаточно хорошо описываются классической физикой.
• В общем случае поведение квантовых объектов описывается уравнением Шрёдингера: $$ i\hbar\frac{\partial}{\partial t}\Psi(\mathbf{r},t)=\hat{H}\Psi(\mathbf{r},t) $$
  • Здесь $\Psi$ — волновая функция системы,
  • $\mathbf{r}$ — обобщённые координаты,
  • $i, \hbar$ — константы (мнимая единица и приведённая постоянная Планка),
  • $\hat{H}$ — Гамильтониан, также известный как оператор энергии системы.
    • Как и в классической механике, он состоит из кинетической и потенциальной энергии.
    • Которые, в свою очередь, являются функциями от положения и импульса объекта.
    • Однако в квантовой механике положение и импульс суть фукнции, не точки, а кинетическая и потенциальная энергии суть операторы, преобразующие их. Явно всё это записывается так: $$ \begin{align} \hat{H} & = \hat{T}+\hat{V} \\ \hat{V} & = V(\mathbf{r}, t) \\ \hat{T} & = \frac{\mathbf{\hat{p}\cdot\hat{p}}}{2m} = -\frac{\hbar^2}{2m}\nabla^2 \\ \end{align} $$
      • Здесь $V$ есть классическая функция потенциальной энергии,
      • $\mathbf{\hat{p}} = i\hbar\nabla$ есть оператор импульса.

• Полезно отметить следующее:
  • Это уравнение является линейным: для любых волновых функций $\Psi_1, \Psi_2$, являющихся его решениями, их линейная комбинация $\Psi_3=\alpha\Psi_1+\beta\Psi_2$ также будет решением.
  • Это значит, что множество решений уравнения образует линейное пространство.
  • Для статических систем (систем, у которых потенциальная энергия не зависит от времени) принято выбирать базисом этого пространства собственные векторы гамильтониана, и называть базисные векторы* чистыми квантовыми состояниями, а не-базисные векторы* — смешанными квантовыми состояниями.
  • * Если точнее, одному и тому же (чистому или смешанному) квантовому состоянию соответствует не конкретный вектор, а все векторы, отличающиеся друг от друга на константу: если функция $\Psi$ соответствует некоторому квантовому состоянию, то все функции $\lambda\Psi$ соответствуют ему же.

• Нас интересует стационарный электрон в атоме водорода. Это очень удобно.
  • Во-первых, удобна стационарность. Несложно показать (см., например), что если оператор потенциальной энергии системы не зависит от времени, то:
    • Волновую функцию чистых состояний можно разложить в зависящую и не зависящую от времени части: $\Psi(\mathbf{r},t) = \tau(t)\psi(\mathbf{r})$;
    • А нестационарное уравнение Шрёдингера можно (для чистых состояний) свести к стационарному: $$\hat{H}\psi(\mathbf{r}) = E\psi(\mathbf{r})$$
    • После чего получить все остальные решения (т.е. смешанные состояния) как линейные комбинации чистых.
  • Во-вторых, ядро атома водорода — одинокий протон — создаёт простой и понятный потенциал, из которого мы можем извлечь явное выражение для потенциальной энергии. Кинетическая энергия всегда устроена одинаково, и имея выражения для кинетической и потенциальной энергий, мы можем в явном виде записать Гамильтониан.
    • Строго говоря, атом водорода отличается от одинокого электрона в потенциале протона: в атоме водорода также есть собственно протон с собственной кинетической энергией.
      • К счастью, есть финт ушами под названием приведённая масса, позволяющий нам свести задачу о взаимном расположении двух тел к задаче о расположении единственного тела, и тем самым завернуть кинетическую энергию протона в кинетическую энергию электрона.
      • Вообще-то этот финт ушами совершенно не квантовый, но почему-то работает.
      • Возможно, он работает потому, что эффект его очень мал: приведённая масса электрона почти не отличается от просто массы электрона. Ну да ладно.

• Теперь мы можем, наконец, записать в явном виде уравнение, описывающее электрон в атоме водорода. В сферических координатах оно выглядит так: $$ \left(-\frac{\hbar^2}{2\mu}\nabla^2-\frac{e^2}{4\pi \epsilon_0r}\right)\psi(r,\theta,\varphi)=E\psi(r,\theta,\varphi) $$
  • Здесь $\mu=\frac{m_eM}{m_e+M}$ — приведённая масса,
  • $E$ — переменная; собственно решениями являются пары $(\psi, E)$,
  • $e$ и $\epsilon_0$ — константы (элементарный заряд и электрическая постоянная).
• Решать мы его не будем — это, пожалуй, выходит за рамки нашего краткого введения. За подробностями можно обратиться, к примеру, сюда.
• Нам достаточно знать его решение. Решение его — это семейство функций с тремя целочисленными параметрами $n,l,m$, имеющее общий вид: $$ % This is _not_ copypasted from wiki, I so swear. \psi_{n l m}(r, \theta, \varphi) = \sqrt{\left( \frac2{na^*_0} \right)^3 \frac{(n-l-1)!}{2n(n+l)!}} e^{-\rho/2} \rho^l L^{2l+1}_{n-l-1}(\rho) Y^m_l(\theta, \varphi) $$
  • Здесь $\rho = \frac{2r}{na^*_0}$,
  • $a^*_0 = \frac{4 \pi \epsilon_0 \hbar^2}{\mu e^2}$ — константа,
  • $L$ — обобщённые многочлены Лагерра,
  • $Y$ — сферические гармоники.
• Параметры $n, l, m$ также называются квантовыми числами и могут принимать следующие значения:
  • $n = 1, 2, 3, \dots$
  • $l = 0, 1, 2, \dots, n-1$
  • $m = -l, \dots, +l$

• Стоит иметь в виду, что это уравнение не описывает структуру орбиталей в совершенстве.
  • Например, если заменить классический оператор импульса $\hat{T} = \frac {\mathbf{\hat{p}}^2}{2m}$ на релятивистский $\hat{T} = \sqrt{\mathbf{\hat{p}}^2c^2+m^2c^4}-mc^2$ (и проделать ещё немного изысканий), можно увидеть тонкую структуру;
  • А если посчитать потенциальную энергию не в электрическом, а в электромагнитном поле, становится видна сверхтонкая.
• Но нам и так хватит.

Теперь, когда мы знаем, что именно мы моделируем, можно приступить к расчётам.

import numpy
import scipy.special
import scipy.misc
from scipy.misc import factorial as fact
import os

import npy2cube

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

ANGSTROM = 10**-11

def w(n,l,m,d,s=30):
    # Generates 30×30×30 space of coord-triplets, 
    # with d being boundaries of the space.
    
    # input-d is measured in angstroms. There's probably no vice in
    # measuring it in angsroms throughout the calculations?
    # d = d*ANGSTROM
    
    # This generates precisely 30 points (on every axis),
    # starting at -d and ending at -d. Between 30 points there are
    # 29 intervals, so step of the grid is equal to 2d/(30-1).
    x,y,z = numpy.mgrid[-d:d:s*1j,-d:d:s*1j,-d:d:s*1j]

    # Spherical transform: computes spherical coords from Cartesian.
    # Formulas are fairly self-explanatory.
    # Unfortunately, more intuitive explanations
    # require instructional diagrams, 
    # which are not text-comment-compatible.
    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)
    
    # Yes, this is the value of a_0.
    # For some reason, it was set to 1.0 before? Wtf?
    a0 = 0.529 # * ANGSTROM
    
    # Not-quite-radius
    # Useful for scaling
    rho = lambda x,y,z,n: 2*r(x,y,z)/n/a0*2**n
    
    # Here, we define/compute the wavefunction as function of position.
    # Firstly, radial dependence: root with factorials, power part,
    # exponential and laguerre parts. In fact, we'll skip factorials, 
    # as they are position-independent, and do not change the overall 
    # form of wavefunction.
    # 2**-n is scaling factor: we want our orbital displayed side-by-side, 
    # and care more about shapes than comparative sizes.
    R = lambda rho,n,l: \
        ((2/n/a0)**1.5 * 0 + 1) * \
        ((fact(n-l-1) / (2*n*fact(n+l)))**0.5 * 0 + 1) * \
        rho**l * numpy.exp(-rho/2) * \
        scipy.special.genlaguerre(n-l-1,2*l+1)(rho)
    # Secondly, angular dependence; namely, spherical harmonics.
    WF = lambda rho,theta,phi,n,l,m: \
        R(rho,n,l) * scipy.special.sph_harm(m,l,phi,theta)
    # Optionally thirdly, squared modulus of wavefunction.
    absWF = lambda rho,theta,phi,n,l,m: \
        numpy.absolute(WF(rho,theta,phi,n,l,m))**2
    
    return absWF(rho(x,y,z,n),theta(x,y,z),phi(x,y,z),n,l,m)

Переберём несколько значений квантовых чисел и сохраним модели функции для каждой тройки.

# Optimal boxes to fit hydrogen orbitals vary from orbital to orbital;
# lacking insight to decide box sizes individually, we arbitrarily designate
# 2-angstrom box as universal. (We actually experimented some and found 
# a good box size)
d = 2
_step = 2.*d/(30-1)

# Actually, npy2cube discards imaginary part of the wavefunction.
# Let's see how it goes, I guess?
#for i in range(30):
#    n = numpy.random.randint(1, 30+1)
#    l = numpy.random.randint(0, n)
#    m = numpy.random.randint(0, l+1)
for n in range(1, 4+1):
    for l in range(n):
        for m in range(l+1):
            grid = w(n,l,m,d,s=60)
            name = '%s_%s_%s' % (n,l,m)
            npy2cube.npy2cube(grid, 
                              start = (-d, -d, -d), 
                              step = (_step, _step, _step), 
                              cube_f = "cube/"+name+'.cube')

Подключаем PyMOL, в которым мы будем рассматривать полученные орбитали.

# pymol launching
import __main__
import time
from xmlrpclib import ServerProxy
#cmd = ServerProxy(uri="http://localhost:9123/RPC2")

#__main__.pymol_argv = [ 'pymol', '-c', '-s "pymol.log']

import pymol
pymol.finish_launching(['pymol', '-c'])
 
from pymol import cmd

from IPython.display import Image

# Borrowed from http://kodomo.cmm.msu.su/~sapsan/v2/terms/term8/PyMolPractice1.html.
default_image = './temp_pymol_img.png'
def prepare_image(width=300, height=300, sleep=2, filename=default_image):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png(filename)
    time.sleep(sleep)

Обнаруживается, что:

• При работе с pymol'ом на сервере можно загружать файлы .cube, но невозможно визуализировать объём, поскольку для этого нужно графическое окно;
• При работе с pymol'ом на персональной машине можно визуализировать объёмы, но невозможно загружать .cube: для этого нужно выставить правильные флаги при сборке pymol'а, а сборка — это сложно.

Конфуз, однако.
Можно, впрочем, загрузить файлы на сервере, сохранить их в более удобочитаемом формате, загрузить формат на персональной машине и потом отрисовать.
Так и сделаем.
На сервере исполняем это:

cmd.reinitialize()
for _name in os.listdir("./cube"):
    name = _name[:-5]
    cmd.load("cube/{}.cube".format(name), "cube")
    cmd.save("pse/{}.pse".format(name))
    time.sleep(.1)
    cmd.delete("all")

На другой машине исполняем это:

cmd.reinitialize()
# Откуда взялись эти числа?
# В принципе, их производит volume_ramp_new, 
# если дать ему на вход названия цветов.
# Но в данном случае мы взяли пример вывода volume_ramp_new,
# и долго и упорно обрабатывали его напильником.
cmd.volume_ramp_new('ramp007', [\
      0.0002,  1.00, 0.00, 0.00, 0.0005, \
      0.0010,  1.00, 1.00, 0.00, 0.005, \
      0.0025,  0.00, 1.00, 0.00, 0.01, \
      0.005,   0.00, 1.00, 1.00, 0.02, \
      0.01,    0.00, 0.00, 1.00, 0.01, \
      0.015,   1.00, 0.00, 1.00, 0.005, \
    ])
cmd.set("ray_volume")

for _name in os.listdir("./pse"):
    name = _name[:-4]
    cmd.load("pse/{}.pse".format(name))
    cmd.turn("x", -30)
    cmd.turn("y", 30)
    cmd.volume("vol", "cube", ramp="ramp007")
    cmd.ray()
    cmd.png("png/{}.png".format(name))
    time.sleep(1)
    cmd.delete("all")

Получаем красивые картинки.
Красивые картинки нарисованы не в правильном масштабе, а так, чтобы все орюбитали помещались примерно в одно и то же пространство.

Image("png_1/2_0_0.png")

Image("png_1/3_2_0.png")

Image("png_2/3_2_0.png")

Image("png_2/4_2_1.png")

По соседству считаем орбитали с помощью ORCA. Это весело.

• Во-первых, orca по-своему нумерует орбитали. Скорее всего, по уровню энергии.
• Во-вторых, считать орбитали водорода после пятой она отказывается.
  • Возможно, где-то есть опция, задающая число орбиталей, но найти её в семисотстраничном мануале — задача не из лёгких.
  • Посему мы считаем орбитали не водорода, а неона. (Неона с зарядодм 9 и одним электроном.) У него из коробки считается 14 орбиталей, нам для поверхностного взгляда хватит.
    • Собственно, для нашей задачи — построения высоких орбиталей водорода — orca не предназнчена. Она хорошо умеет считать «привычные» атомные орбитали и собирать из них орбитали молекулярные. n=3,4 для атома водорода не является ни тем, ни другим.

Входной файл:

%%writefile h.in

! 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);
 MO("H-5.cube",5,0);
 MO("H-6.cube",6,0);
 MO("H-7.cube",7,0);
 MO("H-8.cube",8,0);
 MO("H-9.cube",9,0);
 MO("H-10.cube",10,0);
 MO("H-11.cube",11,0);
 MO("H-12.cube",12,0);
 MO("H-13.cube",13,0);
end
* xyz 9 2
 Ne 0 0 0
*

Overwriting h.in

Код bash:

%%sh
rm -r orca
rm -r ./cube-orca
mkdir orca
mkdir cube-orca
mv h.in orca
cd orca
orca h.in > h.out
cd ..
mv orca/h.in ./
mv orca/*.cube ./cube-orca/

Повторяем пляски с бубном.

cmd.reinitialize()
for _name in os.listdir("./cube-orca"):
    name = _name[:-5]
    cmd.load("cube-orca/{}.cube".format(name), "cube")
    cmd.save("pse-orca/{}.pse".format(name))
    time.sleep(.1)
    cmd.delete("all")

Получаем новые красивые картинки.

Image("png-orca/H-1.png")

Image("png-orca/H-4.png")

Image("png-orca/H-8.png")

Image("png-orca/H-13.png")

Эти картинки, скажем так, отдалённо напоминают орбитали, которые у нас вышли раньше. Что, наверное, служит уроком в неточности симуляций.
С другой стороны, подсчёты с помощью orca драматически быстрее точных подсчётов. Что, в сущности, логично: нельзя получить и точность, и скорость, но можно обменять одно на другое.