[gmx-users] free energy binding problem

[Daniele Veclani]

Dear Useres

I am studying an antibiotic in water (box with antibiotic and 1000
molecules of water), I am currently doing a calculation of free
energy binding. My problem is that, for some values of lambda, I have this
error:

Step 1447079, time 2894.16 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 437.814673, max 1855.867310 (between atoms 3031 and 3033)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   3001   3002   53.3    0.1000   0.1000      0.1000
   3017   3018   40.1    0.2721   0.1090      0.1090
   3023   3024   90.0    0.1100   5.3888      0.1100
   3023   3025   90.0    0.1100   3.4338      0.1100
   3026   3027   90.0    0.1000   0.1903      0.1000
   3031   3032   90.0    0.1105   5.8869      0.1090
   3031   3033   90.0    0.3034 204.2554      0.1100

Back Off! I just backed up step1447079b.pdb to ./#step1447079b.pdb.1#

Back Off! I just backed up step1447079c.pdb to ./#step1447079c.pdb.1#
Wrote pdb files with previous and current coordinates

WARNING: Listed nonbonded interaction between particles 3014 and 3024
at distance 2.258 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.


My .MDP file is:

 Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 2500000   ; 1 ns
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.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 coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 300
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 6
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state      = 6
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       ; linear interpolation of Coulomb (none
in this case)
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = QUI  ; name of moleculetype to decouple
couple-lambda0           = vdw      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case
only vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = h-bonds  ; we only have C-H bonds here
; 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

Is it a symptom of poor stability of the system (not good energy
minimization, or not good nvt npt)?

Best regards
Daniele Veclani