How does verlet neighbour searching list set rlist for dynamic simulations

I am running dynamics simulations using Gromacs 2021.5. In my MDP file, I use Verlet for Neighbour searching. I understand from the manual that the buffer is set automatically and rlist is ignored.

However, when I run my simulations the rlist is always set to 1nm as per my log files unless I assign the rlist a different value.
This is what my log file looks like without setting rlist
Input Parameters:
integrator = md
tinit = 0
dt = 0.001
nsteps = 3000000
init-step = 0
simulation-part = 1
mts = false
comm-mode = Linear
nstcomm = 100
bd-fric = 0
ld-seed = -8388813
emtol = 10
emstep = 0.01
niter = 20
fcstep = 0
nstcgsteep = 1000
nbfgscorr = 10
rtpi = 0.05
nstxout = 2000
nstvout = 2000
nstfout = 0
nstlog = 2000
nstcalcenergy = 100
nstenergy = 2000
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
periodic-molecules = false
verlet-buffer-tolerance = 0.005
rlist = 1
coulombtype = PME
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1
epsilon-r = 1
epsilon-rf = inf
vdw-type = PME
vdw-modifier = Potential-shift
rvdw-switch = 0
rvdw = 1
DispCorr = No
table-extension = 1
fourierspacing = 0.12
fourier-nx = 96
fourier-ny = 96
fourier-nz = 96
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 0
epsilon-surface = 0
tcoupl = V-rescale
nsttcouple = 10
nh-chain-length = 0
print-nose-hoover-chain-variables = false
pcoupl = Berendsen
pcoupltype = Isotropic
nstpcouple = 10
tau-p = 1
compressibility (3x3):
compressibility[ 0]={ 4.50000e-05, 0.00000e+00, 0.00000e+00}
compressibility[ 1]={ 0.00000e+00, 4.50000e-05, 0.00000e+00}
compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 4.50000e-05}
ref-p (3x3):
ref-p[ 0]={ 1.00000e+02, 0.00000e+00, 0.00000e+00}
ref-p[ 1]={ 0.00000e+00, 1.00000e+02, 0.00000e+00}
ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 1.00000e+02}
refcoord-scaling = No
posres-com (3):
posres-com[0]= 0.00000e+00
posres-com[1]= 0.00000e+00
posres-com[2]= 0.00000e+00
posres-comB (3):
posres-comB[0]= 0.00000e+00
posres-comB[1]= 0.00000e+00
posres-comB[2]= 0.00000e+00
QMMM = false
qm-opts:
ngQM = 0
constraint-algorithm = Lincs
continuation = false
Shake-SOR = false
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
nwall = 0
wall-type = 9-3
wall-r-linpot = -1
wall-atomtype[0] = -1
wall-atomtype[1] = -1
wall-density[0] = 0
wall-density[1] = 0
wall-ewald-zfac = 3
pull = false
awh = false
rotation = false
interactiveMD = false
disre = No
disre-weighting = Conservative
disre-mixed = false
dr-fc = 1000
dr-tau = 0
nstdisreout = 100
orire-fc = 0
orire-tau = 0
nstorireout = 100
free-energy = no
cos-acceleration = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
simulated-tempering = false
swapcoords = no
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
applied-forces:
electric-field:
x:
E0 = 0
omega = 0
t0 = 0
sigma = 0
y:
E0 = 0
omega = 0
t0 = 0
sigma = 0
z:
E0 = 0
omega = 0
t0 = 0
sigma = 0
density-guided-simulation:
active = false
group = protein
similarity-measure = inner-product
atom-spreading-weight = unity
force-constant = 1e+09
gaussian-transform-spreading-width = 0.2
gaussian-transform-spreading-range-in-multiples-of-width = 4
reference-density-filename = reference.mrc
nst = 1
normalize-densities = true
adaptive-force-scaling = false
adaptive-force-scaling-time-constant = 4
shift-vector =
transformation-matrix =
grpopts:
nrdf: 26797
ref-t: 300
tau-t: 0.1
annealing: No
annealing-npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp-flags[ 0]: 0

Changing nstlist from 10 to 100, rlist from 1 to 1

1 GPU selected for this run.
Mapping of GPU IDs to the 1 GPU task in the 1 rank on this node:
PP:0
PP tasks will do (non-perturbed) short-ranged and most bonded interactions on the GPU
PP task will update and constrain coordinates on the CPU
Using 1 MPI thread

Non-default thread affinity set, disabling internal thread affinity

Using 20 OpenMP threads

System total charge: -0.000
Will do PME sum in reciprocal space for LJ dispersion interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Using a Gaussian width (1/beta) of 0.298423 nm for LJ Ewald
Will do PME sum in reciprocal space for electrostatic interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Using a Gaussian width (1/beta) of 0.320163 nm for Ewald
Potential shift: LJ r^-12: -1.000e+00 r^-6: -1.000e-03, Ewald -1.000e-05
Initialized non-bonded Ewald tables, spacing: 9.33e-04 size: 1073

