<html xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Title" content="">
<meta name="Keywords" content="">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Wingdings;
panose-1:5 0 0 0 0 0 0 0 0 0;}
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:12.0pt;
font-family:"Times New Roman";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
span.apple-converted-space
{mso-style-name:apple-converted-space;}
span.EmailStyle18
{mso-style-type:personal-reply;
font-family:Calibri;
color:windowtext;}
span.msoIns
{mso-style-type:export-only;
mso-style-name:"";
text-decoration:underline;
color:teal;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style>
</head>
<body bgcolor="white" lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri">So, the updated data can be seen most clearly in mdebin.c, in the function upd_mdebin, right after:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri"> <b>if</b> ((md->fp_dhdl || md->dhc) && bDoDHDL)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri">The components for dH/dL are in enerd->term[F_DVDL+i]), where i=0,1,2,3,4, and you can see by reading the dhdl output what the order is.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri">The potential energy differences are calculated from: from enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0]<o:p></o:p></span></p>
<div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">If you are only interested in the neighboring differences, the can be calculated from enerd->enerpart_lambda[i+1]- enerd->enerpart_lambda[i] (as the [0] component cancels out).<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">So the question is where in the code are they “in step” with the timestep?<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">When one gets to ExpandedEnsembleDynamics in md.c, then they are all updated. However, if you are using velocity verlet, then they are not finished updating until after sum_dhdl(enerd,
state->lambda, ir->fepvals) is called. If you only acess these arrays within the ExpandedEnsembleDynamics (which it sounds like you would be doing, since AIM deals with transitions between states), then you should be fine.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Is this helpful?<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Best,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">~~~~~~~~~~~~~~~~<o:p></o:p></span></p>
<div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Michael Shirts<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Associate Professor<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><a href="mailto:michael.shirts@colorado.edu">michael.shirts@colorado.edu</a><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black"><a href="http://www.colorado.edu/lab/shirtsgroup/"><span style="text-decoration:none">http://www.colorado.edu/lab/shirtsgroup/</span></a><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Phone: (303) 735-7860<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Office: JSCBB C123<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">Department of Chemical and Biological Engineering<o:p></o:p></span></p>
</div>
</div>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri;color:black">University of Colorado Boulder</span><span style="font-size:11.0pt;font-family:Calibri"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;font-family:Calibri"><o:p> </o:p></span></p>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b><span style="font-family:Calibri;color:black">From: </span>
</b><span style="font-family:Calibri;color:black"><gromacs.org_gmx-developers-bounces@maillist.sys.kth.se> on behalf of "Mirabzadeh, Christopher (mira2978@vandals.uidaho.edu)" <mira2978@vandals.uidaho.edu><br>
<b>Reply-To: </b>"gmx-developers@gromacs.org" <gmx-developers@gromacs.org><br>
<b>Date: </b>Thursday, April 20, 2017 at 5:35 PM<br>
<b>To: </b>"gmx-developers@gromacs.org" <gmx-developers@gromacs.org><br>
<b>Subject: </b>Re: [gmx-developers] energy difference for lambda in free energy calculations<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal">Fair enough. Specifically I want the energy difference between lambda steps in order to incorporate the adaptive integration method (AIM) inside of GROMACS.
<o:p></o:p></p>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><a href="http://redmine.gromacs.org/issues/1849">http://redmine.gromacs.org/issues/1849</a><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">I need the de term as part of the acceptance criteria for the Monte Carlo moves between lambda steps. I’ve been working on this for a while now and thought I had things working correctly last year until I dissected the alchemical-analysis
code and found that I wasn’t properly accounting for the individual delta lambda terms. AIM uses the same integral as TI in it’s acceptance criteria. AIM also uses the difference between neighboring lambda simulations, thus the acceptance criteria is:<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">exp(-de + df) <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">where I’m calculating df as:<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<div>
<p class="MsoNormal"> for (j=0; j < efptNR; j++)<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> {<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> if (fep->separate_dvdl[j])<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> {<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> // integrate from fep_state to lamtrial<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> df += 0.5<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> * (fep->all_lambda[j][lamtrial] - fep->all_lambda[j][fep_state])<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> * (dfhist->sum_dvdl[j][lamtrial] + dfhist->sum_dvdl[j][fep_state])<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> * (1.0/expand->mc_temp*BOLTZ);<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> }<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"> }<o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Based on comparisons with the output of TI, I believe I’m calculating df correctly but my histogram isn’t flat and my estimate for the free energy is low. I’m trying to make sure I’m using the correct value for de within my acceptance criteria.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Is that better?<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">Cheers,<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<div>
<p class="MsoNormal">-ChrisM<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt">Mirabzadeh, Christopher <br>
Graduate Research Assistant/Physics Instructor<br>
Department of Physics<br>
University of Idaho<br>
Moscow, Id<br>
(509)339-5647<br>
<br>
<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<p class="MsoNormal">On Apr 20, 2017, at 4:23 PM, Michael R Shirts <<a href="mailto:Michael.Shirts@Colorado.EDU">Michael.Shirts@Colorado.EDU</a>> wrote:<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div style="margin-left:.5in">
<p class="MsoNormal" style="text-indent:-.25in;background:white"><span style="font-family:Wingdings">Ø</span><span style="font-size:7.0pt"> <span class="apple-converted-space"> </span></span>I actually want to do calculations between lambda steps inside of
GROMACS code, before the end of the simulation.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">“calculations” is a bit vague, as is “between lambda steps”. Can you be more specific?</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Best,</span><o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">~~~~~~~~~~~~~~~~</span><o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Michael Shirts</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Associate Professor</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"><a href="mailto:michael.shirts@colorado.edu"><span style="color:purple">michael.shirts@colorado.edu</span></a></span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"><a href="http://www.colorado.edu/lab/shirtsgroup/"><span style="color:purple;text-decoration:none">http://www.colorado.edu/lab/shirtsgroup/</span></a></span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Phone: (303) 735-7860</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Office: JSCBB C123</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Department of Chemical and Biological Engineering</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">University of Colorado Boulder</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<div>
<p class="MsoNormal" style="background:white"><b><span style="font-family:Calibri">From:<span class="apple-converted-space"> </span></span></b><span style="font-family:Calibri"><<a href="mailto:gromacs.org_gmx-developers-bounces@maillist.sys.kth.se">gromacs.org_gmx-developers-bounces@maillist.sys.kth.se</a>>
on behalf of "Mirabzadeh, Christopher (<a href="mailto:mira2978@vandals.uidaho.edu">mira2978@vandals.uidaho.edu</a>)" <<a href="mailto:mira2978@vandals.uidaho.edu">mira2978@vandals.uidaho.edu</a>><br>
<b>Reply-To:<span class="apple-converted-space"> </span></b>"<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>" <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
<b>Date:<span class="apple-converted-space"> </span></b>Thursday, April 20, 2017 at 5:19 PM<br>
<b>To:<span class="apple-converted-space"> </span></b>"<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>" <<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a>><br>
<b>Subject:<span class="apple-converted-space"> </span></b>Re: [gmx-developers] energy difference for lambda in free energy calculations</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal" style="background:white">Hello. Thanks!<span class="apple-converted-space"> </span><o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white">I actually want to do calculations between lambda steps inside of GROMACS code, before the end of the simulation.<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white">-ChrisM<o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt;background:white">Mirabzadeh, Christopher <br>
Graduate Research Assistant/Physics Instructor<br>
Department of Physics<br>
University of Idaho<br>
Moscow, Id<br>
(509)339-5647<br>
<br>
<br>
<o:p></o:p></p>
</div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
<div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<p class="MsoNormal" style="background:white">On Apr 20, 2017, at 4:14 PM, Michael R Shirts <<a href="mailto:Michael.Shirts@Colorado.EDU"><span style="color:purple">Michael.Shirts@Colorado.EDU</span></a>> wrote:<o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Hi, Chris-</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">So, the dhdl.xvg output prints the potential energy difference between neighboring simulations (for BAR, MBAR, etc.), and the dh/dl values (for TI). So that file
should have all the information you need for any method – you can see the<span class="apple-converted-space"> </span><a href="https://github.com/MobleyLab/alchemical-analysis/"><span style="color:#954F72">https://github.com/MobleyLab/alchemical-analysis/</span></a><span class="apple-converted-space"> </span>scripts
for how this information is used.</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">The ROUTE to getting the values right is rather complicated for a number of reasons. So before answering in more detail-</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
</div>
<div style="margin-left:.5in">
<div>
<p class="MsoNormal" style="text-indent:-.25in;background:white"><span style="font-size:11.0pt;font-family:Calibri">1.</span><span style="font-size:7.0pt"> <span class="apple-converted-space"> </span></span><span style="font-size:11.0pt;font-family:Calibri">If
you want to just USE this information, it’s all output by GROMACS already.</span><o:p></o:p></p>
</div>
</div>
<div style="margin-left:.5in">
<div>
<p class="MsoNormal" style="text-indent:-.25in;background:white"><span style="font-size:11.0pt;font-family:Calibri">2.</span><span style="font-size:7.0pt"> <span class="apple-converted-space"> </span></span><span style="font-size:11.0pt;font-family:Calibri">If
you want to MODIFY what GROMACS is outputting, can you be more specific what you want to do?</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">~~~~~~~~~~~~~~~~</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Michael Shirts</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Associate Professor</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"><a href="mailto:michael.shirts@colorado.edu"><span style="color:purple">michael.shirts@colorado.edu</span></a></span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"><a href="http://www.colorado.edu/lab/shirtsgroup/"><span style="color:purple;text-decoration:none">http://www.colorado.edu/lab/shirtsgroup/</span></a></span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Phone: (303) 735-7860</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Office: JSCBB C123</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">Department of Chemical and Biological Engineering</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri">University of Colorado Boulder</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:11.0pt;font-family:Calibri"> </span><o:p></o:p></p>
</div>
</div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<div>
<div>
<p class="MsoNormal" style="background:white"><b><span style="font-family:Calibri">From:<span class="apple-converted-space"> </span></span></b><span style="font-family:Calibri"><<a href="mailto:gromacs.org_gmx-developers-bounces@maillist.sys.kth.se"><span style="color:#954F72">gromacs.org_gmx-developers-bounces@maillist.sys.kth.se</span></a>>
on behalf of "Mirabzadeh, Christopher (<a href="mailto:mira2978@vandals.uidaho.edu"><span style="color:#954F72">mira2978@vandals.uidaho.edu</span></a>)" <<a href="mailto:mira2978@vandals.uidaho.edu"><span style="color:#954F72">mira2978@vandals.uidaho.edu</span></a>><br>
<b>Reply-To:<span class="apple-converted-space"> </span></b>"<a href="mailto:gmx-developers@gromacs.org"><span style="color:#954F72">gmx-developers@gromacs.org</span></a>" <<a href="mailto:gmx-developers@gromacs.org"><span style="color:#954F72">gmx-developers@gromacs.org</span></a>><br>
<b>Date:<span class="apple-converted-space"> </span></b>Thursday, April 20, 2017 at 5:01 PM<br>
<b>To:<span class="apple-converted-space"> </span></b>"<a href="mailto:gmx-developers@gromacs.org"><span style="color:#954F72">gmx-developers@gromacs.org</span></a>" <<a href="mailto:gmx-developers@gromacs.org"><span style="color:#954F72">gmx-developers@gromacs.org</span></a>><br>
<b>Subject:<span class="apple-converted-space"> </span></b>[gmx-developers] energy difference for lambda in free energy calculations</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white">Hello all,<span class="apple-converted-space"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">What I’m trying to do is get the energy difference between neighboring expanded ensemble simulations, i.e U(lambda_trial, x) - U(lambda_current, x). As such, I am trying to better understand the de term used in
expanded ensemble calculations (de = weighted_lamee[lamtrial] - weighted_lamee[fep_state]). <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">It is defined in forcerec.h as:<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">typedef struct {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> real term[F_NRE]; /* The energies for all different interaction types */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> gmx_grppairener_t grpp;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> double dvdl_lin[efptNR]; /* Contributions to dvdl with linear lam-dependence */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> double dvdl_nonlin[efptNR]; /* Idem, but non-linear dependence */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> int n_lambda;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> int fep_state; /*current fep state -- just for printing */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"><span style="color:#FF2600"> double *enerpart_lambda; /* Partial energy for lambda and flambda[] */</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> real foreign_term[F_NRE]; /* alternate array for storing foreign lambda energies */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">} gmx_enerdata_t;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">/* The idea is that dvdl terms with linear lambda dependence will be added<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> * automatically to enerpart_lambda. Terms with non-linear lambda dependence<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> * should explicitly determine the energies at foreign lambda points<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> * when n_lambda > 0.<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">I see that this term is used in the xvg output inside of mdebin.c and <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> /* BAR + thermodynamic integration values */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> if ((md->fp_dhdl || md->dhc) && bDoDHDL)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> for (i = 0; i < enerd->n_lambda-1; i++)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> /* zero for simulated tempering */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> md->dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">……<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">Also, in the replica exchange code (repl_ex.c) we find this term again;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> if (re->type == ereLAMBDA || re->type == ereTL)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> bDLambda = TRUE;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> /* lambda differences. */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> minus the energy of the jth simulation in the jth Hamiltonian */<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> for (i = 0; i < re->nrepl; i++)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> for (j = 0; j < re->nrepl; j++)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> re->de[i][j] = 0;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> }<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> }<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> for (i = 0; i < re->nrepl; i++)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> re->de[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> ……<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">I’m trying to understand how this term is used in Thermodynamic integration (Pymbar for instance). Also, I’ve noticed that the enerd->enerpart_lambda[0] term is only updated every nstenergy step. Could anyone explain
what the 0 term is and how this enerpart_lambda[i]-enerpart_lambda[0] term would be used in TI? <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">The bottom line is that I need the correct value for the energy difference between neighboring lambda simulations. Does the de term do that for me? <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">If de is the correct term, is delta lambda part of the calculation? The reason I ask this question is from not understanding the output of debug information in force.c, specifically;<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> if (debug)<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> {<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> fepvals->all_lambda[j][i], efpt_names[j],<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> dlam, enerd->dvdl_lin[j]);<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> }<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">Does this mean that the energy difference is (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]) + dlam*dvdl_lin[i]? I’m very confused.<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">Thanks,<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white">-ChrisM<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<p class="MsoNormal" style="margin-bottom:12.0pt;background:white;background-position:initial initial;background-repeat:initial initial">
Mirabzadeh, Christopher <br>
Graduate Research Assistant/Physics Instructor<br>
Department of Physics<br>
University of Idaho<br>
Moscow, Id<br>
(509)339-5647<br>
<br>
<br>
<br>
<o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<p class="MsoNormal" style="background:white"><span style="font-size:9.0pt;font-family:Helvetica;background:white">--<span class="apple-converted-space"> </span></span><span style="font-size:9.0pt;font-family:Helvetica"><br>
<span style="background:white">Gromacs Developers mailing list</span><br>
<br>
<span style="background:white">* Please search the archive at<span class="apple-converted-space"> </span></span></span><a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List"><span style="font-size:9.0pt;font-family:Helvetica;color:#954F72;background:white">http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</span></a><span class="apple-converted-space"><span style="font-size:9.0pt;font-family:Helvetica;background:white"> </span></span><span style="font-size:9.0pt;font-family:Helvetica;background:white">before
posting!</span><span style="font-size:9.0pt;font-family:Helvetica"><br>
<br>
<span style="background:white">* Can't post? Read<span class="apple-converted-space"> </span></span></span><a href="http://www.gromacs.org/Support/Mailing_Lists"><span style="font-size:9.0pt;font-family:Helvetica;color:#954F72;background:white">http://www.gromacs.org/Support/Mailing_Lists</span></a><span style="font-size:9.0pt;font-family:Helvetica"><br>
<br>
<span style="background:white">* For (un)subscribe requests visit</span><br>
</span><a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers"><span style="font-size:9.0pt;font-family:Helvetica;color:#954F72;background:white">https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</span></a><span class="apple-converted-space"><span style="font-size:9.0pt;font-family:Helvetica;background:white"> </span></span><span style="font-size:9.0pt;font-family:Helvetica;background:white">or
send a mail to<span class="apple-converted-space"> </span></span><a href="mailto:gmx-developers-request@gromacs.org"><span style="font-size:9.0pt;font-family:Helvetica;color:#954F72;background:white">gmx-developers-request@gromacs.org</span></a><span style="font-size:9.0pt;font-family:Helvetica;background:white">.</span><o:p></o:p></p>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal" style="background:white"> <o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><span style="font-size:9.0pt;font-family:Helvetica;background:white">--<span class="apple-converted-space"> </span></span><span style="font-size:9.0pt;font-family:Helvetica"><br>
<span style="background:white">Gromacs Developers mailing list</span><br>
<br>
<span style="background:white">* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List">
http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</a> before posting!</span><br>
<br>
<span style="background:white">* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists">
http://www.gromacs.org/Support/Mailing_Lists</a></span><br>
<br>
<span style="background:white">* For (un)subscribe requests visit</span><br>
<span style="background:white"><a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers">https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</a> or send a mail to
<a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.</span></span><o:p></o:p></p>
</div>
</blockquote>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</div>
</body>
</html>