No domain decomposition while using restraints

GROMACS version:version 2021.3-MODIFIED
GROMACS modification: Yes

Hi, I have a small system with 87 residues of a protein. And I want to maintain the tertiary structure of the protein, so I came up with some restraints which I am supplying in the topology file under the [bonds] directive. It is as follows

   137  1188     6    3.253    8000
   271  1125     6    3.043    8000
    82  1004     6    1.386    8000
   319   979     6    2.600    8000
    82   923     6    1.309    8000
     7   764     6    2.059    8000
   257   662     6    2.415    8000
   319   829     6    1.809    8000
   151   371     6    1.859    8000
   371   500     6    1.752    8000
   462   764     6    1.440    8000
   923  1049     6    1.770    8000
  1004  1188     6    2.215    8000

Initially, I have a water padding of 1 nm (total number of atoms ~ 60k) and without the above restraints the simulation runs fine, but when I add the restraints, I got the following error. So, I am pretty sure that my type 6 restraints are causing the problem.

Fatal error:
There is no domain decomposition for 96 ranks that is compatible with the
given box and a minimum cell size of 4.21705 nm
Change the number of ranks or mdrun option -rdd or -dds
Look in the log file for details on the domain decomposition

I have searched the mailing list and found that if you increase the water padding (or make the box larger) you can avoid this error. So I tried going for 60-angstrom padding, but still have the same error and my total number of atoms reach close to 700k.

Is there a workaround for this problem? I do not want to simulate a bigger system because of my limited computational resources?

I am attaching my log file for the NVT run.

nvt.log (15.5 KB)

I found a similar post in the mailing archive, a solution to this problem.
Here is the post:

A real fix would then be adding
i!=F_HARMONIC &&
to the if statement on line 272 of src/gmxlib/mshift.c
and recompiling Gromacs, such that the harmonic potentials
are not put into the graph.

But I couldn’t find this line In the current GROMACS version (2021.3) mshift.cpp, As this answer was written 15 years ago, I assume the code changed slightly. Could anyone help regarding where should I change the code and recompile it?

The code has changed massively in 15 years, and the fix proposed above may even pre-date domain decomposition (which I believe was released in 2008).

Hello, I was wondering if you finally managed to overcome this problem somehow, @vasi786

Sincerely,

Agustin

@agustin.zavala, Hello. I did not solve the problem but found a work around it. Instead of doing restraints using [bonds] directive, I did pulling with rate of pulling set to zero. I checked by running a simulation, and the tertiary structure of the protein was maintained during the run.
See the code below

pull-group1-name     = R1
pull-group2-name     = R2
pull-coord1-groups   = 1 2
pull-coord1-type     = umbrella
pull-coord1-geometry = distance
pull-coord1-dim      = Y Y Y
pull-coord1-k        = 3000
pull-coord1-init     = 1.2130461160248225
pull-coord1-rate     = 0

Thank you @vasi786 .
I thought about doing the same. Except I want to use an upper and lower limit… so a Piecewise linear/harmonic restraint, or, if I want to use the pull code, pull-coord1-type = flat-bottom and pull-coord2-type = flat-bottom-high. But so far I haven’t been able to, when I set that type of restraint gromacs complains about the decomposition of the system, the length of these restraint “bonds” is much larger than the portions the simulation is decomposed in, and if I use the pull code, for some reason gromacs complains that the energy of the system is non-finite and stops… and I’m not sure why that is, maybe I’m setting up the pull code incorrectly

Hello,

update to this issue, I finally got help to install gromacs 2023 interlaced with PLUMED (Max Bonomi - PLUMED - Research - Institut Pasteur) on out cluster and used this software to impose upper and lower wall limits for a set of distances between different pairs of atoms in my system, as they did here: Karami et al, Computational and biochemical analysis of type IV pilus dynamics and stability 10.1016/j.str.2021.07.008

It was quite simple to use, efficient, and to install it I got help from the developpers, right away. I understand PLUMED also offers many other features as well.

But how is the parallel performance then? As PLUMED needs to gather all data to a single rank, performance is going to be severely affected when running on 96 MPI ranks.

When fewer ranks are needs, the GROMACS restrains or pull code will work fine.

Well, I can say that I get between 30 and 45 ns/day for a ~300000 particle system, running 1 MPI process with 24 OpenMP threads on 24 cores and 1 GPU, using GROMACS 2022.5-plumed_2.9.0. I think it varies because GPU could be A40, A100 or V100 (depending on what’s available on our cluster), and also I think the Performance is different when I have multiple simulations running on the same node. These conditions were chosen based on the resources we had available most of the time in our cluster, and our needs, 3 systems, ~300.000 particles each, for 1000ns, by triplicate.

I can’t find the Performance for this system without the PLUMED restraints, I think perhaps at some point I reached 60-70ns/day, but I could be mistaken, as I’ve rebuilt and modified our system several times.

But I imagine you’re right, yes, if one could use the GROMACS restraints or pull code, they should work faster, I suppose, and be simpler than installing PLUMED alongside GROMACS.

But as I explained above, I just couldn’t get them to work as I needed… it’s possible I kept making a certain mistake all the time, but I tried many ways to make either of those two strategies work to maintain upper and lower wall distance restraints between 4 proteins in my system (between atoms not quite close to each other) and just couldn’t manage…

PS: and yes, I’m using a single MPI rank per simulation, we didn’t count with dedicated resources for our project and thought of perhaps running even a larger number of systems, so we aimed for throughput more than Performance.

But if you would use 1 MPI rank also with the pull code there cannot be issues with the domain decomposition as there is no domain decomposition.

The problem with dd was for the distance restraints or bond restraints, not for the pull code, if I remember correctly.
But yeah, perhaps having changed our strategy from using several MPI to 1 MPI per simulation (which happened after we solved this issue by installing gromacs with plumed), we should try again the pull code, or the built-in gromacs restraints, and if there’s no issue then I would expect them to be faster than using gromacs with plumed. I should try this when I get a chance.

But my original plan was to use several MPI ranks, I couldn’t make the pull code or the restraints work, so I found this other option and I posted it here in case anybody stumbled upon the same issue as me when trying to use multiple MPIs and these types of forces.

It’s also possible the issue happened because we were using an older Gromacs version before, and for the new Gromacs+Plumed installation we upgraded to a newer version… I’m not sure…