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

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

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

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

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

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

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

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

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

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

${\displaystyle 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)}$.

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

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

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

Это важный шаг алгоритма, и это то, что делает алгоритм Ванга-Ландау немарковским: случайный процесс теперь зависит от истории процесса. Таким образом, когда в следующий раз будет предложено состояние с энергией ${\displaystyle E_{i}}$, это состояние будет отклонено с большей вероятностью; в этом смысле, алгоритм принуждает систему посещать все состояния с одинаковой частотой. Как следствие, гистограмма ${\displaystyle H(E)}$ становится все более и более плоской. Хотя, эта равномерность зависит от того, насколько посчитанная энтропия близка к точной энтропии, что зависит от ${\displaystyle f}$. Для улучшения приближения точной энтропии (и, таким образом, равномерности гистограммы), ${\displaystyle f}$ уменьшается после ${\displaystyle M}$ шагов генерации-принятия:

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

Через некоторое время было показано, что при изменении ${\displaystyle f}$постоянным делением на два алгоритм может не сходиться[5]. Небольшая модификация метода Ванга-Ландау позволяет избежать этого: производится деление не на два, а на ${\displaystyle t}$, при чем ${\displaystyle t}$ пропорционально шагу моделирования.

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

Ниже показан пример кода на Python, в котором предполагается симметричность функции распределения ${\displaystyle g}$:

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

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