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

Сер. 4. 2009. Вып. 2

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

ФИЗИКА

УДК 51-72, 519.245, 539.182

М. А. Вознесенский, П. Н. Воронцов-Вельяминов, А. П. Любарцев

РАСЧЁТЫ РАВНОВЕСНЫХ СВОЙСТВ КВАНТОВЫХ СИСТЕМ С КУЛОНОВСКИМ ВЗАИМОДЕЙСТВИЕМ

МЕТОДОМ МОНТЕ-КАРЛО В РАСШИРЕННОМ АНСАМБЛЕ *)

Введение. Метод Монте-Карло интегралов по траекториям, основанный на фейн-мановской трактовке квантовой механики [1], является одним из наиболее надёжных способов изучения квантовых систем при конечных температурах. Однако при его реализации возникают определённые трудности, связанные с поведением квантовых систем при низких температурах. Одной из таких трудностей является проблема знаков для систем фермионов. Волновая функция, описывающая систему нескольких фер-мионов, меняет знак (антисимметрична) при перестановке любых двух частиц. Соответственно, матрица плотности, с помощью которой описываются свойства квантовых систем, также должна быть антисимметричной относительно перестановок частиц. Это приводит к появлению отрицательных вкладов в матрицу плотности, и как следствие, в статистическую сумму и в выражения для средних. При понижении температуры разность положительных и отрицательных вкладов экспоненциально уменьшается и становится трудно отделимой от статистического шума.

Существует достаточно большое число работ, посвященных попыткам решить указанную проблему тем или иным способом. Основные идеи попыток решения проблемы знака в рамках методов Монте-Карло заключаются в сокращении основной доли знаков не стохастически, а на этапе вычисления весов траекторий. Например, в работе [2] было предложено симметризовать (антисимметризовать) каждую высокотемпературную матрицу плотности, входящую в разложение. Это приводит к набору незамкнутых полимеров с гармоническими связями, на концы которого наложен особого вида потенциал - детерминант из таких связей. Этот подход, в силу геометрических особенностей конфигурационного пространства, даже позволяет решить проблему знаков для одномерных систем со взаимодействием, однако для систем в пространствах с большей размерностью при низких температурах этот метод оказывается мало эффективен. Существуют и более изощрённые способы, так например, идея Ньюмана и Куки [3] основана на предложении заменить каждую вершину траектории целым набором состояний, тогда ещё на этапе вычисления веса такой траектории будет происходить существенная доля сокращения знаков. Данный подход является численно точным, однако расширить его на системы, содержащие больше двух частиц, пока не удалось.

В настоящей работе предлагается метод ослабления проблемы знаков, основанный на точном вычислении отношений статистических сумм циклов [4]. Метод апробирован

*) Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант № 05-02-17428, а так же Шведской Королевской академии наук

© М. А. Вознесенский, П. Н. Воронцов-Вельяминов, А. П. Любарцев, 2009

на системе двух невзаимодействующих фермионов в гармоническом поле и применён к системам двух фермионов с кулоновским отталкиванием в гармоническом и кулонов-ском внешнем поле.

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

ставлении

