[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:35:08 CEST 2017


Fair enough. Specifically I want the energy difference between lambda steps in order to incorporate the adaptive integration method (AIM) inside of GROMACS.

http://redmine.gromacs.org/issues/1849


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:

exp(-de + df)

where I’m calculating df as:

                for (j=0; j < efptNR; j++)
                {
                    if (fep->separate_dvdl[j])
                    {
                        // integrate from fep_state to lamtrial
                        df += 0.5
                           * (fep->all_lambda[j][lamtrial] - fep->all_lambda[j][fep_state])
                           * (dfhist->sum_dvdl[j][lamtrial] + dfhist->sum_dvdl[j][fep_state])
                           * (1.0/expand->mc_temp*BOLTZ);
                    }
                }

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.

Is that better?

Cheers,

-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:23 PM, Michael R Shirts <Michael.Shirts at Colorado.EDU<mailto:Michael.Shirts at Colorado.EDU>> wrote:

>  I actually want to do calculations between lambda steps inside of GROMACS code, before the end of the simulation.

“calculations” is a bit vague, as is “between lambda steps”.  Can you be more specific?

Best,
~~~~~~~~~~~~~~~~
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:19 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: Re: [gmx-developers] energy difference for lambda in free energy calculations

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>.

--
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/b722c041/attachment-0001.html>


More information about the gromacs.org_gmx-developers mailing list