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