GROMACS version: 2020.1
GROMACS modification: No
I am trying to use virtual sites for my protein in the CHARMM36m forcefield to achieve a 4 fs timestep with constraints = h_bonds. According to the Gromacs documentation (ref 1 below):
“Taking as a guideline that with a Verlet (leap-frog) integration scheme a minimum of 5 numerical integration steps should be performed per period of a harmonic oscillation in order to integrate it with reasonable accuracy,”
However, when I run my system through grompp I get the following note:
NOTE 1 [file topol.top, line 49]:
The bond in molecule-type Protein_chain_A between atoms 290 CD and 291
OE1 has an estimated oscillational period of 2.2e-02 ps, which is less
than 10 times the time step of 4.0e-03 ps.
Maybe you forgot to change the constraints mdp option.
My question is, am I safe to follow the guideline that states 5 numerical integration steps per period is reasonably accurate (as is the case in my system), or should I consider changing all bonds to constraints? I am somewhat reluctant to use constraints = all_bonds based on comments on previous posts (ref 2):
“If the force field was designed to only have bonds to H constrained, that’s how it should be used. If it was designed for all bonds to be constrained, that’s how it should be used. This is important because if you change this setting, you’re introducing preventable errors in the simulation.”
Additionally, I have seen (ref 3 and my own tests) that constraints = all_bonds results in a loss in accuracy in addition to an increase in computational expense.
- Removing fastest degrees of freedom — GROMACS 2021.2 documentation
- Bond oscillations versus timestep
- Miscellaneous — GROMACS 2019.2 documentation
Can anyone provide any recommendations to continue? Thank you.