Subject: Re: CCL:Thermodynamic integration

From: <[email-masked]>
Date: Tue, 11 Mar 1997 11:53:05 +0100 (NFT)

Hi, there:

> I am trying to perform a molecular dynamics calculation of the free 
> energy for the mutation of a water molecule into a methane molecule 
> in solution. I am using the thermodynamic integration technique. 
> Although I have been trying for quite a while, gradually finding 
> the bugs in my code, I have still not been able to reproduce the 
> literature values for this reaction.
> 
> I am starting to suspect that I may not have a correct understanding
> of how this technique is used. While the many references I have 
> consulted talk of 'accumulating the partial derivative dU/dlambda'
> they do not describe exactly what this partial derviative is in 
> terms of the positions, velocites etc that can be accumulated each 
> time step.

First thing: what you want to calculate is:

                        1
      delta G = integral  < dU(x,lambda) / dlambda >       dlambda
                        0                           lambda

In the simulation, the integrand comes from a discrete number of points
and your integral is approximated via a standard rule, such as Simpson,
or trapezoidal:
                   N
      delta G ~ sum  { < dU(x,lambda) / dlambda >            +
                   0                             lambda(k)  
                       < dU(x,lambda) / dlambda >            }  x
                                                 lambda(k+1)
                     { lambda(k+1) - lambda(k) } / 2

The choice of N, the number of windows, can be critical for `mutation'-
type free energy calculation...  My experience is that it is often more 
efficient to use few windows (say less than 30) with extensive sampling
than the contrary. Be sure, though, that the integrand varies smoothly.

Now, for that kind of simulation,  I assume that you `grow' bonds, con-
comittently with increasing van der Waals radii. This deserves particu-
lar attention, especially if bonds are shaken. What you should evaluate
here is the contribution to the free energy of the constraints changing
the bond lengths from lambda = 0 to 1. Simply, the latter can be calcu-
lated from:

      dUconstraint / dlambda = {(Fij . rij) / |rij|} x (drij / dlambda)

where Fij is the force dur to the constraint and rij is the constrained
bond length.

One last,  trivial thing.  For the above mutation, I guess you decouple
the process into subprocesses, dedicated to electrostatic and internal+
van der Waals parameters... It is always safer to do so, especially for 
`absolute' free energies, than mutating everything simultaneously - un-
less species at lambda = 0 and 1 are topologically identical (say K+ ->
Na+).

Hope this helps...

Chris
_______________________________________________________________________

Chris Chipot                         
Laboratoire de Chimie Theorique          Phone:      (33) 3-83-91-25-96
U.R.A. C.N.R.S. No 510                   Fax:        (33) 3-83-91-25-30
Universite Henri Poincare - Nancy I       
B.P. 239
54506 Vandoeuvre-les-Nancy Cedex         E-mail: [email-masked]
_______________________________________________________________________