4 fs timestep vibrational frequency integration note - charmm36m forcefield

User discussions
,
1

GROMACS version: 2020.1
GROMACS modification: No

Hi,

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.

refs:

  1. Removing fastest degrees of freedom — GROMACS 2021.2 documentation
  2. Bond oscillations versus timestep
  3. Miscellaneous — GROMACS 2019.2 documentation

Can anyone provide any recommendations to continue? Thank you.

2

I would also appreciate any references that anyone can point me to that discuss this numerical error or why a minimum value of 10 iterations is chosen by grompp but 5 by the documentation? Is there some quantification of the % error somewhere?

3

Hi!

Why would you want to use 4 fs timestep?

This is a timestep of a CG system! Constraing all bonds would get you only to 2 fs. May be angles could get you a bit further, but why?

Just as gromacs advises you - your vibrations are too fast for such a big timestep. Don’t do it!

V

4

This approach is quite common, especially in AMBER force fields, via virtualization of hydrogens or hydrogen mass repartitioning. It’s not common with CHARMM, perhaps for this reason, but it can be done and studies have been published using this approach.

This isn’t true; normal usage of the CHARMM force field constrains only bonds to H and still uses a 2-fs time step, and bonds that are not designed to be constrained should not be as this is simply incorrect use of the force field. Angles should generally never be constrained in common biomolecular force fields. That approach generally isn’t numerically stable.

1 Like