GROMACS version:2018.3
GROMACS modification: No
Hello
I have plotted the energy for one of my systems (solution of 9 protein chains, 168 aminoacids each, waters, 150mM NaCl in Martini3 FF, simulated with Gromacs 2018.3). The conserved energy, which supposed to be total energy exchange between the barostat and the system, linear drift due to num.rounding errors and approximations has huge drifts which apparently depend on the length of simulated interval (first 5 intervals were 3.5ns and last one 8ns). Can anyone advise why I’m generating these drifts and what might be wrong with my protocol (link to files attached)? Thanks.
You files are not accessible.
But some linear drift is normal and usually not problematic. There are many sources of drift. One is the pair-list buffering. This drift is controlled by the verlet-buffer-tolerance parameter. The default value is 0.005 kJ/mol/ps per atom. So for a 200 000 atom system you get max 1000 kJ/mol drift in a ps, or a million kJ/mol in a ns. This might seem like a lot, but per atom per ps this is negligble and does not cause artifacts.
Hello, guys.
I’ll post my question here because I’m experiencing some huge drift in the Conserved Energy. Even if my system is not under the Martini force field, I think that our problems have the same solution through improving lincs-order and lincs-iter parameters.
First, I detail my system (~50k atoms): protein+small ligand in water using ff19SB/GAFF2/OPC3 with hydrogen mass repartitioning (HMR) built with AmberTool25 and converted to the GROMACS format (CMAP per amino acid are correct, just like as in CHARMM-GUI). I performe a 1us per replica.
In my first tests (dt = 4 fs), the system was exploding or presenting the cuda error #700 after 100-300 ns. I was using 4 and 1 to lincs-order and lincs-iter, respectively, PR barostat with tau-p 5 ps, and v-rescale thermostat with tau-t 1 ps. In such a system, the drift was around 1.69e-04 kJ/mol/ps/atom.
In a similar system (with only the protein) without HMR (and using dt = 2 fs), the conserved energy drift was really small, resulting in a drift of 1.16e-05 kJ/mol/ps/atom.
It seems that I solved the crash problem by using lincs-order and lincs-iter to 6 and 2, respectively, as suggested by GROMACS (interestingly, the Martini force field recommends 8 and 2 for these parameters). However, I was wondering if I could obtain a smoother drift as observed in the system with 2 fs. So I found that using the C-rescale barostat with tau-p of 6 ps presented a drift of 5.83e-05 kJ/mol/ps/atom and a value of tau-p of 8 ps presented a drift of 6.11e-06 kJ/mol/ps/atom.
Would it be better to use tau-p equals 8 instead of 6 (maintaining lincs-order=6 and lincs-iter=2) just to obtain a reduced drift in the simulation?
And would a smaller change in the conserved energy generate more precise/stable simulations?
Best,
I think you should use much larger tau-p. I see no reason for having it this small. Something like 20 ps.
Low energy drift is not a goal per se. Also note that there drift contributions that can have opposite sign, so smaller total drift is not necessarily better. 1e-4 kJ/mol/ps/atom is very small. There is no reason for pushing for even smaller values in typical simulations, unless you want to run NVE.