# MD: How to estimate number of timesteps to converge a property

If I want to equilibrate a NPT system until a property reaches its equilibrium value (let’s say the volume V), what is the common practice to have a rough estimate of the total simulation time required?

Of course you could just do a simulation and check if the value is equilibrated. But I would like to know if there is another method to estimate the number of timesteps needed prior to the calculation. That would also allow to double check and be able to assert that there is something wrong if a property takes much too long to converge.

I was thinking the fluctuation dissipation theorem might come in handy in estimating the rate at which my volume is going to relax to its equilibrium value, but I struggle to find a proper formulation.

• While this says nothing about the average case behavior, it's easy to show that in the worse case, there is essentially no such estimate. For example you simulate water just below the freezing point, where the initial condition is a liquid. This system only equilibrates when the supercooled water freezes, but this takes forever. But I doubt if any method, even including direct simulation, will consistently show that the equilibration does take forever. Jul 16, 2021 at 14:03
• This is a good point to make, but if there's a reliable way to estimate the equilibration time, it could in principle estimate "forever" in the case of the supercooled water. I don't think the fact that some systems may not equilibrate on their own in a finite or reasonnable time shows that there is no way to predict it, right? Jul 16, 2021 at 17:24
• @BarbaudJulien I don't think a nanosecond simulation can always estimate the equilibration time, even if "forever" is a valid answer. Think of another simulation of water, but now slightly above the freezing point. A nanosecond simulation cannot spot qualitative differences between the two simulations, as supercooled water is AFAIK not qualitatively different from non-supercooled water on the ns timescale. Yet non-supercooled water equilibrates in nanoseconds, if the initial state is liquid. Jul 16, 2021 at 20:08

tl;dr: simulate for greater than N*period of oscillation

To understand how I arrived at this answer consider that properties can be written as a function of the partition function Q.

For example, the average value of the volume will be

$$\langle V\rangle =\pm {\frac {\partial \ln Q}{\partial \beta P}}.}$$

The partition function for N identical classical particles can be evaluated as $$Q = \frac{1}{N! \, h^{3N}} \int d^N \textbf{q} \int d^N \textbf{p} \ \exp \Big[ {-\beta \sum \limits_{i=1}^N H(\textbf{q}, \textbf{p})} \Big]$$

where $$h$$ is Planck's constant, the bolded $$\textbf{q}$$ and $$\textbf{p}$$ are the 3-D position and momentum vectors, respectively, and $$\beta$$ is $$\frac{1}{k_B T}$$, where $$k_B$$ is the Boltzmann constant and $$T$$ is the temperature.

This is a 6N dimensional function, which is basically impossible to evaluate for any meaningful system.

Nevertheless, we can still learn a bit by considering a model system.

Consider a classical simple 1-D harmonic oscillator (HO) with mass $$m$$ and spring constant $$k$$. The hamiltonian for this system is $$H(x, p) = \frac{p^2}{2m} + \frac12 k x^2$$.

We an analytically solve this for this : $$Q = \frac{1}{N!}\frac{kT}{\hbar \omega}$$

But we can also numerically solve for this, which would require that we at least need to sample over the period of the oscillation in $$x$$ and $$p$$.  Phase space can get pretty tricky, tho. In some cases you can get stuck in a region of phase space which might skew your results. So my formula above isn't a guarantee, more like a lower bound.

• +1 it reminds me of the "Nyquist condition". Jul 16, 2021 at 19:22
• The problem is that, while the vibrations of individual molecules may be relatively uncoupled from each other, their translations and rotations are very strongly coupled, which yield very soft modes. This is a problem even if these soft modes are harmonic, in which case your formula more or less holds, but the period is the oscillation period of the softest mode, which itself grows as a function of N; it's even worse here since the softest modes are extremely anharmonic. On the other hand, a MD step actually samples all modes, so your factor N may be an overestimation. Jul 16, 2021 at 19:53
• @wzkchem5 good points. Also this doesn't answer the question of how to calculate equilibria time. Instead it attempts to answer the question assuming you are already in equilibrium in order to calculate equilibrium properties. I don't know of any good rule of thumb to determine this a priori. Jul 16, 2021 at 20:37

