Re: AMBER: Steered MD & Jarzynski's equality

From: wang <c00jsw00.nchc.org.tw>
Date: Thu, 03 May 2007 16:06:43 +0800

Dear all ,
I want to use distance constraint ( Ca atom of 121th residue and Ca of
251th residue ) .
Could you tell me how to do that ? I tried NOE module . But I still
didn't success .
thank you
Wang

David Mobley 提到:

> It's probably also worth pointing out at this point that the Jarzynski
> equality requires (on point 4) a large number (n) of simulations
> because the most important values for the average are those at the
> tails of the distribution. How large is hard to say -- but do some
> testing at the beginning to see whether this is a price you can afford
> to pay!
>
> (See the recent Jarzynski paper on convergence of exponentially
> averaged work values). The Crooks approach is more efficient, but
> requires doing both forward and reverse simulations for every
> transformation. Not sure whether this is easy (possible?) to do with
> the current code, as I'm not an SMD guy... You'd need fewer
> realizations then.
>
> David
>
>
> On 5/2/07, Gustavo Seabra <gustavo.seabra.gmail.com> wrote:
>
>> Hi James,
>>
>> What your code is doing, it seems to me, is averaging over all work
>> values from *one* run, which has no meaning. To use the Jarzynski
>> relationship, what you want is to, *at each point*, average a large
>> number of work values. So, what you would do is:
>>
>> 1. Equilibrate the system, using a restraint to keep the coordinate of
>> interest in the initial value. This will give you an equilibrated
>> initial ensemble. (Now, if only we could define "equilibrate"... :-) )
>>
>> 2. From that equilibrated ensemble, pick a large number of structures
>> (let's call it 'n'). These will be your initial structures for the next
>> step, the pulling using the smd module.
>>
>> 3. (And this is the part you seem to be missing.) For each one of the
>> 'n' structures you've picked, run an *independent* smd simulation,
>> slowly pulling the coordinate from the initial to the final value.
>> Notice that this will give you 'n' dist_vs_t files. Now the fun
>> begins...
>>
>> 4. At every value of the reaction coordinate (1st column), average the
>> work values from all the 'n' different simulations, using the Jarzynski
>> formula. This will give you a (relative) free energy value at every
>> value of the reaction coordinate, your "PMF" if you will.
>>
>> I hope that makes it a bit clearer for you.
>>
>> Good luck,
>>
>> Gustavo.
>>
>>
>>
>> James W wrote:
>> > Dear Prof. Roitberg ,
>> > I used the amber smd module to run a simulation . This is my output :
>> > ____________________________________________________________________
>> > x0(t) x force work
>> > 57.00000 57.73576 -1471.51842 0.00000
>> > 57.00093 57.34648 -691.10171 -1.00562
>> > .....
>> > ____________________________________________________________________
>> > And this is my code for com[uting pmf :
>> > ____________________________________________________________________
>> > program pmf
>> > real :: coor(100000),p(100000) ,work(100000)
>> > beta = 1.04652438
>> > open(unit=10,file='dist_vs_t')
>> > do i=1,100000,1
>> > read(10,*) a,b,c,d
>> > coor(i) = a
>> > work(i) = d
>> > end do
>> > ! PMF
>> > do i=1,100000,1
>> > pt = exp ( - work(i)/beta ) / i
>> > p(i) = - beta * log (pt)
>> > write(*,*) coor(i) , p(i)
>> > end do
>> > end program
>> > ____________________________________________________________________
>> > But the result is very curious , could you tell me that my code is
>> right ?
>> > thanks
>> >
>> > wang
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > On Wed, 02 May 2007 09:26:30 +0200, Adrian Roitberg wrote
>> >
>> >> James W wrote:
>> >>
>> >>> Dear all ,
>> >>> #the reference : (http://en.wikipedia.org/wiki/Jarzynski_equality)
>> >>> I confused the Jarzynski's equality . The equality is shown:
>> >>> exp (-F/KT) = < exp (-W/KT) > .
>> >>> &
>> >>> F = < W > + ......
>> >>> I used SMD module of CHARMM and obtained the reaction coordinates &
>> >>> work ,like this :
>> >>> ______________________________________________________
>> >>> 57.00000 57.73576 -1471.51842 0.00000
>> >>> 57.00093 57.34648 -691.10171 -1.00562
>> >>> (RC) (WORK)
>> >>> *RC : reaction coordinates
>> >>> I wanted to obtain the " < W > " form my data . My method is :
>> >>> 1. S = sum { exp (-W(i) /kT )} , i = 1 ....N
>> >>> 2. p(i)= exp (-W(i) /kT ) / S
>> >>> 3. <W(i)> = W(i)*p(i) / S
>> >>>
>> >>> Could you tell me that my method is right ?
>> >>>
>> >> James,
>> >> I am not sure I fully understand your question, because you do not
>> >> mention what the subindex i is in your formulas.
>> >> Now, the formula from wikipedia is basically right, except for a
>> couple
>> >> of comments.
>> >>
>> >> First, it must be clear that F is really \Delta F. This is
>> designed for
>> >> free energy differences and NOT absolute free energies. Second, the
>> >> exponential average of work values is over a set of runs, ALL
>> >> starting from an equilibrated ensemble at a certain value of the
>> >> coordinates (s) to be changed. The average then involves MANY runs.
>> >>
>> >> The second line F = < W > + ...... is a cumulant expansion of
>> >> exponential average. I do not recommend using using this unless you
>> >> really know what you are doing.
>> >>
>> >> Your formulas on how to implement this are not quite right.
>> >> If you are thinking about the index i meaning time steps on one
>> >> simulation, then it is not right.
>> >> The subindex i should refer to the same distance and DIFFERENT
>> simulations.
>> >>
>> >> Once that is done:
>> >>
>> >> 1. S = ( sum { exp (-W(i) /kT )}/N) , i = 1 ....N
>> >> 2. \Delta F = - kT ln(S)
>> >>
>> >> I hope this makes things a bit clearer.
>> >>
>> >> Adrian
>> >>
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to amber.scripps.edu
>> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>


-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun May 06 2007 - 06:07:15 PDT
Custom Search