src/programs/mdrun/repl_ex.cpp

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