- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: James W <c00jsw00.nchc.org.tw>

Date: Wed, 2 May 2007 17:44:28 +0800

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
*

Date: Wed, 2 May 2007 17:44:28 +0800

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

============================================================================

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

Custom Search