src/programs/mdrun/repl_ex.cpp

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