AMBER Archive (2002)

Subject: question on Gibbs calculation

From: X. Tan (xjtan_at_u.washington.edu)
Date: Tue Mar 12 2002 - 15:50:24 CST

Dear all,

I'd like to perform the free energy perturbation for the uneqivent charged
system. Before I start the calculation on my system, I carried out the
gibbs calculation, using CH3COOH------->CH3COO- as a test system,
designing the following cycle:

CH3COOH(vacum)----3----->CH3COO-(vacum)
        | |
        |(detaG1) |(detaG2)
        | |
CH3COOH(aq)-------4----->CH3COO-(aq)

then calculated the free energy changes between 3 and 4 processes. Since
the free solvation energies of CH3COOH and CH3COO- are known, in princle,
detaG1-detaG2=detaG3-detaG4 but according my calculation, detaG3-detaG4 is
about 61 kcal/mol for non-PME calculation(both decoupling method and
general gibbs were tried, the result are same) , detaG3-detaG4 is aobut 88
kcal/mol, when taken PME method into account, while the expt. value,
detaG1-detaG2 is 70-72 kcal/mol.

In my calculation, I used RESP charge, add 10 A' TIP3P water box, using
Parm99 force field. the following are some parameters for the calculations
in water(similar parameters were used for vacum phase and non PME
calculation).

 -----------------------------------
 Min.in
 &cntrl imin = 1, maxcyc =2000, nsnb = 25,
  ntc = 1, ntf = 1,
  cut =10.0, ntb = 1, scee = 1.2,
 &end
 -------------------------------------
 md0.in
 &cntrl
  imin = 0, ntx = 1,
  nstlim = 7000, nsnb = 20,
  tempi=0.0, temp0 = 300.0,
  ntt = 1, npscal = 1, ntb = 2, ntp = 1, taup = 0.01,
  ntc = 2, ntf = 2, cut =10.0, scee =1.2,
 &end
  28.9118210 29.1471220 28.1085790 90.0000000 90.0000000 90.0000000
  30 30 30 4 0 0 0
  0.00001
---------------------------------------
 md1.in
 &cntrl
    imin = 0, irest = 1, ntx = 7,
    nstlim =10000, nsnb = 20,
    tempi=0.0, temp0 = 300.0,
    ntt = 1, npscal = 1, ntb = 2, ntp = 1, taup = 0.01,
    ntc = 2, ntf = 2, cut =10.0, scee =1.2,
    ntt = 1, npscal = 1, ntb = 2, ntp = 1, taup = 0.01,
    ntc = 2, ntf = 2, cut =10.0, scee =1.2,
 &end
 24.7464390 24.9478397 24.0589216 90.0000000 90.0000000 90.0000000
 25 25 24 4 0 0
 0.00001
 ----------------------------------------
 gibbs.in (using PME method)
 &cntrl
    timlim = 0.1000000E+09 , irest = 0 , ibelly = 0 ,
    ichdna = 0 , ipol = 0 , i3bod = 0 ,
    iewald = 1 ,
    ntx = 7 , ntxo = 1 , ig = 71277 ,
    tempi = 300.0000 , heat = 0.0000000E+00 ,
    ntb = 2 , iftres = 0 , boxx(1)= 0.0000000E+00 ,
    boxx(2)= 0.0000000E+00 , boxx(3)= 0.0000000E+00 , betar =
    0.0000000E+00 , ibxrd = 0 ,
    nrun = 21 , ntt = 5 , temp0 = 300.0000 ,
    dtemp = 999.0000 , tautp = 0.2000000 , tauts = 0.2000000 ,
    isolvp = 0 , nsel = 0 , dtuse = 0.5000000 ,
    ntp = 1 , npscal = 0 , pres0 = 1.000000 ,
    comp = 44.60000 , taup = 0.6000000 ,
    ndfmin = 0 , ntcm = 0 , nscm = -1 ,
    istay = 0 , nstay = 0 , isvat = 1 ,
    nstlim = 100000 , init = 4 , ntu = 1 ,
    t = 0.0000000E+00 , dt = 0.1000000E-02 , vlimit = 20.00000 ,
    ivemax = 0 ,
    ntc = 3 , ntcc = 0 , nconp = 0 ,
    tol = 0.1000000E-07 , tolr2 = 0.1000000E-03 , ncorc = 1 ,
    ishkfl = 6 , itimth = 0 , jfastw = 0 ,
    ntf = 3 , ntid = 0 , ntn = 3 ,
    ntnb = 1 , nsnb = 25 , idiel = 1 ,
    inbper = 0 , ielper = 0 , imgslt = 0 ,
    idsx0 = 400 , itrslu = 1 , ioleps = 0 ,
    intprt = 5 , itip = 0 ,
    cut = 10.00000 , scnb = 2.000000 , scee = 1.200000 ,
    dielc = 1.000000 , cut2nd = 0.0000000E+00 , cutprt = 0.0000000E+00 ,
    ntpr = 200000 , ntwx = -1 , ntwv = -1 ,
    ntwe = -1 , ntwxm = 0 , ntwvm = 0 ,
    ntwem = 0 , ntpp = 0 , ioutfm = 0 ,
    isande = 1 , iperat = 0 , iatcmp = 0 ,
    ntatdp = 0 , icmpdr = 0 , ncmpdr = 0 ,
    ntwprt = 0 , ntwpr0 = 0 ,
    ntr = 0 , nrc = 0 , ntrx = 1 ,
    taur = 0.1000000 , intr = 0 , ibigm = 0 ,
    nmrmax = 0 , iwtmax = 0 , isftrp = 0 ,
    rwell = 4.000000 ,
    iftime = 0 , ctimt = 40.00000 , almda = 1.000000 ,
    almdel = 0.5000000E-01 , isldyn = -3 , idifrg = 1 ,
    nstmeq = 25000 , nstmul = 75000 , ndmpmc = 0 ,
    idwide = 0 , ibndlm = 0 ,
    iavslp = 8 , iavslm = 2 , islp = 0 ,
    corrsl = 0.8000000 , amxmov = 0.1000000E-01 ,
    iavdel = -1 , iavdem = 2 , amxdel = 0.1000000 ,
    almdl0 = 0.1000000E-02 , dlmin = 0.1000000E-04 , dlmax =
0.1000000E-01,
    amxrst = 0.9999999E-01 , norsts = 0 , ntsd = 0 ,
    almstp = -1.000000 ,
    nstpe = 25000 , nstpa = 75000 , dte = 0.1000000E-02 ,
    dta = 0.1000000E-02 ,
    ivcap = 2 , natcap = 0 , fcap = 1.500000 ,
 &end
 24.6036942 24.8039333 23.9201427 90.0000000 90.0000000 90.0000000
 25 25 24 4 1.0 0.0
 0.00001
 ----------------------------------------------------------------

In general, the difference between the experimental values and calculation
should be within 2 kcal/mol. I don't know why there is such a large gap in
my case. I think the calculations were converged, for I also try 101
windows, it gave the same results as 21 windows. So, I really don't know
why the results looked like this and how to deal with this problem, or is
there something wrong in my calculations. There are so few papers
on this kind of calculation! Could anyone give me some suggestion?
Thanks in advance.

Xiaojian Tan