Energy drift in NVE

GROMACS version: GROMACS/2021.3
GROMACS modification: No

I’ve encountered an issue with NVE for a system consisting of a single type of organic molecule (C24H10N2O4), which seems to be well equilibrated before running NVE. See the included figures for plots of the total system energy, pressure and temperature as a function of simulation length, where the different colored curves corresponds to the different parts of my simulation protocol. The system looks well equilibrated before NVE in my opinion, but as soon as I remove the thermostat/barostat the energy, pressure and temperature drifts.

I’ve tried decreasing the timestep in the NVE simulation to 0.05 fs, as well as run a longer NVE simulation without any luck. I also added an extra round of NVT (nvt_intermediate) before NVE as I read in another post on this forum that that would be beneficial.

Does anyone know what might be causing this behavior? Any tips or feedback on what I might try next is greatly appreciated!

Here is my simulation protocol in more detail:
Timestep: 0.5fs
ref_t: 400K

  1. em (Energy minimization): steep
  2. nvt: tcoupl=v-rescale, tau-t=1, pcoupl=no, generate velocities at 400K
  3. npt_high_p: tcoupl=berendsen, tau-t=2.0, pcoupl=berendsen, tau-p=1, ref-p: 2000 bar (to fix cavitation in the system)
  4. npt: tcoupl=v-rescale, tau-t=2.0, pcoupl=Parinello-Rahman, tau-p=5, ref-p: 1 bar
  5. nvt_intermediate: v-rescale, tau-t=1, pcoupl=no
  6. nve_prod: tcoupl=no, pcoupl=no

Plot of energy during simulation

Plot of pressure during simulation

Plot of temperature during simulation

there are numerous things in an MD simulation that can contribute to the energy drift. Maybe the question in your case is from where the largest contribution to the drift comes. I would check the following:

  • to start with, how does your observed drift compare to published drift settings, e.g. Fig. 8 in “A flexible algorithm for calculating pair interactions in SIMD architectures” (Pall, Hess) 2013. Does extending the Verlet buffer size reduce your drift?
  • if you constrain bonds, the LINCS settings can be refined for higher energy conservation, see .mdp section in the user guide
  • using double precision can help, whereas a too small time step might even be counter-productive
  • tweaking nonbonded parameters (vdW/PME) might help