Año 2 · Módulo 16: Física estadística: microestados, equilibrio y transiciones · Transferencia e integración

159.1. El problema computacional de observar un equilibrio

El modelo de Ising permite describir un material mediante variables muy simples. En cada sitio de una red hay un espín que vale +1+ 1 o −1- 1 y las interacciones favorecen o penalizan configuraciones según sus orientaciones. Esa economía de componentes no se traduce en facilidad de cálculo. Una red cuadrada de lado LL contiene N=L2N = L^{2} espines y admite 2N2^{N} configuraciones. Con L=32L = 32, enumerarlas todas significaría recorrer 210242^{1024} estados, una tarea computacionalmente inabordable. Sin embargo, las observaciones físicas que buscamos suelen ser pocas (energía por sitio, magnetización, correlaciones y fluctuaciones).

En el conjunto canónico, a temperatura TT, la probabilidad de una configuración s=(s1,…,sN)s = \left( s_{1},\ldots,s_{N} \right) viene dada por

π(s)=e−βH(s)Z,Z=∑se−βH(s),β=1kBT.\pi(s) = \frac{e^{- \beta H(s)}}{Z},\quad\quad Z = \sum_{s}^{}e^{- \beta H(s)},\quad\quad\beta = \frac{1}{k_{B}T}.

Si quisiéramos calcular exactamente ⟨A⟩=∑sA(s)π(s)\langle A\rangle = \sum_{s}^{}A(s)\pi(s), la suma recorrería todos los microestados. El método de Monte Carlo propone sustituirla por un promedio de configuraciones generadas con frecuencias compatibles con π\pi. La dificultad estriba en producirlas sin conocer ZZ, precisamente la suma que no podemos evaluar. El algoritmo de Metropolis resuelve esa dificultad con cocientes de probabilidades, en los que la función de partición se cancela al comparar dos estados.

La idea va más allá del modelo magnético y aparece también al muestrear conformaciones moleculares, distribuciones posteriores bayesianas y configuraciones de redes con interacciones. Esas aplicaciones no siempre tienen una interpretación térmica, pero conservan el principio de diseñar una cadena de Markov cuya distribución estacionaria sea una distribución objetivo especificada. El nombre histórico procede del trabajo de Metropolis, Arianna y Marshall Rosenbluth, Augusta Teller y Edward Teller, publicado en 1953. La generalización de Hastings, en 1970, permitió trabajar con propuestas asimétricas.

Conviene fijar desde el comienzo una distinción conceptual. Una secuencia de actualizaciones de Metropolis es un procedimiento de muestreo. Sus pasos no corresponden automáticamente a segundos, a choques moleculares ni a la cinética real del material. Podemos obtener buenos promedios de equilibrio con una dinámica artificial. También podemos construir un simulador que parezca físicamente plausible y que, por un error de implementación, no muestree la distribución canónica. El criterio de corrección debe ser explícito y exige que la cadena conserve la distribución buscada, pueda explorar su espacio de estados y produzca estimaciones cuya incertidumbre podamos evaluar.

159.2. La regla de aceptación y su justificación probabilística

Una cadena de Markov cambia de configuración mediante una regla que depende solamente del estado actual y de la propuesta de transición. Supongamos que, estando en ss, elegimos una configuración candidata s′s\prime con probabilidad q(s′∣s)q(s\prime \mid s). Luego decidimos si aceptamos el cambio mediante α(s,s′)\alpha(s,s\prime). Si se rechaza, conservamos ss; esta permanencia, que a veces omiten las descripciones informales, es parte de la transición probabilística del algoritmo.

Para una propuesta simétrica, q(s′∣s)=q(s∣s′)q(s\prime \mid s) = q(s \mid s\prime), Metropolis acepta con

α(s,s′)=min⁡{1,e−β[H(s′)−H(s)]}.\boxed{\alpha(s,s\prime) = \min\{ 1,e^{- \beta\left\lbrack H(s\prime) - H(s) \right\rbrack}\}.}

Cuando el cambio disminuye la energía, ΔH=H(s′)−H(s)≤0\Delta H = H(s\prime) - H(s) \leq 0, se acepta siempre. Si la energía aumenta, se acepta con probabilidad e−βΔHe^{- \beta\Delta H}, positiva a temperatura positiva. Este segundo caso es indispensable, porque un algoritmo que aceptara únicamente descensos de energía buscaría mínimos locales y dejaría de representar un equilibrio térmico, en el que las fluctuaciones sostienen ocasionalmente estados más costosos.

