src/programs/mdrun/repl_ex.c

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