src/kernel/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");
 164     check_multi_int(fplog,ms,ir->eI,"the integrator");
 165     check_multi_large_int(fplog,ms,ir->init_step+ir->nsteps,"init_step+nsteps");
 166     check_multi_large_int(fplog,ms,(ir->init_step+nst-1)/nst,
 167                           "first exchange step: init_step/-replex");
 168     check_multi_int(fplog,ms,ir->etc,"the temperature coupling");
 169     check_multi_int(fplog,ms,ir->opts.ngtc,
 170                     "the number of temperature coupling groups");
 171     check_multi_int(fplog,ms,ir->epc,"the pressure coupling");
 172     check_multi_int(fplog,ms,ir->efep,"free energy");
 173     check_multi_int(fplog,ms,ir->fepvals->n_lambda,"number of lambda states");
 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     fprintf(fplog,"Order After Exchange: ");
 731     for (i=0;i<n;i++)
 732     {
 733         fprintf(fplog,"%d ",allswaps[i]);
 734     }
 735     fprintf(fplog,"\n\n");
 736 }
 737
 738 static void print_prob(FILE *fplog,const char *leg,int n,real *prob)
 739 {
 740     int  i;
 741     char buf[8];
 742
 743     fprintf(fplog,"Repl %2s ",leg);
 744     for(i=1; i<n; i++)
 745     {
 746         if (prob[i] >= 0)
 747         {
 748             sprintf(buf,"%4.2f",prob[i]);
 749             fprintf(fplog,"  %3s",buf[0]=='1' ? "1.0" : buf+1);
 750         }
 751         else
 752         {
 753             fprintf(fplog,"     ");
 754         }
 755     }
 756     fprintf(fplog,"\n");
 757 }
 758
 759 static void print_count(FILE *fplog,const char *leg,int n,int *count)
 760 {
 761     int i;
 762
 763     fprintf(fplog,"Repl %2s ",leg);
 764     for(i=1; i<n; i++)
 765     {
 766         fprintf(fplog," %4d",count[i]);
 767     }
 768     fprintf(fplog,"\n");
 769 }
 770
 771 static real calc_delta(FILE *fplog, gmx_bool bPrint, struct gmx_repl_ex *re, int a, int b, int ap, int bp) {
 772
 773     real ediff,dpV,delta=0;
 774     real *Epot = re->Epot;
 775     real *Vol = re->Vol;
 776     real **de = re->de;
 777     real *beta = re->beta;
 778
 779     /* Two cases; we are permuted and not.  In all cases, setting ap = a and bp = b will reduce
 780        to the non permuted case */
 781
 782     switch (re->type)
 783     {
 784     case ereTEMP:
 785         /*
 786          * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
 787          */
 788         ediff = Epot[b] - Epot[a];
 789         delta = -(beta[bp] - beta[ap])*ediff;
 790         break;
 791     case ereLAMBDA:
 792         /* two cases:  when we are permuted, and not.  */
 793         /* non-permuted:
 794            ediff =  E_new - E_old
 795                  =  [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
 796                  =  [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
 797                  =  de[b][a] + de[a][b] */
 798         /* permuted:
 799            ediff =  E_new - E_old
 800                  =  [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
 801                  =  [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
 802                  =  [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)]
 803                  =  [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)]
 804                  =  (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b])    */
 805         ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]);
 806         delta = ediff*beta[a]; /* assume all same temperature in this case */
 807         break;
 808     case ereTL:
 809         /* not permuted:  */
 810         /* delta =  reduced E_new - reduced E_old
 811                  =  [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
 812                  =  [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
 813                  =  [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
 814                     [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
 815                  =  [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
 816                     beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
 817                  =  beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
 818         /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
 819         /* permuted (big breath!) */
 820         /*   delta =  reduced E_new - reduced E_old
 821                  =  [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
 822                  =  [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
 823                  =  [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
 824                     - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
 825                     - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
 826                  =  [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
 827                     [(beta_ap H_ap(x_b)  - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
 828                     + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
 829                  =  [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
 830                     [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
 831                     + beta_pb (H_a(x_a) - H_b(x_b))  - beta_ap (H_a(x_a) - H_b(x_b))
 832                  =  ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b]  - beta_bp de[bp][b])
 833                  + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b))  */
 834         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]);
 835         break;
 836     default:
 837         gmx_incons("Unknown replica exchange quantity");
 838     }
 839     if (bPrint)
 840     {
 841         fprintf(fplog,"Repl %d <-> %d  dE_term = %10.3e (kT)\n",a,b,delta);
 842     }
 843     if (re->bNPT)
 844     {
 845         /* revist the calculation for 5.0.  Might be some improvements. */
 846         dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
 847         if (bPrint)
 848         {
 849             fprintf(fplog,"  dpV = %10.3e  d = %10.3e\nb",dpV,delta + dpV);
 850         }
 851         delta += dpV;
 852     }
 853     return delta;
 854 }
 855
 856 static void
 857 test_for_replica_exchange(FILE *fplog,
 858                           const gmx_multisim_t *ms,
 859                           struct gmx_repl_ex *re,
 860                           gmx_enerdata_t *enerd,
 861                           real vol,
 862                           gmx_large_int_t step,
 863                           real time)
 864 {
 865     int  m,i,j,a,b,ap,bp,i0,i1,tmp;
 866     real ediff=0,delta=0,dpV=0;
 867     gmx_bool bPrint,bMultiEx;
 868     gmx_bool *bEx = re->bEx;
 869     real *prob = re->prob;
 870     int *pind = re->destinations;
 871     gmx_bool bEpot=FALSE;
 872     gmx_bool bDLambda=FALSE;
 873     gmx_bool bVol=FALSE;
 874
 875     bMultiEx = (re->nex > 1);  /* multiple exchanges at each state */
 876     fprintf(fplog,"Replica exchange at step " gmx_large_int_pfmt " time %g\n",step,time);
 877
 878     if (re->bNPT)
 879     {
 880         for (i=0;i<re->nrepl;i++)
 881         {
 882             re->Vol[i] = 0;
 883         }
 884         bVol = TRUE;
 885         re->Vol[re->repl]  = vol;
 886     }
 887     if ((re->type == ereTEMP || re->type == ereTL))
 888     {
 889         for (i=0;i<re->nrepl;i++)
 890         {
 891             re->Epot[i] = 0;
 892         }
 893         bEpot = TRUE;
 894         re->Epot[re->repl] = enerd->term[F_EPOT];
 895         /* temperatures of different states*/
 896         for (i=0;i<re->nrepl;i++)
 897         {
 898             re->beta[i] = 1.0/(re->q[ereTEMP][i]*BOLTZ);
 899         }
 900     }
 901     else
 902     {
 903         for (i=0;i<re->nrepl;i++)
 904         {
 905             re->beta[i] = 1.0/(re->temp*BOLTZ);  /* we have a single temperature */
 906         }
 907     }
 908     if (re->type == ereLAMBDA || re->type == ereTL)
 909     {
 910         bDLambda = TRUE;
 911         /* lambda differences. */
 912         /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
 913            minus the energy of the jth simulation in the jth Hamiltonian */
 914         for (i=0;i<re->nrepl;i++)
 915         {
 916             for (j=0;j<re->nrepl;j++)
 917             {
 918                 re->de[i][j] = 0;
 919             }
 920         }
 921         for (i=0;i<re->nrepl;i++)
 922         {
 923             re->de[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);
 924         }
 925     }
 926
 927     /* now actually do the communication */
 928     if (bVol)
 929     {
 930         gmx_sum_sim(re->nrepl,re->Vol,ms);
 931     }
 932     if (bEpot)
 933     {
 934         gmx_sum_sim(re->nrepl,re->Epot,ms);
 935     }
 936     if (bDLambda)
 937     {
 938         for (i=0;i<re->nrepl;i++)
 939         {
 940             gmx_sum_sim(re->nrepl,re->de[i],ms);
 941         }
 942     }
 943
 944     /* make a duplicate set of indices for shuffling */
 945     for(i=0;i<re->nrepl;i++)
 946     {
 947         pind[i] = re->ind[i];
 948     }
 949
 950     if (bMultiEx)
 951     {
 952         /* multiple random switch exchange */
 953         for (i=0;i<re->nex;i++)
 954         {
 955             /* randomly select a pair  */
 956             /* find out which state it is from, and what label that state currently has */
 957             i0 = (int)(re->nrepl*rando(&(re->seed)));
 958             i1 = (int)(re->nrepl*rando(&(re->seed)));
 959             if (i0==i1)
 960             {
 961                 i--;
 962                 continue;  /* got the same pair, back up and do it again */
 963             }
 964
 965             a = re->ind[i0];
 966             b = re->ind[i1];
 967             ap = pind[i0];
 968             bp = pind[i1];
 969
 970             bPrint = FALSE; /* too noisy */
 971             delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); /* calculate the energy difference */
 972
 973             /* we actually only use the first space, since there are actually many switches between pairs. */
 974
 975             if (delta <= 0)
 976             {
 977                 /* accepted */
 978                 prob[0] = 1;
 979                 bEx[0] = TRUE;
 980             }
 981             else
 982             {
 983                 if (delta > PROBABILITYCUTOFF)
 984                 {
 985                     prob[0] = 0;
 986                 }
 987                 else
 988                 {
 989                     prob[0] = exp(-delta);
 990                 }
 991                 /* roll a number to determine if accepted */
 992                 bEx[0] = (rando(&(re->seed)) < prob[0]);
 993             }
 994             re->prob_sum[0] += prob[0];
 995
 996             if (bEx[0])
 997             {
 998                 /* swap the states */
 999                 tmp = pind[i0];
1000                 pind[i0] = pind[i1];
1001                 pind[i1] = tmp;
1002             }
1003         }
1004         re->nattempt[0]++;  /* keep track of total permutation trials here */
1005         print_allswitchind(fplog,re->nrepl,re->ind,pind,re->allswaps,re->tmpswap);
1006     }
1007     else
1008     {
1009         /* standard nearest neighbor replica exchange */
1010         m = (step / re->nst) % 2;
1011         for(i=1; i<re->nrepl; i++)
1012         {
1013             a = re->ind[i-1];
1014             b = re->ind[i];
1015
1016             bPrint = (re->repl==a || re->repl==b);
1017             if (i % 2 == m)
1018             {
1019                 delta = calc_delta(fplog,bPrint,re,a,b,a,b);
1020                 if (delta <= 0) {
1021                     /* accepted */
1022                     prob[i] = 1;
1023                     bEx[i] = TRUE;
1024                 }
1025                 else
1026                 {
1027                     if (delta > PROBABILITYCUTOFF)
1028                     {
1029                         prob[i] = 0;
1030                     }
1031                     else
1032                     {
1033                         prob[i] = exp(-delta);
1034                     }
1035                     /* roll a number to determine if accepted */
1036                     bEx[i] = (rando(&(re->seed)) < prob[i]);
1037                 }
1038                 re->prob_sum[i] += prob[i];
1039
1040                 if (bEx[i])
1041                 {
1042                     /* swap these two */
1043                     tmp = pind[i-1];
1044                     pind[i-1] = pind[i];
1045                     pind[i] = tmp;
1046                     re->nexchange[i]++;  /* statistics for back compatibility */
1047                 }
1048             }
1049             else
1050             {
1051                 prob[i] = -1;
1052                 bEx[i] = FALSE;
1053             }
1054         }
1055         /* print some statistics */
1056         print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
1057         print_prob(fplog,"pr",re->nrepl,prob);
1058         fprintf(fplog,"\n");
1059         re->nattempt[m]++;
1060     }
1061
1062     /* record which moves were made and accepted */
1063     for (i=0;i<re->nrepl;i++)
1064     {
1065         re->nmoves[re->ind[i]][pind[i]] +=1;
1066         re->nmoves[pind[i]][re->ind[i]] +=1;
1067     }
1068 }
1069
1070 static void write_debug_x(t_state *state)
1071 {
1072     int i;
1073
1074     if (debug)
1075     {
1076         for(i=0; i<state->natoms; i+=10)
1077         {
1078             fprintf(debug,"dx %5d %10.5f %10.5f %10.5f\n",i,state->x[i][XX],state->x[i][YY],state->x[i][ZZ]);
1079         }
1080     }
1081 }
1082
1083 static void
1084 cyclic_decomposition(FILE *fplog,
1085                      const int *destinations,
1086                      int **cyclic,
1087                      gmx_bool *incycle,
1088                      const int nrepl,
1089                      int *nswap)
1090 {
1091
1092     int i,j,c,p;
1093     int maxlen = 1;
1094     for (i=0;i<nrepl;i++)
1095     {
1096         incycle[i] = FALSE;
1097     }
1098     for (i=0;i<nrepl;i++)  /* one cycle for each replica */
1099     {
1100         if (incycle[i])
1101         {
1102             cyclic[i][0] = -1;
1103             continue;
1104         }
1105         cyclic[i][0] = i;
1106         incycle[i] = TRUE;
1107         c = 1;
1108         p = i;
1109         for (j=0;j<nrepl;j++)  /* potentially all cycles are part, but we will break first */
1110         {
1111             p = destinations[p]; /* start permuting */
1112             if (p==i)
1113             {
1114                 cyclic[i][c] = -1;
1115                 if (c > maxlen)
1116                 {
1117                     maxlen = c;
1118                 }
1119                 break; /* we've reached the original element, the cycle is complete, and we marked the end. */
1120             }
1121             else
1122             {
1123                 cyclic[i][c] = p;  /* each permutation gives a new member of the cycle */
1124                 incycle[p] = TRUE;
1125                 c++;
1126             }
1127         }
1128     }
1129     *nswap = maxlen - 1;
1130
1131     if (debug)
1132     {
1133         for (i=0;i<nrepl;i++)
1134         {
1135             fprintf(debug,"Cycle %d:",i);
1136             for (j=0;j<nrepl;j++)
1137             {
1138                 if (cyclic[i][j] < 0)
1139                 {
1140                     break;
1141                 }
1142                 fprintf(debug,"%2d",cyclic[i][j]);
1143             }
1144             fprintf(debug,"\n");
1145         }
1146         fflush(debug);
1147     }
1148 }
1149
1150 static void
1151 compute_exchange_order(FILE *fplog,
1152                        int **cyclic,
1153                        int **order,
1154                        const int nrepl,
1155                        const int maxswap)
1156 {
1157     int i,j;
1158
1159     for (j=0;j<maxswap;j++)
1160     {
1161         for (i=0;i<nrepl;i++)
1162         {
1163             if (cyclic[i][j+1] >= 0)
1164             {
1165                 order[cyclic[i][j+1]][j] = cyclic[i][j];
1166                 order[cyclic[i][j]][j] = cyclic[i][j+1];
1167             }
1168         }
1169         for (i=0;i<nrepl;i++)
1170         {
1171             if (order[i][j] < 0)
1172             {
1173                 order[i][j] = i; /* if it's not exchanging, it should stay this round*/
1174             }
1175         }
1176     }
1177
1178     if (debug)
1179     {
1180         fprintf(fplog,"Replica Exchange Order\n");
1181         for (i=0;i<nrepl;i++)
1182         {
1183             fprintf(fplog,"Replica %d:",i);
1184             for (j=0;j<maxswap;j++)
1185             {
1186                 if (order[i][j] < 0) break;
1187                 fprintf(debug,"%2d",order[i][j]);
1188             }
1189             fprintf(fplog,"\n");
1190         }
1191         fflush(fplog);
1192     }
1193 }
1194
1195 static void
1196 prepare_to_do_exchange(FILE *fplog,
1197                        const int *destinations,
1198                        const int replica_id,
1199                        const int nrepl,
1200                        int *maxswap,
1201                        int **order,
1202                        int **cyclic,
1203                        int *incycle,
1204                        gmx_bool *bThisReplicaExchanged)
1205 {
1206     int i,j;
1207     /* Hold the cyclic decomposition of the (multiple) replica
1208      * exchange. */
1209     gmx_bool bAnyReplicaExchanged = FALSE;
1210     *bThisReplicaExchanged = FALSE;
1211
1212     for (i = 0; i < nrepl; i++)
1213     {
1214         if (destinations[i] != i) {
1215             /* only mark as exchanged if the index has been shuffled */
1216             bAnyReplicaExchanged = TRUE;
1217             break;
1218         }
1219     }
1220     if (bAnyReplicaExchanged)
1221     {
1222         /* reinitialize the placeholder arrays */
1223         for (i = 0; i < nrepl; i++)
1224         {
1225             for (j = 0; j < nrepl; j++)
1226             {
1227                 cyclic[i][j] = -1;
1228                 order[i][j] = -1;
1229             }
1230         }
1231
1232         /* Identify the cyclic decomposition of the permutation (very
1233          * fast if neighbor replica exchange). */
1234         cyclic_decomposition(fplog,destinations,cyclic,incycle,nrepl,maxswap);
1235
1236         /* Now translate the decomposition into a replica exchange
1237          * order at each step. */
1238         compute_exchange_order(fplog,cyclic,order,nrepl,*maxswap);
1239
1240         /* Did this replica do any exchange at any point? */
1241         for (j = 0; j < *maxswap; j++)
1242         {
1243             if (replica_id != order[replica_id][j])
1244             {
1245                 *bThisReplicaExchanged = TRUE;
1246                 break;
1247             }
1248         }
1249     }
1250 }
1251
1252 gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
1253                           t_state *state,gmx_enerdata_t *enerd,
1254                           t_state *state_local,gmx_large_int_t step,real time)
1255 {
1256     int i,j;
1257     int replica_id = 0;
1258     int exchange_partner;
1259     int maxswap = 0;
1260     /* Number of rounds of exchanges needed to deal with any multiple
1261      * exchanges. */
1262     /* Where each replica ends up after the exchange attempt(s). */
1263     /* The order in which multiple exchanges will occur. */
1264     gmx_bool bThisReplicaExchanged = FALSE;
1265
1266     if (MASTER(cr))
1267     {
1268         replica_id  = re->repl;
1269         test_for_replica_exchange(fplog,cr->ms,re,enerd,det(state_local->box),step,time);
1270         prepare_to_do_exchange(fplog,re->destinations,replica_id,re->nrepl,&maxswap,
1271                                re->order,re->cyclic,re->incycle,&bThisReplicaExchanged);
1272     }
1273     /* Do intra-simulation broadcast so all processors belonging to
1274      * each simulation know whether they need to participate in
1275      * collecting the state. Otherwise, they might as well get on with
1276      * the next thing to do. */
1277     if (PAR(cr))
1278     {
1279 #ifdef GMX_MPI
1280         MPI_Bcast(&bThisReplicaExchanged,sizeof(gmx_bool),MPI_BYTE,MASTERRANK(cr),
1281                   cr->mpi_comm_mygroup);
1282 #endif
1283     }
1284
1285     if (bThisReplicaExchanged)
1286     {
1287         /* Exchange the states */
1288
1289         if (PAR(cr))
1290         {
1291             /* Collect the global state on the master node */
1292             if (DOMAINDECOMP(cr))
1293             {
1294                 dd_collect_state(cr->dd,state_local,state);
1295             }
1296             else
1297             {
1298                 pd_collect_state(cr,state);
1299             }
1300         }
1301
1302         if (MASTER(cr))
1303         {
1304             /* There will be only one swap cycle with standard replica
1305              * exchange, but there may be multiple swap cycles if we
1306              * allow multiple swaps. */
1307             for (j = 0; j < maxswap; j++)
1308             {
1309                 exchange_partner = re->order[replica_id][j];
1310
1311                 if (exchange_partner != replica_id)
1312                 {
1313                     /* Exchange the global states between the master nodes */
1314                     if (debug)
1315                     {
1316                         fprintf(debug,"Exchanging %d with %d\n",replica_id,exchange_partner);
1317                     }
1318                     exchange_state(cr->ms,exchange_partner,state);
1319                 }
1320             }
1321             /* For temperature-type replica exchange, we need to scale
1322              * the velocities. */
1323             if (re->type == ereTEMP || re->type == ereTL)
1324             {
1325                 scale_velocities(state,sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));
1326             }
1327
1328         }
1329
1330         /* With domain decomposition the global state is distributed later */
1331         if (!DOMAINDECOMP(cr))
1332         {
1333             /* Copy the global state to the local state data structure */
1334             copy_state_nonatomdata(state,state_local);
1335
1336             if (PAR(cr))
1337             {
1338                 bcast_state(cr,state,FALSE);
1339             }
1340         }
1341     }
1342
1343     return bThisReplicaExchanged;
1344 }
1345
1346 void print_replica_exchange_statistics(FILE *fplog,struct gmx_repl_ex *re)
1347 {
1348     int  i;
1349
1350     fprintf(fplog,"\nReplica exchange statistics\n");
1351
1352     if (re->nex == 0)
1353     {
1354         fprintf(fplog,"Repl  %d attempts, %d odd, %d even\n",
1355                 re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
1356
1357         fprintf(fplog,"Repl  average probabilities:\n");
1358         for(i=1; i<re->nrepl; i++)
1359         {
1360             if (re->nattempt[i%2] == 0)
1361             {
1362                 re->prob[i] = 0;
1363             }
1364             else
1365             {
1366                 re->prob[i] =  re->prob_sum[i]/re->nattempt[i%2];
1367             }
1368         }
1369         print_ind(fplog,"",re->nrepl,re->ind,NULL);
1370         print_prob(fplog,"",re->nrepl,re->prob);
1371
1372         fprintf(fplog,"Repl  number of exchanges:\n");
1373         print_ind(fplog,"",re->nrepl,re->ind,NULL);
1374         print_count(fplog,"",re->nrepl,re->nexchange);
1375
1376         fprintf(fplog,"Repl  average number of exchanges:\n");
1377         for(i=1; i<re->nrepl; i++)
1378         {
1379             if (re->nattempt[i%2] == 0)
1380             {
1381                 re->prob[i] = 0;
1382             }
1383             else
1384             {
1385                 re->prob[i] =  ((real)re->nexchange[i])/re->nattempt[i%2];
1386             }
1387         }
1388         print_ind(fplog,"",re->nrepl,re->ind,NULL);
1389         print_prob(fplog,"",re->nrepl,re->prob);
1390
1391         fprintf(fplog,"\n");
1392     }
1393     /* print the transition matrix */
1394     print_transition_matrix(fplog,"",re->nrepl,re->nmoves,re->nattempt);
1395 }