[gmx-developers] energy difference for lambda in free energy calculations
Mirabzadeh, Christopher (mira2978@vandals.uidaho.edu)
mira2978 at vandals.uidaho.edu
Fri Apr 21 01:19:30 CEST 2017
Hello. Thanks!
I actually want to do calculations between lambda steps inside of GROMACS code, before the end of the simulation.
-ChrisM
Mirabzadeh, Christopher
Graduate Research Assistant/Physics Instructor
Department of Physics
University of Idaho
Moscow, Id
(509)339-5647
On Apr 20, 2017, at 4:14 PM, Michael R Shirts <Michael.Shirts at Colorado.EDU<mailto:Michael.Shirts at Colorado.EDU>> wrote:
Hi, Chris-
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 https://github.com/MobleyLab/alchemical-analysis/ scripts for how this information is used.
The ROUTE to getting the values right is rather complicated for a number of reasons. So before answering in more detail-
1. If you want to just USE this information, it’s all output by GROMACS already.
2. If you want to MODIFY what GROMACS is outputting, can you be more specific what you want to do?
~~~~~~~~~~~~~~~~
Michael Shirts
Associate Professor
michael.shirts at colorado.edu<mailto:michael.shirts at colorado.edu>
http://www.colorado.edu/lab/shirtsgroup/
Phone: (303) 735-7860
Office: JSCBB C123
Department of Chemical and Biological Engineering
University of Colorado Boulder
From: <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se<mailto:gromacs.org_gmx-developers-bounces at maillist.sys.kth.se>> on behalf of "Mirabzadeh, Christopher (mira2978 at vandals.uidaho.edu<mailto:mira2978 at vandals.uidaho.edu>)" <mira2978 at vandals.uidaho.edu<mailto:mira2978 at vandals.uidaho.edu>>
Reply-To: "gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>" <gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>>
Date: Thursday, April 20, 2017 at 5:01 PM
To: "gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>" <gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>>
Subject: [gmx-developers] energy difference for lambda in free energy calculations
Hello all,
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]).
It is defined in forcerec.h as:
typedef struct {
real term[F_NRE]; /* The energies for all different interaction types */
gmx_grppairener_t grpp;
double dvdl_lin[efptNR]; /* Contributions to dvdl with linear lam-dependence */
double dvdl_nonlin[efptNR]; /* Idem, but non-linear dependence */
int n_lambda;
int fep_state; /*current fep state -- just for printing */
double *enerpart_lambda; /* Partial energy for lambda and flambda[] */
real foreign_term[F_NRE]; /* alternate array for storing foreign lambda energies */
gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */
} gmx_enerdata_t;
/* The idea is that dvdl terms with linear lambda dependence will be added
* automatically to enerpart_lambda. Terms with non-linear lambda dependence
* should explicitly determine the energies at foreign lambda points
* when n_lambda > 0.
*/
I see that this term is used in the xvg output inside of mdebin.c and
/* BAR + thermodynamic integration values */
if ((md->fp_dhdl || md->dhc) && bDoDHDL)
{
for (i = 0; i < enerd->n_lambda-1; i++)
{
/* zero for simulated tempering */
md->dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];
……
Also, in the replica exchange code (repl_ex.c) we find this term again;
if (re->type == ereLAMBDA || re->type == ereTL)
{
bDLambda = TRUE;
/* lambda differences. */
/* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
minus the energy of the jth simulation in the jth Hamiltonian */
for (i = 0; i < re->nrepl; i++)
{
for (j = 0; j < re->nrepl; j++)
{
re->de[i][j] = 0;
}
}
for (i = 0; i < re->nrepl; i++)
{
re->de[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);
……
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?
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?
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;
if (debug)
{
fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
fepvals->all_lambda[j][i], efpt_names[j],
(enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),
dlam, enerd->dvdl_lin[j]);
}
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.
Thanks,
-ChrisM
Mirabzadeh, Christopher
Graduate Research Assistant/Physics Instructor
Department of Physics
University of Idaho
Moscow, Id
(509)339-5647
--
Gromacs Developers mailing list
* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers or send a mail to gmx-developers-request at gromacs.org<mailto:gmx-developers-request at gromacs.org>.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170420/7edff23e/attachment-0001.html>
More information about the gromacs.org_gmx-developers
mailing list