Restarting from checkpoint no longer reinitializes WL weights.

Fixes a problem where mdrun was reinitializing the initial Wang-Landau
delta for expanded ensemble simulations, because the flag turning it
off was stored the expanded ensemble data structure (not saved in cpt)
instead of the df_history structure (is saved in checkpoint).  In the
process, some moderate encasulaton of the df_history structure and
the expanded ensemble methods.

Fixes #1350

Change-Id: I13492a7a9773fcb417fcd0ee106d851d9838ce25

diff --git a/include/mdrun.h b/include/mdrun.h
index a11fa78ac34145c2bdf2f912be8dc3e1e9e941e7..98ff760bc6b2ce821f5e598523a26bd24fa8d006 100644 (file)
--- a/include/mdrun.h
+++ b/include/mdrun.h
@@ -177,15 +177,18 @@ gmx_integrator_t do_tpi;
 
 void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit);
 
+GMX_LIBMD_EXPORT
+void init_expanded_ensemble(gmx_bool bStateFromCP, t_inputrec *ir, gmx_rng_t *mcrng, df_history_t *dfhist);
+
 GMX_LIBMD_EXPORT
 int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *enerd,
-                             t_state *state, t_extmass *MassQ, df_history_t *dfhist,
+                             t_state *state, t_extmass *MassQ, int fep_state, df_history_t *dfhist,
                              gmx_large_int_t step, gmx_rng_t mcrng,
                              rvec *v, t_mdatoms *mdatoms);
 
 GMX_LIBMD_EXPORT
 void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *expand, t_simtemp *simtemp, df_history_t *dfhist,
-                               int nlam, int frequency, gmx_large_int_t step);
+                               int fep_state, int frequency, gmx_large_int_t step);
 
 GMX_LIBMD_EXPORT
 void get_mc_state(gmx_rng_t rng, t_state *state);

diff --git a/include/typedefs.h b/include/typedefs.h
index c23aeb7f76754abb4659f3104ccc520d81d7a56f..a0853c7b74c1c5b9f12c204057a3b785c1844fcd 100644 (file)
--- a/include/typedefs.h
+++ b/include/typedefs.h
@@ -152,7 +152,9 @@ void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength);
 GMX_LIBGMX_EXPORT
 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda);
 GMX_LIBGMX_EXPORT
-void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta);
+void init_df_history(df_history_t *dfhist, int nlambda);
+GMX_LIBGMX_EXPORT
+void done_df_history(df_history_t *dfhist);
 GMX_LIBGMX_EXPORT
 void copy_df_history(df_history_t * df_dest, df_history_t *df_source);
 

diff --git a/include/types/inputrec.h b/include/types/inputrec.h
index 0652e66ecf16ff2a1aa97cf656efa78e3f6713dc..d3ab9218cf8d367293763a7d2cecc9562ca337b7 100644 (file)
--- a/include/types/inputrec.h
+++ b/include/types/inputrec.h
@@ -190,7 +190,7 @@ typedef struct {
     real     wl_ratio;            /* ratio between largest and smallest number for freezing the weights */
     real     init_wl_delta;       /* starting delta for wang-landau */
     gmx_bool bWLoneovert;         /* use one over t convergece for wang-landau when the delta get sufficiently small */
-    gmx_bool bInit_weights;       /* did we initialize the weights? */
+    gmx_bool bInit_weights;       /* did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic */
     real     mc_temp;             /* To override the main temperature, or define it if it's not defined */
     real    *init_lambda_weights; /* user-specified initial weights to start with  */
 } t_expanded;

diff --git a/include/types/state.h b/include/types/state.h
index fbad796183db83fd43d1d9e648549e757e196e25..2a0e975a79dd7ccf23b2cd5c5943bfe47663cb46 100644 (file)
--- a/include/types/state.h
+++ b/include/types/state.h
@@ -124,7 +124,7 @@ typedef struct
 {
     int      nlambda;        /* total number of lambda states - for history*/
 
-    gmx_bool bEquil;         /* reached equilibration */
+    gmx_bool bEquil;         /* Have we reached equilibration */
     int     *n_at_lam;       /* number of points observed at each lambda */
     real    *wl_histo;       /* histogram for WL flatness determination */
     real     wl_delta;       /* current wang-landau delta */
@@ -141,6 +141,7 @@ typedef struct
 
     real   **Tij;            /* transition matrix */
     real   **Tij_empirical;  /* Empirical transition matrix */
+
 } df_history_t;
 
 typedef struct

