Algoritmo de Wang y Landau

El algoritmo de Wang y Landau , propuesto por Fugao Wang y David P. Landau , ^{[ 1 ] [ 2 ] [ 3 ]} es un método de Monte Carlo diseñado para estimar la densidad de estados de un sistema. El método realiza un paseo aleatorio no markoviano para construir la densidad de estados recorriendo rápidamente todo el espectro de energía disponible. El algoritmo de Wang y Landau es un método importante para obtener la densidad de estados necesaria para realizar una simulación multicanónica .

El algoritmo de Wang-Landau puede aplicarse a cualquier sistema caracterizado por una función de coste (o energía). Por ejemplo, se ha aplicado a la solución de integrales numéricas ^{[ 4 ]} y al plegamiento de proteínas. ^{[ 5 ] [ 6 ]} El muestreo de Wang-Landau está relacionado con el algoritmo de metadinámica . ^{[ 7 ]}

El algoritmo de Wang y Landau se utiliza para obtener una estimación de la densidad de estados de un sistema caracterizado por una función de coste. Emplea un proceso estocástico no markoviano que converge asintóticamente a un conjunto multicanónico ^{[ 1 ]} (es decir, a un algoritmo de Metropolis-Hastings con una distribución de muestreo inversa a la densidad de estados). La principal consecuencia es que esta distribución de muestreo conduce a una simulación donde las barreras de energía son invisibles. Esto significa que el algoritmo visita todos los estados accesibles (favorables y menos favorables) mucho más rápido que un algoritmo de Metropolis ^{[ 8 ] .}

Consideremos un sistema definido en un espacio de fases.${\displaystyle \Omega }$y una función de coste,${\displaystyle E}$(la energía), limitada en un espectro${\displaystyle E\in \Gamma =[E_{\min },E_{\max }]}$.${\displaystyle E}$tiene una densidad de estados asociada,${\displaystyle \rho (E)}$, que debe ser estimado. El estimador es${\displaystyle {\hat {\rho }}(E)\equiv \exp(S(E))}$. Debido a que el algoritmo de Wang y Landau funciona en espectros discretos, ^{[ 1 ]} el espectro${\displaystyle \Gamma }$está dividido en${\displaystyle N}$valores discretos con una diferencia entre ellos de${\displaystyle \Delta }$, de tal manera que

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

Dado este espectro discreto, el algoritmo se inicializa de la siguiente manera:

• estableciendo a cero todas las entradas de la entropía microcanónica,${\displaystyle S(E_{i})=0\ \ i=1,2,...,N}$
• inicializando${\displaystyle f=1}$y
• inicializando el sistema aleatoriamente, introduciendo una configuración aleatoria.${\displaystyle {\boldsymbol {r}}\in \Omega }$.

El algoritmo realiza entonces una simulación de conjunto multicanónico : ^{[ 1 ]} un paseo aleatorio de Metropolis-Hastings en el espacio de fases del sistema con una distribución de probabilidad dada por${\displaystyle P({\boldsymbol {r}})=1/{\hat {\rho }}(E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))}$y una probabilidad de proponer un nuevo estado dada por una distribución de probabilidad${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$Un histograma${\displaystyle H(E)}$Se almacena la información sobre las energías visitadas. Al igual que en el algoritmo de Metropolis-Hastings, se realiza un paso de aceptación de propuestas, que consiste en (véase la descripción general del algoritmo de Metropolis-Hastings ):

1. proponiendo un estado${\displaystyle {\boldsymbol {r}}'\in \Omega }$según la distribución de propuestas arbitrarias${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$
2. aceptar/rechazar el estado propuesto según

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

dónde${\displaystyle S=S(E({\boldsymbol {r}}))}$y${\displaystyle S'=S(E({\boldsymbol {r}}'))}$.

Después de cada paso de aceptación de la propuesta, el sistema transita a algún valor.${\displaystyle E_{i}}$,${\displaystyle H(E_{i})}$se incrementa en uno y se realiza la siguiente actualización:

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

Este es el paso crucial del algoritmo, y es lo que hace que el algoritmo de Wang y Landau no sea markoviano: el proceso estocástico ahora depende de la historia del proceso. Por lo tanto, la próxima vez que haya una propuesta a un estado con esa energía en particular${\displaystyle E_{i}}$, ahora es más probable que se rechace esa propuesta; en este sentido, el algoritmo obliga al sistema a visitar todo el espectro por igual. ^{[ 1 ]} La consecuencia es que el histograma${\displaystyle H(E)}$es cada vez más plano. Sin embargo, esta planitud depende de cuán bien se aproxime la entropía calculada a la entropía exacta, que naturalmente depende del valor de f. ^{[ 9 ]} Para aproximar cada vez mejor la entropía exacta (y por lo tanto la planitud del histograma ), f se disminuye después de M pasos de propuesta-aceptación:

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