Using shifted Lennard-Jones, switch between 0 and 1 nm
Generated table with 1000 data points for 1-4 COUL.
Tabscale = 500 points/nm
Generated table with 1000 data points for 1-4 LJ6.
Tabscale = 500 points/nm
Generated table with 1000 data points for 1-4 LJ12.
Tabscale = 500 points/nm

Using GPU 8x8 nonbonded short-range kernels

Using a 8x8 pair-list setup:
updated every 100 steps, buffer 0.000 nm, rlist 1.000 nm
At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be:
updated every 100 steps, buffer 0.000 nm, rlist 1.000 nm
Removing pbc first time

*This is what my log file with setting rlist=1.5 look like

Input Parameters:
integrator = md
tinit = 0
dt = 0.001
nsteps = 3500000
init-step = 0
simulation-part = 1
mts = false
comm-mode = Linear
nstcomm = 100
bd-fric = 0
ld-seed = 1945558523
emtol = 10
emstep = 0.01
niter = 20
fcstep = 0
nstcgsteep = 1000
nbfgscorr = 10
rtpi = 0.05
nstxout = 5000
nstvout = 5000
nstfout = 0
nstlog = 5000
nstcalcenergy = 100
nstenergy = 5000
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
periodic-molecules = false
verlet-buffer-tolerance = 0.005
rlist = 1.5
coulombtype = PME
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1.5
epsilon-r = 1
epsilon-rf = inf
vdw-type = PME
vdw-modifier = Potential-shift
rvdw-switch = 0
rvdw = 1.5
DispCorr = EnerPres
table-extension = 1
fourierspacing = 0.12
fourier-nx = 84
fourier-ny = 84
fourier-nz = 84
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 0
epsilon-surface = 0
tcoupl = V-rescale
nsttcouple = 10
nh-chain-length = 0
print-nose-hoover-chain-variables = false
pcoupl = Berendsen
pcoupltype = Isotropic
nstpcouple = 10
tau-p = 1
compressibility (3x3):
compressibility[ 0]={ 4.50000e-05, 0.00000e+00, 0.00000e+00}
compressibility[ 1]={ 0.00000e+00, 4.50000e-05, 0.00000e+00}
compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 4.50000e-05}
ref-p (3x3):
ref-p[ 0]={ 1.00000e+02, 0.00000e+00, 0.00000e+00}
ref-p[ 1]={ 0.00000e+00, 1.00000e+02, 0.00000e+00}
ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 1.00000e+02}
refcoord-scaling = No
posres-com (3):
posres-com[0]= 0.00000e+00
posres-com[1]= 0.00000e+00
posres-com[2]= 0.00000e+00
posres-comB (3):
posres-comB[0]= 0.00000e+00
posres-comB[1]= 0.00000e+00
posres-comB[2]= 0.00000e+00
QMMM = false
qm-opts:
ngQM = 0
constraint-algorithm = Lincs
continuation = false
Shake-SOR = false
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
nwall = 0
wall-type = 9-3
wall-r-linpot = -1
wall-atomtype[0] = -1
wall-atomtype[1] = -1
wall-density[0] = 0
wall-density[1] = 0
wall-ewald-zfac = 3
pull = false
awh = false
rotation = false
interactiveMD = false
disre = No
disre-weighting = Conservative
disre-mixed = false
dr-fc = 1000
dr-tau = 0
nstdisreout = 100
orire-fc = 0
orire-tau = 0
nstorireout = 100
free-energy = no
cos-acceleration = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
simulated-tempering = false
swapcoords = no
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
applied-forces:
electric-field:
x:
E0 = 0
omega = 0
t0 = 0
sigma = 0
y:
E0 = 0
omega = 0
t0 = 0
sigma = 0
z:
E0 = 0
omega = 0
t0 = 0
sigma = 0
density-guided-simulation:
active = false
group = protein
similarity-measure = inner-product
atom-spreading-weight = unity
force-constant = 1e+09
gaussian-transform-spreading-width = 0.2
gaussian-transform-spreading-range-in-multiples-of-width = 4
reference-density-filename = reference.mrc
nst = 1
normalize-densities = true
adaptive-force-scaling = false
adaptive-force-scaling-time-constant = 4
shift-vector =
transformation-matrix =
grpopts:
nrdf: 23112
ref-t: 300
tau-t: 0.1
annealing: No
annealing-npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp-flags[ 0]: 0

Changing nstlist from 10 to 100, rlist from 1.5 to 1.5

