# How to tell convergence in Metropolis Monte Carlo simulation?

I just started learning MC simulation and I wrote my own code for it based on Lennard-Jones potential.

Per cycle:

1. Move particle $$i$$
2. Apply boundary condition
3. Calculate change in energy, $$\Delta E$$ due to this move
4. If $$\Delta E < 0$$, accept this move, elseif rand$$<\exp(-\Delta E/kT)$$, accept this move.
5. Steps $$1-4$$ are looped through a system of $$216$$ methane particles. Then record the energy of the system.

People at school told me that since it is stochastic, I will not be able to tell convergence by looking at how the energy changes over cycle. I should be looking at acceptance rate - I could not tell anything from this either.

I adjusted my displacement step for a system at $$300K$$ to get an acceptance rate of around $$50\%$$. But I am still getting fluctuation that does not look minimal to me whatsoever. I ran my simulation for $$5000$$ and $$10000$$ cycles but both look just the same.

• One test is to run multiple chains with different random seeds, and see if you get approximately the same distribution with each chain. Mar 22, 2021 at 2:47
• Hi Andrew, in the context of molecular simulation, what do chain and seed refer to? Mar 22, 2021 at 2:48
• A chain is a specific instance of the MC sampler. So, "running two chains" means run the MC algorithm twice. The "seed" is the initial state of the pseudo-random number generator used by the algorithm. Essentially the seed defines a specific realization of the stochastic process. If your run has converged, then the distribution shouldn't depend on the seed (the algorithm should have forgotten its initial state and converged to the true distribution). Running multiple chains with different seeds directly tests whether this criterion (ie: the answer shouldn't depend on the seed) is met or not. Mar 22, 2021 at 2:55
• As a caveat my experience is in statistical sampling of Bayesian posterior distributions for parameter estimation (Markov Chain Monte Carlo), and I don't know whether or not there are better diagnostics for a molecular dynamics simulations. Mar 22, 2021 at 3:11
• Perhaps a better fit for the Computational Science Stack Exchange? Mar 22, 2021 at 5:05

I would start from the very beginning. What the Monte Carlo (MC) algorithm is supposed to do, and what do we mean by the convergence of a simulation.

Metropolis MC (MMC) is one algorithm of the family of MC algorithms, all characterized by the use of random variables. The goal of the specific MMC you are using is to generate a sequence of atomic configurations (positions) with the property that asymptotically (in the limit of a huge number of steps) the frequency of each configuration approaches the probability the same configuration would have in an equilibrium canonical ensemble.

As a consequence of this basic requirement of MMC chains, averages of physical observables over the set of configurations obtained in a simulation will converge towards the canonical equilibrium averages. It should also be clear that, due to the underlying random sampling of configurations, such a convergence is a non-monotonic, probabilistic convergence. Therefore, the instantaneous values of the potential energy of the $$k-$$th configuration ($$U_k$$) will keep fluctuating forever around the asymptotic value with a non-zero variance. Such a non-zero variance is an unavoidable feature of the physical system: in the canonical ensemble, the kinetic, potential, and total energies fluctuate around their average values. Moreover, such fluctuations are connected to physical quantities. For instance, fluctuations of total energy are connected to the constant-volume specific heat of the system.

The running average of these quantities shows different behavior. The quantity $$P_i = \frac{1}{i}\sum_{k=1}^i U_k,$$ at the end of a long run ($$i \rightarrow \infty$$) will give the average potential energy over the generated configurations. For longer and longer runs, it will approach the canonical average of the potential energy.

It is on the approach of the running average to its asymptotic value that it is meaningful to speak about convergence.

There are some consequences of this state of things that are usually not appreciated by beginners in the field. The most important is that, while the asymptotic value (i.e., the canonical average) of a running average is unique, the convergence to this asymptotic value depends on the starting value. The less likely is the starting configuration as measured by its canonical probability, the longer the run is required to achieve a given threshold of accuracy. This is the reason for not starting the averages at the beginning of the simulation. In practice, all simulations are split into one or more equilibration runs, followed by the runs where averages are evaluated (statistics runs).

Deciding where the equilibration part of a simulation is finished requires some training. Still, as a general indication, a look at the sequence of instantaneous values should fluctuate around a central stationary value without indications of global trends of variation. For example, the first plot included in the question can be considered an example of an equilibrated sample. Starting an average from any of the corresponding configurations is probably safe. Then, usual techniques to evaluate the statistical uncertainty of the averages should be used to decide when a required threshold of accuracy has been achieved.

In principle, that's all. But what about the acceptance ratio?

Well, that is a related but different story. The ratio between accepted moves and attempted moves is an indirect measure of the efficiency of the MMC chain and, in particular, of the choice of the Markovian transition matrix. In most implementations, the acceptance ratio is related to the maximum change of the coordinates (maximum displacement) allowed at each step. There are two trivial limiting cases:

• a zero maximum displacement would result in a useless simulation where all configurations are accepted (100% acceptance) just because nothing change;
• a large maximum displacement would correspond to completely random sampling. In the case of an atomic system, this is highly inefficient: in almost all the cases, the new configuration would imply short interatomic distances, then huge potential energy, and consequently, almost all the moves will be rejected (0% acceptance).

Therefore, one has to find a non-zero displacement providing a non-negligible acceptance.

The optimal acceptance depends on the precise meaning of efficiency. The often-cited value of 50% is a traditional number coming from the pioneering studies, where there was no careful evaluation of the computational efficiency of the algorithm. More accurate studies have shown that a lower acceptance ratio (say between 10% and 20%) is better since it corresponds to a larger diffusion of the system in the configuration space in a given number of steps. Such larger diffusion, in turn, implies a smaller statistical error. However, the increase of efficiency between acceptances of 50% and 15%, although measurable, is smaller than a factor 2 for typical simple liquids. Therefore, using the traditional 50% is not a huge penalty. In case of doubts, evaluating the root mean square displacement from the initial configuration is a good way to estimate the effect of different choices.

It should be clear that the link between acceptance ratio and convergence is only indirect and related to the efficiency of the sampling.

• That is everything I needed, thank you! Mar 22, 2021 at 19:16