[gmx-users] My FEP so far (1)
pronk at cbr.su.se
Tue Sep 28 09:40:52 CEST 2010
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:
0.2 0.0 0.4
0.4 0.2 0.6
0.6 0.4 0.8
0.8 0.6 1.0
(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.
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
> 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...
More information about the gromacs.org_gmx-users