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

Question

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.

Answer 1

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

${\displaystyle \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.

Answer 2

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.

Answer 3

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.