<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hi,<div><br></div><div>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. </div><div><br></div><div>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:</div><div><br></div><div>init_lambda<span class="Apple-tab-span" style="white-space:pre"> </span>foreign_lambda</div><div>0.0<span class="Apple-tab-span" style="white-space:pre"> </span> 0.2</div><div>0.2<span class="Apple-tab-span" style="white-space:pre"> </span>0.0 0.4</div><div>0.4<span class="Apple-tab-span" style="white-space:pre"> </span>0.2 0.6</div><div>0.6<span class="Apple-tab-span" style="white-space:pre"> </span>0.4 0.8</div><div>0.8<span class="Apple-tab-span" style="white-space:pre"> </span>0.6 1.0</div><div>1.0<span class="Apple-tab-span" style="white-space:pre"> </span>0.8</div><div><br></div><div>(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.</div><div><br></div><div>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). </div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>Sander</div><div><br></div><div><br><div><div>On Sep 28, 2010, at 02:47 , TJ Mustard wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite">
<div><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Hello all,</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">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</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">First if you are using 4.5+ I found that it has arguments that can be set in its mdp files:</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">; Free energy control stuff<br>
free-energy = yes ; = no<br>
init-lambda = 0.00 ; = 0<br>
delta-lambda = 0 ; = 0<br>
foreign_lambda = ; =<br>
sc-alpha = 0.5 ; = 0<br>
sc-power = 1.0 ; = 0<br>
sc-sigma = 0.3 ; = 0.3<br>
nstdhdl = 10 ; = 10<br>
separate-dhdl-file = yes ; = yes<br>
dhdl-derivatives = yes ; = yes<br>
dh_hist_size = 0 ; = 0<br>
dh_hist_spacing = 0.1 ; = 0.1<br>
couple-moltype = RNA_chain_A ; =<br>
couple-lambda0 = vdw-q ; = vdw-q<br>
couple-lambda1 = none ; = vdw-q<br>
couple-intramol = no ; = no<span></span></div><p style="margin: 0px;"> </p><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">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 <a href="http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial">http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial</a>) 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.</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">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.)</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">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.</div><p style="margin: 0px;"> </p><div style="margin-top: 0px; margin-right: 0px; margin-bottom: 0px; margin-left: 0px; ">Thank you,</div><p style="margin: 0px;"><span></span></p><p style="margin: 0px;"> </p><p style="font-family: monospace; white-space: nowrap; margin: 5px 0px 5px 0px;">TJ Mustard<br>
Email: <a href="mailto:mustardt@onid.orst.edu">mustardt@onid.orst.edu</a></p>
</div>
-- <br>gmx-users mailing list <a href="mailto:gmx-users@gromacs.org">gmx-users@gromacs.org</a><br><a href="http://lists.gromacs.org/mailman/listinfo/gmx-users">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!<br>Please don't post (un)subscribe requests to the list. Use the <br>www interface or send it to gmx-users-request@gromacs.org.<br>Can't post? Read http://www.gromacs.org/Support/Mailing_Lists</blockquote></div><br></div></body></html>