src/programs/mdrun/repl_ex.c

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