src/gromacs/mdlib/expanded.c

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