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

[Mirabzadeh, Christopher (mira2978@vandals.uidaho.edu)]

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



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