How do I determine equilibration in a Monte Carlo $NVT $ simulation?

Question

I'm running an NVT (constant number of particles, volume and temperature) Monte Carlo simulation (Metropolis algorithm) of particles in two dimensions interacting via Lennard-Jonse potential $$U = 4\left(\frac{1}{r^{12}} - \frac{1}{r^6}\right),$$ (in reduced units). Boundary conditions are periodic.

From this simulation I'm calculating the instantaneous pressure and potential energy. In the first steps the system is not in equilibrium, so I need to start averaging after the system is in equilibrium.

I'm starting my simulation from a random configuration.

My question: even after the system has reached equilibrium, it fluctuates around this equilibrium. These fluctuations may be large for large temperatures. So how do I know that I have reached equilibrium?

Here are some examples of the curve:

$$\uparrow$$ The energy Vs. simulation step, for a high temperature (warmer color is higher density)

$$\uparrow$$The energy Vs. simulation step, for a low temperature (warmer color is higher density)

$$\uparrow$$The energy Vs. simulation step, for a high temperature, only for low densities. in this graph it's harder to tell if we reached equilibrium (warmer color is higher density)

Answer 1

I see that this question has been bumped to the home page again. I recently answered a similar question at https://math.stackexchange.com/a/2920136/575517 so I'll give the essence of it here, in case anyone finds it helpful.

This question "How do I know that the simulation run has reached equilibrium?" is often glossed over and left as a "rule of thumb". As discussed in other answers, usually one requires that an "equilibration" or "burn-in" time should be at least as long as the correlation time $\tau_A$ of the variable $A$ of interest or, even better, all variables of interest (taking the longest $\tau$). One discards that part of the trajectory, and starts accumulating averages thereafter.

Problems here are that one needs an estimate of $\tau_A$ beforehand, and also that this argument is loosely based on linear response theory: namely that the relaxation of a slightly perturbed state to equilibrium occurs on a timescale given by $\tau_A$, which is a property of the equilibrium time correlation function. There's no guarantee that the relaxation from an arbitrarily prepared initial configuration will follow this law, even though $\tau_A$ might provide a reasonable guide.

However I'm aware of at least one paper where an attempt has been made, by John Chodera, to tackle it objectively: https://doi.org/10.1101/021659 which was also published in J Chem Theo Comp, 12, 1799 (2016).

I won't try to reproduce the mathematics here, but the basic idea is to use the procedure for estimating statistical errors in correlated sequences of data - which involves estimating the correlation time (or the statistical inefficiency, which is the spacing between effectively independent samples) - and applying it to the interval $(t_0,t_{\text{max}})$ which covers the period between the (proposed) end of the equilibration period, $t_0$, and the end of the whole dataset, $t_{\text{max}}$. This calculation behaves in a predictable way if the dataset is at equilibrium: the fluctuations $\langle\delta A_{\Delta t}^2\rangle$ of a finite-time-average $$ \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) $$ depend in a known way on the averaging interval $\Delta t$. The time origins $t$ are chosen within the interval of interest, $(t_0,t_{\text{max}}-\Delta t)$; the fluctuations are calculated from all the periods $\Delta t$ lying within that interval. The method systematically carries out this calculation as a function of $t_0$, reducing it from an initial high value towards the start of the run. At some point, it is expected that deviations from the expected behaviour are seen. That point is taken to be the end of the equilibration period.

Anyway, reading that paper should be helpful in clarifying this issue. The author also provides a piece of Python software to implement the calculation automatically, so it may be of practical use as well.

Answer 2

As many comments say, there is not a single and best answer, each one uses a different method. The solution that you found is a good one, but how do you define when the equilibrium has been reached?

In order to do that you need check the last values of the simulation (Energy, pressure, etc.), so you choose a set of previous configurations that you'll check:

$$N = 10$$

And with the parameters that define your equilibrium you compute the mean value and the standard deviation: $$\langle P\rangle, \langle \Delta P^2\rangle$$

Those values shouldn't change too much after some steps. if you store some mean values and their variance you will see that the mean converges to the value of the system and the variance over the temporal mean values at each step will go to zero. $$Var(\langle P\rangle_i , \langle P\rangle_{i-1},...,\langle P\rangle_{i-n})\rightarrow0$$

Thus what you need to choose as a parameter for the equilibrium are how many steps you consider for the mean value and how many mean values you use to compute the variance.

P.D: The fluctuations after the system has reached the equilibrium are normal and also that fluctuations increases with temperature.

Answer 3

You can't know for sure.

In the simulation, the system is constantly changing chaotically and in general it is possible that after some period of relative constancy of macroscopic variables, some radical change such as phase transition will occur.

However, if the system is very simple, such as hard sphere gas, one can calculate theoretically how the equilibrium state should look like for given constraints such as volume and temperature. Then, if the values derived from the simulation approach those theoretical values, we can say that most probably we are observing a transition to equilibrium and nothing surprising is going to happen.

If we do not know the values of variables in thermodynamic equilibrium, we can't say the system is approaching it.

However, there is also the idea of restricted thermodynamic equilibrium, defined only for some limited time period that we are interested in. The system is then in equilibrium, if estimation of average kinetic energy is uniform throughout the system's region of space and if its observed fluctuation in time and space is consistent with Boltzmann's distribution, and if the same is true for other physical quantities of interest (average pressure, net potential energy).

Answer 4

I think it is totally normal to have greater fluctuations of the thermodynamic functions when increasing temperature. What I am not so sure about is whether that belongs to a realistic situation or not...

If you know the shape of the resulting particle distribution function, you can guess if your system has reached the equilibrium by means of the computation of its average, making the distribution function calculation at each n time-steps (300 time-steps has been a good n in my Monte Carlo simulations).

You have to take into account that a Monte Carlo simulation is already supposed to be simulating an equilibrium system, so you should not expect to observe any dynamic phenomena when using it (you cannot observe the appearance of bubbles in the bulk of a mixture, for instance).

If you are setting the system for a large range of temperatures, changing the displacement parameter could be a good idea since the acceptation rate depends on the ratio (new energy-old energy)/T.

Answer 5

I found an answer in the book Monte Carlo simulations in statistical physics - an introduction (by K.Binder, D.W.Hermann), page 35.

To determine equilibration we need to run the simulation a few times, let's say $n_{run}$ times. we define the average $<>_T$ as an average after $t$ steps of the simulation: $$\langle A(t)\rangle_T = \frac{1}{n_{run}} \sum_{l=1}^{n_{run}}{A(t,l)}$$ Where A is some physical property, the energy or the pressure for example. $A(t,n)$ is the value of this property in time (simulation step) t of the n'th simulation. Now we define the "non-linear relaxation function": $$\phi_A(k)=\frac{\langle A(k)\rangle_T - \langle A(\infty)\rangle_T}{\langle A(0)\rangle_T - \langle A(\infty)\rangle_T}$$ Where $\langle A(\infty)\rangle_T$ is the average of A in the last simulation step, for $n_{run}$ simulations. This relaxation function goes to zero after some step t. the relaxation time of this function is: $$\tau_A=\int_0^{\infty}\phi(t)dt$$ still not sure why this is the relaxation time, I will be happy if someone could explane it to me So if this is the relaxation time, than the time for which the system reaches equilibrium must be much larger than this time: $$t_{simulation} \gg \tau_A$$ since different properties may have different relaxation times, we must consider the slowest relaxation time to get the system equilibrated.