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