Алгоритм Ванга и Ландау

Алгоритм Ванга и Ландау , предложенный Фугао Вангом и Дэвидом П. Ландау , ^[1] представляет собой метод Монте-Карло , предназначенный для оценки плотности состояний системы. Метод выполняет немарковское случайное блуждание для построения плотности состояний, быстро посещая весь доступный энергетический спектр. Алгоритм Ванга и Ландау — важный метод получения плотности состояний, необходимой для выполнения мультиканонического моделирования .

Алгоритм Ванга – Ландау можно применить к любой системе, которая характеризуется функцией стоимости (или энергии). Например, его применяли для решения числовых интегралов ^[2] и сворачивания белков. ^{[3] [4]} Выборка Ванга – Ландау связана с алгоритмом метадинамики . ^[5]

Обзор

Алгоритм Ванга и Ландау используется для получения оценки плотности состояний системы, характеризуемой функцией стоимости. Он использует немарковский случайный процесс , который асимптотически сходится к мультиканоническому ансамблю . ^[1] (То есть к алгоритму Метрополиса – Гастингса с распределением выборки, обратным плотности состояний) Основным следствием является то, что такое распределение выборки приводит к моделированию, в котором энергетические барьеры невидимы. Это означает, что алгоритм посещает все доступные состояния (благоприятные и менее благоприятные) гораздо быстрее, чем алгоритм Метрополиса. ^[6]

Алгоритм

Рассмотрим систему, определенную в фазовом пространстве , и функцию стоимости E (например, энергию), ограниченную спектром , которому соответствует плотность состояний , которую необходимо оценить. Оценщик есть . _ Поскольку алгоритм Ванга и Ландау работает в дискретных спектрах, ^[1] спектр делится на N дискретных значений с разницей между ними , таких что ${\displaystyle \Omega }$ ${\displaystyle E\in \Gamma =[E_{\min },E_{\max }]}$ ${\displaystyle \rho (E)}$ ${\displaystyle {\hat {\rho }}(E)\equiv \exp(S(E))}$ ${\displaystyle \Gamma }$ ${\displaystyle \Delta }$

${\displaystyle N={\frac {E_{\max }-E_{\min }}{\Delta }}}$.

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

• обнуление всех записей микроканонической энтропии, ${\displaystyle S(E_{i})=0\ \ i=1,2,...,N}$
• инициализация и ${\displaystyle f=1}$
• случайная инициализация системы путем ввода случайной конфигурации . ${\displaystyle {\boldsymbol {r}}\in \Omega }$