Posteriormente se demostró que actualizar f dividiendo constantemente por dos puede conducir a errores de saturación. ^{[ 9 ]} Una pequeña modificación al método de Wang y Landau para evitar este problema es utilizar el factor f proporcional a${\displaystyle 1/t}$, dónde${\displaystyle t}$es proporcional al número de pasos de la simulación. ^{[ 9 ]}

Queremos obtener la densidad de estados (DOS) para el potencial del oscilador armónico .

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

La DOS analítica viene dada por,

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

Al realizar la última integral obtenemos

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

En general, la densidad de estados (DOS) para un oscilador armónico multidimensional vendrá dada por alguna potencia de E , cuyo exponente será una función de la dimensión del sistema.

Por lo tanto, podemos usar un potencial de oscilador armónico simple para probar la precisión del algoritmo de Wang-Landau porque ya conocemos la forma analítica de la densidad de estados. Por consiguiente, comparamos la densidad de estados estimada.${\displaystyle {\hat {\rho }}(E)}$obtenido mediante el algoritmo de Wang-Landau con${\displaystyle \rho (E)}$.

A continuación se muestra un código de ejemplo del algoritmo de Wang-Landau en Python , donde asumimos que se utiliza una distribución de propuesta simétrica g:

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

El código considera un "sistema", que es el sistema subyacente que se está estudiando.

energía_actual = sistema.configuración_aleatoria ( ) # Una configuración inicial aleatoriaMientras f > épsilon : system.propos_configuration () # Se propone una configuración propuesta proposed_energy = system.proposed_energy ( ) # Se calcula la energía de la configuración propuestaif random () < exp ( entropy [ current_energy ] - entropy [ proposed_energy ]): # Si se acepta, actualiza la energía y el sistema: current_energy = proposed_energy system.accept_proposed_configuration () else : # Si se rechaza system.reject_proposed_configuration ( )H [ energía_actual ] += 1 entropía [ energía_actual ] += fif is_flat ( H ): # is_flat comprueba si el histograma es plano (por ejemplo, 95% de planitud) H [:] = 0 f *= 0.5 # Refinar el parámetro f

El algoritmo de Wang-Landau se presta a la paralelización (utilizando múltiples núcleos de CPU y/o GPU para mejorar el rendimiento de la simulación) de tres maneras clave. Primero, el dominio de energía puede dividirse en varios subdominios de energía más pequeños, con un caminante de Wang-Landau independiente que recopila estadísticas en cada subdominio. ^{[ 2 ]} (En este caso, la DOS global se recupera periódicamente a partir de la DOS local). Estos subdominios de energía pueden elegirse para que sean no uniformes en tamaño, para asegurar una mayor concentración de caminantes en regiones donde la densidad de estados es más difícil de converger. ^{[ 10 ]} Segundo, en un dominio de energía dado, múltiples caminantes de Wang-Landau independientes pueden recopilar estadísticas que luego se combinan. ^{[ 2 ] [ 11 ]} Tercero, en el caso de que se utilicen múltiples subdominios de energía, una superposición entre subdominios de energía facilita el intercambio de réplicas , es decir , el intercambio de caminantes de un subdominio de energía a otro, lo que puede mejorar la exploración del espacio de configuración. ^{[ 12 ] [ 13 ] [ 14 ]}

La dinámica molecular (DM) suele ser preferible a la de Monte Carlo (MC), por lo que es deseable contar con un algoritmo de DM que incorpore la idea básica de WL para el muestreo de energía plana. Dicho algoritmo es la Dinámica Molecular de Temperatura Estadística (STMD), desarrollada ^{[ 15 ]} por Jaegil Kim et al. en la Universidad de Boston.

Se dio un primer paso esencial con el algoritmo de Monte Carlo de temperatura estadística (STMC). WLMC requiere un aumento extenso en el número de intervalos de energía con el tamaño del sistema, causado por trabajar directamente con la densidad de estados. STMC se centra en una cantidad intensiva, la temperatura estadística,${\displaystyle T(E)=1/(dS(E)/dE)}$, donde E es la energía potencial. Cuando se combina con la relación,${\displaystyle \Omega (E)=e^{S(E)}}$, donde nos instalamos${\displaystyle k_{B}=1}$, la regla WL para actualizar la densidad de estados proporciona la regla para actualizar la temperatura estadística discretizada,

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

dónde${\displaystyle \alpha _{j\pm 1}=1/(1\mp \delta f{\tilde {T}}_{j\pm 1}),\delta f=(\ln f/2\Delta E),\Delta E}$es el tamaño del contenedor de energía, y${\displaystyle {\tilde {T}}}$denota la estimación móvil. Definimos f como en, ^{[ 1 ]} un factor >1 que multiplica la estimación de la DOS para el i-ésimo intervalo de energía cuando el sistema visita una energía en ese intervalo.