p (г г'; в)

exp ( -в^

А

может быть использована для вычисления средних [5]:

(А) = (г (Р))-1 У * I (г ехр (-ей

где г = (г(1),)) - вектор из координат всех частиц, т(к) - координата к-й частицы (^-мерный вектор в случае ^-мерного пространства), г (в) - каноническая статистическая сумма

Z (в)

dr (r

exp (-вН

в = (kT) 1 - обратная температура в энергетических единицах, k - постоянная Больцмана. n

С помощью очевидного тождества exp (-ря) = ^exp (-РЯ/n, где n - любое натуральное число, матрица плотности р и статистическая сумма Z могут быть выражены в следующем виде:

/п П

dr2... dr^^[ ^Гі exp (-в#/nj Гі+1^ ■' i=1

Z (в)

dr1... drn П Гі |exp (-в#/

где в уравнении (1) г1 = г, rn+1 ее г', а в уравнении (2) rn+1 = г1. Величина p (гі, гі+1; в/n) =

(1) (2)

гі+1) является матрицей плотно-

Гі+1) ,

exp \—РН/nj

сти, соответствующей более высокой температуре nT, или обратной температуре P/n.

Имея ввиду, что оператор Н есть сумма операторов кинетической и потенциальной энергий, Н = K + V, высокотемпературную матрицу плотности можно записать в виде

exp (-Р (К + /n^ = exp ^—piK/n^ exp (-pV/n^ + O ((p/n)3) .

Отсюда приходим к следующей аппроксимации для высокотемпературной матрицы плотности:

r exp (-Р (k + У) /nj ri+i^ « ^ exp (-рк/nj exp (-pV/nj ri+^ =

exp (-pK/nnj rj+^ exp ( РV (ri+1)) =

4 -dN/2

f— \27lP^2 )

I mn < <2 в т/ / \

ЄХР(“2р^(Гї“Гї+і) '/(r*+l)

(З)

r

r

r

Здесь в последнем равенстве использовано выражение для матрицы плотности свободной частицы [1]. Подставляя (3) в соотношение (2), получаем выражение для статистической суммы в п-вершинном приближении:

/ \ —ndN/2 гг / 'п 0 п

^ (Р) = (^) / *1- / еХР ^ (Г< “ Гт)2 “ п £{Гг) I ’ (4)

где по-прежнему

Гп+1 = Г1. (5)

В рамках п-вершинного приближения возникает квантово-классический изоморфизм: каждая квантовая частица представляется в виде формальной классической системы - «траектории» - набора из «вершин», соединённых последовательно между собой гармоническими пружинками с жёсткостью к = шп/$Ь?. Как видно из формулы (4), жёсткость пружинок прямо пропорциональна числу вершин п и обратно пропорциональна обратной температуре р. Потенциал V (г) включает в себя как взаимодействие со внешним полем, так и межчастичное взаимодействие.

Рассмотрим теперь систему N тождественных частиц со спином (фермионов или бозонов). Выражая матрицу плотности как антисимметричную или симметричную, в зависимости от статистики, сумму по всем N! перестановкам Р координат частиц в матрице плотности р(в) системы N различимых частиц, мы получаем следующее выражение для статистической суммы [6]:

£(Л,я) (р) = ±_^2^1Р] ( с1х^..Ях^р^ , ...Р (х^ ; р), (6)

! {Р}

где х(к) = (т(к), а(й)) обозначает совокупность пространственных координат и спиновой переменной соответствующую к-й частице, соответственно / dx(k) = ^с(ь) /

=1 соответствует симметричному случаю (бозоны), £, = —1 - антисимметричному (фермионы), [Р] - чётность перестановки.

Матрица плотности в уравнении (6) в случае, когда гамильтониан не зависит от спина, может быть представлена как произведение координатной и спиновой частей:

р(0) = р{0) ^(1),...г^),р ^(1)^ , ...Р ^^ р8 (р), (7)

р8 (Р) = 5 (а(1),Р (а(1))) •... • 5 (o(N),Р (o(N^ . (8)

Подставляя уравнения (7) и (8) в уравнение (6) приходим к следующему выражению для статистической суммы:

(Р) = (Р)> (9)

где

- спиновая, а

{Р}

к (р ) = ^ р8 (р )

v 1 ^'с(1) ,...с(ЯК ^ '

гр) (Р) = / dr(1)...dr(N )р(с°) (г(1) ,...г^ ),р(т(1)) ,...р(т^ ^ ; р) (10)

соответственно координатная часть статистической суммы, отвечающие перестановке Р.

Матрица плотности р(в) в выражении (10) может быть представлена в виде разложения (1) по высокотемпературным матрицам плотности, для которых в свою очередь будет справедлива аппроксимация (3). В итоге мы получим выражение для статистической суммы, аналогичное выражению (4), но, как не трудно понять, условие (5) в этом случае будет выглядеть так:

Гп+1 = Р (Г1) = (р (г(1)) , Р (г(2) ) ,...,Р (г^ ))) ,

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

Перепишем статистическую сумму (9) как сумму по классам перестановок [6]. В этом случае ^[р] ^ ^[о],

K (Р) ^ K (О) = (2в + 1)ЕV с(о),

где С (О) - число циклов длинной V в данном классе С, так что ^у С (О) есть полное число циклов в классе О ^1 < Х^=1 CV (О) < . Так как значение интеграла (10)

по всем пространственным переменным зависит только от циклической структуры перестановки, то г^^')0}) (Р) = гО) (Р). В свете этого статистическая сумма будет выглядеть следующим образом:

^(А,5) (Р) = м м 2\о) (Р). (п)

о

N

где М (О) = N!/ П (CV- число элементов в классе [6].

V=1

Конкретизируем выражение для статистической суммы в случае системы двух тождественных фермионов в термостате:

44) (Р) = \ [к (I2) V) (Р) - к (2) г(2) (Р)] . (12)

Здесь два класса перестановок обозначены как 12 (две отдельные частицы) и 2 (цикл из двух частиц), М (12) = М (2) = 1, г(12) (Р) обозначает статистическую сумму системы из двух различимых частиц при обратной температуре Р, а г(2) (Р) - соответственно статистическую сумму для цикла из двух частиц. В случае спина, равного нулю (в = 0), спиновая часть статистической суммы К (О) в (11) будет равна 1 для всех классов перестановок О. Если же спин равен половине (в = 1/2), то для N = 2 К (12) = 4, а К (2) = 2.

Метод расширенных ансамблей. Данный метод, позволяющий в ходе одного расчёта получить отношение статистических сумм, был предложен в работе [7]. Пусть Н (д) есть конфигурационная часть классического гамильтониана, описывающего рассматриваемую систему. Введём приведённый гамильтониан Н = —РН (д). Пусть Но - гамильтониан системы, для которой статистическая сумма может быть рассчитана точно. Построим цепочку гамильтонианов

Но, Н1, ...Нм = Н,

в которой с изменением индекса т от 0 до М гамильтониан ко постепенно превращается в км, например, по линейному закону:

кт — (1 ^т) к0 + ^ткМ , (13)

где 0 — ^о < ^1 < ••• < ^м — 1- Для каждого т имеем канонический ансамбль с конфигурационным интегралом

2т — ^ ехр (кт) •

Составим расширенный и модифицированный ансамбль со статистической суммой:

м

2 — ^ %т ехр (-Пт), (14)

т=0

где Пт - так называемые балансировочные параметры. Расширение ансамбля означает введение суммирования по дополнительному индексу т, причём \т в частных случаях может совпадать с обычными термодинамическими параметрами: температурой, давлением, числом частиц и т. п. Модификация означает введение весового множителя ехр(-пт)- Каждый из канонических ансамблей с номером т становится подансамб-лем расширенного ансамбля (14). В построенном расширенном ансамбле организуется МК-блуждание в соответствии со стандартным алгоритмом Метрополиса. При этом реализуют два типа шагов: изменение конфигурации системы внутри подансамбля (\т не меняется) и переходы между подансамблями (изменение \т при неизменных значениях координат). Вероятность перехода 1 ^ 2 определяется стандартным выражением алгоритма Метрополиса

р1^2 — тіп (і,ехр ((к(2) — — (к(1) — , (15)

причём для шагов первого типа п(1) — П(2) •

В ходе описанной МК-процедуры подсчитывается число раз пт, когда система побывала в подансамбле с индексом т. В результате получаем оценку вероятности состояния характеризуемого значением Хт: рт ^ пт/пМс, где пМс - полное число МК-шагов.

С другой стороны, исходя из (14), имеем рт — 2т ехр(—пт) /2, и, следовательно,

пт 2т / , ч

--- — ехр(-Г|т +щ).

пк 2к

Таким образом, можно получить отношение статистических сумм для любых двух под-ансамблей из расширенного ансамбля (14). Поскольку для подансамбля с индексом «ноль» статистическая сумма известна, можно получить значения статистической суммы для любого подансамбля, в частности исследуемого, обозначенного индексом М:

пм

— А)----ехР (-Л' + Лм) •

по

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

предполагает проведение нескольких предварительных прогонов с целью оптимизации значений балансировочных параметров, однако эту процедуру оказывается возможным автоматизировать при помощи алгоритма Ванга-Ландау [8].

Алгоритм Ванга-Ландау. Данный алгоритм [8] является самонастраивающейся процедурой организации марковского процесса, реализующего равномерное блуждание по пространству энергий в генерируемой выборке.

Блуждание по алгоритму Метрополиса приводит к канонической функции распределения p (E) = Q (E) exp (—E/квT), где Q(E) - плотность числа состояний. Если больцмановский фактор exp (—E/квT) заменить весом 1/Q (E), естественный вес Q (E) будет полностью скомпенсирован, и распределение по энергиям p (E) будет равномерным. Однако плотность числа состояний a priori неизвестна.

Для реализации алгоритма рассматриваемый промежуток энергий Emin ^ E ^ Emax делится на конечное число интервалов Nb. Вводятся два массива размером Nb каждый. Один предназначается для весов Qm, второй - для чисел посещений nm. Все начальные значения Qm полагаются равными единице, значения nm - нулю. После этого начинается блуждание по области конфигурационного пространства, отвечающей выбранному диапазону энергий.

Пусть Ei и E2 есть энергии исходного и пробного состояний, соответственно. Каждое из этих значений попадает в какой-то интервал, и каждому из них соответствует какой-то индекс. Для реализации веса 1/Q(E) вероятность перехода следует выбрать как

Независимо от того, состоялся переход или нет, значение Пт, соответствующее текущему значению энергии, умножается на некоторое ] - модифицирующий фактор, т. е. 0.т заменяется на ]Пт, а счётчик посещений пт увеличивается на 1.

Чтобы избежать работы с большими числами удобно, ввести функцию Б (Е) —

— ІпП (Е). Тогда вместо массива для 0,т вводится массив сумматоров для Бт, и вероятность перехода переписывается как

Модификация значения Бт сведётся к добавлению ДБ — 1п ]. Начальное значение модифицирующего фактора может быть взято, например, f — е — 2, 7182^ [8].

С фиксированным значением ] выполняется серия шагов такого размера, чтобы гистограмма пт стала «плоской». После выполнения этой серии модифицирующий фактор уменьшается по некоторому закону, например: / —>■ а/7 (АБ —*■ АБ/2). Возможно использовать любую функцию, монотонно убывающую к единице. Счётчики посещений пт обнуляются, и запускается новая серия. После её выполнения, описанные выше действия повторяются. Процедура продолжается до тех пор, пока модифицирующий фактор ДБ не станет меньше некоторого заранее определённого значения. Очевидно, что модифицирующий фактор действует как параметр, определяющий точность полученных значений плотностей состояний. После этого выполняется несколько серий с фиксированными значениями величин Бт. Поскольку невозможно получить идеально равномерное блуждание, плоской гистограммой считают такую, в которой значения пт для любого т отличаются от среднего значения {пт) не более чем на несколько процентов.

P (Ei ^ E2) = min (X exp (Sm(Ei) — Sm(E2})) •

(І6)

5. Алгоритм Ванга-Ландау в методе расширенных ансамблей. В данной работе алгоритм Ванга-Ландау используется для настройки балансировочных параметров Пт в методе расширенных ансамблей. Следуя методу расширенных ансамблей, мы строим цепочку гамильтонианов, например, по закону (13). Вводим два массива размером М +1, один для весов Бт, второй для чисел посещения пт. В начале расчёта все Бт и пт берутся равными нулю.

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

Шаги, осуществляющие попытку перехода между подансамблями, реализуется следующим образом. С вероятностью 0,5 выбирается «направление» - в сторону увеличения или уменьшения индекса подансамбля на 1. В том случае, если мы находимся в крайнем подансамбле с индексом 0 или М, то достоверно совершается попытка перейти в единственный соседний подансамбль с индексом 1 или М — 1, соответственно. Исходя из (16), вероятность перехода определяется выражением

рт^т±1 = (1, ехр (Ьт±1 кт $т±1 + ^т)) • (17)

Вне зависимости от того, состоялся переход или нет, в счётчик пт, соответствующий текущему подансамблю, добавляется единица, а значение Бт увеличивается на величину

. Таким образом, последним действием мы уменьшаем вес текущего подансамбля, а, следовательно, вероятность попасть в него в следующий раз. Сравнивая (15) и (17), не трудно убедиться, что балансировочные параметры цт в точности совпадают с параметрами Бт алгоритма Ванга-Ландау.

Так как каждый из крайних подансамблей (т равно 0 или М) граничит только с одним подансамблем (т = 1 или М — 1, соответственно), а все остальные с двумя (с индексами т± 1), очевидно, что при попадании в крайний подансамбль к значениям элементов массивов пт и Бт необходимо добавлять удвоенные добавки.

Изменение конформации производится при неизменном значении Хт. Тип шага выбирается случайно, но между ними сохраняется определённая пропорция - на один шаг, осуществляющий смену подансамбля, приходится п шагов изменения конформации (п - число вершин).

Описанным образом выполняется серия, обычно порядка 106 + 107 шагов. Такого количества вполне достаточно, по крайней мере, для рассматриваемых в этой работе систем, для формирования равномерного распределения чисел посещения подансам-блей пт.

С использованием текущих значений Бт и пт вычисляется статистическая сумма исследуемой системы с индексом М для оценки погрешности и достаточности числа шагов. Далее, величина Д$ домножается на константу а < 1, и запускается новая серия, тем самым начинается более тонкая настройка балансировочных параметров.

Применение к простейшим задачам квантовой статистики. Описанная совокупность методов была всесторонне протестирована нами на системах из 2-х и 3-х ферми-частиц в гармоническом внешнем поле [9]. Мы рассмотрели систему, состоящую из двух невзаимодействующих частиц, для которой существует аналитическое решение. Расширяя ансамбль по перестановке, иначе говоря, переходя от системы двух различимых частиц к циклу из двух частиц, нам удалось воспроизвести отношения статистических сумм, соответствующих двум имеющимся классам перестановок - полученные результаты хорошо совпадают с точными значениями. В серии расчётов для различных температур удалось получить хорошие результаты вплоть до обратной

температуры Ь = 10. При таких температурах значения энергии выходят на уровень основного состояния.

Были также проведены расчёты для систем двух и трёх квантовых частиц с взаимодействием в гармоническом внешнем поле. В этом случае каждый класс перестановки рассматривался отдельно. В качестве опорной системы была взята система без взаимодействия. Расширение ансамбля сводилось к включению кулоновского отталкивания между частицами. Энергия основного состояния системы двух тождественных бесспи-новых частиц с кулоновским отталкиванием в одномерном гармоническом поле получилась равной Е = 2, 87 [9] (энергия системы двух невзаимодействующих тождественных бесспиновых частиц Е = 2). Для двух бесспиновых фермионов в трёхмерном гармоническом поле было получено Е = 4, 73 (энергия системы без взаимодействия Е = 4).

Полученное нами значение энергии для системы тождественных бесспиновых частиц с кулоновским отталкиванием в одномерном кулоновском поле (Е = 2, 87 [9]) оказалось несколько выше значения Е = 2, 67, полученного в [10] другим методом. Данные в [10] вызывают больше доверия, так как фейнмановские континуальные интегралы были взяты там численно в 2048-вершинном приближении. В предыдущей работе [9] мы использовали малое число вершин. В настоящей работе нам удалось улучшить результат моделирования.

Аппроксимация для кулоновского потенциала. В моделировании кулонов-ских систем имеется принципиальная сложность - точка сингулярности притягивающего кулоновского потенциала. Во-первых, в машинном счёте существует максимальное представимое значение в пределах установленной точности, следовательно, значение потенциала в точке сингулярности должно быть заменено конечным значением. Во-вторых, в методе МК-интегралов по траекториям экспонента, стоящая под знаком интеграла в статистической сумме (4), является весом в процедуре существенной выборки. Соответственно возникает эффект «затягивания», траектория будет стремиться локализоваться в малой окрестности точки притяжения, и результаты моделирования будут сильно искажены. Для нейтрализации описанной проблемы сингулярности необходимо сгладить потенциал в особой точке и подобрать число вершин, соответствующее сглаживанию. Так как область в окрестности нуля заметно влияет на результаты моделирования, сглаживание необходимо осуществлять так, чтобы минимизировать расхождение с истинным потенциалом. При выборе же числа вершин необходимо исходить из следующих соображений. С одной стороны, при недостаточном числе вершин мы будем наблюдать эффект затягивания траектории, с другой стороны - чрезмерное увеличение числа вершин приводит к увеличению необходимого объема вычислений. Поиск оптимального режима является очень кропотливой задачей.

В данной работе для устранения проблемы сингулярности мы воспользовались предложенной в работе Кола-де-Рэдта [11] идеей заменить кулоновский потенциал эффективным потенциалом, полученным в результате применения неравенства Дженсона и последующей свёртки истинного потенциала с пропагатером свободной частицы. Слагаемое в показателе экспоненты в уравнении (4), описывающее потенциальную энергию кулоновского взаимодействия двух заряженных частиц, заменится на интегральное выражение:

-Г2д) (ГУ.Г23.Г1Й+1).Г2Й+1)Д).

3 = 1 3 = 1

где гг(п+1) = rг1, а

I (т^,Г2з,Т1и+1),Г20'+!),^ =

п ,

Т ((|)1/2 |(гу -Г2^Ф + (г1(і+1) -Г20-+1))с^ф|)

с£ф- 4 7

0

|(гу - Т2з )tgф + (Гі0-+і) - г2(.,'+1)) СІ&ф|

(приведено в безразмерных переменных г = (тй/Й)1/2 г, Ь = кав, ® = те4/к3, т - масса частиц, е - заряд электрона).

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

IV ( сЦ, сі;+і, 8, ^

«л1/2 ( (п\1/2 (п\1/2

ь) 1 (ь) ‘‘-(ь) ^«’0

Ь) ■\Ь.

Ч2 егї Г(d2tg2ф + 2(1■ й]+1 соя 0 + d2j+1 ctg2ф)1/2 с^ф

0

ф + 2^ йз+1 соя 0 + й2+1С^2ф)

1/2

где и ]-+1 - длины векторов Гу — Г2у и гу+1 — Г2у+1, а 0 - угол между ними.

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

Применение к кулоновским системам. Рассмотрим систему двух тождественных взаимодействующих частиц, подчиняющихся статистике Ферми во внешнем поле. Перепишем статистическую сумму (12) в следующем виде (полагая при этом спин частиц равным 0):

гЪл) = -

Z 0

12

2

70

Z12

7огі

“ 2 Ж,

где ноликами обозначены соответствующие статсуммы классов перестановки частиц в опорной системе.

В качестве опорной системы мы использовали систему невзаимодействующих частиц в гармоническом поле, для которой известны точные формулы для статистической суммы и энергии. Статистическая сумма N различимых частиц в ]-мерном пространстве в этом случае даётся формулой:

Z

(о)

N

Z

1«

^1)

N

— dN

Для неразличимых частиц в случае отсутствия взаимодействия между частица-

Г7(°)

ми статистическая сумма класса перестановки ) распадается на произведение статистических сумм циклов. Статистическая сумма цикла длиной V выглядит следующим образом [7]:

Zv (Ь) = ^ (уЬ)

— d

Мы рассмотрели, во-первых, систему двух тождественных частиц с кулоновским отталкиванием в гармоническом внешнем поле с целью уточнить наши результаты,

полученные ранее. Во-вторых, были произведены расчёты для системы двух частиц в кулоновском внешнем поле. Центральный заряд был взят равным 2 (атом гелия). Мы построили две отличающиеся условиями замыкания цепочки гамильтонианов (13)

2 п т 2 п

Кг = -^ЕЕ(ГК_Г‘(''+1)) - - Е Е У ) ’

к=1 г=1 к=1 г=1

где Гк(^+1) = Гк1 для случая двух различимых частиц, гц^+1) = Г21, ^(N+1) = Гц для цикла из двух частиц. Число вершин п = 100. Для системы в гармоническом поле

V (Гкг, Кп) = ^|гн - Г2*| , |г1(*+1) “ г2(*+1) | , 9, ^ ,

в кулоновском внешнем поле

V (Гкг, ^ш) = (1 — ^ш) Г1 +

+

УГ - Ы , |г!(*+1) -Г2(*+1)| ’9’ _ Ш (|Г^1 ’ |Гй(;+1)| >0> ^

МК-блуждание в расширенном ансамбле было организовано следующим образом. Для изменения конформации использовались простые шаги: смещение одной случайно выбранной вершины - каждая компонента вектора смещения равномерно выбиралась из отрезка [—]х,]х]. Таким образом, вектор смещения равномерно выбирался из кубика с ребром 2]х. Попытка смены подансамбля предпринималась в среднем один раз на 100 шагов изменения конформации. Начальное значение величины Д$ бралось равным 0,01. В конце серии шагов величина Д$ домножалась на а = 0, 975. Такой выбор параметров был продиктован следующими соображениями. Использование больших значений Д$ приводило к очень большим флуктуациям результатов расчёта и выливалось в систематические ошибки, использование меньших значений приводило к ухудшению сходимости. Кроме уменьшения добавки Д$, в конце каждой серии масштаб шага ]х корректировался для поддержания процента удачных переходов в пределах 50 — 70 %.

В случае системы в гармоническом внешнем поле статистические суммы опорной и исследуемой систем отличались всего на 2 порядка, поэтому расширенный ансамбль был составлен из 10 подансамблей. Длина расчётов составляла 2000 серий по 1 000 000 шагов.

На рис. 1 представлены зависимости от обратной температуры статистических сумм для систем двух различимых частиц, цикла из двух частиц и тождественных частиц с нулевым спином (в = 0). При достаточно низких температурах вклады, соответствующие возбуждённым состояниям, становятся пренебрежимо малы по сравнению с вкладом, соответствующим основному состоянию. Исходя из этого, оставляя в статистической сумме только одно слагаемое, при помощи соотношения Гиббса-Гельмгольца Е = —д 1п Л/дв, мы можем получить значение энергии основного состояния. Пунктирной линией показана прямая, проведённая по методу наименьших квадратов через точки, соответствующие логарифму статистической суммы системы тождественных частиц, взятому со знаком минус. Значение энергии, полученное этим способом, оказывается равным 2, 70 ± 0,02, что гораздо лучше согласуется с данными, полученными в [10] (Е = 2, 67), чем результаты, полученные нами ранее в [9]. После публикации работы [10], её авторами был получен уточнённый результат для энергии Е = 2, 71, что практически совпадает в пределах погрешности с нашим результатом.

Рис. 1. Зависимость статистических сумм от обратной температуры: для системы двух различимых частиц (квадраты), цикла из двух частиц (круги) и системы тождественных частиц со спином равным нулю (треугольники) с кулоновским отталкиванием в одномерном гармоническом поле; пунктирная прямая проведена по методу наименьших квадратов через точки, соответствующие системе тождественных частиц; значение энергии, полученное по наклону, Е = 2,70 ±0,02

10-2 10- 3

Б 10- 4

£

и 10- 5 § 10- 6 1ю- . |ю- 8

н

и

10-

10-

'-а

4 5 6 7

обратная температура Ь

1

8

3

Рис. 2. Нормированное число посещений подансамблей для перехода от невзаимодействующих различимых частиц в гармоническом поле к частицам с кулоновским отталкиванием в кулоновском поле:

обратная температура Ь = 8; горизонталь показывает идеальное распределение

°,

н

е

я 0

е

о

о

по 0, о 0,

5

я ^ 0, е о н

3 0, а

ЕС

о

р

я 0, м р о н

0,

0475

0470

0465

0460

0455

0450

0445

0,0 0,2 0,4 0,6 0,8 1,0 X

В расчётах для трёхмерной системы двух тождественных частиц с кулоновским отталкиванием в гармоническом поле значение энергии получилось равным 4,55 ± 0,05, что также несколько меньше значения, полученного нами ранее (Е = 4, 73 [9]). Метод, использованный в [10], не пригоден для исследования трёхмерных систем, поэтому нам было не с чем сравнивать результат.

В случае включения кулоновского поля значение статистической суммы существенно возрастает по отношению к статистической сумме невзаимодействующих частиц в гармоническом поле. Уже при Ь = 4 различие между ними достигает 10 порядков. Это существенное препятствие, тем не менее, алгоритм Ванга-Ландау оказывается достаточно мощным методом, чтобы, увеличив число подансамблей до 22, нам удалось получить хорошие результаты и в этом случае. Длина расчётов составляла 2000 серий по 5 000 000 шагов.

На рис. 2 представлена гистограмма посещений подансамблей для перехода от невзаимодействующих различимых частиц в гармоническом поле (Хо = 0) к частицам с кулоновским отталкиванием в кулоновском поле (Хм = 1). Отличие полученных оценок вероятностей посещения от идеального равномерного распределения не более чем на 3 %.

На рис. 3 приводятся зависимости статистических сумм подансамблей для двух имеющихся классов перестановки от обратной температуры. Из 22 кривых демонстрируются 7. Кружками обозначены кривые, соответствующие классу двух различимых

1012

109

106

103

кГ

я 1

кГ 10- 3

10- 6

10- 9

10- 12

0,857

0,714

0,571

0,429

0,286

0,143

0

2 3 4 5 6 7

обратная температура Ь

8 9

Рис. 3. Зависимости статистических сумм подансамблей для классов перестановок от обратной температуры:

класс 12 — кружки, класс 2 — кресты; переход от двух невзаимодействующих частиц в гармоническом поле (Яо = 0) к системе двух частиц с кулоновским отталкиванием в кулоновском поле (Ям = 1); значения Я для промежуточных подансамблей указаны на рисунке

0

1

обратная температура Ь

Рис. 4. Зависимости статистических сумм от обратной температуры:

для системы двух различных частиц (треугольники), цикла из двух частиц (перевёрнутые треугольники), системы тождественных частиц со спином, равным нулю (круги) с кулоновским отталкиванием в трёхмерном кулоновском поле заряда 2 (атом гелия); пунктирными линиями показаны асимптоты, соответствующие энергиям Е = —2,904 основного состояния и Е = = —2,175 триплетного состояния атома гелия [12]

частиц, крестами - классу цикла из двух частиц. Каждая пара кривых отмечена соответствующим значением параметра Ят.

Далее, на рис. 4 представлены зависимости от обратной температуры статистических сумм систем двух различимых частиц, цикла из двух частиц и системы тождественных частиц без спина. Пунктирными линиями показаны асимптоты, соответствующие значениям энергий Е = -2, 904 основного состояния и Е = -2,175 триплетного состояния атома гелия, полученных путём квантовомеханических расчётов [12]. Как видно из рисунка, точки хорошо ложатся на асимптоты.

На рис. 5 изображена зависимость от температуры взятого со знаком минус логарифма статистической суммы системы тождественных частиц со спином равным 0 с кулоновским отталкиванием в кулоновском поле (триплетное состояние атома гелия) и системы со спином равным 1/2 (основное состояние атома гелия). Прямые линии проведены по методу наименьших квадратов. Из рис. 4 ясно видно, что статистические суммы при обратной температуре меньше Ь = 5 содержат вклады, соответствующие возбуждённым состояниям и не могут использованы для построения асимптоты. По наклону этих прямых можно определить значения энергий, которые получились равными -2, 3 ± 0,1 для в = 0 (Е = -2,175 [12]) и -2, 83 ± 0,05 для в = 1/2 (Е = -2, 904 [12]), что хорошо согласуется с результатами квантовомеханических расчётов.

Рис. 5. Зависимость логарифма статистической суммы:

для системы тождественных частиц со спином равным 0 (триплетное состояние атома гелия — кружки) и со спином равным 1/2 (основное состояние — квадраты) с кулоновским отталкиванием в куло-новском поле от обратной температуры; прямые линии проведены по методу наименьших квадратов

N

£

обратная температура b

Заключение. В настоящей работе метод расширенных ансамблей с настройкой балансировочных параметров по алгоритму Ванга-Ландау был использован для исследования равновесных свойств малых квантовых кулоновских систем. В начале были проведены расчёты для системы двух частиц с кулоновским отталкиванием в гармоническом поле. Полученные результаты очень хорошо согласуются с результатами численного интегрирования конфигурационных интегралов. Далее, было проведено моделирование системы двух взаимодействующих частиц в кулоновском поле (атом гелия). Полученные значения энергий хорошо согласуются с экспериментальными данными. Совокупность методов, предложенная в данной статье, может быть применена для исследования систем из большего числа квантовых частиц. С увеличением числа частиц, число классов перестановки тоже возрастает. Но поскольку расчёт для каждого класса проводится независимо от остальных, в подобных расчётах представляется весьма перспективным использование вычислительных grid-систем. Все расчёты были проведены на ПК Пентиум 4,3 ГГц. Время расчёта для одного значения температуры составляло от нескольких часов до суток.

Литература

1. Фейнман Р., Хиббс А. Квантовая механика и интегралы по траекториям / пер. с англ. М., 1968.

2. Takahashi M., Imada M. Monte Carlo calculation of quantum systems // J. Phys. Soc. Jpn. 1984. Vol. 53. N 3. P. 963-974.

3. Newman W. H., Kuki A. Improved method for path integral Monte Carlo integration in fermionic systems // J. Chem. Phys. 1992. Vol. 92. N 2. P. 1409-1417.

4. Lubartsev A. P., Vorontsov-Velyaminov P. N. Path integral Monte Carlo method in quantum statistics for a system of N identical fermions // Phys. Rev. (A). 1993. Vol. 48. N 6. P. 4075-4083.

5. Фейнман Р. Статистическая механика. Курс лекций / пер. с англ. М., 1975.

6. Vorontsov-Velyaminov P. N., Nesvit M. O., Gorbunov R. I. Bead-Fourier path-integral Monte Carlo method applied to systems of indentical particles // Phys. Rev. (E). 1997. Vol. 55. N 2. P. 1979-1997.

7. Lyubartsev A. P., Martsinovski A. A., Shevkunov S. V., Vorontsov-Velyaminov P. N. New approach to Monte Carlo calculation of the free energy: mеthod of expanded ensembles // J. Chem. Phys. 1992. Vol. 96. N 3. P. 1776-1783.

8. Wang F., Landau D. P. Efficient, multiple-range random walk algorithm to calculate the density of states // Phys. Rev. Lett. 2001. Vol. 86. N 10. P. 2050-2053.

9. Вознесенский М. А., Воронцов-Вельяминов П. Н., Любарцев А. П. Расчёты равновесных свойств квантовых систем методом Монте-Карло в расширенном ансамбле // Вычислит. методы и программир. 2008. Т. 9. C. 170-183.

10. Поляков Е. А., Воронцов-Вельяминов П. Н. Квантовый газ во внешнем поле при конечных температурах. Точное выражение для плотности и возбуждённых состояний // Там же. 2007. Т. 8. C. 334-351.

11. Kole J. S., De Raedt H. Quantum Monte Carlo method for attractive Coulomb potentials // Phys. Rev. (E). 2001. Vol. 64. P. 016704-(1)-016704-(6).

12. Pekeris C. L. I1 S, 21 S and 23S states of H~ and of He // Phys. Rev. 1962. Vol. 126. N 4. P. 1470-1476.

Принято к публикации 14 октября 2008 г.