src/programs/mdrun/repl_ex.cpp

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