diff --git a/src/gmxlib/checkpoint.c b/src/gmxlib/checkpoint.c
index 53e49f6376b5c2676c5fd6882419e80c0cd88efe..ae1e316f458eb1d113a1ff4a1f46bf5ce4f30e1f 100644 (file)
--- a/src/gmxlib/checkpoint.c
+++ b/src/gmxlib/checkpoint.c
@@ -156,7 +156,7 @@ enum {
 /* free energy history variable names  */
 const char *edfh_names[edfhNR] =
 {
-    "bEquilibrated", "N_at_state", "Wang-Landau_Histogram", "Wang-Landau-delta", "Weights", "Free Energies", "minvar", "variance",
+    "bEquilibrated", "N_at_state", "Wang-Landau Histogram", "Wang-Landau Delta", "Weights", "Free Energies", "minvar", "variance",
     "accumulated_plus", "accumulated_minus", "accumulated_plus_2",  "accumulated_minus_2", "Tij", "Tij_empirical"
 };
 
@@ -1013,6 +1013,12 @@ static int do_cpt_state(XDR *xd, gmx_bool bRead,
         rng_p  = NULL;
         rngi_p = NULL;
     }
+
+    if (bRead) /* we need to allocate space for dfhist if we are reading */
+    {
+        init_df_history(&state->dfhist,state->dfhist.nlambda);
+    }
+
     /* We want the MC_RNG the same across all the notes for now -- lambda MC is global */
 
     sflags = state->flags;
@@ -2395,7 +2401,6 @@ void list_checkpoint(const char *fn, FILE *out)
 
     if (ret == 0)
     {
-        init_df_history(&state.dfhist, state.dfhist.nlambda, 0); /* reinitialize state with correct sizes */
         ret = do_cpt_df_hist(gmx_fio_getxdr(fp), TRUE,
                              flags_dfh, &state.dfhist, out);
     }

diff --git a/src/gmxlib/typedefs.c b/src/gmxlib/typedefs.c
index a928755687887f3b8ee5ff4e316b4ec88093ab06..2500c29b43dd45fb18bad761c4e8db2fd1587a29 100644 (file)
--- a/src/gmxlib/typedefs.c
+++ b/src/gmxlib/typedefs.c
@@ -624,7 +624,7 @@ void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainle
     zero_history(&state->hist);
     zero_ekinstate(&state->ekinstate);
     init_energyhistory(&state->enerhist);
-    init_df_history(&state->dfhist, nlambda, 0);
+    init_df_history(&state->dfhist,nlambda);
     state->ddp_count       = 0;
     state->ddp_count_cg_gl = 0;
     state->cg_gl           = NULL;
@@ -884,41 +884,43 @@ real max_cutoff(real cutoff1, real cutoff2)
     }
 }
 
