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 * Copyright (c) 2011,2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/gmxlib/network.h"
50 #include "gromacs/math/units.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdlib/main.h"
53 #include "gromacs/mdtypes/commrec.h"
54 #include "gromacs/mdtypes/forcerec.h" // only for gmx_enerdata_t
55 #include "gromacs/mdtypes/inputrec.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/mdtypes/state.h"
58 #include "gromacs/random/threefry.h"
59 #include "gromacs/random/uniformintdistribution.h"
60 #include "gromacs/random/uniformrealdistribution.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/pleasecite.h"
63 #include "gromacs/utility/smalloc.h"
66 #define PROBABILITYCUTOFF 100
67 /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
69 //! Rank in the multisimulaiton
70 #define MSRANK(ms, nodeid) (nodeid)
73 ereTEMP
, ereLAMBDA
, ereENDSINGLE
, ereTL
, ereNR
75 const char *erename
[ereNR
] = { "temperature", "lambda", "end_single_marker", "temperature and lambda"};
76 /* end_single_marker merely notes the end of single variable replica exchange. All types higher than
77 it are multiple replica exchange methods */
78 /* Eventually, should add 'pressure', 'temperature and pressure', 'lambda_and_pressure', 'temperature_lambda_pressure'?;
79 Let's wait until we feel better about the pressure control methods giving exact ensembles. Right now, we assume constant pressure */
81 typedef struct gmx_repl_ex
83 int repl
; /* replica ID */
84 int nrepl
; /* total number of replica */
85 real temp
; /* temperature */
86 int type
; /* replica exchange type from ere enum */
87 real
**q
; /* quantity, e.g. temperature or lambda; first index is ere, second index is replica ID */
88 gmx_bool bNPT
; /* use constant pressure and temperature */
89 real
*pres
; /* replica pressures */
90 int *ind
; /* replica indices */
91 int *allswaps
; /* used for keeping track of all the replica swaps */
92 int nst
; /* replica exchange interval (number of steps) */
93 int nex
; /* number of exchanges per interval */
94 int seed
; /* random seed */
95 int nattempt
[2]; /* number of even and odd replica change attempts */
96 real
*prob_sum
; /* sum of probabilities */
97 int **nmoves
; /* number of moves between replicas i and j */
98 int *nexchange
; /* i-th element of the array is the number of exchanges between replica i-1 and i */
100 /* these are helper arrays for replica exchange; allocated here so they
101 don't have to be allocated each time */
109 /* helper arrays to hold the quantities that are exchanged */
118 static gmx_bool
repl_quantity(const gmx_multisim_t
*ms
,
119 struct gmx_repl_ex
*re
, int ere
, real q
)
125 snew(qall
, ms
->nsim
);
127 gmx_sum_sim(ms
->nsim
, qall
, ms
);
130 for (s
= 1; s
< ms
->nsim
; s
++)
132 if (qall
[s
] != qall
[0])
140 /* Set the replica exchange type and quantities */
143 snew(re
->q
[ere
], re
->nrepl
);
144 for (s
= 0; s
< ms
->nsim
; s
++)
146 re
->q
[ere
][s
] = qall
[s
];
154 init_replica_exchange(FILE *fplog
,
155 const gmx_multisim_t
*ms
,
156 int numAtomsInSystem
,
157 const t_inputrec
*ir
,
158 const ReplicaExchangeParameters
&replExParams
)
162 struct gmx_repl_ex
*re
;
164 gmx_bool bLambda
= FALSE
;
166 fprintf(fplog
, "\nInitializing Replica Exchange\n");
168 if (!isMultiSim(ms
) || ms
->nsim
== 1)
170 gmx_fatal(FARGS
, "Nothing to exchange with only one replica, maybe you forgot to set the -multidir option of mdrun?");
172 if (!EI_DYNAMICS(ir
->eI
))
174 gmx_fatal(FARGS
, "Replica exchange is only supported by dynamical simulations");
175 /* Note that PAR(cr) is defined by cr->nnodes > 1, which is
176 * distinct from isMultiSim(ms). A multi-simulation only runs
177 * with real MPI parallelism, but this does not imply PAR(cr)
180 * Since we are using a dynamical integrator, the only
181 * decomposition is DD, so PAR(cr) and DOMAINDECOMP(cr) are
182 * synonymous. The only way for cr->nnodes > 1 to be true is
183 * if we are using DD. */
189 re
->nrepl
= ms
->nsim
;
190 snew(re
->q
, ereENDSINGLE
);
192 fprintf(fplog
, "Repl There are %d replicas:\n", re
->nrepl
);
194 /* We only check that the number of atoms in the systms match.
195 * This, of course, do not guarantee that the systems are the same,
196 * but it does guarantee that we can perform replica exchange.
198 check_multi_int(fplog
, ms
, numAtomsInSystem
, "the number of atoms", FALSE
);
199 check_multi_int(fplog
, ms
, ir
->eI
, "the integrator", FALSE
);
200 check_multi_int64(fplog
, ms
, ir
->init_step
+ir
->nsteps
, "init_step+nsteps", FALSE
);
201 const int nst
= replExParams
.exchangeInterval
;
202 check_multi_int64(fplog
, ms
, (ir
->init_step
+nst
-1)/nst
,
203 "first exchange step: init_step/-replex", FALSE
);
204 check_multi_int(fplog
, ms
, ir
->etc
, "the temperature coupling", FALSE
);
205 check_multi_int(fplog
, ms
, ir
->opts
.ngtc
,
206 "the number of temperature coupling groups", FALSE
);
207 check_multi_int(fplog
, ms
, ir
->epc
, "the pressure coupling", FALSE
);
208 check_multi_int(fplog
, ms
, ir
->efep
, "free energy", FALSE
);
209 check_multi_int(fplog
, ms
, ir
->fepvals
->n_lambda
, "number of lambda states", FALSE
);
211 re
->temp
= ir
->opts
.ref_t
[0];
212 for (i
= 1; (i
< ir
->opts
.ngtc
); i
++)
214 if (ir
->opts
.ref_t
[i
] != re
->temp
)
216 fprintf(fplog
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
217 fprintf(stderr
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
222 bTemp
= repl_quantity(ms
, re
, ereTEMP
, re
->temp
);
223 if (ir
->efep
!= efepNO
)
225 bLambda
= repl_quantity(ms
, re
, ereLAMBDA
, (real
)ir
->fepvals
->init_fep_state
);
227 if (re
->type
== -1) /* nothing was assigned */
229 gmx_fatal(FARGS
, "The properties of the %d systems are all the same, there is nothing to exchange", re
->nrepl
);
231 if (bLambda
&& bTemp
)
238 please_cite(fplog
, "Sugita1999a");
239 if (ir
->epc
!= epcNO
)
242 fprintf(fplog
, "Repl Using Constant Pressure REMD.\n");
243 please_cite(fplog
, "Okabe2001a");
245 if (ir
->etc
== etcBERENDSEN
)
247 gmx_fatal(FARGS
, "REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
248 ETCOUPLTYPE(ir
->etc
), ETCOUPLTYPE(etcVRESCALE
));
253 if (ir
->fepvals
->delta_lambda
!= 0) /* check this? */
255 gmx_fatal(FARGS
, "delta_lambda is not zero");
260 snew(re
->pres
, re
->nrepl
);
261 if (ir
->epct
== epctSURFACETENSION
)
263 pres
= ir
->ref_p
[ZZ
][ZZ
];
269 for (i
= 0; i
< DIM
; i
++)
271 if (ir
->compress
[i
][i
] != 0)
273 pres
+= ir
->ref_p
[i
][i
];
279 re
->pres
[re
->repl
] = pres
;
280 gmx_sum_sim(re
->nrepl
, re
->pres
, ms
);
283 /* Make an index for increasing replica order */
284 /* only makes sense if one or the other is varying, not both!
285 if both are varying, we trust the order the person gave. */
286 snew(re
->ind
, re
->nrepl
);
287 for (i
= 0; i
< re
->nrepl
; i
++)
292 if (re
->type
< ereENDSINGLE
)
295 for (i
= 0; i
< re
->nrepl
; i
++)
297 for (j
= i
+1; j
< re
->nrepl
; j
++)
299 if (re
->q
[re
->type
][re
->ind
[j
]] < re
->q
[re
->type
][re
->ind
[i
]])
301 /* Unordered replicas are supposed to work, but there
302 * is still an issues somewhere.
303 * Note that at this point still re->ind[i]=i.
305 gmx_fatal(FARGS
, "Replicas with indices %d < %d have %ss %g > %g, please order your replicas on increasing %s",
308 re
->q
[re
->type
][i
], re
->q
[re
->type
][j
],
312 re
->ind
[i
] = re
->ind
[j
];
315 else if (re
->q
[re
->type
][re
->ind
[j
]] == re
->q
[re
->type
][re
->ind
[i
]])
317 gmx_fatal(FARGS
, "Two replicas have identical %ss", erename
[re
->type
]);
323 /* keep track of all the swaps, starting with the initial placement. */
324 snew(re
->allswaps
, re
->nrepl
);
325 for (i
= 0; i
< re
->nrepl
; i
++)
327 re
->allswaps
[i
] = re
->ind
[i
];
333 fprintf(fplog
, "\nReplica exchange in temperature\n");
334 for (i
= 0; i
< re
->nrepl
; i
++)
336 fprintf(fplog
, " %5.1f", re
->q
[re
->type
][re
->ind
[i
]]);
338 fprintf(fplog
, "\n");
341 fprintf(fplog
, "\nReplica exchange in lambda\n");
342 for (i
= 0; i
< re
->nrepl
; i
++)
344 fprintf(fplog
, " %3d", (int)re
->q
[re
->type
][re
->ind
[i
]]);
346 fprintf(fplog
, "\n");
349 fprintf(fplog
, "\nReplica exchange in temperature and lambda state\n");
350 for (i
= 0; i
< re
->nrepl
; i
++)
352 fprintf(fplog
, " %5.1f", re
->q
[ereTEMP
][re
->ind
[i
]]);
354 fprintf(fplog
, "\n");
355 for (i
= 0; i
< re
->nrepl
; i
++)
357 fprintf(fplog
, " %5d", (int)re
->q
[ereLAMBDA
][re
->ind
[i
]]);
359 fprintf(fplog
, "\n");
362 gmx_incons("Unknown replica exchange quantity");
366 fprintf(fplog
, "\nRepl p");
367 for (i
= 0; i
< re
->nrepl
; i
++)
369 fprintf(fplog
, " %5.2f", re
->pres
[re
->ind
[i
]]);
372 for (i
= 0; i
< re
->nrepl
; i
++)
374 if ((i
> 0) && (re
->pres
[re
->ind
[i
]] < re
->pres
[re
->ind
[i
-1]]))
376 fprintf(fplog
, "\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
377 fprintf(stderr
, "\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
382 if (replExParams
.randomSeed
== -1)
386 re
->seed
= static_cast<int>(gmx::makeRandomSeed());
392 gmx_sumi_sim(1, &(re
->seed
), ms
);
396 re
->seed
= replExParams
.randomSeed
;
398 fprintf(fplog
, "\nReplica exchange interval: %d\n", re
->nst
);
399 fprintf(fplog
, "\nReplica random seed: %d\n", re
->seed
);
404 snew(re
->prob_sum
, re
->nrepl
);
405 snew(re
->nexchange
, re
->nrepl
);
406 snew(re
->nmoves
, re
->nrepl
);
407 for (i
= 0; i
< re
->nrepl
; i
++)
409 snew(re
->nmoves
[i
], re
->nrepl
);
411 fprintf(fplog
, "Replica exchange information below: ex and x = exchange, pr = probability\n");
413 /* generate space for the helper functions so we don't have to snew each time */
415 snew(re
->destinations
, re
->nrepl
);
416 snew(re
->incycle
, re
->nrepl
);
417 snew(re
->tmpswap
, re
->nrepl
);
418 snew(re
->cyclic
, re
->nrepl
);
419 snew(re
->order
, re
->nrepl
);
420 for (i
= 0; i
< re
->nrepl
; i
++)
422 snew(re
->cyclic
[i
], re
->nrepl
+1);
423 snew(re
->order
[i
], re
->nrepl
);
425 /* allocate space for the functions storing the data for the replicas */
426 /* not all of these arrays needed in all cases, but they don't take
427 up much space, since the max size is nrepl**2 */
428 snew(re
->prob
, re
->nrepl
);
429 snew(re
->bEx
, re
->nrepl
);
430 snew(re
->beta
, re
->nrepl
);
431 snew(re
->Vol
, re
->nrepl
);
432 snew(re
->Epot
, re
->nrepl
);
433 snew(re
->de
, re
->nrepl
);
434 for (i
= 0; i
< re
->nrepl
; i
++)
436 snew(re
->de
[i
], re
->nrepl
);
438 re
->nex
= replExParams
.numExchanges
;
442 static void exchange_reals(const gmx_multisim_t gmx_unused
*ms
, int gmx_unused b
, real
*v
, int n
)
452 MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
453 buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
454 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
459 MPI_Isend(v
, n
*sizeof(real
), MPI_BYTE
, MSRANK(ms
, b
), 0,
460 ms
->mpi_comm_masters
, &mpi_req
);
461 MPI_Recv(buf
, n
*sizeof(real
), MPI_BYTE
, MSRANK(ms
, b
), 0,
462 ms
->mpi_comm_masters
, MPI_STATUS_IGNORE
);
463 MPI_Wait(&mpi_req
, MPI_STATUS_IGNORE
);
466 for (i
= 0; i
< n
; i
++)
475 static void exchange_doubles(const gmx_multisim_t gmx_unused
*ms
, int gmx_unused b
, double *v
, int n
)
485 MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
486 buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
487 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
492 MPI_Isend(v
, n
*sizeof(double), MPI_BYTE
, MSRANK(ms
, b
), 0,
493 ms
->mpi_comm_masters
, &mpi_req
);
494 MPI_Recv(buf
, n
*sizeof(double), MPI_BYTE
, MSRANK(ms
, b
), 0,
495 ms
->mpi_comm_masters
, MPI_STATUS_IGNORE
);
496 MPI_Wait(&mpi_req
, MPI_STATUS_IGNORE
);
499 for (i
= 0; i
< n
; i
++)
507 static void exchange_rvecs(const gmx_multisim_t gmx_unused
*ms
, int gmx_unused b
, rvec
*v
, int n
)
517 MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
518 buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
519 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
524 MPI_Isend(v
[0], n
*sizeof(rvec
), MPI_BYTE
, MSRANK(ms
, b
), 0,
525 ms
->mpi_comm_masters
, &mpi_req
);
526 MPI_Recv(buf
[0], n
*sizeof(rvec
), MPI_BYTE
, MSRANK(ms
, b
), 0,
527 ms
->mpi_comm_masters
, MPI_STATUS_IGNORE
);
528 MPI_Wait(&mpi_req
, MPI_STATUS_IGNORE
);
531 for (i
= 0; i
< n
; i
++)
533 copy_rvec(buf
[i
], v
[i
]);
539 static void exchange_state(const gmx_multisim_t
*ms
, int b
, t_state
*state
)
541 /* When t_state changes, this code should be updated. */
543 ngtc
= state
->ngtc
* state
->nhchainlength
;
544 nnhpres
= state
->nnhpres
* state
->nhchainlength
;
545 exchange_rvecs(ms
, b
, state
->box
, DIM
);
546 exchange_rvecs(ms
, b
, state
->box_rel
, DIM
);
547 exchange_rvecs(ms
, b
, state
->boxv
, DIM
);
548 exchange_reals(ms
, b
, &(state
->veta
), 1);
549 exchange_reals(ms
, b
, &(state
->vol0
), 1);
550 exchange_rvecs(ms
, b
, state
->svir_prev
, DIM
);
551 exchange_rvecs(ms
, b
, state
->fvir_prev
, DIM
);
552 exchange_rvecs(ms
, b
, state
->pres_prev
, DIM
);
553 exchange_doubles(ms
, b
, state
->nosehoover_xi
.data(), ngtc
);
554 exchange_doubles(ms
, b
, state
->nosehoover_vxi
.data(), ngtc
);
555 exchange_doubles(ms
, b
, state
->nhpres_xi
.data(), nnhpres
);
556 exchange_doubles(ms
, b
, state
->nhpres_vxi
.data(), nnhpres
);
557 exchange_doubles(ms
, b
, state
->therm_integral
.data(), state
->ngtc
);
558 exchange_doubles(ms
, b
, &state
->baros_integral
, 1);
559 exchange_rvecs(ms
, b
, as_rvec_array(state
->x
.data()), state
->natoms
);
560 exchange_rvecs(ms
, b
, as_rvec_array(state
->v
.data()), state
->natoms
);
563 static void copy_state_serial(const t_state
*src
, t_state
*dest
)
567 /* Currently the local state is always a pointer to the global
568 * in serial, so we should never end up here.
569 * TODO: Implement a (trivial) t_state copy once converted to C++.
571 GMX_RELEASE_ASSERT(false, "State copying is currently not implemented in replica exchange");
575 static void scale_velocities(t_state
*state
, real fac
)
579 if (as_rvec_array(state
->v
.data()))
581 for (i
= 0; i
< state
->natoms
; i
++)
583 svmul(fac
, state
->v
[i
], state
->v
[i
]);
588 static void print_transition_matrix(FILE *fplog
, int n
, int **nmoves
, int *nattempt
)
593 ntot
= nattempt
[0] + nattempt
[1];
594 fprintf(fplog
, "\n");
595 fprintf(fplog
, "Repl");
596 for (i
= 0; i
< n
; i
++)
598 fprintf(fplog
, " "); /* put the title closer to the center */
600 fprintf(fplog
, "Empirical Transition Matrix\n");
602 fprintf(fplog
, "Repl");
603 for (i
= 0; i
< n
; i
++)
605 fprintf(fplog
, "%8d", (i
+1));
607 fprintf(fplog
, "\n");
609 for (i
= 0; i
< n
; i
++)
611 fprintf(fplog
, "Repl");
612 for (j
= 0; j
< n
; j
++)
615 if (nmoves
[i
][j
] > 0)
617 Tprint
= nmoves
[i
][j
]/(2.0*ntot
);
619 fprintf(fplog
, "%8.4f", Tprint
);
621 fprintf(fplog
, "%3d\n", i
);
625 static void print_ind(FILE *fplog
, const char *leg
, int n
, int *ind
, gmx_bool
*bEx
)
629 fprintf(fplog
, "Repl %2s %2d", leg
, ind
[0]);
630 for (i
= 1; i
< n
; i
++)
632 fprintf(fplog
, " %c %2d", (bEx
!= nullptr && bEx
[i
]) ? 'x' : ' ', ind
[i
]);
634 fprintf(fplog
, "\n");
637 static void print_allswitchind(FILE *fplog
, int n
, int *pind
, int *allswaps
, int *tmpswap
)
641 for (i
= 0; i
< n
; i
++)
643 tmpswap
[i
] = allswaps
[i
];
645 for (i
= 0; i
< n
; i
++)
647 allswaps
[i
] = tmpswap
[pind
[i
]];
650 fprintf(fplog
, "\nAccepted Exchanges: ");
651 for (i
= 0; i
< n
; i
++)
653 fprintf(fplog
, "%d ", pind
[i
]);
655 fprintf(fplog
, "\n");
657 /* the "Order After Exchange" is the state label corresponding to the configuration that
658 started in state listed in order, i.e.
663 configuration starting in simulation 3 is now in simulation 0,
664 configuration starting in simulation 0 is now in simulation 1,
665 configuration starting in simulation 1 is now in simulation 2,
666 configuration starting in simulation 2 is now in simulation 3
668 fprintf(fplog
, "Order After Exchange: ");
669 for (i
= 0; i
< n
; i
++)
671 fprintf(fplog
, "%d ", allswaps
[i
]);
673 fprintf(fplog
, "\n\n");
676 static void print_prob(FILE *fplog
, const char *leg
, int n
, real
*prob
)
681 fprintf(fplog
, "Repl %2s ", leg
);
682 for (i
= 1; i
< n
; i
++)
686 sprintf(buf
, "%4.2f", prob
[i
]);
687 fprintf(fplog
, " %3s", buf
[0] == '1' ? "1.0" : buf
+1);
694 fprintf(fplog
, "\n");
697 static void print_count(FILE *fplog
, const char *leg
, int n
, int *count
)
701 fprintf(fplog
, "Repl %2s ", leg
);
702 for (i
= 1; i
< n
; i
++)
704 fprintf(fplog
, " %4d", count
[i
]);
706 fprintf(fplog
, "\n");
709 static real
calc_delta(FILE *fplog
, gmx_bool bPrint
, struct gmx_repl_ex
*re
, int a
, int b
, int ap
, int bp
)
712 real ediff
, dpV
, delta
= 0;
713 real
*Epot
= re
->Epot
;
716 real
*beta
= re
->beta
;
718 /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce
719 to the non permuted case */
725 * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
727 ediff
= Epot
[b
] - Epot
[a
];
728 delta
= -(beta
[bp
] - beta
[ap
])*ediff
;
731 /* two cases: when we are permuted, and not. */
733 ediff = E_new - E_old
734 = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
735 = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
736 = de[b][a] + de[a][b] */
739 ediff = E_new - E_old
740 = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
741 = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
742 = [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)]
743 = [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)]
744 = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]) */
745 /* but, in the current code implementation, we flip configurations, not indices . . .
746 So let's examine that.
747 = [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)]
748 = [H_b(x_ap) - H_a(x_ap)] + [H_a(x_bp) - H_b(x_pb)]
749 = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp]
750 So, if we exchange b<=> bp and a<=> ap, we return to the same result.
751 So the simple solution is to flip the
752 position of perturbed and original indices in the tests.
755 ediff
= (de
[bp
][a
] - de
[ap
][a
]) + (de
[ap
][b
] - de
[bp
][b
]);
756 delta
= ediff
*beta
[a
]; /* assume all same temperature in this case */
760 /* delta = reduced E_new - reduced E_old
761 = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
762 = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
763 = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
764 [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
765 = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
766 beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
767 = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
768 /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
769 /* permuted (big breath!) */
770 /* delta = reduced E_new - reduced E_old
771 = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
772 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
773 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
774 - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
775 - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
776 = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
777 [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
778 + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
779 = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
780 [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
781 + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b))
782 = ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b] - beta_bp de[bp][b])
783 + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */
784 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
]);
787 gmx_incons("Unknown replica exchange quantity");
791 fprintf(fplog
, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a
, b
, delta
);
795 /* revist the calculation for 5.0. Might be some improvements. */
796 dpV
= (beta
[ap
]*re
->pres
[ap
]-beta
[bp
]*re
->pres
[bp
])*(Vol
[b
]-Vol
[a
])/PRESFAC
;
799 fprintf(fplog
, " dpV = %10.3e d = %10.3e\n", dpV
, delta
+ dpV
);
807 test_for_replica_exchange(FILE *fplog
,
808 const gmx_multisim_t
*ms
,
809 struct gmx_repl_ex
*re
,
810 const gmx_enerdata_t
*enerd
,
815 int m
, i
, j
, a
, b
, ap
, bp
, i0
, i1
, tmp
;
817 gmx_bool bPrint
, bMultiEx
;
818 gmx_bool
*bEx
= re
->bEx
;
819 real
*prob
= re
->prob
;
820 int *pind
= re
->destinations
; /* permuted index */
821 gmx_bool bEpot
= FALSE
;
822 gmx_bool bDLambda
= FALSE
;
823 gmx_bool bVol
= FALSE
;
824 gmx::ThreeFry2x64
<64> rng(re
->seed
, gmx::RandomDomain::ReplicaExchange
);
825 gmx::UniformRealDistribution
<real
> uniformRealDist
;
826 gmx::UniformIntDistribution
<int> uniformNreplDist(0, re
->nrepl
-1);
828 bMultiEx
= (re
->nex
> 1); /* multiple exchanges at each state */
829 fprintf(fplog
, "Replica exchange at step %" GMX_PRId64
" time %.5f\n", step
, time
);
833 for (i
= 0; i
< re
->nrepl
; i
++)
838 re
->Vol
[re
->repl
] = vol
;
840 if ((re
->type
== ereTEMP
|| re
->type
== ereTL
))
842 for (i
= 0; i
< re
->nrepl
; i
++)
847 re
->Epot
[re
->repl
] = enerd
->term
[F_EPOT
];
848 /* temperatures of different states*/
849 for (i
= 0; i
< re
->nrepl
; i
++)
851 re
->beta
[i
] = 1.0/(re
->q
[ereTEMP
][i
]*BOLTZ
);
856 for (i
= 0; i
< re
->nrepl
; i
++)
858 re
->beta
[i
] = 1.0/(re
->temp
*BOLTZ
); /* we have a single temperature */
861 if (re
->type
== ereLAMBDA
|| re
->type
== ereTL
)
864 /* lambda differences. */
865 /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
866 minus the energy of the jth simulation in the jth Hamiltonian */
867 for (i
= 0; i
< re
->nrepl
; i
++)
869 for (j
= 0; j
< re
->nrepl
; j
++)
874 for (i
= 0; i
< re
->nrepl
; i
++)
876 re
->de
[i
][re
->repl
] = (enerd
->enerpart_lambda
[(int)re
->q
[ereLAMBDA
][i
]+1]-enerd
->enerpart_lambda
[0]);
880 /* now actually do the communication */
883 gmx_sum_sim(re
->nrepl
, re
->Vol
, ms
);
887 gmx_sum_sim(re
->nrepl
, re
->Epot
, ms
);
891 for (i
= 0; i
< re
->nrepl
; i
++)
893 gmx_sum_sim(re
->nrepl
, re
->de
[i
], ms
);
897 /* make a duplicate set of indices for shuffling */
898 for (i
= 0; i
< re
->nrepl
; i
++)
900 pind
[i
] = re
->ind
[i
];
903 rng
.restart( step
, 0 );
907 /* multiple random switch exchange */
911 for (i
= 0; i
< re
->nex
+ nself
; i
++)
913 // For now this is superfluous, but just in case we ever add more
914 // calls in different branches it is safer to always reset the distribution.
915 uniformNreplDist
.reset();
917 /* randomly select a pair */
918 /* in theory, could reduce this by identifying only which switches had a nonneglibible
919 probability of occurring (log p > -100) and only operate on those switches */
920 /* find out which state it is from, and what label that state currently has. Likely
921 more work that useful. */
922 i0
= uniformNreplDist(rng
);
923 i1
= uniformNreplDist(rng
);
927 continue; /* self-exchange, back up and do it again */
930 a
= re
->ind
[i0
]; /* what are the indices of these states? */
935 bPrint
= FALSE
; /* too noisy */
936 /* calculate the energy difference */
937 /* if the code changes to flip the STATES, rather than the configurations,
938 use the commented version of the code */
939 /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */
940 delta
= calc_delta(fplog
, bPrint
, re
, ap
, bp
, a
, b
);
942 /* we actually only use the first space in the prob and bEx array,
943 since there are actually many switches between pairs. */
953 if (delta
> PROBABILITYCUTOFF
)
959 prob
[0] = exp(-delta
);
961 // roll a number to determine if accepted. For now it is superfluous to
962 // reset, but just in case we ever add more calls in different branches
963 // it is safer to always reset the distribution.
964 uniformRealDist
.reset();
965 bEx
[0] = uniformRealDist(rng
) < prob
[0];
967 re
->prob_sum
[0] += prob
[0];
971 /* swap the states */
977 re
->nattempt
[0]++; /* keep track of total permutation trials here */
978 print_allswitchind(fplog
, re
->nrepl
, pind
, re
->allswaps
, re
->tmpswap
);
982 /* standard nearest neighbor replica exchange */
984 m
= (step
/ re
->nst
) % 2;
985 for (i
= 1; i
< re
->nrepl
; i
++)
990 bPrint
= (re
->repl
== a
|| re
->repl
== b
);
993 delta
= calc_delta(fplog
, bPrint
, re
, a
, b
, a
, b
);
1002 if (delta
> PROBABILITYCUTOFF
)
1008 prob
[i
] = exp(-delta
);
1010 // roll a number to determine if accepted. For now it is superfluous to
1011 // reset, but just in case we ever add more calls in different branches
1012 // it is safer to always reset the distribution.
1013 uniformRealDist
.reset();
1014 bEx
[i
] = uniformRealDist(rng
) < prob
[i
];
1016 re
->prob_sum
[i
] += prob
[i
];
1020 /* swap these two */
1022 pind
[i
-1] = pind
[i
];
1024 re
->nexchange
[i
]++; /* statistics for back compatibility */
1033 /* print some statistics */
1034 print_ind(fplog
, "ex", re
->nrepl
, re
->ind
, bEx
);
1035 print_prob(fplog
, "pr", re
->nrepl
, prob
);
1036 fprintf(fplog
, "\n");
1040 /* record which moves were made and accepted */
1041 for (i
= 0; i
< re
->nrepl
; i
++)
1043 re
->nmoves
[re
->ind
[i
]][pind
[i
]] += 1;
1044 re
->nmoves
[pind
[i
]][re
->ind
[i
]] += 1;
1046 fflush(fplog
); /* make sure we can see what the last exchange was */
1050 cyclic_decomposition(const int *destinations
,
1059 for (i
= 0; i
< nrepl
; i
++)
1063 for (i
= 0; i
< nrepl
; i
++) /* one cycle for each replica */
1074 for (j
= 0; j
< nrepl
; j
++) /* potentially all cycles are part, but we will break first */
1076 p
= destinations
[p
]; /* start permuting */
1084 break; /* we've reached the original element, the cycle is complete, and we marked the end. */
1088 cyclic
[i
][c
] = p
; /* each permutation gives a new member of the cycle */
1094 *nswap
= maxlen
- 1;
1098 for (i
= 0; i
< nrepl
; i
++)
1100 fprintf(debug
, "Cycle %d:", i
);
1101 for (j
= 0; j
< nrepl
; j
++)
1103 if (cyclic
[i
][j
] < 0)
1107 fprintf(debug
, "%2d", cyclic
[i
][j
]);
1109 fprintf(debug
, "\n");
1116 compute_exchange_order(int **cyclic
,
1123 for (j
= 0; j
< maxswap
; j
++)
1125 for (i
= 0; i
< nrepl
; i
++)
1127 if (cyclic
[i
][j
+1] >= 0)
1129 order
[cyclic
[i
][j
+1]][j
] = cyclic
[i
][j
];
1130 order
[cyclic
[i
][j
]][j
] = cyclic
[i
][j
+1];
1133 for (i
= 0; i
< nrepl
; i
++)
1135 if (order
[i
][j
] < 0)
1137 order
[i
][j
] = i
; /* if it's not exchanging, it should stay this round*/
1144 fprintf(debug
, "Replica Exchange Order\n");
1145 for (i
= 0; i
< nrepl
; i
++)
1147 fprintf(debug
, "Replica %d:", i
);
1148 for (j
= 0; j
< maxswap
; j
++)
1150 if (order
[i
][j
] < 0)
1154 fprintf(debug
, "%2d", order
[i
][j
]);
1156 fprintf(debug
, "\n");
1163 prepare_to_do_exchange(struct gmx_repl_ex
*re
,
1164 const int replica_id
,
1166 gmx_bool
*bThisReplicaExchanged
)
1169 /* Hold the cyclic decomposition of the (multiple) replica
1171 gmx_bool bAnyReplicaExchanged
= FALSE
;
1172 *bThisReplicaExchanged
= FALSE
;
1174 for (i
= 0; i
< re
->nrepl
; i
++)
1176 if (re
->destinations
[i
] != re
->ind
[i
])
1178 /* only mark as exchanged if the index has been shuffled */
1179 bAnyReplicaExchanged
= TRUE
;
1183 if (bAnyReplicaExchanged
)
1185 /* reinitialize the placeholder arrays */
1186 for (i
= 0; i
< re
->nrepl
; i
++)
1188 for (j
= 0; j
< re
->nrepl
; j
++)
1190 re
->cyclic
[i
][j
] = -1;
1191 re
->order
[i
][j
] = -1;
1195 /* Identify the cyclic decomposition of the permutation (very
1196 * fast if neighbor replica exchange). */
1197 cyclic_decomposition(re
->destinations
, re
->cyclic
, re
->incycle
, re
->nrepl
, maxswap
);
1199 /* Now translate the decomposition into a replica exchange
1200 * order at each step. */
1201 compute_exchange_order(re
->cyclic
, re
->order
, re
->nrepl
, *maxswap
);
1203 /* Did this replica do any exchange at any point? */
1204 for (j
= 0; j
< *maxswap
; j
++)
1206 if (replica_id
!= re
->order
[replica_id
][j
])
1208 *bThisReplicaExchanged
= TRUE
;
1215 gmx_bool
replica_exchange(FILE *fplog
, const t_commrec
*cr
,
1216 const gmx_multisim_t
*ms
, struct gmx_repl_ex
*re
,
1217 t_state
*state
, const gmx_enerdata_t
*enerd
,
1218 t_state
*state_local
, gmx_int64_t step
, real time
)
1222 int exchange_partner
;
1224 /* Number of rounds of exchanges needed to deal with any multiple
1226 /* Where each replica ends up after the exchange attempt(s). */
1227 /* The order in which multiple exchanges will occur. */
1228 gmx_bool bThisReplicaExchanged
= FALSE
;
1232 replica_id
= re
->repl
;
1233 test_for_replica_exchange(fplog
, ms
, re
, enerd
, det(state_local
->box
), step
, time
);
1234 prepare_to_do_exchange(re
, replica_id
, &maxswap
, &bThisReplicaExchanged
);
1236 /* Do intra-simulation broadcast so all processors belonging to
1237 * each simulation know whether they need to participate in
1238 * collecting the state. Otherwise, they might as well get on with
1239 * the next thing to do. */
1240 if (DOMAINDECOMP(cr
))
1243 MPI_Bcast(&bThisReplicaExchanged
, sizeof(gmx_bool
), MPI_BYTE
, MASTERRANK(cr
),
1244 cr
->mpi_comm_mygroup
);
1248 if (bThisReplicaExchanged
)
1250 /* Exchange the states */
1251 /* Collect the global state on the master node */
1252 if (DOMAINDECOMP(cr
))
1254 dd_collect_state(cr
->dd
, state_local
, state
);
1258 copy_state_serial(state_local
, state
);
1263 /* There will be only one swap cycle with standard replica
1264 * exchange, but there may be multiple swap cycles if we
1265 * allow multiple swaps. */
1267 for (j
= 0; j
< maxswap
; j
++)
1269 exchange_partner
= re
->order
[replica_id
][j
];
1271 if (exchange_partner
!= replica_id
)
1273 /* Exchange the global states between the master nodes */
1276 fprintf(debug
, "Exchanging %d with %d\n", replica_id
, exchange_partner
);
1278 exchange_state(ms
, exchange_partner
, state
);
1281 /* For temperature-type replica exchange, we need to scale
1282 * the velocities. */
1283 if (re
->type
== ereTEMP
|| re
->type
== ereTL
)
1285 scale_velocities(state
, sqrt(re
->q
[ereTEMP
][replica_id
]/re
->q
[ereTEMP
][re
->destinations
[replica_id
]]));
1290 /* With domain decomposition the global state is distributed later */
1291 if (!DOMAINDECOMP(cr
))
1293 /* Copy the global state to the local state data structure */
1294 copy_state_serial(state
, state_local
);
1298 return bThisReplicaExchanged
;
1301 void print_replica_exchange_statistics(FILE *fplog
, struct gmx_repl_ex
*re
)
1305 fprintf(fplog
, "\nReplica exchange statistics\n");
1309 fprintf(fplog
, "Repl %d attempts, %d odd, %d even\n",
1310 re
->nattempt
[0]+re
->nattempt
[1], re
->nattempt
[1], re
->nattempt
[0]);
1312 fprintf(fplog
, "Repl average probabilities:\n");
1313 for (i
= 1; i
< re
->nrepl
; i
++)
1315 if (re
->nattempt
[i
%2] == 0)
1321 re
->prob
[i
] = re
->prob_sum
[i
]/re
->nattempt
[i
%2];
1324 print_ind(fplog
, "", re
->nrepl
, re
->ind
, nullptr);
1325 print_prob(fplog
, "", re
->nrepl
, re
->prob
);
1327 fprintf(fplog
, "Repl number of exchanges:\n");
1328 print_ind(fplog
, "", re
->nrepl
, re
->ind
, nullptr);
1329 print_count(fplog
, "", re
->nrepl
, re
->nexchange
);
1331 fprintf(fplog
, "Repl average number of exchanges:\n");
1332 for (i
= 1; i
< re
->nrepl
; i
++)
1334 if (re
->nattempt
[i
%2] == 0)
1340 re
->prob
[i
] = ((real
)re
->nexchange
[i
])/re
->nattempt
[i
%2];
1343 print_ind(fplog
, "", re
->nrepl
, re
->ind
, nullptr);
1344 print_prob(fplog
, "", re
->nrepl
, re
->prob
);
1346 fprintf(fplog
, "\n");
1348 /* print the transition matrix */
1349 print_transition_matrix(fplog
, re
->nrepl
, re
->nmoves
, re
->nattempt
);