Cutoff-scheme and intermolecular interactions

User discussions
,
1

GROMACS version: 2020 and above
GROMACS modification: No

Hi everyone,
thank you for reaching out. I am writing to deepen into two main topics.

  1. I would like to ask whether the “group” cutoff-scheme will be will be reimplemented, as it is around years (I think) it is not available anymore. In particular, I think that would be absolutely interesting in decoupling a solute from the environment via a lambda parameter.

Here I provide an example. At the very moment, if going for thermodynamic integration to compute free energy of solvation, either I choose to put couple-intramol = no, for which case I noticed I have to set rlist to a value which is larger than the maximum distance between two atoms of the solute, along with verlet-buffer-tolerance = -1 so to keep that value, or otherwise I face this very problem:

“Fatal error:
There are perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.2 nm, which is not supported. This can happen because the system is
unstable or because intra-molecular interactions at long distances are
excluded. If the latter is the case, you can try to increase nstlist or rlist
to avoid this.The error is likely triggered by the use of couple-intramol=no
and the maximal distance in the decoupled molecule exceeding rlist.”

I really made sure the system is equilibrated, and in fact my error is simply removed if I respect precisely what I stated. Otherwise, an alternative to switch off this very error is putting couple-intramol = yes, and rlist can be finally set to the value desired (which is de facto rcoulomb, if using PME), and also verlet-buffer-tolerance is free to be set again to the default value. But as you know, in that case one should redo another thermodynamic integration procedure on the solute in vacuum.

If using gromacs 2019, one can (I verified) switch off the error with this very possibility: using the finally available cutoff-scheme = group - instead of verlet -, and putting energygrps = solute along with energygrp_excl = solute solute . This very alternative is by far the preferred one if the solute is kept frozen for some reasons/approximations one is willing to do (see manual of 24.02 version, last sentence in the description of energygrp_excl).

  1. With regard to the very first alternative I stated, I would like to ask whether there is simply a slow down in the trajectory propagation, or if putting rlist to values like 2.5 nm and fixing it via verlet-buffer-tolerance = -1 also has a negative physical reverberation on the behavior of the system. Specifically, I ask this for cases in which rlist is fixed to values larger than the typical 1.0 and above one usually finds: should I be worried of this first alternative, in this very case?

Thanking you for your time, I wish you a lovely day.

Bests,
Jacopo

2

The energygroup exclusion solution is not a solution. You would run into exactly the same exclusion issue, but maybe there was no check for this in GROMACS 2019.

You can increase rlist, this only has effects on performance, not on the results. But I would choose to use allParameters.mdp couple-intramol = yes. Larger molecules tend to collapse when decoupled from solvent which can maks sampling (much) more difficult.