Расчет возбужденных состояний квантовых систем методом Монте-Карло Текст научной статьи по специальности «Физика»

УДК 530.145

И. А. Семенихин, Ю. И. Ожигов

РАСЧЕТ ВОЗБУЖДЕННЫХ СОСТОЯНИЙ КВАНТОВЫХ СИСТЕМ МЕТОДОМ МОНТЕ-КАРЛО1

(кафедра квантовой информатики факультета ВМиК)

Введение. В последние два десятилетия произошел значительный прогресс в области моделирования сложных физических систем с помощью компьютера [1]. Компьютерное моделирование, или "компьютерный эксперимент", играет сегодня в физике не меньшую роль, чем обычные эксперименты. Тем не менее на сегодняшний день имеется большой класс задач, которые не могут быть эффективно рассчитаны на компьютере. Среди них такая чрезвычайно важная задача, как расчет многочастичных квантовых систем, имеющая массу приложений в области физики наноструктур, квантовой химии, биологии и др. В случае прямого счета данная задача требует экспоненциально больших вычислительных ресурсов, и это сильно сдерживает прогресс в вышеперечисленных областях. Одним из возможных путей ее решения является создание квантового компьютера [2, 3] — устройства, работающего по квантово-механическим законам и тем самым способного эффективно моделировать квантовые системы. В последние годы в этой области был достигнут значительный прогресс, в особенности в сфере реализации квантового компьютера на ионных ловушках [4]. Тем не менее вопрос о возможности получить масштабируемый квантовый компьютер в ближайшие несколько десятков лет все еще остается открытым. Поэтому важнейшей задачей является создание эффективных математических моделей и алгоритмов, способных как можно более полно и точно описывать сложные квантово-механические системы, не требуя при этом экспоненциально большого вычислительного ресурса.

В области моделирования электронной структуры "из первых принципов", когда учитываются квантовые корреляции, наиболее популярными сегодня являются расчеты с использованием функционала плотности [5]. Данный метод, являющийся достаточно точным и быстрым, тем не менее имеет ряд существенных ограничений [6]. В частности, отсутствие полного математического описания так называемого обменно-корреляционного функционала энергии (exchange-correlation energy functional), который в общем случае заранее не известен, может привести к большим погрешностям в вычислениях, а иногда и к качественно неправильному результату. Альтернативным и в то же время взаимодополняющим методом в том случае, когда важна точность результатов, является квантовый метод Монте-Карло (КМК) [7]. Данный метод явно учитывает квантовые корреляции и при этом эффективно масштабируется как третья степень от числа частиц в системе. Для небольших систем метод КМК дает точность, достигающую химической точности (около 1 ккал на моль или 0,04 eV на молекулу), что недостижимо при расчетах в рамках теории среднего поля, в частности в методе Хартри-Фока. Такая точность вполне достаточна для расчета межатомных сил и химических реакций.

Существенной чертой метода КМК является возможность легкого распараллеливания вычислений, что ввиду имеющихся на сегодня тенденций развития вычислительной техники может стать решающим преимуществом метода КМК над другими методами расчетов. Здесь следует отметить, что уже сейчас квантовые методы Монте-Карло позволяют рассчитывать системы, содержащие до тысячи электронов [7].

Квантовый метод Монте-Карло является общим названием для целого ряда методов расчета, использующих случайные величины. Простейшим из КМК является вариационный метод Монте-Карло (variational Monte Carlo) (ВМК). В данном методе расчета, основанном на вариационном принципе квантовой механики, интегрирование пробной волновой функции проводят с помощью метода Монте-Карло. Для системы, содержащей N электронов, где N 1, необходимо рассчитывать ЗА-мерные интегралы, и метод Монте-Карло оказывается намного более эффективным, чем детерминированные подходы. Главным недостатком ВМК является сильная зависимость точности получаемого результата от выбора пробной волновой функции. Данный недостаток отсутствует в другом КМК-методе, так называемом диффузионном методе Монте-Карло (diffusion Monte Carlo) (ДМК), который будет подробно рассмотрен далее. Существует еще ряд методов КМК, например метод вспомогательного поля

1 Работа поддержана грантом INTAS 04-77-7289 и фондом NIX Foundation — грант F793/8-05.

