[gmx-developers] energy difference for lambda in free energy calculations

Michael R Shirts Michael.Shirts at Colorado.EDU
Fri Apr 21 01:15:07 CEST 2017

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>
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> on behalf of "Mirabzadeh, Christopher (mira2978 at vandals.uidaho.edu)" <mira2978 at vandals.uidaho.edu>
Reply-To: "gmx-developers at gromacs.org" <gmx-developers at gromacs.org>
Date: Thursday, April 20, 2017 at 5:01 PM
To: "gmx-developers at gromacs.org" <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.



Mirabzadeh, Christopher
Graduate Research Assistant/Physics Instructor
Department of Physics
University of Idaho
Moscow, Id

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170420/a14f9608/attachment-0001.html>

More information about the gromacs.org_gmx-developers mailing list