Algoritmo de Wang-Landau

El algoritmo Wang-Landau es una extensión del método de Monte Carlo, propuesto por Fugao Wang y David P. Landau, que permite calcular la densidad de estados (en inglés density of states o DOS) de un sistema dado sin tener ningún conocimiento previo sobre ella. Se dice entonces que la densidad de estados se calcula on the fly, o sea durante la simulación. El algoritmo es independiente de la temperatura (un problema común con otros algoritmos Monte Carlo) y es fácil de implementar.^[1]

Este algoritmo fue pensado inicialmente para muestrear el espacio de configuraciones de ciertos sistemas que no podían ser tratados mediante el algoritmo de Metropolis. El algoritmo de Metropolis realiza un muestreo usando el factor de Boltzmann para aceptar o descartar configuraciones. Metropolis funciona bien si todas las configuraciones posibles del sistema se encuentran dentro de un rango relativamente angosto de energías. Si por el contrario la separación entre las configuraciones de menor energía y las de más alta energía es muy grande el algoritmo de Metropolis tenderá a quedarse atrapado en algunos de los mínimos locales del sistema ya que el factor de Boltzmann hará que la probabilidad de escapar sea muy baja. A temperaturas altas el algoritmo de Metropolis efectuará un muestreo aleatorio de las configuraciones de alta energía, las de baja energía no serán muestreadas adecuadamente. A temperaturas bajas Metropolis explora las configuraciones de baja energía pero existe el riesgo de quedar atrapado en un mínimo local.

Para evitar los problemas anteriores Wang y Landau propusieron el algoritmo que lleva su nombre y que nos permite realizar un muestreo uniforme del espacio de configuraciones de un sistema dado. Fue aplicado inicialmente al estudio de los cristales de espín, pero su uso se extendió a otras áreas.

En principio, el algoritmo Wang-Landau, se aplica a cualquier sistema que sea descrito por una función de costo o energía. Ha sido aplicado a la solución de integrales^[2] y al plegamiento de proteínas.^[3]

Algoritmo

De acuerdo a la versión inicial del algoritmo Wang-Landau se deberían conocer los estados de máxima y mínima energía del sistema para establecer de qué tamaño será el arreglo donde será almacenada la densidad de estados (DOS). Sin embargo este requisito no es indispensable si en nuestro programa empleamos arreglos dinámicos de memoria ya que en este caso podemos modificar el tamaño del arreglo durante la simulación computacional.

Cabe aclarar que la energía no es el único parámetro que podemos emplear en este algoritmo. Podemos utilizar otra propiedad, por ejemplo en proteínas se utiliza la distancia extremo-extremo o bien el radio de giro.

Supongamos que tenemos los valores máximos y mínimos de la energía y deseamos una precisión Δ. El número necesario de bins o elementos de nuestro arreglo sería:

N={\frac {E_{\mathrm {max} }-E_{\mathrm {min} }}{\Delta }},

Se debe tener cuidado a la hora de elegir la precisión de nuestra simulación ya que el tiempo de computo aumenta rápidamente cuando incrementamos el número de bins.

Para comenzar la simulación inicializamos el arreglo de la DOS g(E) a cero. Se emplea un arreglo auxiliar que guarda un registro de cuando un bin de energía ha sido visitado añadiendo 1 a esta entrada. Este arreglo H(E) también es inicializado, al principio de la simulación, a cero.

Se empieza con una cierta configuración del sistema y entonces se propone un movimiento de Monte Carlo, el movimiento se acepta si se cumple la ecuación:

