Re: [AMBER] TI in Amber and the TUTORIAL A9

From: Marek Maly <>
Date: Tue, 06 Dec 2011 00:38:28 +0100

Hello Thomas,
thanks a lot for your comments !

I would be very grateful if you could check that my ones and
in case of for me the most important point #1 also provide
some your explanation or reference info.

Thank you very much in advance for your time !

Dne Mon, 05 Dec 2011 13:05:41 +0100 <> napsal/-a:

> Hi,
>> #1
>> I would like to know how from the MD run in some fixed value of lambda
>> (L)
>> is possible
>> to get <dV/dL>. I would naturally assume that two MD runs (first with L
>> and the second with L + dL) has to
>> be done from which <V(L)> and <V(L+dL)> are obtained and then <dV/dL>
>> might be estimated as
>> ( <V(L+dL)> - <V(L)> )/dL or as the ( <V(L+dL,t) - V(L,t)> )/dL where t
>> is the simulation time.
> the lambda-dependent potential can be written down, so you can work out
> analytical derivatives and calculate them at a single lambda-value.

In this case (which you are talking about) you probably taking in account
(if speaking about V(L) or dV/dL) just
the L-weighted interactions between "unique" and the "common" atoms it
means just this particular contribution
to the total potential energy of the whole simulated system, where the
"unique" atoms are vdw coupled/decoupled or charged/uncharged
and the "common" atoms are all the other atoms which are not going to be
changed during TI. Am I right ?

Then I fully agree with you that once you have the analytical formula for
V(L) (potential for "unique"-"common" atom pairs)
you have in turn also the analytical formula for dV/dL(L) and so you can
evaluate it easily during
simulation with the given L value.

But in my opinion also the rest of the potential energy which comes mainly
 from "common"-"common" atom pairs is
influenced (perturbed) by introduction of the L-weighted interactions
between "unique"-"common" pairs.
I think it is clear that for example "disappearing" of some molecule X
 from the solvent brings solvent
molecules from neighbourhood of X in closer distances as the space which
was previously occupied by
molecule X is now more and more free and hence "common"-"common"
interaction (here solvent-solvent) changes here dramatically
I think.

So I see V(L) in two different contexts.

1) As the L-weighted interaction potential between "unique"-"common" atoms
which has analytical formula
and so also analytical formula for dV/dL.

2) As the total potential energy of the whole system which is composed of
"unique"-"common", "common"-"common", but also
  "unique"-"unique" atom pairs interactions. Unfortunately for example the
dependence of "common"-"common" part of the total potential
  energy on L is nontrivial and might to be pretty hard to express it
analytically as well as it's derivation.
  That was the reason for which I was thinking exclusively about the
numerical estimation of the dV/dL (for the given L value) from
  the relevant MD run considering here V as the TOTAL potential energy
which dependents "non-trivially" as a whole on L.

  From your response is it clear to me, that the values of dV/dL(L) are
computed just with respect to the analytically
  well defined part ("unique"-"common" interacations) of the total
potential energy in another words derivated V(L) has here meaning 1).

  So if I am not wrong you claim that this equation

  G1-G0 = integral_0_1<dU(L)/dL> (1)

  derived directly from the free energy definition is possible write
  also like this:

  G1-G0 = integral_0_1_<dV(L)/dL> (2)

  where U(L) is the total potential energy and V(L) is the potential energy
just from "unique"-"common" interactions
  which is possible to express analytically.

  If yes, may you please explain why (2) is valid in general or to advice
the proper reference where this
  probably trivial ( unfortunately not for me :)) ) equivalence is derived

>> #2
>> I also did not understand fully (as it is not explained anywhere) why
>> for
>> the MD run with the given fixed L value the two simultaneous
>> sander simulations has to be done for any TI calculation or it.s stage
>> (one with the prmtop_V0 file and the second with prmtop_V1 file)
>> and how the data obtained from these two simultaneous simulations are
>> combined at the end.
>> Could you please put some light on this ? (e.g. using as the example the
>> study from your tutorial A9 for
>> the explanation what is going on in each of that two MD simultaneous
>> threads during the all three stages)
> For implementation details, please see our 2007 paper (Steinbrecher,
> Mobley, Case, JCP). If you want to know why the Amber run is conducted in
> the way it is, the best way to start may be the code, as the
> implementation is quite straightforward. Check out thermo_int.f, where
> the
> actual force and energy mixing is done.

Thanks for the reference and name of the relevant file ! I will check it.

