src/gromacs/mdlib/expanded.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "expanded.h"
  39
  40 #include <cmath>
  41 #include <cstdio>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/domdec/domdec.h"
  46 #include "gromacs/fileio/confio.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/xtcio.h"
  49 #include "gromacs/gmxlib/network.h"
  50 #include "gromacs/gmxlib/nrnb.h"
  51 #include "gromacs/listed_forces/disre.h"
  52 #include "gromacs/listed_forces/orires.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/mdlib/calcmu.h"
  57 #include "gromacs/mdlib/constr.h"
  58 #include "gromacs/mdlib/force.h"
  59 #include "gromacs/mdlib/update.h"
  60 #include "gromacs/mdtypes/enerdata.h"
  61 #include "gromacs/mdtypes/forcerec.h"
  62 #include "gromacs/mdtypes/inputrec.h"
  63 #include "gromacs/mdtypes/md_enums.h"
  64 #include "gromacs/mdtypes/mdatom.h"
  65 #include "gromacs/mdtypes/state.h"
  66 #include "gromacs/random/threefry.h"
  67 #include "gromacs/random/uniformrealdistribution.h"
  68 #include "gromacs/timing/wallcycle.h"
  69 #include "gromacs/utility/fatalerror.h"
  70 #include "gromacs/utility/gmxmpi.h"
  71 #include "gromacs/utility/smalloc.h"
  72
  73 #include "expanded_internal.h"
  74
  75 static void init_df_history_weights(df_history_t* dfhist, const t_expanded* expand, int nlim)
  76 {
  77     int i;
  78     dfhist->wl_delta = expand->init_wl_delta;
  79     for (i = 0; i < nlim; i++)
  80     {
  81         dfhist->sum_weights[i] = expand->init_lambda_weights[i];
  82         dfhist->sum_dg[i]      = expand->init_lambda_weights[i];
  83     }
  84 }
  85
  86 /* Eventually should contain all the functions needed to initialize expanded ensemble
  87    before the md loop starts */
  88 void init_expanded_ensemble(gmx_bool bStateFromCP, const t_inputrec* ir, df_history_t* dfhist)
  89 {
  90     if (!bStateFromCP)
  91     {
  92         init_df_history_weights(dfhist, ir->expandedvals, ir->fepvals->n_lambda);
  93     }
  94 }
  95
  96 static void GenerateGibbsProbabilities(const real* ene, double* p_k, double* pks, int minfep, int maxfep)
  97 {
  98
  99     int  i;
 100     real maxene;
 101
 102     *pks   = 0.0;
 103     maxene = ene[minfep];
 104     /* find the maximum value */
 105     for (i = minfep; i <= maxfep; i++)
 106     {
 107         if (ene[i] > maxene)
 108         {
 109             maxene = ene[i];
 110         }
 111     }
 112     /* find the denominator */
 113     for (i = minfep; i <= maxfep; i++)
 114     {
 115         *pks += std::exp(ene[i] - maxene);
 116     }
 117     /*numerators*/
 118     for (i = minfep; i <= maxfep; i++)
 119     {
 120         p_k[i] = std::exp(ene[i] - maxene) / *pks;
 121     }
 122 }
 123
 124 static void
 125 GenerateWeightedGibbsProbabilities(const real* ene, double* p_k, double* pks, int nlim, real* nvals, real delta)
 126 {
 127
 128     int   i;
 129     real  maxene;
 130     real* nene;
 131     *pks = 0.0;
 132
 133     snew(nene, nlim);
 134     for (i = 0; i < nlim; i++)
 135     {
 136         if (nvals[i] == 0)
 137         {
 138             /* add the delta, since we need to make sure it's greater than zero, and
 139                we need a non-arbitrary number? */
 140             nene[i] = ene[i] + std::log(nvals[i] + delta);
 141         }
 142         else
 143         {
 144             nene[i] = ene[i] + std::log(nvals[i]);
 145         }
 146     }
 147
 148     /* find the maximum value */
 149     maxene = nene[0];
 150     for (i = 0; i < nlim; i++)
 151     {
 152         if (nene[i] > maxene)
 153         {
 154             maxene = nene[i];
 155         }
 156     }
 157
 158     /* subtract off the maximum, avoiding overflow */
 159     for (i = 0; i < nlim; i++)
 160     {
 161         nene[i] -= maxene;
 162     }
 163
 164     /* find the denominator */
 165     for (i = 0; i < nlim; i++)
 166     {
 167         *pks += std::exp(nene[i]);
 168     }
 169
 170     /*numerators*/
 171     for (i = 0; i < nlim; i++)
 172     {
 173         p_k[i] = std::exp(nene[i]) / *pks;
 174     }
 175     sfree(nene);
 176 }
 177
 178 static int FindMinimum(const real* min_metric, int N)
 179 {
 180
 181     real min_val;
 182     int  min_nval, nval;
 183
 184     min_nval = 0;
 185     min_val  = min_metric[0];
 186
 187     for (nval = 0; nval < N; nval++)
 188     {
 189         if (min_metric[nval] < min_val)
 190         {
 191             min_val  = min_metric[nval];
 192             min_nval = nval;
 193         }
 194     }
 195     return min_nval;
 196 }
 197
 198 static gmx_bool CheckHistogramRatios(int nhisto, const real* histo, real ratio)
 199 {
 200
 201     int      i;
 202     real     nmean;
 203     gmx_bool bIfFlat;
 204
 205     nmean = 0;
 206     for (i = 0; i < nhisto; i++)
 207     {
 208         nmean += histo[i];
 209     }
 210
 211     if (nmean == 0)
 212     {
 213         /* no samples! is bad!*/
 214         bIfFlat = FALSE;
 215         return bIfFlat;
 216     }
 217     nmean /= static_cast<real>(nhisto);
 218
 219     bIfFlat = TRUE;
 220     for (i = 0; i < nhisto; i++)
 221     {
 222         /* make sure that all points are in the ratio < x <  1/ratio range  */
 223         if (!((histo[i] / nmean < 1.0 / ratio) && (histo[i] / nmean > ratio)))
 224         {
 225             bIfFlat = FALSE;
 226             break;
 227         }
 228     }
 229     return bIfFlat;
 230 }
 231
 232 static gmx_bool CheckIfDoneEquilibrating(int nlim, const t_expanded* expand, const df_history_t* dfhist, int64_t step)
 233 {
 234
 235     int      i, totalsamples;
 236     gmx_bool bDoneEquilibrating = TRUE;
 237     gmx_bool bIfFlat;
 238
 239     /* If we are doing slow growth to get initial values, we haven't finished equilibrating */
 240     if (expand->lmc_forced_nstart > 0)
 241     {
 242         for (i = 0; i < nlim; i++)
 243         {
 244             if (dfhist->n_at_lam[i]
 245                 < expand->lmc_forced_nstart) /* we are still doing the initial sweep, so we're
 246                                                 definitely not done equilibrating*/
 247             {
 248                 bDoneEquilibrating = FALSE;
 249                 break;
 250             }
 251         }
 252     }
 253     else
 254     {
 255         /* assume we have equilibrated the weights, then check to see if any of the conditions are not met */
 256         bDoneEquilibrating = TRUE;
 257
 258         /* calculate the total number of samples */
 259         switch (expand->elmceq)
 260         {
 261             case elmceqNO:
 262                 /* We have not equilibrated, and won't, ever. */
 263                 bDoneEquilibrating = FALSE;
 264                 break;
 265             case elmceqYES:
 266                 /* we have equilibrated -- we're done */
 267                 bDoneEquilibrating = TRUE;
 268                 break;
 269             case elmceqSTEPS:
 270                 /* first, check if we are equilibrating by steps, if we're still under */
 271                 if (step < expand->equil_steps)
 272                 {
 273                     bDoneEquilibrating = FALSE;
 274                 }
 275                 break;
 276             case elmceqSAMPLES:
 277                 totalsamples = 0;
 278                 for (i = 0; i < nlim; i++)
 279                 {
 280                     totalsamples += dfhist->n_at_lam[i];
 281                 }
 282                 if (totalsamples < expand->equil_samples)
 283                 {
 284                     bDoneEquilibrating = FALSE;
 285                 }
 286                 break;
 287             case elmceqNUMATLAM:
 288                 for (i = 0; i < nlim; i++)
 289                 {
 290                     if (dfhist->n_at_lam[i]
 291                         < expand->equil_n_at_lam) /* we are still doing the initial sweep, so we're
 292                                                      definitely not 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
 318                        to convert to 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                 break;
 334             default: bDoneEquilibrating = TRUE; break;
 335         }
 336     }
 337     return bDoneEquilibrating;
 338 }
 339
 340 static gmx_bool UpdateWeights(int           nlim,
 341                               t_expanded*   expand,
 342                               df_history_t* dfhist,
 343                               int           fep_state,
 344                               const real*   scaled_lamee,
 345                               const real*   weighted_lamee,
 346                               int64_t       step)
 347 {
 348     gmx_bool bSufficientSamples;
 349     real     acceptanceWeight;
 350     int      i;
 351     int      min_nvalm, min_nvalp, maxc;
 352     real     omega_m1_0, omega_p1_0;
 353     real     zero_sum_weights;
 354     real *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array,
 355             *dwp_array, *dwm_array;
 356     real    clam_varm, clam_varp, clam_osum, clam_weightsm, clam_weightsp, clam_minvar;
 357     real *  lam_variance, *lam_dg;
 358     double* p_k;
 359     double  pks = 0;
 360
 361     /* Future potential todos for this function (see #3848):
 362      *  - Update the names in the dhist structure to be clearer. Not done for now since this
 363      *    a bugfix update and we are mininizing other code changes.
 364      *  - Modularize the code some more.
 365      *  - potentially merge with accelerated weight histogram functionality, since it's very similar.
 366      */
 367     /*  if we have equilibrated the expanded ensemble weights, we are not updating them, so exit now */
 368     if (dfhist->bEquil)
 369     {
 370         return FALSE;
 371     }
 372
 373     if (CheckIfDoneEquilibrating(nlim, expand, dfhist, step))
 374     {
 375         dfhist->bEquil = TRUE;
 376         /* zero out the visited states so we know how many equilibrated states we have
 377            from here on out.*/
 378         for (i = 0; i < nlim; i++)
 379         {
 380             dfhist->n_at_lam[i] = 0;
 381         }
 382         return TRUE;
 383     }
 384
 385     /* If we reached this far, we have not equilibrated yet, keep on
 386        going resetting the weights */
 387
 388     if (EWL(expand->elamstats))
 389     {
 390         if (expand->elamstats == elamstatsWL) /* Using standard Wang-Landau for weight updates */
 391         {
 392             dfhist->sum_weights[fep_state] -= dfhist->wl_delta;
 393             dfhist->wl_histo[fep_state] += 1.0;
 394         }
 395         else if (expand->elamstats == elamstatsWWL)
 396         /* Using weighted Wang-Landau for weight updates.
 397          * Very closly equivalent to accelerated weight histogram approach
 398          * applied to expanded ensemble. */
 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] += static_cast<real>(p_k[i]);
 407             }
 408
 409             /* then increment weights (uses count) */
 410             pks = 0.0;
 411             GenerateWeightedGibbsProbabilities(
 412                     weighted_lamee, p_k, &pks, nlim, dfhist->wl_histo, dfhist->wl_delta);
 413
 414             for (i = 0; i < nlim; i++)
 415             {
 416                 dfhist->sum_weights[i] -= dfhist->wl_delta * static_cast<real>(p_k[i]);
 417             }
 418             /* Alternate definition, using logarithms. Shouldn't make very much difference! */
 419             /*
 420                real di;
 421                for (i=0;i<nlim;i++)
 422                {
 423                 di = (real)1.0 + dfhist->wl_delta*(real)p_k[i];
 424                 dfhist->sum_weights[i] -= log(di);
 425                }
 426              */
 427             sfree(p_k);
 428         }
 429
 430         zero_sum_weights = dfhist->sum_weights[0];
 431         for (i = 0; i < nlim; i++)
 432         {
 433             dfhist->sum_weights[i] -= zero_sum_weights;
 434         }
 435     }
 436
 437     if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMETROPOLIS
 438         || expand->elamstats == elamstatsMINVAR)
 439     {
 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 values of the free energy differences and the
 456          * variance in their estimates between nearby lambdas. We will
 457          * only actually update 2 of these, the state we are currently
 458          * at and the one we end up moving to
 459          */
 460
 461         for (i = 0; i < nlim - 1; i++)
 462         { /* only through the second to last */
 463             lam_dg[i] = dfhist->sum_dg[i + 1] - dfhist->sum_dg[i];
 464             lam_variance[i] =
 465                     gmx::square(dfhist->sum_variance[i + 1]) - gmx::square(dfhist->sum_variance[i]);
 466         }
 467
 468         /* accumulate running averages of thermodynamic averages for Bennett Acceptance Ratio-based
 469          * estimates of the free energy .
 470          * Rather than peforming self-consistent estimation of the free energies at each step,
 471          * we keep track of an array of possible different free energies (cnvals),
 472          * and we self-consistently choose the best one. The one that leads to a free energy estimate
 473          * that is closest to itself is the best estimate of the free energy.  It is essentially a
 474          * parallellized version of self-consistent iteration.  maxc is the number of these constants. */
 475
 476         for (int nval = 0; nval < maxc; nval++)
 477         {
 478             const real cnval = static_cast<real>(nval - expand->c_range);
 479
 480             /* Compute acceptance criterion weight to the state below this one for use in averages.
 481              * Note we do not have to have just moved from that state to use this free energy
 482              * estimate; these are essentially "virtual" moves. */
 483
 484             if (fep_state > 0)
 485             {
 486                 const auto lambdaEnergyDifference =
 487                         cnval - (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
 488                 acceptanceWeight =
 489                         gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
 490                 dfhist->accum_m[fep_state][nval] += acceptanceWeight;
 491                 dfhist->accum_m2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
 492             }
 493
 494             // Compute acceptance criterion weight to transition to the next state
 495             if (fep_state < nlim - 1)
 496             {
 497                 const auto lambdaEnergyDifference =
 498                         -cnval + (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
 499                 acceptanceWeight =
 500                         gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
 501                 dfhist->accum_p[fep_state][nval] += acceptanceWeight;
 502                 dfhist->accum_p2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
 503             }
 504
 505             /* Determination of Metropolis transition and Barker transition weights */
 506
 507             int numObservationsCurrentState = dfhist->n_at_lam[fep_state];
 508             /* determine the number of observations above and below the current state */
 509             int numObservationsLowerState = 0;
 510             if (fep_state > 0)
 511             {
 512                 numObservationsLowerState = dfhist->n_at_lam[fep_state - 1];
 513             }
 514             int numObservationsHigherState = 0;
 515             if (fep_state < nlim - 1)
 516             {
 517                 numObservationsHigherState = dfhist->n_at_lam[fep_state + 1];
 518             }
 519
 520             /* Calculate the biases for each expanded ensemble state that minimize the total
 521              * variance, as implemented in Martinez-Veracoechea and Escobedo,
 522              * J. Phys. Chem. B 2008, 112, 8120-8128
 523              *
 524              * The variance associated with the free energy estimate between two states i and j
 525              * is calculated as
 526              *     Var(i,j) = {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} / numObservations(i->j)
 527              *              + {avg[xi(j->i)^2] / avg[xi(j->i)]^2 - 1} / numObservations(j->i)
 528              * where xi(i->j) is the acceptance factor / weight associated with moving from state i to j
 529              * As we are calculating the acceptance factor to the neighbors every time we're visiting
 530              * a state, numObservations(i->j) == numObservations(i) and numObservations(j->i) == numObservations(j)
 531              */
 532
 533             /* Accumulation of acceptance weight averages between the current state and the
 534              * states +1 (p1) and -1 (m1), averaged at current state (0)
 535              */
 536             real avgAcceptanceCurrentToLower  = 0;
 537             real avgAcceptanceCurrentToHigher = 0;
 538             /* Accumulation of acceptance weight averages quantities between states 0
 539              *  and states +1 and -1, squared
 540              */
 541             real avgAcceptanceCurrentToLowerSquared  = 0;
 542             real avgAcceptanceCurrentToHigherSquared = 0;
 543             /* Accumulation of free energy quantities from lower state (m1) to current state (0) and squared */
 544             real avgAcceptanceLowerToCurrent        = 0;
 545             real avgAcceptanceLowerToCurrentSquared = 0;
 546             /* Accumulation of free energy quantities from upper state (p1) to current state (0) and squared */
 547             real avgAcceptanceHigherToCurrent        = 0;
 548             real avgAcceptanceHigherToCurrentSquared = 0;
 549
 550             if (numObservationsCurrentState > 0)
 551             {
 552                 avgAcceptanceCurrentToLower = dfhist->accum_m[fep_state][nval] / numObservationsCurrentState;
 553                 avgAcceptanceCurrentToHigher =
 554                         dfhist->accum_p[fep_state][nval] / numObservationsCurrentState;
 555                 avgAcceptanceCurrentToLowerSquared =
 556                         dfhist->accum_m2[fep_state][nval] / numObservationsCurrentState;
 557                 avgAcceptanceCurrentToHigherSquared =
 558                         dfhist->accum_p2[fep_state][nval] / numObservationsCurrentState;
 559             }
 560
 561             if ((fep_state > 0) && (numObservationsLowerState > 0))
 562             {
 563                 avgAcceptanceLowerToCurrent =
 564                         dfhist->accum_p[fep_state - 1][nval] / numObservationsLowerState;
 565                 avgAcceptanceLowerToCurrentSquared =
 566                         dfhist->accum_p2[fep_state - 1][nval] / numObservationsLowerState;
 567             }
 568
 569             if ((fep_state < nlim - 1) && (numObservationsHigherState > 0))
 570             {
 571                 avgAcceptanceHigherToCurrent =
 572                         dfhist->accum_m[fep_state + 1][nval] / numObservationsHigherState;
 573                 avgAcceptanceHigherToCurrentSquared =
 574                         dfhist->accum_m2[fep_state + 1][nval] / numObservationsHigherState;
 575             }
 576             /* These are accumulation of positive values (see definition of acceptance functions
 577              * above), or of squares of positive values.
 578              * We're taking this for granted in the following calculation, so make sure
 579              * here that nothing weird happened. Although technically all values should be positive,
 580              * because of floating point precisions, they might be numerically zero. */
 581             GMX_RELEASE_ASSERT(
 582                     avgAcceptanceCurrentToLower >= 0 && avgAcceptanceCurrentToLowerSquared >= 0
 583                             && avgAcceptanceCurrentToHigher >= 0
 584                             && avgAcceptanceCurrentToHigherSquared >= 0 && avgAcceptanceLowerToCurrent >= 0
 585                             && avgAcceptanceLowerToCurrentSquared >= 0 && avgAcceptanceHigherToCurrent >= 0
 586                             && avgAcceptanceHigherToCurrentSquared >= 0,
 587                     "By definition, the acceptance factors should all be nonnegative.");
 588
 589             real varianceCurrentToLower   = 0;
 590             real varianceCurrentToHigher  = 0;
 591             real weightDifferenceToLower  = 0;
 592             real weightDifferenceToHigher = 0;
 593             real varianceToLower          = 0;
 594             real varianceToHigher         = 0;
 595
 596             if (fep_state > 0)
 597             {
 598                 if (numObservationsCurrentState > 0)
 599                 {
 600                     /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
 601                      *
 602                      * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
 603                      * acceptances are all positive!), and hence
 604                      *     {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0  for  avg[xi(i->j)] -> 0
 605                      * We're catching that case explicitly to avoid numerical
 606                      * problems dividing by zero when the overlap between states is small (#3304)
 607                      */
 608                     if (avgAcceptanceCurrentToLower > 0)
 609                     {
 610                         varianceCurrentToLower =
 611                                 avgAcceptanceCurrentToLowerSquared
 612                                         / (avgAcceptanceCurrentToLower * avgAcceptanceCurrentToLower)
 613                                 - 1.0;
 614                     }
 615                     if (numObservationsLowerState > 0)
 616                     {
 617                         /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
 618                          *
 619                          * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
 620                          * acceptances are all positive!), and hence
 621                          *     {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0  for  avg[xi(i->j)] -> 0
 622                          * We're catching that case explicitly to avoid numerical
 623                          * problems dividing by zero when the overlap between states is small (#3304)
 624                          */
 625                         real varianceLowerToCurrent = 0;
 626                         if (avgAcceptanceLowerToCurrent > 0)
 627                         {
 628                             varianceLowerToCurrent =
 629                                     avgAcceptanceLowerToCurrentSquared
 630                                             / (avgAcceptanceLowerToCurrent * avgAcceptanceLowerToCurrent)
 631                                     - 1.0;
 632                         }
 633                         /* Free energy difference to the state one state lower */
 634                         /* if these either of these quantities are zero, the energies are */
 635                         /* way too large for the dynamic range.  We need an alternate guesstimate */
 636                         if ((avgAcceptanceCurrentToLower == 0) || (avgAcceptanceLowerToCurrent == 0))
 637                         {
 638                             weightDifferenceToLower =
 639                                     (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
 640                         }
 641                         else
 642                         {
 643                             weightDifferenceToLower = (std::log(avgAcceptanceCurrentToLower)
 644                                                        - std::log(avgAcceptanceLowerToCurrent))
 645                                                       + cnval;
 646                         }
 647                         /* Variance of the free energy difference to the one state lower */
 648                         varianceToLower =
 649                                 (1.0 / numObservationsCurrentState) * (varianceCurrentToLower)
 650                                 + (1.0 / numObservationsLowerState) * (varianceLowerToCurrent);
 651                     }
 652                 }
 653             }
 654
 655             if (fep_state < nlim - 1)
 656             {
 657                 if (numObservationsCurrentState > 0)
 658                 {
 659                     /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
 660                      *
 661                      * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
 662                      * acceptances are all positive!), and hence
 663                      *     {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0  for  avg[xi(i->j)] -> 0
 664                      * We're catching that case explicitly to avoid numerical
 665                      * problems dividing by zero when the overlap between states is small (#3304)
 666                      */
 667
 668                     if (avgAcceptanceCurrentToHigher < 0)
 669                     {
 670                         varianceCurrentToHigher =
 671                                 avgAcceptanceCurrentToHigherSquared
 672                                         / (avgAcceptanceCurrentToHigher * avgAcceptanceCurrentToHigher)
 673                                 - 1.0;
 674                     }
 675                     if (numObservationsHigherState > 0)
 676                     {
 677                         /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
 678                          *
 679                          * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
 680                          * acceptances are all positive!), and hence
 681                          *     {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0  for  avg[xi(i->j)] -> 0
 682                          * We're catching that case explicitly to avoid numerical
 683                          * problems dividing by zero when the overlap between states is small (#3304)
 684                          */
 685                         real varianceHigherToCurrent = 0;
 686                         if (avgAcceptanceHigherToCurrent > 0)
 687                         {
 688                             varianceHigherToCurrent =
 689                                     avgAcceptanceHigherToCurrentSquared
 690                                             / (avgAcceptanceHigherToCurrent * avgAcceptanceHigherToCurrent)
 691                                     - 1.0;
 692                         }
 693                         /* Free energy difference to the state one state higher */
 694                         /* if these either of these quantities are zero, the energies are */
 695                         /* way too large for the dynamic range.  We need an alternate guesstimate */
 696                         if ((avgAcceptanceHigherToCurrent == 0) || (avgAcceptanceCurrentToHigher == 0))
 697                         {
 698                             weightDifferenceToHigher =
 699                                     (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
 700                         }
 701                         else
 702                         {
 703                             weightDifferenceToHigher = (std::log(avgAcceptanceHigherToCurrent)
 704                                                         - std::log(avgAcceptanceCurrentToHigher))
 705                                                        + cnval;
 706                         }
 707                         /* Variance of the free energy difference to the one state higher */
 708                         varianceToHigher =
 709                                 (1.0 / numObservationsHigherState) * (varianceHigherToCurrent)
 710                                 + (1.0 / numObservationsCurrentState) * (varianceCurrentToHigher);
 711                     }
 712                 }
 713             }
 714
 715             if (numObservationsCurrentState > 0)
 716             {
 717                 omegam_array[nval] = varianceCurrentToLower;
 718             }
 719             else
 720             {
 721                 omegam_array[nval] = 0;
 722             }
 723             weightsm_array[nval] = weightDifferenceToLower;
 724             varm_array[nval]     = varianceToLower;
 725             if (numObservationsLowerState > 0)
 726             {
 727                 dwm_array[nval] =
 728                         fabs((cnval + std::log((1.0 * numObservationsCurrentState) / numObservationsLowerState))
 729                              - lam_dg[fep_state - 1]);
 730             }
 731             else
 732             {
 733                 dwm_array[nval] = std::fabs(cnval - lam_dg[fep_state - 1]);
 734             }
 735
 736             if (numObservationsCurrentState > 0)
 737             {
 738                 omegap_array[nval] = varianceCurrentToHigher;
 739             }
 740             else
 741             {
 742                 omegap_array[nval] = 0;
 743             }
 744             weightsp_array[nval] = weightDifferenceToHigher;
 745             varp_array[nval]     = varianceToHigher;
 746             if ((numObservationsHigherState > 0) && (numObservationsCurrentState > 0))
 747             {
 748                 dwp_array[nval] =
 749                         fabs((cnval + std::log((1.0 * numObservationsHigherState) / numObservationsCurrentState))
 750                              - lam_dg[fep_state]);
 751             }
 752             else
 753             {
 754                 dwp_array[nval] = std::fabs(cnval - lam_dg[fep_state]);
 755             }
 756         }
 757
 758         /* find the free energy estimate closest to the guessed weight's value */
 759
 760         min_nvalm     = FindMinimum(dwm_array, maxc);
 761         omega_m1_0    = omegam_array[min_nvalm];
 762         clam_weightsm = weightsm_array[min_nvalm];
 763         clam_varm     = varm_array[min_nvalm];
 764
 765         min_nvalp     = FindMinimum(dwp_array, maxc);
 766         omega_p1_0    = omegap_array[min_nvalp];
 767         clam_weightsp = weightsp_array[min_nvalp];
 768         clam_varp     = varp_array[min_nvalp];
 769
 770         clam_osum   = omega_m1_0 + omega_p1_0;
 771         clam_minvar = 0;
 772         if (clam_osum > 0)
 773         {
 774             clam_minvar = 0.5 * std::log(clam_osum);
 775         }
 776
 777         if (fep_state > 0)
 778         {
 779             lam_dg[fep_state - 1]       = clam_weightsm;
 780             lam_variance[fep_state - 1] = clam_varm;
 781         }
 782
 783         if (fep_state < nlim - 1)
 784         {
 785             lam_dg[fep_state]       = clam_weightsp;
 786             lam_variance[fep_state] = clam_varp;
 787         }
 788
 789         if (expand->elamstats == elamstatsMINVAR)
 790         {
 791             bSufficientSamples = TRUE;
 792             /* make sure the number of samples in each state are all
 793              * past a user-specified threshold
 794              */
 795             for (i = 0; i < nlim; i++)
 796             {
 797                 if (dfhist->n_at_lam[i] < expand->minvarmin)
 798                 {
 799                     bSufficientSamples = FALSE;
 800                 }
 801             }
 802             if (bSufficientSamples)
 803             {
 804                 dfhist->sum_minvar[fep_state] = clam_minvar;
 805                 if (fep_state == 0)
 806                 {
 807                     for (i = 0; i < nlim; i++)
 808                     {
 809                         dfhist->sum_minvar[i] += (expand->minvar_const - clam_minvar);
 810                     }
 811                     expand->minvar_const          = clam_minvar;
 812                     dfhist->sum_minvar[fep_state] = 0.0;
 813                 }
 814                 else
 815                 {
 816                     dfhist->sum_minvar[fep_state] -= expand->minvar_const;
 817                 }
 818             }
 819         }
 820
 821         /* we need to rezero minvar now, since it could change at fep_state = 0 */
 822         dfhist->sum_dg[0]       = 0.0;
 823         dfhist->sum_variance[0] = 0.0;
 824         dfhist->sum_weights[0]  = dfhist->sum_dg[0] + dfhist->sum_minvar[0]; /* should be zero */
 825
 826         for (i = 1; i < nlim; i++)
 827         {
 828             dfhist->sum_dg[i] = lam_dg[i - 1] + dfhist->sum_dg[i - 1];
 829             dfhist->sum_variance[i] =
 830                     std::sqrt(lam_variance[i - 1] + gmx::square(dfhist->sum_variance[i - 1]));
 831             dfhist->sum_weights[i] = dfhist->sum_dg[i] + dfhist->sum_minvar[i];
 832         }
 833
 834         sfree(lam_dg);
 835         sfree(lam_variance);
 836
 837         sfree(omegam_array);
 838         sfree(weightsm_array);
 839         sfree(varm_array);
 840         sfree(dwm_array);
 841
 842         sfree(omegap_array);
 843         sfree(weightsp_array);
 844         sfree(varp_array);
 845         sfree(dwp_array);
 846     }
 847     return FALSE;
 848 }
 849
 850 static int ChooseNewLambda(int               nlim,
 851                            const t_expanded* expand,
 852                            df_history_t*     dfhist,
 853                            int               fep_state,
 854                            const real*       weighted_lamee,
 855                            double*           p_k,
 856                            int64_t           seed,
 857                            int64_t           step)
 858 {
 859     /* Choose new lambda value, and update transition matrix */
 860
 861     int                  i, ifep, minfep, maxfep, lamnew, lamtrial, starting_fep_state;
 862     real                 r1, r2, de, trialprob, tprob = 0;
 863     double *             propose, *accept, *remainder;
 864     double               pks;
 865     real                 pnorm;
 866     gmx::ThreeFry2x64<0> rng(
 867             seed, gmx::RandomDomain::ExpandedEnsemble); // We only draw once, so zero bits internal counter is fine
 868     gmx::UniformRealDistribution<real> dist;
 869
 870     starting_fep_state = fep_state;
 871     lamnew             = fep_state; /* so that there is a default setting -- stays the same */
 872
 873     if (!EWL(expand->elamstats)) /* ignore equilibrating the weights if using WL */
 874     {
 875         if ((expand->lmc_forced_nstart > 0) && (dfhist->n_at_lam[nlim - 1] <= expand->lmc_forced_nstart))
 876         {
 877             /* Use a marching method to run through the lambdas and get preliminary free energy data,
 878                before starting 'free' sampling.  We start free sampling when we have enough at each lambda */
 879
 880             /* if we have enough at this lambda, move on to the next one */
 881
 882             if (dfhist->n_at_lam[fep_state] == expand->lmc_forced_nstart)
 883             {
 884                 lamnew = fep_state + 1;
 885                 if (lamnew == nlim) /* whoops, stepped too far! */
 886                 {
 887                     lamnew -= 1;
 888                 }
 889             }
 890             else
 891             {
 892                 lamnew = fep_state;
 893             }
 894             return lamnew;
 895         }
 896     }
 897
 898     snew(propose, nlim);
 899     snew(accept, nlim);
 900     snew(remainder, nlim);
 901
 902     for (i = 0; i < expand->lmc_repeats; i++)
 903     {
 904         rng.restart(step, i);
 905         dist.reset();
 906
 907         for (ifep = 0; ifep < nlim; ifep++)
 908         {
 909             propose[ifep] = 0;
 910             accept[ifep]  = 0;
 911         }
 912
 913         if ((expand->elmcmove == elmcmoveGIBBS) || (expand->elmcmove == elmcmoveMETGIBBS))
 914         {
 915             /* use the Gibbs sampler, with restricted range */
 916             if (expand->gibbsdeltalam < 0)
 917             {
 918                 minfep = 0;
 919                 maxfep = nlim - 1;
 920             }
 921             else
 922             {
 923                 minfep = fep_state - expand->gibbsdeltalam;
 924                 maxfep = fep_state + expand->gibbsdeltalam;
 925                 if (minfep < 0)
 926                 {
 927                     minfep = 0;
 928                 }
 929                 if (maxfep > nlim - 1)
 930                 {
 931                     maxfep = nlim - 1;
 932                 }
 933             }
 934
 935             GenerateGibbsProbabilities(weighted_lamee, p_k, &pks, minfep, maxfep);
 936
 937             if (expand->elmcmove == elmcmoveGIBBS)
 938             {
 939                 for (ifep = minfep; ifep <= maxfep; ifep++)
 940                 {
 941                     propose[ifep] = p_k[ifep];
 942                     accept[ifep]  = 1.0;
 943                 }
 944                 /* Gibbs sampling */
 945                 r1 = dist(rng);
 946                 for (lamnew = minfep; lamnew <= maxfep; lamnew++)
 947                 {
 948                     if (r1 <= p_k[lamnew])
 949                     {
 950                         break;
 951                     }
 952                     r1 -= p_k[lamnew];
 953                 }
 954             }
 955             else if (expand->elmcmove == elmcmoveMETGIBBS)
 956             {
 957
 958                 /* Metropolized Gibbs sampling */
 959                 for (ifep = minfep; ifep <= maxfep; ifep++)
 960                 {
 961                     remainder[ifep] = 1 - p_k[ifep];
 962                 }
 963
 964                 /* find the proposal probabilities */
 965
 966                 if (remainder[fep_state] == 0)
 967                 {
 968                     /* only the current state has any probability */
 969                     /* we have to stay at the current state */
 970                     lamnew = fep_state;
 971                 }
 972                 else
 973                 {
 974                     for (ifep = minfep; ifep <= maxfep; ifep++)
 975                     {
 976                         if (ifep != fep_state)
 977                         {
 978                             propose[ifep] = p_k[ifep] / remainder[fep_state];
 979                         }
 980                         else
 981                         {
 982                             propose[ifep] = 0;
 983                         }
 984                     }
 985
 986                     r1 = dist(rng);
 987                     for (lamtrial = minfep; lamtrial <= maxfep; lamtrial++)
 988                     {
 989                         pnorm = p_k[lamtrial] / remainder[fep_state];
 990                         if (lamtrial != fep_state)
 991                         {
 992                             if (r1 <= pnorm)
 993                             {
 994                                 break;
 995                             }
 996                             r1 -= pnorm;
 997                         }
 998                     }
 999
1000                     /* we have now selected lamtrial according to p(lamtrial)/1-p(fep_state) */
1001                     tprob = 1.0;
1002                     /* trial probability is min{1,\frac{1 - p(old)}{1-p(new)} MRS 1/8/2008 */
1003                     trialprob = (remainder[fep_state]) / (remainder[lamtrial]);
1004                     if (trialprob < tprob)
1005                     {
1006                         tprob = trialprob;
1007                     }
1008                     r2 = dist(rng);
1009                     if (r2 < tprob)
1010                     {
1011                         lamnew = lamtrial;
1012                     }
1013                     else
1014                     {
1015                         lamnew = fep_state;
1016                     }
1017                 }
1018
1019                 /* now figure out the acceptance probability for each */
1020                 for (ifep = minfep; ifep <= maxfep; ifep++)
1021                 {
1022                     tprob = 1.0;
1023                     if (remainder[ifep] != 0)
1024                     {
1025                         trialprob = (remainder[fep_state]) / (remainder[ifep]);
1026                     }
1027                     else
1028                     {
1029                         trialprob = 1.0; /* this state is the only choice! */
1030                     }
1031                     if (trialprob < tprob)
1032                     {
1033                         tprob = trialprob;
1034                     }
1035                     /* probability for fep_state=0, but that's fine, it's never proposed! */
1036                     accept[ifep] = tprob;
1037                 }
1038             }
1039
1040             if (lamnew > maxfep)
1041             {
1042                 /* it's possible some rounding is failing */
1043                 if (gmx_within_tol(remainder[fep_state], 0, 50 * GMX_DOUBLE_EPS))
1044                 {
1045                     /* numerical rounding error -- no state other than the original has weight */
1046                     lamnew = fep_state;
1047                 }
1048                 else
1049                 {
1050                     /* probably not a numerical issue */
1051                     int   loc    = 0;
1052                     int   nerror = 200 + (maxfep - minfep + 1) * 60;
1053                     char* errorstr;
1054                     snew(errorstr, nerror);
1055                     /* if its greater than maxfep, then something went wrong -- probably underflow
1056                        in the calculation of sum weights. Generated detailed info for failure */
1057                     loc += sprintf(
1058                             errorstr,
1059                             "Something wrong in choosing new lambda state with a Gibbs move -- "
1060                             "probably underflow in weight determination.\nDenominator is: "
1061                             "%3d%17.10e\n  i                dE        numerator          weights\n",
1062                             0,
1063                             pks);
1064                     for (ifep = minfep; ifep <= maxfep; ifep++)
1065                     {
1066                         loc += sprintf(&errorstr[loc],
1067                                        "%3d %17.10e%17.10e%17.10e\n",
1068                                        ifep,
1069                                        weighted_lamee[ifep],
1070                                        p_k[ifep],
1071                                        dfhist->sum_weights[ifep]);
1072                     }
1073                     gmx_fatal(FARGS, "%s", errorstr);
1074                 }
1075             }
1076         }
1077         else if ((expand->elmcmove == elmcmoveMETROPOLIS) || (expand->elmcmove == elmcmoveBARKER))
1078         {
1079             /* use the metropolis sampler with trial +/- 1 */
1080             r1 = dist(rng);
1081             if (r1 < 0.5)
1082             {
1083                 if (fep_state == 0)
1084                 {
1085                     lamtrial = fep_state;
1086                 }
1087                 else
1088                 {
1089                     lamtrial = fep_state - 1;
1090                 }
1091             }
1092             else
1093             {
1094                 if (fep_state == nlim - 1)
1095                 {
1096                     lamtrial = fep_state;
1097                 }
1098                 else
1099                 {
1100                     lamtrial = fep_state + 1;
1101                 }
1102             }
1103
1104             de = weighted_lamee[lamtrial] - weighted_lamee[fep_state];
1105             if (expand->elmcmove == elmcmoveMETROPOLIS)
1106             {
1107                 tprob = 1.0;
1108                 if (de < 0)
1109                 {
1110                     tprob = std::exp(de);
1111                 }
1112                 propose[fep_state] = 0;
1113                 propose[lamtrial]  = 1.0; /* note that this overwrites the above line if fep_state = ntrial, which only occurs at the ends */
1114                 accept[fep_state] =
1115                         1.0; /* doesn't actually matter, never proposed unless fep_state = ntrial, in which case it's 1.0 anyway */
1116                 accept[lamtrial] = tprob;
1117             }
1118             else if (expand->elmcmove == elmcmoveBARKER)
1119             {
1120                 if (de > 0) /* Numerically stable version */
1121                 {
1122                     tprob = 1.0 / (1.0 + std::exp(-de));
1123                 }
1124                 else if (de < 0)
1125                 {
1126                     tprob = std::exp(de) / (std::exp(de) + 1.0);
1127                 }
1128                 propose[fep_state] = (1 - tprob);
1129                 propose[lamtrial] +=
1130                         tprob; /* we add, to account for the fact that at the end, they might be the same point */
1131                 accept[fep_state] = 1.0;
1132                 accept[lamtrial]  = 1.0;
1133             }
1134
1135             r2 = dist(rng);
1136             if (r2 < tprob)
1137             {
1138                 lamnew = lamtrial;
1139             }
1140             else
1141             {
1142                 lamnew = fep_state;
1143             }
1144         }
1145
1146         for (ifep = 0; ifep < nlim; ifep++)
1147         {
1148             dfhist->Tij[fep_state][ifep] += propose[ifep] * accept[ifep];
1149             dfhist->Tij[fep_state][fep_state] += propose[ifep] * (1.0 - accept[ifep]);
1150         }
1151         fep_state = lamnew;
1152     }
1153
1154     dfhist->Tij_empirical[starting_fep_state][lamnew] += 1.0;
1155
1156     sfree(propose);
1157     sfree(accept);
1158     sfree(remainder);
1159
1160     return lamnew;
1161 }
1162
1163 /* print out the weights to the log, along with current state */
1164 void PrintFreeEnergyInfoToFile(FILE*               outfile,
1165                                const t_lambda*     fep,
1166                                const t_expanded*   expand,
1167                                const t_simtemp*    simtemp,
1168                                const df_history_t* dfhist,
1169                                int                 fep_state,
1170                                int                 frequency,
1171                                int64_t             step)
1172 {
1173     int         nlim, i, ifep, jfep;
1174     real        dw, dg, dv, Tprint;
1175     const char* print_names[efptNR] = { " FEPL", "MassL", "CoulL",   " VdwL",
1176                                         "BondL", "RestT", "Temp.(K)" };
1177     gmx_bool    bSimTemp            = FALSE;
1178
1179     nlim = fep->n_lambda;
1180     if (simtemp != nullptr)
1181     {
1182         bSimTemp = TRUE;
1183     }
1184
1185     if (step % frequency == 0)
1186     {
1187         fprintf(outfile, "             MC-lambda information\n");
1188         if (EWL(expand->elamstats) && (!(dfhist->bEquil)))
1189         {
1190             fprintf(outfile, "  Wang-Landau incrementor is: %11.5g\n", dfhist->wl_delta);
1191         }
1192         fprintf(outfile, "  N");
1193         for (i = 0; i < efptNR; i++)
1194         {
1195             if (fep->separate_dvdl[i])
1196             {
1197                 fprintf(outfile, "%7s", print_names[i]);
1198             }
1199             else if ((i == efptTEMPERATURE) && bSimTemp)
1200             {
1201                 fprintf(outfile, "%10s", print_names[i]); /* more space for temperature formats */
1202             }
1203         }
1204         fprintf(outfile, "    Count   ");
1205         if (expand->elamstats == elamstatsMINVAR)
1206         {
1207             fprintf(outfile, "W(in kT)   G(in kT)  dG(in kT)  dV(in kT)\n");
1208         }
1209         else
1210         {
1211             fprintf(outfile, "G(in kT)  dG(in kT)\n");
1212         }
1213         for (ifep = 0; ifep < nlim; ifep++)
1214         {
1215             if (ifep == nlim - 1)
1216             {
1217                 dw = 0.0;
1218                 dg = 0.0;
1219                 dv = 0.0;
1220             }
1221             else
1222             {
1223                 dw = dfhist->sum_weights[ifep + 1] - dfhist->sum_weights[ifep];
1224                 dg = dfhist->sum_dg[ifep + 1] - dfhist->sum_dg[ifep];
1225                 dv = std::sqrt(gmx::square(dfhist->sum_variance[ifep + 1])
1226                                - gmx::square(dfhist->sum_variance[ifep]));
1227             }
1228             fprintf(outfile, "%3d", (ifep + 1));
1229             for (i = 0; i < efptNR; i++)
1230             {
1231                 if (fep->separate_dvdl[i])
1232                 {
1233                     fprintf(outfile, "%7.3f", fep->all_lambda[i][ifep]);
1234                 }
1235                 else if (i == efptTEMPERATURE && bSimTemp)
1236                 {
1237                     fprintf(outfile, "%9.3f", simtemp->temperatures[ifep]);
1238                 }
1239             }
1240             if (EWL(expand->elamstats)
1241                 && (!(dfhist->bEquil))) /* if performing WL and still haven't equilibrated */
1242             {
1243                 if (expand->elamstats == elamstatsWL)
1244                 {
1245                     fprintf(outfile, " %8d", static_cast<int>(dfhist->wl_histo[ifep]));
1246                 }
1247                 else
1248                 {
1249                     fprintf(outfile, " %8.3f", dfhist->wl_histo[ifep]);
1250                 }
1251             }
1252             else /* we have equilibrated weights */
1253             {
1254                 fprintf(outfile, " %8d", dfhist->n_at_lam[ifep]);
1255             }
1256             if (expand->elamstats == elamstatsMINVAR)
1257             {
1258                 fprintf(outfile,
1259                         " %10.5f %10.5f %10.5f %10.5f",
1260                         dfhist->sum_weights[ifep],
1261                         dfhist->sum_dg[ifep],
1262                         dg,
1263                         dv);
1264             }
1265             else
1266             {
1267                 fprintf(outfile, " %10.5f %10.5f", dfhist->sum_weights[ifep], dw);
1268             }
1269             if (ifep == fep_state)
1270             {
1271                 fprintf(outfile, " <<\n");
1272             }
1273             else
1274             {
1275                 fprintf(outfile, "   \n");
1276             }
1277         }
1278         fprintf(outfile, "\n");
1279
1280         if ((step % expand->nstTij == 0) && (expand->nstTij > 0) && (step > 0))
1281         {
1282             fprintf(outfile, "                     Transition Matrix\n");
1283             for (ifep = 0; ifep < nlim; ifep++)
1284             {
1285                 fprintf(outfile, "%12d", (ifep + 1));
1286             }
1287             fprintf(outfile, "\n");
1288             for (ifep = 0; ifep < nlim; ifep++)
1289             {
1290                 for (jfep = 0; jfep < nlim; jfep++)
1291                 {
1292                     if (dfhist->n_at_lam[ifep] > 0)
1293                     {
1294                         if (expand->bSymmetrizedTMatrix)
1295                         {
1296                             Tprint = (dfhist->Tij[ifep][jfep] + dfhist->Tij[jfep][ifep])
1297                                      / (dfhist->n_at_lam[ifep] + dfhist->n_at_lam[jfep]);
1298                         }
1299                         else
1300                         {
1301                             Tprint = (dfhist->Tij[ifep][jfep]) / (dfhist->n_at_lam[ifep]);
1302                         }
1303                     }
1304                     else
1305                     {
1306                         Tprint = 0.0;
1307                     }
1308                     fprintf(outfile, "%12.8f", Tprint);
1309                 }
1310                 fprintf(outfile, "%3d\n", (ifep + 1));
1311             }
1312
1313             fprintf(outfile, "                  Empirical Transition Matrix\n");
1314             for (ifep = 0; ifep < nlim; ifep++)
1315             {
1316                 fprintf(outfile, "%12d", (ifep + 1));
1317             }
1318             fprintf(outfile, "\n");
1319             for (ifep = 0; ifep < nlim; ifep++)
1320             {
1321                 for (jfep = 0; jfep < nlim; jfep++)
1322                 {
1323                     if (dfhist->n_at_lam[ifep] > 0)
1324                     {
1325                         if (expand->bSymmetrizedTMatrix)
1326                         {
1327                             Tprint = (dfhist->Tij_empirical[ifep][jfep] + dfhist->Tij_empirical[jfep][ifep])
1328                                      / (dfhist->n_at_lam[ifep] + dfhist->n_at_lam[jfep]);
1329                         }
1330                         else
1331                         {
1332                             Tprint = dfhist->Tij_empirical[ifep][jfep] / (dfhist->n_at_lam[ifep]);
1333                         }
1334                     }
1335                     else
1336                     {
1337                         Tprint = 0.0;
1338                     }
1339                     fprintf(outfile, "%12.8f", Tprint);
1340                 }
1341                 fprintf(outfile, "%3d\n", (ifep + 1));
1342             }
1343         }
1344     }
1345 }
1346
1347 int ExpandedEnsembleDynamics(FILE*                 log,
1348                              t_inputrec*           ir,
1349                              const gmx_enerdata_t* enerd,
1350                              t_state*              state,
1351                              t_extmass*            MassQ,
1352                              int                   fep_state,
1353                              df_history_t*         dfhist,
1354                              int64_t               step,
1355                              rvec*                 v,
1356                              const t_mdatoms*      mdatoms)
1357 /* Note that the state variable is only needed for simulated tempering, not
1358    Hamiltonian expanded ensemble.  May be able to remove it after integrator refactoring. */
1359 {
1360     real *      pfep_lamee, *scaled_lamee, *weighted_lamee;
1361     double*     p_k;
1362     int         i, nlim, lamnew, totalsamples;
1363     real        oneovert, maxscaled = 0, maxweighted = 0;
1364     t_expanded* expand;
1365     t_simtemp*  simtemp;
1366     gmx_bool    bIfReset, bSwitchtoOneOverT, bDoneEquilibrating = FALSE;
1367
1368     expand  = ir->expandedvals;
1369     simtemp = ir->simtempvals;
1370     nlim    = ir->fepvals->n_lambda;
1371
1372     snew(scaled_lamee, nlim);
1373     snew(weighted_lamee, nlim);
1374     snew(pfep_lamee, nlim);
1375     snew(p_k, nlim);
1376
1377     /* update the count at the current lambda*/
1378     dfhist->n_at_lam[fep_state]++;
1379
1380     /* need to calculate the PV term somewhere, but not needed here? Not until there's a lambda
1381        state that's pressure controlled.*/
1382     /*
1383        pVTerm = 0;
1384        where does this PV term go?
1385        for (i=0;i<nlim;i++)
1386        {
1387        fep_lamee[i] += pVTerm;
1388        }
1389      */
1390
1391     /* determine the minimum value to avoid overflow.  Probably a better way to do this */
1392     /* we don't need to include the pressure term, since the volume is the same between the two.
1393        is there some term we are neglecting, however? */
1394
1395     if (ir->efep != efepNO)
1396     {
1397         for (i = 0; i < nlim; i++)
1398         {
1399             if (ir->bSimTemp)
1400             {
1401                 /* Note -- this assumes no mass changes, since kinetic energy is not added  . . . */
1402                 scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (simtemp->temperatures[i] * BOLTZ)
1403                                   + enerd->term[F_EPOT]
1404                                             * (1.0 / (simtemp->temperatures[i])
1405                                                - 1.0 / (simtemp->temperatures[fep_state]))
1406                                             / BOLTZ;
1407             }
1408             else
1409             {
1410                 scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (expand->mc_temp * BOLTZ);
1411                 /* mc_temp is currently set to the system reft unless otherwise defined */
1412             }
1413
1414             /* save these energies for printing, so they don't get overwritten by the next step */
1415             /* they aren't overwritten in the non-free energy case, but we always print with these
1416                for simplicity */
1417         }
1418     }
1419     else
1420     {
1421         if (ir->bSimTemp)
1422         {
1423             for (i = 0; i < nlim; i++)
1424             {
1425                 scaled_lamee[i] =
1426                         enerd->term[F_EPOT]
1427                         * (1.0 / simtemp->temperatures[i] - 1.0 / simtemp->temperatures[fep_state]) / BOLTZ;
1428             }
1429         }
1430     }
1431
1432     for (i = 0; i < nlim; i++)
1433     {
1434         pfep_lamee[i] = scaled_lamee[i];
1435
1436         weighted_lamee[i] = dfhist->sum_weights[i] - scaled_lamee[i];
1437         if (i == 0)
1438         {
1439             maxscaled   = scaled_lamee[i];
1440             maxweighted = weighted_lamee[i];
1441         }
1442         else
1443         {
1444             if (scaled_lamee[i] > maxscaled)
1445             {
1446                 maxscaled = scaled_lamee[i];
1447             }
1448             if (weighted_lamee[i] > maxweighted)
1449             {
1450                 maxweighted = weighted_lamee[i];
1451             }
1452         }
1453     }
1454
1455     for (i = 0; i < nlim; i++)
1456     {
1457         scaled_lamee[i] -= maxscaled;
1458         weighted_lamee[i] -= maxweighted;
1459     }
1460
1461     /* update weights - we decide whether or not to actually do this inside */
1462
1463     bDoneEquilibrating =
1464             UpdateWeights(nlim, expand, dfhist, fep_state, scaled_lamee, weighted_lamee, step);
1465     if (bDoneEquilibrating)
1466     {
1467         if (log)
1468         {
1469             fprintf(log,
1470                     "\nStep %" PRId64 ": Weights have equilibrated, using criteria: %s\n",
1471                     step,
1472                     elmceq_names[expand->elmceq]);
1473         }
1474     }
1475
1476     lamnew = ChooseNewLambda(
1477             nlim, expand, dfhist, fep_state, weighted_lamee, p_k, ir->expandedvals->lmc_seed, step);
1478     /* if using simulated tempering, we need to adjust the temperatures */
1479     if (ir->bSimTemp && (lamnew != fep_state)) /* only need to change the temperatures if we change the state */
1480     {
1481         int   i, j, n, d;
1482         real* buf_ngtc;
1483         real  told;
1484         int   nstart, nend, gt;
1485
1486         snew(buf_ngtc, ir->opts.ngtc);
1487
1488         for (i = 0; i < ir->opts.ngtc; i++)
1489         {
1490             if (ir->opts.ref_t[i] > 0)
1491             {
1492                 told              = ir->opts.ref_t[i];
1493                 ir->opts.ref_t[i] = simtemp->temperatures[lamnew];
1494                 buf_ngtc[i]       = std::sqrt(ir->opts.ref_t[i] / told); /* using the buffer as temperature scaling */
1495             }
1496         }
1497
1498         /* we don't need to manipulate the ekind information, as it isn't due to be reset until the next step anyway */
1499
1500         nstart = 0;
1501         nend   = mdatoms->homenr;
1502         for (n = nstart; n < nend; n++)
1503         {
1504             gt = 0;
1505             if (mdatoms->cTC)
1506             {
1507                 gt = mdatoms->cTC[n];
1508             }
1509             for (d = 0; d < DIM; d++)
1510             {
1511                 v[n][d] *= buf_ngtc[gt];
1512             }
1513         }
1514
1515         if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir))
1516         {
1517             /* we need to recalculate the masses if the temperature has changed */
1518             init_npt_masses(ir, state, MassQ, FALSE);
1519             for (i = 0; i < state->nnhpres; i++)
1520             {
1521                 for (j = 0; j < ir->opts.nhchainlength; j++)
1522                 {
1523                     state->nhpres_vxi[i + j] *= buf_ngtc[i];
1524                 }
1525             }
1526             for (i = 0; i < ir->opts.ngtc; i++)
1527             {
1528                 for (j = 0; j < ir->opts.nhchainlength; j++)
1529                 {
1530                     state->nosehoover_vxi[i + j] *= buf_ngtc[i];
1531                 }
1532             }
1533         }
1534         sfree(buf_ngtc);
1535     }
1536
1537     /* now check on the Wang-Landau updating critera */
1538
1539     if (EWL(expand->elamstats))
1540     {
1541         bSwitchtoOneOverT = FALSE;
1542         if (expand->bWLoneovert)
1543         {
1544             totalsamples = 0;
1545             for (i = 0; i < nlim; i++)
1546             {
1547                 totalsamples += dfhist->n_at_lam[i];
1548             }
1549             oneovert = (1.0 * nlim) / totalsamples;
1550             /* oneovert has decreasd by a bit since last time, so we actually make sure its within one of this number */
1551             /* switch to 1/t incrementing when wl_delta has decreased at least once, and wl_delta is now less than 1/t */
1552             if ((dfhist->wl_delta <= ((totalsamples) / (totalsamples - 1.00001)) * oneovert)
1553                 && (dfhist->wl_delta < expand->init_wl_delta))
1554             {
1555                 bSwitchtoOneOverT = TRUE;
1556             }
1557         }
1558         if (bSwitchtoOneOverT)
1559         {
1560             dfhist->wl_delta =
1561                     oneovert; /* now we reduce by this each time, instead of only at flatness */
1562         }
1563         else
1564         {
1565             bIfReset = CheckHistogramRatios(nlim, dfhist->wl_histo, expand->wl_ratio);
1566             if (bIfReset)
1567             {
1568                 for (i = 0; i < nlim; i++)
1569                 {
1570                     dfhist->wl_histo[i] = 0;
1571                 }
1572                 dfhist->wl_delta *= expand->wl_scale;
1573                 if (log)
1574                 {
1575                     fprintf(log, "\nStep %d: weights are now:", static_cast<int>(step));
1576                     for (i = 0; i < nlim; i++)
1577                     {
1578                         fprintf(log, " %.5f", dfhist->sum_weights[i]);
1579                     }
1580                     fprintf(log, "\n");
1581                 }
1582             }
1583         }
1584     }
1585     sfree(pfep_lamee);
1586     sfree(scaled_lamee);
1587     sfree(weighted_lamee);
1588     sfree(p_k);
1589
1590     return lamnew;
1591 }