Initializing Domain Decomposition on 4 ranks
Dynamic load balancing: auto
Minimum cell size due to atom displacement: 0.482 nm
Initial maximum distances in bonded interactions:
two-body bonded interactions: 0.404 nm, LJ-14, atoms 5258 5261
multi-body bonded interactions: 0.404 nm, Ryckaert-Bell., atoms 5258 5261
Minimum cell size due to bonded interactions: 0.444 nm
Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.218 nm
Estimated maximum distance required for P-LINCS: 0.218 nm
Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
Using 0 separate PME ranks
Optimizing the DD grid for 4 cells with a minimum initial size of 0.603 nm
The maximum allowed number of cells is: X 16 Y 16 Z 16
Domain decomposition grid 4 x 1 x 1, separate PME ranks 0
PME domain decomposition: 4 x 1 x 1
Domain decomposition rank 0, coordinates 0 0 0

The initial number of communication pulses is: X 1
The initial domain decomposition cell size is: X 2.50 nm

The maximum allowed distance for atoms involved in interactions is:
non-bonded interactions 1.500 nm
(the following are initial values, they could change due to box deformation)
two-body bonded interactions (-rdd) 1.500 nm
multi-body bonded interactions (-rdd) 1.500 nm
atoms separated by up to 5 constraints (-rcon) 2.500 nm

When dynamic load balancing gets turned on, these settings will change to:
The maximum number of communication pulses is: X 1
The minimum size for domain decomposition cells is 1.500 nm
The requested allowed shrink of DD cells (option -dds) is: 0.80
The allowed shrink of domain decomposition cells is: X 0.60
The maximum allowed distance for atoms involved in interactions is:
non-bonded interactions 1.500 nm
two-body bonded interactions (-rdd) 1.500 nm
multi-body bonded interactions (-rdd) 1.500 nm
atoms separated by up to 5 constraints (-rcon) 1.500 nm

On host r2i7n8 1 GPU selected for this run.
Mapping of GPU IDs to the 4 GPU tasks in the 4 ranks on this node:
PP:0,PP:0,PP:0,PP:0
PP tasks will do (non-perturbed) short-ranged and most bonded interactions on the GPU
PP task will update and constrain coordinates on the CPU
Using 4 MPI processes

Non-default thread affinity set, disabling internal thread affinity

Using 5 OpenMP threads per MPI process

System total charge: -0.000
Will do PME sum in reciprocal space for LJ dispersion interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Using a Gaussian width (1/beta) of 0.447634 nm for LJ Ewald
Will do PME sum in reciprocal space for electrostatic interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Using a Gaussian width (1/beta) of 0.480244 nm for Ewald
Potential shift: LJ r^-12: -7.707e-03 r^-6: -8.779e-05, Ewald -6.667e-06
Initialized non-bonded Ewald tables, spacing: 1.14e-03 size: 1314

Using shifted Lennard-Jones, switch between 0 and 1.5 nm
Generated table with 1250 data points for 1-4 COUL.
Tabscale = 500 points/nm
Generated table with 1250 data points for 1-4 LJ6.
Tabscale = 500 points/nm
Generated table with 1250 data points for 1-4 LJ12.
Tabscale = 500 points/nm
Long Range LJ corr.: 3.0097e-15

Using GPU 8x8 nonbonded short-range kernels

Using a 8x8 pair-list setup:
updated every 100 steps, buffer 0.000 nm, rlist 1.500 nm
At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be:
updated every 100 steps, buffer 0.000 nm, rlist 1.500 nm
Removing pbc first time

Initializing Parallel LINear Constraint Solver

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
B. Hess
P-LINCS: A Parallel Linear Constraint Solver for molecular simulation
J. Chem. Theory Comput. 4 (2008) pp. 116-122
-------- -------- — Thank You — -------- --------

The number of constraints is 6900
There are constraints between atoms in different decomposition domains,
will communicate selected coordinates each lincs iteration

Linking all bonded interactions to atoms

Intra-simulation communication will occur every 10 steps.

I would like to understand how this works. How is the rlist determined ??

Thank you.

Under the verlet buffer scheme, the important parameter is verlet-buffer-tolerance.

This parameter is the acceptable level of error/drift due to particles moving such that they should be in another particle short-range interaction list, but are not (because neighbor search happen every nstlist step, not every step).

These borderline particles are included in the neighbor list, by taking all pairs within r_{list} = r_{cutoff} + r_{buffer}, with r_{buffer} some positive value. Under NVT condition, one can estimate the diffusion of the borderline particles , so that r_{buffer} can be set to include them and thus satisfy the error tolerance.

Additionally, at a given allowed error level, you can increase nstlist if you increase the buffer (and conversely). That gives a free parameter that can be tuned for performance. This is what it done when the rlist/nstlist are set automatically.

If you really want to set nstlist/rlist exactly (which is not recommended in the general case), you need to set verlet-buffer-tolerance to -1.

See this part of the reference manual, this part of the user manual, and this publication.

1 Like