p<\min \left\{1,{\frac {g(E)}{g(E')}}\right\},

donde p es un número uniforme en el intervalo [0, 1), E es la energía de la configuración actual y E′ la energía de la configuración propuesta.

Si el movimiento es rechazado, el histograma H(E) es incrementado en la energía actual y la DOS es multiplicada por el factor constante f:

g(E)\rightarrow g(E)\times f,

usualmente se elige f = e = 2,72 (número de Euler). Si la movida es aceptada se actualizan igualmente ambos arreglos pero en el bin de energía propuesto. Cuando el histograma de visitas H(E) es suficientemente plano lo reinicializamos a cero y el factor f se actualiza de la siguiente manera:

f_{i+1}\rightarrow {\sqrt {f_{i}}},

Si tratamos de repetir este esquema pronto notaremos que el programa marca un error. El problema radica en que al multiplicar continuamente cada entrada de la densidad de estados por f obtendremos números extremadamente grandes que no podrán ser almacenados por ningún ordenador. Por esta razón se elige generalmente trabajar en una escala logarítmica reemplazando multiplicaciones con sumas, divisiones con restas, etc. de la siguiente manera:

\log[g(E)]\rightarrow \log[g(E)]+\log[f],

Nótese que la elección f = e es fortuita ya que log(e)=1 en la base natural de los logaritmos.

El factor de modificación se actualiza como:

\log(f_{i+1})\rightarrow (1/2)\log(f_{i}),

Para entender por qué el algoritmo Wang-Landau funciona miremos primeramente la distribución de probabilidad en el algoritmo metrópolis que está dada por:

P(E)=P_{\text{Metropolis}}g(E)=e^{-\beta E}g(E),

en el algoritmo de Wang-Landau la distribución de probabilidad es:

P(E)=P_{\text{WangLandau}}g(E)={\frac {1}{g(E)}}g(E)=1,

de modo que esta probabilidad es constante para todo E.

Lo anterior nos indica que con el algoritmo Wang-Landau realizamos un muestreo uniforme en el espacio de configuraciones del sistema bajo estudio.

Observaciones

La elección del factor de modificación clásica (tomando la raíz cuadrada del factor en la simulación anterior) puede conducirnos al así llamado «error de saturación».^[4] El cual consiste en que después de un cierto número de actualizaciones del factor de modificación el error de la simulación alcanza un valor constante y por tanto la precisión de nuestros cálculos no puede ser ya mejorada. Para evitar este problema se propuso el algoritmo $1/t$ , donde $t$ es el tiempo de simulación. Este algoritmo es una variante del algoritmo de Wang-Landau clásico.

Otro punto importante respecto al algoritmo Wang-Landau es que es completamente independiente de la temperatura a diferencia de las simulaciones multicanónicas, metadinámica y simulaciones de temperaturas en paralelo.

Sistema de prueba

Tomemos el caso de una partícula en un potencial armónico en una dimensión (1D):

E(x)=x^{2},

la densidad de estados de este sistema puede ser calculada analíticamente y está dada por:

g(E)=\int \delta (E(x)-E_{0})dx=\int \delta (x^{2}-E_{0})dx,

el cálculo de esta integral arroja:

g(E)\sim E^{-1/2},

en general, la DOS para un oscilador armónico multidimensional será alguna potencia de E cuyo exponente dependerá de la dimensión del sistema.

Ejemplo de un código

Tomemos como un ejemplo el siguiente código:

real[real] log_g;     // esta matriz asociativa almacena el logaritmo de la densidad de estados
uint[real] H;         // esto almacena el histograma

real E_old = the_system.energy();    // esto almacena la energía anterior.
real F = 1;                          // el factor de modificación.

while(F > epsilon) {
  the_system.evolve();               // propone una nueva configuración al sistema.
  real E_new = the_system.energy();  // computa la energía del sistema actual.

  if (log_g[E_new] < log_g[E_old])   // comprueba si el supuesto es aceptado.
    E_old = E_new;
  else if (random() < exp(log_g[E_old] - log_g[E_new]))
    E_old = E_new;
  else
    the_system.reject_and_undo();

  H[E_old]++;                 // actualiza el histograma y la densidad de estados.
  log_g[E_old] += F;

  if (is_flat(H)) {           // cuando el histograma es plano,,
    H[] = 0;                  //   reinícialo y reduce el factor de modificación.
    F *= 0.5;
  }
}

return log_g;

Referencias

↑ Wang, Fugao and Landau, D. P. (marzo de 2001). «Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States». Phys. Rev. Lett. (en inglés) (American Physical Society) 86 (10): 2050-2053. doi:10.1103/PhysRevLett.86.2050.
↑ R. E. Belardinelli and S. Manzi and V. D. Pereyra (diciembre de 2008). «Analysis of the convergence of the 1∕t and Wang-Landau algorithms in the calculation of multidimensional integrals». Phys. Rev. E (en inglés) (American Physical Society) 78: 067701. doi:10.1103/PhysRevE.78.067701.
↑ P. Ojeda and M. Garcia and A. Londono and N.Y. Chen (febrero de 2009). «Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States». Biophys. Jour. (en inglés) (Biophysical Society) 96 (3): 1076-1082. doi:10.1529/biophysj.107.125369.
↑ Belardinelli, R. E. and Pereyra, V. D. (noviembre de 2007). «Wang-Landau algorithm: A theoretical analysis of the saturation of the error». Jour. Chem. Phys. (en inglés) 127 (18): 184105. doi:10.1063/1.2803061.

Datos: Q6522961

[1] Wang, Fugao and Landau, D. P. (marzo de 2001). «Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States». Phys. Rev. Lett. (en inglés) (American Physical Society) 86 (10): 2050-2053. doi:10.1103/PhysRevLett.86.2050.

[2] R. E. Belardinelli and S. Manzi and V. D. Pereyra (diciembre de 2008). «Analysis of the convergence of the 1∕t and Wang-Landau algorithms in the calculation of multidimensional integrals». Phys. Rev. E (en inglés) (American Physical Society) 78: 067701. doi:10.1103/PhysRevE.78.067701.

[3] P. Ojeda and M. Garcia and A. Londono and N.Y. Chen (febrero de 2009). «Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States». Biophys. Jour. (en inglés) (Biophysical Society) 96 (3): 1076-1082. doi:10.1529/biophysj.107.125369.

[4] Belardinelli, R. E. and Pereyra, V. D. (noviembre de 2007). «Wang-Landau algorithm: A theoretical analysis of the saturation of the error». Jour. Chem. Phys. (en inglés) 127 (18): 184105. doi:10.1063/1.2803061.

[1]

[2]

[3]

[4]