[gmx-users] Periodic molecules and minimum bonded cut-off

[Mark Abraham]

On 8/03/2012 12:55 PM, Andrew DeYoung wrote:
> Hi,
>
> I am running an NVT simulation involving a periodic molecule, which has
> bonds that cross the boundaries of the simulation box.  So, I am using
> periodic_molecules = yes in my .mdp file.  I am also using the following
> cut-offs in my .mdp file:
>
> vdwtype = cut-off
> coulombtype = PME
> rlist = 1.65
> rvdw = 1.65
> rcoulomb = 1.65
>
> grompp runs fine, with no errors.  However, mdrun gives me the following
> note message at the beginning of the run (the simulation runs and finishes
> fine):
>
> ---
> NOTE: Periodic molecules: can not easily determine the required minimum
> bonded cut-off, using half the non-bonded cut-off
> ---
>
> If you have time, do you know if this message means that the non-bonded
> cut-off (is this referring to rvdw and rcoulomb?) which is used in the
> simulation is actually different than the cut-off I specify in the .mdp
> file?

No, it is referring to the need to determine the minimum size of the DD 
cells, which depends on a whole range of factors including the largest 
diameter of bonded interactions (usually dihedrals). Periodic molecules 
make determining that largest diameter a pain, so you're being warned 
that the algorithm is falling back on a heuristic for approximating that 
minimum size. Your cut-offs are not being changed - none of them are a 
"minimum bonded cut-off".

>    In other words, given the .mdp parameters above and the note message,
> are these parameters being used?
>
> vdwtype = cut-off
> coulombtype = PME
> rlist = 1.65
> rvdw = 1.65
> rcoulomb = 1.65
>
> Or are these actually being used (one-half of what I specified) by the code?
>
> vdwtype = cut-off
> coulombtype = PME
> rlist = 0.825
> rvdw = 0.825
> rcoulomb = 0.825
>
> I have checked the .log file resulting from mdrun.  In the "Input
> Parameters" section, it lists (among many others):
>
> vdwtype              = Cut-off
> coulombtype          = PME
> rlist                = 1.65
> rlistlong            = 1.65
> rvdw                 = 1.65
> rcoulomb             = 1.65
>
> whereas in the following section in the .log file, it says "Minimum cell
> size due to bonded interactions: 0.825 nm."

As it promised above, it's using half of the nonbonded cut-off here.

>   Later in the .log file, it says
> "Cut-off's:   NS: 1.65   Coulomb: 1.65   LJ: 1.65".  Near the end of the
> .log file, it says:
>
> ---
> The maximum allowed distance for charge groups involved in interactions is:
>                   non-bonded interactions           1.650 nm
>              two-body bonded interactions  (-rdd)   1.650 nm
>            multi-body bonded interactions  (-rdd)   0.849 nm
>    atoms separated by up to 5 constraints  (-rcon)  0.849 nm
> ---
>
> My question is, why does it use 0.849 nm for the "multi-body bonded
> interactions" and "atoms separated by up to 5 constraints," rather than
> 1.650 nm?

The heuristic says so. There's been a scaling of the resulting (1.65/2) 
to suit some other aspect of the algorithm.

>    Is there any way that I can force these terms to be 1.650 nm
> instead?  (Or should I?)  Also, if you have time, what is a multi-body
> bonded interaction?  I cannot seem to find this in the manual.

One that involves more than two bodies, so angles, dihedrals, etc.

The hints with -rdd and -rcon refer to mdrun command line arguments of 
those names. You should not get involved unless you know your bonded 
interactions are so large that their diameter might not fit. However, 
even five colinear C-C bonds will fit in 0.849 nm. If there was a 
problem, mdrun will crash and report missing bonded interactions.

Mark