src/mdlib/expanded.c

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