-extern void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta)
+void init_df_history(df_history_t *dfhist, int nlambda)
 {
     int i;
 
-    dfhist->bEquil   = 0;
     dfhist->nlambda  = nlambda;
-    dfhist->wl_delta = wl_delta;
-    snew(dfhist->sum_weights, dfhist->nlambda);
-    snew(dfhist->sum_dg, dfhist->nlambda);
-    snew(dfhist->sum_minvar, dfhist->nlambda);
-    snew(dfhist->sum_variance, dfhist->nlambda);
-    snew(dfhist->n_at_lam, dfhist->nlambda);
-    snew(dfhist->wl_histo, dfhist->nlambda);
-
-    /* allocate transition matrices here */
-    snew(dfhist->Tij, dfhist->nlambda);
-    snew(dfhist->Tij_empirical, dfhist->nlambda);
+    dfhist->bEquil   = 0;
+    dfhist->wl_delta = 0;
 
-    for (i = 0; i < dfhist->nlambda; i++)
+    if (nlambda > 0)
     {
-        snew(dfhist->Tij[i], dfhist->nlambda);
-        snew(dfhist->Tij_empirical[i], dfhist->nlambda);
-    }
+        snew(dfhist->sum_weights, dfhist->nlambda);
+        snew(dfhist->sum_dg, dfhist->nlambda);
+        snew(dfhist->sum_minvar, dfhist->nlambda);
+        snew(dfhist->sum_variance, dfhist->nlambda);
+        snew(dfhist->n_at_lam, dfhist->nlambda);
+        snew(dfhist->wl_histo, dfhist->nlambda);
 
-    snew(dfhist->accum_p, dfhist->nlambda);
-    snew(dfhist->accum_m, dfhist->nlambda);
-    snew(dfhist->accum_p2, dfhist->nlambda);
-    snew(dfhist->accum_m2, dfhist->nlambda);
+        /* allocate transition matrices here */
+        snew(dfhist->Tij, dfhist->nlambda);
+        snew(dfhist->Tij_empirical, dfhist->nlambda);
 
-    for (i = 0; i < dfhist->nlambda; i++)
-    {
-        snew((dfhist->accum_p)[i], dfhist->nlambda);
-        snew((dfhist->accum_m)[i], dfhist->nlambda);
-        snew((dfhist->accum_p2)[i], dfhist->nlambda);
-        snew((dfhist->accum_m2)[i], dfhist->nlambda);
+        /* allocate accumulators for various transition matrix
+           free energy methods here */
+        snew(dfhist->accum_p, dfhist->nlambda);
+        snew(dfhist->accum_m, dfhist->nlambda);
+        snew(dfhist->accum_p2, dfhist->nlambda);
+        snew(dfhist->accum_m2, dfhist->nlambda);
+
+        for (i = 0; i < dfhist->nlambda; i++)
+        {
+            snew(dfhist->Tij[i], dfhist->nlambda);
+            snew(dfhist->Tij_empirical[i], dfhist->nlambda);
+            snew((dfhist->accum_p)[i], dfhist->nlambda);
+            snew((dfhist->accum_m)[i], dfhist->nlambda);
+            snew((dfhist->accum_p2)[i], dfhist->nlambda);
+            snew((dfhist->accum_m2)[i], dfhist->nlambda);
+        }
     }
 }
 
@@ -926,9 +928,12 @@ extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
 {
     int i, j;
 
-    init_df_history(df_dest, df_source->nlambda, df_source->wl_delta);
+    /* Currently, there should not be any difference in nlambda between the two,
+       but this is included for completeness for potential later functionality */
     df_dest->nlambda = df_source->nlambda;
     df_dest->bEquil  = df_source->bEquil;
+    df_dest->wl_delta = df_source->wl_delta;
+
     for (i = 0; i < df_dest->nlambda; i++)
     {
         df_dest->sum_weights[i]  = df_source->sum_weights[i];
@@ -937,18 +942,46 @@ extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
         df_dest->sum_variance[i] = df_source->sum_variance[i];
         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
         df_dest->wl_histo[i]     = df_source->wl_histo[i];
-        df_dest->accum_p[i]      = df_source->accum_p[i];
-        df_dest->accum_m[i]      = df_source->accum_m[i];
-        df_dest->accum_p2[i]     = df_source->accum_p2[i];
-        df_dest->accum_m2[i]     = df_source->accum_m2[i];
     }
 
     for (i = 0; i < df_dest->nlambda; i++)
     {
         for (j = 0; j < df_dest->nlambda; j++)
         {
+            df_dest->accum_p[i][j]      = df_source->accum_p[i][j];
+            df_dest->accum_m[i][j]      = df_source->accum_m[i][j];
+            df_dest->accum_p2[i][j]     = df_source->accum_p2[i][j];
+            df_dest->accum_m2[i][j]     = df_source->accum_m2[i][j];
             df_dest->Tij[i][j]            = df_source->Tij[i][j];
             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
         }
     }
 }
