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