>> I can only guess that this approach allow for simultaneous vdw
>> coupling/decoupling of the different atoms
>> like in our case in second stage when benzene H6 is decoupled and phenol
>> O1,H6 atoms are coupled.
> yes, that is true, you can simultaneously have atoms appear and
> disappear.
>> However this should be also eventually possible to divide into 2
>> separate
>> stages ((i)only H6 disappearing and then (ii) only O1,H6 appearing)
>> with separate dG contributions which together should get the dG of the
>> simultaneous H6 disappearing and O1,H6 appearing or not ?
> You could do that if you wanted with the currect code. However, I do not
> see why you would want to break it down into these two steps. Removing
> and
> adding a groups at the same time seems a better way to do this for me,
> since it doesn't introduce yet another non-physical state quite different
> from the end states.

Yes I agree with you. It was just to check if I understood well whats
going on here...

>> Regarding to uncharging of benzene H6 atom stage or the last stage where
>> O1,H6 phenol atoms are charged I really have no
>> explanation for two simultaneous sander simulations during the given
>> lambda run.
> one prmtop has charges on the atoms the other hasn't. Each is represented
> as its own groupsander process. This could have been done in a single
> sander simulation as well, but neither the simulation speed nor the
> ease-of-use would have benefited from such an option in my opinion.

I will for sure check the relevant article and the code as you advised but
for this
moment I still do not see the necessity of two simultaneous sander runs in
case on
just one way TI processes (i.e. only charging or only uncharging or only
vdw coupling or only vdw decoupling).

Here in principle just one PRMTOP should be enough if I am able to
provide the rest information
(definition of "unique" atoms and type of the TI change so to define V0
and V1 state) as a part of *.in file. For this
just one mask (like SCMASK) and another parameter which define one of the
4 possible TI changes
should be enough. Is still mystery for me what the second sander thread is
calculating in such one way cases
but maybe after reading that article and checking that *.f file everything
will be clear ( hopefully :)) ).

>> #3
>> I also do not understand why do you have three columns of the graphs
>> here
>> ( )
>> labeled as Step X - Complex, Step X - Water, Step X - Complex-Water
>> where
>> X is going from 1 to 3.
> Well, as the column labels say, one gives the values for the
> transformation in complex, one in water and the last one the difference
> between both.
>> #4
>> Do you think that TI might be possible to use also for the calculation
>> of
>> the absolute free energy of binding
>> of the small lignad (L) to the surface of bigger receptor (R) for
>> example
>> using this scenario.
>> L bound to R in given water-box.
>> L and the R separated in the same water-box sufficiently.
>> (Of course assuming the identical positions of all atoms except L in the
>> initial files for STARTING and FINAL system)
>> So I assume that during the whole TI procedure divided into relevant
>> decoupling/coupling (El.,vdw.) stages
>> L continuously "disappear" from it's original position and "appear" on
>> it's final position
>> sufficiently far from the R in the same water box.
>> Might be this approach OK in your opinion ? If not what is wrong ? How
>> to
>> improve it if TI could be really useful here ?
> While you could do things this way, there is extensive literature on
> doing
> absolute dG-TI (cited e.g. in our papers). The principle is simple, but
> there are numerous pitfalls and details that people have worked out, so
> make sure you consult these before attempting to setup a simulation.

OK thanks

>> BTW do you recommend in such case which I described above to use TI
>> (however maybe with much more complicated protocol than I suggested)
>> or is better and also more common to use here Umbrella sampling ?
> Depending on what your reaction coordinate for the dissociation looks
> like, US may be a better option than TI to obtain the binding energetics.
> I can't say anything definite here, since it really is very system and
> structure dependent.

OK, since I am interested about the small L bound weakly to the surface of
R (unfortunately not too much flat
considering dimensions of L) the reaction coordinate should be some kind
of distance L from R but that.s why I like the idea of TI
where one absolutely don.t care about any physical
dissociation/association pathway instead just clearly define
initial and the final state of the system and then transform the initial
to the final state using sequence of several elementary
"appear"/"disappear" processes.

   Thanks again for all !

       Best wishes,


> Kind Regards,
> Thomas
> Dr. Thomas Steinbrecher
> formerly at the
> BioMaps Institute
> Rutgers University
> 610 Taylor Rd.
> Piscataway, NJ 08854
> _______________________________________________
> AMBER mailing list
> __________ Informace od ESET NOD32 Antivirus, verze databaze 6684
> (20111205) __________
> Tuto zpravu proveril ESET NOD32 Antivirus.

Tato zpráva byla vytvořena převratným poštovním klientem Opery:
AMBER mailing list
Received on Mon Dec 05 2011 - 16:00:03 PST
Custom Search