Ask.Cyberinfrastructure

Implementing full geometry constraints in GROMACS simulations

Does anyone have experience with molecular dynamics simulations in GROMACS with geometry constraints? I know that GROMACS implements LINCS, SHAKE, and SETTLE, in order to constrain chemical bond lengths, angles, etc., and each of these has different capabilities/advantages. That being said, I rarely run constrained simulations.

I am trying to constrain all bond lengths and angles, meanwhile using my own table files (custom force field), and hopefully run on more than 1 CPU.

I run into trouble because the use of custom table files forces me to alter the “charge group” of atoms in molecules, which causes SHAKE to fail on >1 CPU. But LINCS can not constrain angles. My current workaround has been to make all bond angles very stiff by using arbitrarily large angle force constants. This allows the “rigid” simulation to run on multiple CPUs using LINCS, but I won’t generate proper energetics any more.

Hi Ryan,

I haven’t worked with GROMACS before, so my assistance may be limited. Which version are you using and how are you using parallelization? Is this document helpful? http://www.gromacs.org/Documentation/Acceleration_and_parallelization

If you are using GROMACS 4.6, you might be able to try using MPI and/or OpenMP. Apologies if you are already doing this and still running into this problem.

I am using GROMACS 2018.4. I am indeed already using parallelization with MPI. I don’t think it is the MPI that is the problem, per se.
The problem is caused by the fact that the algorithm needed for constraining bond lengths and angles (SHAKE) does not seem to work on more than one CPU, if custom table files are used. There is a constraint algorithm that works find with multiple CPUs and custom table files (LINCS), but this algorithm can not constrain bond angles, only lengths.

As far as I understand, Gromacs can constraint both bonds and angles with LINCS using the constraints = all-angles option (http://manual.gromacs.org/documentation/2018-current/user-guide/mdp-options.html#bonds ). As stated in the manual:

Convert all bonds to constraints and all angles to bond-constraints.

The shake algorithm is only included for reference in Gromacs, and therefore is only used for testing.