Project

General

Profile

Bug #1372

Energy drift - comparison of double and single precision

Added by guillaume chevrot over 6 years ago. Updated over 5 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
documentation
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I compared the total energy of 2 simulations:
lysozyme in water / NVE ensemble / single precision / Gromacs 4.6.3
lysozyme in water / NVE ensemble / double precision / Gromacs 4.6.3

I observe a constant drift in energy in the case of the single precision simulation (see file total_energy_lys.pdf)

I enclose the ‘mdout.mdp’ file so you can check the mdp options I used. Note: in this example, I used the velocity Verlet integrator. I also tested the leap-frog integrator and I observed the same drift in energy.

I also join the log, edr and tpr files.

Guillaume

PS: In the literature, this problem has already been mentioned: http://jcp.aip.org/resource/1/jcpsa6/v126/i4/p046101_s1

mdout.mdp (11 KB) mdout.mdp guillaume chevrot, 10/28/2013 07:34 PM
total_energy_lys.pdf (32.5 KB) total_energy_lys.pdf guillaume chevrot, 10/28/2013 07:34 PM
md_nve_lysspce_0_10ns.tpr (2.02 MB) md_nve_lysspce_0_10ns.tpr guillaume chevrot, 10/28/2013 07:34 PM
md_nve_lysspce_0_10ns.log (3.3 MB) md_nve_lysspce_0_10ns.log guillaume chevrot, 10/28/2013 07:34 PM
md_nve_lysspce_0_10ns.edr (3.09 MB) md_nve_lysspce_0_10ns.edr guillaume chevrot, 10/28/2013 07:35 PM
lys.top (552 KB) lys.top guillaume chevrot, 10/29/2013 07:11 PM
npt.gro (3.38 MB) npt.gro guillaume chevrot, 10/29/2013 07:12 PM
Water.tar.bz2 (281 KB) Water.tar.bz2 guillaume chevrot, 10/29/2013 07:12 PM

Related issues

Related to GROMACS - Task #1063: update documentation to correct "single" precisionClosed12/06/2012

History

#1 Updated by Mark Abraham over 6 years ago

  • Status changed from New to Feedback wanted

Thanks for the report. This seems to be the same system discussed on gmx-users already (https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-October/084733.html).

The complexity of this system is a bit high to be a useful diagnostic of where a problem might lie. SETTLE is known to drift acceptably under some conditions, but the point is moot because Guillame's literature comparisons are with water-only simulations and his buffer is so long that it is off the scale of that reference (and his 1fs time step is half of theirs). If Guillame can show excess drift in a water-only simulation that may be directly compared would be useful.

Also useful would be the .top + input .gro, e.g. so that if the issue is with the time step or the verlet-buffer-list setup we can see that.

Guillame's compiler was gcc 4.1.x, which could itself be a problem (unlikely, given the below).

Over 100ps with gcc 4.8 using a GPU, I observe with Guillame's .tpr

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Bond                        2630.98        120    294.859    708.723  (kJ/mol)
Angle                       3729.16         87    209.644    522.779  (kJ/mol)
Proper Dih.                 5378.84        6.1     47.195    24.3992  (kJ/mol)
Improper Dih.                260.84        2.4    20.2846    15.0637  (kJ/mol)
LJ-14                        1862.6        2.6     38.728    14.3807  (kJ/mol)
Coulomb-14                  16162.4         20    117.834    75.1824  (kJ/mol)
LJ (SR)                      147780        150    980.677    915.677  (kJ/mol)
Disper. corr.              -6674.06       -nan 5.73499e-06 -1.05932e-06  (kJ/mol)
Coulomb (SR)                -968839        460    1471.89   -2913.21  (kJ/mol)
Coul. recip.                4302.05        5.7    61.9105    -27.377  (kJ/mol)
Potential                   -793406         99    516.228   -664.376  (kJ/mol)
Kinetic En.                  128458         37    479.866    166.286  (kJ/mol)
Total Energy                -664949         70    144.078   -498.095  (kJ/mol)
Temperature                 295.202      0.086    1.10275    0.38212  (K)

This is comparable with Guillame's observed drift in total energy. Clearly something is wrong, but we can't tell what yet.

#2 Updated by guillaume chevrot about 6 years ago

Thanks for the test with another compiler.

I enclose the .top + .gro input files of the lysozyme simulation.

I can show excess drift in a water-only simulation of 1 ns with gromacs 4.6.1 (same drift with velocity Verlet or leap-frog integrator): see the tar file Water.tar.bz2.
Note: I only observe the drift in the case of single precision simulation (no drift in the case of double precision simulation).

#3 Updated by Mark Abraham almost 6 years ago

  • Category set to mdrun
  • Assignee changed from Michael Shirts to Mark Abraham

guillaume chevrot wrote:

I compared the total energy of 2 simulations:
lysozyme in water / NVE ensemble / single precision / Gromacs 4.6.3
lysozyme in water / NVE ensemble / double precision / Gromacs 4.6.3

I observe a constant drift in energy in the case of the single precision simulation (see file total_energy_lys.pdf)

Drift is normal and expected. The amount of drift acceptable for a simulation is an open question. There's little or no value in incurring expense to ensure the drift is tiny if the other approximations in the model are more significant - in particular, if sampling remains unconverged and the lack of convergence is partly caused by the cost of taking small integration steps.

I enclose the ‘mdout.mdp’ file so you can check the mdp options I used. Note: in this example, I used the velocity Verlet integrator. I also tested the leap-frog integrator and I observed the same drift in energy.

Guillame's simulation setup used 1fs time steps without constraints, which would not normally be considered a good idea. grompp even issues a note that the there is a protein bond vibration that has an expected period of 9fs. So who knows what conservation to expect? :)

