Alchemical - non-bonded table-extension warning

User discussions
1

GROMACS version:5.1.4 and 2018
GROMACS modification: No

Hi,

When performing an alchemical calculation of a relatively large molecule, for instance a large hydrocarbon (say icosane Nb Carbons = 20) in water, I get multiple warnings concerning non-bonded interactions (see below)

WARNING: Listed nonbonded interaction between particles x and y
at distance 2.047 which is larger than the table limit 2.000 nm.

This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.

IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.

The x and y above can be either both atoms of the hydrocarbon or atoms of the hydrocarbon and water.

The force field I used was OPLS-aa and I tried 2 different versions of gromacs (5.1.4 and 2018).
I am running MD (NPT) and I get no LINCS warnings.

These warnings do not occur if I carry out a “normal” MD (i.e., no decoupling).

For smaller hydrocarbons, even with with 1,4 interactions I did not get such a warning (e.g., dodecane Nb carbons = 12).

The simulations for different lambdas end and it is possible to calculate the free energy afterwards.

Although the warning clearly points out a direct relation with the decoupling approach I am not sure what it means exactly and its consequences.

Increasing table-extension “solves” the problem. However, I do not understand the nature of the problem, and I do not want to simply increase table-extension without understanding the source of it.

Any hint on this would be highly appreciated.

Thanks in advance

Nuno

2

Given the large size of the molecule being transformed, you are probably getting high-energy clashes among the atoms in the chain. Please post your .mdp file to further troubleshoot. Transforming very large molecules often suffers from such problems, and more often, great difficulty in converging the calculations.

3

Hi,

Thanks a lot for your reply. I actually found out, in the mean time, that by switching “couple-intramol” to “yes” apparently solves the problem.

couple-intramol

no

All intra-molecular non-bonded interactions for “moleculetype” are replaced by exclusions and explicit pair interactions. In this manner the decoupled state of the molecule corresponds to the proper vacuum state without periodicity effects.

yes

The intra-molecular Van der Waals and Coulomb interactions are also turned on/off. This can be useful for partitioning free-energies of relatively large molecules, where the intra-molecular non-bonded interactions might lead to kinetically trapped vacuum conformations. The 1-4 pair interactions are not turned off.

However, I am not sure I understand the “on/off” above and, therefore, the differences between the “YES” and “NO” options to their full extent. The fact that I got warnings regarding interactions between atoms from the solute and the solvent and not just between the solute is also confusing. These also disappeared now.

This is something I am still trying to figure out. Any help on this would be really appreciated.

My .mdp is given below:

; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 5000000
nstcomm = 100

; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 0

; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = verlet
nstlist = 20
ns_type = grid
pbc = xyz
rlist = 1.0

; Electrostatics
coulombtype = PME
rcoulomb = 1.0

; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 0.9
rvdw = 1.0

; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres

; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12

; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-06
epsilon_surface = 0

; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps = system
tau_t = 1.0
ref_t = 298

; Pressure coupling is on for NPT
Pcoupl = Parrinello-Rahman
tau_p = 2.0
compressibility = 4.5e-05
ref_p = 1.0

; Free energy control stuff
free_energy = yes
init_lambda_state = 0
delta_lambda = 0
calc_lambda_neighbors = 1

; lambdas 20 points
coul_lambdas = 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw_lambdas = 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0

bonded_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
mass_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
temperature_lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

; Options for the decoupling
sc-alpha = 0.5
sc-coul = no
sc-power = 1
sc-sigma = 0.3
couple-moltype = MOL
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
;couple-intramol = yes
;nstdhdl = 10

; Do not generate velocities
gen_vel = no

; options for bonds
constraints = h-bonds

; Type of constraint algorithm
constraint-algorithm = lincs

; Constrain the starting configuration
; since we are continuing from NPT
continuation = yes

; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12

Thanks once again.

Best wishes.

Nuno

4

Setting couple-intramol to yes will soften the intramolecular nonbonded interactions and allow the calculation to proceed, but depending on why you’re performing the calculation, you may or may not want that setting to be active. If you’re looking for an absolute free energy of solvation, couple-intramol should be no, otherwise you’ll need to do more calculations (corresponding gas-phase transformation with couple-intramol = yes).

5