Los detalles se dan en la Ref. ^{[ 15 ]} Con una suposición inicial para${\displaystyle T(E)}$y el rango restringido a estar entre${\displaystyle T_{L}}$y${\displaystyle T_{U}}$, la simulación procede como en WLMC, con diferencias numéricas significativas. Una interpolación de${\displaystyle {\tilde {T}}(E)}$proporciona una expresión continua de la estimación${\displaystyle S(E)}$al integrar su inversa, lo que permite el uso de intervalos de energía más grandes que en WL. Diferentes valores de${\displaystyle S(E)}$están disponibles dentro del mismo intervalo de energía al evaluar la probabilidad de aceptación. Cuando las fluctuaciones del histograma son menores al 20% de la media,${\displaystyle f}$ se reduce según${\displaystyle f\rightarrow {\sqrt {f}}}$.

STMC se comparó con WL para el modelo de Ising y el líquido de Lennard-Jones. Al aumentar el tamaño del intervalo de energía, STMC obtiene los mismos resultados en un rango considerable, mientras que el rendimiento de WL se deteriora rápidamente. STMD puede utilizar valores iniciales más pequeños de${\displaystyle f_{d}=f-1}$para una convergencia más rápida. En resumen, STMC necesita menos pasos para obtener la misma calidad de resultados.

Ahora consideremos el resultado principal, STMD. Se basa en la observación de que en una simulación MD estándar a temperatura${\displaystyle T_{0}}$con fuerzas derivadas de la energía potencial${\displaystyle E([x])}$, dónde${\displaystyle [x]}$denota todas las posiciones, el peso de muestreo para una configuración es${\displaystyle e^{-E([x])/T_{0}}}$Además, si las fuerzas se derivan de una función${\displaystyle W(E)}$, el peso del muestreo es${\displaystyle e^{-W(E([x]))/T_{0}}}$.

Para el muestreo de energía plana, sea el potencial efectivo${\displaystyle T_{0}S(E)}$- dinámica molecular entrópica. Entonces el peso es${\displaystyle e^{-S(E)}}$. Dado que la densidad de estados es${\displaystyle e^{+S(E)}}$Su producto ofrece un muestreo de energía plano.

Las fuerzas se calculan como

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

dónde${\displaystyle F^{0}}$denota la fuerza usual derivada de la energía potencial. Escalando las fuerzas usuales por el factor${\displaystyle (T_{0}/T(E))}$produce un muestreo de energía plano.

STMD comienza con un algoritmo MD ordinario a velocidad constante.${\displaystyle T_{0}}$y V. Las fuerzas se escalan como se indica, y la temperatura estadística se actualiza en cada paso de tiempo, utilizando el mismo procedimiento que en STMC. A medida que la simulación converge al muestreo de energía plano, la estimación en curso${\displaystyle {\tilde {T}}(E)}$converge a la verdadera${\displaystyle T(E)}$Los detalles técnicos, incluidos los pasos para acelerar la convergencia, se describen en ^{[ 15 ]} y ^{[ 16 ] .}

En STMD${\displaystyle T_{0}}$se denomina temperatura cinética ya que controla las velocidades como de costumbre, pero no entra en el muestreo configuracional, lo cual es inusual. Por lo tanto, STMD puede sondear bajas energías con partículas rápidas. Cualquier promedio canónico se puede calcular con reponderación, pero la temperatura estadística,${\displaystyle T(E)}$, está disponible de inmediato sin análisis adicional. Es extremadamente valioso para estudiar transiciones de fase. En nanosistemas finitos${\displaystyle T(E)}$tiene una característica correspondiente a cada "transición de subfase". Para una transición suficientemente fuerte, una construcción de área igual en un bucle S en${\displaystyle 1/T(E)}$da la temperatura de transición.

STMD ha sido perfeccionado por el grupo BU, ^{[ 16 ]} y aplicado a varios sistemas por ellos y otros. D. Stelter reconoció que, a pesar de nuestro énfasis en trabajar con cantidades intensivas,${\displaystyle \ln(f)}$es extenso. Sin embargo${\displaystyle \delta f=(\ln(f)/2\Delta E)}$es intensivo y el procedimiento${\displaystyle f\rightarrow {\sqrt {f}}}$basado en la planitud del histograma se reemplaza por corte${\displaystyle \delta f}$a la mitad cada número fijo de pasos de tiempo. Este simple cambio hace que STMD sea totalmente intensivo y mejora sustancialmente el rendimiento para sistemas grandes. ^{[ 16 ]} Además, el valor final del intensivo${\displaystyle \delta f}$es una constante que determina la magnitud del error en la convergencia${\displaystyle T(E)}$y es independiente del tamaño del sistema. STMD se implementa en LAMMPS como fix stmd.

La simulación STMD es particularmente útil para las transiciones de fase. Es imposible obtener información de equilibrio con una simulación canónica, ya que se requiere sobreenfriamiento o sobrecalentamiento para que se produzca la transición. Sin embargo, una simulación STMD obtiene un muestreo de energía uniforme con una progresión natural de calentamiento y enfriamiento, sin quedar atrapado en el estado de baja o alta energía. Más recientemente, se ha aplicado a la transición fluido/gel ^{[ 16 ]} en nanopartículas recubiertas de lípidos.

El intercambio de réplicas STMD ^{[ 17 ]} también ha sido presentado por el grupo BU.