+
+void done_df_history(df_history_t *dfhist)
+{
+    int i;
+
+    if (dfhist->nlambda > 0)
+    {
+        sfree(dfhist->n_at_lam);
+        sfree(dfhist->wl_histo);
+        sfree(dfhist->sum_weights);
+        sfree(dfhist->sum_dg);
+        sfree(dfhist->sum_minvar);
+        sfree(dfhist->sum_variance);
+
+        for (i=0;i<dfhist->nlambda;i++)
+        {
+            sfree(dfhist->Tij[i]);
+            sfree(dfhist->Tij_empirical[i]);
+            sfree(dfhist->accum_p[i]);
+            sfree(dfhist->accum_m[i]);
+            sfree(dfhist->accum_p2[i]);
+            sfree(dfhist->accum_m2[i]);
+        }
+    }
+    dfhist->bEquil   = 0;
+    dfhist->nlambda  = 0;
+    dfhist->wl_delta = 0;
+}

diff --git a/src/kernel/md.c b/src/kernel/md.c
index de36746cbe44aeada7f201a028b610b38413a01f..81fca912a8746a7b6f297836762a0e5912f7ad86 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -181,7 +181,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     int               nchkpt  = 1;
     gmx_localtop_t   *top;
     t_mdebin         *mdebin = NULL;
-    df_history_t      df_history;
     t_state          *state    = NULL;
     rvec             *f_global = NULL;
     int               n_xtc    = -1;
@@ -330,14 +329,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         snew(f, top_global->natoms);
     }
 
-    /* lambda Monte carlo random number generator  */
-    if (ir->bExpanded)
-    {
-        mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
-    }
-    /* copy the state into df_history */
-    copy_df_history(&df_history, &state_global->dfhist);
-
     /* Kinetic energy data */
     snew(ekind, 1);
     init_ekindata(fplog, top_global, &(ir->opts), ekind);
@@ -462,6 +453,11 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         bStateFromCP = FALSE;
     }
 
+    if (ir->bExpanded)
+    {
+        init_expanded_ensemble(bStateFromCP,ir,&mcrng,&state->dfhist);
+    }
+
     if (MASTER(cr))
     {
         if (bStateFromCP)
@@ -858,7 +854,8 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
             bDoFEP       = (do_per_step(step, nstfep) && (ir->efep != efepNO));
-            bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded) && (ir->bExpanded) && (step > 0));
+            bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
+                            && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt));
         }
 
         if (bSimAnn)
@@ -1415,7 +1412,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                statistics, but if performing simulated tempering, we
                do update the velocities and the tau_t. */
 
-            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, &df_history, step, mcrng, state->v, mdatoms);
+            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, mcrng, state->v, mdatoms);
+            /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
+            copy_df_history(&state_global->dfhist,&state->dfhist);
         }
         /* ################## START TRAJECTORY OUTPUT ################# */
 
@@ -1494,14 +1493,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                         state_global->ekinstate.bUpToDate = TRUE;
                     }
                     update_energyhistory(&state_global->enerhist, mdebin);
-                    if (ir->efep != efepNO || ir->bSimTemp)
-                    {
-                        state_global->fep_state = state->fep_state; /* MRS: seems kludgy. The code should be
-                                                                       structured so this isn't necessary.
-                                                                       Note this reassignment is only necessary
-                                                                       for single threads.*/
-                        copy_df_history(&state_global->dfhist, &df_history);
-                    }
                 }
             }
             write_traj(fplog, cr, outf, mdof_flags, top_global,
@@ -2004,7 +1995,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             {
                 /* only needed if doing expanded ensemble */
                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
-                                          &df_history, state->fep_state, ir->nstlog, step);
+                                          &state_global->dfhist, state->fep_state, ir->nstlog, step);
             }
             if (!(bStartingFromCpt && (EI_VV(ir->eI))))
             {
@@ -2043,13 +2034,11 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         if (bDoExpanded)
         {
-            /* Have to do this part after outputting the logfile and the edr file */
+            /* Have to do this part _after_ outputting the logfile and the edr file */
+            /* Gets written into the state at the beginning of next loop*/
             state->fep_state = lamnew;
-            for (i = 0; i < efptNR; i++)
-            {
-                state_global->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
-            }
         }
+
         /* Remaining runtime */
         if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning)
         {

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index a9fd5bd6c1334d354bc4a99fd464fb1d542e1915..a8ad2bdf18a6fe01c8c6fa1be9ef7a2ce90f8a7c 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -1518,7 +1518,6 @@ static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weight
     parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
     if (nweights == 0)
     {
-        expand->bInit_weights = FALSE;
         snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
     }
     else if (nweights != fep->n_lambda)
@@ -1526,10 +1525,6 @@ static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weight
         gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
                   nweights, fep->n_lambda);
     }
