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



Небольшое введение


Как нам любит рассказывать всякая научная фантастика и псевдонаучная ерунда, типа фильма «Секрет», законы микромира, очень сильно отличаются от привычных нам, классических.
В мире квантовой механики правит всем вероятность, задаваемая волновой функцией $\psi$ (интересующиеся деталями могут заглянуть, например, в пост «Мюонный катализ с точки зрения квантовой химии. Часть I: обычный водород vs. мюонный водород»).
Из вероятностных свойств квантмеха растут ноги всяких забавных вещей, типа котов Шрёдингера, принципов неопределённости Гейзенберга и неравенств Белла.

Но все эти картинки со всякими орбиталями электрона не давали ответа на вопрос «как же электрон летит в пространстве». Чтобы прояснить эту ситуацию, физики потратили много времени, но так и не справились с этим. Зато Дэвид Бом (известный многим по эффекту Ааронова — Бома) окончательно создал один из формализмов квантовой механики (имени себя), в котором всё-таки есть траектории, по которым движется квантовая частица. И, в отличие от фейнмановских интегралов по траекториям, эта траектория у каждой частицы ровно одна. Это свойство принципиально позволяет отследить движение частиц, и сравнить движение классических и квантовых частиц, чем мы и займёмся в этой заметке.
не только формализм
Собственно, сам формализм никому особо не интересен, но из этого формализма можно построить одну из интерпретаций квантовой механики, которую, из-за кажущийся простоты классической механики любят некоторые фрики (не многие, ибо въехать в это дело не очень просто).
Эту (как и другие) интерпретацию мы обсуждать не будем.


Классические и квантовые траектории


Мы будем рассматривать достаточно скучную систему: один электрон в поле нескольких протонов. Об этой системе, а также о классической и квантовой механике можно почитать в первой и второй частях постов «Мюонный катализ с точки зрения квантовой химии».

Классическая задача о движении частицы в некотором потенциале даётся вторым законом Ньютона:

$m \ddot{x} = F$


где m — масса частицы, x — координата, F — сила, действующая на частицу, а $\ddot{x} = \frac{d^2 x}{dt^2}$ — вторая производная координаты частицы по времени, или ускорение. Если в системе действуют только потенциальные силы, то силу можно выразить через новую сущность, потенциальную энергию V как

$F = - \frac{dV}{dx}$


В нашем случае, электрона в поле нескольких протонов,

где электрон взаимодействует с каждым из протонов по закону Кулона

$V(R) = - k e^2/R$

, где k — коэффициент, равный 1 в атомных единицах, e — заряд электрона, а R — расстояние от электрона до протона.
В этом случае суммарный потенциал, действующий на электрон будет равен

$V = \sum_{n=1}^N V_n(R_n) = - \sum_{n=1}^N \frac{k e^2}{R_n} \ ,$


где индекс n нумерует протоны (всего протонов N штук), а Rn — расстояние от электрона до n-го протона.

Численно решить диффур, который представляет из себя второй закон Ньютона, задача халтурная, главное задать начальные положение и скорость. Если электрон будет лететь слишком быстро, то он вырвется за пределы притяжения протона(ов) и улетит на бесконечность, а если энергии будет чуть-чуть, то он будет вечно бултыхаться туды-сюды в поле одного из ядер, никогда не посещая другие.
Лучистое трение
Если бы мы учли лучистое трение, которое возникает из-за того, что при движении с ускорением, электрон будет часть своей энергии движения отдавать электромагнитному полю, излучая его куда-то, то электрон в итоге скатился бы на ядро за какое-то время.

Так что происходит в классике, мы знаем.

А что же будет в бомовской динамике?
В этом случае частица тоже будет двигаться по второму закону Ньютона с потенциалом $V = V_\mathrm{C} + V_\mathrm{Q}$, где $V_\mathrm{C}$ — классический потенциал из обычного закона Ньютона, который в нашем случае имеет вид, данный выше.
Т.е. помимо классического потенциала, на неё будет действовать ещё и другая сущность: квантовый потенциал $V_\mathrm{Q}$, имеющий (в 1D случае) вид

