Русская Википедия:Алгоритм Ванга — Ландау

Содержание

Текст

Текст

Текст

Текст

Текст

Текст

Текст

Текст

Текст

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

Алгоритм Ванга-Ландау может быть применен к любой системе, которая характеризуется некоторым параметром (например, энергией, объемом и др.). К примеру, он может быть использован для численного интегрирования и моделирования белков.

Описание

Алгоритм Ванга-Ландау является реализацией метода энтропического моделирования, в котором изучается плотность состояний с помощью блуждания в пространстве энергий с равновероятным посещением всех энергетических состояний. Алгоритм решает проблему подбора подходящих вероятностей перехода для получения требуемого при энтропическом моделировании равномерного посещения энергетических состояний и, следовательно, позволяет получить плотность состояний <math>\Omega(E)</math>.

Алгоритм

Рассмотрим систему в фазовом пространстве <math>\Omega</math> и энергию <math>E</math>, изменение которой ограничено диапазоном <math>E_{min} \leq E \leq E_{max}</math>. Пусть рассматриваемая система имеет плотность вероятности <math>\rho(E) \sim \exp(S(E))</math>, которую нам требуется посчитать. Поскольку алгоритм Ванга и Ландау работает с дискретным спектром энергии, диапазон <math>E_{min} \leq E \leq E_{max}</math> разбивается на конечное число (<math>M</math>) равных отрезков («ящиков»), размер которых равен <math>\Delta</math>. Таким образом:

<math>M = \frac{E_{max}-E_{min}}{\Delta}</math>.

С учетом этого дискретного спектра, алгоритм имеет следующие начальные условия:

• <math>S(E_i) = 0, i = 1, 2, ..., M</math>
• <math>f = 1</math> (используется далее, как добавка к энтропии <math>S(E_i)</math> после каждого принятого шага)
• Выбирается случайная конфигурация системы <math>\vec{r} \in \Omega</math>

Алгоритм выполняет моделирование в мультиканоническом ансамбле: случайное блуждание Метрополиса-Гастингса по фазовому пространству <math>\Omega</math> системы с распределением вероятности <math>P(\vec{r}) = 1/\rho(E(\vec{r})) = \exp(-S(E(\vec{r})))</math> и вероятностью генерации нового состояния, данной распределением вероятности <math>g(\vec{r} \rightarrow \vec{r}')</math>, которое выбирается произвольно (обычно любое состояние может быть сгенерировано с равной вероятностью). В процессе моделирования, посещение каждого «ящика» записывается в гистограмму <math>H(E)</math> (то есть, значение <math>H(E)</math> увеличивается на единицу). Как и в алгоритме Метрополиса-Гастингса, генерация и принятие нового состояния выполняется следующим образом:

1. генерация нового состояния <math>\vec{r}' \in \Omega</math> согласно распределению вероятности <math>g(\vec{r} \rightarrow \vec{r}')</math>
2. принятие/отклонение нового состояния, производится следующим образом:

Если энтропия нового состояния меньше текущего, то оно сразу принимается. Если же энтропия увеличилась, то новое состояние принимается с вероятностью: <math>A = e^{S - S'}\frac{g(\vec{r}'\rightarrow \vec{r})}{g(\vec{r}\rightarrow \vec{r}')}</math>, где <math>S = S(E(\vec{r}))</math> и <math>S' = S(E(\vec{r}'))</math>.

То есть, общая формула выглядит следующим образом:

<math>A(\vec{r}\rightarrow \vec{r}') = \min\left(1,e^{S - S'}\frac{g(\vec{r}'\rightarrow \vec{r})}{g(\vec{r}\rightarrow \vec{r}')}\right)</math>.

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

После каждого шага генерации-принятия система переходит в некоторое состояние <math>E_i</math>, значение <math>H(E_i)</math> увеличивается на единицу, а также выполняется следующее изменение:

<math>S(E_i) \leftarrow S(E_i) + f</math>

Это важный шаг алгоритма, и это то, что делает алгоритм Ванга-Ландау немарковским: случайный процесс теперь зависит от истории процесса. Таким образом, когда в следующий раз будет предложено состояние с энергией <math>E_i</math>, это состояние будет отклонено с большей вероятностью; в этом смысле, алгоритм принуждает систему посещать все состояния с одинаковой частотой. Как следствие, гистограмма <math>H(E)</math> становится все более и более плоской. Хотя, эта равномерность зависит от того, насколько посчитанная энтропия близка к точной энтропии, что зависит от <math>f</math>. Для улучшения приближения точной энтропии (и, таким образом, равномерности гистограммы), <math>f</math> уменьшается после <math>M</math> шагов генерации-принятия:

<math>f \leftarrow f/2</math>

Через некоторое время было показано, что при изменении <math>f</math>постоянным делением на два алгоритм может не сходиться. Небольшая модификация метода Ванга-Ландау позволяет избежать этого: производится деление не на два, а на <math>t</math>, при чем <math>t</math> пропорционально шагу моделирования.

В результате использования этого алгоритма происходит автоматическая настройка весов вероятности перехода, которые одновременно определяют плотности состояний. По окончании расчета вычисляется массив <math>\Omega(E) = \exp S(E)</math>и нормируется на единицу.

Пример кода

Ниже показан пример кода на Python, в котором предполагается симметричность функции распределения <math>g</math>:

<math>g(\boldsymbol{x}'\rightarrow \boldsymbol{x})=g(\boldsymbol{x}\rightarrow \boldsymbol{x}')</math>

currentEnergy = system.randomConfiguration() # случайная генерация начального состояния системы
while (f > epsilon):
    system.proposeConfiguration() # генерация новой конфигурации
    proposedEnergy = system.proposedEnergy() # вычисление энергии нового состояния

    if (random() < exp(entropy[currentEnergy]-entropy[proposedEnergy])):
        # если принято, обновляем энергию и систему
        currentEnergy = proposedEnergy
        system.acceptProposedConfiguration()
    else:
        # если отклонено
        system.rejectProposedConfiguration()
    
    H[currentEnergy] += 1
    entropy[currentEnergy] += f
    
    if (isFlat(H)): # isFlat проверяет достаточно ли гладкая гистограмма (например, 95%)
        H[:] = 0
        f *= 0.5 # refine the f parameter