¿Cómo determino el equilibrio en una simulación Monte Carlo NVTNVTNVT?

Veo que esta pregunta ha sido trasladada a la página de inicio nuevamente. Recientemente respondí una pregunta similar en https://math.stackexchange.com/a/2920136/575517 , así que daré la esencia aquí, en caso de que alguien lo encuentre útil.

Esta pregunta "¿Cómo sé que la ejecución de la simulación ha alcanzado el equilibrio?" a menudo se pasa por alto y se deja como una "regla general". Como se discutió en otras respuestas, generalmente se requiere que un tiempo de "equilibrado" o "quemado" sea al menos tan largo como el tiempo de correlación$\tau_A$de la variable$A$de interés o, mejor aún, todas las variables de interés (tomando la mayor$\tau$). Uno descarta esa parte de la trayectoria y comienza a acumular promedios a partir de entonces.

Los problemas aquí son que uno necesita una estimación de$\tau_A$de antemano, y también que este argumento se basa libremente en la teoría de la respuesta lineal : es decir, que la relajación de un estado ligeramente perturbado al equilibrio ocurre en una escala de tiempo dada por$\tau_A$, que es una propiedad de la función de correlación de tiempo de equilibrio. No hay garantía de que la relajación de una configuración inicial arbitrariamente preparada siga esta ley, aunque$\tau_A$podría proporcionar una guía razonable.

Sin embargo, conozco al menos un artículo en el que John Chodera intentó abordarlo de manera objetiva: https://doi.org/10.1101/021659 , que también se publicó en J Chem Theo Comp, 12, 1799 (2016) .

No intentaré reproducir las matemáticas aquí, pero la idea básica es usar el procedimiento para estimar errores estadísticos en secuencias de datos correlacionados, lo que implica estimar el tiempo de correlación (o la ineficiencia estadística, que es el espacio entre muestras efectivamente independientes). ) - y aplicándolo al intervalo$(t_0,t_{\text{max}})$que cubre el período entre el final (propuesto) del período de equilibrio,$t_0$, y el final de todo el conjunto de datos,$t_{\text{max}}$. Este cálculo se comporta de forma predecible si el conjunto de datos está en equilibrio: las fluctuaciones$\langle\delta A_{\Delta t}^2\rangle$de un promedio de tiempo finito

$$\delta A_{\Delta t} = A_{\Delta t}-\langle A\rangle, \qquad A_{\Delta t} = \frac{1}{\Delta t} \int_t^{t+\Delta t} dt \, A(t)$$ depender de manera conocida del intervalo de promediación$\Delta t$. Los orígenes del tiempo$t$se eligen dentro del intervalo de interés,$(t_0,t_{\text{max}}-\Delta t)$; las fluctuaciones se calculan a partir de todos los periodos$\Delta t$que se encuentra dentro de ese intervalo. El método realiza sistemáticamente este cálculo en función de$t_0$, reduciéndolo desde un valor alto inicial hacia el inicio de la ejecución. En algún momento, se espera que se observen desviaciones del comportamiento esperado. Ese punto se considera el final del período de equilibrio.

De todos modos, leer ese documento debería ser útil para aclarar este problema. El autor también proporciona una pieza de software de Python para implementar el cálculo automáticamente, por lo que también puede ser de uso práctico.

Como dicen muchos comentarios, no hay una única y mejor respuesta, cada uno usa un método diferente. La solución que encontraste es buena, pero ¿cómo defines cuándo se ha alcanzado el equilibrio?

Para ello necesitas comprobar los últimos valores de la simulación (Energía, presión, etc.), por lo que eliges un conjunto de configuraciones previas que comprobarás:

$$N = 10$$

Y con los parámetros que definen tu equilibrio calculas el valor medio y la desviación estándar:

$$\langle P\rangle, \langle \Delta P^2\rangle$$

Esos valores no deberían cambiar demasiado después de algunos pasos. si almacena algunos valores medios y su varianza, verá que la media converge al valor del sistema y la varianza sobre los valores medios temporales en cada paso irá a cero.

$$Var(\langle P\rangle_i , \langle P\rangle_{i-1},...,\langle P\rangle_{i-n})\rightarrow0$$

