[gmx-users] Free energy perturbation - table extension and blowing system

[Krzysztof Makuch]

Hi,
I want to calculate 'solvation' energy of lutein in bilayer and water by
Free Energy Perturbation. When I'm running a normal MD, or some kind of
steering dynamics,my system works just fine. After addition of FEP
parameters, em runs, but instantly  returns series of warnings:

WARNING: Listed nonbonded interaction between particles 6461 and 6482
at distance 2.545 which is larger than the table limit 2.200 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 same thing repeats for several pairs of atoms, only in one
transmembrane molecule. As I said, everything is fine without FEP
parameters. For this instance I'm using the following parameters adapted
stright from the Bevanlab tutorial:

; Run control
integrator               = steep
nsteps                   = 5000
; EM criteria and other stuff
emtol                    = 100
emstep                   = 0.01
niter                    = 20
nbfgscorr                = 10
; Output control
nstlog                   = 1
nstenergy                = 1
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme            = verlet
nstlist                  = 1
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.2
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
; van der Waals
vdwtype                  = cutoff
vdw-modifier             = potential-switch
rvdw-switch              = 1.0
rvdw                     = 1.2
; 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 and pressure coupling are off during EM
tcoupl                   = no
pcoupl                   = no
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = 1
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state        0    1    2    3    4    5    6    7    8
9    10   11   12   13   14   15   16   17   18   19   20
vdw_lambdas              = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_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 0.00
; We are not transforming any bonded or restrained interactions
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 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 0.00
; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)
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 0.00
; Not doing simulated temperting here
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 0.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = lutein
couple-lambda0           = vdw
couple-lambda1           = none
couple-intramol          = no
nstdhdl                  = 10
; No velocities during EM
gen_vel                  = no
; options for bonds
constraints              = h-bonds
; Type of constraint algorithm
constraint-algorithm     = lincs
; Do not constrain the starting configuration
continuation             = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12


What I've tried by now and didn't work:
- I've tried to increase table extension, but it didn't do any good.
- I've tried to change the charge groups, which are now 1 atom = 1 group
for this molecule
- I was trying to solve this with gromacs 4, 5 and 16'. Warnings were
similar.
If I try to run further FEP steps the system with FEP parameters blows up.
MD without FEP simulates doesn't cause any problems.
My guess is I'm missing something inside FEP parameters. I'm using GROMACS
for several years by now and for the first time I can't solve the problem
by myself. I'd be glad for any tips You could give.
Best,
Krzysiek
-- 
Jagiellonian University
Department of Computational Biophysics and Bioinformatics
tel.1: (12) 664 61 49
tel.2: (48) 664 086 049