La regla satisface una condición suficiente para que π\pi sea estacionaria, el balance detallado,

π(s)q(s′∣s)α(s,s′)=π(s′)q(s∣s′)α(s′,s).\pi(s)q(s\prime \mid s)\alpha(s,s\prime) = \pi(s\prime)q(s \mid s\prime)\alpha(s\prime,s).

La igualdad se verifica separando dos casos. Si H(s′)≥H(s)H(s\prime) \geq H(s), la transición hacia el estado más energético incorpora e−βΔHe^{- \beta\Delta H}; la transición inversa se acepta siempre. Como π(s′)/π(s)=e−βΔH\pi(s\prime)/\pi(s) = e^{- \beta\Delta H}, los flujos en ambas direcciones coinciden. Si H(s′)<H(s)H(s\prime) < H(s), basta intercambiar sus papeles. La suma de los flujos entrantes reproduce entonces la probabilidad de cada estado.

El balance detallado no es la única manera de construir una cadena correcta. Existe la condición más general de balance global y hay muestreadores no reversibles que usan esa libertad. Tampoco la estacionariedad equivale a convergencia desde cualquier estado inicial. En un espacio finito se necesitan condiciones como irreducibilidad y aperiodicidad para obtener la distribución estacionaria única y converger hacia ella. En el Ising con actualizaciones aleatorias de un único espín, T>0T > 0, campo finito y conectividad habitual, cualquier configuración es alcanzable mediante una secuencia de inversiones aceptadas con probabilidad positiva. Los rechazos producen auto-transiciones que en general hacen aperiódica la cadena; aunque algunas configuraciones acepten siempre su siguiente propuesta, en el caso ordinario hay auto-transiciones en otras configuraciones accesibles que bastan para evitar una periodicidad global.

Cuando las propuestas no son simétricas, el cociente correcto es

α(s,s′)=min⁡{1,π(s′)q(s∣s′)π(s)q(s′∣s)}.\alpha(s,s\prime) = \min\left\{ 1,\frac{\pi(s\prime)q(s \mid s\prime)}{\pi(s)q(s\prime \mid s)} \right\}.

Ignorar los factores qq en una propuesta sesgada altera la distribución estacionaria. Esta es la diferencia sustantiva entre la forma elemental de Metropolis y el esquema de Metropolis–Hastings. La demostración de validez probabilística no garantiza rapidez de exploración; una cadena perfectamente correcta puede necesitar un número enorme de actualizaciones para atravesar una barrera de energía libre.

159.3. Un cambio local en la red de Ising

Consideremos el Hamiltoniano ferromagnético con vecinos próximos,

H(s)=−J∑⟨i,j⟩sisj−h∑isi,J>0,si∈{−1,+1}.H(s) = - J\sum_{\langle i,j\rangle}^{}s_{i}s_{j} - h\sum_{i}^{}s_{i},\quad\quad J > 0,\quad s_{i} \in \{ - 1, + 1\}.

El primer término suma cada enlace físico una sola vez. El segundo acopla todos los espines a un campo externo hh, expresado aquí en unidades de energía. Cambiar el signo de sis_{i} conserva todos los demás espines. Por ello, aunque calcular HH desde cero requiera recorrer toda la red, la diferencia energética de una inversión depende solo del entorno inmediato y vale

ΔHi=2si(J∑j∈∂isj+h).\boxed{\Delta H_{i} = 2s_{i}\left( J\sum_{j \in \partial i}^{}s_{j} + h \right).}

La fórmula se obtiene observando que cada producto −Jsisj- Js_{i}s_{j} cambia de signo y aporta 2Jsisj2Js_{i}s_{j} a la diferencia; ocurre lo mismo con el término de campo. Para una red cuadrada interior o periódica, cada sitio tiene cuatro vecinos, de modo que la suma vecinal puede valer −4,−2,0,2,4- 4, - 2,0,2,4. Esta estructura discreta permite incluso precalcular probabilidades de aceptación cuando JJ, hh y TT permanecen constantes.

