Custom Search
----
I'm pasting part of the mdout file from Sander run:
Free energy options:
icfe = 1, ifsc = 1, klambda = 1
clambda = 0.0100, scalpha = 0.5000, scbeta = 12.0000
-----
Further at the end of the run where it says the averages for DV/DL over
xxxx steps (pasted below)
DV/DL, AVERAGES OVER 500000 STEPS
NSTEP = 500000 TIME(PS) = 550.000 TEMP(K) = 0.01 PRESS =
0.0
Etot = 0.0576 EKtot = 0.0576 EPtot =
82.0805
BOND = 0.0000 ANGLE = 0.0000 DIHED =
0.0000
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-7.1591
EELEC = 0.0000 EHBOND = 0.0000 RESTRAINT =
89.2397
EAMBER (non-restraint) = -7.1591
DV/DL = 82.9328
Ewald error estimate: 0.0000E+00
The restraint energy seems as much as the calculated DV/DL. The information
is for lambda = 0.01. I'm not sure if it has taken out the restraint energy
from DV/DL.
I'm pasting (below) the lambda vs DV/DL values for all the lambda values I
tried.
lambda DV/DL
0.01 82.9328
0.02 38.1868
0.03 22.6341
0.04 14.8417
0.05 10.4916
0.10 1.6545
0.15 -4.3380
0.20 -8.1609
0.25 -4.5246
0.30 1.3307
0.35 6.2843
0.40 9.3272
0.45 9.6986
0.50 6.6347
0.55 5.5444
0.60 4.1688
0.65 3.5596
0.70 3.0633
0.75 2.9437
0.80 2.6626
0.85 2.2263
0.90 2.2796
0.95 1.9180
0.96 1.8828
0.97 2.0306
0.98 1.9645
0.99 1.6034
The input file for the V0 state:
NPT production
&cntrl
imin = 0, ntx = 5, irest = 1,
ntpr = 10000, ntwr = 100000, ntwx = 10000,
ntf = 1, ntc = 1,
ntb = 2, cut = 10.0,
nstlim = 500000, dt = 0.001,
temp0 = 300.0, ntt = 3, gamma_ln = 2,
ntr = 1, restraint_wt = 50.0,
restraintmask = ':1'
ntp = 1, pres0 = 1.0, taup = 2.0,
icfe=1, clambda = 0.01,
ifsc=1, dvdl_norest = 1,
crgmask=':K+.K+',
scmask=':K+.K+',
&end
The mdin file for production run for state V1:
NPT production
&cntrl
imin = 0, ntx = 5, irest = 1,
ntpr = 10000, ntwr = 100000, ntwx = 10000,
ntf = 1, ntc = 1,
ntb = 2, cut = 10.0,
nstlim = 500000, dt = 0.001,
temp0 = 300.0, ntt = 3, gamma_ln = 2,
ntr = 1, restraint_wt = 50.0,
restraintmask = ':1'
ntp = 1, pres0 = 1.0, taup = 2.0,
icfe=1, clambda = 0.01,
ifsc=1, dvdl_norest = 1,
crgmask=':Cl-.Cl-',
scmask=':Cl-.Cl-',
&end
Can you suggest something based on my inputs why the calculated DV/DL is so
high in the beginning of lambda = 0.01, 0.02..
When I tried before without specifying dvdl_norest = 1, the endpoint
singularity was observed close to lambda =0 and lambda = 1. But now it
appears to be the case when lambda is close to 0.
Thanks for your time.
Regards
Sai
On Thu, Jan 19, 2012 at 1:22 PM, Sai Kumar Ramadugu <sramadugu.gmail.com>wrote:
> Dr Steinbrecher,
> The ultimate goal in this project is to mutate the charged aminoacids in
> wild type protein to glycine or alanine and carry on the TI to find
> relative free energy differences.
> When I was looking through literature, there were very few papers that
> talk about mutating a charged residue to neutral residue.
> One of the papers that I have come across is title "Accurate Estimates of
> Free Energy Changes in Charge Mutations", JCTC, 6, 1884.
>
> The authors suggest that when charged species are mutated, due to the
> Ewald treatment of electrostatics, there is a self-energy term that appears
> and one has to account for it and they suggested a method to do dual
> mutations.
> In their work they have shown it work for Glu to Gly or Ala and mutation
> of K+ and Cl- to neutral dummy atoms.
> I'm trying to follow that protocol for my charged mutations. When they did
> for KCl in water to neutral atoms, they were suggesting that the charged
> species in a concurrent dual mutation should be kept at a certain distance.
> I quote from the paper.... "In order to guarantee the independence of the
> two concurrent mutations, it is necessary to position the counterion at an
> adequate distance from the protein in order to minimize the electrostatic
> interaction between the protein and the ion. The Bjerrum length is defined
> as the distance at which the Coulomb interaction between two monovalent
> ions in a uniform dielectric is equal to the thermal energy, kBT".
>
> This was reason for me to set the position restraints for K+ and Cl-. In
> the initial file they are separated by 15 A in my case.
>
>
>
> On Thu, Jan 19, 2012 at 4:49 AM, <steinbrt.rci.rutgers.edu> wrote:
>
>> Hi,
>>
>> > When I observed the dV/dL averages for 500ps for lambda = 0.99, the
>> > restraint energy is -177.4749 kcal/mol which is the almost all the dV/dL
>> > for my system.
>> > The reason I had the restraints was to make sure that the K+ and Cl-
>> dont
>> > come close during the dual mutation.
>> >
>> > Is there a better way to approach this?
>>
>> I am not sure what you ultimately want to do, so it is difficult to
>> distinguish better and worse ways to do it here. Decoupling a restrained
>> molecule from its surroundings is a different thing from decoupling a
>> non-restrained one, with resulting differences in free energy. Both ways
>> make sense in terms of different calculations.
>>
>> If you want the two ions to be restrained in V0 only, then your setup is
>> good to go, but shows (unavoidable?) convergence problems at high lambdas.
>> If you want the restraints to not count for dvdl, set dvdl_norest, or
>>
>
> I want the restraints for both V0 and V1 but since V1 was only with water,
> I was not sure how to apply the position restraints for the corresponding
> atoms in V1 state.
> Is this assumption valid?
>
> I think your suggestion of setting dvdl_norest can be helpful in my case.
> I'll set up some calculations and report here.
>
>
> Thanks for your help.
>
> Regards
> Sai
>
>
>
>
>
>> define the restraint as a distance between K and Cl instead of cartesian
>> ones. Or, remove the restraints altogether (which, indeed may also give
>> you poor sampling as K and Cl can stick together).
>>
>> None of these things is wrong per se, they just answer different
>> questions.
>>
>> Kind Regards,
>>
>> Thomas
>>
>> Dr. Thomas Steinbrecher
>> formerly at the
>> BioMaps Institute
>> Rutgers University
>> 610 Taylor Rd.
>> Piscataway, NJ 08854
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jan 24 2012 - 15:00:03 PST