Energy drift in NVE

[elindgren]

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

[Carsten]

Hi,
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
  Best,
  Carsten

[elindgren]

Thank you for your detailed reply! I apologize for not replying very quickly. I’ve been trying some of your suggestions.

Comparison of drift when increasing verlet buffer

I made the buffer larger (I think) by setting the following cutoff scheme:

; Neighborlist
cutoff-scheme   = Verlet
nstlist         = 40 ; Frequency to update the neighborlist
verlet-buffer-tolerance     = 0.0005 ; Decrease by a factor of 10

Gromacs reports the following drift in the total energy:

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -25475.3      22000    45867.5     158876  (kJ/mol)

Statistics over 2000001 steps [ 0.0000 through 1000.0000 ps ], 1 data sets
All statistics are over 20001 points

This is for 1000ps, with 59880 atoms in the simulation box. Converting this to drift per atom in kBT/ns (@400K) I obtain a drift of 4*10^4 kBT/ns per atom, which is at least four orders of magnitude higher than the ones in Fig. 8 in the paper you mentioned. The drift for my other runs (e.g. npt equilibration) are much more in line with the numbers in the paper though.

Heuristically it does seem that modifying the verlet buffer improved the drift a little bit though; comparing the figures below with the ones in my previous post the total energy drifts slightly less.

Are my settings for the verlet buffer reasonable? Can I change them further to hopefully decrease the drift? Given the effect they’ve had though it seems unlikely that this is the only source of the drift.

Lincs and constraining bonds

I’m not constraining any bonds explicitly, so I guess changing the LINCS settings will not help me?

Double precision Gromacs

I unfortunately don’t have a double-precision version of Gromacs compiled on my HPC clusters. If nothing else works I can try to gain access to one.

Nonbonded parameters, vdW/PME

I haven’t been able to run with PME, since I get a warning from Gromacs when trying to enable it. Could this perhaps be the source of the issue? See below:

NOTE 1 [file topol.top, line 11]:
  System has non-zero total charge: 0.149549
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.

I don’t see how my system could be charged, since I have neutral perylene molecules and no extra charges. I’ve always gotten this warning when using these .gro-files though, so perhaps it could be something up with them? I’ve generated the .gro-files from .pdb-files using the LigParGen webserver, which parametrize them to OPLS compatible files.

Should I change force field, or try to add charges to neutralize the system?

Energy and temperature when running with increased Verlet buffer

Total energy

e

Temperature

[hess]

You surely want to use PME.
That you system doesn’t have integer charge indicates that there likely is an issue with you topology.

Why are you not using constraints on bonds involving hydrogens? With most force fields you have to use such constraints for correctness.
If you have forcefield that is parametrized without constraints on hydrogen atoms, you will likely need to compiler GROMACS in double precision to get good energy conservation with the smaller time step required.

[elindgren]

Thank you for your reply!

You surely want to use PME.
That you system doesn’t have integer charge indicates that there likely is an issue with you topology.

Alright! I checked the charges, and there seems to be a slight charge associated with each molecule of 0.0001, since the partial charges doesn’t add up. Here is part of my .itp-file:

;
; GENERATED BY LigParGen Server
; Jorgensen Lab @ Yale University
;
...
MOL                   3
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass
     1   opls_800      1    MOL   H00      1     0.5273     1.0080
     2   opls_801      1    MOL   N01      1    -1.0719    14.0070
     3   opls_802      1    MOL   C02      1     0.6600    12.0110
     4   opls_803      1    MOL   O03      1    -0.4108    15.9990
     5   opls_804      1    MOL   C04      1    -0.1503    12.0110
     6   opls_805      1    MOL   C05      1    -0.0462    12.0110
     7   opls_806      1    MOL   H06      1     0.1918     1.0080
     8   opls_807      1    MOL   C07      1    -0.1448    12.0110
     9   opls_808      1    MOL   H08      1     0.1616     1.0080
    10   opls_809      1    MOL   C09      1     0.0127    12.0110
    11   opls_810      1    MOL   C0A      1     0.0129    12.0110
    12   opls_811      1    MOL   C0B      1    -0.1450    12.0110
    13   opls_812      1    MOL   H0C      1     0.1617     1.0080
    14   opls_813      1    MOL   C0D      1    -0.0458    12.0110
    15   opls_814      1    MOL   H0E      1     0.1918     1.0080
    16   opls_815      1    MOL   C0F      1    -0.1499    12.0110
    17   opls_816      1    MOL   C0G      1     0.6595    12.0110
    18   opls_817      1    MOL   O0H      1    -0.4107    15.9990
    19   opls_818      1    MOL   N0I      1    -1.0693    14.0070
    20   opls_819      1    MOL   H0J      1     0.5275     1.0080
    21   opls_820      1    MOL   C0K      1     0.6596    12.0110
    22   opls_821      1    MOL   O0M      1    -0.4109    15.9990
    23   opls_822      1    MOL   C0N      1    -0.1485    12.0110
    24   opls_823      1    MOL   C0O      1    -0.0484    12.0110
    25   opls_824      1    MOL   H0P      1     0.1918     1.0080
    26   opls_825      1    MOL   C0Q      1    -0.1445    12.0110
    27   opls_826      1    MOL   H0R      1     0.1617     1.0080
    28   opls_827      1    MOL   C0S      1     0.0139    12.0110
    29   opls_828      1    MOL   C0T      1    -0.0376    12.0110
    30   opls_829      1    MOL   C0U      1     0.0310    12.0110
    31   opls_830      1    MOL   C0V      1     0.0134    12.0110
    32   opls_831      1    MOL   C0W      1    -0.1452    12.0110
    33   opls_832      1    MOL   H0X      2     0.1616     1.0080
    34   opls_833      1    MOL   C0Y      2    -0.0491    12.0110
    35   opls_834      1    MOL   H0Z      2     0.1918     1.0080
    36   opls_835      1    MOL   C10      2    -0.1489    12.0110
    37   opls_836      1    MOL   C11      2     0.6611    12.0110
    38   opls_837      1    MOL   O12      2    -0.4100    15.9990
    39   opls_838      1    MOL   C13      2     0.0326    12.0110
    40   opls_839      1    MOL   C14      2    -0.0374    12.0110