Imaginemos una región con todos los espines alineados en +1+ 1, campo nulo y J=1J = 1. Invertir un espín del interior rompe cuatro enlaces favorables y produce ΔH=8\Delta H = 8. A temperatura kBT=3k_{B}T = 3, la aceptación de ese movimiento es exp⁡(−8/3)\exp( - 8/3), aproximadamente 0.069480.06948. Una inversión que recupere un espín desalineado de ese mismo entorno disminuye la energía y se acepta siempre. De esta asimetría surgen dominios persistentes, aunque la cadena conserve la posibilidad de escapar de ellos a temperatura positiva.

La condición de frontera interviene directamente en el cálculo. Con fronteras periódicas, el vecino derecho del último sitio de una fila es el primero. Con bordes abiertos, los sitios periféricos tienen menos vecinos y la energía incluye menos enlaces; no se debe fingir que tienen cuatro vecinos asignándoles espines inexistentes. En redes muy pequeñas, incluso el modo de contabilizar enlaces periódicos necesita cuidado. En una red de lado dos, sumar enlaces «derecha» y «abajo» por sitio puede describir un multigrafo con enlaces repetidos entre las mismas parejas. Para comparar resultados con la convención de una red ordinaria conviene elegir L≥3L \geq 3 o especificar explícitamente ese detalle.

También importa la convención de unidades. Si fijamos J=kB=1J = k_{B} = 1, la temperatura del código es TT medida en unidades J/kBJ/k_{B}. La energía, la magnetización y el campo siguen teniendo normalizaciones distintas. Un error frecuente consiste en multiplicar dos veces por β\beta o usar la fórmula de aceptación correspondiente a h=0h = 0 cuando el Hamiltoniano sí incluye el campo. Ambos errores pueden producir gráficos visualmente razonables, así que cada implementación debería verificarse contra recuentos exactos en sistemas pequeños.

Los movimientos de energía negativa se aceptan siempre. Invertir un espín con cuatro vecinos alineados cuesta 8J y se acepta un 6,95 % a T = 3.
Figura 159.1. Con propuestas simétricas, los descensos se aceptan y los ascensos se aceptan con el factor de Boltzmann. A campo nulo, invertir un espín alineado con sus cuatro vecinos cambia cuatro enlaces de −J a +J: ΔH = 8J. Con T = 3 y J = kB = 1, la aceptación es exp(−8/3). Ampliar figura ↗

159.4. Construir una cadena reproducible sin confundir intentos y barridos

Para actualizar la red se selecciona al azar un sitio, se calcula ΔHi\Delta H_{i}, se genera un número uniforme u∈[0,1)u \in \lbrack 0,1) y se invierte el espín si u<min⁡(1,e−βΔHi)u < \min\left( 1,e^{- \beta\Delta H_{i}} \right). Una vez aceptada la propuesta, puede actualizarse la energía almacenada con H←H+ΔHiH \leftarrow H + \Delta H_{i}; si se rechaza, ni los espines ni la energía cambian. El generador de números pseudoaleatorios necesita un estado inicial reproducible, pero la semilla no sustituye una evaluación de errores estadísticos.

El esquema siguiente describe el algoritmo con suficiente precisión.

inicializar la red s y calcular H(s)
repetir para cada barrido:
    repetir N veces:
        elegir un sitio i uniformemente, con reemplazo
        calcular delta_H del cambio s_i -> -s_i
        si delta_H <= 0 o U(0,1) < exp(-beta * delta_H):
            invertir s_i y actualizar H
    registrar H, M y otros observables en los barridos seleccionados

Llamamos intento a cada propuesta de inversión y barrido a NN intentos. Al elegir sitios con reemplazo, un barrido no garantiza haber visitado todos los sitios una sola vez. Aun así, es una unidad de coste computacional útil para comparar sistemas de distinto tamaño. Otra implementación puede visitar los sitios en orden fijo o en una permutación aleatoria. Una composición de actualizaciones por sitios en orden fijo puede preservar la distribución objetivo si cada subactualización la preserva, aunque la composición completa no tenga por qué satisfacer balance detallado ni reproducir exactamente la cadena de propuestas uniformes aleatorias. Si la validación depende de las probabilidades de transición del algoritmo elemental, estas variantes deben declararse por separado.

Conviene distinguir el coste algorítmico de la escala física. Mil barridos equivalen aproximadamente a mil propuestas por espín, sin correspondencia directa con unidades de tiempo real. Las reglas de aceptación son una dinámica artificial diseñada para el muestreo canónico. Incluso cuando una regla de inversión local se usa como aproximación a una dinámica física de relajación, establecer la correspondencia temporal requiere hipótesis adicionales sobre tasas y mecanismos microscópicos.