(auxiliary-field Monte Carlo, см., обзор [8]), применяющийся для расчетов модельных гамильтонианов, в частности модели Хаббарда. Необходимо также упомянуть метод КМК, основанный на интегралах по траекториям (path integral Monte Carlo, см., например, обзор [9]), позволяющий проводить расчеты при ненулевой температуре и при этом не требующий нахождения конкретных возбужденных состояний и точных волновых функций. Данный метод применяют, как правило, к системам, претерпевающим фазовые переходы, в частности он используется для моделирования жидкого гелия.

Следует отметить, что наиболее распространенные и имеющие наибольшее значение КМК-методы такие, как ВМК и ДМК, были разработаны для исследования основного состояния квантовых систем, и имеется ряд трудностей, мешающих их применению для расчета возбужденных состояний квантовых систем [7]. В нашей работе мы остановимся на исследовании диффузионного метода Монте-Карло и рассмотрим простой подход, позволяющий в рамках данного метода проводить расчет первых возбужденных состояний.

Описание алгоритма. Метод ДМК основан на решении нестационарного уравнения Шрёдинге-ра в мнимом времени (1). При этом уравнение Шрёдингера переходит в уравнение диффузии, решение которого при больших значениях времени сходится к основному состоянию соответствующего стационарного уравнения Шрёдингера:

- t) = t) + (V(R) - Е,)ф(К, t) = Нф(П, t). (1)

Здесь ^(R, i) — волновая функция системы из iV-частиц, R и V2 — соответственно координата и оператор Лапласа в ЗА^мерном пространстве, V(R) — оператор потенциальной энергии. Величина Es служит для задания начала отсчета энергии. Уравнение (1) написано в атомных единицах, которые мы будем использовать везде далее.

Как уже говорилось, уравнение (1) характеризуется тем, что произвольное начальное состояние ф(R, i = 0) сходится при t —> оо к основному состоянию соответствующей квантово-механическоей системы.

Для практического применения алгоритма используется выборка по значимости (importance sampling), которая дает возможность работать с сильно меняющимися в пространстве потенциалами [10, 11]. Для этого вводится функция /(R, t) = t), где ^y(R) — пробная волновая функция

(trial or guiding wave function). Подставляя /(R, t) в (1), получаем

- R, t) = t/>T(R)#V>t(R)~7(R, t) = - V/CR, t) + VMR)/(R, t)] + (El(R) - Es)f(R, t) =

= [T + D + ÉL]f(R,t), (2)

где ü£)(R) — ЗА^мерная дрейфовая скорость, определяемая как

vd(R) = фт(К)~1^фт(R) = Vin \фт(р.)\ , (3)

El — локальная энергия:

ЕЬ(R) = фт(R)"1 (-^V2 + V(R)) фт(В.). (4)

Уравнение (2) (уравнение Фоккера-Планка) совпадает с уравнением диффузии ансамбля классических частиц в SiV-мерном пространстве. Функция /(R, t) описывает распределение этих ЗА^мерных псевдочастиц. Поведение ансамбля определяется функцией Грина:

/(R, t + Ai) = J G(R,R', Ai)/(R',i) <ZR. (5)

В пределе бесконечно малого времени эволюцию системы можно разбить на три процесса: свободная диффузия с коэффициентом 1/2, соответствующая кинетической энергии Т в (2), смещение под действием силы D и процесс рождения/уничтожения псевдочастиц, определяемый величиной ЕЭти стохастические процессы моделируются методом Монте-Карло.

Вычисления проходят следующим образом: вначале распределяют псевдочастицы в ЗА^мерном пространстве с плотностью |^т|2) затем ансамбль псевдочастиц эволюционирует в соответствии с

уравнением (5). После периода установления, во время которого вклад возбужденных состояний в состояние системы уменьшается по экспоненте, псевдочастицы распределяются по пространству с плотностью вероятности, пропорциональной /(Л, = £), где — функция основного состояния. Далее мы можем начать накапливать статистику по нужным нам величинам, например по энергии основного состояния, равной

где Ri — координата г-ж псевдочастицы.

Выборка по значимости (2) решает одновременно несколько проблем. Во-первых, она избавляет нас от сильных флуктуаций локальной энергии £x(R), которые связаны с расходимостью кулоновского потенциала. При правильном выборе пробной волновой функции t/>y(R) значения локальной энергии будут близки к энергии основного состояния и практически постоянны в разных точках ЗА-мерного пространства, минимизируя дисперсию Е.

