phase average uncertainty in molecular dynamics

I was wondering if I could find some scaling relation of the uncertainity of phase averages of thermodynamic quantities with the number of atoms $$N_a$$, the number of timesteps $$N_t$$ and size of the timesteps $$dt$$ in a molecular dynamics (MD) simulation.

For example, I have learned that energy fluctuations $$\Delta E$$ in the canonical ensemble scale with $$N_a$$ as: $$$$\frac{\Delta E}{E} \sim \frac{1}{\sqrt{N_a}}$$$$ If I understand correctly, $$\Delta E$$ does not depend on the simulation length. Instead, larger simulation lengths for a fixed $$N_a$$ will only decrease our uncertainity on $$\Delta E$$ but not the actual value of $$\Delta E$$. However, I'm confused because I understand that in MD one can enlarge $$N_a$$ to decrease uncertainty on phase averages for the same simulation length. I believe this is desired as we can process more atoms in parallel, but not more timesteps in parallel. Is there some general relation like the one above where both the number of atoms and the simulation length is also present?

To think about statistical mechanics, one must very occasionally¹ think about statistics.

So imagine I toss a coin six times, and you toss a coin six hundred times. My coin comes up heads three times, and your coin comes up heads three hundred times.

You would be much more confident that your coin is fair (that is, it has a 50-50 chance of coming up heads) than mine. In fact (using a frequentist approach) you would be about ten times more confident.

And yet, the standard deviation of your sample and mine (using the population formula) are exactly the same! (Don't take my word for it -- calculate it.)

What differs between my sample and yours is that the standard error of the estimator of the mean is much smaller in your sample and mine.

Similarly, if I sample a molecular dynamics system for 10 ns of time, and then for 1000 ns of time, the standard deviation of energy measured along the trajectory should not change. (Up to sampling error. Variance, and thus standard deviation, converges much more slowly than the mean.) However, the standard error of estimated mean energy does become 1/10th of its former value, varying inversely with sample size.

You should then be asking how to calculate the sample size along an MD trajectory. After all, let's say I run a 10 ns simulation and record data every 5 fs instead of every 20 fs. I hope your intuition is sharp enough to recognise I don't actually have four times as much data!

The correct answer is that the "sample size" of a simulation, for each observable, is the simulation time divided by the observable's correlation time. Sampling the simulation more often than the correlation time does not increase the sample size. For further exploration, see Grossfield's paper on quantifying uncertainty in biomolecular simulations.

Note also that just as there are correlation times, there are also correlation lengths. Most MD simulations are large enough that correlation lengths can be ignored, but occasionally some observables are affected enough by spatial correlation that finite size corrections are needed. Two notable examples are the power spectra of transverse vibrations of lipid bilayer membranes, and the Yeh-Hummer correction for estimating diffusivity from mean-squared displacement.

¹this is a playful, extreme understatement

• Thank you very much for your answer. Now considering we are sampling data not more often than the correlation time, would increasing the number of atoms used in our ensemble decrease our uncertainty on the average value of the energy for the same simulation length and sample size?
– n1ps
Jan 24 at 17:17
• In the provided reference I have seen that the standard error of the mean of some observable $f$ is $SE(f)=\sigma_f/\sqrt{N_f^{ind}}$ where $\sigma_f$ is the variance of the observable and $N_f^{ind}$ is the number of independent samples. So if I understand correctly, for an ensemble with more atoms, $\sigma_f$ will decrease and we won't need as many $N_f^{ind}$ to keep same SE.
– n1ps
Jan 24 at 17:33
• this is EXTREMELY dependent on your exact observable. "A particle's energy" is not uniquely defined (you can arbitrarily re-partition interaction energy between any of the interacting particles, and their mediating field, without changing the observed dynamics). The average energy per particle is a property of the configuration, not of individual particles, and so it is observed once per configuration and not per particle. Jan 30 at 0:58
• But suppose you are modelling a box of flexible water molecules. The H-O-H angle of each molecule is well-defined and simulating 2x water molecules, for the same amount of time, should indeed yield ~ 0.7x standard error in estimated angle (with simulation size and time far exceeding correlation scales). But: simulating 2x water molecules will usually be >2x slower and the tradeoff is not worthwhile, in terms of computational cost. Jan 30 at 1:03
• Thanks a lot again, this solves the questions I had regarding this topic.
– n1ps
Jan 30 at 12:17