Mantener la energía de forma incremental es eficiente, pero basta una diferencia ΔHi\Delta H_{i} mal calculada para que se acumule un error. Una prueba automática consiste en escoger configuraciones aleatorias, calcular HH, invertir un espín, recalcular el Hamiltoniano completo y comprobar que H(s′)−H(s)H(s\prime) - H(s) coincide exactamente con la expresión local. Durante una ejecución prolongada es útil recalcular ocasionalmente HH para verificar que la variable almacenada no se haya desviado. También debe probarse que invertir dos veces el mismo espín restaura el estado y la energía originales.

En la evaluación científica se documentan tamaño, geometría, límites, acoplamientos, temperatura, campo, algoritmo, definición de barrido, generador y semilla. Ese registro vale más que mostrar solo la imagen final de un dominio, que puede ser llamativa y aun así proceder de un cálculo que no haya alcanzado la distribución objetivo.

159.5. Observables, fluctuaciones y consistencia termodinámica

Las configuraciones producidas deben transformarse luego en observables definidos de antemano. Para N=L2N = L^{2}, la magnetización total es M=∑isiM = \sum_{i}^{}s_{i}, la magnetización por espín es m=M/Nm = M/N y la energía por espín es e=H/Ne = H/N. La media muestral de una magnitud AA, medida en nn configuraciones después del período inicial, es ⟨A⟩^=n−1∑t=1nA(st)\widehat{\langle A\rangle} = n^{- 1}\sum_{t = 1}^{n}A\left( s_{t} \right). Si la cadena ha alcanzado el régimen estacionario y es ergódica, esa media converge al valor esperado canónico cuando aumenta la longitud de la muestra.

Algunas identidades independientes del algoritmo permiten hacer pruebas cruzadas. Con campo hh y acoplamiento −hM- hM, la susceptibilidad por espín es

χN=βN(⟨M2⟩−⟨M⟩2).\chi_{N} = \frac{\beta}{N}\left( \langle M^{2}\rangle - \langle M\rangle^{2} \right).

Si el Hamiltoniano no depende explícitamente de la temperatura, la capacidad calorífica a campo fijo satisface

CN=kBβ2(⟨H2⟩−⟨H⟩2).C_{N} = k_{B}\beta^{2}\left( \langle H^{2}\rangle - \langle H\rangle^{2} \right).

Esta CNC_{N} es la capacidad total; para expresarla por espín dividimos por NN. Ambas relaciones requieren distinguir promedios de cuadrados de cuadrados de promedios. Sustituir ⟨M2⟩\langle M^{2}\rangle por ⟨∣M∣⟩2\langle|M|\rangle^{2} en la fórmula de susceptibilidad cambiaría el observable termodinámico, aunque en la práctica se use la magnetización absoluta para localizar regiones ordenadas de una simulación finita.

En campo nulo, la simetría global si↦−sis_{i} \mapsto - s_{i} implica ⟨m⟩=0\langle m\rangle = 0 para cualquier red finita correctamente equilibrada. A temperaturas bajas puede costar observar cambios de signo porque las configuraciones de magnetización positiva y negativa están separadas por barreras. Una trayectoria que permanezca todo el tiempo en una orientación producirá ⟨m⟩^≠0\widehat{\langle m\rangle} \neq 0 y esa cifra puede confundirse con la magnetización espontánea. En cambio, ⟨∣m∣⟩\langle|m|\rangle, ⟨m2⟩\langle m^{2}\rangle y el histograma PN(m)P_{N}(m) permiten observar orden finito sin exigir que las inversiones globales ocurran frecuentemente. No deben usarse, sin embargo, para ocultar un mezclado incompleto; hay que informar si se exploraron los dos sectores de magnetización.

La distribución de magnetización contiene información que una media pierde. Un histograma de dos picos simétricos puede producir promedio cero, aunque casi todas las configuraciones tengan magnetización de módulo grande. La razón de momentos de cuarto orden que define el cumulante de Binder, U4=1−⟨m4⟩/(3⟨m2⟩2)U_{4} = 1 - \langle m^{4}\rangle/\left( 3\langle m^{2}\rangle^{2} \right), da una medida adicional de la forma de esa distribución. Sus cruces al variar LL pueden ser útiles cerca de la temperatura crítica, bajo controles de tamaño, frontera e incertidumbre. Ninguno de estos indicadores remedia un muestreo insuficiente.