Во-вторых, использование выборки по значимости позволяет работать с волновыми функциями, меняющими знак в разных точках пространства и соответственно имеющими нули. Этим свойством обладают все многоэлектронные волновые функции, и применение выборки по значимости помогает в некотором смысле обойти так называемую проблему знака. Поскольку поверхность нулей волновой функции ^(R, i) полностью определяется поверхностью нулей пробной волновой функции t/>y(R), данный подход называют приближением фиксированной поверхности нулей (fixed-node approximation) [12-15]. В процессе временной эволюции распределение /(R, t) сходится к /(R, t —> оо) = t/>y(R)t/>(R), где ф(И) соответствует состоянию с минимальной энергией, имеющему такое же расположение нулей, как ф(И). Полученная энергия будет точной, если поверхность нулей ф(И) совпадает с нулями собственной функции Н. В общем случае можно показать [7], что ошибка в вычислении средних значений операторов, коммутирующих с гамильтонианом, составит второй порядок малости от ошибки в определении поверхности нулей пробной волновой функции.

Метод ДМК в приближении фиксированной поверхности нулей позволяет также рассчитывать возбужденные состояния. Если расположение нулей пробной волновой функции ф(И) такое же, как у возбужденного состояния Н, то алгоритм ДМК даст точное собственное значение энергии этого состояния. Главная проблема заключается в нахождении поверхности нулей возбужденных состояний. Решение данной проблемы на сегодняшний день не найдено. В нашей работе мы предлагаем алгоритм, позволяющий приближенно находить энергию возбужденных состояний без знания расположения нулей соответствующей собственной функции.

Предположим, что имеется система из А взаимодействующих частиц в 3-мерном пространстве. Это могут быть как бозоны, так и фермионы. Требуется найти первый возбужденный уровень этой системы. Формально данная задача может быть представлена как задача об одной ЗА-мерной частице. Если мы добавим к нашей системе еще одну такую же ЗА-мерную частицу, то при условии антисимметрии волновой функции этих двух частиц и отсутствии вырождения, вторая частица займет первый возбужденный уровень системы. Таким образом, энергия двух таких невзаимодействующих ЗА-мерных "ферми-частиц" будет равна сумме энергий основного и первого возбужденного состояний нашей первоначальной системы. Зная энергию основного состояния, мы легко получим энергию, соответствующую первому возбужденному состоянию. Аналогично можно найти и более высокие уровни энергии. Платой за простоту данного алгоритма является увеличение размерности задачи. Для нахождения к-го уровня энергии мы должны будем увеличить размерность задачи как минимум в к раз. Тем не менее для нахождения нескольких первых возбужденных уровней энергии данный метод может быть полезен, поскольку он не требует от нас знания поверхности нулей волновой функции возбужденного состояния.

Исследуем работу вышеизложенного алгоритма на примере расчета первых возбужденных состояний нескольких простых систем. В качестве первого примера возьмем гармонический осциллятор. Соответствующий гамильтониан запишется как

Для расчета первого возбужденного состояния (6) нужно будет найти основной уровень системы с гамильтонианом