-    else
-    {
-        expand->bInit_weights = TRUE;
-    }
     if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
     {
         expand->nstexpanded = fep->nstdhdl;

diff --git a/src/mdlib/domdec.c b/src/mdlib/domdec.c
index 046ef5133003f2cba1cb9953aba6e107a0dfe4cf..d78ed38a1d97d86ceadb0bea7a994b6b4255290c 100644 (file)
--- a/src/mdlib/domdec.c
+++ b/src/mdlib/domdec.c
@@ -1551,7 +1551,6 @@ void dd_collect_state(gmx_domdec_t *dd,
         copy_mat(state_local->fvir_prev, state->fvir_prev);
         copy_mat(state_local->pres_prev, state->pres_prev);
 
-
         for (i = 0; i < state_local->ngtc; i++)
         {
             for (j = 0; j < nh; j++)
@@ -1809,6 +1808,34 @@ static void dd_distribute_vec(gmx_domdec_t *dd, t_block *cgs, rvec *v, rvec *lv)
     }
 }
 
+static void dd_distribute_dfhist(gmx_domdec_t *dd, df_history_t *dfhist)
+{
+    int i;
+    dd_bcast(dd, sizeof(int), &dfhist->bEquil);
+    dd_bcast(dd, sizeof(int), &dfhist->nlambda);
+    dd_bcast(dd, sizeof(real), &dfhist->wl_delta);
+
+    if (dfhist->nlambda > 0)
+    {
+        int nlam = dfhist->nlambda;
+        dd_bcast(dd, sizeof(int)*nlam, dfhist->n_at_lam);
+        dd_bcast(dd, sizeof(real)*nlam, dfhist->wl_histo);
+        dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_weights);
+        dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_dg);
+        dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_minvar);
+        dd_bcast(dd, sizeof(real)*nlam, dfhist->sum_variance);
+
+        for (i = 0; i<nlam; i++) {
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p[i]);
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m[i]);
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_p2[i]);
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->accum_m2[i]);
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij[i]);
+            dd_bcast(dd, sizeof(real)*nlam, dfhist->Tij_empirical[i]);
+        }
+    }
+}
+
 static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
                                 t_state *state, t_state *state_local,
                                 rvec **f)
@@ -1831,6 +1858,7 @@ static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
         copy_mat(state->boxv, state_local->boxv);
         copy_mat(state->svir_prev, state_local->svir_prev);
         copy_mat(state->fvir_prev, state_local->fvir_prev);
+        copy_df_history(&state_local->dfhist,&state->dfhist);
         for (i = 0; i < state_local->ngtc; i++)
         {
             for (j = 0; j < nh; j++)
@@ -1864,6 +1892,9 @@ static void dd_distribute_state(gmx_domdec_t *dd, t_block *cgs,
     dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_xi);
     dd_bcast(dd, ((state_local->nnhpres*nh)*sizeof(double)), state_local->nhpres_vxi);
 
+    /* communicate df_history -- required for restarting from checkpoint */
+    dd_distribute_dfhist(dd,&state_local->dfhist);
+
     if (dd->nat_home > state_local->nalloc)
     {
         dd_realloc_state(state_local, f, dd->nat_home);

diff --git a/src/mdlib/expanded.c b/src/mdlib/expanded.c
index 672ebcdbdb660c3726d841992aa37fd272739fd6..a7bc879cb85d692051b47df19f1961811fb2893f 100644 (file)
--- a/src/mdlib/expanded.c
+++ b/src/mdlib/expanded.c
@@ -92,7 +92,30 @@
 #include "tmpi.h"
 #endif
 
-void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int minfep, int maxfep)
+
+static void init_df_history_weights(df_history_t *dfhist, t_expanded *expand, int nlim)
+{
+    int i;
+    dfhist->wl_delta = expand->init_wl_delta;
+    for (i = 0; i < nlim; i++)
+    {
+        dfhist->sum_weights[i] = expand->init_lambda_weights[i];
+        dfhist->sum_dg[i]      = expand->init_lambda_weights[i];
+    }
+}
+
+/* Eventually should contain all the functions needed to initialize expanded ensemble
+   before the md loop starts */
+extern void init_expanded_ensemble(gmx_bool bStateFromCP, t_inputrec *ir, gmx_rng_t *mcrng, df_history_t *dfhist)
+{
+    *mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
+    if (!bStateFromCP)
+    {
+        init_df_history_weights(dfhist,ir->expandedvals,ir->fepvals->n_lambda);
+    }
+}
+
+static void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int minfep, int maxfep)
 {
 
     int  i;
@@ -120,7 +143,7 @@ void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int minfep,
     }
 }
 