$V_\mathrm{Q} = - \frac{\hbar^2}{2mA} \frac{ d^2A}{ dx^2}$


где A — это амплитуда (модуль) волновой функции $A = |\psi|$ ($\psi = A \exp(i \varphi)$, где $\varphi$ — фаза волновой функции).
Значит, чтобы получить уравнение движения квантовой частицы, нам всё равно нужно знать что-то о волновой функции.
Про скрытые параметры
Формализм Бома является теорией со скрытыми параметрами. Но поскольку скрытый параметр (волновая функция) — нелокален, то результаты вычислений в этом формализме всё равно удовлетворяют вышеупомянутым неравенствам Белла.


В случае одного протона нам известно (см., например, тут) точное выражение волновой функции электрона в основном (1s) состоянии [в атомных единицах]:

$\psi(R) = \exp(-R)$


О нормировке и единицах
В формуле для квантового потенциала, нормировка числителя сократится со знаменателем, поэтому мы о ней париться и не будем.
В аргументе экспоненты, на самом деле, стоит не $R$, а $R/a_0$, где $a_0$ — это боровский радиус (0.529 Å). Но, поскольку мы юзаем атомные единицы, где $a_0 = 1$, эту единицу длины мы можем позволить себе не писать. Поподробнее об этом можно почитать тут.

В случае нескольких протонов, в рамках метода молекулярных орбиталей как комбинаций атомных орбиталей (МО ЛКАО, см. тут) основное состояние с достаточной степенью точности будет даваться суммой 1s-орбиталей каждого из атомов:

$\psi \approx \sum_{n=1}^N \psi_n(R_n) = \sum_{n=1}^N \exp(-R_n)$


Теперь, чтобы узнать квантовый потенциал, нужно всего-то воспользоваться этим выражением.
Ну<s>д</s>жные вычисления
Функция $\psi$ как сумма 1s орбиталей — действительная, поэтому $A = \psi$.
Поскольку электрон может двигаться в трёх измерениях, нужно одномерную производную $A_{xx}'' = \frac{d^2 A}{dx^2}$ заменить на её трёхмерное обобщение: $\Delta A = A_{xx}'' + A_{yy}'' + A_{zz}''$. Оператор $\Delta$ можно представить как квадрат оператора набла: $\Delta = \nabla ^2 $. Также можно представить расстояние $R_n$ как $R_n = \sqrt{\mathbf{R}_n^2}$, где $\mathbf{R}_n$ — радиус-вектор электрона относительно n-го протона.

Тогда

$\Delta A = \nabla^2 \psi = \sum_{n=1}^N \nabla^2 \psi_n(R_n)$


Первая производная считается легко:

$\nabla \psi_n(R_n) = \nabla \exp(-R_n) = \exp(-R_n) \cdot (-1) \cdot \frac{1}{2 \underbrace{\sqrt{\mathbf{R}_n^2}}_{R_n}} \cdot 2 \mathbf{R}_n = -\exp(-R_n) \cdot \frac{\mathbf{R}_n}{R_n}$


Вторая же производная уже несколько сложнее:

$\nabla (\nabla \exp(-R_n)) = -\frac{\mathbf{R}_n}{R_n}\nabla \exp(-R_n) -\exp(-R_n) \nabla \frac{\mathbf{R}_n}{R_n} = \exp(-R_n) - \frac{2\exp(-R_n) }{R_n}$


