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
Received on Sun May 06 2007 - 06:07:07 PDT