I reran Guillame's 10ns 51k atom lysozyme sim, using pre-4.6.6 git version 53bfa7d under several configurations, the drift in Total Energy from 1ns to 10ns is:

Setup                                           Drift (kJ/mol)
dt=0.0005, constraints=none (ie. twice nsteps)  -67487.3 (which is -0.0589 k_B T / ns / atom)
dt=0.001,  constraints=none                     -40081.5
dt=0.001,  constraints=lincs, lincs-iter=1      -39620.8
dt=0.001,  constraints=lincs, lincs-iter=2      -40751.5

This seems consistent with past experience of SETTLE in single precision.

I also join the log, edr and tpr files.

Guillaume

PS: In the literature, this problem has already been mentioned: http://jcp.aip.org/resource/1/jcpsa6/v126/i4/p046101_s1

As discussed on gmx-users, that publication used GROMACS 3.3.1, and its observations are inapplicable to GROMACS 4.6.

I can show excess drift in a water-only simulation of 1 ns with gromacs 4.6.1 (same drift with velocity Verlet or leap-frog integrator): see the tar file Water.tar.bz2.
Note: I only observe the drift in the case of single precision simulation (no drift in the case of double precision simulation).

For example, the drift in Guillame's water.tar.bz2 with 511 waters is -153.708 kJ/mol over 1ns. Figure 1 of the paper he cites is for 901 molecules of water, so we can get -153 / 511 * 901 = -270 kJ/mol/ns. That's clearly better than it used to be. Guillame used lincs-iter=1, which is documented (http://manual.gromacs.org/online/mdp_opt.html#bond and in the PDF manual) to be suboptimal for NVE. I re-ran his simulation with lincs-iter=2 (with pre-4.6.6 git version 53bfa7d), and the drift over 1ns was reduced only slightly, to -145.904. Using SHAKE with tolerance 1e-6 was -148.192.

Note also that -153.708 is -0.0404 k_B T / ns / atom with a Verlet buffer of 0.2nm, which is quite consistent with (say) Figure 8 of http://dx.doi.org/10.1016/j.cpc.2013.06.003. So things seem to be working as intended. Whether that underlying drift is acceptable for observing things from a given simulation is up to the person doing it, but in many cases for other observables and other ensembles, I think the drift shown will prove acceptable.

#4 Updated by Berk Hess almost 6 years ago

Also note that with SETTLE in single precision the energy drift goes up with smaller time steps, as also noticed in the reference. With a 2 fs time step the drift will be smaller.
I think the current accuracy for single precision in Gromacs is enough for the vast majority of applications. Double precision can be used when this is not the case. This is not to say that we should try to improve the energy conservation is single (or better: mixed) precision whenever possible.

#5 Updated by Mark Abraham almost 6 years ago

NB https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2007-January/002028.html also states that the source of drift discussed in Guillame's paper reference has been resolved for a long time.

#6 Updated by Konrad Hinsen over 5 years ago

I think there is a bug hiding behind this discussion, though it's in the Gromacs manual rather than in the software.

As Mark says correctly, "Drift is normal and expected. The amount of drift acceptable for a simulation is an open question." In practice, that means that each scientist planning a simulation with Gromacs must decide where his/her priorities are. However, they can make an informed decision only if they are aware of the problem and understand its implications.

The single/double precision choice is discussed in the Gromacs manual, but only in the appendix "Technical Details". Moreover, the discussion there is highly technical and doesn't really address the impact on the final trajectories.

I don't agree with the Gromacs manual's atttitude that the tradeoff between stability/precision on one hand and speed/sampling quality on the other hand is a technical detail. It's something that anyone doing non-trivial simulations should be aware of, and a choice that any Gromacs user should report in the papers describing the simulations.

#7 Updated by Mark Abraham over 5 years ago

Konrad Hinsen wrote:

I think there is a bug hiding behind this discussion, though it's in the Gromacs manual rather than in the software.

As Mark says correctly, "Drift is normal and expected. The amount of drift acceptable for a simulation is an open question." In practice, that means that each scientist planning a simulation with Gromacs must decide where his/her priorities are. However, they can make an informed decision only if they are aware of the problem and understand its implications.

Indeed.

The single/double precision choice is discussed in the Gromacs manual, but only in the appendix "Technical Details". Moreover, the discussion there is highly technical and doesn't really address the impact on the final trajectories.

I don't agree with the Gromacs manual's atttitude that the tradeoff between stability/precision on one hand and speed/sampling quality on the other hand is a technical detail. It's something that anyone doing non-trivial simulations should be aware of, and a choice that any Gromacs user should report in the papers describing the simulations.

From 5.0, we plan to refer to the two implementations as mixed and double precision, since the most critical parts (of virial computation) are always done in double precision. But spending some quality time with the documentation and sed is still to be done ;-) I plan also to introduce a GROMACS "Usage Guide" to address a number of these kinds of practical issues in a setting distinct from the reference manual, and I will incorporate your suggestions there. Thanks!

#8 Updated by Rossen Apostolov over 5 years ago

  • Status changed from Feedback wanted to Closed

#9 Updated by Mark Abraham over 5 years ago

  • Category changed from mdrun to documentation
  • Status changed from Closed to In Progress
  • Target version set to 5.0

I still want to fix the documentation aspect

#10 Updated by Rossen Apostolov over 5 years ago

  • Related to Task #1063: update documentation to correct "single" precision added

#11 Updated by Erik Lindahl over 5 years ago

  • Status changed from In Progress to Fix uploaded

Documentation fix is at https://gerrit.gromacs.org/#/c/3719/

#12 Updated by Erik Lindahl over 5 years ago

  • Status changed from Fix uploaded to Resolved

#13 Updated by Erik Lindahl over 5 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF