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