[gmx-users] Need Help in Running Thermodynamic Integration (TI) Simulations

[Yew Mun Yip]

Hi, I recently tried running TI simulations with a ligand solvated in water
using the CHARMM force field.

When I was running equilibration, everything was fine and the thermodynamic
properties (density, volume, pressure, etc) stabilized.

Then I started my TI simulations.

When I ran till init-lambda = 2, the following error message occurred.

step XXX: Water molecule starting at atom XXXX can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates.

Although visualizing the PDBs generated with VMD does show that the water
molecule have very large bond lengths between the atoms, I'm more concerned
with what could be the problem with my .mdp file that may have lead to this
error?

This is my TI's .mdp file:

; we'll use the sd integrator with 100000 time steps (200ps)
integrator               = sd
nsteps                   = 100000
dt                       = 0.002
nstenergy                = 1000
nstlog                   = 5000
; turn off trajectory writing
nstxout                  = 0
nstvout                  = 0
; We use the old group scheme because as of writing, the Verlet scheme
; does not support free energy calculations with coupled molecules.
cutoff-scheme            = group
; cut-offs at 1.0nm
rlist                    = 1.0
dispcorr                 = EnerPres
vdw-type                 = cut-off
rvdw                     = 1.0
; Coulomb interactions
coulombtype              = pme
rcoulomb                 = 1.0
fourierspacing           = 0.12
; Constraints
constraints              = all-bonds
; set temperature to 300K
tcoupl                   = v-rescale
tc-grps                  = system
tau-t                    = 0.2
ref-t                    = 300
; set pressure to 1 bar with a thermostat that gives a correct
; thermodynamic ensemble
pcoupl                   = parrinello-rahman
ref-p                    = 1
compressibility          = 4.5e-5
tau-p                    = 5
; and set the free energy parameters
free-energy              = yes
couple-moltype           = LIG
; these 'soft-core' parameters make sure we never get overlapping
; charges as lambda goes to 0
sc-power                 = 1
sc-sigma                 = 0.3
sc-alpha                 = 1.0
; we still want the molecule to interact with itself at lambda=0
couple-intramol          = no
couple-lambda1           = vdwq
couple-lambda0           = none
init-lambda-state        = 2
; These are the lambda states at which we simulate
; for separate LJ and Coulomb decoupling, use
fep-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.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
vdw-lambdas              = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

Any advice is greatly appreciated. =)

-- 
*Yip Yew Mun (Mr)* | PhD Research Scholar | Division of Chemistry
& Biological Chemistry
School of Physical & Mathematical Sciences | Nanyang
Technological University | Singapore 639798
Tel: (+65) 97967803 | Email: yipy0010 at e.ntu.edu.sg | GMT+8h