# In a Monte Carlo $NVT$simulation how do I determine equilibration?

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(\frac{1}{r^{12}} - \frac{1}{r^6})$, 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)

• This is a good question, but the answer in the scientific literature for when a system has equilibriated is "when you feel like it". Indeed you can easily find people publishing results of new phenomena, when in fact they are merely artefacts of nonequilibrium. Even when you have reached a "steady state", you are not guaranteed to have found a representative state of the free energy minimum. An accidentally supercooled liquid (e.g. through finite size effects) serves as an example. – alarge Jan 19 '16 at 21:27
• interesting other post with the same question but in a theoretical context Functions and Length Scales – user46925 Jan 22 '16 at 16:04
• That question is different than mine, because my equilibrium has fluctuations... – Adi Ro Jan 22 '16 at 20:23
• It's all in my 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(1r12−1r6)U=4(1r12−1r6), in reduced units). boundary conditions are periodic." I don't understand what are saturation parameters. the graphs show 500000 runs – Adi Ro Jan 23 '16 at 15:06
• @igael saturation is when thermal equilibrium is reached. In a MonteCarlo simulation we try to estimate an integral: $\int{e^{-\frac{u}{kT}}udu}$ where the integral is over all possible states. so we create a set of systems that represent the most likely configurations (the thermal equilibrium configurations). we give each configuration a probability of $e^{-\frac{u}{kT}}$ to reflect the exponential term. so "saturation" is reached when we have simulated enough configurations of the system, so that we can say that our finite sum reflects the infinite integral. I'm asking when is it reached – Adi Ro Jan 23 '16 at 17:40

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://www.biorxiv.org/content/early/2015/07/04/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.

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: $$<P>, <\Delta P^2>$$

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(<P>_i , <P>_{i-1},...,<P>_{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.

• I'm not sure that I understand you correctly. when you say "mean value" ($<P>$) do you mean an average over N steps of the simulation? what do you mean when you say "temporal"? there is no time involved... – Adi Ro Jan 21 '16 at 18:22
• I mean you do an average on some values of <P> during the last N steps (this N is something you decide, the last 10 steps or 20) and this average value needs to converge. So you save some average values to check again the variance that should be a small number (again another decision that you should do, 1e-3? 1e-2? 1e-5?). I know it's chaotic the fact that you do more statistics with some statistical data... – Jose M Gonzalez Jan 23 '16 at 13:32
• And I apology because I wrote medium value when I meant mean value and thanks for the edition – Jose M Gonzalez Jan 23 '16 at 13:37
• Ok now I understand you. It seams logical that the variance should go to zero, but what if doesn't have to? Because we have temperature maybe in equilibrium, the system fluctuates with constant larger than zero variance? I't even logical to think that the variance is not constant! It's thermal fluctuations after all – Adi Ro Jan 25 '16 at 19:01

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: $$<A(t)>_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{<A(k)>_T - <A(\infty)>_T}{<A(0)>_T - <A(\infty)>_T}$$ Where $<A(\infty)>_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.

• With respect to your question: The typical relaxation behaviour is something exponential like $\phi=\exp{-t/\tau}$. You see that you can infer $\tau$ by integration. The problem (and Im sure Mr. Binder points that out) is: How to choose $A$? You could simply search for the largest relaxation times, but this will turn out not to be practical. So you will have to find the timescale on which the processes relevant for your problem have relaxated. – Bort Feb 9 '16 at 15:37

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.

I hope this may help you, Susana (Engineering Physics student, Polythectic University of Catalonia)

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).