В этой статье я подробно опишу как численно смоделировать простейшую квантовую систему — систему с двумя состояниями. Ценность этой симуляции не только в ней самой, но также в той базе, которой она является для любой квантовомеханической симуляции.

Система с двумя состояниями — простейшая квантомеханическая система и основа квантовой механики. Частными случаями такой системы являются кубит, спин 1/2, поляризация фотона, кот Шрёдингера и многое другое. Теория таких систем служит базой и для более сложных систем вроде спиновых цепочек. Поэтому очень естественно начать моделирование именно с системы с двумя состояниями.

Максимально подробно и на весьма элементарном уровне теория такой системы изложена в 7 главе третьего тома (выпуски 8 и 9) "Фейнмановских лекций по физике". Там приведены уравнения, которым подчиняется эволюция состояния такой системы, равно как и их аналитическое решение. В данной же статье я покажу, как решить эти же уравнения численно. Хотя это может выглядеть стрельбой из пушки по воробьям, но база, которую мы заложим, позволит в дальнейшем перейти к более сложным системам, где аналитическое решение затруднительно.

Использовать для расчётов я буду (современный) Фортран, так что можно рассматривать эту статью и как своеобразное введение в программирование на Фортране. Это по-своему красивый язык, паритета с которым в ряде аспектов не имеет ни один другой язык программирования.

Уравнение Шрёдингера для системы с двумя состояниями.

Квантовомеханическое состояние системы интересующего нас вида есть просто пара комплексных чисел,, где— это амплитуда вероятности обнаружить систему в-м базисном состоянии. (Например, для кота Шрёдингера базисными состояниями могут быть "кот жив" и "кот мёртв".) Удобно их объединить в вектор-столбец

$\psi = \begin{pmatrix} C_1 \\ C_2 \end{pmatrix}.$

Поскольку полная вероятность оказаться в каком-либо состоянии должна быть единица, то должно выполняться свойство нормировки

$|C_1|^2+|C_2|^2=1. \qquad(*)$

Эволюция (почти) любой квантовомеханической системы задаётся гамильтонианом, то есть комплексной матрицей, в нашем случае размера 2 на 2:

$H=\begin{pmatrix} H_{11} & H_{12} \\ H_{21} & H_{22} \end{pmatrix},$

которая должна быть эрмитовой, что сводится в нашем случае к тому, что $H_{11}$ и $H_{22}$ должны быть вещественными, а $H_{12}$ и $H_{21}$ — комплексно-сопряжёнными друг к другу.

Тогда главное квантомеханическое уравнение — уравнение Шрёдингера — выглядит так:

$i \hbar \frac {d \psi} {dt} = H \psi,$

где-- время,-- мнимая единица, $\hbar$ -- постоянная Планка, а справа стоит матричное произведениеи $\psi$ .

Именно это уравнение мы и будем решать.

Фортран: подготовка окружения и проекта

Для работы потребуется установить gfortran, Fortran Package Manager, DISLIN и LAPACK с BLAS. Для Debian это делается следующим образом. gfortran, LAPACK и BLAS ставим из штатного репозитория:

# apt install gfortran liblapack-dev openblas-dev

Затем качаем FPM в виде единого *.F90 файла и компилируем его:

$ gfortran -o fpm fpm-*.F90

Полученный бинарник для удобства стоит поместить куда-нибудь в $PATH (в идеале куда-то вроде $HOME/bin).

DISLIN берём с официального сайта в виде deb-пакета и устанавливаем:

# apt install dislin-*.linux.i586_64.deb

Наконец мы готовы создать проект:

$ fpm new tss
<...>
$ cd tss
$ fpm run
+ mkdir -p build/dependencies
[  0%]                        tss.f90
[ 25%]                        tss.f90  done.
[ 25%]                       libtss.a
[ 50%]                       libtss.a  done.
[ 50%]                       main.f90
[ 75%]                       main.f90  done.
[ 75%]                            tss
[100%]                            tss  done.
[100%] Project compiled successfully.
 Hello, tss!

Как решать уравнение Шрёдингера

Уравнение Шрёдингера похоже на уравнение типа уравнения теплопроводности, но оно очень хитрое. Если попробовать применить что-то вроде обычного метода Эйлера, то получим такую схему:

$\psi(t_{n+1})=\left(1-\frac {i \Delta t} \hbar H \right) \psi(t_n)$

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

