src/gromacs/mdlib/expanded.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include <math.h>
  38 #include <stdio.h>
  39
  40 #include "gromacs/legacyheaders/typedefs.h"
  41 #include "gromacs/utility/smalloc.h"
  42 #include "gromacs/legacyheaders/names.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "gromacs/legacyheaders/txtdump.h"
  45 #include "gromacs/legacyheaders/chargegroup.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/legacyheaders/nrnb.h"
  48 #include "gromacs/legacyheaders/mdrun.h"
  49 #include "gromacs/legacyheaders/update.h"
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/legacyheaders/mdatoms.h"
  52 #include "gromacs/legacyheaders/force.h"
  53 #include "gromacs/legacyheaders/pme.h"
  54 #include "gromacs/legacyheaders/disre.h"
  55 #include "gromacs/legacyheaders/orires.h"
  56 #include "gromacs/legacyheaders/network.h"
  57 #include "gromacs/legacyheaders/calcmu.h"
  58 #include "gromacs/legacyheaders/constr.h"
  59 #include "gromacs/random/random.h"
  60 #include "gromacs/legacyheaders/domdec.h"
  61 #include "gromacs/legacyheaders/macros.h"
  62
  63 #include "gromacs/fileio/confio.h"
  64 #include "gromacs/fileio/gmxfio.h"
  65 #include "gromacs/fileio/trnio.h"
  66 #include "gromacs/fileio/xtcio.h"
  67 #include "gromacs/timing/wallcycle.h"
  68 #include "gromacs/utility/fatalerror.h"
  69 #include "gromacs/utility/gmxmpi.h"
  70
  71 static void init_df_history_weights(df_history_t *dfhist, t_expanded *expand, int nlim)
  72 {
  73     int i;
  74     dfhist->wl_delta = expand->init_wl_delta;
  75     for (i = 0; i < nlim; i++)
  76     {
  77         dfhist->sum_weights[i] = expand->init_lambda_weights[i];
  78         dfhist->sum_dg[i]      = expand->init_lambda_weights[i];
  79     }
  80 }
  81
  82 /* Eventually should contain all the functions needed to initialize expanded ensemble
  83    before the md loop starts */
  84 extern void init_expanded_ensemble(gmx_bool bStateFromCP, t_inputrec *ir, df_history_t *dfhist)
  85 {
  86     if (!bStateFromCP)
  87     {
  88         init_df_history_weights(dfhist, ir->expandedvals, ir->fepvals->n_lambda);
  89     }
  90 }
  91
  92 static void GenerateGibbsProbabilities(real *ene, double *p_k, double *pks, int minfep, int maxfep)
  93 {
  94
  95     int  i;
  96     real maxene;
  97
  98     *pks   = 0.0;
  99     maxene = ene[minfep];
 100     /* find the maximum value */
 101     for (i = minfep; i <= maxfep; i++)
 102     {
 103         if (ene[i] > maxene)
 104         {
 105             maxene = ene[i];
 106         }
 107     }
 108     /* find the denominator */
 109     for (i = minfep; i <= maxfep; i++)
 110     {
 111         *pks += exp(ene[i]-maxene);
 112     }
 113     /*numerators*/
 114     for (i = minfep; i <= maxfep; i++)
 115     {
 116         p_k[i] = exp(ene[i]-maxene) / *pks;
 117     }
 118 }
 119
 120 static void GenerateWeightedGibbsProbabilities(real *ene, double *p_k, double *pks, int nlim, real *nvals, real delta)
 121 {
 122
 123     int   i;
 124     real  maxene;
 125     real *nene;
 126     *pks = 0.0;
 127
 128     snew(nene, nlim);
 129     for (i = 0; i < nlim; i++)
 130     {
 131         if (nvals[i] == 0)
 132         {
 133             /* add the delta, since we need to make sure it's greater than zero, and
 134                we need a non-arbitrary number? */
 135             nene[i] = ene[i] + log(nvals[i]+delta);
 136         }
 137         else
 138         {
 139             nene[i] = ene[i] + log(nvals[i]);
 140         }
 141     }
 142
 143     /* find the maximum value */
 144     maxene = nene[0];
 145     for (i = 0; i < nlim; i++)
 146     {
 147         if (nene[i] > maxene)
 148         {
 149             maxene = nene[i];
 150         }
 151     }
 152
 153     /* subtract off the maximum, avoiding overflow */
 154     for (i = 0; i < nlim; i++)
 155     {
 156         nene[i] -= maxene;
 157     }
 158
 159     /* find the denominator */
 160     for (i = 0; i < nlim; i++)
 161     {
 162         *pks += exp(nene[i]);
 163     }
 164
 165     /*numerators*/
 166     for (i = 0; i < nlim; i++)
 167     {
 168         p_k[i] = exp(nene[i]) / *pks;
 169     }
 170     sfree(nene);
 171 }
 172
 173 real do_logsum(int N, real *a_n)
 174 {
 175
 176     /*     RETURN VALUE */
 177     /* log(\sum_{i=0}^(N-1) exp[a_n]) */
 178     real maxarg;
 179     real sum;
 180     int  i;
 181     real logsum;
 182     /*     compute maximum argument to exp(.) */
 183
 184     maxarg = a_n[0];
 185     for (i = 1; i < N; i++)
 186     {
 187         maxarg = max(maxarg, a_n[i]);
 188     }
 189
 190     /* compute sum of exp(a_n - maxarg) */
 191     sum = 0.0;
 192     for (i = 0; i < N; i++)
 193     {
 194         sum = sum + exp(a_n[i] - maxarg);
 195     }
 196
 197     /*     compute log sum */
 198     logsum = log(sum) + maxarg;
 199     return logsum;
 200 }
 201
 202 int FindMinimum(real *min_metric, int N)
 203 {
 204
 205     real min_val;
 206     int  min_nval, nval;
 207
 208     min_nval = 0;
 209     min_val  = min_metric[0];
 210
 211     for (nval = 0; nval < N; nval++)
 212     {
 213         if (min_metric[nval] < min_val)
 214         {
 215             min_val  = min_metric[nval];
 216             min_nval = nval;
 217         }
 218     }
 219     return min_nval;
 220 }
 221
 222 static gmx_bool CheckHistogramRatios(int nhisto, real *histo, real ratio)
 223 {
 224
 225     int      i;
 226     real     nmean;
 227     gmx_bool bIfFlat;
 228
 229     nmean = 0;
 230     for (i = 0; i < nhisto; i++)
 231     {
 232         nmean += histo[i];
 233     }
 234
 235     if (nmean == 0)
 236     {
 237         /* no samples! is bad!*/
 238         bIfFlat = FALSE;
 239         return bIfFlat;
 240     }
 241     nmean /= (real)nhisto;
 242
 243     bIfFlat = TRUE;
 244     for (i = 0; i < nhisto; i++)
 245     {
 246         /* make sure that all points are in the ratio < x <  1/ratio range  */
 247         if (!((histo[i]/nmean < 1.0/ratio) && (histo[i]/nmean > ratio)))
 248         {
 249             bIfFlat = FALSE;
 250             break;
 251         }
 252     }
 253     return bIfFlat;
 254 }
 255
 256 static gmx_bool CheckIfDoneEquilibrating(int nlim, t_expanded *expand, df_history_t *dfhist, gmx_int64_t step)
 257 {
 258
 259     int      i, totalsamples;
 260     gmx_bool bDoneEquilibrating = TRUE;
 261     gmx_bool bIfFlat;
 262
 263     /* assume we have equilibrated the weights, then check to see if any of the conditions are not met */
 264
 265     /* calculate the total number of samples */
 266     switch (expand->elmceq)
 267     {
 268         case elmceqNO:
 269             /* We have not equilibrated, and won't, ever. */
 270             return FALSE;
 271         case elmceqYES:
 272             /* we have equilibrated -- we're done */
 273             return TRUE;
 274         case elmceqSTEPS:
 275             /* first, check if we are equilibrating by steps, if we're still under */
 276             if (step < expand->equil_steps)
 277             {
 278                 bDoneEquilibrating = FALSE;
 279             }
 280             break;
 281         case elmceqSAMPLES:
 282             totalsamples = 0;
 283             for (i = 0; i < nlim; i++)
 284             {
 285                 totalsamples += dfhist->n_at_lam[i];
 286             }
 287             if (totalsamples < expand->equil_samples)
 288             {
 289                 bDoneEquilibrating = FALSE;
 290             }
 291             break;
 292         case elmceqNUMATLAM:
 293             for (i = 0; i < nlim; i++)
 294             {
 295                 if (dfhist->n_at_lam[i] < expand->equil_n_at_lam) /* we are still doing the initial sweep, so we're definitely not
 296                                                                      done equilibrating*/
 297                 {
 298                     bDoneEquilibrating  = FALSE;
 299                     break;
 300                 }
 301             }
 302             break;
 303         case elmceqWLDELTA:
 304             if (EWL(expand->elamstats)) /* This check is in readir as well, but
 305                                            just to be sure */
 306             {
 307                 if (dfhist->wl_delta > expand->equil_wl_delta)
 308                 {
 309                     bDoneEquilibrating = FALSE;
 310                 }
 311             }
 312             break;
 313         case elmceqRATIO:
 314             /* we can use the flatness as a judge of good weights, as long as
 315                we're not doing minvar, or Wang-Landau.
 316                But turn off for now until we figure out exactly how we do this.
 317              */
 318
 319             if (!(EWL(expand->elamstats) || expand->elamstats == elamstatsMINVAR))
 320             {
 321                 /* we want to use flatness -avoiding- the forced-through samples.  Plus, we need to convert to
 322                    floats for this histogram function. */
 323
 324                 real *modhisto;
 325                 snew(modhisto, nlim);
 326                 for (i = 0; i < nlim; i++)
 327                 {
 328                     modhisto[i] = 1.0*(dfhist->n_at_lam[i]-expand->lmc_forced_nstart);
 329                 }
 330                 bIfFlat = CheckHistogramRatios(nlim, modhisto, expand->equil_ratio);
 331                 sfree(modhisto);
 332                 if (!bIfFlat)
 333                 {
 334                     bDoneEquilibrating = FALSE;
 335                 }
 336             }
 337         default:
 338             bDoneEquilibrating = TRUE;
 339     }
 340     /* one last case to go though, if we are doing slow growth to get initial values, we haven't finished equilibrating */
 341
 342     if (expand->lmc_forced_nstart > 0)
 343     {
 344         for (i = 0; i < nlim; i++)
 345         {
 346             if (dfhist->n_at_lam[i] < expand->lmc_forced_nstart) /* we are still doing the initial sweep, so we're definitely not
 347                                                                     done equilibrating*/
 348             {
 349                 bDoneEquilibrating = FALSE;
 350                 break;
 351             }
 352         }
 353     }
 354     return bDoneEquilibrating;
 355 }
 356
 357 static gmx_bool UpdateWeights(int nlim, t_expanded *expand, df_history_t *dfhist,
 358                               int fep_state, real *scaled_lamee, real *weighted_lamee, gmx_int64_t step)
 359 {
 360     real     maxdiff = 0.000000001;
 361     gmx_bool bSufficientSamples;
 362     int      i, k, n, nz, indexi, indexk, min_n, max_n, totali;
 363     int      n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
 364     real     omega_m1_0, omega_p1_m1, omega_m1_p1, omega_p1_0, clam_osum;
 365     real     de, de_function, dr, denom, maxdr;
 366     real     min_val, cnval, zero_sum_weights;
 367     real    *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array, *dwp_array, *dwm_array;
 368     real     clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar;
 369     real    *lam_weights, *lam_minvar_corr, *lam_variance, *lam_dg;
 370     double  *p_k;
 371     double   pks = 0;
 372     real    *numweighted_lamee, *logfrac;
 373     int     *nonzero;
 374     real     chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1;
 375
 376     /* if we have equilibrated the weights, exit now */
 377     if (dfhist->bEquil)
 378     {
 379         return FALSE;
 380     }
 381
 382     if (CheckIfDoneEquilibrating(nlim, expand, dfhist, step))
 383     {
 384         dfhist->bEquil = TRUE;
 385         /* zero out the visited states so we know how many equilibrated states we have
 386            from here on out.*/
 387         for (i = 0; i < nlim; i++)
 388         {
 389             dfhist->n_at_lam[i] = 0;
 390         }
 391         return TRUE;
 392     }
 393
 394     /* If we reached this far, we have not equilibrated yet, keep on
 395        going resetting the weights */
 396
 397     if (EWL(expand->elamstats))
 398     {
 399         if (expand->elamstats == elamstatsWL)  /* Standard Wang-Landau */
 400         {
 401             dfhist->sum_weights[fep_state] -= dfhist->wl_delta;
 402             dfhist->wl_histo[fep_state]    += 1.0;
 403         }
 404         else if (expand->elamstats == elamstatsWWL) /* Weighted Wang-Landau */
 405         {
 406             snew(p_k, nlim);
 407
 408             /* first increment count */
 409             GenerateGibbsProbabilities(weighted_lamee, p_k, &pks, 0, nlim-1);
 410             for (i = 0; i < nlim; i++)
 411             {
 412                 dfhist->wl_histo[i] += (real)p_k[i];
 413             }
 414
 415             /* then increment weights (uses count) */
 416             pks = 0.0;
 417             GenerateWeightedGibbsProbabilities(weighted_lamee, p_k, &pks, nlim, dfhist->wl_histo, dfhist->wl_delta);
 418
 419             for (i = 0; i < nlim; i++)
 420             {
 421                 dfhist->sum_weights[i] -= dfhist->wl_delta*(real)p_k[i];
 422             }
 423             /* Alternate definition, using logarithms. Shouldn't make very much difference! */
 424             /*
 425                real di;
 426                for (i=0;i<nlim;i++)
 427                {
 428                 di = (real)1.0 + dfhist->wl_delta*(real)p_k[i];
 429                 dfhist->sum_weights[i] -= log(di);
 430                }
 431              */
 432             sfree(p_k);
 433         }
 434
 435         zero_sum_weights =  dfhist->sum_weights[0];
 436         for (i = 0; i < nlim; i++)
 437         {
 438             dfhist->sum_weights[i] -= zero_sum_weights;
 439         }
 440     }
 441
 442     if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMETROPOLIS || expand->elamstats == elamstatsMINVAR)
 443     {
 444
 445         de_function = 0;  /* to get rid of warnings, but this value will not be used because of the logic */
 446         maxc        = 2*expand->c_range+1;
 447
 448         snew(lam_dg, nlim);
 449         snew(lam_variance, nlim);
 450
 451         snew(omegap_array, maxc);
 452         snew(weightsp_array, maxc);
 453         snew(varp_array, maxc);
 454         snew(dwp_array, maxc);
 455
 456         snew(omegam_array, maxc);
 457         snew(weightsm_array, maxc);
 458         snew(varm_array, maxc);
 459         snew(dwm_array, maxc);
 460
 461         /* unpack the current lambdas -- we will only update 2 of these */
 462
 463         for (i = 0; i < nlim-1; i++)
 464         {   /* only through the second to last */
 465             lam_dg[i]       = dfhist->sum_dg[i+1] - dfhist->sum_dg[i];
 466             lam_variance[i] = pow(dfhist->sum_variance[i+1], 2) - pow(dfhist->sum_variance[i], 2);
 467         }
 468
 469         /* accumulate running averages */
 470         for (nval = 0; nval < maxc; nval++)
 471         {
 472             /* constants for later use */
 473             cnval = (real)(nval-expand->c_range);
 474             /* actually, should be able to rewrite it w/o exponential, for better numerical stability */
 475             if (fep_state > 0)
 476             {
 477                 de = exp(cnval - (scaled_lamee[fep_state]-scaled_lamee[fep_state-1]));
 478                 if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
 479                 {
 480                     de_function = 1.0/(1.0+de);
 481                 }
 482                 else if (expand->elamstats == elamstatsMETROPOLIS)
 483                 {
 484                     if (de < 1.0)
 485                     {
 486                         de_function = 1.0;
 487                     }
 488                     else
 489                     {
 490                         de_function = 1.0/de;
 491                     }
 492                 }
 493                 dfhist->accum_m[fep_state][nval]  += de_function;
 494                 dfhist->accum_m2[fep_state][nval] += de_function*de_function;
 495             }
 496
 497             if (fep_state < nlim-1)
 498             {
 499                 de = exp(-cnval + (scaled_lamee[fep_state+1]-scaled_lamee[fep_state]));
 500                 if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
 501                 {
 502                     de_function = 1.0/(1.0+de);
 503                 }
 504                 else if (expand->elamstats == elamstatsMETROPOLIS)
 505                 {
 506                     if (de < 1.0)
 507                     {
 508                         de_function = 1.0;
 509                     }
 510                     else
 511                     {
 512                         de_function = 1.0/de;
 513                     }
 514                 }
 515                 dfhist->accum_p[fep_state][nval]  += de_function;
 516                 dfhist->accum_p2[fep_state][nval] += de_function*de_function;
 517             }
 518
 519             /* Metropolis transition and Barker transition (unoptimized Bennett) acceptance weight determination */
 520
 521             n0  = dfhist->n_at_lam[fep_state];
 522             if (fep_state > 0)
 523             {
 524                 nm1 = dfhist->n_at_lam[fep_state-1];
 525             }
 526             else
 527             {
 528                 nm1 = 0;
 529             }
 530             if (fep_state < nlim-1)
 531             {
 532                 np1 = dfhist->n_at_lam[fep_state+1];
 533             }
 534             else
 535             {
 536                 np1 = 0;
 537             }
 538
 539             /* logic SHOULD keep these all set correctly whatever the logic, but apparently it can't figure it out. */
 540             chi_m1_0 = chi_p1_0 = chi_m2_0 = chi_p2_0 = chi_p1_m1 = chi_p2_m1 = chi_m1_p1 = chi_m2_p1 = 0;
 541
 542             if (n0 > 0)
 543             {
 544                 chi_m1_0 = dfhist->accum_m[fep_state][nval]/n0;
 545                 chi_p1_0 = dfhist->accum_p[fep_state][nval]/n0;
 546                 chi_m2_0 = dfhist->accum_m2[fep_state][nval]/n0;
 547                 chi_p2_0 = dfhist->accum_p2[fep_state][nval]/n0;
 548             }
 549
 550             if ((fep_state > 0 ) && (nm1 > 0))
 551             {
 552                 chi_p1_m1 = dfhist->accum_p[fep_state-1][nval]/nm1;
 553                 chi_p2_m1 = dfhist->accum_p2[fep_state-1][nval]/nm1;
 554             }
 555
 556             if ((fep_state < nlim-1) && (np1 > 0))
 557             {
 558                 chi_m1_p1 = dfhist->accum_m[fep_state+1][nval]/np1;
 559                 chi_m2_p1 = dfhist->accum_m2[fep_state+1][nval]/np1;
 560             }
 561
 562             omega_m1_0    = 0;
 563             omega_p1_0    = 0;
 564             clam_weightsm = 0;
 565             clam_weightsp = 0;
 566             clam_varm     = 0;
 567             clam_varp     = 0;
 568
 569             if (fep_state > 0)
 570             {
 571                 if (n0 > 0)
 572                 {
 573                     omega_m1_0 = chi_m2_0/(chi_m1_0*chi_m1_0) - 1.0;
 574                 }
 575                 if (nm1 > 0)
 576                 {
 577                     omega_p1_m1 = chi_p2_m1/(chi_p1_m1*chi_p1_m1) - 1.0;
 578                 }
 579                 if ((n0 > 0) && (nm1 > 0))
 580                 {
 581                     clam_weightsm = (log(chi_m1_0) - log(chi_p1_m1)) + cnval;
 582                     clam_varm     = (1.0/n0)*(omega_m1_0) + (1.0/nm1)*(omega_p1_m1);
 583                 }
 584             }
 585
 586             if (fep_state < nlim-1)
 587             {
 588                 if (n0 > 0)
 589                 {
 590                     omega_p1_0 = chi_p2_0/(chi_p1_0*chi_p1_0) - 1.0;
 591                 }
 592                 if (np1 > 0)
 593                 {
 594                     omega_m1_p1 = chi_m2_p1/(chi_m1_p1*chi_m1_p1) - 1.0;
 595                 }
 596                 if ((n0 > 0) && (np1 > 0))
 597                 {
 598                     clam_weightsp = (log(chi_m1_p1) - log(chi_p1_0)) + cnval;
 599                     clam_varp     = (1.0/np1)*(omega_m1_p1) + (1.0/n0)*(omega_p1_0);
 600                 }
 601             }
 602
 603             if (n0 > 0)
 604             {
 605                 omegam_array[nval]             = omega_m1_0;
 606             }
 607             else
 608             {
 609                 omegam_array[nval]             = 0;
 610             }
 611             weightsm_array[nval]           = clam_weightsm;
 612             varm_array[nval]               = clam_varm;
 613             if (nm1 > 0)
 614             {
 615                 dwm_array[nval]  = fabs( (cnval + log((1.0*n0)/nm1)) - lam_dg[fep_state-1] );
 616             }
 617             else
 618             {
 619                 dwm_array[nval]  = fabs( cnval - lam_dg[fep_state-1] );
 620             }
 621
 622             if (n0 > 0)
 623             {
 624                 omegap_array[nval]             = omega_p1_0;
 625             }
 626             else
 627             {
 628                 omegap_array[nval]             = 0;
 629             }
 630             weightsp_array[nval]           = clam_weightsp;
 631             varp_array[nval]               = clam_varp;
 632             if ((np1 > 0) && (n0 > 0))
 633             {
 634                 dwp_array[nval]  = fabs( (cnval + log((1.0*np1)/n0)) - lam_dg[fep_state] );
 635             }
 636             else
 637             {
 638                 dwp_array[nval]  = fabs( cnval - lam_dg[fep_state] );
 639             }
 640
 641         }
 642
 643         /* find the C's closest to the old weights value */
 644
 645         min_nvalm     = FindMinimum(dwm_array, maxc);
 646         omega_m1_0    = omegam_array[min_nvalm];
 647         clam_weightsm = weightsm_array[min_nvalm];
 648         clam_varm     = varm_array[min_nvalm];
 649
 650         min_nvalp     = FindMinimum(dwp_array, maxc);
 651         omega_p1_0    = omegap_array[min_nvalp];
 652         clam_weightsp = weightsp_array[min_nvalp];
 653         clam_varp     = varp_array[min_nvalp];
 654
 655         clam_osum   = omega_m1_0 + omega_p1_0;
 656         clam_minvar = 0;
 657         if (clam_osum > 0)
 658         {
 659             clam_minvar = 0.5*log(clam_osum);
 660         }
 661
 662         if (fep_state > 0)
 663         {
 664             lam_dg[fep_state-1]       = clam_weightsm;
 665             lam_variance[fep_state-1] = clam_varm;
 666         }
 667
 668         if (fep_state < nlim-1)
 669         {
 670             lam_dg[fep_state]       = clam_weightsp;
 671             lam_variance[fep_state] = clam_varp;
 672         }
 673
 674         if (expand->elamstats == elamstatsMINVAR)
 675         {
 676             bSufficientSamples = TRUE;
 677             /* make sure they are all past a threshold */
 678             for (i = 0; i < nlim; i++)
 679             {
 680                 if (dfhist->n_at_lam[i] < expand->minvarmin)
 681                 {
 682                     bSufficientSamples = FALSE;
 683                 }
 684             }
 685             if (bSufficientSamples)
 686             {
 687                 dfhist->sum_minvar[fep_state] = clam_minvar;
 688                 if (fep_state == 0)
 689                 {
 690                     for (i = 0; i < nlim; i++)
 691                     {
 692                         dfhist->sum_minvar[i] += (expand->minvar_const-clam_minvar);
 693                     }
 694                     expand->minvar_const          = clam_minvar;
 695                     dfhist->sum_minvar[fep_state] = 0.0;
 696                 }
 697                 else
 698                 {
 699                     dfhist->sum_minvar[fep_state] -= expand->minvar_const;
 700                 }
 701             }
 702         }
 703
 704         /* we need to rezero minvar now, since it could change at fep_state = 0 */
 705         dfhist->sum_dg[0]       = 0.0;
 706         dfhist->sum_variance[0] = 0.0;
 707         dfhist->sum_weights[0]  = dfhist->sum_dg[0] + dfhist->sum_minvar[0]; /* should be zero */
 708
 709         for (i = 1; i < nlim; i++)
 710         {
 711             dfhist->sum_dg[i]       = lam_dg[i-1] + dfhist->sum_dg[i-1];
 712             dfhist->sum_variance[i] = sqrt(lam_variance[i-1] + pow(dfhist->sum_variance[i-1], 2));
 713             dfhist->sum_weights[i]  = dfhist->sum_dg[i] + dfhist->sum_minvar[i];
 714         }
 715
 716         sfree(lam_dg);
 717         sfree(lam_variance);
 718
 719         sfree(omegam_array);
 720         sfree(weightsm_array);
 721         sfree(varm_array);
 722         sfree(dwm_array);
 723
 724         sfree(omegap_array);
 725         sfree(weightsp_array);
 726         sfree(varp_array);
 727         sfree(dwp_array);
 728     }
 729     return FALSE;
 730 }
 731
 732 static int ChooseNewLambda(int nlim, t_expanded *expand, df_history_t *dfhist, int fep_state, real *weighted_lamee, double *p_k,
 733                            gmx_int64_t seed, gmx_int64_t step)
 734 {
 735     /* Choose new lambda value, and update transition matrix */
 736
 737     int      i, ifep, jfep, minfep, maxfep, lamnew, lamtrial, starting_fep_state;
 738     real     r1, r2, de_old, de_new, de, trialprob, tprob = 0;
 739     real   **Tij;
 740     double  *propose, *accept, *remainder;
 741     double   pks;
 742     real     sum, pnorm;
 743     gmx_bool bRestricted;
 744
 745     starting_fep_state = fep_state;
 746     lamnew             = fep_state; /* so that there is a default setting -- stays the same */
 747
 748     if (!EWL(expand->elamstats))    /* ignore equilibrating the weights if using WL */
 749     {
 750         if ((expand->lmc_forced_nstart > 0) && (dfhist->n_at_lam[nlim-1] <= expand->lmc_forced_nstart))
 751         {
 752             /* Use a marching method to run through the lambdas and get preliminary free energy data,
 753                before starting 'free' sampling.  We start free sampling when we have enough at each lambda */
 754
 755             /* if we have enough at this lambda, move on to the next one */
 756
 757             if (dfhist->n_at_lam[fep_state] == expand->lmc_forced_nstart)
 758             {
 759                 lamnew = fep_state+1;
 760                 if (lamnew == nlim)  /* whoops, stepped too far! */
 761                 {
 762                     lamnew -= 1;
 763                 }
 764             }
 765             else
 766             {
 767                 lamnew = fep_state;
 768             }
 769             return lamnew;
 770         }
 771     }
 772
 773     snew(propose, nlim);
 774     snew(accept, nlim);
 775     snew(remainder, nlim);
 776
 777     for (i = 0; i < expand->lmc_repeats; i++)
 778     {
 779         double rnd[2];
 780
 781         gmx_rng_cycle_2uniform(step, i, seed, RND_SEED_EXPANDED, rnd);
 782
 783         for (ifep = 0; ifep < nlim; ifep++)
 784         {
 785             propose[ifep] = 0;
 786             accept[ifep]  = 0;
 787         }
 788
 789         if ((expand->elmcmove == elmcmoveGIBBS) || (expand->elmcmove == elmcmoveMETGIBBS))
 790         {
 791             bRestricted = TRUE;
 792             /* use the Gibbs sampler, with restricted range */
 793             if (expand->gibbsdeltalam < 0)
 794             {
 795                 minfep      = 0;
 796                 maxfep      = nlim-1;
 797                 bRestricted = FALSE;
 798             }
 799             else
 800             {
 801                 minfep = fep_state - expand->gibbsdeltalam;
 802                 maxfep = fep_state + expand->gibbsdeltalam;
 803                 if (minfep < 0)
 804                 {
 805                     minfep = 0;
 806                 }
 807                 if (maxfep > nlim-1)
 808                 {
 809                     maxfep = nlim-1;
 810                 }
 811             }
 812
 813             GenerateGibbsProbabilities(weighted_lamee, p_k, &pks, minfep, maxfep);
 814
 815             if (expand->elmcmove == elmcmoveGIBBS)
 816             {
 817                 for (ifep = minfep; ifep <= maxfep; ifep++)
 818                 {
 819                     propose[ifep] = p_k[ifep];
 820                     accept[ifep]  = 1.0;
 821                 }
 822                 /* Gibbs sampling */
 823                 r1 = rnd[0];
 824                 for (lamnew = minfep; lamnew <= maxfep; lamnew++)
 825                 {
 826                     if (r1 <= p_k[lamnew])
 827                     {
 828                         break;
 829                     }
 830                     r1 -= p_k[lamnew];
 831                 }
 832             }
 833             else if (expand->elmcmove == elmcmoveMETGIBBS)
 834             {
 835
 836                 /* Metropolized Gibbs sampling */
 837                 for (ifep = minfep; ifep <= maxfep; ifep++)
 838                 {
 839                     remainder[ifep] = 1 - p_k[ifep];
 840                 }
 841
 842                 /* find the proposal probabilities */
 843
 844                 if (remainder[fep_state] == 0)
 845                 {
 846                     /* only the current state has any probability */
 847                     /* we have to stay at the current state */
 848                     lamnew = fep_state;
 849                 }
 850                 else
 851                 {
 852                     for (ifep = minfep; ifep <= maxfep; ifep++)
 853                     {
 854                         if (ifep != fep_state)
 855                         {
 856                             propose[ifep] = p_k[ifep]/remainder[fep_state];
 857                         }
 858                         else
 859                         {
 860                             propose[ifep] = 0;
 861                         }
 862                     }
 863
 864                     r1 = rnd[0];
 865                     for (lamtrial = minfep; lamtrial <= maxfep; lamtrial++)
 866                     {
 867                         pnorm = p_k[lamtrial]/remainder[fep_state];
 868                         if (lamtrial != fep_state)
 869                         {
 870                             if (r1 <= pnorm)
 871                             {
 872                                 break;
 873                             }
 874                             r1 -= pnorm;
 875                         }
 876                     }
 877
 878                     /* we have now selected lamtrial according to p(lamtrial)/1-p(fep_state) */
 879                     tprob = 1.0;
 880                     /* trial probability is min{1,\frac{1 - p(old)}{1-p(new)} MRS 1/8/2008 */
 881                     trialprob = (remainder[fep_state])/(remainder[lamtrial]);
 882                     if (trialprob < tprob)
 883                     {
 884                         tprob = trialprob;
 885                     }
 886                     r2 = rnd[1];
 887                     if (r2 < tprob)
 888                     {
 889                         lamnew = lamtrial;
 890                     }
 891                     else
 892                     {
 893                         lamnew = fep_state;
 894                     }
 895                 }
 896
 897                 /* now figure out the acceptance probability for each */
 898                 for (ifep = minfep; ifep <= maxfep; ifep++)
 899                 {
 900                     tprob = 1.0;
 901                     if (remainder[ifep] != 0)
 902                     {
 903                         trialprob = (remainder[fep_state])/(remainder[ifep]);
 904                     }
 905                     else
 906                     {
 907                         trialprob = 1.0; /* this state is the only choice! */
 908                     }
 909                     if (trialprob < tprob)
 910                     {
 911                         tprob = trialprob;
 912                     }
 913                     /* probability for fep_state=0, but that's fine, it's never proposed! */
 914                     accept[ifep] = tprob;
 915                 }
 916             }
 917
 918             if (lamnew > maxfep)
 919             {
 920                 /* it's possible some rounding is failing */
 921                 if (gmx_within_tol(remainder[fep_state], 0, 50*GMX_DOUBLE_EPS))
 922                 {
 923                     /* numerical rounding error -- no state other than the original has weight */
 924                     lamnew = fep_state;
 925                 }
 926                 else
 927                 {
 928                     /* probably not a numerical issue */
 929                     int   loc    = 0;
 930                     int   nerror = 200+(maxfep-minfep+1)*60;
 931                     char *errorstr;
 932                     snew(errorstr, nerror);
 933                     /* if its greater than maxfep, then something went wrong -- probably underflow in the calculation
 934                        of sum weights. Generated detailed info for failure */
 935                     loc += sprintf(errorstr, "Something wrong in choosing new lambda state with a Gibbs move -- probably underflow in weight determination.\nDenominator is: %3d%17.10e\n  i                dE        numerator          weights\n", 0, pks);
 936                     for (ifep = minfep; ifep <= maxfep; ifep++)
 937                     {
 938                         loc += sprintf(&errorstr[loc], "%3d %17.10e%17.10e%17.10e\n", ifep, weighted_lamee[ifep], p_k[ifep], dfhist->sum_weights[ifep]);
 939                     }
 940                     gmx_fatal(FARGS, errorstr);
 941                 }
 942             }
 943         }
 944         else if ((expand->elmcmove == elmcmoveMETROPOLIS) || (expand->elmcmove == elmcmoveBARKER))
 945         {
 946             /* use the metropolis sampler with trial +/- 1 */
 947             r1 = rnd[0];
 948             if (r1 < 0.5)
 949             {
 950                 if (fep_state == 0)
 951                 {
 952                     lamtrial = fep_state;
 953                 }
 954                 else
 955                 {
 956                     lamtrial = fep_state-1;
 957                 }
 958             }
 959             else
 960             {
 961                 if (fep_state == nlim-1)
 962                 {
 963                     lamtrial = fep_state;
 964                 }
 965                 else
 966                 {
 967                     lamtrial = fep_state+1;
 968                 }
 969             }
 970
 971             de = weighted_lamee[lamtrial] - weighted_lamee[fep_state];
 972             if (expand->elmcmove == elmcmoveMETROPOLIS)
 973             {
 974                 tprob     = 1.0;
 975                 trialprob = exp(de);
 976                 if (trialprob < tprob)
 977                 {
 978                     tprob = trialprob;
 979                 }
 980                 propose[fep_state] = 0;
 981                 propose[lamtrial]  = 1.0; /* note that this overwrites the above line if fep_state = ntrial, which only occurs at the ends */
 982                 accept[fep_state]  = 1.0; /* doesn't actually matter, never proposed unless fep_state = ntrial, in which case it's 1.0 anyway */
 983                 accept[lamtrial]   = tprob;
 984
 985             }
 986             else if (expand->elmcmove == elmcmoveBARKER)
 987             {
 988                 tprob = 1.0/(1.0+exp(-de));
 989
 990                 propose[fep_state] = (1-tprob);
 991                 propose[lamtrial] += tprob; /* we add, to account for the fact that at the end, they might be the same point */
 992                 accept[fep_state]  = 1.0;
 993                 accept[lamtrial]   = 1.0;
 994             }
 995
 996             r2 = rnd[1];
 997             if (r2 < tprob)
 998             {
 999                 lamnew = lamtrial;
1000             }
1001             else
1002             {
1003                 lamnew = fep_state;
1004             }
1005         }
1006
1007         for (ifep = 0; ifep < nlim; ifep++)
1008         {
1009             dfhist->Tij[fep_state][ifep]      += propose[ifep]*accept[ifep];
1010             dfhist->Tij[fep_state][fep_state] += propose[ifep]*(1.0-accept[ifep]);
1011         }
1012         fep_state = lamnew;
1013     }
1014
1015     dfhist->Tij_empirical[starting_fep_state][lamnew] += 1.0;
1016
1017     sfree(propose);
1018     sfree(accept);
1019     sfree(remainder);
1020
1021     return lamnew;
1022 }
1023
1024 /* print out the weights to the log, along with current state */
1025 extern void PrintFreeEnergyInfoToFile(FILE *outfile, t_lambda *fep, t_expanded *expand, t_simtemp *simtemp, df_history_t *dfhist,
1026                                       int fep_state, int frequency, gmx_int64_t step)
1027 {
1028     int         nlim, i, ifep, jfep;
1029     real        dw, dg, dv, dm, Tprint;
1030     real       *temps;
1031     const char *print_names[efptNR] = {" FEPL", "MassL", "CoulL", " VdwL", "BondL", "RestT", "Temp.(K)"};
1032     gmx_bool    bSimTemp            = FALSE;
1033
1034     nlim = fep->n_lambda;
1035     if (simtemp != NULL)
1036     {
1037         bSimTemp = TRUE;
1038     }
1039
1040     if (mod(step, frequency) == 0)
1041     {
1042         fprintf(outfile, "             MC-lambda information\n");
1043         if (EWL(expand->elamstats) && (!(dfhist->bEquil)))
1044         {
1045             fprintf(outfile, "  Wang-Landau incrementor is: %11.5g\n", dfhist->wl_delta);
1046         }
1047         fprintf(outfile, "  N");
1048         for (i = 0; i < efptNR; i++)
1049         {
1050             if (fep->separate_dvdl[i])
1051             {
1052                 fprintf(outfile, "%7s", print_names[i]);
1053             }
1054             else if ((i == efptTEMPERATURE) && bSimTemp)
1055             {
1056                 fprintf(outfile, "%10s", print_names[i]); /* more space for temperature formats */
1057             }
1058         }
1059         fprintf(outfile, "    Count   ");
1060         if (expand->elamstats == elamstatsMINVAR)
1061         {
1062             fprintf(outfile, "W(in kT)   G(in kT)  dG(in kT)  dV(in kT)\n");
1063         }
1064         else
1065         {
1066             fprintf(outfile, "G(in kT)  dG(in kT)\n");
1067         }
1068         for (ifep = 0; ifep < nlim; ifep++)
1069         {
1070             if (ifep == nlim-1)
1071             {
1072                 dw = 0.0;
1073                 dg = 0.0;
1074                 dv = 0.0;
1075                 dm = 0.0;
1076             }
1077             else
1078             {
1079                 dw = dfhist->sum_weights[ifep+1] - dfhist->sum_weights[ifep];
1080                 dg = dfhist->sum_dg[ifep+1] - dfhist->sum_dg[ifep];
1081                 dv = sqrt(pow(dfhist->sum_variance[ifep+1], 2) - pow(dfhist->sum_variance[ifep], 2));
1082                 dm = dfhist->sum_minvar[ifep+1] - dfhist->sum_minvar[ifep];
1083
1084             }
1085             fprintf(outfile, "%3d", (ifep+1));
1086             for (i = 0; i < efptNR; i++)
1087             {
1088                 if (fep->separate_dvdl[i])
1089                 {
1090                     fprintf(outfile, "%7.3f", fep->all_lambda[i][ifep]);
1091                 }
1092                 else if (i == efptTEMPERATURE && bSimTemp)
1093                 {
1094                     fprintf(outfile, "%9.3f", simtemp->temperatures[ifep]);
1095                 }
1096             }
1097             if (EWL(expand->elamstats) && (!(dfhist->bEquil)))  /* if performing WL and still haven't equilibrated */
1098             {
1099                 if (expand->elamstats == elamstatsWL)
1100                 {
1101                     fprintf(outfile, " %8d", (int)dfhist->wl_histo[ifep]);
1102                 }
1103                 else
1104                 {
1105                     fprintf(outfile, " %8.3f", dfhist->wl_histo[ifep]);
1106                 }
1107             }
1108             else   /* we have equilibrated weights */
1109             {
1110                 fprintf(outfile, " %8d", dfhist->n_at_lam[ifep]);
1111             }
1112             if (expand->elamstats == elamstatsMINVAR)
1113             {
1114                 fprintf(outfile, " %10.5f %10.5f %10.5f %10.5f", dfhist->sum_weights[ifep], dfhist->sum_dg[ifep], dg, dv);
1115             }
1116             else
1117             {
1118                 fprintf(outfile, " %10.5f %10.5f", dfhist->sum_weights[ifep], dw);
1119             }
1120             if (ifep == fep_state)
1121             {
1122                 fprintf(outfile, " <<\n");
1123             }
1124             else
1125             {
1126                 fprintf(outfile, "   \n");
1127             }
1128         }
1129         fprintf(outfile, "\n");
1130
1131         if ((mod(step, expand->nstTij) == 0) && (expand->nstTij > 0) && (step > 0))
1132         {
1133             fprintf(outfile, "                     Transition Matrix\n");
1134             for (ifep = 0; ifep < nlim; ifep++)
1135             {
1136                 fprintf(outfile, "%12d", (ifep+1));
1137             }
1138             fprintf(outfile, "\n");
1139             for (ifep = 0; ifep < nlim; ifep++)
1140             {
1141                 for (jfep = 0; jfep < nlim; jfep++)
1142                 {
1143                     if (dfhist->n_at_lam[ifep] > 0)
1144                     {
1145                         if (expand->bSymmetrizedTMatrix)
1146                         {
1147                             Tprint = (dfhist->Tij[ifep][jfep]+dfhist->Tij[jfep][ifep])/(dfhist->n_at_lam[ifep]+dfhist->n_at_lam[jfep]);
1148                         }
1149                         else
1150                         {
1151                             Tprint = (dfhist->Tij[ifep][jfep])/(dfhist->n_at_lam[ifep]);
1152                         }
1153                     }
1154                     else
1155                     {
1156                         Tprint = 0.0;
1157                     }
1158                     fprintf(outfile, "%12.8f", Tprint);
1159                 }
1160                 fprintf(outfile, "%3d\n", (ifep+1));
1161             }
1162
1163             fprintf(outfile, "                  Empirical Transition Matrix\n");
1164             for (ifep = 0; ifep < nlim; ifep++)
1165             {
1166                 fprintf(outfile, "%12d", (ifep+1));
1167             }
1168             fprintf(outfile, "\n");
1169             for (ifep = 0; ifep < nlim; ifep++)
1170             {
1171                 for (jfep = 0; jfep < nlim; jfep++)
1172                 {
1173                     if (dfhist->n_at_lam[ifep] > 0)
1174                     {
1175                         if (expand->bSymmetrizedTMatrix)
1176                         {
1177                             Tprint = (dfhist->Tij_empirical[ifep][jfep]+dfhist->Tij_empirical[jfep][ifep])/(dfhist->n_at_lam[ifep]+dfhist->n_at_lam[jfep]);
1178                         }
1179                         else
1180                         {
1181                             Tprint = dfhist->Tij_empirical[ifep][jfep]/(dfhist->n_at_lam[ifep]);
1182                         }
1183                     }
1184                     else
1185                     {
1186                         Tprint = 0.0;
1187                     }
1188                     fprintf(outfile, "%12.8f", Tprint);
1189                 }
1190                 fprintf(outfile, "%3d\n", (ifep+1));
1191             }
1192         }
1193     }
1194 }
1195
1196 extern int ExpandedEnsembleDynamics(FILE *log, t_inputrec *ir, gmx_enerdata_t *enerd,
1197                                     t_state *state, t_extmass *MassQ, int fep_state, df_history_t *dfhist,
1198                                     gmx_int64_t step,
1199                                     rvec *v, t_mdatoms *mdatoms)
1200 /* Note that the state variable is only needed for simulated tempering, not
1201    Hamiltonian expanded ensemble.  May be able to remove it after integrator refactoring. */
1202 {
1203     real       *pfep_lamee, *scaled_lamee, *weighted_lamee;
1204     double     *p_k;
1205     int         i, nlim, lamnew, totalsamples;
1206     real        oneovert, maxscaled = 0, maxweighted = 0;
1207     t_expanded *expand;
1208     t_simtemp  *simtemp;
1209     double     *temperature_lambdas;
1210     gmx_bool    bIfReset, bSwitchtoOneOverT, bDoneEquilibrating = FALSE;
1211
1212     expand  = ir->expandedvals;
1213     simtemp = ir->simtempvals;
1214     nlim    = ir->fepvals->n_lambda;
1215
1216     snew(scaled_lamee, nlim);
1217     snew(weighted_lamee, nlim);
1218     snew(pfep_lamee, nlim);
1219     snew(p_k, nlim);
1220
1221     /* update the count at the current lambda*/
1222     dfhist->n_at_lam[fep_state]++;
1223
1224     /* need to calculate the PV term somewhere, but not needed here? Not until there's a lambda state that's
1225        pressure controlled.*/
1226     /*
1227        pVTerm = 0;
1228        where does this PV term go?
1229        for (i=0;i<nlim;i++)
1230        {
1231        fep_lamee[i] += pVTerm;
1232        }
1233      */
1234
1235     /* determine the minimum value to avoid overflow.  Probably a better way to do this */
1236     /* we don't need to include the pressure term, since the volume is the same between the two.
1237        is there some term we are neglecting, however? */
1238
1239     if (ir->efep != efepNO)
1240     {
1241         for (i = 0; i < nlim; i++)
1242         {
1243             if (ir->bSimTemp)
1244             {
1245                 /* Note -- this assumes no mass changes, since kinetic energy is not added  . . . */
1246                 scaled_lamee[i] = (enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0])/(simtemp->temperatures[i]*BOLTZ)
1247                     + enerd->term[F_EPOT]*(1.0/(simtemp->temperatures[i])- 1.0/(simtemp->temperatures[fep_state]))/BOLTZ;
1248             }
1249             else
1250             {
1251                 scaled_lamee[i] = (enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0])/(expand->mc_temp*BOLTZ);
1252                 /* mc_temp is currently set to the system reft unless otherwise defined */
1253             }
1254
1255             /* save these energies for printing, so they don't get overwritten by the next step */
1256             /* they aren't overwritten in the non-free energy case, but we always print with these
1257                for simplicity */
1258         }
1259     }
1260     else
1261     {
1262         if (ir->bSimTemp)
1263         {
1264             for (i = 0; i < nlim; i++)
1265             {
1266                 scaled_lamee[i] = enerd->term[F_EPOT]*(1.0/simtemp->temperatures[i] - 1.0/simtemp->temperatures[fep_state])/BOLTZ;
1267             }
1268         }
1269     }
1270
1271     for (i = 0; i < nlim; i++)
1272     {
1273         pfep_lamee[i] = scaled_lamee[i];
1274
1275         weighted_lamee[i] = dfhist->sum_weights[i] - scaled_lamee[i];
1276         if (i == 0)
1277         {
1278             maxscaled   = scaled_lamee[i];
1279             maxweighted = weighted_lamee[i];
1280         }
1281         else
1282         {
1283             if (scaled_lamee[i] > maxscaled)
1284             {
1285                 maxscaled = scaled_lamee[i];
1286             }
1287             if (weighted_lamee[i] > maxweighted)
1288             {
1289                 maxweighted = weighted_lamee[i];
1290             }
1291         }
1292     }
1293
1294     for (i = 0; i < nlim; i++)
1295     {
1296         scaled_lamee[i]   -= maxscaled;
1297         weighted_lamee[i] -= maxweighted;
1298     }
1299
1300     /* update weights - we decide whether or not to actually do this inside */
1301
1302     bDoneEquilibrating = UpdateWeights(nlim, expand, dfhist, fep_state, scaled_lamee, weighted_lamee, step);
1303     if (bDoneEquilibrating)
1304     {
1305         if (log)
1306         {
1307             fprintf(log, "\nStep %d: Weights have equilibrated, using criteria: %s\n", (int)step, elmceq_names[expand->elmceq]);
1308         }
1309     }
1310
1311     lamnew = ChooseNewLambda(nlim, expand, dfhist, fep_state, weighted_lamee, p_k,
1312                              ir->expandedvals->lmc_seed, step);
1313     /* if using simulated tempering, we need to adjust the temperatures */
1314     if (ir->bSimTemp && (lamnew != fep_state)) /* only need to change the temperatures if we change the state */
1315     {
1316         int   i, j, n, d;
1317         real *buf_ngtc;
1318         real  told;
1319         int   nstart, nend, gt;
1320
1321         snew(buf_ngtc, ir->opts.ngtc);
1322
1323         for (i = 0; i < ir->opts.ngtc; i++)
1324         {
1325             if (ir->opts.ref_t[i] > 0)
1326             {
1327                 told              = ir->opts.ref_t[i];
1328                 ir->opts.ref_t[i] =  simtemp->temperatures[lamnew];
1329                 buf_ngtc[i]       = sqrt(ir->opts.ref_t[i]/told); /* using the buffer as temperature scaling */
1330             }
1331         }
1332
1333         /* we don't need to manipulate the ekind information, as it isn't due to be reset until the next step anyway */
1334
1335         nstart = 0;
1336         nend   = mdatoms->homenr;
1337         for (n = nstart; n < nend; n++)
1338         {
1339             gt = 0;
1340             if (mdatoms->cTC)
1341             {
1342                 gt = mdatoms->cTC[n];
1343             }
1344             for (d = 0; d < DIM; d++)
1345             {
1346                 v[n][d] *= buf_ngtc[gt];
1347             }
1348         }
1349
1350         if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir))
1351         {
1352             /* we need to recalculate the masses if the temperature has changed */
1353             init_npt_masses(ir, state, MassQ, FALSE);
1354             for (i = 0; i < state->nnhpres; i++)
1355             {
1356                 for (j = 0; j < ir->opts.nhchainlength; j++)
1357                 {
1358                     state->nhpres_vxi[i+j] *= buf_ngtc[i];
1359                 }
1360             }
1361             for (i = 0; i < ir->opts.ngtc; i++)
1362             {
1363                 for (j = 0; j < ir->opts.nhchainlength; j++)
1364                 {
1365                     state->nosehoover_vxi[i+j] *= buf_ngtc[i];
1366                 }
1367             }
1368         }
1369         sfree(buf_ngtc);
1370     }
1371
1372     /* now check on the Wang-Landau updating critera */
1373
1374     if (EWL(expand->elamstats))
1375     {
1376         bSwitchtoOneOverT = FALSE;
1377         if (expand->bWLoneovert)
1378         {
1379             totalsamples = 0;
1380             for (i = 0; i < nlim; i++)
1381             {
1382                 totalsamples += dfhist->n_at_lam[i];
1383             }
1384             oneovert = (1.0*nlim)/totalsamples;
1385             /* oneovert has decreasd by a bit since last time, so we actually make sure its within one of this number */
1386             /* switch to 1/t incrementing when wl_delta has decreased at least once, and wl_delta is now less than 1/t */
1387             if ((dfhist->wl_delta <= ((totalsamples)/(totalsamples-1.00001))*oneovert) &&
1388                 (dfhist->wl_delta < expand->init_wl_delta))
1389             {
1390                 bSwitchtoOneOverT = TRUE;
1391             }
1392         }
1393         if (bSwitchtoOneOverT)
1394         {
1395             dfhist->wl_delta = oneovert; /* now we reduce by this each time, instead of only at flatness */
1396         }
1397         else
1398         {
1399             bIfReset = CheckHistogramRatios(nlim, dfhist->wl_histo, expand->wl_ratio);
1400             if (bIfReset)
1401             {
1402                 for (i = 0; i < nlim; i++)
1403                 {
1404                     dfhist->wl_histo[i] = 0;
1405                 }
1406                 dfhist->wl_delta *= expand->wl_scale;
1407                 if (log)
1408                 {
1409                     fprintf(log, "\nStep %d: weights are now:", (int)step);
1410                     for (i = 0; i < nlim; i++)
1411                     {
1412                         fprintf(log, " %.5f", dfhist->sum_weights[i]);
1413                     }
1414                     fprintf(log, "\n");
1415                 }
1416             }
1417         }
1418     }
1419     sfree(pfep_lamee);
1420     sfree(scaled_lamee);
1421     sfree(weighted_lamee);
1422     sfree(p_k);
1423
1424     return lamnew;
1425 }