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