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