$\psi(t_{n+1})=\left(1+\frac {i \Delta t} \hbar H \right)^{-1} \psi(t_n)$

Такая схема стабильна, но она не сохраняет нормировку, что тоже нехорошо. Оказывается можно взять половинку от одной схемы и половинку от другой и получить стабильную схему, сохраняющую нормировку:

$\psi(t_{n+1})=\left(1+\frac {i \Delta t} {2 \hbar} H \right)^{-1} \left(1-\frac {i \Delta t} {2 \hbar} H \right) \psi(t_n)$

Это выражение можно упростить:

$\left(1+\frac {i \Delta t} {2 \hbar} H \right)^{-1} \left(1-\frac {i \Delta t} {2 \hbar} H \right) = \left(1+\frac {i \Delta t} {2 \hbar} H \right)^{-1} \left(-1-\frac {i \Delta t} {2 \hbar} H + 2 \right) =$ $= -1 + 2\left(1+\frac {i \Delta t} {2 \hbar} H \right)^{-1},$

и схема получается следующей.

Находим $\chi$ такое, что $\left(1+\frac {i \Delta t} {2 \hbar} H \right)\chi = \psi(t_n)$ .
Находим $\psi(t_{n+1}) = 2 \chi - \psi(t_n)$ .

Кодируем расчёт

Открываем src/tss.f90 и приступаем к написанию кода. Объявим несколько констант. Глобальные константы в Фортран объявляются в первой части модуля, то есть до строчки contains. Будем пользоваться системой единиц, в которой $\hbar = 1$ . Начнём со вспомогательной константы, амплитуды перехода между двумя состояниями:

real, parameter :: A = 1.0

Тут всё просто: real — это вещественный тип одинарной точности, parameter — так в Фортране называются и обозначаются константы, A — имя константы (используем обозначение из фейнмановских лекций), 1.0 — значение константы.

Используя эту константу, создадим константный же массив размера 2 на 2, который будет содержать гамильтониан:

complex, dimension(2, 2), parameter :: H = reshape([0.0, -A, -A, 0.0], [2, 2])

То, что написано слева от знака присваивания =, опять же просто: complex — комплексное одинарной точности, dimension задаёт размерность массива, parameter — константа, H — имя константы. Справа посложнее. Дело в том, что в Фортране нет готовой возможности записывать значение многомерного массива, только одномерного — перечислением элементов в квадратных скобках. Поэтому чтобы задать двухмерный массив мы сначала задаём одномерный массив с нужным числом элементов ( $2 \cdot 2 = 4$ ) и пользуемся встроенной функцией reshape, которая преобразует массив в массив такого же размера, но иной формы — в данном случае в двумерный массив два на два. Важно: в Фортране многомерные массивы хранятся не в построчном формате, а в поколоночном. То есть порядок перечисления элементовздесь такой: $H_{11}$ , $H_{21}$ , $H_{12}$ , $H_{22}$ .

Нам также понадобится единичная матрица, поэтому добавим представляющий её массив:

complex, dimension(2, 2), parameter :: I = reshape([1.0, 0.0, 0.0, 1.0], [2, 2])

Начальное состояние:

complex, dimension(2), parameter :: Psi_0 = [1.0, 0.0]

Наконец, объявим целочисленную константу количества шагов, которые мы будем вычислять и веществ��нную константу $\Delta t$ :

integer, parameter :: M = 100
real, parameter :: DeltaT = 0.1

Теперь у нас есть всё, чтобы вычислить $\left(1+\frac {i \Delta t} {2 \hbar} H \right)$ . Обозначим эту величину через. Поскольку она зависит только от констант, её тоже можно сделать константой:

  complex, dimension(2, 2), parameter :: U = I + (0.0, 0.5) * DeltaT * H

Обратите внимание, что там, где в формуле была единица, мы использовали единичную матрицу. Именно это подразумевается в таких формулах: все числовые константы там, где это необходимо, неявно умножаются на единичную матрицу. В Фортране же операции с константами — поэлементные, то есть если просто прибавить константу к матрице, она добавится ко всем элементам, что совсем не то же самое. Поэтому при переводе физических формул на Фортран следует не забывать вставлять в нужные места единичную матрицу.

Следующий этап — вычисление $\chi$ . Напишем заготовку соответствующей функции:

