How to tell convergence in Metropolis Monte Carlo simulation?

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.