src/programs/mdrun/repl_ex.c

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