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