# Reproducibility in molecular dynamics trajectories

Let's say I have two molecular dynamics (MD) simulators, A and B, which implement the following interface: $$\left\{\vec{x}_i(t = 0), T_i, t \right\} \rightarrow \left\{\vec{x}_i(t), \vec{F}_i(t) \right\}$$ where $\vec{x}_i$ is the position of atom/ion $i$, $T_i$ is the type of atom/ion $i$ (e.g. hydrogen or carbon), and $\vec{F}_i$ is the force on atom/ion $i$. Given a set of times, they can produce discrete trajectories: $$\mathcal{T} = \left\{ \left(t, \vec{x}_i(t), \vec{F}_i(t)\right) : t \in (dt) \mathbb{N} \right\}$$ Unless A and B are identical, there are going to be differences in the trajectories. Given a pair of trajectories, $(\mathcal{T}_A, \mathcal{T}_B)$, how should I decide if A and B are in agreement (for that particular trajectory)? or what would I have to do to convince you that $B$ was reproducing the behavior of $A$?

To give some context, B is an altered version of A, using a slightly different Hamiltonian, so they have identical interfaces. I understand that it may be useful to sample the neighborhood of initial conditions around $\vec{x}(t=0)$, but producing a trajectory is computationally intensive, so using a small number of trajectories is preferable. The regime is ab initio molecular dynamics using plane-wave DFT. I'm particularly interested in the motion of adatoms on surfaces.

• Do you have a metric of interest such as the diffusion constant of your adatom? It would be easier to answer if you gave us more specifics on how A and B differ (temperature, barostat, thermostat, initial conditions, basis sets, functionals, etc.) – Deathbreath Aug 1 '13 at 15:46
• @Deathbreath, I don't have a specific metric, and am trying to stay away from complex metrics unless there are a set that the community looks to and trusts. The difference is most like differing basis sets. – Max Hutchinson Aug 1 '13 at 15:57
• If the initial conditions are the same (same coordinates and momenta), then you can actually do RMSDs for each step to assess divergence of the paths. – Deathbreath Aug 1 '13 at 17:47
• @Deathbreath, the Lyapunov instability, mentioned in an answer, could cause problems for RMSD across multiple steps. Did you mean for a single step? What RMSD would you classify as agreement/reproduction? – Max Hutchinson Aug 1 '13 at 21:37

Any numerical differences between A and B will become exponentially larger with time (i.e. the Lyapunov instability, as discussed in Frenkel and Smit). Even a small difference due to basis set size could result in dramatic differences in the trajectory over time. So I'm not sure that a comparison between individual trajectories will be meaningful. It may be more relevant to compare some equilibrium quantity as Deathbreath suggested, such as average energy or diffusion constant.

But if you really wanted to check whether the trajectories were the same, one thing you could try is a time-reversal test. Sometimes as a test people will run a simulation for a fixed period of time, and then run it backwards (i.e. reverse the directions of all momenta). If the integrator is time-reversible (not all are), then you should get back to where you started (or at least close, depending on machine precision and the length of the run). Maybe you could try a similar strategy: run a simulation forwards with A. Then run it reverse, once with A and once with B. If method A and method B are identical, then the final coordinates of these two runs should be more or less equally close to the starting coordinates. This could be measured by the mean square displacement of the atoms, for example.

• I suspect the condition 'should be more or less equally close to the starting coordinates' is only true in a statistical sense; either trajectory could get 'lucky' and end up uncharacteristically close to the starting point, which would throw off the comparison. – Max Hutchinson Aug 5 '13 at 16:25
• Yes, that's certainly possible. I suppose that one might have to do several trials to get statistically convincing results. – Max Radin Aug 15 '13 at 1:29

First thought is to use A and B on standard validation test cases. I'm not sure what's available for validation on MD (googling "molecular dynamics validation" turned up a lot), but in CFD there's plenty of databases.

If you need to validate something very specific (it sounds like you are), then you'll need to compile some statistics to convince me that A and B are reproducing the same behavior. You'd do this by sampling an input space of initial conditions given to both methods. There are some sophisticated sampling methods available in DAKOTA that I've used; things like polynomial chaos expansion or the various design-of-experiments/uncertainty quantification capabilities can dramatically reduce the number of runs you have to do to compile meaningful statistics.

• I'm currently interested in validating for a particular initial condition, not rigorously over a wide range of conditions (that may come later). I've updated the question to reflect this. – Max Hutchinson Aug 1 '13 at 15:52

Reproducibility of the trajectory is generally impossible; even different compiler versions on the same code implementing an integrator will not reproduce each other's first step - chaotic divergence follows. Reproducibility of the resulting thermodynamic observables is a much more relevant and achievable criterion.

One simple approach is to recognize that the distribution of (say) potential energies over the trajectory is known a priori, and varies in a known way with (say) temperature. A statistically significant observation that the same (and the expected) temperature dependence of the distribution is observed with the two integrators is an excellent start to demonstrating their equivalence. See http://dx.doi.org/10.1021/ct300688p (article also available on arxiv.org; code on https://simtk.org/home/checkensemble).

• The difference in the two schemes that I am interested in comparing are Hamiltonians, not integrators. It doesn't seem like these methods take into account differences in the definitions of the energy of the particles, even the kinetic energy. I'm really looking for an answer that takes $x_i(t)$ and $F_i(t)$ as the only inputs. – Max Hutchinson Oct 10 '13 at 13:38
• Your original question is phrased in terms of "simulators," which could be more clearly phrased if the integrator is constant and only the Hamiltonian varies. – mabraham Oct 10 '13 at 19:17
• Still, the integrator plus Hamiltonian sample an ensemble - your question is that of whether each Hamiltonian with the same integrator samples the same ensemble. Reproducing the trajectory is meaningless, even if possible, unless you can assert that the trajectory samples an ensemble with the correct properties. Shirts's method does exactly this. – mabraham Oct 10 '13 at 19:28
• The primary issue is that my question is about trajectories (or dynamics), not about ensembles (or thermodynamics). I think solutions should try not to assume the same integrator or Hamiltonian, so I think the simulator name fits best. – Max Hutchinson Oct 10 '13 at 22:41
• What purpose does your trajectory serve? – mabraham Oct 10 '13 at 22:56