2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
56 #define PROBABILITYCUTOFF 100
57 /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
59 enum { ereTEMP, ereLAMBDA, ereENDSINGLE ,ereTL, ereNR };
60 const char *erename[ereNR] = { "temperature", "lambda", "end_single_marker", "temperature and lambda"};
61 /* end_single_marker merely notes the end of single variable replica exchange. All types higher than
62 it are multiple replica exchange methods */
63 /* Eventually, should add 'pressure', 'temperature and pressure', 'lambda_and_pressure', 'temperature_lambda_pressure'?;
64 Let's wait until we feel better about the pressure control methods giving exact ensembles. Right now, we assume constant pressure */
66 typedef struct gmx_repl_ex
85 /* these are helper arrays for replica exchange; allocated here so they
86 don't have to be allocated each time */
94 /* helper arrays to hold the quantities that are exchanged */
103 static gmx_bool repl_quantity(FILE *fplog,const gmx_multisim_t *ms,
104 struct gmx_repl_ex *re,int ere,real q)
112 gmx_sum_sim(ms->nsim,qall,ms);
115 for (s=1; s<ms->nsim; s++)
117 if (qall[s] != qall[0])
125 /* Set the replica exchange type and quantities */
128 snew(re->q[ere],re->nrepl);
129 for(s=0; s<ms->nsim; s++)
131 re->q[ere][s] = qall[s];
138 gmx_repl_ex_t init_replica_exchange(FILE *fplog,
139 const gmx_multisim_t *ms,
140 const t_state *state,
141 const t_inputrec *ir,
142 int nst, int nex, int init_seed)
146 struct gmx_repl_ex *re;
148 gmx_bool bLambda=FALSE;
150 fprintf(fplog,"\nInitializing Replica Exchange\n");
152 if (ms == NULL || ms->nsim == 1)
154 gmx_fatal(FARGS,"Nothing to exchange with only one replica, maybe you forgot to set the -multi option of mdrun?");
160 re->nrepl = ms->nsim;
161 snew(re->q,ereENDSINGLE);
163 fprintf(fplog,"Repl There are %d replicas:\n",re->nrepl);
165 check_multi_int(fplog,ms,state->natoms,"the number of atoms",FALSE);
166 check_multi_int(fplog,ms,ir->eI,"the integrator",FALSE);
167 check_multi_large_int(fplog,ms,ir->init_step+ir->nsteps,"init_step+nsteps",FALSE);
168 check_multi_large_int(fplog,ms,(ir->init_step+nst-1)/nst,
169 "first exchange step: init_step/-replex",FALSE);
170 check_multi_int(fplog,ms,ir->etc,"the temperature coupling",FALSE);
171 check_multi_int(fplog,ms,ir->opts.ngtc,
172 "the number of temperature coupling groups",FALSE);
173 check_multi_int(fplog,ms,ir->epc,"the pressure coupling",FALSE);
174 check_multi_int(fplog,ms,ir->efep,"free energy",FALSE);
175 check_multi_int(fplog,ms,ir->fepvals->n_lambda,"number of lambda states",FALSE);
177 re->temp = ir->opts.ref_t[0];
178 for(i=1; (i<ir->opts.ngtc); i++)
180 if (ir->opts.ref_t[i] != re->temp)
182 fprintf(fplog,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
183 fprintf(stderr,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
188 bTemp = repl_quantity(fplog,ms,re,ereTEMP,re->temp);
189 if (ir->efep != efepNO)
191 bLambda = repl_quantity(fplog,ms,re,ereLAMBDA,(real)ir->fepvals->init_fep_state);
193 if (re->type == -1) /* nothing was assigned */
195 gmx_fatal(FARGS,"The properties of the %d systems are all the same, there is nothing to exchange",re->nrepl);
197 if (bLambda && bTemp) {
203 please_cite(fplog,"Sugita1999a");
204 if (ir->epc != epcNO)
207 fprintf(fplog,"Repl Using Constant Pressure REMD.\n");
208 please_cite(fplog,"Okabe2001a");
210 if (ir->etc == etcBERENDSEN)
212 gmx_fatal(FARGS,"REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
213 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
217 if (ir->fepvals->delta_lambda != 0) /* check this? */
219 gmx_fatal(FARGS,"delta_lambda is not zero");
224 snew(re->pres,re->nrepl);
225 if (ir->epct == epctSURFACETENSION)
227 pres = ir->ref_p[ZZ][ZZ];
235 if (ir->compress[i][i] != 0)
237 pres += ir->ref_p[i][i];
243 re->pres[re->repl] = pres;
244 gmx_sum_sim(re->nrepl,re->pres,ms);
247 /* Make an index for increasing replica order */
248 /* only makes sense if one or the other is varying, not both!
249 if both are varying, we trust the order the person gave. */
250 snew(re->ind,re->nrepl);
251 for(i=0; i<re->nrepl; i++)
256 if (re->type<ereENDSINGLE) {
258 for(i=0; i<re->nrepl; i++)
260 for(j=i+1; j<re->nrepl; j++)
262 if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]])
265 re->ind[i] = re->ind[j];
268 else if (re->q[re->type][re->ind[j]] == re->q[re->type][re->ind[i]])
270 gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]);
276 /* keep track of all the swaps, starting with the initial placement. */
277 snew(re->allswaps,re->nrepl);
278 for(i=0; i<re->nrepl; i++)
280 re->allswaps[i] = re->ind[i];
286 fprintf(fplog,"\nReplica exchange in temperature\n");
287 for(i=0; i<re->nrepl; i++)
289 fprintf(fplog," %5.1f",re->q[re->type][re->ind[i]]);
294 fprintf(fplog,"\nReplica exchange in lambda\n");
295 for(i=0; i<re->nrepl; i++)
297 fprintf(fplog," %3d",(int)re->q[re->type][re->ind[i]]);
302 fprintf(fplog,"\nReplica exchange in temperature and lambda state\n");
303 for(i=0; i<re->nrepl; i++)
305 fprintf(fplog," %5.1f",re->q[ereTEMP][re->ind[i]]);
308 for(i=0; i<re->nrepl; i++)
310 fprintf(fplog," %5d",(int)re->q[ereLAMBDA][re->ind[i]]);
315 gmx_incons("Unknown replica exchange quantity");
319 fprintf(fplog,"\nRepl p");
320 for(i=0; i<re->nrepl; i++)
322 fprintf(fplog," %5.2f",re->pres[re->ind[i]]);
325 for(i=0; i<re->nrepl; i++)
327 if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
329 fprintf(fplog,"\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
330 fprintf(stderr,"\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
339 re->seed = make_seed();
345 gmx_sumi_sim(1,&(re->seed),ms);
349 re->seed = init_seed;
351 fprintf(fplog,"\nReplica exchange interval: %d\n",re->nst);
352 fprintf(fplog,"\nReplica random seed: %d\n",re->seed);
357 snew(re->prob_sum,re->nrepl);
358 snew(re->nexchange,re->nrepl);
359 snew(re->nmoves,re->nrepl);
360 for (i=0;i<re->nrepl;i++)
362 snew(re->nmoves[i],re->nrepl);
364 fprintf(fplog,"Replica exchange information below: x=exchange, pr=probability\n");
366 /* generate space for the helper functions so we don't have to snew each time */
368 snew(re->destinations,re->nrepl);
369 snew(re->incycle,re->nrepl);
370 snew(re->tmpswap,re->nrepl);
371 snew(re->cyclic,re->nrepl);
372 snew(re->order,re->nrepl);
373 for (i=0;i<re->nrepl;i++)
375 snew(re->cyclic[i],re->nrepl);
376 snew(re->order[i],re->nrepl);
378 /* allocate space for the functions storing the data for the replicas */
379 /* not all of these arrays needed in all cases, but they don't take
380 up much space, since the max size is nrepl**2 */
381 snew(re->prob,re->nrepl);
382 snew(re->bEx,re->nrepl);
383 snew(re->beta,re->nrepl);
384 snew(re->Vol,re->nrepl);
385 snew(re->Epot,re->nrepl);
386 snew(re->de,re->nrepl);
387 for (i=0;i<re->nrepl;i++)
389 snew(re->de[i],re->nrepl);
395 static void exchange_reals(const gmx_multisim_t *ms,int b,real *v,int n)
405 MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
406 buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
407 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
412 MPI_Isend(v,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
413 ms->mpi_comm_masters,&mpi_req);
414 MPI_Recv(buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
415 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
416 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
428 static void exchange_ints(const gmx_multisim_t *ms,int b,int *v,int n)
437 MPI_Sendrecv(v, n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
438 buf,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
439 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
444 MPI_Isend(v,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
445 ms->mpi_comm_masters,&mpi_req);
446 MPI_Recv(buf,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
447 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
448 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
459 static void exchange_doubles(const gmx_multisim_t *ms,int b,double *v,int n)
469 MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
470 buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
471 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
476 MPI_Isend(v,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
477 ms->mpi_comm_masters,&mpi_req);
478 MPI_Recv(buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
479 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
480 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
491 static void exchange_rvecs(const gmx_multisim_t *ms,int b,rvec *v,int n)
501 MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
502 buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
503 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
508 MPI_Isend(v[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
509 ms->mpi_comm_masters,&mpi_req);
510 MPI_Recv(buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
511 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
512 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
517 copy_rvec(buf[i],v[i]);
523 static void exchange_state(const gmx_multisim_t *ms,int b,t_state *state)
525 /* When t_state changes, this code should be updated. */
527 ngtc = state->ngtc * state->nhchainlength;
528 nnhpres = state->nnhpres* state->nhchainlength;
529 exchange_rvecs(ms,b,state->box,DIM);
530 exchange_rvecs(ms,b,state->box_rel,DIM);
531 exchange_rvecs(ms,b,state->boxv,DIM);
532 exchange_reals(ms,b,&(state->veta),1);
533 exchange_reals(ms,b,&(state->vol0),1);
534 exchange_rvecs(ms,b,state->svir_prev,DIM);
535 exchange_rvecs(ms,b,state->fvir_prev,DIM);
536 exchange_rvecs(ms,b,state->pres_prev,DIM);
537 exchange_doubles(ms,b,state->nosehoover_xi,ngtc);
538 exchange_doubles(ms,b,state->nosehoover_vxi,ngtc);
539 exchange_doubles(ms,b,state->nhpres_xi,nnhpres);
540 exchange_doubles(ms,b,state->nhpres_vxi,nnhpres);
541 exchange_doubles(ms,b,state->therm_integral,state->ngtc);
542 exchange_rvecs(ms,b,state->x,state->natoms);
543 exchange_rvecs(ms,b,state->v,state->natoms);
544 exchange_rvecs(ms,b,state->sd_X,state->natoms);
547 static void copy_rvecs(rvec *s,rvec *d,int n)
555 copy_rvec(s[i],d[i]);
560 static void copy_doubles(const double *s,double *d,int n)
573 static void copy_reals(const real *s,real *d,int n)
586 static void copy_ints(const int *s,int *d,int n)
599 #define scopy_rvecs(v,n) copy_rvecs(state->v,state_local->v,n);
600 #define scopy_doubles(v,n) copy_doubles(state->v,state_local->v,n);
601 #define scopy_reals(v,n) copy_reals(state->v,state_local->v,n);
602 #define scopy_ints(v,n) copy_ints(state->v,state_local->v,n);
604 static void copy_state_nonatomdata(t_state *state,t_state *state_local)
606 /* When t_state changes, this code should be updated. */
608 ngtc = state->ngtc * state->nhchainlength;
609 nnhpres = state->nnhpres* state->nhchainlength;
610 scopy_rvecs(box,DIM);
611 scopy_rvecs(box_rel,DIM);
612 scopy_rvecs(boxv,DIM);
613 state_local->veta = state->veta;
614 state_local->vol0 = state->vol0;
615 scopy_rvecs(svir_prev,DIM);
616 scopy_rvecs(fvir_prev,DIM);
617 scopy_rvecs(pres_prev,DIM);
618 scopy_doubles(nosehoover_xi,ngtc);
619 scopy_doubles(nosehoover_vxi,ngtc);
620 scopy_doubles(nhpres_xi,nnhpres);
621 scopy_doubles(nhpres_vxi,nnhpres);
622 scopy_doubles(therm_integral,state->ngtc);
623 scopy_rvecs(x,state->natoms);
624 scopy_rvecs(v,state->natoms);
625 scopy_rvecs(sd_X,state->natoms);
626 copy_ints(&(state->fep_state),&(state_local->fep_state),1);
627 scopy_reals(lambda,efptNR);
630 static void scale_velocities(t_state *state,real fac)
636 for(i=0; i<state->natoms; i++)
638 svmul(fac,state->v[i],state->v[i]);
643 static void pd_collect_state(const t_commrec *cr,t_state *state)
649 fprintf(debug,"Collecting state before exchange\n");
651 shift = cr->nnodes - cr->npmenodes - 1;
652 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->x,NULL,shift,NULL);
655 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->v,NULL,shift,NULL);
659 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->sd_X,NULL,shift,NULL);
663 static void print_transition_matrix(FILE *fplog,const char *leg,int n,int **nmoves, int *nattempt)
668 ntot = nattempt[0] + nattempt[1];
670 fprintf(fplog,"Repl");
673 fprintf(fplog," "); /* put the title closer to the center */
675 fprintf(fplog,"Empirical Transition Matrix\n");
677 fprintf(fplog,"Repl");
680 fprintf(fplog,"%8d",(i+1));
686 fprintf(fplog,"Repl");
690 if (nmoves[i][j] > 0)
692 Tprint = nmoves[i][j]/(2.0*ntot);
694 fprintf(fplog,"%8.4f",Tprint);
696 fprintf(fplog,"%3d\n",i);
700 static void print_ind(FILE *fplog,const char *leg,int n,int *ind,gmx_bool *bEx)
704 fprintf(fplog,"Repl %2s %2d",leg,ind[0]);
707 fprintf(fplog," %c %2d",(bEx!=0 && bEx[i]) ? 'x' : ' ',ind[i]);
712 static void print_allswitchind(FILE *fplog,int n,int *ind,int *pind, int *allswaps, int *tmpswap)
718 tmpswap[i] = allswaps[i];
722 allswaps[i] = tmpswap[pind[i]];
725 fprintf(fplog,"\nAccepted Exchanges: ");
728 fprintf(fplog,"%d ",pind[i]);
732 /* the "Order After Exchange" is the state label corresponding to the configuration that
733 started in state listed in order, i.e.
738 configuration starting in simulation 3 is now in simulation 0,
739 configuration starting in simulation 0 is now in simulation 1,
740 configuration starting in simulation 1 is now in simulation 2,
741 configuration starting in simulation 2 is now in simulation 3
743 fprintf(fplog,"Order After Exchange: ");
746 fprintf(fplog,"%d ",allswaps[i]);
748 fprintf(fplog,"\n\n");
751 static void print_prob(FILE *fplog,const char *leg,int n,real *prob)
756 fprintf(fplog,"Repl %2s ",leg);
761 sprintf(buf,"%4.2f",prob[i]);
762 fprintf(fplog," %3s",buf[0]=='1' ? "1.0" : buf+1);
772 static void print_count(FILE *fplog,const char *leg,int n,int *count)
776 fprintf(fplog,"Repl %2s ",leg);
779 fprintf(fplog," %4d",count[i]);
784 static real calc_delta(FILE *fplog, gmx_bool bPrint, struct gmx_repl_ex *re, int a, int b, int ap, int bp) {
786 real ediff,dpV,delta=0;
787 real *Epot = re->Epot;
790 real *beta = re->beta;
792 /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce
793 to the non permuted case */
799 * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
801 ediff = Epot[b] - Epot[a];
802 delta = -(beta[bp] - beta[ap])*ediff;
805 /* two cases: when we are permuted, and not. */
807 ediff = E_new - E_old
808 = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
809 = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
810 = de[b][a] + de[a][b] */
813 ediff = E_new - E_old
814 = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
815 = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
816 = [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)]
817 = [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)]
818 = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]) */
819 /* but, in the current code implementation, we flip configurations, not indices . . .
820 So let's examine that.
821 = [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)]
822 = [H_b(x_ap) - H_a(x_ap)] + [H_a(x_bp) - H_b(x_pb)]
823 = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp]
824 So, if we exchange b<=> bp and a<=> ap, we return to the same result.
825 So the simple solution is to flip the
826 position of perturbed and original indices in the tests.
829 ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]);
830 delta = ediff*beta[a]; /* assume all same temperature in this case */
834 /* delta = reduced E_new - reduced E_old
835 = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
836 = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
837 = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
838 [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
839 = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
840 beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
841 = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
842 /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
843 /* permuted (big breath!) */
844 /* delta = reduced E_new - reduced E_old
845 = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
846 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
847 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
848 - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
849 - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
850 = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
851 [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
852 + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
853 = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
854 [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
855 + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b))
856 = ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b] - beta_bp de[bp][b])
857 + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */
858 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]);
861 gmx_incons("Unknown replica exchange quantity");
865 fprintf(fplog,"Repl %d <-> %d dE_term = %10.3e (kT)\n",a,b,delta);
869 /* revist the calculation for 5.0. Might be some improvements. */
870 dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
873 fprintf(fplog," dpV = %10.3e d = %10.3e\nb",dpV,delta + dpV);
881 test_for_replica_exchange(FILE *fplog,
882 const gmx_multisim_t *ms,
883 struct gmx_repl_ex *re,
884 gmx_enerdata_t *enerd,
886 gmx_large_int_t step,
889 int m,i,j,a,b,ap,bp,i0,i1,tmp;
890 real ediff=0,delta=0,dpV=0;
891 gmx_bool bPrint,bMultiEx;
892 gmx_bool *bEx = re->bEx;
893 real *prob = re->prob;
894 int *pind = re->destinations; /* permuted index */
895 gmx_bool bEpot=FALSE;
896 gmx_bool bDLambda=FALSE;
899 bMultiEx = (re->nex > 1); /* multiple exchanges at each state */
900 fprintf(fplog,"Replica exchange at step " gmx_large_int_pfmt " time %g\n",step,time);
904 for (i=0;i<re->nrepl;i++)
909 re->Vol[re->repl] = vol;
911 if ((re->type == ereTEMP || re->type == ereTL))
913 for (i=0;i<re->nrepl;i++)
918 re->Epot[re->repl] = enerd->term[F_EPOT];
919 /* temperatures of different states*/
920 for (i=0;i<re->nrepl;i++)
922 re->beta[i] = 1.0/(re->q[ereTEMP][i]*BOLTZ);
927 for (i=0;i<re->nrepl;i++)
929 re->beta[i] = 1.0/(re->temp*BOLTZ); /* we have a single temperature */
932 if (re->type == ereLAMBDA || re->type == ereTL)
935 /* lambda differences. */
936 /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
937 minus the energy of the jth simulation in the jth Hamiltonian */
938 for (i=0;i<re->nrepl;i++)
940 for (j=0;j<re->nrepl;j++)
945 for (i=0;i<re->nrepl;i++)
947 re->de[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);
951 /* now actually do the communication */
954 gmx_sum_sim(re->nrepl,re->Vol,ms);
958 gmx_sum_sim(re->nrepl,re->Epot,ms);
962 for (i=0;i<re->nrepl;i++)
964 gmx_sum_sim(re->nrepl,re->de[i],ms);
968 /* make a duplicate set of indices for shuffling */
969 for(i=0;i<re->nrepl;i++)
971 pind[i] = re->ind[i];
976 /* multiple random switch exchange */
977 for (i=0;i<re->nex;i++)
979 /* randomly select a pair */
980 /* in theory, could reduce this by identifying only which switches had a nonneglibible
981 probability of occurring (log p > -100) and only operate on those switches */
982 /* find out which state it is from, and what label that state currently has. Likely
983 more work that useful. */
984 i0 = (int)(re->nrepl*rando(&(re->seed)));
985 i1 = (int)(re->nrepl*rando(&(re->seed)));
989 continue; /* self-exchange, back up and do it again */
992 a = re->ind[i0]; /* what are the indices of these states? */
997 bPrint = FALSE; /* too noisy */
998 /* calculate the energy difference */
999 /* if the code changes to flip the STATES, rather than the configurations,
1000 use the commented version of the code */
1001 /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */
1002 delta = calc_delta(fplog,bPrint,re,ap,bp,a,b);
1004 /* we actually only use the first space in the prob and bEx array,
1005 since there are actually many switches between pairs. */
1015 if (delta > PROBABILITYCUTOFF)
1021 prob[0] = exp(-delta);
1023 /* roll a number to determine if accepted */
1024 bEx[0] = (rando(&(re->seed)) < prob[0]);
1026 re->prob_sum[0] += prob[0];
1030 /* swap the states */
1032 pind[i0] = pind[i1];
1036 re->nattempt[0]++; /* keep track of total permutation trials here */
1037 print_allswitchind(fplog,re->nrepl,re->ind,pind,re->allswaps,re->tmpswap);
1041 /* standard nearest neighbor replica exchange */
1042 m = (step / re->nst) % 2;
1043 for(i=1; i<re->nrepl; i++)
1048 bPrint = (re->repl==a || re->repl==b);
1051 delta = calc_delta(fplog,bPrint,re,a,b,a,b);
1059 if (delta > PROBABILITYCUTOFF)
1065 prob[i] = exp(-delta);
1067 /* roll a number to determine if accepted */
1068 bEx[i] = (rando(&(re->seed)) < prob[i]);
1070 re->prob_sum[i] += prob[i];
1074 /* swap these two */
1076 pind[i-1] = pind[i];
1078 re->nexchange[i]++; /* statistics for back compatibility */
1087 /* print some statistics */
1088 print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
1089 print_prob(fplog,"pr",re->nrepl,prob);
1090 fprintf(fplog,"\n");
1094 /* record which moves were made and accepted */
1095 for (i=0;i<re->nrepl;i++)
1097 re->nmoves[re->ind[i]][pind[i]] +=1;
1098 re->nmoves[pind[i]][re->ind[i]] +=1;
1100 fflush(fplog); /* make sure we can see what the last exchange was */
1103 static void write_debug_x(t_state *state)
1109 for(i=0; i<state->natoms; i+=10)
1111 fprintf(debug,"dx %5d %10.5f %10.5f %10.5f\n",i,state->x[i][XX],state->x[i][YY],state->x[i][ZZ]);
1117 cyclic_decomposition(FILE *fplog,
1118 const int *destinations,
1127 for (i=0;i<nrepl;i++)
1131 for (i=0;i<nrepl;i++) /* one cycle for each replica */
1142 for (j=0;j<nrepl;j++) /* potentially all cycles are part, but we will break first */
1144 p = destinations[p]; /* start permuting */
1152 break; /* we've reached the original element, the cycle is complete, and we marked the end. */
1156 cyclic[i][c] = p; /* each permutation gives a new member of the cycle */
1162 *nswap = maxlen - 1;
1166 for (i=0;i<nrepl;i++)
1168 fprintf(debug,"Cycle %d:",i);
1169 for (j=0;j<nrepl;j++)
1171 if (cyclic[i][j] < 0)
1175 fprintf(debug,"%2d",cyclic[i][j]);
1177 fprintf(debug,"\n");
1184 compute_exchange_order(FILE *fplog,
1192 for (j=0;j<maxswap;j++)
1194 for (i=0;i<nrepl;i++)
1196 if (cyclic[i][j+1] >= 0)
1198 order[cyclic[i][j+1]][j] = cyclic[i][j];
1199 order[cyclic[i][j]][j] = cyclic[i][j+1];
1202 for (i=0;i<nrepl;i++)
1204 if (order[i][j] < 0)
1206 order[i][j] = i; /* if it's not exchanging, it should stay this round*/
1213 fprintf(fplog,"Replica Exchange Order\n");
1214 for (i=0;i<nrepl;i++)
1216 fprintf(fplog,"Replica %d:",i);
1217 for (j=0;j<maxswap;j++)
1219 if (order[i][j] < 0) break;
1220 fprintf(debug,"%2d",order[i][j]);
1222 fprintf(fplog,"\n");
1229 prepare_to_do_exchange(FILE *fplog,
1230 const int *destinations,
1231 const int replica_id,
1237 gmx_bool *bThisReplicaExchanged)
1240 /* Hold the cyclic decomposition of the (multiple) replica
1242 gmx_bool bAnyReplicaExchanged = FALSE;
1243 *bThisReplicaExchanged = FALSE;
1245 for (i = 0; i < nrepl; i++)
1247 if (destinations[i] != i) {
1248 /* only mark as exchanged if the index has been shuffled */
1249 bAnyReplicaExchanged = TRUE;
1253 if (bAnyReplicaExchanged)
1255 /* reinitialize the placeholder arrays */
1256 for (i = 0; i < nrepl; i++)
1258 for (j = 0; j < nrepl; j++)
1265 /* Identify the cyclic decomposition of the permutation (very
1266 * fast if neighbor replica exchange). */
1267 cyclic_decomposition(fplog,destinations,cyclic,incycle,nrepl,maxswap);
1269 /* Now translate the decomposition into a replica exchange
1270 * order at each step. */
1271 compute_exchange_order(fplog,cyclic,order,nrepl,*maxswap);
1273 /* Did this replica do any exchange at any point? */
1274 for (j = 0; j < *maxswap; j++)
1276 if (replica_id != order[replica_id][j])
1278 *bThisReplicaExchanged = TRUE;
1285 gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
1286 t_state *state,gmx_enerdata_t *enerd,
1287 t_state *state_local,gmx_large_int_t step,real time)
1291 int exchange_partner;
1293 /* Number of rounds of exchanges needed to deal with any multiple
1295 /* Where each replica ends up after the exchange attempt(s). */
1296 /* The order in which multiple exchanges will occur. */
1297 gmx_bool bThisReplicaExchanged = FALSE;
1301 replica_id = re->repl;
1302 test_for_replica_exchange(fplog,cr->ms,re,enerd,det(state_local->box),step,time);
1303 prepare_to_do_exchange(fplog,re->destinations,replica_id,re->nrepl,&maxswap,
1304 re->order,re->cyclic,re->incycle,&bThisReplicaExchanged);
1306 /* Do intra-simulation broadcast so all processors belonging to
1307 * each simulation know whether they need to participate in
1308 * collecting the state. Otherwise, they might as well get on with
1309 * the next thing to do. */
1313 MPI_Bcast(&bThisReplicaExchanged,sizeof(gmx_bool),MPI_BYTE,MASTERRANK(cr),
1314 cr->mpi_comm_mygroup);
1318 if (bThisReplicaExchanged)
1320 /* Exchange the states */
1324 /* Collect the global state on the master node */
1325 if (DOMAINDECOMP(cr))
1327 dd_collect_state(cr->dd,state_local,state);
1331 pd_collect_state(cr,state);
1337 /* There will be only one swap cycle with standard replica
1338 * exchange, but there may be multiple swap cycles if we
1339 * allow multiple swaps. */
1341 for (j = 0; j < maxswap; j++)
1343 exchange_partner = re->order[replica_id][j];
1345 if (exchange_partner != replica_id)
1347 /* Exchange the global states between the master nodes */
1350 fprintf(debug,"Exchanging %d with %d\n",replica_id,exchange_partner);
1352 exchange_state(cr->ms,exchange_partner,state);
1355 /* For temperature-type replica exchange, we need to scale
1356 * the velocities. */
1357 if (re->type == ereTEMP || re->type == ereTL)
1359 scale_velocities(state,sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));
1364 /* With domain decomposition the global state is distributed later */
1365 if (!DOMAINDECOMP(cr))
1367 /* Copy the global state to the local state data structure */
1368 copy_state_nonatomdata(state,state_local);
1372 bcast_state(cr,state,FALSE);
1377 return bThisReplicaExchanged;
1380 void print_replica_exchange_statistics(FILE *fplog,struct gmx_repl_ex *re)
1384 fprintf(fplog,"\nReplica exchange statistics\n");
1388 fprintf(fplog,"Repl %d attempts, %d odd, %d even\n",
1389 re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
1391 fprintf(fplog,"Repl average probabilities:\n");
1392 for(i=1; i<re->nrepl; i++)
1394 if (re->nattempt[i%2] == 0)
1400 re->prob[i] = re->prob_sum[i]/re->nattempt[i%2];
1403 print_ind(fplog,"",re->nrepl,re->ind,NULL);
1404 print_prob(fplog,"",re->nrepl,re->prob);
1406 fprintf(fplog,"Repl number of exchanges:\n");
1407 print_ind(fplog,"",re->nrepl,re->ind,NULL);
1408 print_count(fplog,"",re->nrepl,re->nexchange);
1410 fprintf(fplog,"Repl average number of exchanges:\n");
1411 for(i=1; i<re->nrepl; i++)
1413 if (re->nattempt[i%2] == 0)
1419 re->prob[i] = ((real)re->nexchange[i])/re->nattempt[i%2];
1422 print_ind(fplog,"",re->nrepl,re->ind,NULL);
1423 print_prob(fplog,"",re->nrepl,re->prob);
1425 fprintf(fplog,"\n");
1427 /* print the transition matrix */
1428 print_transition_matrix(fplog,"",re->nrepl,re->nmoves,re->nattempt);