159.6. Equilibrar una cadena y reconocer cuándo todavía no ha llegado

La distribución canónica es estacionaria para la transición diseñada, pero una cadena iniciada en un estado particular no comienza necesariamente en equilibrio. Si todos los espines están alineados y la temperatura es alta, la magnetización tardará cierto tiempo en aproximarse a su régimen fluctuante. Si se parte de una red aleatoria a baja temperatura, aparecerán y crecerán dominios antes de que los observables se estabilicen. Descartar una fracción inicial de mediciones es razonable; elegir esa fracción por una regla universal de «tantos barridos para cualquier modelo» no lo es.

El período de adaptación o burn-in debe fijarse según el tamaño, la temperatura, la geometría y la lentitud de los modos colectivos. Un diagnóstico elemental compara varias cadenas iniciadas en estados contrastantes (todos los espines positivos, todos negativos y una configuración aleatoria). Si después de un intervalo suficiente coinciden las distribuciones de energía y magnetización, aumenta la confianza en la equilibración. La coincidencia de medias no asegura convergencia cuando varias cadenas han quedado atrapadas en estados metastables equivalentes, o cuando los eventos lentos relevantes casi nunca aparecen. Por eso también se inspeccionan histogramas, series temporales y cambios de signo cuando estos son posibles.

Una curva plana de energía no prueba equilibrio. A baja temperatura, tanto una fase positiva como una negativa pueden tener prácticamente la misma energía, mientras una cadena permanece confinada durante toda la simulación. Si el campo favorece una de ellas, una trayectoria atrapada en la otra puede mostrar una energía aparentemente estable y un sesgo considerable en observables sensibles al campo. Un promedio constante entre bloques consecutivos puede deberse simplemente a que todos los bloques contienen el mismo estado retenido.

También debe distinguirse entre equilibrar después de un cambio de temperatura y reutilizar configuraciones anteriores. Un barrido térmico en el que cada temperatura comienza desde la última configuración de la temperatura precedente puede ahorrar tiempo, pero induce dependencia entre etapas. Cerca de una transición, esa estrategia puede generar una rama aparente de histéresis si la velocidad de cambio supera la capacidad de relajación. Para estimar cantidades de equilibrio en cada temperatura se necesita demostrar equilibración independiente o usar métodos diseñados para explorar correctamente la familia de distribuciones.

En la práctica se prueban longitudes crecientes de burn-in, se comparan semillas independientes y se examina cuánto cambian las estimaciones al ampliar el intervalo descartado. Las decisiones deben registrarse antes de seleccionar la versión del resultado más atractiva. Cuando la cadena no mezcla entre sectores, corresponde informar esa limitación y modificar el método; no tiene sentido prolongar indefinidamente una ejecución a la espera de que el problema desaparezca por azar. Para redes muy pequeñas, la enumeración exacta permite distinguir un error de inicialización de un error de implementación.

159.7. Muestras dependientes y precisión de un promedio de Monte Carlo

Aunque cada intento de Metropolis use números aleatorios independientes, las configuraciones sucesivas son dependientes. Si el cambio propuesto solo afecta a un espín, dos configuraciones vecinas difieren, como máximo, en una de sus NN variables. La magnetización, la energía y las correlaciones espaciales retienen memoria durante muchos barridos. Tratar cien mil observaciones consecutivas como cien mil muestras independientes produce barras de error demasiado estrechas, a veces por factores grandes cerca de una transición.

Para un observable estacionario AtA_{t}, definamos la autocorrelación normalizada ρA(k)=Cov(At,At+k)/Var(At)\rho_{A}(k) = Cov\left( A_{t},A_{t + k} \right)/Var\left( A_{t} \right). La escala integrada, con la convención habitual, es

τint,A=12+∑k=1∞ρA(k).\tau_{int,A} = \frac{1}{2} + \sum_{k = 1}^{\infty}\rho_{A}(k).

Cuando la suma converge y se cumplen las condiciones de un teorema central del límite para cadenas de Markov, la varianza aproximada del promedio de nn mediciones consecutivas es