, { d 1 d ж 1 ж о , , , ,

+ Т + (7)

при условии антисимметрии волновой функции ф(х±, х2). Пробную волновую функцию выберем в виде

Фт(х 1,ж2) = (хг (8)

Первый множитель в уравнении (8) определяет антисимметрию волновой функции. Для простоты мы взяли просто разность координат. Второй множитель мы выбрали в виде экспоненты для более быстрой сходимости метода ДМК, поскольку волновая функция основного состояния гармонического осциллятора пропорциональна ехр( — х2/2). Как легко заметить, мы получили точный вид искомой волновой функции уравнения (7), которая является антисимметризованным произведением основного и первого возбужденного состояний гармонического осциллятора. В данном случае метод ДМК сразу даст точное значение энергии, равное сумме энергий основного и первого возбужденного состояний: Е= 1/2 + (1/2 + 1) = 2.

Совершенно аналогичную картину мы получим для второго возбужденного состояния, если возьмем антисимметричную пробную волновую функцию в виде

фт(х 1,ж2,ж3) = (Ж1 - х2){х1 - х3)(х2 - Х3)е-{х1+х2+х^/2. (9)

Данная функция, как легко показать, является антисимметризованным произведением первых трех истинных волновых функций гармонического осциллятора. Метод ДМК даст точное значение энергии и в этом случае: Е = 1/2 + (1/2 + 1) + (1/2 + 2) = 4,5. Этот процесс можно продолжить для произвольного возбужденного состояния гармонического осциллятора.

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

В качестве следующего примера проведем расчет первого возбужденного состояния атома водорода:

Н(г)ф(г)= (10)

Для расчета возбужденного состояния (10) необходимо найти основное состояние системы, имеющей следующий гамильтониан:

ВД(гь г2) = (Я(п) + Н(т2))ф(тъ г2) = (~\V\ - - 1 - -) ф(тъ г2), (11)

V 2 2 r-i r2J

где волновая функция ф(ri,r2) должна быть антисимметричной. Аналогично предыдущему случаю выберем пробную волновую функцию (trial wave function) в виде произведения функций основного состояния (10) и простейшей антисимметричной части:

фт(г1,г2) = (г1-г2)е-^+г^. (12)

Экспоненциальный множитель мы опять же выбрали из соображений лучшей сходимости метода ДМК. Множитель, задающий расположение нулей, мы взяли в виде (г 1 — г2), а не в виде |г 1 — г2|, поскольку истинная поверхность нулей должна быть (ЗА — 1)-мерной поверхностью, в нашем случае 5-мерной. Пробная волновая функция (12) уже не является собственной функцией гамильтониана (11), так же как и поверхность нулей (12) уже не отвечает расположению нулей истинной волновой функции основного состояния гамильтониана (11). Тем не менее, если поверхность нулей (12) достаточно близка к истинной, то с помощью данной функции можно получить достаточно точные результаты, поскольку ошибка при вычислении значения энергии будет иметь второй порядок ошибки в определении расположения нулей фт. Как мы уже говорили, значение первого возбужденного уровня энергии

(10) есть разность энергий основного состояния (11) и основного состояния (10). На рис. 1 показаны результаты расчета энергии первого возбужденного состояния атома водорода методом ДМК для разной величины шага АЬ. В данном случае мы используем несколько модифицированный нами метод ДМК второго порядка [16], поэтому при уменьшении шага значение энергии Е(АЬ) должно квадратично сходится к точному значению Е\ = —1/4.

Рис. 1. Полученная методом ДМК величина энергии первого возбужденного состояния атома Н как функция шага по времени At. Статистическая погрешность значений энергии меньше величины символа. Горизонтальная линия отвечает истинному значению энергии Ei = —1/4

В качестве еще одного примера проведем расчет первого возбужденного уровня атома Не:

НФ(гг, г2) = -- 1 - 1 + ф(гъ r2). (13)

Для этого введем гамильтониан Н\.

Нг = Я(П, г2) + Я(гз, г4) = - Jv? - iv2 - iv2 - iv2 - 1 - 1 - 1 - 1 + J-. (14)

2 2 2 2 ri r2 r-i rA r34

Отметим, что электроны 1 и 2 никак не взаимодействуют с электронами 3 и 4, хотя и взаимодействуют между собой. Это является выражением того, что электроны 1 и 2 составляют первую 6-мерную ферми-частицу, а электроны 3 и 4 — вторую. И в соответствии с нашим методом расчета эти две частицы не взаимодействуют. Выберем пробную волновую функцию в виде

ЫГ1,Г2,ГЗ,Г4) = (Д12 - Д34) -е 2^Г'Х(Г1,Г2)Х(ГЗ,Г4), Х(Г,Г') = ео,5|р-р'|/(1+о,2|р-р'|)- (15)

Величина R42 является радиус-вектором первой 6-мерной ферми-частицы, и R23 — второй. Антисимметрию по перестановкам этих двух частиц мы ввели самым примитивным способом с помощью разности величин этих векторов. Функцию х(г5 г') мы выбрали из соображений лучшей сходимости метода ДМК, она позволяет функции (15) удовлетворить так называемым касповым условиям (cusp conditions). Результаты расчета приведены на рис. 2. Видно, что хотя фт была выбрана не самым оптимальным образом, алгоритм сходится к значению, близкому истинному. Этот результат можно значительно улучшить, если выбрать пробную волновую функцию в виде

Фт(Г1, г2, г3, г4) =

Ф\ (ri)<^o (гг) <Мгз)<Мг4)

