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