Var(Aˉn)≃2τint,AnVarπ(A).Var\left( {\bar{A}}_{n} \right) \simeq \frac{2\tau_{int,A}}{n}{Var}_{\pi}(A).

En ese régimen podemos interpretar nef≃n/(2τint,A)n_{ef} \simeq n/\left( 2\tau_{int,A} \right) como un tamaño efectivo de muestra. Si, por ejemplo, τint,A=25\tau_{int,A} = 25 barridos y observamos 50 00050\, 000 valores, disponemos de una precisión comparable, aproximadamente, a la de mil observaciones independientes con igual varianza. La cifra es ilustrativa; distintos observables tienen distintas autocorrelaciones y una estimación deficiente de la cola de ρ\rho puede subestimar seriamente la incertidumbre.

Para intervalos de confianza, una estrategia operacional consiste en agrupar mediciones en bloques de longitud mucho mayor que la escala de correlación y estudiar la estabilidad de la varianza entre bloques. También sirven los estimadores espectrales, el remuestreo por bloques y las réplicas independientes de cadenas. El bootstrap ordinario de observaciones individuales destruye la estructura de dependencia y, sin corrección, no estima bien la precisión. Aumentar el número de puntos del gráfico no mejora la precisión si no aumenta también el tiempo total de exploración efectivo.

La incertidumbre estadística se distingue del sesgo por falta de equilibración. Un intervalo calculado sobre una cadena atrapada puede ser estrechísimo, porque dentro del pozo metastable las fluctuaciones son pequeñas, aunque la media esté lejos del equilibrio global. Ninguna fórmula de autocorrelación aplicada únicamente a los datos retenidos descubre estados que nunca se visitaron. Por eso la inferencia requiere tanto un análisis interno de dependencia como controles externos, entre ellos inicios contrastantes, ejecuciones más largas y la comparación con resultados analíticos cuando existen.

El término ralentización crítica designa el crecimiento de escalas de relajación al acercarse a ciertos puntos críticos. Para las actualizaciones locales, regiones correlacionadas cada vez mayores tienden a reorganizarse lentamente. En ese caso el número efectivo de configuraciones puede crecer mucho menos que el número de barridos. Registrar datos con más frecuencia no lo remedia, porque las mediciones adicionales apenas añaden información; conviene cambiar la estrategia de actualización o invertir más coste en trayectorias realmente exploratorias.

La elección de la ventana usada al sumar autocorrelaciones también modifica el error comunicado. Wolff (2004; errata de 2007) mostró que truncar demasiado pronto la suma puede subestimar la varianza, mientras que incluir retardos muy largos añade ruido al estimador. La estimación debe examinar la estabilidad del resultado frente a la ventana elegida.

159.8. Caso calculado con nueve espines y una comprobación exacta

Una red de 3×33 \times 3 sitios con fronteras periódicas contiene solo 29=5122^{9} = 512 configuraciones. Es lo bastante pequeña para sumar el conjunto canónico sin aproximaciones de Monte Carlo y suficientemente grande para representar una red cuadrada donde cada sitio tiene cuatro vecinos distintos. Elegimos J=kB=1J = k_{B} = 1, T=3T = 3 y h=0.2h = 0.2. El campo positivo rompe la simetría entre magnetizaciones positivas y negativas. Para cada configuración se calcula la energía contando un enlace hacia la derecha y otro hacia abajo por sitio; así hay exactamente 18 enlaces en esta red.

El cálculo exhaustivo, con pesos proporcionales a e−H/3e^{- H/3}, produce ⟨m⟩=0.331915\langle m\rangle = 0.331915 y ⟨∣m∣⟩=0.719730\langle|m|\rangle = 0.719730. Aunque el campo favorece una orientación, el sistema pequeño todavía visita configuraciones de magnetización negativa y la media firmada es mucho menor que la media absoluta. El segundo momento vale ⟨m2⟩=0.606902\langle m^{2}\rangle = 0.606902. Con la fórmula de fluctuaciones, la susceptibilidad canónica por sitio es χ=1.490205\chi = 1.490205 en las unidades escogidas. Estos valores provienen del recuento exacto y no dependen de una semilla.