I’m using the LigParGen server to generate .itp and .gro-files like this from .pdb-files. I get the pdb-files from the ATB, for instance this one.
I guess removing/changing the partial charge by modifying the .itp-file is not a very smart option. Would relaxing the structures myself via DFT be an option? Or is the topology just plain weird?

Why are you not using constraints on bonds involving hydrogens? With most force fields you have to use such constraints for correctness.
If you have forcefield that is parametrized without constraints on hydrogen atoms, you will likely need to compiler GROMACS in double precision to get good energy conservation with the smaller time step required.

I’m new to Gromacs, so I didn’t know that you had to constrain bonds yourself, I thought that was handled by the force field. How do I know if my force field is bonded or not? I’ve looked around but I haven’t been able to discern if OPLS is bonded or not.

Assuming OPLS is bonded, if I want to use bond constraints, would I just add

constraints = all-bonds
constraint-algorithm = LINCS
continuation = no

and run the simulation?

[hess]

Just change the charge of one atom by 0.0001 to neutralize the molecule.

You should use constraints=h-bonds.
Then you can use a time step of 2 fs and get reasonable energy conservation.

[elindgren]

Alright, I’ll try that!

[elindgren]

Using PME works much better! The drift isn’t nearly as bad, but there is still a marked slope to the total energy and temperature when I run NVE (see figures below).

Furthermore, the temperature also jumps up by about 15K when I switch from NVT to NVE (see the third figure). Is this common, and I just should run my simulation for longer to get it to stabilize?

Energy during NVE simulation

Temperature during NVE simulation

Temperature jump when going from NVT to NVE

[hess]

The will always be a slope due to numerical precision in all algorithms involved. Often verlet-buffer-tolerance gives the largest drift. You can lower that and see if your energy conservation improves. But a negative drift is usually caused by the constraint algorithm. I assume you are using LINCS with default settings? You can also try double precision.

Is the temperature at step 0 in the NVT ensemble not also 15 K lower than the target temperature?

[elindgren]

Alright, I see. I’m currently testing lowering the verlet-buffer-tolerance to see if that makes an improvement. I checked the log-file from the nve-run and Gromacs recommended me to run with a buffer-tolerance of ~6e5, so I’m currently trying 5e-5.

I’m using LINCS with mostly default settings (I think), but I increased the lincs-iter to 2 in accordance with what I read on in the .mdp-docs.

Here are my constraint settings:

; Constraints
constraints             = h-bonds
constraint-algorithm    = LINCS
continuation            = no
lincs-iter              = 2 ; 1 is sufficient, but for NVE use 2

[hess]

You mean 6e-5? 5e-5 is not a significant difference.

My guess is that your drift comes from the constraints and this can’t be improved more in single precision.

Why do you want to use and NVE ensemble?

[elindgren]

I was unclear. I was initially using the default buffer-tolerance of 5e-3, which Gromacs recommended me to lower to 6e-5. To be on the safe side, I’ve lowered it to 5e-5.

My guess is that your drift comes from the constraints and this can’t be improved more in single precision.

Hm I see, I’ll consider double precision then.

Why do you want to use and NVE ensemble?

From what I’ve understood NVE is the “gold standard” when it comes to MD simulations, since there is no barostat/thermostat. I’m studying the structural and dynamical properties of a system of molecules, where I’m particularly interested in things such as dynamical structure factors, density of states etc. I thought the thermostat/barostat would affect the dynamics, and thus I want to run without them. But maybe I’m mistaken?

For example, one issue I’ve been having is that the gmx_dos module seems to give me an offset in the DOS, as well as it not being properly normalized to 3N, even though I’m running with a very short dump frequency (1fs). My guess was that this might have something to do with the thermostat/barostat, and hence I wanted to try it. (It’s possible that this offset is yet another artifact of not running with PME as we discussed earlier, but I’ve yet had time to try this.)

[hess]

Some physicists think anything but NVE is wrong when looking at dynamics, which theoretically might be correct. But there are many approximations and other factors that affect MD results that the effect of NVT versus NVE is completely negligible in nearly all applications.

Just use NVT, or even NPT.

[elindgren]

Hm I see. I’ll read up on how much NVE/NVT affects the result in practice then. Thank you for your help!