src/programs/mdrun/repl_ex.c

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