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 typedef struct gmx_repl_ex {
  55   int  repl;
  56   int  nrepl;
  57   real temp;
  58   int  type;
  59   real *q;
  60   gmx_bool bNPT;
  61   real *pres;
  62   int  *ind;
  63   int  nst;
  64   int  seed;
  65   int  nattempt[2];
  66   real *prob_sum;
  67   int  *nexchange;
  68 } t_gmx_repl_ex;
  69
  70 enum { ereTEMP, ereLAMBDA, ereNR };
  71 const char *erename[ereNR] = { "temperature", "lambda" };
  72
  73 static void repl_quantity(FILE *fplog,const gmx_multisim_t *ms,
  74                           struct gmx_repl_ex *re,int ere,real q)
  75 {
  76   real *qall;
  77   gmx_bool bDiff;
  78   int  s;
  79
  80   snew(qall,ms->nsim);
  81   qall[re->repl] = q;
  82   gmx_sum_sim(ms->nsim,qall,ms);
  83
  84   bDiff = FALSE;
  85   for(s=1; s<ms->nsim; s++)
  86     if (qall[s] != qall[0])
  87       bDiff = TRUE;
  88
  89   if (bDiff) {
  90     if (re->type >= 0 && re->type < ereNR) {
  91       gmx_fatal(FARGS,"For replica exchange both %s and %s differ",
  92                 erename[re->type],erename[ere]);
  93     } else {
  94       /* Set the replica exchange type and quantities */
  95       re->type = ere;
  96       snew(re->q,re->nrepl);
  97       for(s=0; s<ms->nsim; s++)
  98         re->q[s] = qall[s];
  99     }
 100   }
 101
 102   sfree(qall);
 103 }
 104
 105 gmx_repl_ex_t init_replica_exchange(FILE *fplog,
 106                                     const gmx_multisim_t *ms,
 107                                     const t_state *state,
 108                                     const t_inputrec *ir,
 109                                     int nst,int init_seed)
 110 {
 111   real temp,pres;
 112   int  i,j,k;
 113   struct gmx_repl_ex *re;
 114
 115   fprintf(fplog,"\nInitializing Replica Exchange\n");
 116
 117   if (ms == NULL || ms->nsim == 1)
 118     gmx_fatal(FARGS,"Nothing to exchange with only one replica, maybe you forgot to set the -multi option of mdrun?");
 119
 120   snew(re,1);
 121
 122   re->repl     = ms->sim;
 123   re->nrepl    = ms->nsim;
 124
 125   fprintf(fplog,"Repl  There are %d replicas:\n",re->nrepl);
 126
 127   check_multi_int(fplog,ms,state->natoms,"the number of atoms");
 128   check_multi_int(fplog,ms,ir->eI,"the integrator");
 129   check_multi_int(fplog,ms,ir->init_step+ir->nsteps,"init_step+nsteps");
 130   check_multi_int(fplog,ms,(ir->init_step+nst-1)/nst,
 131                   "first exchange step: init_step/-replex");
 132   check_multi_int(fplog,ms,ir->etc,"the temperature coupling");
 133   check_multi_int(fplog,ms,ir->opts.ngtc,
 134                   "the number of temperature coupling groups");
 135   check_multi_int(fplog,ms,ir->epc,"the pressure coupling");
 136   check_multi_int(fplog,ms,ir->efep,"free energy");
 137
 138   re->temp = ir->opts.ref_t[0];
 139   for(i=1; (i<ir->opts.ngtc); i++) {
 140     if (ir->opts.ref_t[i] != re->temp) {
 141       fprintf(fplog,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 142       fprintf(stderr,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
 143     }
 144   }
 145
 146     re->type = -1;
 147     for(i=0; i<ereNR; i++)
 148     {
 149         switch (i)
 150         {
 151         case ereTEMP:
 152             repl_quantity(fplog,ms,re,i,re->temp);
 153             break;
 154         case ereLAMBDA:
 155             if (ir->efep != efepNO)
 156             {
 157                 repl_quantity(fplog,ms,re,i,ir->init_lambda);
 158             }
 159             break;
 160         default:
 161             gmx_incons("Unknown replica exchange quantity");
 162         }
 163     }
 164     if (re->type == -1)
 165     {
 166         gmx_fatal(FARGS,"The properties of the %d systems are all the same, there is nothing to exchange",re->nrepl);
 167     }
 168
 169     switch (re->type)
 170     {
 171     case ereTEMP:
 172         please_cite(fplog,"Hukushima96a");
 173         if (ir->epc != epcNO)
 174         {
 175             re->bNPT = TRUE;
 176             fprintf(fplog,"Repl  Using Constant Pressure REMD.\n");
 177             please_cite(fplog,"Okabe2001a");
 178         }
 179         if (ir->etc == etcBERENDSEN)
 180         {
 181             gmx_fatal(FARGS,"REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
 182                       ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
 183         }
 184         break;
 185     case ereLAMBDA:
 186         if (ir->delta_lambda != 0)
 187         {
 188             gmx_fatal(FARGS,"delta_lambda is not zero");
 189         }
 190         break;
 191     }
 192
 193   if (re->bNPT) {
 194     snew(re->pres,re->nrepl);
 195     if (ir->epct == epctSURFACETENSION) {
 196       pres = ir->ref_p[ZZ][ZZ];
 197     } else {
 198       pres = 0;
 199       j = 0;
 200       for(i=0; i<DIM; i++)
 201         if (ir->compress[i][i] != 0) {
 202           pres += ir->ref_p[i][i];
 203           j++;
 204         }
 205       pres /= j;
 206     }
 207     re->pres[re->repl] = pres;
 208     gmx_sum_sim(re->nrepl,re->pres,ms);
 209   }
 210
 211   snew(re->ind,re->nrepl);
 212   /* Make an index for increasing temperature order */
 213   for(i=0; i<re->nrepl; i++)
 214     re->ind[i] = i;
 215   for(i=0; i<re->nrepl; i++) {
 216     for(j=i+1; j<re->nrepl; j++) {
 217       if (re->q[re->ind[j]] < re->q[re->ind[i]]) {
 218         k = re->ind[i];
 219         re->ind[i] = re->ind[j];
 220         re->ind[j] = k;
 221       } else if (re->q[re->ind[j]] == re->q[re->ind[i]]) {
 222         gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]);
 223       }
 224     }
 225   }
 226   fprintf(fplog,"Repl   ");
 227   for(i=0; i<re->nrepl; i++)
 228     fprintf(fplog," %3d  ",re->ind[i]);
 229   switch (re->type) {
 230   case ereTEMP:
 231     fprintf(fplog,"\nRepl  T");
 232     for(i=0; i<re->nrepl; i++)
 233       fprintf(fplog," %5.1f",re->q[re->ind[i]]);
 234     break;
 235   case ereLAMBDA:
 236     fprintf(fplog,"\nRepl  l");
 237     for(i=0; i<re->nrepl; i++)
 238       fprintf(fplog," %5.3f",re->q[re->ind[i]]);
 239     break;
 240   default:
 241     gmx_incons("Unknown replica exchange quantity");
 242   }
 243   if (re->bNPT) {
 244     fprintf(fplog,"\nRepl  p");
 245     for(i=0; i<re->nrepl; i++)
 246     {
 247       fprintf(fplog," %5.2f",re->pres[re->ind[i]]);
 248     }
 249
 250     for(i=0; i<re->nrepl; i++)
 251     {
 252       if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
 253       {
 254         gmx_fatal(FARGS,"The reference pressure decreases with increasing temperature");
 255       }
 256     }
 257   }
 258   fprintf(fplog,"\nRepl  ");
 259
 260   re->nst = nst;
 261   if (init_seed == -1) {
 262     if (MASTERSIM(ms))
 263       re->seed = make_seed();
 264     else
 265       re->seed = 0;
 266     gmx_sumi_sim(1,&(re->seed),ms);
 267   } else {
 268     re->seed = init_seed;
 269   }
 270   fprintf(fplog,"\nRepl  exchange interval: %d\n",re->nst);
 271   fprintf(fplog,"\nRepl  random seed: %d\n",re->seed);
 272
 273   re->nattempt[0] = 0;
 274   re->nattempt[1] = 0;
 275   snew(re->prob_sum,re->nrepl);
 276   snew(re->nexchange,re->nrepl);
 277
 278   fprintf(fplog,
 279           "Repl  below: x=exchange, pr=probability\n");
 280
 281   return re;
 282 }
 283
 284 static void exchange_reals(const gmx_multisim_t *ms,int b,real *v,int n)
 285 {
 286   real *buf;
 287   int  i;
 288
 289   if (v) {
 290     snew(buf,n);
 291 #ifdef GMX_MPI
 292     /*
 293     MPI_Sendrecv(v,  n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 294                  buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 295                  ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 296     */
 297     {
 298       MPI_Request mpi_req;
 299
 300       MPI_Isend(v,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 301                 ms->mpi_comm_masters,&mpi_req);
 302       MPI_Recv(buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
 303                ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 304       MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
 305     }
 306 #endif
 307     for(i=0; i<n; i++)
 308       v[i] = buf[i];
 309     sfree(buf);
 310   }
 311 }
 312
 313 static void exchange_doubles(const gmx_multisim_t *ms,int b,double *v,int n)
 314 {
 315   double *buf;
 316   int  i;
 317
 318   if (v) {
 319     snew(buf,n);
 320 #ifdef GMX_MPI
 321     /*
 322     MPI_Sendrecv(v,  n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 323                  buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 324                  ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 325     */
 326     {
 327       MPI_Request mpi_req;
 328
 329       MPI_Isend(v,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 330                 ms->mpi_comm_masters,&mpi_req);
 331       MPI_Recv(buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
 332                ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 333       MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
 334     }
 335 #endif
 336     for(i=0; i<n; i++)
 337       v[i] = buf[i];
 338     sfree(buf);
 339   }
 340 }
 341
 342 static void exchange_rvecs(const gmx_multisim_t *ms,int b,rvec *v,int n)
 343 {
 344   rvec *buf;
 345   int  i;
 346
 347   if (v) {
 348     snew(buf,n);
 349 #ifdef GMX_MPI
 350     /*
 351     MPI_Sendrecv(v[0],  n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 352                  buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 353                  ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 354     */
 355     {
 356       MPI_Request mpi_req;
 357
 358       MPI_Isend(v[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 359                 ms->mpi_comm_masters,&mpi_req);
 360       MPI_Recv(buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
 361                ms->mpi_comm_masters,MPI_STATUS_IGNORE);
 362       MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
 363     }
 364 #endif
 365     for(i=0; i<n; i++)
 366       copy_rvec(buf[i],v[i]);
 367     sfree(buf);
 368   }
 369 }
 370
 371 static void exchange_state(const gmx_multisim_t *ms,int b,t_state *state)
 372 {
 373   /* When t_state changes, this code should be updated. */
 374   int ngtc,nnhpres;
 375   ngtc = state->ngtc * state->nhchainlength;
 376   nnhpres = state->nnhpres* state->nhchainlength;
 377   exchange_rvecs(ms,b,state->box,DIM);
 378   exchange_rvecs(ms,b,state->box_rel,DIM);
 379   exchange_rvecs(ms,b,state->boxv,DIM);
 380   exchange_reals(ms,b,&(state->veta),1);
 381   exchange_reals(ms,b,&(state->vol0),1);
 382   exchange_rvecs(ms,b,state->svir_prev,DIM);
 383   exchange_rvecs(ms,b,state->fvir_prev,DIM);
 384   exchange_rvecs(ms,b,state->pres_prev,DIM);
 385   exchange_doubles(ms,b,state->nosehoover_xi,ngtc);
 386   exchange_doubles(ms,b,state->nosehoover_vxi,ngtc);
 387   exchange_doubles(ms,b,state->nhpres_xi,nnhpres);
 388   exchange_doubles(ms,b,state->nhpres_vxi,nnhpres);
 389   exchange_doubles(ms,b,state->therm_integral,state->ngtc);
 390   exchange_rvecs(ms,b,state->x,state->natoms);
 391   exchange_rvecs(ms,b,state->v,state->natoms);
 392   exchange_rvecs(ms,b,state->sd_X,state->natoms);
 393 }
 394
 395 static void copy_rvecs(rvec *s,rvec *d,int n)
 396 {
 397     int i;
 398
 399     if (d != NULL)
 400     {
 401         for(i=0; i<n; i++)
 402         {
 403             copy_rvec(s[i],d[i]);
 404         }
 405     }
 406 }
 407
 408 static void copy_doubles(const double *s,double *d,int n)
 409 {
 410     int i;
 411
 412     if (d != NULL)
 413     {
 414         for(i=0; i<n; i++)
 415         {
 416             d[i] = s[i];
 417         }
 418     }
 419 }
 420
 421 #define scopy_rvecs(v,n)   copy_rvecs(state->v,state_local->v,n);
 422 #define scopy_doubles(v,n) copy_doubles(state->v,state_local->v,n);
 423
 424 static void copy_state_nonatomdata(t_state *state,t_state *state_local)
 425 {
 426   /* When t_state changes, this code should be updated. */
 427   int ngtc,nnhpres;
 428   ngtc = state->ngtc * state->nhchainlength;
 429   nnhpres = state->nnhpres* state->nhchainlength;
 430   scopy_rvecs(box,DIM);
 431   scopy_rvecs(box_rel,DIM);
 432   scopy_rvecs(boxv,DIM);
 433   state_local->veta = state->veta;
 434   state_local->vol0 = state->vol0;
 435   scopy_rvecs(svir_prev,DIM);
 436   scopy_rvecs(fvir_prev,DIM);
 437   scopy_rvecs(pres_prev,DIM);
 438   scopy_doubles(nosehoover_xi,ngtc);
 439   scopy_doubles(nosehoover_vxi,ngtc);
 440   scopy_doubles(nhpres_xi,nnhpres);
 441   scopy_doubles(nhpres_vxi,nnhpres);
 442   scopy_doubles(therm_integral,state->ngtc);
 443   scopy_rvecs(x,state->natoms);
 444   scopy_rvecs(v,state->natoms);
 445   scopy_rvecs(sd_X,state->natoms);
 446 }
 447
 448 static void scale_velocities(t_state *state,real fac)
 449 {
 450   int i;
 451
 452   if (state->v)
 453     for(i=0; i<state->natoms; i++)
 454       svmul(fac,state->v[i],state->v[i]);
 455 }
 456
 457 static void pd_collect_state(const t_commrec *cr,t_state *state)
 458 {
 459   int shift;
 460
 461   if (debug)
 462     fprintf(debug,"Collecting state before exchange\n");
 463 shift = cr->nnodes - cr->npmenodes - 1;
 464   move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->x,NULL,shift,NULL);
 465   if (state->v)
 466     move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->v,NULL,shift,NULL);
 467   if (state->sd_X)
 468     move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->sd_X,NULL,shift,NULL);
 469 }
 470
 471 static void print_ind(FILE *fplog,const char *leg,int n,int *ind,gmx_bool *bEx)
 472 {
 473   int i;
 474
 475   fprintf(fplog,"Repl %2s %2d",leg,ind[0]);
 476   for(i=1; i<n; i++) {
 477     fprintf(fplog," %c %2d",(bEx!=0 && bEx[i]) ? 'x' : ' ',ind[i]);
 478   }
 479   fprintf(fplog,"\n");
 480 }
 481
 482 static void print_prob(FILE *fplog,const char *leg,int n,real *prob)
 483 {
 484   int  i;
 485   char buf[8];
 486
 487   fprintf(fplog,"Repl %2s ",leg);
 488   for(i=1; i<n; i++) {
 489     if (prob[i] >= 0) {
 490       sprintf(buf,"%4.2f",prob[i]);
 491       fprintf(fplog,"  %3s",buf[0]=='1' ? "1.0" : buf+1);
 492     } else {
 493       fprintf(fplog,"     ");
 494     }
 495   }
 496   fprintf(fplog,"\n");
 497 }
 498
 499 static void print_count(FILE *fplog,const char *leg,int n,int *count)
 500 {
 501   int i;
 502
 503   fprintf(fplog,"Repl %2s ",leg);
 504   for(i=1; i<n; i++) {
 505     fprintf(fplog," %4d",count[i]);
 506   }
 507   fprintf(fplog,"\n");
 508 }
 509
 510 static int get_replica_exchange(FILE *fplog,const gmx_multisim_t *ms,
 511                                 struct gmx_repl_ex *re,real *ener,real vol,
 512                                 int step,real time)
 513 {
 514   int  m,i,a,b;
 515   real *Epot=NULL,*Vol=NULL,*dvdl=NULL,*prob;
 516   real ediff=0,delta=0,dpV=0,betaA=0,betaB=0;
 517   gmx_bool *bEx,bPrint;
 518   int  exchange;
 519
 520   fprintf(fplog,"Replica exchange at step %d time %g\n",step,time);
 521
 522   switch (re->type) {
 523   case ereTEMP:
 524     snew(Epot,re->nrepl);
 525     snew(Vol,re->nrepl);
 526     Epot[re->repl] = ener[F_EPOT];
 527     Vol[re->repl]  = vol;
 528     gmx_sum_sim(re->nrepl,Epot,ms);
 529     gmx_sum_sim(re->nrepl,Vol,ms);
 530     break;
 531   case ereLAMBDA:
 532     snew(dvdl,re->nrepl);
 533     dvdl[re->repl] = ener[F_DVDL];
 534     gmx_sum_sim(re->nrepl,dvdl,ms);
 535     break;
 536   }
 537
 538   snew(bEx,re->nrepl);
 539   snew(prob,re->nrepl);
 540
 541   exchange = -1;
 542   m = (step / re->nst) % 2;
 543   for(i=1; i<re->nrepl; i++) {
 544     a = re->ind[i-1];
 545     b = re->ind[i];
 546     bPrint = (re->repl==a || re->repl==b);
 547     if (i % 2 == m) {
 548       switch (re->type) {
 549       case ereTEMP:
 550         /* Use equations from:
 551          * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
 552          */
 553         ediff = Epot[b] - Epot[a];
 554         betaA = 1.0/(re->q[a]*BOLTZ);
 555         betaB = 1.0/(re->q[b]*BOLTZ);
 556         delta = (betaA - betaB)*ediff;
 557         break;
 558       case ereLAMBDA:
 559         /* Here we exchange based on a linear extrapolation of dV/dlambda.
 560          * We would like to have the real energies
 561          * from foreign lambda calculations.
 562          */
 563         ediff = (dvdl[a] - dvdl[b])*(re->q[b] - re->q[a]);
 564         delta = ediff/(BOLTZ*re->temp);
 565         break;
 566       default:
 567         gmx_incons("Unknown replica exchange quantity");
 568       }
 569       if (bPrint)
 570         fprintf(fplog,"Repl %d <-> %d  dE = %10.3e",a,b,delta);
 571       if (re->bNPT) {
 572         dpV = (betaA*re->pres[a]-betaB*re->pres[b])*(Vol[b]-Vol[a])/PRESFAC;
 573         if (bPrint)
 574           fprintf(fplog,"  dpV = %10.3e  d = %10.3e",dpV,delta + dpV);
 575         delta += dpV;
 576       }
 577       if (bPrint)
 578         fprintf(fplog,"\n");
 579       if (delta <= 0) {
 580         prob[i] = 1;
 581         bEx[i] = TRUE;
 582       } else {
 583         if (delta > 100)
 584           prob[i] = 0;
 585         else
 586           prob[i] = exp(-delta);
 587         bEx[i] = (rando(&(re->seed)) < prob[i]);
 588       }
 589       re->prob_sum[i] += prob[i];
 590       if (bEx[i]) {
 591         if (a == re->repl) {
 592           exchange = b;
 593         } else if (b == re->repl) {
 594           exchange = a;
 595         }
 596         re->nexchange[i]++;
 597       }
 598     } else {
 599       prob[i] = -1;
 600       bEx[i] = FALSE;
 601     }
 602   }
 603   print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
 604   print_prob(fplog,"pr",re->nrepl,prob);
 605   fprintf(fplog,"\n");
 606
 607   sfree(bEx);
 608   sfree(prob);
 609   sfree(Epot);
 610   sfree(Vol);
 611   sfree(dvdl);
 612
 613   re->nattempt[m]++;
 614
 615   return exchange;
 616 }
 617
 618 static void write_debug_x(t_state *state)
 619 {
 620   int i;
 621
 622   if (debug) {
 623     for(i=0; i<state->natoms; i+=10)
 624       fprintf(debug,"dx %5d %10.5f %10.5f %10.5f\n",i,state->x[i][XX],state->x[i][YY],state->x[i][ZZ]);
 625   }
 626 }
 627
 628 gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
 629                       t_state *state,real *ener,
 630                       t_state *state_local,
 631                       int step,real time)
 632 {
 633     gmx_multisim_t *ms;
 634     int  exchange=-1,shift;
 635     gmx_bool bExchanged=FALSE;
 636
 637     ms = cr->ms;
 638
 639     if (MASTER(cr))
 640     {
 641         exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),
 642                                         step,time);
 643         bExchanged = (exchange >= 0);
 644     }
 645
 646     if (PAR(cr))
 647     {
 648 #ifdef GMX_MPI
 649         MPI_Bcast(&bExchanged,sizeof(gmx_bool),MPI_BYTE,MASTERRANK(cr),
 650                   cr->mpi_comm_mygroup);
 651 #endif
 652     }
 653
 654     if (bExchanged)
 655     {
 656         /* Exchange the states */
 657
 658         if (PAR(cr))
 659         {
 660             /* Collect the global state on the master node */
 661             if (DOMAINDECOMP(cr))
 662             {
 663                 dd_collect_state(cr->dd,state_local,state);
 664             }
 665             else
 666             {
 667                 pd_collect_state(cr,state);
 668             }
 669         }
 670
 671         if (MASTER(cr))
 672         {
 673             /* Exchange the global states between the master nodes */
 674             if (debug)
 675             {
 676                 fprintf(debug,"Exchanging %d with %d\n",ms->sim,exchange);
 677             }
 678             exchange_state(ms,exchange,state);
 679
 680             if (re->type == ereTEMP)
 681             {
 682                 scale_velocities(state,sqrt(re->q[ms->sim]/re->q[exchange]));
 683             }
 684         }
 685
 686         /* With domain decomposition the global state is distributed later */
 687         if (!DOMAINDECOMP(cr))
 688         {
 689             /* Copy the global state to the local state data structure */
 690             copy_state_nonatomdata(state,state_local);
 691
 692             if (PAR(cr))
 693             {
 694                 bcast_state(cr,state,FALSE);
 695             }
 696         }
 697     }
 698
 699     return bExchanged;
 700 }
 701
 702 void print_replica_exchange_statistics(FILE *fplog,struct gmx_repl_ex *re)
 703 {
 704   real *prob;
 705   int  i;
 706
 707   fprintf(fplog,"\nReplica exchange statistics\n");
 708   fprintf(fplog,"Repl  %d attempts, %d odd, %d even\n",
 709           re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
 710
 711   snew(prob,re->nrepl);
 712
 713   fprintf(fplog,"Repl  average probabilities:\n");
 714   for(i=1; i<re->nrepl; i++) {
 715     if (re->nattempt[i%2] == 0)
 716       prob[i] = 0;
 717     else
 718       prob[i] =  re->prob_sum[i]/re->nattempt[i%2];
 719   }
 720   print_ind(fplog,"",re->nrepl,re->ind,NULL);
 721   print_prob(fplog,"",re->nrepl,prob);
 722
 723   fprintf(fplog,"Repl  number of exchanges:\n");
 724   print_ind(fplog,"",re->nrepl,re->ind,NULL);
 725   print_count(fplog,"",re->nrepl,re->nexchange);
 726
 727   fprintf(fplog,"Repl  average number of exchanges:\n");
 728   for(i=1; i<re->nrepl; i++) {
 729     if (re->nattempt[i%2] == 0)
 730       prob[i] = 0;
 731     else
 732       prob[i] =  ((real)re->nexchange[i])/re->nattempt[i%2];
 733   }
 734   print_ind(fplog,"",re->nrepl,re->ind,NULL);
 735   print_prob(fplog,"",re->nrepl,prob);
 736
 737   sfree(prob);
 738
 739   fprintf(fplog,"\n");
 740 }