function calc_Chi(Psi) result(Chi)
  complex, dimension(2), intent(in) :: Psi
  complex, dimension(2) :: Chi
end function

Функции и процедуры в Фортране необходимо размещать во второй части модуля, то есть после строки contains. Здесь calc_Chi — это название функции, Psi — имя единственного аргумента. Конструкция result(Chi) задаёт имя, которое будет использовано внутри функции чтобы вернуть результат. Далее идёт описание типов аргументов и результата. Тут нам уже всё знакомо, кроме intent(in), который указывает, что соотвествующий параметр является именно входным, а не выходным (out) или двунаправленным (inout).

Для вычисления $\chi$ необходимо решить матричное линейное уравнение. Это делается с помощью библиотеки LAPACK. Открываем в документации раздел Driver Routines: Linear Equations и читаем. Понимаем, что нам нужна функция CGESV: она предназначена для работы с комплексными числами одинарной точности и для матрицы общего вида — как раз то, что нам нужно. Читаем документацию этой функции и приступаем к заполнению аргументов.

Первый параметр — количество уравнений, то есть строк матрицы коэффициентов, то есть в нашем случае 2. Второй параметр — сколько уравнений за раз мы хотим решить. Нам нужно решить только одно, значит 1. Дальше идёт матрица коэффициентов, у нас это U. Но просто так передавать U нельзя, так как это константа, а соответствующий параметр является и входным и выходным (то есть inout). Поэтому копируем её в массив-переменную A и передаём A.

Потом опять количество строк матрицы коэффициентов (на случай если матрица содержит лишние строки), то есть опять 2. Затем мало интересный нам выходной параметр под названием IPIV, какой-то целочисленный массив размера — в нашем случае — 2. Наконец, правая часть. У нас это Psi, но тут нужен двумерный массив, хотя и с единичным размером второго измерения. Кроме того, в этот же массив будут записано решение, то есть Chi, поэтому делаем так. Объявляем вспомогательный массив B:

complex, dimension(2, 1) :: B

Затем нам надо загнать в него Psi. Фортран позволяет сделать это одной командой, одним присваиванием:

B(:, 1) = Psi

Единичка тут потому, что по умолчанию нумерация в массивах начинается с 1. Двоеточие означает "все возможные значения". Таким образом, вся эта запись в левой части присваивания переводится как "первый столбец B".

В конце аналогичным образом мы передаём результат в Chi:

Chi = B(:, 1)

Последние два аргумента — это количество строк B, то есть 2, и выходной аргумент — статус. Его надо сравнить с нулём, чтобы понять, не произошла ли ошибка. В итоге получаем следующее:

function calc_Chi(Psi) result(Chi)
  complex, dimension(2), intent(in) :: Psi
  complex, dimension(2) :: Chi

  external :: cgesv
  complex, dimension(2, 2) :: A
  integer, dimension(2) :: ipiv
  complex, dimension(2, 1) :: B
  integer :: info

  A = U
  B(:, 1) = Psi
  call cgesv(2, 1, A, 2, ipiv, B, 2, info)
  if (info /= 0) then
    error stop "cgesv error"
  end if
  Chi = B(:, 1)
end function

Строка external :: cgesv сообщает компилятору, что где-то во внешней библиотеке есть процедура (или функция) с таким названием. Остальное думаю вполне понятно.

Чтобы скомпилировать получившийся код, надо подправить проектный файл, fpm.toml. Во-первых, разрешим использование external процедур, для чего з��меним implicit-external = false на implicit-external = true. Во-вторых, слинкуем программу с LAPACK и BLAS, добавив link=["lapack", "blas"] в секцию [build]. После этого можно убедиться с помощью fpm build что всё собирается.

Выполним вычисление M значений Psi, увеличивая каждый раз время на DeltaT и собирая результат в два массива X и Y, где X(I) будет хранить время I-го шага, а Y(I) — квадрат модуля амплитуды вероятности пребывания в первом базисном состоянии. В дальнейшем мы используем собранные данные для построения графика.

subroutine run
  complex, dimension(2) :: Psi
  integer :: I
  real, dimension(0:M) :: X, Y
  real :: T

  T = 0.0
  Psi = Psi_0
  X(0) = T
  Y(0) = abs(Psi(1)) ** 2
  do I = 1, M
    T = T + DeltaT
    call evol(Psi)
    X(I) = T
    Y(I) = abs(Psi(1)) ** 2
  end do
