[gmx-users] My FEP so far (1)

[Sander Pronk]

Hi,

The procedure you're describing is correct (and I believe g_analyze can do the integration for you). However, with 4.5 there's a new way to calculate free energies, using Bennett' s Acceptance Ratio. 

That method relies on the energy differences between the state the simulation is at (init_lambda in this case) and the state you want to calculate the free energy difference for, as a set of instantaneous values, written to 'dhdl.xvg'. You can enable calculating the energies w.r.t.  those states with the 'foreign lambda' mdp parameter. For example, if you'd like 6 lambda points between 0 and 1, your mdp parameters would look like this:

init_lambda		foreign_lambda
0.0				       0.2
0.2				0.0 0.4
0.4				0.2 0.6
0.6				0.4 0.8
0.8				0.6 1.0
1.0				0.8

(so for each simulation at a specific lambda point there are two foreign lambdas, one with the 'previous' lambda value and one with the 'next' lambda value). BTW, 6 lambda points is probably not enough.

You could also set all foreign_lambda values to all lambda values you're simulating at: '0.0 0.2 0.4 0.6 0.8 1.0' for all simulations - the method would work but you'll lose some disk space (and maybe some performance when there are a lot of lambda points). 

You can then use g_bar to calculate the full free energy difference with error estimates (for example, g_bar -f lambda=*/dhdl.xvg) . This method has much better statistical (and numerical) properties than integrating dH/dl.

I hope all this makes sense. We're working on getting an updated free energy tutorial - I'll announce it here when we're done.

Sander


On Sep 28, 2010, at 02:47 , TJ Mustard wrote:

> Hello all,
>  
> As there is very little out on the web for FEP in gromacs that is up to date I thought I would throw my two cents out there for approval of my techniques (hopefully) and to help others that find this subject hard to get a significant amount of information on. Of course there are papers to be read but the specific settings and how to go about them are
>  
> First if you are using 4.5+ I found that it has arguments that can be set in its mdp files:
>  
> ; Free energy control stuff
> free-energy              = yes  ; = no
> init-lambda              = 0.00         ; = 0
> delta-lambda             = 0    ; = 0
> foreign_lambda           =      ; =
> sc-alpha                 = 0.5  ; = 0
> sc-power                 = 1.0  ; = 0
> sc-sigma                 = 0.3  ; = 0.3
> nstdhdl                  = 10   ; = 10
> separate-dhdl-file       = yes  ; = yes
> dhdl-derivatives         = yes  ; = yes
> dh_hist_size             = 0    ; = 0
> dh_hist_spacing          = 0.1  ; = 0.1
> couple-moltype           = RNA_chain_A  ; =
> couple-lambda0           = vdw-q        ; = vdw-q
> couple-lambda1           = none ; = vdw-q
> couple-intramol          = no   ;       = no
>  
>  
> On this run I set the RNA_chain_A to be perturbation. I also set the change to be all interactions to none. To keep artifacts to a minimum I set the "sc" settings as above. (found these in the FEP tutorial online here http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial) I would then make several identical runs with different init_lambda values. The md or sd log file would output the average dVpot/dLambda and I would put these values in a spreadsheet vs their respective lambda values. I would then use Simpson's Rule or other to approximate the area under the curve. With this technique and lambda steps of 0.05 I calculated the FEP for the hydration of methane to be ~2.5 kcal/mol. Not perfect but in the ballpark according to the tutorial above.
>  
> I found that the addition of this argument set in 4.5+ allowed for easier FEP. (ie not having to manually edit the state b for each atom in question.)
>  
> If anyone else has anything so add or correct please do so. I found this very difficult to understand as the tutorial is for a different force field and my future systems were going to be much larger.
>  
> Thank you,
>  
> TJ Mustard
> Email: mustardt at onid.orst.edu
> 
> -- 
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100928/c5b07972/attachment.html>