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