end subroutine

subroutine evol(Psi)
  complex, dimension(2), intent(inout) :: Psi

  Psi = 2.0 * calc_Chi(Psi) - Psi
end subroutine

В этом коде два новшества. Первое — это массивы, нумерация элементов в которых начинается не с 1, а с нуля: dimension(0:M). Тут M — это максимальный номер элемента, так что всего в этих массивахэлемент. Второе — цикл по I от 1 до M, что выглядит как do I = 1, M.

Осталось вывести красиво результат, для чего нам пригодится библиотека DISLIN. Чтобы воспользоваться ей, включим в проект соответствующий модуль:

$ ln -s /usr/local/dislin/gf/dislin.f90 src

и добавим "dislin" в список подключаемых библиотек в fpm.toml.

Читаем код примера построения графиков, выкидываем лишние украшательства, читаем документацию процедур GRAF и CURVE и пишем такой код в конце нашей процедуры run:

call metafl('CONS')
call scrmod('REVERS')
call disini()
call pagera()
call complx()
call axspos(450, 1800)
call axslen(2200, 1200)
call color('FORE')
call graf(0.0, X(M), 0.0, 1.0, 0.0, 1.0, 0.0, 0.2)
call curve(X, Y, M + 1)
call disfin()

В начало модуля надо дописать use dislin.

Осталось только выкинуть процедуру say_hello и заменить её на run в src/tss.f90 и в app/main.f90.

Запускаем с помощью fpm run и получаем красивейший результат:

Отчётливая косинусоида, при том, что в программе нет ни одного вызова тригонометрических или показательных функций! Вероятность, как ей и положено, плавно плещется туда-сюда между базисными состояниями.

Полный код src/tss.f90:

module tss
  use dislin
  implicit none (type, external)
  private

  public :: run

  real, parameter :: A = 1.0
  complex, dimension(2, 2), parameter :: H = reshape([0.0, -A, -A, 0.0], [2, 2])
  complex, dimension(2), parameter :: Psi_0 = [1.0, 0.0]
  integer, parameter :: M = 100
  real, parameter :: DeltaT = 0.1

  complex, dimension(2, 2), parameter :: I = reshape([1.0, 0.0, 0.0, 1.0], [2, 2])
  complex, dimension(2, 2), parameter :: U = I + (0.0, 0.5) * DeltaT * H
contains
  subroutine run
    complex, dimension(2) :: Psi
    integer :: I
    real, dimension(0:M) :: X, Y
    real :: T

    T = 0.0
    Psi = Psi_0
    X(0) = T
    Y(0) = abs(Psi(1)) ** 2
    do I = 1, M
      T = T + DeltaT
      call evol(Psi)
      X(I) = T
      Y(I) = abs(Psi(1)) ** 2
    end do

    call metafl('CONS')
    call scrmod('REVERS')
    call disini()
    call pagera()
    call complx()
    call axspos(450, 1800)
    call axslen(2200, 1200)
    call color('FORE')
    call graf(0.0, X(M), 0.0, 1.0, 0.0, 1.0, 0.0, 0.2)
    call curve(X, Y, M + 1)
    call disfin()
  end subroutine

  subroutine evol(Psi)
    complex, dimension(2), intent(inout) :: Psi

    Psi = 2.0 * calc_Chi(Psi) - Psi
  end subroutine

  function calc_Chi(Psi) result(Chi)
    complex, dimension(2), intent(in) :: Psi
    complex, dimension(2) :: Chi

    external :: cgesv
    complex, dimension(2, 2) :: A
    integer, dimension(2) :: ipiv
    complex, dimension(2, 1) :: B
    integer :: info

    A = U
    B(:, 1) = Psi
    call cgesv(2, 1, A, 2, ipiv, B, 2, info)
    if (info /= 0) then
      error stop "cgesv error"
    end if
    Chi = B(:, 1)
  end function
end module tss

Литература:

Фейнман, Лейтон, Сэндс. Фейнмановские лекции по физике. Т. 3 (вып. 8, 9).
Scherer. Computational Physics. Simulation of Classical and Quantum Systems.
Koonin, Meredith. Computational Physics. Fortran Version.

Квантовая физика на компьютере. Часть 1. Система с двумя состояниями