-void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *pks, int nlim, real *nvals, real delta)
+static void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *pks, int nlim, real *nvals, real delta)
 {
 
     int   i;
@@ -362,7 +385,7 @@ static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist
 {
     real     maxdiff = 0.000000001;
     gmx_bool bSufficientSamples;
-    int      i, k, n, nz, indexi, indexk, min_n, max_n, nlam, totali;
+    int      i, k, n, nz, indexi, indexk, min_n, max_n, totali;
     int      n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
     real     omega_m1_0, omega_p1_m1, omega_m1_p1, omega_p1_0, clam_osum;
     real     de, de_function, dr, denom, maxdr;
@@ -1022,7 +1045,7 @@ static int ChooseNewLambda(FILE *log, int nlim, t_expanded *expand, df_history_t
 
 /* print out the weights to the log, along with current state */
 extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *expand, t_simtemp *simtemp, df_history_t *dfhist,
-                                      int nlam, int frequency, gmx_large_int_t step)
+                                      int fep_state, int frequency, gmx_large_int_t step)
 {
     int         nlim, i, ifep, jfep;
     real        dw, dg, dv, dm, Tprint;
@@ -1116,7 +1139,7 @@ extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *
             {
                 fprintf(outfile, " %10.5f %10.5f", dfhist->sum_weights[ifep], dw);
             }
-            if (ifep == nlam)
+            if (ifep == fep_state)
             {
                 fprintf(outfile, " <<\n");
             }
@@ -1203,13 +1226,15 @@ extern void set_mc_state(gmx_rng_t rng, t_state *state)
 }
 
 extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *enerd,
-                                    t_state *state, t_extmass *MassQ, df_history_t *dfhist,
+                                    t_state *state, t_extmass *MassQ, int fep_state, df_history_t *dfhist,
                                     gmx_large_int_t step, gmx_rng_t mcrng,
                                     rvec *v, t_mdatoms *mdatoms)
+/* Note that the state variable is only needed for simulated tempering, not
+   Hamiltonian expanded ensemble.  May be able to remove it after integrator refactoring. */
 {
     real       *pfep_lamee, *scaled_lamee, *weighted_lamee;
     double     *p_k;
-    int         i, nlam, nlim, lamnew, totalsamples;
+    int         i, nlim, lamnew, totalsamples;
     real        oneovert, maxscaled = 0, maxweighted = 0;
     t_expanded *expand;
     t_simtemp  *simtemp;
@@ -1219,26 +1244,14 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
     expand  = ir->expandedvals;
     simtemp = ir->simtempvals;
     nlim    = ir->fepvals->n_lambda;
-    nlam    = state->fep_state;
 
     snew(scaled_lamee, nlim);
     snew(weighted_lamee, nlim);
     snew(pfep_lamee, nlim);
     snew(p_k, nlim);
 
-    if (expand->bInit_weights)                    /* if initialized weights, we need to fill them in */
-    {
-        dfhist->wl_delta = expand->init_wl_delta; /* MRS -- this would fit better somewhere else? */
-        for (i = 0; i < nlim; i++)
-        {
-            dfhist->sum_weights[i] = expand->init_lambda_weights[i];
-            dfhist->sum_dg[i]      = expand->init_lambda_weights[i];
-        }
-        expand->bInit_weights = FALSE;
-    }
-
     /* update the count at the current lambda*/
-    dfhist->n_at_lam[nlam]++;
+    dfhist->n_at_lam[fep_state]++;
 
     /* need to calculate the PV term somewhere, but not needed here? Not until there's a lambda state that's
        pressure controlled.*/
@@ -1263,7 +1276,7 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
             {
                 /* Note -- this assumes no mass changes, since kinetic energy is not added  . . . */
                 scaled_lamee[i] = (enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0])/(simtemp->temperatures[i]*BOLTZ)
-                    + enerd->term[F_EPOT]*(1.0/(simtemp->temperatures[i])- 1.0/(simtemp->temperatures[nlam]))/BOLTZ;
+                    + enerd->term[F_EPOT]*(1.0/(simtemp->temperatures[i])- 1.0/(simtemp->temperatures[fep_state]))/BOLTZ;
             }
             else
             {
@@ -1282,7 +1295,7 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
         {
             for (i = 0; i < nlim; i++)
             {
-                scaled_lamee[i] = enerd->term[F_EPOT]*(1.0/simtemp->temperatures[i] - 1.0/simtemp->temperatures[nlam])/BOLTZ;
+                scaled_lamee[i] = enerd->term[F_EPOT]*(1.0/simtemp->temperatures[i] - 1.0/simtemp->temperatures[fep_state])/BOLTZ;
             }
         }
     }