Implementamos luego un muestreador local que selecciona sitios con reemplazo y aplica la regla de aceptación. Para hacerlo reproducible se fija la semilla 2026102120261021, se comienza con los nueve espines positivos, se ejecutan 60 00060\, 000 barridos y se descartan los primeros 10 00010\, 000. Las 50 00050\, 000 configuraciones restantes dan ⟨m⟩^=0.331618\widehat{\langle m\rangle} = 0.331618, ⟨∣m∣⟩^=0.717111\widehat{\langle|m|\rangle} = 0.717111 y ⟨m2⟩^=0.604002\widehat{\langle m^{2}\rangle} = 0.604002. La fracción de propuestas aceptadas, calculada sobre todos los intentos, es aproximadamente 0.305470.30547. El pequeño desacuerdo con la enumeración es compatible con un error de muestreo, pero esa comparación aislada no equivale a un intervalo de confianza ni valida la convergencia.

La comparación incluye pruebas de naturaleza diferente. La energía calculada mediante actualizaciones locales debe coincidir con la recalculada desde todas las parejas de vecinos. La distribución observada debe aproximarse a la canónica al prolongar el muestreo. La respuesta medida mediante fluctuaciones tiene que ser compatible con la que se obtiene derivando la función de partición exacta. En una red pequeña estas comprobaciones pueden hacerse por separado, sin acudir a una figura que «parezca» correcta.

El protocolo también permite detectar errores más sutiles. Si una implementación cuenta dos veces los enlaces, cambia el factor de Boltzmann y desplaza sistemáticamente todas las medias. Si se usa e−ΔH/Te^{- \Delta H/T} sin fijar kB=1k_{B} = 1, las unidades quedan ocultas. Si se acepta siempre un cambio de energía positiva, el algoritmo calentará artificialmente la muestra. La enumeración exacta detecta esos defectos aunque la cadena haya recorrido millones de barridos y sus gráficos sean muy suaves.

El script de verificación adjunto al capítulo permite reconstruir estos datos y conserva los parámetros, la convención de enlaces, la semilla y los valores exactos. Una ejecución corta como esta no demuestra un mezclado perfecto. El caso es una referencia controlada para auditar un método antes de usarlo en redes grandes, donde la solución exacta deja de estar disponible.

Las frecuencias de Metropolis se aproximan a las probabilidades de los 512 estados. La media acumulada fluctúa alrededor de la media exacta 0,331915.
Figura 159.2. Enumeración exhaustiva y una ejecución nueva con NumPy PCG64, semilla 20261021, 60 000 barridos y descarte de los primeros 10 000. Cada barrido contiene nueve propuestas aleatorias con reemplazo; se conservan los rechazos. El error estándar se estima con 50 bloques de 1000 barridos y no convierte el tiempo algorítmico en tiempo físico. Ampliar figura ↗

159.9. Cuándo Metropolis deja de ser eficiente

En una red ferromagnética grande, cerca de la temperatura crítica, aparecen fluctuaciones colectivas cuya longitud se aproxima al tamaño del sistema. Cambiar un espín por vez resulta poco eficaz para reorganizar regiones correlacionadas, porque muchas propuestas alteran detalles locales sin producir un estado estadísticamente muy diferente. La cadena puede cumplir rigurosamente el balance detallado y, sin embargo, conservar memoria durante miles de barridos. Esa lentitud es una propiedad dinámica del algoritmo; no invalida la distribución estacionaria, pero sí modifica el coste necesario para estimarla.

Por debajo de la temperatura crítica y con campo próximo a cero aparece otro obstáculo, porque la inversión de la magnetización global puede exigir la formación temporal de dominios separados por interfaces. El coste de una interfaz crece, de manera aproximada, con su área; en dos dimensiones, con una longitud proporcional a LL. La probabilidad de cruzar ese tipo de barrera puede hacerse extremadamente pequeña en muestreo local. El problema se agrava cuando hay coexistencia de fases de distinta densidad o estructura, como ocurre en ciertas transiciones de primer orden. Un histograma con un valle profundo entre dos picos indica que hay regiones importantes del espacio de estados difíciles de conectar.

Los algoritmos de clúster desarrollados por Swendsen y Wang y por Wolff proponen actualizar grupos correlacionados de espines. Para el Ising ferromagnético en condiciones apropiadas reducen sustancialmente la ralentización crítica. No son un reemplazo universal, ya que dependen de la estructura de interacciones y deben preservar la distribución objetivo. El «tiempo» de un movimiento de clúster es, con más claridad todavía, un tiempo computacional sin interpretación como duración física. Si el objetivo fuera estudiar nucleación real o cinética microscópica, escoger un algoritmo solo porque mezcla mejor destruiría precisamente las propiedades temporales que se intentan representar.