Затем алгоритм выполняет мультиканоническое ансамблевое моделирование: ^[1] случайное блуждание Метрополиса – Гастингса в фазовом пространстве системы с распределением вероятностей, заданным и вероятностью предложения нового состояния, заданной распределением вероятностей . Сохраняется гистограмма посещенных энергий. Как и в алгоритме Метрополиса–Гастингса, выполняется этап принятия предложения и заключается в (см. обзор алгоритма Метрополиса–Гастингса ): ${\displaystyle P({\boldsymbol {r}})=1/{\hat {\rho }}(E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))}$ ${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$ ${\displaystyle H(E)}$

1. предложение штата в соответствии с произвольным распределением предложений ${\displaystyle {\boldsymbol {r}}'\in \Omega }$ ${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$
2. принять/отклонить предложенное состояние в соответствии с

${\displaystyle A({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')=\min \left(1,e^{S-S'}{\frac {g({\boldsymbol {r}}'\rightarrow {\boldsymbol {r}})}{g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}}\right)}$

где и . ${\displaystyle S=S(E({\boldsymbol {r}}))}$ ${\displaystyle S'=S(E({\boldsymbol {r}}'))}$

После каждого шага принятия предложения система переходит к некоторому значению , увеличивается на единицу и выполняется следующее обновление: ${\displaystyle E_{i}}$ ${\displaystyle H(E_{i})}$

${\displaystyle S(E_{i})\leftarrow S(E_{i})+f}$.

Это решающий шаг алгоритма, и именно он делает алгоритм Ванга и Ландау немарковским: теперь случайный процесс зависит от истории процесса. Следовательно, в следующий раз, когда государству с такой же энергией поступит предложение , оно, скорее всего, будет отклонено; в этом смысле алгоритм заставляет систему одинаково посещать весь спектр. ^[1] В результате гистограмма становится все более плоской. Однако эта неравномерность зависит от того, насколько хорошо вычисленная энтропия приближена к точной энтропии, которая, естественно, зависит от значения f. ^[7] Чтобы все лучше и лучше аппроксимировать точную энтропию (и, следовательно, плоскостность гистограммы), f уменьшается после M шагов принятия предложения: ${\displaystyle E_{i}}$ ${\displaystyle H(E)}$

${\displaystyle f\leftarrow f/2}$.

Позже было показано, что обновление f путем постоянного деления на два может привести к ошибкам насыщения. ^[7] Небольшая модификация метода Ванга и Ландау, позволяющая избежать этой проблемы, заключается в использовании коэффициента f , пропорционального , где пропорционально количеству шагов моделирования. ^[7] ${\displaystyle 1/t}$ ${\displaystyle t}$

Тестовая система

Мы хотим получить DOS для потенциала гармонического осциллятора .

${\displaystyle E(x)=x^{2},\,}$

Аналитический DOS определяется выражением:

${\displaystyle \rho (E)=\int \delta (E(x)-E)\,dx=\int \delta (x^{2}-E)\,dx,}$

выполнив последний интеграл, получим

${\displaystyle \rho (E)\propto {\begin{cases}E^{-1/2},{\text{for }}x\in \mathbb {R} ^{1}\\{\text{const}},{\text{for }}x\in \mathbb {R} ^{2}\\E^{1/2},{\text{for }}x\in \mathbb {R} ^{3}\\\end{cases}}}$

В общем, DOS для многомерного гармонического осциллятора будет задаваться некоторой степенью E , показатель степени будет функцией размерности системы.

Следовательно, мы можем использовать простой потенциал гармонического осциллятора для проверки точности алгоритма Ванга – Ландау, поскольку мы уже знаем аналитическую форму плотности состояний. Поэтому мы сравниваем оценку плотности состояний, полученную алгоритмом Ванга–Ландау, с . ${\displaystyle {\hat {\rho }}(E)}$ ${\displaystyle \rho (E)}$

Образец кода

Ниже приведен пример кода алгоритма Ванга-Ландау на Python , где мы предполагаем, что используется симметричное распределение предложений g:

${\displaystyle g({\boldsymbol {x}}'\rightarrow {\boldsymbol {x}})=g({\boldsymbol {x}}\rightarrow {\boldsymbol {x}}')}$

Код рассматривает «систему», которая является базовой изучаемой системой.

текущаяЭнергия  =  система . Случайная конфигурация ()  # Случайная начальная конфигурацияв то время как  f  >  эпсилон :  система . OfferConfiguration ()  # Предлагается предлагаемая конфигурация предложенная  Энергия  =  система . предложенная энергия ()  # Вычисленная энергия предлагаемой конфигурации if  random ()  <  exp ( entropy [ currentEnergy ]  -  entropy [ предложенная энергия ]):  # Если принято, обновите энергию и систему:  currentEnergy  =  предложенная энергетическая  система . AcceptProposeConfiguration ()  else :  # Если  система отклонена . отклонить предложенную конфигурацию () H [ currentEnergy ]  +=  1  энтропия [ currentEnergy ]  +=  f if  isFlat ( H ):  # isFlat проверяет, является ли гистограмма плоской (например, плоскостность 95%)  H [:]  =  0  f  *=  0,5  # Уточняем параметр f

Молекулярная динамика Ванга и Ландау: статистическая температурная молекулярная динамика (STMD)

Молекулярная динамика (МД) обычно предпочтительнее метода Монте-Карло (МК), поэтому желательно иметь алгоритм МД, включающий основную идею WL для выборки с плоской энергией. Этот алгоритм — статистическая температурная молекулярная динамика (STMD), разработанный ^[8] Джагилом Кимом и др. в Бостонском университете.

Важный первый шаг был сделан с помощью алгоритма статистической температуры Монте-Карло (STMC). WLMC требует значительного увеличения количества энергетических ячеек в зависимости от размера системы, что вызвано непосредственной работой с плотностью состояний. STMC сосредоточен на интенсивной величине, статистической температуре, где E - потенциальная энергия. В сочетании с соотношением , где мы установили , правило WL для обновления плотности состояний дает правило для обновления дискретизированной статистической температуры, ${\displaystyle T(E)=1/(dS(E)/dE)}$ ${\displaystyle \Omega (E)=e^{S(E)}}$ ${\displaystyle k_{B}=1}$

${\displaystyle {\tilde {T}}_{j\pm 1}^{\prime }=\alpha _{j\pm 1}{\tilde {T}}_{j\pm 1},}$

где где - размер элемента энергии и обозначает текущую оценку. Мы определяем f как в ^[1] коэффициент >1, который умножает оценку DOS для i-го интервала энергии, когда система посещает энергию в этом интервале. ${\displaystyle \alpha _{j\pm 1}=1/(1\mp \delta f{\tilde {T}}_{j\pm 1}),\delta f=(lnf/2\Delta E),\Delta E}$ ${\displaystyle {\tilde {T}}}$

Подробности приведены в Ref. ^[8] При первоначальном предположении для и диапазоне, ограниченном между и , моделирование происходит так же, как в WLMC, со значительными численными различиями. Интерполяция дает континуальное выражение оценки после интегрирования ее обратного значения, что позволяет использовать более крупные элементы разрешения по энергии, чем в WL. При оценке вероятности принятия в пределах одного и того же энергетического интервала доступны разные значения . Когда колебания гистограммы составляют менее 20% от среднего значения, значение уменьшается в соответствии с . ${\displaystyle T(E)}$ ${\displaystyle T_{L}}$ ${\displaystyle T_{U}}$ ${\displaystyle {\tilde {T}}(E)}$ ${\displaystyle S(E)}$ ${\displaystyle S(E)}$ ${\displaystyle f}$ ${\displaystyle f\rightarrow {\sqrt {f}}}$

STMC сравнивали с WL для модели Изинга и жидкости Леннарда-Джонса. При увеличении размера энергетического бункера STMC получает те же результаты в значительном диапазоне, в то время как производительность WL быстро ухудшается. STMD может использовать меньшие начальные значения для более быстрой сходимости. В целом, STMC требуется меньше шагов для получения результатов того же качества. ${\displaystyle f_{d}=f-1}$

Теперь рассмотрим основной результат — STMD. Он основан на наблюдении, что в стандартном МД-моделировании при температуре с силами, полученными из потенциальной энергии , где обозначает все положения, вес выборки для конфигурации равен . Более того, если силы получены из функции , вес выборки равен . ${\displaystyle T_{0}}$ ${\displaystyle E([x])}$ ${\displaystyle [x]}$ ${\displaystyle e^{-E([x])/T_{0}}}$ ${\displaystyle W(E)}$ ${\displaystyle e^{-W(E([x]))/T_{0}}}$

Для плоской выборки энергии пусть эффективный потенциал будет энтропийной молекулярной динамикой. Тогда вес . Поскольку плотность состояний равна , их произведение дает плоскую выборку энергии. ${\displaystyle T_{0}S(E)}$ ${\displaystyle e^{-S(E)}}$ ${\displaystyle e^{+S(E)}}$

Силы рассчитываются как

${\displaystyle F=(-d/dx)T_{0}S(E)=T_{0}(dS/dE)(-d/dx)E([x])=(T_{0}/T(E))F^{0}}$

где обозначает обычную силу, полученную из потенциальной энергии. Масштабирование обычных сил с помощью коэффициента дает плоскую выборку энергии. ${\displaystyle F^{0}}$ ${\displaystyle (T_{0}/T(E))}$

STMD начинается с обычного алгоритма MD при константе и V. Силы масштабируются, как указано, а статистическая температура обновляется на каждом временном шаге с использованием той же процедуры, что и в STMC. Поскольку симуляция сходится к выборке плоской энергии, текущая оценка сходится к истинному значению . Технические детали, включая меры по ускорению конвергенции, описаны в ^[8] и. ^[9] ${\displaystyle T_{0}}$ ${\displaystyle {\tilde {T}}(E)}$ ${\displaystyle T(E)}$

В STMD температура называется кинетической, так как она контролирует скорости как обычно, но не входит в конфигурационную выборку, что необычно. Таким образом, STMD может исследовать низкие энергии с помощью быстрых частиц. Любое каноническое среднее значение можно рассчитать с помощью повторного взвешивания, но статистическую температуру можно получить сразу же, без дополнительного анализа. Это чрезвычайно ценно для изучения фазовых переходов. В конечных наносистемах имеется особенность, соответствующая каждому «субфазовому переходу». Для достаточно сильного перехода равновеликая конструкция на S-петле дает температуру перехода. ${\displaystyle T_{0}}$ ${\displaystyle T(E)}$ ${\displaystyle T(E)}$ ${\displaystyle 1/T(E)}$

STMD был усовершенствован группой BU ^[9] и применен ими и другими к нескольким системам. Д. Стелтер признавал, что, несмотря на наш упор на работу с интенсивными объемами, она является обширной. Однако это трудоемкий процесс, и процедура, основанная на неравномерности гистограммы, заменяется сокращением пополам через каждое фиксированное количество временных шагов. Это простое изменение делает STMD полностью интенсивным и существенно повышает производительность больших систем. ^[9] Кроме того, конечное значение интенсивности является константой, определяющей величину ошибки в конвергентной системе и не зависящей от размера системы. STMD реализован в LAMMPS как исправление stmd. ${\displaystyle ln(f)}$ ${\displaystyle \delta f=(ln(f)/2\Delta E)}$ ${\displaystyle f\rightarrow {\sqrt {f}}}$ ${\displaystyle \delta f}$ ${\displaystyle \delta f}$ ${\displaystyle T(E)}$

STMD особенно полезен для фазовых переходов. Информацию о равновесии невозможно получить с помощью канонического моделирования, поскольку для возникновения перехода необходимо переохлаждение или перегрев. Однако прогон STMD позволяет получить плоскую выборку энергии с естественным развитием нагрева и охлаждения, не попадая в состояние низкой или высокой энергии. Совсем недавно он был применен к переходу жидкость/гель ^[9] в наночастицах, обернутых липидами.

Обмен репликами STMD ^[10] также был представлен группой BU.

Рекомендации

1. Ван, Фугао и Ландау, ДП (март 2001 г.). «Эффективный многодиапазонный алгоритм случайного блуждания для расчета плотности состояний». Физ. Преподобный Летт . 86 (10): 2050–2053. arXiv : cond-mat/0011174 . Бибкод : 2001PhRvL..86.2050W. doi :10.1103/PhysRevLett.86.2050. PMID  11289852. S2CID  2941153.
2. RE Белардинелли, С. Манзи и В.Д. Перейра (декабрь 2008 г.). «Анализ сходимости алгоритмов 1/t и Ванга – Ландау при вычислении многомерных интегралов». Физ. Преподобный Е. 78 (6): 067701. arXiv : 0806.0268 . Бибкод : 2008PhRvE..78f7701B. doi : 10.1103/PhysRevE.78.067701. PMID  19256982. S2CID  8645288.
3. П. Охеда, М. Гарсия, А. Лондоно и Нью-Йорк Чен (февраль 2009 г.). «Моделирование белков в клетках методом Монте-Карло: влияние заключения на стабильность промежуточных состояний». Биофиз. Дж . 96 (3): 1076–1082. arXiv : 0711.0916 . Бибкод : 2009BpJ....96.1076O. дои : 10.1529/biophysj.107.125369. ПМЦ 2716574 . ПМИД  18849410. 
4. П. Охеда и М. Гарсия (июль 2010 г.). «Вызванное электрическим полем разрушение нативной конформации белка бета-листа и образование альфа-спиральной структуры». Биофиз. Дж . 99 (2): 595–599. Бибкод : 2010BpJ....99..595O. дои : 10.1016/j.bpj.2010.04.040. ПМК 2905109 . ПМИД  20643079. 
5. Кристоф Юнгханс, Дэнни Перес и Томас Фогель. «Молекулярная динамика в мультиканоническом ансамбле: эквивалентность выборки Ванга – Ландау, статистическая температурная молекулярная динамика и метадинамика». Журнал химической теории и вычислений 10.5 (2014): 1843–1847. дои : 10.1021/ct500077d
6. Берг, Б.; Нейхаус, Т. (1992). «Мультиканонический ансамбль: новый подход к моделированию фазовых переходов первого рода». Письма о физических отзывах . 68 (1): 9–12. arXiv : hep-lat/9202004 . Бибкод : 1992PhRvL..68....9B. doi :10.1103/PhysRevLett.68.9. PMID  10045099. S2CID  19478641.
7. Белардинелли, Р.Э. и Перейра, В.Д. (2007). «Алгоритм Ванга – Ландау: теоретический анализ насыщения ошибки». Журнал химической физики . 127 (18): 184105. arXiv : cond-mat/0702414 . Бибкод : 2007JChPh.127r4105B. дои : 10.1063/1.2803061. PMID  18020628. S2CID  25162388.
8. Ким, Джегил; Штрауб, Джон и Киз, Том (август 2006 г.). «Статистически-температурный метод Монте-Карло и алгоритмы молекулярной динамики». Физ. Преподобный Летт . 97 (5): 50601–50604. doi : 10.1103/PhysRevLett.97.050601.
9. Стелтер, Дэвид и Киз, Том (2019). «Моделирование фазового равновесия жидкость/гель в липидных везикулах». Мягкая материя . 15 : 8102–8112. дои : 10.1039/c9sm00854c.
10. Ким, Джегил; Штрауб, Джон и Киз, Том (апрель 2012 г.). «Алгоритм статистико-температурной молекулярной динамики репликационного обмена». Журнал физической химии Б. 116 : 8646–8653. дои : 10.1021/jp300366j.