@@ -1318,7 +1331,7 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
 
     /* update weights - we decide whether or not to actually do this inside */
 
-    bDoneEquilibrating = UpdateWeights(nlim, expand, dfhist, nlam, scaled_lamee, weighted_lamee, step);
+    bDoneEquilibrating = UpdateWeights(nlim, expand, dfhist, fep_state, scaled_lamee, weighted_lamee, step);
     if (bDoneEquilibrating)
     {
         if (log)
@@ -1327,9 +1340,9 @@ extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *e
         }
     }
 
-    lamnew = ChooseNewLambda(log, nlim, expand, dfhist, nlam, weighted_lamee, p_k, mcrng);
+    lamnew = ChooseNewLambda(log, nlim, expand, dfhist, fep_state, weighted_lamee, p_k, mcrng);
     /* if using simulated tempering, we need to adjust the temperatures */
-    if (ir->bSimTemp && (lamnew != nlam)) /* only need to change the temperatures if we change the state */
+    if (ir->bSimTemp && (lamnew != fep_state)) /* only need to change the temperatures if we change the state */
     {
         int   i, j, n, d;
         real *buf_ngtc;

diff --git a/src/mdlib/init.c b/src/mdlib/init.c
index e957cf46c2b6ae1ec6f7fc9cdab54aa5a2756abc..50c0a2d83e8aa00cf75c58a18be977e1e4a00247 100644 (file)
--- a/src/mdlib/init.c
+++ b/src/mdlib/init.c
@@ -190,7 +190,7 @@ void set_state_entries(t_state *state, const t_inputrec *ir, int nnodes)
     init_ekinstate(&state->ekinstate, ir);
 
     init_energyhistory(&state->enerhist);
-    init_df_history(&state->dfhist, ir->fepvals->n_lambda, ir->expandedvals->init_wl_delta);
+    init_df_history(&state->dfhist, ir->fepvals->n_lambda);
 }
 
 

diff --git a/src/mdlib/md_support.c b/src/mdlib/md_support.c
index 9ff6d07a65a1399a73ec67d1c4dd2925cc0b4db4..b1b73bab5992810ef85239c30e1a10975fd4f23e 100644 (file)
--- a/src/mdlib/md_support.c
+++ b/src/mdlib/md_support.c
@@ -582,6 +582,16 @@ void set_current_lambdas(gmx_large_int_t step, t_lambda *fepvals, gmx_bool bReru
                 }
             }
         }
+        else
+        {
+            if (state->fep_state > 0) {
+                state_global->fep_state = state->fep_state; /* state->fep is the one updated by bExpanded */
+                for (i = 0; i < efptNR; i++)
+                {
+                    state_global->lambda[i] = fepvals->all_lambda[i][state_global->fep_state];
+                }
+            }
+        }
     }
     for (i = 0; i < efptNR; i++)
     {