src/gromacs/mdrun/replicaexchange.cpp

   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-2004, The GROMACS development team.
   6  * Copyright (c) 2011,2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 /*! \internal \file
  39  *
  40  * \brief Implements the replica exchange routines.
  41  *
  42  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  43  * \author Mark Abraham <mark.j.abraham@gmail.com>
  44  * \ingroup module_mdrun
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "replicaexchange.h"
  49
  50 #include "config.h"
  51
  52 #include <cmath>
  53
  54 #include <random>
  55
  56 #include "gromacs/domdec/collect.h"
  57 #include "gromacs/gmxlib/network.h"
  58 #include "gromacs/math/units.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/mdrun/multisim.h"
  61 #include "gromacs/mdtypes/commrec.h"
  62 #include "gromacs/mdtypes/enerdata.h"
  63 #include "gromacs/mdtypes/inputrec.h"
  64 #include "gromacs/mdtypes/md_enums.h"
  65 #include "gromacs/mdtypes/state.h"
  66 #include "gromacs/random/threefry.h"
  67 #include "gromacs/random/uniformintdistribution.h"
  68 #include "gromacs/random/uniformrealdistribution.h"
  69 #include "gromacs/utility/fatalerror.h"
  70 #include "gromacs/utility/pleasecite.h"
  71 #include "gromacs/utility/smalloc.h"
  72
  73 //! Helps cut off probability values.
  74 constexpr int c_probabilityCutoff = 100;
  75
  76 /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
  77
  78 //! Rank in the multisimulation
  79 #define MSRANK(ms, nodeid)  (nodeid)
  80
  81 //! Enum for replica exchange flavours
  82 enum {
  83     ereTEMP, ereLAMBDA, ereENDSINGLE, ereTL, ereNR
  84 };
  85 /*! \brief Strings describing replica exchange flavours.
  86  *
  87  *  end_single_marker merely notes the end of single variable replica
  88  *  exchange. All types higher than it are multiple replica exchange
  89  *  methods.
  90  *
  91  * Eventually, should add 'pressure', 'temperature and pressure',
  92  *  'lambda_and_pressure', 'temperature_lambda_pressure'?; Let's wait
  93  *  until we feel better about the pressure control methods giving
  94  *  exact ensembles.  Right now, we assume constant pressure */
  95 static const char *erename[ereNR] = { "temperature", "lambda", "end_single_marker", "temperature and lambda"};
  96
  97 //! Working data for replica exchange.
  98 struct gmx_repl_ex
  99 {
 100     //! Replica ID
 101     int       repl;
 102     //! Total number of replica
 103     int       nrepl;
 104     //! Temperature
 105     real      temp;
 106     //! Replica exchange type from ere enum
 107     int       type;
 108     //! Quantity, e.g. temperature or lambda; first index is ere, second index is replica ID
 109     real    **q;
 110     //! Use constant pressure and temperature
 111     gmx_bool  bNPT;
 112     //! Replica pressures
 113     real     *pres;
 114     //! Replica indices
 115     int      *ind;
 116     //! Used for keeping track of all the replica swaps
 117     int      *allswaps;
 118     //! Replica exchange interval (number of steps)
 119     int       nst;
 120     //! Number of exchanges per interval
 121     int       nex;
 122     //! Random seed
 123     int       seed;
 124     //! Number of even and odd replica change attempts
 125     int       nattempt[2];
 126     //! Sum of probabilities
 127     real     *prob_sum;
 128     //! Number of moves between replicas i and j
 129     int     **nmoves;
 130     //! i-th element of the array is the number of exchanges between replica i-1 and i
 131     int      *nexchange;
 132
 133     /*! \brief Helper arrays for replica exchange; allocated here
 134      * so they don't have to be allocated each time */
 135     //! \{
 136     int      *destinations;
 137     int     **cyclic;
 138     int     **order;
 139     int      *tmpswap;
 140     gmx_bool *incycle;
 141     gmx_bool *bEx;
 142     //! \}
 143
 144     //! Helper arrays to hold the quantities that are exchanged.
 145     //! \{
 146     real  *prob;
 147     real  *Epot;
 148     real  *beta;
 149     real  *Vol;
 150     real **de;
 151     //! \}
 152 };
 153
 154 // TODO We should add Doxygen here some time.
 155 //! \cond
 156
 157 static gmx_bool repl_quantity(const gmx_multisim_t *ms,
 158                               struct gmx_repl_ex *re, int ere, real q)
 159 {
 160     real    *qall;
 161     gmx_bool bDiff;
 162     int      s;
 163
 164     snew(qall, ms->nsim);
 165     qall[re->repl] = q;
 166     gmx_sum_sim(ms->nsim, qall, ms);
 167
 168     bDiff = FALSE;
 169     for (s = 1; s < ms->nsim; s++)
 170     {
 171         if (qall[s] != qall[0])
 172         {
 173             bDiff = TRUE;
 174         }
 175     }
 176
 177     if (bDiff)
 178     {
 179         /* Set the replica exchange type and quantities */
 180         re->type = ere;
 181
 182         snew(re->q[ere], re->nrepl);
 183         for (s = 0; s < ms->nsim; s++)
 184         {
 185             re->q[ere][s] = qall[s];
 186         }
 187     }
 188     sfree(qall);
 189     return bDiff;
 190 }
 191
 192 gmx_repl_ex_t
 193 init_replica_exchange(FILE                            *fplog,
 194                       const gmx_multisim_t            *ms,
 195                       int                              numAtomsInSystem,
 196                       const t_inputrec                *ir,
 197                       const ReplicaExchangeParameters &replExParams)
 198 {
 199     real                pres;
 200     int                 i, j;
 201     struct gmx_repl_ex *re;
 202     gmx_bool            bTemp;
 203     gmx_bool            bLambda = FALSE;
 204
 205     fprintf(fplog, "\nInitializing Replica Exchange\n");
 206
 207     if (!isMultiSim(ms) || ms->nsim == 1)
 208     {
 209         gmx_fatal(FARGS, "Nothing to exchange with only one replica, maybe you forgot to set the -multidir option of mdrun?");
 210     }
 211     if (!EI_DYNAMICS(ir->eI))
 212     {
 213         gmx_fatal(FARGS, "Replica exchange is only supported by dynamical simulations");
 214         /* Note that PAR(cr) is defined by cr->nnodes > 1, which is
 215          * distinct from isMultiSim(ms). A multi-simulation only runs
 216          * with real MPI parallelism, but this does not imply PAR(cr)
 217          * is true!
 218          *
 219          * Since we are using a dynamical integrator, the only
 220          * decomposition is DD, so PAR(cr) and DOMAINDECOMP(cr) are
 221          * synonymous. The only way for cr->nnodes > 1 to be true is
 222          * if we are using DD. */
 223     }
 224
 225     snew(re, 1);
 226
 227     re->repl     = ms->sim;
 228     re->nrepl    = ms->nsim;
 229     snew(re->q, ereENDSINGLE);
 230
 231     fprintf(fplog, "Repl  There are %d replicas:\n", re->nrepl);
 232
 233     /* We only check that the number of atoms in the systms match.
 234      * This, of course, do not guarantee that the systems are the same,
 235      * but it does guarantee that we can perform replica exchange.
 236      */
 237     check_multi_int(fplog, ms, numAtomsInSystem, "the number of atoms", FALSE);
 238     check_multi_int(fplog, ms, ir->eI, "the integrator", FALSE);
 239     check_multi_int64(fplog, ms, ir->init_step+ir->nsteps, "init_step+nsteps", FALSE);
 240     const int nst = replExParams.exchangeInterval;
 241     check_multi_int64(fplog, ms, (ir->init_step+nst-1)/nst,
 242                       "first exchange step: init_step/-replex", FALSE);
 243     check_multi_int(fplog, ms, ir->etc, "the temperature coupling", FALSE);
 244     check_multi_int(fplog, ms, ir->opts.ngtc,
 245                     "the number of temperature coupling groups", FALSE);
 246     check_multi_int(fplog, ms, ir->epc, "the pressure coupling", FALSE);
 247     check_multi_int(fplog, ms, ir->efep, "free energy", FALSE);
 248     check_multi_int(fplog, ms, ir->fepvals->n_lambda, "number of lambda states", FALSE);
 249
 250     re->temp = ir->opts.ref_t[0];
 251     for (i = 1; (i < ir->opts.ngtc); i++)
 252     {
 253         if (ir->opts.ref_t[i] != re->temp)
 254         {
 255             fprintf(fplog, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 256             fprintf(stderr, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 257         }
 258     }
 259
 260     re->type = -1;
 261     bTemp    = repl_quantity(ms, re, ereTEMP, re->temp);
 262     if (ir->efep != efepNO)
 263     {
 264         bLambda = repl_quantity(ms, re, ereLAMBDA, static_cast<real>(ir->fepvals->init_fep_state));
 265     }
 266     if (re->type == -1)  /* nothing was assigned */
 267     {
 268         gmx_fatal(FARGS, "The properties of the %d systems are all the same, there is nothing to exchange", re->nrepl);
 269     }
 270     if (bLambda && bTemp)
 271     {
 272         re->type = ereTL;
 273     }
 274
 275     if (bTemp)
 276     {
 277         please_cite(fplog, "Sugita1999a");
 278         if (ir->epc != epcNO)
 279         {
 280             re->bNPT = TRUE;
 281             fprintf(fplog, "Repl  Using Constant Pressure REMD.\n");
 282             please_cite(fplog, "Okabe2001a");
 283         }
 284         if (ir->etc == etcBERENDSEN)
 285         {
 286             gmx_fatal(FARGS, "REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
 287                       ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
 288         }
 289     }
 290     if (bLambda)
 291     {
 292         if (ir->fepvals->delta_lambda != 0)   /* check this? */
 293         {
 294             gmx_fatal(FARGS, "delta_lambda is not zero");
 295         }
 296     }
 297     if (re->bNPT)
 298     {
 299         snew(re->pres, re->nrepl);
 300         if (ir->epct == epctSURFACETENSION)
 301         {
 302             pres = ir->ref_p[ZZ][ZZ];
 303         }
 304         else
 305         {
 306             pres = 0;
 307             j    = 0;
 308             for (i = 0; i < DIM; i++)
 309             {
 310                 if (ir->compress[i][i] != 0)
 311                 {
 312                     pres += ir->ref_p[i][i];
 313                     j++;
 314                 }
 315             }
 316             pres /= j;
 317         }
 318         re->pres[re->repl] = pres;
 319         gmx_sum_sim(re->nrepl, re->pres, ms);
 320     }
 321
 322     /* Make an index for increasing replica order */
 323     /* only makes sense if one or the other is varying, not both!
 324        if both are varying, we trust the order the person gave. */
 325     snew(re->ind, re->nrepl);
 326     for (i = 0; i < re->nrepl; i++)
 327     {
 328         re->ind[i] = i;
 329     }
 330
 331     if (re->type < ereENDSINGLE)
 332     {
 333
 334         for (i = 0; i < re->nrepl; i++)
 335         {
 336             for (j = i+1; j < re->nrepl; j++)
 337             {
 338                 if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]])
 339                 {
 340                     /* Unordered replicas are supposed to work, but there
 341                      * is still an issues somewhere.
 342                      * Note that at this point still re->ind[i]=i.
 343                      */
 344                     gmx_fatal(FARGS, "Replicas with indices %d < %d have %ss %g > %g, please order your replicas on increasing %s",
 345                               i, j,
 346                               erename[re->type],
 347                               re->q[re->type][i], re->q[re->type][j],
 348                               erename[re->type]);
 349                 }
 350                 else if (re->q[re->type][re->ind[j]] == re->q[re->type][re->ind[i]])
 351                 {
 352                     gmx_fatal(FARGS, "Two replicas have identical %ss", erename[re->type]);
 353                 }
 354             }
 355         }
 356     }
 357
 358     /* keep track of all the swaps, starting with the initial placement. */
 359     snew(re->allswaps, re->nrepl);
 360     for (i = 0; i < re->nrepl; i++)
 361     {
 362         re->allswaps[i] = re->ind[i];
 363     }
 364
 365     switch (re->type)
 366     {
 367         case ereTEMP:
 368             fprintf(fplog, "\nReplica exchange in temperature\n");
 369             for (i = 0; i < re->nrepl; i++)
 370             {
 371                 fprintf(fplog, " %5.1f", re->q[re->type][re->ind[i]]);
 372             }
 373             fprintf(fplog, "\n");
 374             break;
 375         case ereLAMBDA:
 376             fprintf(fplog, "\nReplica exchange in lambda\n");
 377             for (i = 0; i < re->nrepl; i++)
 378             {
 379                 fprintf(fplog, " %3d", static_cast<int>(re->q[re->type][re->ind[i]]));
 380             }
 381             fprintf(fplog, "\n");
 382             break;
 383         case ereTL:
 384             fprintf(fplog, "\nReplica exchange in temperature and lambda state\n");
 385             for (i = 0; i < re->nrepl; i++)
 386             {
 387                 fprintf(fplog, " %5.1f", re->q[ereTEMP][re->ind[i]]);
 388             }
 389             fprintf(fplog, "\n");
 390             for (i = 0; i < re->nrepl; i++)
 391             {
 392                 fprintf(fplog, " %5d", static_cast<int>(re->q[ereLAMBDA][re->ind[i]]));
 393             }
 394             fprintf(fplog, "\n");
 395             break;
 396         default:
 397             gmx_incons("Unknown replica exchange quantity");
 398     }
 399     if (re->bNPT)
 400     {
 401         fprintf(fplog, "\nRepl  p");
 402         for (i = 0; i < re->nrepl; i++)
 403         {
 404             fprintf(fplog, " %5.2f", re->pres[re->ind[i]]);
 405         }
 406
 407         for (i = 0; i < re->nrepl; i++)
 408         {
 409             if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
 410             {
 411                 fprintf(fplog, "\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
 412                 fprintf(stderr, "\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
 413             }
 414         }
 415     }
 416     re->nst = nst;
 417     if (replExParams.randomSeed == -1)
 418     {
 419         if (isMasterSim(ms))
 420         {
 421             re->seed = static_cast<int>(gmx::makeRandomSeed());
 422         }
 423         else
 424         {
 425             re->seed = 0;
 426         }
 427         gmx_sumi_sim(1, &(re->seed), ms);
 428     }
 429     else
 430     {
 431         re->seed = replExParams.randomSeed;
 432     }
 433     fprintf(fplog, "\nReplica exchange interval: %d\n", re->nst);
 434     fprintf(fplog, "\nReplica random seed: %d\n", re->seed);
 435
 436     re->nattempt[0] = 0;
 437     re->nattempt[1] = 0;
 438
 439     snew(re->prob_sum, re->nrepl);
 440     snew(re->nexchange, re->nrepl);
 441     snew(re->nmoves, re->nrepl);
 442     for (i = 0; i < re->nrepl; i++)
 443     {
 444         snew(re->nmoves[i], re->nrepl);
 445     }
 446     fprintf(fplog, "Replica exchange information below: ex and x = exchange, pr = probability\n");
 447
 448     /* generate space for the helper functions so we don't have to snew each time */
 449
 450     snew(re->destinations, re->nrepl);
 451     snew(re->incycle, re->nrepl);
 452     snew(re->tmpswap, re->nrepl);
 453     snew(re->cyclic, re->nrepl);
 454     snew(re->order, re->nrepl);
 455     for (i = 0; i < re->nrepl; i++)
 456     {
 457         snew(re->cyclic[i], re->nrepl+1);
 458         snew(re->order[i], re->nrepl);
 459     }
 460     /* allocate space for the functions storing the data for the replicas */
 461     /* not all of these arrays needed in all cases, but they don't take
 462        up much space, since the max size is nrepl**2 */
 463     snew(re->prob, re->nrepl);
 464     snew(re->bEx, re->nrepl);
 465     snew(re->beta, re->nrepl);
 466     snew(re->Vol, re->nrepl);
 467     snew(re->Epot, re->nrepl);
 468     snew(re->de, re->nrepl);
 469     for (i = 0; i < re->nrepl; i++)
 470     {
 471         snew(re->de[i], re->nrepl);
 472     }
 473     re->nex = replExParams.numExchanges;
 474     return re;
 475 }
 476
 477 static void exchange_reals(const gmx_multisim_t gmx_unused *ms, int gmx_unused b, real *v, int n)
 478 {
 479     real *buf;
 480     int   i;
 481
 482     if (v)
 483     {
 484         snew(buf, n);
 485 #if GMX_MPI
 486         /*
 487            MPI_Sendrecv(v,  n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 488            buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 489            ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 490          */
 491         {
 492             MPI_Request mpi_req;
 493
 494             MPI_Isend(v, n*sizeof(real), MPI_BYTE, MSRANK(ms, b), 0,
 495                       ms->mpi_comm_masters, &mpi_req);
 496             MPI_Recv(buf, n*sizeof(real), MPI_BYTE, MSRANK(ms, b), 0,
 497                      ms->mpi_comm_masters, MPI_STATUS_IGNORE);
 498             MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
 499         }
 500 #endif
 501         for (i = 0; i < n; i++)
 502         {
 503             v[i] = buf[i];
 504         }
 505         sfree(buf);
 506     }
 507 }
 508
 509
 510 static void exchange_doubles(const gmx_multisim_t gmx_unused *ms, int gmx_unused b, double *v, int n)
 511 {
 512     double *buf;
 513     int     i;
 514
 515     if (v)
 516     {
 517         snew(buf, n);
 518 #if GMX_MPI
 519         /*
 520            MPI_Sendrecv(v,  n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 521            buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 522            ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 523          */
 524         {
 525             MPI_Request mpi_req;
 526
 527             MPI_Isend(v, n*sizeof(double), MPI_BYTE, MSRANK(ms, b), 0,
 528                       ms->mpi_comm_masters, &mpi_req);
 529             MPI_Recv(buf, n*sizeof(double), MPI_BYTE, MSRANK(ms, b), 0,
 530                      ms->mpi_comm_masters, MPI_STATUS_IGNORE);
 531             MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
 532         }
 533 #endif
 534         for (i = 0; i < n; i++)
 535         {
 536             v[i] = buf[i];
 537         }
 538         sfree(buf);
 539     }
 540 }
 541
 542 static void exchange_rvecs(const gmx_multisim_t gmx_unused *ms, int gmx_unused b, rvec *v, int n)
 543 {
 544     rvec *buf;
 545     int   i;
 546
 547     if (v)
 548     {
 549         snew(buf, n);
 550 #if GMX_MPI
 551         /*
 552            MPI_Sendrecv(v[0],  n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 553            buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 554            ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 555          */
 556         {
 557             MPI_Request mpi_req;
 558
 559             MPI_Isend(v[0], n*sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0,
 560                       ms->mpi_comm_masters, &mpi_req);
 561             MPI_Recv(buf[0], n*sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0,
 562                      ms->mpi_comm_masters, MPI_STATUS_IGNORE);
 563             MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
 564         }
 565 #endif
 566         for (i = 0; i < n; i++)
 567         {
 568             copy_rvec(buf[i], v[i]);
 569         }
 570         sfree(buf);
 571     }
 572 }
 573
 574 static void exchange_state(const gmx_multisim_t *ms, int b, t_state *state)
 575 {
 576     /* When t_state changes, this code should be updated. */
 577     int ngtc, nnhpres;
 578     ngtc    = state->ngtc * state->nhchainlength;
 579     nnhpres = state->nnhpres* state->nhchainlength;
 580     exchange_rvecs(ms, b, state->box, DIM);
 581     exchange_rvecs(ms, b, state->box_rel, DIM);
 582     exchange_rvecs(ms, b, state->boxv, DIM);
 583     exchange_reals(ms, b, &(state->veta), 1);
 584     exchange_reals(ms, b, &(state->vol0), 1);
 585     exchange_rvecs(ms, b, state->svir_prev, DIM);
 586     exchange_rvecs(ms, b, state->fvir_prev, DIM);
 587     exchange_rvecs(ms, b, state->pres_prev, DIM);
 588     exchange_doubles(ms, b, state->nosehoover_xi.data(), ngtc);
 589     exchange_doubles(ms, b, state->nosehoover_vxi.data(), ngtc);
 590     exchange_doubles(ms, b, state->nhpres_xi.data(), nnhpres);
 591     exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres);
 592     exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc);
 593     exchange_doubles(ms, b, &state->baros_integral, 1);
 594     exchange_rvecs(ms, b, state->x.rvec_array(), state->natoms);
 595     exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms);
 596 }
 597
 598 static void copy_state_serial(const t_state *src, t_state *dest)
 599 {
 600     if (dest != src)
 601     {
 602         /* Currently the local state is always a pointer to the global
 603          * in serial, so we should never end up here.
 604          * TODO: Implement a (trivial) t_state copy once converted to C++.
 605          */
 606         GMX_RELEASE_ASSERT(false, "State copying is currently not implemented in replica exchange");
 607     }
 608 }
 609
 610 static void scale_velocities(gmx::ArrayRef<gmx::RVec> velocities, real fac)
 611 {
 612     for (auto &v : velocities)
 613     {
 614         v *= fac;
 615     }
 616 }
 617
 618 static void print_transition_matrix(FILE *fplog, int n, int **nmoves, const int *nattempt)
 619 {
 620     int   i, j, ntot;
 621     float Tprint;
 622
 623     ntot = nattempt[0] + nattempt[1];
 624     fprintf(fplog, "\n");
 625     fprintf(fplog, "Repl");
 626     for (i = 0; i < n; i++)
 627     {
 628         fprintf(fplog, "    ");  /* put the title closer to the center */
 629     }
 630     fprintf(fplog, "Empirical Transition Matrix\n");
 631
 632     fprintf(fplog, "Repl");
 633     for (i = 0; i < n; i++)
 634     {
 635         fprintf(fplog, "%8d", (i+1));
 636     }
 637     fprintf(fplog, "\n");
 638
 639     for (i = 0; i < n; i++)
 640     {
 641         fprintf(fplog, "Repl");
 642         for (j = 0; j < n; j++)
 643         {
 644             Tprint = 0.0;
 645             if (nmoves[i][j] > 0)
 646             {
 647                 Tprint = nmoves[i][j]/(2.0*ntot);
 648             }
 649             fprintf(fplog, "%8.4f", Tprint);
 650         }
 651         fprintf(fplog, "%3d\n", i);
 652     }
 653 }
 654
 655 static void print_ind(FILE *fplog, const char *leg, int n, int *ind, const gmx_bool *bEx)
 656 {
 657     int i;
 658
 659     fprintf(fplog, "Repl %2s %2d", leg, ind[0]);
 660     for (i = 1; i < n; i++)
 661     {
 662         fprintf(fplog, " %c %2d", (bEx != nullptr && bEx[i]) ? 'x' : ' ', ind[i]);
 663     }
 664     fprintf(fplog, "\n");
 665 }
 666
 667 static void print_allswitchind(FILE *fplog, int n, int *pind, int *allswaps, int *tmpswap)
 668 {
 669     int i;
 670
 671     for (i = 0; i < n; i++)
 672     {
 673         tmpswap[i] = allswaps[i];
 674     }
 675     for (i = 0; i < n; i++)
 676     {
 677         allswaps[i] = tmpswap[pind[i]];
 678     }
 679
 680     fprintf(fplog, "\nAccepted Exchanges:   ");
 681     for (i = 0; i < n; i++)
 682     {
 683         fprintf(fplog, "%d ", pind[i]);
 684     }
 685     fprintf(fplog, "\n");
 686
 687     /* the "Order After Exchange" is the state label corresponding to the configuration that
 688        started in state listed in order, i.e.
 689
 690        3 0 1 2
 691
 692        means that the:
 693        configuration starting in simulation 3 is now in simulation 0,
 694        configuration starting in simulation 0 is now in simulation 1,
 695        configuration starting in simulation 1 is now in simulation 2,
 696        configuration starting in simulation 2 is now in simulation 3
 697      */
 698     fprintf(fplog, "Order After Exchange: ");
 699     for (i = 0; i < n; i++)
 700     {
 701         fprintf(fplog, "%d ", allswaps[i]);
 702     }
 703     fprintf(fplog, "\n\n");
 704 }
 705
 706 static void print_prob(FILE *fplog, const char *leg, int n, real *prob)
 707 {
 708     int  i;
 709     char buf[8];
 710
 711     fprintf(fplog, "Repl %2s ", leg);
 712     for (i = 1; i < n; i++)
 713     {
 714         if (prob[i] >= 0)
 715         {
 716             sprintf(buf, "%4.2f", prob[i]);
 717             fprintf(fplog, "  %3s", buf[0] == '1' ? "1.0" : buf+1);
 718         }
 719         else
 720         {
 721             fprintf(fplog, "     ");
 722         }
 723     }
 724     fprintf(fplog, "\n");
 725 }
 726
 727 static void print_count(FILE *fplog, const char *leg, int n, int *count)
 728 {
 729     int i;
 730
 731     fprintf(fplog, "Repl %2s ", leg);
 732     for (i = 1; i < n; i++)
 733     {
 734         fprintf(fplog, " %4d", count[i]);
 735     }
 736     fprintf(fplog, "\n");
 737 }
 738
 739 static real calc_delta(FILE *fplog, gmx_bool bPrint, struct gmx_repl_ex *re, int a, int b, int ap, int bp)
 740 {
 741
 742     real   ediff, dpV, delta = 0;
 743     real  *Epot = re->Epot;
 744     real  *Vol  = re->Vol;
 745     real **de   = re->de;
 746     real  *beta = re->beta;
 747
 748     /* Two cases; we are permuted and not.  In all cases, setting ap = a and bp = b will reduce
 749        to the non permuted case */
 750
 751     switch (re->type)
 752     {
 753         case ereTEMP:
 754             /*
 755              * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
 756              */
 757             ediff = Epot[b] - Epot[a];
 758             delta = -(beta[bp] - beta[ap])*ediff;
 759             break;
 760         case ereLAMBDA:
 761             /* two cases:  when we are permuted, and not.  */
 762             /* non-permuted:
 763                ediff =  E_new - E_old
 764                      =  [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
 765                      =  [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
 766                      =  de[b][a] + de[a][b] */
 767
 768             /* permuted:
 769                ediff =  E_new - E_old
 770                      =  [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
 771                      =  [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
 772                      =  [H_bp(x_a) - H_a(x_a) + H_a(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_b(x_b) + H_b(x_b) - H_bp(x_b)]
 773                      =  [H_bp(x_a) - H_a(x_a)] - [H_ap(x_a) - H_a(x_a)] + [H_ap(x_b) - H_b(x_b)] - H_bp(x_b) - H_b(x_b)]
 774                      =  (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b])    */
 775             /* but, in the current code implementation, we flip configurations, not indices . . .
 776                So let's examine that.
 777                      =  [H_b(x_ap) - H_a(x_a)] - [H_a(x_ap) - H_a(x_a)] + [H_a(x_bp) - H_b(x_b)] - H_b(x_bp) - H_b(x_b)]
 778                      =  [H_b(x_ap) - H_a(x_ap)]  + [H_a(x_bp) - H_b(x_pb)]
 779                      = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp]
 780                      So, if we exchange b<=> bp and a<=> ap, we return to the same result.
 781                      So the simple solution is to flip the
 782                      position of perturbed and original indices in the tests.
 783              */
 784
 785             ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]);
 786             delta = ediff*beta[a]; /* assume all same temperature in this case */
 787             break;
 788         case ereTL:
 789             /* not permuted:  */
 790             /* delta =  reduced E_new - reduced E_old
 791                      =  [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
 792                      =  [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
 793                      =  [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
 794                         [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
 795                      =  [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
 796                         beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
 797                      =  beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
 798             /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
 799             /* permuted (big breath!) */
 800             /*   delta =  reduced E_new - reduced E_old
 801                      =  [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
 802                      =  [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
 803                      =  [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
 804                         - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
 805                         - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
 806                      =  [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
 807                         [(beta_ap H_ap(x_b)  - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
 808              + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
 809                      =  [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
 810                         [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
 811              + beta_pb (H_a(x_a) - H_b(x_b))  - beta_ap (H_a(x_a) - H_b(x_b))
 812                      =  ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b]  - beta_bp de[bp][b])
 813              + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b))  */
 814             delta = beta[bp]*(de[bp][a] - de[bp][b]) + beta[ap]*(de[ap][b] - de[ap][a]) - (beta[bp]-beta[ap])*(Epot[b]-Epot[a]);
 815             break;
 816         default:
 817             gmx_incons("Unknown replica exchange quantity");
 818     }
 819     if (bPrint)
 820     {
 821         fprintf(fplog, "Repl %d <-> %d  dE_term = %10.3e (kT)\n", a, b, delta);
 822     }
 823     if (re->bNPT)
 824     {
 825         /* revist the calculation for 5.0.  Might be some improvements. */
 826         dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
 827         if (bPrint)
 828         {
 829             fprintf(fplog, "  dpV = %10.3e  d = %10.3e\n", dpV, delta + dpV);
 830         }
 831         delta += dpV;
 832     }
 833     return delta;
 834 }
 835
 836 static void
 837 test_for_replica_exchange(FILE                 *fplog,
 838                           const gmx_multisim_t *ms,
 839                           struct gmx_repl_ex   *re,
 840                           const gmx_enerdata_t *enerd,
 841                           real                  vol,
 842                           int64_t               step,
 843                           real                  time)
 844 {
 845     int                                  m, i, j, a, b, ap, bp, i0, i1, tmp;
 846     real                                 delta = 0;
 847     gmx_bool                             bPrint, bMultiEx;
 848     gmx_bool                            *bEx      = re->bEx;
 849     real                                *prob     = re->prob;
 850     int                                 *pind     = re->destinations; /* permuted index */
 851     gmx_bool                             bEpot    = FALSE;
 852     gmx_bool                             bDLambda = FALSE;
 853     gmx_bool                             bVol     = FALSE;
 854     gmx::ThreeFry2x64<64>                rng(re->seed, gmx::RandomDomain::ReplicaExchange);
 855     gmx::UniformRealDistribution<real>   uniformRealDist;
 856     gmx::UniformIntDistribution<int>     uniformNreplDist(0, re->nrepl-1);
 857
 858     bMultiEx = (re->nex > 1);  /* multiple exchanges at each state */
 859     fprintf(fplog, "Replica exchange at step %" PRId64 " time %.5f\n", step, time);
 860
 861     if (re->bNPT)
 862     {
 863         for (i = 0; i < re->nrepl; i++)
 864         {
 865             re->Vol[i] = 0;
 866         }
 867         bVol               = TRUE;
 868         re->Vol[re->repl]  = vol;
 869     }
 870     if ((re->type == ereTEMP || re->type == ereTL))
 871     {
 872         for (i = 0; i < re->nrepl; i++)
 873         {
 874             re->Epot[i] = 0;
 875         }
 876         bEpot              = TRUE;
 877         re->Epot[re->repl] = enerd->term[F_EPOT];
 878         /* temperatures of different states*/
 879         for (i = 0; i < re->nrepl; i++)
 880         {
 881             re->beta[i] = 1.0/(re->q[ereTEMP][i]*BOLTZ);
 882         }
 883     }
 884     else
 885     {
 886         for (i = 0; i < re->nrepl; i++)
 887         {
 888             re->beta[i] = 1.0/(re->temp*BOLTZ);  /* we have a single temperature */
 889         }
 890     }
 891     if (re->type == ereLAMBDA || re->type == ereTL)
 892     {
 893         bDLambda = TRUE;
 894         /* lambda differences. */
 895         /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
 896            minus the energy of the jth simulation in the jth Hamiltonian */
 897         for (i = 0; i < re->nrepl; i++)
 898         {
 899             for (j = 0; j < re->nrepl; j++)
 900             {
 901                 re->de[i][j] = 0;
 902             }
 903         }
 904         for (i = 0; i < re->nrepl; i++)
 905         {
 906             re->de[i][re->repl] = (enerd->enerpart_lambda[static_cast<int>(re->q[ereLAMBDA][i])+1]-enerd->enerpart_lambda[0]);
 907         }
 908     }
 909
 910     /* now actually do the communication */
 911     if (bVol)
 912     {
 913         gmx_sum_sim(re->nrepl, re->Vol, ms);
 914     }
 915     if (bEpot)
 916     {
 917         gmx_sum_sim(re->nrepl, re->Epot, ms);
 918     }
 919     if (bDLambda)
 920     {
 921         for (i = 0; i < re->nrepl; i++)
 922         {
 923             gmx_sum_sim(re->nrepl, re->de[i], ms);
 924         }
 925     }
 926
 927     /* make a duplicate set of indices for shuffling */
 928     for (i = 0; i < re->nrepl; i++)
 929     {
 930         pind[i] = re->ind[i];
 931     }
 932
 933     rng.restart( step, 0 );
 934
 935     if (bMultiEx)
 936     {
 937         /* multiple random switch exchange */
 938         int nself = 0;
 939
 940
 941         for (i = 0; i < re->nex + nself; i++)
 942         {
 943             // For now this is superfluous, but just in case we ever add more
 944             // calls in different branches it is safer to always reset the distribution.
 945             uniformNreplDist.reset();
 946
 947             /* randomly select a pair  */
 948             /* in theory, could reduce this by identifying only which switches had a nonneglibible
 949                probability of occurring (log p > -100) and only operate on those switches */
 950             /* find out which state it is from, and what label that state currently has. Likely
 951                more work that useful. */
 952             i0 = uniformNreplDist(rng);
 953             i1 = uniformNreplDist(rng);
 954             if (i0 == i1)
 955             {
 956                 nself++;
 957                 continue;  /* self-exchange, back up and do it again */
 958             }
 959
 960             a  = re->ind[i0]; /* what are the indices of these states? */
 961             b  = re->ind[i1];
 962             ap = pind[i0];
 963             bp = pind[i1];
 964
 965             bPrint = FALSE; /* too noisy */
 966             /* calculate the energy difference */
 967             /* if the code changes to flip the STATES, rather than the configurations,
 968                use the commented version of the code */
 969             /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */
 970             delta = calc_delta(fplog, bPrint, re, ap, bp, a, b);
 971
 972             /* we actually only use the first space in the prob and bEx array,
 973                since there are actually many switches between pairs. */
 974
 975             if (delta <= 0)
 976             {
 977                 /* accepted */
 978                 prob[0] = 1;
 979                 bEx[0]  = TRUE;
 980             }
 981             else
 982             {
 983                 if (delta > c_probabilityCutoff)
 984                 {
 985                     prob[0] = 0;
 986                 }
 987                 else
 988                 {
 989                     prob[0] = exp(-delta);
 990                 }
 991                 // roll a number to determine if accepted. For now it is superfluous to
 992                 // reset, but just in case we ever add more calls in different branches
 993                 // it is safer to always reset the distribution.
 994                 uniformRealDist.reset();
 995                 bEx[0] = uniformRealDist(rng) < prob[0];
 996             }
 997             re->prob_sum[0] += prob[0];
 998
 999             if (bEx[0])
1000             {
1001                 /* swap the states */
1002                 tmp      = pind[i0];
1003                 pind[i0] = pind[i1];
1004                 pind[i1] = tmp;
1005             }
1006         }
1007         re->nattempt[0]++;  /* keep track of total permutation trials here */
1008         print_allswitchind(fplog, re->nrepl, pind, re->allswaps, re->tmpswap);
1009     }
1010     else
1011     {
1012         /* standard nearest neighbor replica exchange */
1013
1014         m = (step / re->nst) % 2;
1015         for (i = 1; i < re->nrepl; i++)
1016         {
1017             a = re->ind[i-1];
1018             b = re->ind[i];
1019
1020             bPrint = (re->repl == a || re->repl == b);
1021             if (i % 2 == m)
1022             {
1023                 delta = calc_delta(fplog, bPrint, re, a, b, a, b);
1024                 if (delta <= 0)
1025                 {
1026                     /* accepted */
1027                     prob[i] = 1;
1028                     bEx[i]  = TRUE;
1029                 }
1030                 else
1031                 {
1032                     if (delta > c_probabilityCutoff)
1033                     {
1034                         prob[i] = 0;
1035                     }
1036                     else
1037                     {
1038                         prob[i] = exp(-delta);
1039                     }
1040                     // roll a number to determine if accepted. For now it is superfluous to
1041                     // reset, but just in case we ever add more calls in different branches
1042                     // it is safer to always reset the distribution.
1043                     uniformRealDist.reset();
1044                     bEx[i] = uniformRealDist(rng) < prob[i];
1045                 }
1046                 re->prob_sum[i] += prob[i];
1047
1048                 if (bEx[i])
1049                 {
1050                     /* swap these two */
1051                     tmp       = pind[i-1];
1052                     pind[i-1] = pind[i];
1053                     pind[i]   = tmp;
1054                     re->nexchange[i]++;  /* statistics for back compatibility */
1055                 }
1056             }
1057             else
1058             {
1059                 prob[i] = -1;
1060                 bEx[i]  = FALSE;
1061             }
1062         }
1063         /* print some statistics */
1064         print_ind(fplog, "ex", re->nrepl, re->ind, bEx);
1065         print_prob(fplog, "pr", re->nrepl, prob);
1066         fprintf(fplog, "\n");
1067         re->nattempt[m]++;
1068     }
1069
1070     /* record which moves were made and accepted */
1071     for (i = 0; i < re->nrepl; i++)
1072     {
1073         re->nmoves[re->ind[i]][pind[i]] += 1;
1074         re->nmoves[pind[i]][re->ind[i]] += 1;
1075     }
1076     fflush(fplog); /* make sure we can see what the last exchange was */
1077 }
1078
1079 static void
1080 cyclic_decomposition(const int *destinations,
1081                      int      **cyclic,
1082                      gmx_bool  *incycle,
1083                      const int  nrepl,
1084                      int       *nswap)
1085 {
1086
1087     int i, j, c, p;
1088     int maxlen = 1;
1089     for (i = 0; i < nrepl; i++)
1090     {
1091         incycle[i] = FALSE;
1092     }
1093     for (i = 0; i < nrepl; i++)  /* one cycle for each replica */
1094     {
1095         if (incycle[i])
1096         {
1097             cyclic[i][0] = -1;
1098             continue;
1099         }
1100         cyclic[i][0] = i;
1101         incycle[i]   = TRUE;
1102         c            = 1;
1103         p            = i;
1104         for (j = 0; j < nrepl; j++) /* potentially all cycles are part, but we will break first */
1105         {
1106             p = destinations[p];    /* start permuting */
1107             if (p == i)
1108             {
1109                 cyclic[i][c] = -1;
1110                 if (c > maxlen)
1111                 {
1112                     maxlen = c;
1113                 }
1114                 break; /* we've reached the original element, the cycle is complete, and we marked the end. */
1115             }
1116             else
1117             {
1118                 cyclic[i][c] = p;  /* each permutation gives a new member of the cycle */
1119                 incycle[p]   = TRUE;
1120                 c++;
1121             }
1122         }
1123     }
1124     *nswap = maxlen - 1;
1125
1126     if (debug)
1127     {
1128         for (i = 0; i < nrepl; i++)
1129         {
1130             fprintf(debug, "Cycle %d:", i);
1131             for (j = 0; j < nrepl; j++)
1132             {
1133                 if (cyclic[i][j] < 0)
1134                 {
1135                     break;
1136                 }
1137                 fprintf(debug, "%2d", cyclic[i][j]);
1138             }
1139             fprintf(debug, "\n");
1140         }
1141         fflush(debug);
1142     }
1143 }
1144
1145 static void
1146 compute_exchange_order(int     **cyclic,
1147                        int     **order,
1148                        const int nrepl,
1149                        const int maxswap)
1150 {
1151     int i, j;
1152
1153     for (j = 0; j < maxswap; j++)
1154     {
1155         for (i = 0; i < nrepl; i++)
1156         {
1157             if (cyclic[i][j+1] >= 0)
1158             {
1159                 order[cyclic[i][j+1]][j] = cyclic[i][j];
1160                 order[cyclic[i][j]][j]   = cyclic[i][j+1];
1161             }
1162         }
1163         for (i = 0; i < nrepl; i++)
1164         {
1165             if (order[i][j] < 0)
1166             {
1167                 order[i][j] = i; /* if it's not exchanging, it should stay this round*/
1168             }
1169         }
1170     }
1171
1172     if (debug)
1173     {
1174         fprintf(debug, "Replica Exchange Order\n");
1175         for (i = 0; i < nrepl; i++)
1176         {
1177             fprintf(debug, "Replica %d:", i);
1178             for (j = 0; j < maxswap; j++)
1179             {
1180                 if (order[i][j] < 0)
1181                 {
1182                     break;
1183                 }
1184                 fprintf(debug, "%2d", order[i][j]);
1185             }
1186             fprintf(debug, "\n");
1187         }
1188         fflush(debug);
1189     }
1190 }
1191
1192 static void
1193 prepare_to_do_exchange(struct gmx_repl_ex *re,
1194                        const int           replica_id,
1195                        int                *maxswap,
1196                        gmx_bool           *bThisReplicaExchanged)
1197 {
1198     int i, j;
1199     /* Hold the cyclic decomposition of the (multiple) replica
1200      * exchange. */
1201     gmx_bool bAnyReplicaExchanged = FALSE;
1202     *bThisReplicaExchanged = FALSE;
1203
1204     for (i = 0; i < re->nrepl; i++)
1205     {
1206         if (re->destinations[i] != re->ind[i])
1207         {
1208             /* only mark as exchanged if the index has been shuffled */
1209             bAnyReplicaExchanged = TRUE;
1210             break;
1211         }
1212     }
1213     if (bAnyReplicaExchanged)
1214     {
1215         /* reinitialize the placeholder arrays */
1216         for (i = 0; i < re->nrepl; i++)
1217         {
1218             for (j = 0; j < re->nrepl; j++)
1219             {
1220                 re->cyclic[i][j] = -1;
1221                 re->order[i][j]  = -1;
1222             }
1223         }
1224
1225         /* Identify the cyclic decomposition of the permutation (very
1226          * fast if neighbor replica exchange). */
1227         cyclic_decomposition(re->destinations, re->cyclic, re->incycle, re->nrepl, maxswap);
1228
1229         /* Now translate the decomposition into a replica exchange
1230          * order at each step. */
1231         compute_exchange_order(re->cyclic, re->order, re->nrepl, *maxswap);
1232
1233         /* Did this replica do any exchange at any point? */
1234         for (j = 0; j < *maxswap; j++)
1235         {
1236             if (replica_id != re->order[replica_id][j])
1237             {
1238                 *bThisReplicaExchanged = TRUE;
1239                 break;
1240             }
1241         }
1242     }
1243 }
1244
1245 gmx_bool replica_exchange(FILE *fplog, const t_commrec *cr,
1246                           const gmx_multisim_t *ms, struct gmx_repl_ex *re,
1247                           t_state *state, const gmx_enerdata_t *enerd,
1248                           t_state *state_local, int64_t step, real time)
1249 {
1250     int j;
1251     int replica_id = 0;
1252     int exchange_partner;
1253     int maxswap = 0;
1254     /* Number of rounds of exchanges needed to deal with any multiple
1255      * exchanges. */
1256     /* Where each replica ends up after the exchange attempt(s). */
1257     /* The order in which multiple exchanges will occur. */
1258     gmx_bool bThisReplicaExchanged = FALSE;
1259
1260     if (MASTER(cr))
1261     {
1262         replica_id  = re->repl;
1263         test_for_replica_exchange(fplog, ms, re, enerd, det(state_local->box), step, time);
1264         prepare_to_do_exchange(re, replica_id, &maxswap, &bThisReplicaExchanged);
1265     }
1266     /* Do intra-simulation broadcast so all processors belonging to
1267      * each simulation know whether they need to participate in
1268      * collecting the state. Otherwise, they might as well get on with
1269      * the next thing to do. */
1270     if (DOMAINDECOMP(cr))
1271     {
1272 #if GMX_MPI
1273         MPI_Bcast(&bThisReplicaExchanged, sizeof(gmx_bool), MPI_BYTE, MASTERRANK(cr),
1274                   cr->mpi_comm_mygroup);
1275 #endif
1276     }
1277
1278     if (bThisReplicaExchanged)
1279     {
1280         /* Exchange the states */
1281         /* Collect the global state on the master node */
1282         if (DOMAINDECOMP(cr))
1283         {
1284             dd_collect_state(cr->dd, state_local, state);
1285         }
1286         else
1287         {
1288             copy_state_serial(state_local, state);
1289         }
1290
1291         if (MASTER(cr))
1292         {
1293             /* There will be only one swap cycle with standard replica
1294              * exchange, but there may be multiple swap cycles if we
1295              * allow multiple swaps. */
1296
1297             for (j = 0; j < maxswap; j++)
1298             {
1299                 exchange_partner = re->order[replica_id][j];
1300
1301                 if (exchange_partner != replica_id)
1302                 {
1303                     /* Exchange the global states between the master nodes */
1304                     if (debug)
1305                     {
1306                         fprintf(debug, "Exchanging %d with %d\n", replica_id, exchange_partner);
1307                     }
1308                     exchange_state(ms, exchange_partner, state);
1309                 }
1310             }
1311             /* For temperature-type replica exchange, we need to scale
1312              * the velocities. */
1313             if (re->type == ereTEMP || re->type == ereTL)
1314             {
1315                 scale_velocities(state->v,
1316                                  std::sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));
1317             }
1318
1319         }
1320
1321         /* With domain decomposition the global state is distributed later */
1322         if (!DOMAINDECOMP(cr))
1323         {
1324             /* Copy the global state to the local state data structure */
1325             copy_state_serial(state, state_local);
1326         }
1327     }
1328
1329     return bThisReplicaExchanged;
1330 }
1331
1332 void print_replica_exchange_statistics(FILE *fplog, struct gmx_repl_ex *re)
1333 {
1334     int  i;
1335
1336     fprintf(fplog, "\nReplica exchange statistics\n");
1337
1338     if (re->nex == 0)
1339     {
1340         fprintf(fplog, "Repl  %d attempts, %d odd, %d even\n",
1341                 re->nattempt[0]+re->nattempt[1], re->nattempt[1], re->nattempt[0]);
1342
1343         fprintf(fplog, "Repl  average probabilities:\n");
1344         for (i = 1; i < re->nrepl; i++)
1345         {
1346             if (re->nattempt[i%2] == 0)
1347             {
1348                 re->prob[i] = 0;
1349             }
1350             else
1351             {
1352                 re->prob[i] =  re->prob_sum[i]/re->nattempt[i%2];
1353             }
1354         }
1355         print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1356         print_prob(fplog, "", re->nrepl, re->prob);
1357
1358         fprintf(fplog, "Repl  number of exchanges:\n");
1359         print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1360         print_count(fplog, "", re->nrepl, re->nexchange);
1361
1362         fprintf(fplog, "Repl  average number of exchanges:\n");
1363         for (i = 1; i < re->nrepl; i++)
1364         {
1365             if (re->nattempt[i%2] == 0)
1366             {
1367                 re->prob[i] = 0;
1368             }
1369             else
1370             {
1371                 re->prob[i] =  (static_cast<real>(re->nexchange[i]))/re->nattempt[i%2];
1372             }
1373         }
1374         print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1375         print_prob(fplog, "", re->nrepl, re->prob);
1376
1377         fprintf(fplog, "\n");
1378     }
1379     /* print the transition matrix */
1380     print_transition_matrix(fplog, re->nrepl, re->nmoves, re->nattempt);
1381 }
1382
1383 //! \endcond