Re: AMBER: large RMS fluctuations in turning off VDW interactions by TI

From: Julien Michel <j.michel.soton.ac.uk>
Date: Sun, 01 Jul 2007 23:42:53 +0100

Hi Holly,

Turning off LJ interactions can be problematic in FEP/TI simulations
because the excluded volume of the atoms being turned off effectively
only vanish in the very last steps of the coupling parameter. This is a
well known issue that can be addressed by different techniques like
retracting vanishing atoms in the VdW radius of other atoms, non linear
scaling of the force field terms or a softcore (a modified LJ energy
function). This paper may help you understand the issues and the solutions:

"A Comparison of Non-Bonded Scaling Approaches for Free Energy
Calculations" Pitera J.W.1; van Gunsteren W.F. Molecular
Simulation,Volume 28, Number 1, 1 January 2002 , pp. 45-65(21)

Julien Michel

Holly Freedman wrote:
> Dear AMBER support,
> I have been struggling with a problem I've been having with a thermodynamic integration calculation.
> I have a molecule with a little over 50 atoms, and I am perturbing a small part of it (O-COCH3) to (O-
> CH3). According to advice from AMBER, I am dividng the calculation into four parts: Appearing and
> disappearing atoms (ie the uncoupling of Van Der Waals interactions), and charging and uncharging
> atoms which have appeared or disappeared. I am having trouble with the calculations where I
> decouple Van Der Waals interactions because the RMS fluctuations in DV/DL which I am getting are
> way too high (average over 3 kcal/mol). I am doing this using 9 values of clambda as described in
> the manual.
>
> Here is my infile:
> First equilibration:
> &cntrl
> imin = 0,
> ntx = 1, irest = 0, ntr = 1,
> ntpr = 100, ntwr = 100, ntwx = 5000, ntwe = 5000,
> ntf = 2, ntb = 2, cut = 10,
> nstlim = 500000, dt = .002, ntc = 2,
> tempi = 0.0, temp0 = 298.0, ntt = 1,
> tautp = 2.0, ntp = 1,
> icfe = 1, clambda = 0.33787, klambda = 6,
> &end
> Restraining region
> 10.0
> RES 1 1
> END
> END
>
> Then:
> &cntrl
> imin = 0,
> ntx = 5, irest = 1, ntr = 1, iwrap=1,
> ntpr = 100, ntwr = 100, ntwx = 5000, ntwe = 5000,
> ntf = 2, ntb = 2, cut = 10,
> nstlim = 2500000, dt = .002, ntc = 2,
> tempi = 298.0, temp0 = 298.0, ntt = 1,
> tautp = 1.0, ntp = 1,
> icfe = 1, clambda = 0.33787, klambda = 6,
> &end
> Restraining region
> 10.0
> RES 1 1
> END
> END
>
> My guess was that this was because the frcmod file I used for the dummy atoms wasn't consistent
> with the force field for the analogous real atoms, but I can't see my mistake. Here is a my frcmod:
>
> remark goes here
> MASS
> D2 1.008
> D4 12.01
> D5 16.00 same as O
>
> BOND
> D5-c 648.0 1.214 same as c-o
> c-D4 328.3 1.508 same as c-c3
> D4-D2 337.3 1.092 same as c3-hc
>
> ANGLE
> D5-c-os 76.2 122.43 same as o-c-os
> D5-c-D4 68.0 123.11 same as c3-c-o
> D4-c-os 69.3 111.96 same as c3-c-os
> c-D4-D2 47.2 109.68 same as c-c3-hc
> D2-D4-D2 39.4 108.35 same as hc-c3-hc
>
> DIHE
> X-c -D4-X 0 0.000 0.000 0.000 same as X -c -c3-X
> X-c -c3-X 0 0.000 0.000 0.000 same as X -c -c3-X
>
> IMPROPER
>
> NONBON
> D2 1.0000 0.0000
> D4 1.0000 0.0000
> D5 1.0000 0.0000 ATTN, need revision
>
>
> I was wondering if anyone else has had a similar problem, and discovered what was going wrong?
> I would really appreciate any helpful suggestions anyone has on this.
>
> Thanks in advance,
> Holly Freedman
> --
> Department of Physics, University of Alberta
> Edmonton CANADA
>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>
>

-- 
Julien Michel
http://www.julienmichel.net
Dr J. W. Essex research group
University of Southampton
United Kingdom
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Jul 04 2007 - 06:07:09 PDT
Custom Search