Algoritmo de Wang-Landau

Este artículo o sección sobre matemáticas necesita ser wikificado con un formato acorde a las convenciones de estilo. Por favor, edítalo para que las cumpla. Mientras tanto, no elimines este aviso puesto el 19 de diciembre de 2010. También puedes ayudar wikificando otros artículos.

El algoritmo Wang-Landau es un Método Monte Carlo que permite calcular la densidad de estados (density of states en inglés: 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", osea 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.

El algoritmo Wang-Landau fue propuesto por Fugao Wang y David P. Landau ^[1] es una extensión del muestreo Metropolis Monte Carlo.

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 mas alta energía es muy grande el algoritmo de Metropolis tendera a quedarse atrapado en algunos de los mínimos locales del sistema ya que el factor de Boltzmann hara que la probabilidad de escapar sea muy baja. A temperaturas altas el algoritmo de Metropolis efectuara 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 para estudiar cristales de espín pero su uso se extendido a otras areas.

En principio el algoritmo Wang-Landau se aplica para 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 version inicial del algoritmo Wang-Landau se deberían conocer los estados de maxima y mínima energía del sistema para establecer de que tamano sera el arreglo donde sera 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 tamano 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 numero necesario de bins o elementos de nuestro arreglo seria,

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 numero de bins.

Para comenzar la simulación inicializamos el arreglo de la DOS g(E) a cero. Se emplea un arreglo auxiliar guarda un record de cuando un bin de energía ha sido visitado anadiendo 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

donde p es un numero 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

usualmente se elige f = e = 2.72 (Euler's number). 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,

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,

Notese 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,

Para entender por que el algoritmo Wang-Landau funciona miremos la primeramente distribución de probabilidad en el algortimo metrópolis que es dada por

P(E) = P_Metropolisg(E) = e ^{− βE}g(E),

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

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 numero 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 multicanonicas, metadinamica y simulaciones de temperaturas en paralelo.

Sistema de prueba

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

E(x) = x²,

la densidad de estados de este sistema puede ser calculada analíticamente y es dada por,

el calculo de esta integral arroja

g(E)∼E ^{− 1 / 2},

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

Ejemplo de un código

Tomemos como un ejemplo el siguiente código:

real[real] log_g;     // this associative array stores the log of the density of states
uint[real] H;         // this stores the histogram
 
real E_old = the_system.energy();    // this stores the old energy.
real F = 1;                          // the modification factor.
 
 
while(F > epsilon) {
  the_system.evolve();               // propose a new configuration of the system.
  real E_new = the_system.energy();  // compute the current system energy.
 
  if (log_g[E_new] < log_g[E_old])   // check if the proposal is accepted.
    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]++;                 // update the histogram and density of states.
  log_g[E_old] += F;
 
  if (is_flat(H)) {           // when the histogram is flat,
    H[] = 0;                  //   reset it and reduce the modification factor.
    F *= 0.5;
  }
}
 
 
return log_g;

Referencias