где $-\frac{\mathbf{R}_n}{R_n}\nabla \exp(-R_n) = \exp(-R_n) \cdot \underbrace{ \left(-\frac{\mathbf{R}_n}{R_n} \right)^2}_{1} = \exp(-R_n)$ и
$ -\exp(-R_n) \nabla \frac{\mathbf{R}_n}{R_n} = -\exp(-R_n) \cdot \left( \frac{\overbrace{\nabla \mathbf{R}_n}^{3}}{R_n} - \frac{2\mathbf{R}_n^2}{2 R_n^3} \right) = - \frac{2\exp(-R_n)}{R_n}$.
В итоге остаётся:
$\Delta \psi = \overbrace{\sum_{n=1}^{N} \exp(-R_n)}^{\psi} - \sum_{n=1}^{N} \frac{2\exp(-R_n)}{R_n}$
Поделив всё на $\psi=A$ и домножив на $-\frac{\hbar^2}{2m}$
получим
$V_\mathrm{Q} = -\frac{\hbar^2}{2m} \left( 1 - \sum_{n=1}^N \frac{2 \exp(-R_n)}{R_n}\right)$
Единица при дифференцировании для получения силы исчезнет, поэтому можно смело оставить только второе слагаемое.

В итоге наш квантовый потенциал мы можем записать как

$V_\mathrm{Q} \approx \frac{\hbar^2}{m} \sum_{n=1}^N \frac{ \exp(-R_n)}{R_n}$


и с этим выражением мы можем уже гонять бомовскую динамику электрона в поле многих протонов.

Реализация


Для этого всего безобразия был написан код на питоне, он доступен тут:
Код на питоне
from math import *
import numpy as np

cutoff=5.0e-4
Quantum=True

def dist(r1,r2):
    return np.dot((r1-r2), (r1-r2))

def Vc(r, r0):
    if dist(r, r0)>=cutoff:
        return -1.0/dist(r, r0)
    else:
        return -1.0/cutoff


    
rH=[]
#h1
#rH.append(np.array([ 0.0, 0.0, 0.0]))

#h2
rH.append(np.array([-1.0, 0.0, 0.0]))
rH.append(np.array([+1.0, 0.0, 0.0]))




def Vat(r):
    V=0.0
    for rh in rH:
        V+=Vc(r,rh)
    return V

def PsiA(r):
    psi=0.0
    for rh in rH:
        if dist(r, rh)<1.0e3:
            psi+=np.exp(-dist(r, rh))
    return psi

def Vq(r):
    vq=0.0
    for rh in rH:
        if dist(r, rh)>=cutoff:
            vq-=2.0*np.exp(-dist(r, rh))/dist(r, rh)
        else:
            vq-=2.0*np.exp(-cutoff)/cutoff    

    vq*=(-0.5) # -0.5*hbar**2/me
    return vq

def GradF(F, r):
    grad=np.zeros(3)
    dx=0.1
    for i in range(0,3):
        dr=np.zeros(3)
        dr[i]=dx
        #print(dr)
        #print(F(r+dr)-F(r-dr))
        grad[i]+=(F(r+dr)-F(r-dr))/(2.*dx)
    return grad


    

dt=0.001
tmax=2.0e1
DR=1.0
dx=0.001

MaxR=10.0

t=0.0

cent=np.zeros(3)

Ntrj=30
m=1.0

def GenRvBox(DX):
    return np.random.uniform(-DX,+DX,3)

def GenRvSph(DX):
    r=np.random.uniform(0.0,DX)
    phi=np.random.uniform(0.0,2.0*np.pi)
    theta=np.random.uniform(0.0,np.pi)
    x=r*np.sin(theta)*np.cos(phi)
    y=r*np.sin(theta)*np.sin(phi)
    z=r*np.cos(theta)
    return np.array([x,y,z])

for ntrj in range(0,Ntrj):
    if Quantum:
        outf=open("bmd_%05i" % (ntrj) + ".trj", "w")
    else:
        outf=open("cmd_%05i" % (ntrj) + ".trj", "w")
    
    nat=np.random.randint(0,len(rH))
    r=rH[nat]+GenRvSph(DR)
    rprev=r+GenRvBox(dx)
    outf.write("%15.10f %15.10f %15.10f\n" % tuple(r))
    t=0.0
    while t<=tmax and dist(r,cent)<=MaxR:
        F= -GradF(Vat, r)  
        if Quantum:
            F-= GradF(Vq, r)
        rnew= 2.*r - rprev + (F/m)*dt**2
        rprev=r
        r=rnew
        outf.write("%15.10f %15.10f %15.10f\n" % tuple(r))
    
        t+=dt
    outf.close()

