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