Alchemical - non-bonded table-extension warning

GROMACS user forum
#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