On the mathematical side, there exists no algorithm that estimates the time required to converge a property within any predefined accuracy, however low the accuracy is. The reason is simple. Suppose you perform a series of simulations of a liquid, with successively low temperature. When you are above the freezing point of the liquid, your simulation probably well equilibrates within the nanosecond timescale. But as soon as you reach below the freezing point, since the thermodynamically stable phase is now the solid, your property can only be considered converged when the now supercooled liquid freezes - and this may easily take days, even if you are as much as 10 K below the freezing point! So in order to have a time estimate that is always reliable, you have to be able to predict the freezing point of the liquid to arbitrary precision - this holds even if the estimate is only required to be accurate to $$\pm7$$ orders of magnitude (from nanoseconds to days)! No matter how accurate your estimate of the freezing point is, your estimate of the equilibration time will fail by tens of orders of magnitude, if your simulation temperature happens to be between your estimate of the freezing point and the true freezing point.

The take home message is, any technique that can possibly work, will only work if you are far from a phase transition, and if your initial condition is not in a metastable state. These two conditions should preferably be verified by experimental data of this system and/or related systems.

If these two conditions are already shown to hold, the time required to converge the volume can probably be estimated by comparing to literature values of materials that have similar compressibility as your system. Compressibility describes the sensitivity of volume as a function of pressure, and should correlate well with the time required for the volume to equilibrate, since the primary reason for the volume to fluctuate is the local fluctuations of the pressure.

I wasn't going to answer, because there is no answer, but, anyways...

## First: common practice is...

Common practice is to frequently measure the property you are after, as well as others, such as internal energy, density and make a plot. The property will fluctuate but eventually reach an equilibrium that is easy for the eye to spot on a plot. Or, use the maths as Chodera has done paper here, open access, thanks BioRxiv, and let the computer tell you when you have reached equilibrium.

This is what is done. either manually or algorithmically check that the property is fluctuating about some mean value with a standard deviation you would consider indicative of equilibrium.

This, however, doesn't tell you how long it will take, it only tells you how long it took.

## Second:...

The idea I propose, which is crude, possibly too crude, is a rather simple one...

If you know what something value is, and you know its rate of change, and you know where you want it to be, you can figure out how long it will take. The problem is really knowing the rate of change. If we take a first order expansion of pressure as a (made up) function of time we could say we have

$$\begin{equation} P_{\rm equilibrium}(t) = P_{\rm current}(t_{curr}) + \left(\frac{\partial P}{\partial t}\right)_{t_{curr}} \Delta t \end{equation}$$

It would be simple in theory to estimate the rate of change of pressure by measuring the pressure at two different time steps

$$\begin{equation} \left(\frac{\partial P}{\partial t}\right)_{t_{curr}} \approx \frac{P_2 - P_1}{t_2 -t_1} \end{equation}$$ So you could rearrange for the time

$$\begin{equation} \Delta t = \frac{P_{equil} - P_{curr}}{\frac{P_2 - P_1}{t_2 -t_1}} \end{equation}$$

You could probably replace P for pressure with P for property.

Unfortunately this probably won't work, particularly in the case of pressure, since instantaneous pressure can jump wildly. You would need to calculate the rate of change many times and take the average.

Another reason it likely wouldn't be that accurate would be that the rate of change would probably not be constant. It would be like when a download says 5 minutes remaining, and then in 1 minute it says 10 minutes remaining, and then in another minute, 30 minutes remaining, and then in another minute Download Complete.

Anyways, that is a thought I had.

In the event you were doing an NPT and wanted to guess what the volume would be at equilibrium, you could probably apply chain rule, and come up with a scheme, where there would be a partial derivative that is essentially the compressibility factor.