GROMACS Coordinates and velocity and force trajectories do not match accuracy

Question

I am using GROMACS to perform MD simulations of proteins. So, I tried to record the simulation results nstxout, nstvout, and nstfout with the same step size (nstxout 1, nstvout 1, nstfout 1). Then, I calculated the coordinate (step n) of an arbitrary step of the trajectory (step n) from the velocity of step n - 1, step n, and calculated the velocity of n - 1/2, and used that velocity to calculate the coordinate (n-1) of step n using the velocity of n - 1/2. LeapFrog was used in the simulation.

How can I get the coordinates, velocity, and force to be consistent?

Answer 1

From GROMACS documentation about reproducibility:

The following factors affect the reproducibility of a simulation, and thus its output:

• Precision (mixed / double) with double giving “better” reproducibility. Number of cores, due to different order in which forces are accumulated. For instance (a+b)+c is not necessarily binary identical to a+(b+c) in floating-point arithmetic.
• Type of processors. Even within the same processor family there can be slight differences.
• Optimization level when compiling.
• Optimizations at run time: e.g. the FFTW library that is typically used for fast Fourier transforms determines at startup which version of their algorithms is fastest, and uses that for the remainder of the calculations. Since the speed estimate is not deterministic, the results may vary from run to run.
• Dynamic linking to different versions of shared libraries (e.g. for FFTs)
• Dynamic load balancing, since particles are redistributed to processors based on elapsed wallclock time, which will lead to (a+b)+c != a+(b+c) issues as above
• Number of PME-only ranks (for parallel PME simulations)
• MPI reductions typically do not guarantee the order of the operations, and so the absence of associativity for floating-point arithmetic means the result of a reduction depends on the order actually chosen
• On GPUs, the reduction of e.g. non-bonded forces has a non-deterministic summation order, so any fast implementation is non-reprodudible by design.

I would add that in your case of hard coding everything by hand you may have to look a bit into your barostat and thermostat (although their influence should be found on speed of the atoms). Maybe a sample code with coordinates and speed of 4 atoms would help to give a more precise answer.

You can also try having a look at gromacs source code. I believe OpenMM provides a bit more control and clarity on each step, if for some unfathomable reason you are not doing this only for curiosity sake but because you actually have to make a code that reproduces these results.