Фо {г1)Фо (г2) <^о(гз)<^о(г4)

x(ri,r2)x(r3,r4), (16)

где детерминант составлен из волновых функций фо(г) основного и первого возбужденного ф\(г) состояний водородоподобного атома с зарядом 2 = 2. Результаты расчета приведены на рис. 3. Видно, что теперь алгоритм сходится к истинному значению энергии первого возбужденного уровня.

Таким образом, описанный нами алгоритм позволяет находить возмущенное состояние квантово-механической системы, сводя эту задачу к расчету основного состояния системы вдвое большей размерности. Следует отметить, что данный метод, конечно, не решает проблему нахождения точной поверхности нулей возбужденного состояния, однако его можно использовать в качестве первого приближения, когда мы не имеем информации о расположении нулей истинной волновой функции.

Рис. 2. Полученная методом ДМК величина энергии первого возбужденного состояния атома Не как функция шага по времени At. Статистическая погрешность значений энергии меньше величины символа. Горизонтальная линия отвечает истинному значению энергии Ех = —2,14597404 ...

Рис. 3. Полученная методом ДМК величина энергии первого возбужденного состояния атома Не как функция шага по времени At при использовании пробной волновой функции (16). Статистическая погрешность значений энергии меньше величины символа. Горизонтальная линия отвечает истинному значению энергии

Ег = -2,14597404...

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

СПИСОК ЛИТЕРАТУРЫ

1. Ожигов Ю. И. Кванты амплитуд в моделировании многочастичных систем // Микроэлектроника. 2006. 35. № 1. С. 61-77.

2. Валиев К. А. Квантовые компьютеры: смена парадигмы вычислений // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2005. № 2. С. 3-16.

3. Валиев К. А., Кокин А. А. Квантовые компьютеры: надежды и реальность. М.; Ижевск: НИЦ "Регулярная и хаотическая динамика", 2002.

4. SasuraM., BuzekV. Cold trapped ions as quantum information processors // lanl e-print quant-ph/0112041. 2001.

5. Kohn W., Sham Г.Л. Self-consistent equations including exchange and correlation effects // Phys. Rev. A. 1965. 140. P. A1133-A1138.

6. Parr R. G., Yang W. Density-functional theory of atoms and molecules. N.Y.: Oxford University, 1989.

7. Foulkes W.M.C., Mitas Г., Needs R. J., Rajagopal G. Quantum Monte Carlo simulations of solids // Rev. Mod. Phys. 2001. 73. P. 33-83.

8. Senatore G., March N. H. Recent progress in the field of electron correlation // Rev. Mod. Phys. 1994. 66. P. 445-479.

9. Ceperley D.M. Path integrals in the theory of condensed helium // Rev. Mod. Phys. 1995. 67. P. 279-355.

10. Grimm R. C., Storer R. G. Monte-Carlo solution of Schrodinger's equation //J. Comput. Phys. 1971. 7. P. 134-156.

11. Kalos M.H., Levesque D., Verlet L. Helium at zero temperature with hard-sphere and other forces // Phys. Rev. A. 1974. 9. P. 2178-2195.

12. Anderson J.B. A random-walk simulation of the Schrodinger equation: H^ //J. Chem. Phys. 1975. 63. P. 1499-1503.

13. Anderson J.B. Quantum chemistry by random walk //J. Chem. Phys. 1976. 65. P. 4121-4127.

14. Moskowitz J.W., Schmidt K.E., Lee M.A., Kalos M.H. A new look at correlation energy in atomic and molecular systems. II. The application of the Green's function Monte Carlo method to LiH // J. Chem. Phys. 1982. 77. P. 349-355.

15. Reynolds P.J., Ceperley D.M., Alder B.J., Lester W.A. Fixed-node quantum Monte Carlo for molecules // J. Chem. Phys. 1982. 77. P. 5593-5603.

16. Siu A. Chin. Quadratic diffusion Monte Carlo algorithms for solving atomic many-body problems // Phys. Rev. A. 1990. 42. P. 6991-7005.

Поступила в редакцию 10.02.06