Hay otros caminos. El intercambio entre réplicas a distintas temperaturas puede facilitar que una configuración atraviese regiones difíciles y vuelva al valor de interés. Los métodos multicanónicos procuran aplanar barreras de probabilidad modificando el peso de muestreo y restituyendo después, por reponderación, los observables del conjunto buscado. Esas herramientas introducen nuevas decisiones de diseño y diagnósticos; usarlas sin controlar la reponderación, el soporte explorado o el intercambio entre réplicas puede ocultar sesgos nuevos bajo un método aparentemente sofisticado.

Una tasa de aceptación alta tampoco es sinónimo de muestreo eficiente. Si casi todas las propuestas producen cambios diminutos, la cadena puede avanzar lentamente aunque acepte con frecuencia. Una aceptación baja puede reflejar propuestas demasiado costosas, pero, sin otros datos, no cuantifica la independencia de los estados aceptados. La medida relevante es el coste computacional por estimación con precisión conocida; en este capítulo lo observamos a través de autocorrelaciones y comprobaciones cruzadas, sin fijar como objetivo arbitrario un porcentaje de aceptación.

159.10. De una simulación plausible a una estimación defendible

El algoritmo de Metropolis conecta un Hamiltoniano especificado con los promedios canónicos de sistemas demasiado grandes para enumerarlos. Su logro matemático central es conservar una distribución objetivo mediante movimientos locales cuya aceptación usa solo diferencias de energía. Esa propiedad puede demostrarse sin evaluar la función de partición. La demostración establece corrección estacionaria, pero nada dice sobre si una ejecución de duración finita ha visitado todas las regiones relevantes.

Para evaluar una implementación conviene empezar por un caso pequeño y exactamente enumerable. El ejemplo de nueve espines da valores de referencia para magnetización, energía y fluctuaciones, junto con una prueba directa de la expresión local de ΔH\Delta H. Superada esa etapa, aumentamos el tamaño y estudiamos la dependencia de los resultados respecto de la inicialización, el tiempo descartado y la duración total. En cada ejecución deben quedar definidos los límites de la red, las unidades termodinámicas y la convención de conteo de enlaces.

Las barras de error requieren especial atención porque las muestras no son independientes. La información efectiva depende de la autocorrelación integrada del observable, además del número de configuraciones guardadas. Cuando una trayectoria no cambia de fase, un error estándar pequeño describe únicamente las fluctuaciones del sector explorado. El informe debe reconocer explícitamente ese confinamiento. Si se reemplaza Metropolis por un algoritmo de clúster o un método de temperado, la comparación de equilibrio sigue siendo legítima, pero las escalas temporales ya no describen el mismo proceso algorítmico.

Conviene mantener separados tres enunciados. Uno afirma que el algoritmo tiene la distribución estacionaria correcta; otro, que las configuraciones obtenidas representan adecuadamente esa distribución; el último, que el modelo físico describe el sistema que nos interesa. Cada uno necesita su propia evidencia. El primero descansa en propiedades de la cadena, el segundo en diagnósticos y estimación de incertidumbre, el tercero en el contraste científico entre el Hamiltoniano y los fenómenos observados. Un resultado puede superar la primera prueba y fallar en las otras dos.

El criterio de calidad más exigente es la reproducibilidad interpretativa. Un tercero debe poder identificar qué se simuló, cuáles fueron los observables y cuánto dependen de las decisiones numéricas. La temperatura, el campo, el estado inicial, la regla de propuesta, las semillas y la definición de barrido forman parte del método. También forman parte de él los resultados negativos, como las discrepancias entre cadenas, la lentitud de mezclado y las regiones del espacio de parámetros donde la incertidumbre no se resolvió. Estas limitaciones orientan la elección de simuladores más adecuados.

Al finalizar este capítulo podemos construir un muestreador y comprobarlo contra una referencia exacta. El cierre del módulo usará esas herramientas para una cuestión distinta, la de cómo distinguir coexistencia de equilibrio, retardo cinético e influencia del tamaño finito. Allí una curva con apariencia de histéresis dejará de interpretarse como una respuesta evidente del material y se convertirá en una hipótesis que necesita controles específicos.

Referencias