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: musta...@onid.orst.edu > > -- > gmx-users mailing list gmx-users@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-requ...@gromacs.org. > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- gmx-users mailing list gmx-users@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-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