Por lo tanto, lo que debe elegir como parámetro para el equilibrio es cuántos pasos considera para el valor medio y cuántos valores medios usa para calcular la varianza.

PD: Las fluctuaciones después de que el sistema ha alcanzado el equilibrio son normales y también que las fluctuaciones aumentan con la temperatura.

No puedes saberlo con seguridad.

En la simulación, el sistema cambia constantemente de manera caótica y, en general, es posible que después de algún período de relativa constancia de las variables macroscópicas, ocurra algún cambio radical, como una transición de fase.

Sin embargo, si el sistema es muy simple, como un gas de esfera dura, uno puede calcular teóricamente cómo debería ser el estado de equilibrio para restricciones dadas como el volumen y la temperatura. Entonces, si los valores derivados de la simulación se aproximan a esos valores teóricos, podemos decir que lo más probable es que estemos observando una transición al equilibrio y no va a ocurrir nada sorprendente.

Si no conocemos los valores de las variables en equilibrio termodinámico, no podemos decir que el sistema se está acercando a él.

Sin embargo, también existe la idea de equilibrio termodinámico restringido, definido solo para un período de tiempo limitado que nos interesa. Entonces, el sistema está en equilibrio, si la estimación de la energía cinética promedio es uniforme en toda la región del espacio del sistema y si se observa la fluctuación en el tiempo y el espacio es consistente con la distribución de Boltzmann, y si lo mismo es cierto para otras cantidades físicas de interés (presión promedio, energía potencial neta).

Creo que es totalmente normal tener mayores fluctuaciones de las funciones termodinámicas al aumentar la temperatura. De lo que no estoy tan seguro es si eso pertenece a una situación realista o no...

Si conoce la forma de la función de distribución de partículas resultante, puede adivinar si su sistema ha alcanzado el equilibrio mediante el cálculo de su promedio, haciendo el cálculo de la función de distribución en cada n pasos de tiempo (300 pasos de tiempo ha sido un buena n en mis simulaciones de Monte Carlo).

Debe tener en cuenta que se supone que una simulación de Monte Carlo ya está simulando un sistema de equilibrio, por lo que no debe esperar observar ningún fenómeno dinámico al usarla (no puede observar la aparición de burbujas en el grueso de una mezcla, por instancia).

Si está configurando el sistema para un amplio rango de temperaturas, cambiar el parámetro de desplazamiento podría ser una buena idea ya que la tasa de aceptación depende de la relación (energía nueva-energía antigua)/T.

Encontré una respuesta en el libro Simulaciones de Monte Carlo en física estadística: una introducción (por K.Binder, DWHermann), página 35.

Para determinar el equilibrio, necesitamos ejecutar la simulación varias veces, digamos$n_{run}$veces. definimos el promedio$<>_T$como promedio después$t$pasos de la simulación:

$$\langle A(t)\rangle_T = \frac{1}{n_{run}} \sum_{l=1}^{n_{run}}{A(t,l)}$$ Donde A es alguna propiedad física, la energía o la presión por ejemplo.$A(t,n)$es el valor de esta propiedad en el tiempo (paso de simulación) t de la n-ésima simulación. Ahora definimos la "función de relajación no lineal": $$\phi_A(k)=\frac{\langle A(k)\rangle_T - \langle A(\infty)\rangle_T}{\langle A(0)\rangle_T - \langle A(\infty)\rangle_T}$$ Dónde$\langle A(\infty)\rangle_T$es el promedio de A en el último paso de simulación, para$n_{run}$simulaciones Esta función de relajación llega a cero después de algún paso t. el tiempo de relajación de esta función es: $$\tau_A=\int_0^{\infty}\phi(t)dt$$ Todavía no estoy seguro de por qué este es el tiempo de relajación, me encantaría que alguien pudiera explicarme. Entonces, si este es el tiempo de relajación, entonces el tiempo durante el cual el sistema alcanza el equilibrio debe ser mucho más grande que este tiempo: $$t_{simulation} \gg \tau_A$$ dado que diferentes propiedades pueden tener diferentes tiempos de relajación, debemos considerar el tiempo de relajación más lento para equilibrar el sistema.