ok, thanks. Yes, I am interested in taking the absolute value of the solvation free energy.
If I understand correctly that would amount to the way this calculation used to be performed in older versions of gromacs. The “on/off” then just means lambda will be used in a similar way for the intra-molecular Hamiltonian. Non-bonded intra-molecular interactions would be “on” for lambda = 0 (for my .mdp liquid to gas) and “off” for Lambda = 1 (complete decoupling). Is that it?

I am waiting for some test results using a larger cut-off for the “table-extension”. I am not expecting significant differences since we are talking about very large distances for essentially van der Waals interactions, however, I still do not fully understand the source of the problem, especially because of the solute-water interactions’ warnings.

Thanks once again.

Best

Nuno

6

That’s correct.

7

ok Thanks. I finished a test and apparently increasing “table-extension” to a value large enough to get rid of every warning has no effect on the final free energy.
Thus, although I am still not sure what triggers those warnings I guess one can trust this calculation. Not sure if that will be the case for molecules where electrostatic interactions play a more important role. Probably not.
In that case increasing “table-extension” may be necessary, unless some undesired effects may come from that. Anyway I let this information for anyone who may be interested in using alchemical methods for large molecules. If you know of any inconvenient (other that possibly a slight increase in the time of computation) in using a larger value for “table-extension” I would appreciate if you could pointed those out.

I will keep trying to figure this out to a deeper extent and will post here any conclusions.

Best,

Nuno

8

Hello Justin,

I’m running to the same problem for an AA mutation in a peptide protein binding system.
My .mdp file (taken from http://pmx.mpibpc.mpg.de) is very different from Nuno’s. Also, the distance causing the problem is ridiculously large (1766462.875)!! It happens right after 100 ps.

Interestingly, the same problem happens for a very different system (used in the tutorial: Mutation free energy calculations) if I go beyond 100 ps, (which is the time set in the tutorial).

Any suggestion?
Thank you in advance.
Here is my .mdp file:

;

; VARIOUS PREPROCESSING OPTIONS
title =
include =
define =

; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 100000
;200 ps

ld_seed = -1
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
comm-mode = Linear
; number of steps for center of mass motion removal
nstcomm = 1
nstcalcenergy = 1
; group(s) for center of mass motion removal
comm-grps =

; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells
niter = 0
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 0
nstvout = 0
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 1000
; Output frequency for energies to log file and energy file
nstlog = 10000
nstenergy = 100000
; Output frequency and precision for xtc file
nstxtcout = 10000
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
cutoff-scheme = verlet
nstlist = 20
; ns algorithm (simple or grid)
ns-type = Grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 1.2
domain-decomposition = no

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.1
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-r = 1
; Method for doing Van der Waals
vdw-type = switch
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.1
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no

; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0

; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = v-rescale
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau-t = 0.1
ref-t = 300
; Pressure coupling
Pcoupl = Parrinello-Rahman
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 5
compressibility = 4.6E-5
ref-p = 1
; Random seed for Andersen thermostat
andersen_seed = 815131

; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = no
; Number of time points to use for specifying annealing in each group
annealing_npoints = 2
; List of times at the annealing points for each group
annealing_time = 0 250
; Temp. at each annealing point, for each group.
annealing_temp = 0 300

; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel = no
gen-temp = 300
gen-seed = 173529

; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = yes
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 6
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 2
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials
morse = no

; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =

; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Equal
; Use sqrt of the time averaged times the instantaneous violation
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file
nstdisreout = 100
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file
nstorireout = 100
; Dihedral angle restraints: No, Simple or Ensemble
dihre = No
dihre-fc = 1000
dihre-tau = 0
; Output frequency for dihedral values to energy file
nstdihreout = 100

; Free energy control stuff
free_energy = yes
init_lambda = 0
delta_lambda = 0.00002
sc-alpha = 0
sc-sigma = 0.3

refcoord-scaling = all

; Non-equilibrium MD stuff
acc-grps =
accelerate =
freezegrps =
freezedim =
cos-acceleration = 0

; Electric fields
; Format is number of terms (int) and for all terms an amplitude (real)
; and a phase angle (real)
E-x =
E-xt =
E-y =
E-yt =
E-z =
E-zt =

; User defined thingies
user1-grps =
user2-grps =
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0