Диссертация (1103257), страница 3
Текст из файла (страница 3)
Теоретически оптимальным является алгоритм тепловой бани("heat bath"), в котором конфигурации генерируются с вероятностьюT (x → x0 ) = π(x0 ).(36)A(x → x0 ) = 1,(37)P (x → x0 ) = π(x0 ),(38)При этомто есть все конфигурации принимаются, что обеспечивает наименьшую корреляцию между конфигурациями и наибыстрейшую термализацию.В реальности сгенерировать сразу всю конфигурацию с распределением (38) фактически невозможно, так оно будет иметь весьма сложный вид, а генерация случайных чисел с таким распределением – крайне долгая операция. Однако, если рассмотретьx = (x1 , . . . , xn−1 , xn , xn+1 , . .
. , xN ) и x0 = (x1 , . . . , xn−1 , x0n , xn+1 , . . . , xN ), то можно всевышеизложенные рассуждения применить к xn , считая его конфигурацией с заданнымиграничными условиями; изменение при этом нормировки в (26) очевидно не отразится в(35). ТогдаP (x → x0 ) = P (x → x0 ) = π(x0 ) ∼ e−Sn /~ ,(39)где Sn – часть действия, содержащая все слагаемые, содержащие xn и зависящая от остальных координат как от параметров.
В таком случае процесс генерации состоит в последовательном проходе ("sweep") всех координат с правилом (39), после чего новая конфигурация добавляется в набор.Метод тепловой бани хорошо работает тогда, когда можно эффективно генерироватьраспределение (39). Во многих случаях это невозможно и приходится пользоваться общей схемой. При этом по-прежнему конфигурацию разумно генерировать не полностью, апроходом по временным слоям.
Алгоритм полностью определяется выбором T (x → x0 ), откоторого соответственно зависит производительность. Для большей эффективности пробное распределение должно отвечать двум требованиям: по математическим соображенияоно должно быть близко к π(x0 ), а по вычислительным – допускать быструю генерациюсоответствующих случайных чисел. Исторически первым было предложено использоватьравномерное распределение; оно очевидно удовлетворяет второму требованию и очевидноне удовлетворяет первому.Представляется разумным (и используется в данной работе) пробное распределение,состоящее из кинетической части π(x0 ). Оно имеет вид (см. (39),(23))( 0 )xn−1 +xn+1 2−mxn2T (xn → x0n ) ∼ exp −,a~(40)то есть является гауссовым.
Для генерации гауссова распределения существует крайнеэффективный алгоритм – преобразование Бокса-Мюллера позволяет из двух равномернораспределенных случайных чисел (генератор есть в стандартных библиотеках, но в даннойработе используется и всем рекомендуется "вихрь Мерсенна"из GNU GSL) получить дванезависимых нормально распределенных случайных числа:)pγr ∈ U (0, 1) → r = −2σ 2 ln γrξ1 = r cos ϕ + µ→∈ N (µ, σ 2 ).γϕ ∈ U [0, 1) → ϕ = 2πγϕξ2 = r sin ϕ + µДля пробного распределения (40) вероятность принятия конфигурации имеет вид −V (x0 )a/~ ne0A(xn → xn ) = min 1, −V (xn )a/~ .e(41)(42)Из сравнения (40) и (42) ясно, что при достаточно малом шаге наиболее изменчивая частьраспределения учитывается при генерации нового значения при помощи (40), а вероятность принятия этого значения (42) близка к единице.Так как в вероятность принятия изменения конфигурации в общем случае меньшеединицы, последовательные конфигурации частично совпадают, что не позволяет считатьиз независимыми.
Для уменьшения корреляции можно последовательное примененитьправило (34) несколько раз. Рассматривая для простоты две попытки (далее возможноиндуктивное рассуждение), несложно убедиться, что вероятность переходаRP (x → x00 ) = T (x → x0 )A(x → x0 ) + δ(x − x0 ) 1 − dyT (x → y)A(x → y) ×R× T (x0 → x00 )A(x0 → x00 ) + δ(x0 − x00 ) 1 − dyT (x0 → y)A(x0 → y) ,(43)то есть применение правила (34) фиксированное число раз, удовлетворяет условию детального баланса. Надо отметить, что этому условию не удовлетворяет вероятность переходаRP (x → x00 ) = T (x → x0 )A(x → x0 )δ(x0 − x00 ) + δ(x − x0 ) 1 − dyT (x → y)A(x → y) ×R× T (x0 → x00 )A(x0 → x00 ) + δ(x0 − x00 ) 1 − dyT (x0 → y)A(x0 → y) ,(44)то есть принятие первого "удачного"изменения конфигурации за окончательное, хотя иускоряет вычисления, но не гарантирует правильной термализации.
Также теоретическинеобоснованна генерация до первого удачного значения без ограничения количества попыток, то есть вероятность переходаP (x → x0 ) = RT (x → x0 )A(x → x0 ),dyT (x → y)A(x → y)(45)Впрочем, последние два алгоритма в численных экспериментах вели к правильной термализации (что не противоречит теории, так как условие баланса не является необходимым).Правильный результат вычислений, разумеется, получается при a → 0, кроме того, какуже говорилось, уменьшение шага увеличивает вероятность принятия новых конфигураций. Однако при этом время вычислений увеличивается за счет увеличения числа шагов(так как физическое время N a должно сохраняться), а последовательные конфигурациименьше отличаются, то есть скоррелированы. Далее будет изложен метод, позволяющийбыстрее получать более качественные (с меньшим шагом и меньшей корреляцией) конфигурации.3.2Многоуровневый алгоритмОсновная идеяданного метода состоит в том, что участок конфигурации изменяется непоследовательным проходом его точек, а последовательным приближением с улучшениемкачества траектории на нем[3].
Более грубое приближение генерируется быстро, но позволяет достаточно хорошо оценить статистический вес и принять или отклонить траекториюеще до того, как она полностью сформирована. Будем рассматривать алгоритм в с выделением уровней методом деления пополам, так как он является одновременно наиболееэффективным и простым для понимания и реализации.Выберем число уровней алгоритма Nlevel ≤ log2 N (конкретное значение подбирается для каждой задачи эмпирически, с целью ускорения счета). Случайным образом выберем участок длиной 2Nlevel шагов. В силу периодических граничных условий можнорассмотреть любой такой участок, например s = (x0 , .
. . , x2Nlevel ). Разобьем его на уровни: s0 = (x0 , x2Nlevel ) – граничные точки – 0-й уровень, s1 = (x2Nlevel −1 ) – 1-й, s2 =(x2Nlevel −2 , x2Nlevel −1 +2Nlevel −2 ) – 2-й, и так далее, sNlevael = (x1 , x3 , . . . , x2Nlevel −1 ) – Nlevel -й;на k-m уровне 2k−1 точек. В качестве примера на (Рис. 1) показано разбиение участкатраектории на четыре уровня.20121x0-1-2-2024681012141618nРис. 1: Последовательность генерации участка траектории в многоуровневом алгоритме(Nlevel = 4)Теперь построим "действие уровня"πk (sk ) ≡ πk (s0 , .
. . , sk−1 , sk ), где координаты напредыдущих уровнях являются параметрами. Важно отметить, что данная функция необязана совпадать с действием и может быть выбрана произвольно, исходя и вычислительных соображений. Единственным требованием является совпадение действия последнегоуровня с истинным:πNlevel (sNlevel ) = π(s).(46)Впрочем, оптимальным выбором является именно истинное действие, проинтегрированноепо координатам следующих уровней:Zdsk+1 . . . dsNlevel π(s).πk (sk ) =(47)Далее поочередно, начиная с первого уровня, применяется алгоритм, подобный описанному выше алгоритму Метрополиса-Гастингса.
При помощи пробного распределенияTk (s0k ) = Tk (s00 , . . . , s0k−1 ; sk ; sk+1 , . . . , sNlevel → s00 , . . . , s0k−1 ; s0k ; sk+1 , . . . , sNlevel )генерируются новые значения координат уровня s0k . С вероятностьюAk (s0k ) = A(s00 , . . . , s0k−1 ; sk ; sk+1 , . . . , sNlevel → s00 , . . . , s0k−1 ; s0k ; sk+1 , . . . , sNlevel ) =hiTk (sk )πk (s0k )πk−1 (sk )= min 1, Tk (s0 )πk (sk )πk−1 (s0 )k(48)kэти значения принимаются и генерация продолжается; в противном случае генерация прекращается и за новую принимается конфигурация со старыми координатами на всех уровнях.Таким образом, вероятность перехода для имеет видPk (s0k ) = P (s00 , .
. . , s0k−1 ; sk ; sk+1 , . . . , sNlevel → s00 , . . . , s0k−1 ; s0k ; sk+1 , . . . , sNlevel ) == Tk (s00 , . . . , s0k−1 ; sk ; sk+1 , . . . , sNlevel → s00 , . . . , s0k−1 ; s0k ; sk+1 , . . . , sNlevel ) ×h T (s0 ,...,s ;s0 ;s0 ,...,s0i00000kk−1Nlevel →s0 ,...,sk−1 ;sk ;sk+1 ,...,sNlevel )πk (s0 ,...,sk−1 ,sk )πk−1 (s0 ,...,sk−1 )k+1× min 1, Tk (s0 ,...,s0 ;skk ;sk+1+00000,...,sNlevel →s0 ,...,sk−1 ;sk ;sk+1 ,...,sNlevel )πk (s0 ,...,sk−1 ,sk )πk−1 (s0 ,...,sk−1 )0k−1R+δ(s − s0 ) 1 − dyTk (s00 , . . . , s0k−1 ; sk ; sk+1 , .
. . , sNlevel → s00 , . . . , s0k−1 ; y; sk+1 , . . . , sNlevel )×h T (s0 ,...,s ;y;s0 ,...,s0io0000kk−1Nlevel →s0 ,...,sk−1 ;sk ;sk+1 ,...,sNlevel )πk (s0 ,...,sk−1 ,y)πk−1 (s0 ,...,sk−1 )k+1× min 1, Tk (s0 ,...,s0 ;sk ;sk+1(49),...,sN→s0 ,...,s0 ;s0 ;sk+1 ,...,sN)πk (s0 ,...,sk−1 ,sk )πk−1 (s0 ,...,s0 )0k−1level0k−1klevel0k−1и пробных распределениях Tk (s0k ) удовлетворяет условию детального баланса для данногоуровня:Pk (s0k )πk (sk )πk (s0k )= Pk (sk ).πk−1 (sk−1 )πk−1 (s0k−1 )(50)Учитывая неизменность граничных точек s0 = s00π0 (s0 ) = π0 (s00 ),(51)перемножая уравнения (50) для всех уровнейπ (s0 )1)P1 (s01 ) ππ10 (s= P1 (s1 ) π01 (s01 ) ,(s0 )02)P2 (s02 ) ππ21 (s(s1 )=π (s0 )P2 (s2 ) π12 (s02 ) ,1...πNlevel (sNlevel )(s)level −1 Nlevel −1PNlevel (s0Nlevel ) πNπNlevel (s0N= PNlevel (sNlevel ) πN(s0level −1 Nlevel)level −1),(52)используя очевидную формулу для полной вероятности перехода0P (s ) =NYlevelPk (sk ),(53)k=1и условие (46), получаем, что детальный баланс на каждом уровне ведет к детальномубалансу на всей траектории:P (s0 )π(s) = P (s)π(s0 ),(54)что обеспечивает корректность рассматриваемого алгоритма.За счет отсева конфигураций с низким статистическим весом на ранних этапах и одновременного изменения участка траектории многоуровневый алгоритм дает значительноеувеличение скорости генерации нескоррелированных конфигураций.
Конкретный вид егозадается функциями Tk (s0k ) и πk (sk ). Так как точки одного уровня не являются соседними,то в каждой из них можно использовать пробное распределение (40), разумеется положивв нем xn−1 → xn−2Nlevel −k , xn+1 → xn+2Nlevel −k , a → ak = 2Nlevel −k a. Тогда, выбирая действиеуровня в соответствии с (135), получаем вероятность принятия новых координат уровня −V (s0 )a/~ ke0(55)Ak (sk ) = min 1, −V (s )a/~ ,keгде a – "настоящий"шаг (замены, как в Tk (s0k ) делать не нужно), а потенциал вычисляетсятолько на данном уровне. Алгоритм, в котором действие уровня включает полный шагуровня, также возможен и предложен, в частности в работе [3].
Однако для расмматриваемой задачи о системе многих тел он оказывается значительно мнее производительным,чем выше описанный алгоритм, используемый в данной работе. Наш вариант отличаетсябольшей средней вероятностью принятия конфигурации (acceptance rate) и достаточнобыстрой термализацией в большем диапазоне параметров.3.3Измерения. Среднее значение кинетической энергииФормула (29) позволяет вычислить среднее значение любой величины, представленной ввиде функции координат. Например, несложно измерить hxi, hx2 i, hx4 i, hx0 xn i. Знаниекоррелятора hx0 xn i позволяет узнать энергию низшего возбужденного состояния ∆E =E1 − E0 .
Для этого используется соотношение[1]hx(τ )x(τ + ∆τ )i − hx(τ )ihx(τ + ∆τ )i ∼ exp {−∆E∆τ /~} ,(56)которое на решетке (предполагая, что hxi = 0) переходит вhx0 xn i ∼ exp {−∆Ena/~} .(57)~∆E = − lc−1 ,a(58)Отсюдагде введено обозначение для обратной корреляционной длины, измеренной в шагах решетки:lc−1 =∂ lnhx0 xn i.∂n(59)В задачах квантовой механики несомненно имеет смысл и является важнейшей задачейвычисление энергии основного состояния.
Определение средней потенциальной энергиине представляет проблемы, так как она очевидным образом представляется как функциякоординат и для нее можно применить формулу (29). Среднюю кинетическую энергиюможно вычислить при помощи теоремы вириала:hKi = hxV 0 (x)i/2.(60)Применение теоремы вириала в сложных задачах квантовой механики, например взадаче многих тел, может быть неудобно, а в задачах квантовой теории поля и вовсеневозможно.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