exit()


Мы же только обсудим несколько моментов.
Интегрируется второй закон Ньютона при помощи алгоритма Верле:

$x(t + \Delta t) = 2 x(t) - x(t - \Delta t) + \frac{F(t)}{m} \Delta t^2$


Начальное положение генерируется путём случайного выбора одного из протонов, вокруг него случайно выбирается направле��ие (используя сферические координаты). Чтобы задать начальную скорость, нужно задать ещё одно, предыдущее положение. Оно выбирается ещё при помощи одного маленького случайного вектора.

Включая/выключая квантовый потенциал, мы переходим в квантовый/классический режимы движения.

Ну а дальше, можно построить красивые картинки, используя Gnuplot, для атома водорода

и для молекулы H2+


Как видно, классические траектории (верхние, синие) или очень локализованы, или, если слишком быстро заставить двигаться электрон, убегают от ядер. В квантовом же случае (нижние, розовые), квантовый потенциал позволяет электронам гулять существенно дальше от ядра, а в случае молекулы H2+, позволяет бегать от одного протона к другому, что является косвенной визуализацией химических связей.

Пару слов о построении картинок: чтобы создать эффект неона, каждая траектория рисуется много раз, от тонкой белой, до толстой чёрной, через тени интересующего цвета. Для удобства выбора такой палитры можно, например, воспользоваться сайтом https://www.color-hex.com/
Пример скрипта приведён ниже.
Скрипт для Gnuplot
unset key
set xyplane relative 0

unset box

set view map

set size ratio -1

unset border
unset xtics
unset ytics

set terminal pngcairo size 2160,4096 backgr rgb "black"
set output "tmp.png"

yshift=-5.0
maxiC=29
maxiQ=29
splot \
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 30.0 lc rgb "#030d19" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 18.0 lc rgb "#071b33" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 17.0 lc rgb "#0a294c" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 16.0 lc rgb "#0e3766" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 15.0 lc rgb "#11457f" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 14.0 lc rgb "#155399" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 13.0 lc rgb "#1861b2" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 12.0 lc rgb "#1c6fcc" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 11.0 lc rgb "#1f7de5" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 10.0 lc rgb "#238bff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 9.0 lc rgb "#3896ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 8. lc rgb "#4ea2ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 7. lc rgb "#65adff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 6. lc rgb "#7bb9ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 5. lc rgb "#91c5ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 4. lc rgb "#a7d0ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 3. lc rgb "#bddcff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 2. lc rgb "#d3e7ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 1. lc rgb "#e9f3ff" not,\
for [i=0:maxiC] sprintf("cmd_%05i.trj", i) w l lw 0.5 lc rgb "#ffffff" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 30.0 lc rgb "#190613" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 18.0 lc rgb "#330c27" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 17.0 lc rgb "#4c123b" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 16.0 lc rgb "#66184f" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 15.0 lc rgb "#7f1e63" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 14.0 lc rgb "#992476" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 13.0 lc rgb "#b22a8a" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 12.0 lc rgb "#cc309e" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 11.0 lc rgb "#e536b2" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 10.0 lc rgb "#ff3dc6" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 9.0 lc rgb "#ff50cb" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 8. lc rgb "#ff63d1" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 7. lc rgb "#ff77d7" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 6. lc rgb "#ff8adc" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 5. lc rgb "#ff9ee2" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 4. lc rgb "#ffb1e8" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 3. lc rgb "#ffc4ed" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 2. lc rgb "#ffd8f3" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 1. lc rgb "#ffebf9" not,\
for [i=0:maxiQ] sprintf("bmd_%05i.trj", i) u 1:($2+yshift):3 w l lw 0.5 lc rgb "#ffffff" not



Заключение


Бомовские траектории, хоть и сложны в понимании и вычислении, позволяют рисовать красивые картинки, которые показывают, насколько квантовая механика веселее и богаче классической.

Если есть комментарии, вопросы, предложения: пишите. :)