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
>
> --
> Dr. Adrian E. Roitberg
> Associate Professor
> Quantum Theory Project and Department of Chemistry
>
> University of Florida PHONE 352 392-6972
> P.O. Box 118435 FAX 352 392-8722
> Gainesville, FL 32611-8435 Email adrian.qtp.ufl.edu
>
============================================================================
>
> To announce that there must be no criticism of the president,
> or that we are to stand by the president right or wrong,
> is not only unpatriotic and servile, but is morally treasonable
> to the American public."
> -- Theodore Roosevelt
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
--
National Center High-performance Computing (http://www.nchc.org.tw)
-----------------------------------------------------------------------
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:04 PDT
