# Energy drift in NVE **URL:** \ **Category:** User discussions\ **Tags:** mdrun\ **Created:** [July 6, 2022, 1:45pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423 "2022-07-06T13:45:37Z")\ **Posts on this page:** 14\ **Page:** 1
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [July 6, 2022, 1:45pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/1 "2022-07-06T13:45:37Z")
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](https://gromacs.bioexcel.eu/t/how-to-use-nve-on-simulations/2208/9) 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 ![energy](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/d/d59465ce2c5244ce28f25782142587698b2e625e.png) Plot of pressure during simulation ![pressure](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/a/af600daf8f2f459706ab5b727747dddbcd181a41.png) Plot of temperature during simulation ![temperature](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/d/d128e1b7fc36b136994a08839743177fe32aefc3.png) ---
**Author:** ![Carsten](https://avatars.discourse-cdn.com/v4/letter/c/f14d63/32.png) [@Carsten](https://gromacs.bioexcel.eu/u/Carsten)\ **Post date:** [July 6, 2022, 3:10pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/2 "2022-07-06T15:10:00Z")
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 ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 24, 2022, 11:31am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/3 "2022-08-24T11:31:11Z")
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: ```bash ; 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: ```bash 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: ```bash 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 ![nve_energy](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/9b985b219dff11134ae903c40dffd1ce005607f3.png) e Temperature ![nve_temp](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/97baef26dd6c1f114623f8ce976de813ac9ed404.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 24, 2022, 12:40pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/4 "2022-08-24T12:40:39Z")
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. ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 24, 2022, 1:28pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/5 "2022-08-24T13:28:34Z")
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: ```bash ; ; 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](http://atb.uq.edu.au/molecule.py?molid=23979#panel-information). 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 ```bash constraints = all-bonds constraint-algorithm = LINCS continuation = no ``` and run the simulation? ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 24, 2022, 1:42pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/6 "2022-08-24T13:42:56Z")
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. ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 24, 2022, 1:45pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/7 "2022-08-24T13:45:00Z")
Alright, I’ll try that! ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 29, 2022, 3:23pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/8 "2022-08-29T15:23:25Z")
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 ![energy_drift](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/9141b4a161814b1490914625107a33cbfa8487af.png) Temperature during NVE simulation ![temperature](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/8/8112d3d474c466e3addacfe26cd3775d7b954a0c.png) Temperature jump when going from NVT to NVE ![nvt_to_nve](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/3/34a4990cbe0aae770dcb3a5763dd6f8e1f681b89.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 30, 2022, 7:27am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/9 "2022-08-30T07:27:06Z")
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? ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 30, 2022, 8:40am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/10 "2022-08-30T08:40:15Z")
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: ```bash ; Constraints constraints = h-bonds constraint-algorithm = LINCS continuation = no lincs-iter = 2 ; 1 is sufficient, but for NVE use 2 ``` ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 30, 2022, 8:47am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/11 "2022-08-30T08:47:04Z")
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? ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 30, 2022, 9:01am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/12 "2022-08-30T09:01:43Z")
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.) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [August 30, 2022, 10:14am UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/13 "2022-08-30T10:14:26Z")
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. ---
**Author:** ![elindgren](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@elindgren](https://gromacs.bioexcel.eu/u/elindgren)\ **Post date:** [August 30, 2022, 2:24pm UTC](https://gromacs.bioexcel.eu/t/energy-drift-in-nve/4423/14 "2022-08-30T14:24:13Z")
Hm I see. I’ll read up on how much NVE/NVT affects the result in practice then. Thank you for your help!