Simulations with time step 4 fs using constraints=H-angles

User discussions
1

GROMACS version: Release 2024.1
GROMACS modification: No

Dear Gromacs community

I am investigating how to carry out simulations using a time step equal to 4 fs. In the manual (whose publication date is Feb 28, 2024), I have read that this can be done using two different methods:

A) Applying mass repartitioning [1].
B) Applying virtual sites plus mass repartitioning [2].

The literature suggests that imposing constraints on the angles of hydrogen atoms makes it possible to increase the time step of the simulation [3, 4]. Imposing constraints on angles of hydrogen atoms is supported by GROMACS setting <<constraints=h-angles>>. However, I cannot see this approach mentioned in the documentation as a recommended way to achieve a 4 fs time step. Is there any reason why mass repartitioning/virtual sites are preferred over constraining hydrogen angles?

Thank you very much. With kind regards

[1] From the manual: mass-repartitioning-factor (…) With h-bonds constrained a factor of 3 will usually enable a time step of 4 fs.

[2] From the manual: Disregarding these very fast oscillations of period 13 fs, the next shortest periods are around 20 fs, which will allow a maximum time step of about 4 fs. Removing the bond-angle degrees of freedom from hydrogen atoms can best be done by defining them as virtual interaction sites instead of normal atoms. Whereas a normal atom is connected to the molecule with bonds, angles and dihedrals, a virtual site’s position is calculated from the position of three nearby heavy atoms in a predefined manner (see also sec. Virtual interaction sites (page 435)).For the hydrogens in water and in hydroxyl, sulfhydryl, or amine groups, no degrees of freedom can be removed, because rotational freedom should be preserved. The only other option available to slow down these motions is to increase the mass of the hydrogen atoms at the expense of the mass of the connected heavy atom.

[3] Mazur, A. Hierarchy of Fast Motions in Protein Dynamics. J. Phys. Chem. B 1998, 102, 473–479.

[4] Stocker, U.; Juchli, D.; van Gunsteren, W. F. Increasing the time step and efficiency of molecular dynamics simulations: optimal solutions for equilibrium simulations or structure refinement of large biomolecules. Molecular Simulation 2003, 29, 123–138.

2

You should not constrain all H-angles as that will adversely affect the conformational distribution.

3

Dear Prof. Hess
Thank you very much for your reply. If I am not mistaken, <<constraints=H-angles>> just fixes H-X-H and X-O-H angles, where X is a heavy (non-H) atomic species. May you please explain a bit more why such constraints distort the dynamics more than mass repartitioning or virtual sites?
Thanks again!

4

I thought it constrained all angles involving hydrogens, but I now see that you are right. Where did you find this information?

I don’t know what the vibrational frequency of an angle involving a single hydrogen is. Maybe it’s just ok for a 4 fs time step. But constraining hydrogen angles requires constraining all bonds, which is generally not advisable as this changes the energy landscape. Also constraining will be expensive.

5

Futhermore, the h-angles option in GROMACS also constrains all bonds. This is not advisable with most force fields, as they have been parametrized with flexible bonds between heavy atoms and constraining these can increase energy barriers of dihedrals.

6

Dear Professor Hess
Thank you very much for your kind reply. You asked <<Where did you find this information?>>. Actually we did not find it in any document, but by carrying out some calculations of proteins.
We are also interested in knowing how many parameters in a force field depend on H-angles for a given H atom (this is, for one given H atom, how many H-angles have a contribution to the energy). May you give us a clue of how we may answer this question?
Thank you very much again.

7

All atomistic force field I know of have angle potentials between all atoms connected by two bonds. Then there are dihedral potentials. Some force field use dihedral potentials between all quadruplets connected by three bonds in a line. Some force fields exclude quadruplets with hydrogens on both ends. And the there can be improper dihedral potentials.