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-2019,2020,2021, 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.
40 * \brief Implements the replica exchange routines.
42 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
43 * \author Mark Abraham <mark.j.abraham@gmail.com>
44 * \ingroup module_mdrun
48 #include "gromacs/utility/enumerationhelpers.h"
49 #include "replicaexchange.h"
57 #include "gromacs/domdec/collect.h"
58 #include "gromacs/gmxlib/network.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/mdrunutility/multisim.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/enerdata.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/mdtypes/state.h"
67 #include "gromacs/random/threefry.h"
68 #include "gromacs/random/uniformintdistribution.h"
69 #include "gromacs/random/uniformrealdistribution.h"
70 #include "gromacs/utility/enumerationhelpers.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/pleasecite.h"
73 #include "gromacs/utility/smalloc.h"
75 //! Helps cut off probability values.
76 constexpr int c_probabilityCutoff = 100;
78 /* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
80 //! Rank in the multisimulation
81 #define MSRANK(ms, nodeid) (nodeid)
83 //! Enum for replica exchange flavours
84 enum class ReplicaExchangeType : int
92 /*! \brief Strings describing replica exchange flavours.
94 * end_single_marker merely notes the end of single variable replica
95 * exchange. All types higher than it are multiple replica exchange
98 * Eventually, should add 'pressure', 'temperature and pressure',
99 * 'lambda_and_pressure', 'temperature_lambda_pressure'?; Let's wait
100 * until we feel better about the pressure control methods giving
101 * exact ensembles. Right now, we assume constant pressure */
102 static const char* enumValueToString(ReplicaExchangeType enumValue)
104 constexpr gmx::EnumerationArray<ReplicaExchangeType, const char*> replicateExchangeTypeNames = {
105 "temperature", "lambda", "end_single_marker", "temperature and lambda"
107 return replicateExchangeTypeNames[enumValue];
110 //! Working data for replica exchange.
115 //! Total number of replica
119 //! Replica exchange type from ReplicaExchangeType enum
120 ReplicaExchangeType type;
121 //! Quantity, e.g. temperature or lambda; first index is ere, second index is replica ID
122 gmx::EnumerationArray<ReplicaExchangeType, real*> q;
123 //! Use constant pressure and temperature
125 //! Replica pressures
129 //! Used for keeping track of all the replica swaps
131 //! Replica exchange interval (number of steps)
133 //! Number of exchanges per interval
137 //! Number of even and odd replica change attempts
139 //! Sum of probabilities
141 //! Number of moves between replicas i and j
143 //! i-th element of the array is the number of exchanges between replica i-1 and i
146 /*! \brief Helper arrays for replica exchange; allocated here
147 * so they don't have to be allocated each time */
157 //! Helper arrays to hold the quantities that are exchanged.
167 // TODO We should add Doxygen here some time.
170 static gmx_bool repl_quantity(const gmx_multisim_t* ms, struct gmx_repl_ex* re, ReplicaExchangeType ere, real q)
176 snew(qall, ms->numSimulations_);
178 gmx_sum_sim(ms->numSimulations_, qall, ms);
181 for (s = 1; s < ms->numSimulations_; s++)
183 if (qall[s] != qall[0])
191 /* Set the replica exchange type and quantities */
194 snew(re->q[ere], re->nrepl);
195 for (s = 0; s < ms->numSimulations_; s++)
197 re->q[ere][s] = qall[s];
204 gmx_repl_ex_t init_replica_exchange(FILE* fplog,
205 const gmx_multisim_t* ms,
206 int numAtomsInSystem,
207 const t_inputrec* ir,
208 const ReplicaExchangeParameters& replExParams)
212 struct gmx_repl_ex* re;
214 gmx_bool bLambda = FALSE;
216 fprintf(fplog, "\nInitializing Replica Exchange\n");
218 if (!isMultiSim(ms) || ms->numSimulations_ == 1)
221 "Nothing to exchange with only one replica, maybe you forgot to set the "
222 "-multidir option of mdrun?");
224 if (replExParams.numExchanges < 0)
226 gmx_fatal(FARGS, "Replica exchange number of exchanges needs to be positive");
229 if (!EI_DYNAMICS(ir->eI))
231 gmx_fatal(FARGS, "Replica exchange is only supported by dynamical simulations");
232 /* Note that PAR(cr) is defined by cr->nnodes > 1, which is
233 * distinct from isMultiSim(ms). A multi-simulation only runs
234 * with real MPI parallelism, but this does not imply PAR(cr)
237 * Since we are using a dynamical integrator, the only
238 * decomposition is DD, so PAR(cr) and haveDDAtomOrdering(*cr) are
239 * synonymous. The only way for cr->nnodes > 1 to be true is
240 * if we are using DD. */
245 re->repl = ms->simulationIndex_;
246 re->nrepl = ms->numSimulations_;
248 fprintf(fplog, "Repl There are %d replicas:\n", re->nrepl);
250 /* We only check that the number of atoms in the systms match.
251 * This, of course, do not guarantee that the systems are the same,
252 * but it does guarantee that we can perform replica exchange.
254 check_multi_int(fplog, ms, numAtomsInSystem, "the number of atoms", FALSE);
255 check_multi_int(fplog, ms, static_cast<int>(ir->eI), "the integrator", FALSE);
256 check_multi_int64(fplog, ms, ir->init_step + ir->nsteps, "init_step+nsteps", FALSE);
257 const int nst = replExParams.exchangeInterval;
259 fplog, ms, (ir->init_step + nst - 1) / nst, "first exchange step: init_step/-replex", FALSE);
260 check_multi_int(fplog, ms, static_cast<int>(ir->etc), "the temperature coupling", FALSE);
261 check_multi_int(fplog, ms, ir->opts.ngtc, "the number of temperature coupling groups", FALSE);
262 check_multi_int(fplog, ms, static_cast<int>(ir->epc), "the pressure coupling", FALSE);
263 check_multi_int(fplog, ms, static_cast<int>(ir->efep), "free energy", FALSE);
264 check_multi_int(fplog, ms, ir->fepvals->n_lambda, "number of lambda states", FALSE);
266 re->temp = ir->opts.ref_t[0];
267 for (i = 1; (i < ir->opts.ngtc); i++)
269 if (ir->opts.ref_t[i] != re->temp)
272 "\nWARNING: The temperatures of the different temperature coupling groups are "
273 "not identical\n\n");
275 "\nWARNING: The temperatures of the different temperature coupling groups are "
276 "not identical\n\n");
280 re->type = ReplicaExchangeType::Count;
281 bTemp = repl_quantity(ms, re, ReplicaExchangeType::Temperature, re->temp);
282 if (ir->efep != FreeEnergyPerturbationType::No)
284 bLambda = repl_quantity(
285 ms, re, ReplicaExchangeType::Lambda, static_cast<real>(ir->fepvals->init_fep_state));
287 if (re->type == ReplicaExchangeType::Count) /* nothing was assigned */
290 "The properties of the %d systems are all the same, there is nothing to exchange",
293 if (bLambda && bTemp)
295 re->type = ReplicaExchangeType::TemperatureLambda;
300 please_cite(fplog, "Sugita1999a");
301 if (ir->epc != PressureCoupling::No)
304 fprintf(fplog, "Repl Using Constant Pressure REMD.\n");
305 please_cite(fplog, "Okabe2001a");
307 if (ir->etc == TemperatureCoupling::Berendsen)
310 "REMD with the %s thermostat does not produce correct potential energy "
311 "distributions, consider using the %s thermostat instead",
312 enumValueToString(ir->etc),
313 enumValueToString(TemperatureCoupling::VRescale));
318 if (ir->fepvals->delta_lambda != 0) /* check this? */
320 gmx_fatal(FARGS, "delta_lambda is not zero");
325 snew(re->pres, re->nrepl);
326 if (ir->epct == PressureCouplingType::SurfaceTension)
328 pres = ir->ref_p[ZZ][ZZ];
334 for (i = 0; i < DIM; i++)
336 if (ir->compress[i][i] != 0)
338 pres += ir->ref_p[i][i];
344 re->pres[re->repl] = pres;
345 gmx_sum_sim(re->nrepl, re->pres, ms);
348 /* Make an index for increasing replica order */
349 /* only makes sense if one or the other is varying, not both!
350 if both are varying, we trust the order the person gave. */
351 snew(re->ind, re->nrepl);
352 for (i = 0; i < re->nrepl; i++)
357 if (re->type < ReplicaExchangeType::EndSingle)
360 for (i = 0; i < re->nrepl; i++)
362 for (j = i + 1; j < re->nrepl; j++)
364 if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]])
366 /* Unordered replicas are supposed to work, but there
367 * is still an issues somewhere.
368 * Note that at this point still re->ind[i]=i.
371 "Replicas with indices %d < %d have %ss %g > %g, please order your "
372 "replicas on increasing %s",
375 enumValueToString(re->type),
378 enumValueToString(re->type));
380 else if (re->q[re->type][re->ind[j]] == re->q[re->type][re->ind[i]])
382 gmx_fatal(FARGS, "Two replicas have identical %ss", enumValueToString(re->type));
388 /* keep track of all the swaps, starting with the initial placement. */
389 snew(re->allswaps, re->nrepl);
390 for (i = 0; i < re->nrepl; i++)
392 re->allswaps[i] = re->ind[i];
397 case ReplicaExchangeType::Temperature:
398 fprintf(fplog, "\nReplica exchange in temperature\n");
399 for (i = 0; i < re->nrepl; i++)
401 fprintf(fplog, " %5.1f", re->q[re->type][re->ind[i]]);
403 fprintf(fplog, "\n");
405 case ReplicaExchangeType::Lambda:
406 fprintf(fplog, "\nReplica exchange in lambda\n");
407 for (i = 0; i < re->nrepl; i++)
409 fprintf(fplog, " %3d", static_cast<int>(re->q[re->type][re->ind[i]]));
411 fprintf(fplog, "\n");
413 case ReplicaExchangeType::TemperatureLambda:
414 fprintf(fplog, "\nReplica exchange in temperature and lambda state\n");
415 for (i = 0; i < re->nrepl; i++)
417 fprintf(fplog, " %5.1f", re->q[ReplicaExchangeType::Temperature][re->ind[i]]);
419 fprintf(fplog, "\n");
420 for (i = 0; i < re->nrepl; i++)
422 fprintf(fplog, " %5d", static_cast<int>(re->q[ReplicaExchangeType::Lambda][re->ind[i]]));
424 fprintf(fplog, "\n");
426 default: gmx_incons("Unknown replica exchange quantity");
430 fprintf(fplog, "\nRepl p");
431 for (i = 0; i < re->nrepl; i++)
433 fprintf(fplog, " %5.2f", re->pres[re->ind[i]]);
436 for (i = 0; i < re->nrepl; i++)
438 if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i - 1]]))
441 "\nWARNING: The reference pressures decrease with increasing "
444 "\nWARNING: The reference pressures decrease with increasing "
450 if (replExParams.randomSeed == -1)
454 re->seed = static_cast<int>(gmx::makeRandomSeed());
460 gmx_sumi_sim(1, &(re->seed), ms);
464 re->seed = replExParams.randomSeed;
466 fprintf(fplog, "\nReplica exchange interval: %d\n", re->nst);
467 fprintf(fplog, "\nReplica random seed: %d\n", re->seed);
472 snew(re->prob_sum, re->nrepl);
473 snew(re->nexchange, re->nrepl);
474 snew(re->nmoves, re->nrepl);
475 for (i = 0; i < re->nrepl; i++)
477 snew(re->nmoves[i], re->nrepl);
479 fprintf(fplog, "Replica exchange information below: ex and x = exchange, pr = probability\n");
481 /* generate space for the helper functions so we don't have to snew each time */
483 snew(re->destinations, re->nrepl);
484 snew(re->incycle, re->nrepl);
485 snew(re->tmpswap, re->nrepl);
486 snew(re->cyclic, re->nrepl);
487 snew(re->order, re->nrepl);
488 for (i = 0; i < re->nrepl; i++)
490 snew(re->cyclic[i], re->nrepl + 1);
491 snew(re->order[i], re->nrepl);
493 /* allocate space for the functions storing the data for the replicas */
494 /* not all of these arrays needed in all cases, but they don't take
495 up much space, since the max size is nrepl**2 */
496 snew(re->prob, re->nrepl);
497 snew(re->bEx, re->nrepl);
498 snew(re->beta, re->nrepl);
499 snew(re->Vol, re->nrepl);
500 snew(re->Epot, re->nrepl);
501 snew(re->de, re->nrepl);
502 for (i = 0; i < re->nrepl; i++)
504 snew(re->de[i], re->nrepl);
506 re->nex = replExParams.numExchanges;
510 static void exchange_reals(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, real* v, int n)
520 MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
521 buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
522 ms->mastersComm_,MPI_STATUS_IGNORE);
527 MPI_Isend(v, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, &mpi_req);
528 MPI_Recv(buf, n * sizeof(real), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, MPI_STATUS_IGNORE);
529 MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
532 for (i = 0; i < n; i++)
541 static void exchange_doubles(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, double* v, int n)
551 MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
552 buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
553 ms->mastersComm_,MPI_STATUS_IGNORE);
558 MPI_Isend(v, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, &mpi_req);
559 MPI_Recv(buf, n * sizeof(double), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, MPI_STATUS_IGNORE);
560 MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
563 for (i = 0; i < n; i++)
571 static void exchange_rvecs(const gmx_multisim_t gmx_unused* ms, int gmx_unused b, rvec* v, int n)
581 MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
582 buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
583 ms->mastersComm_,MPI_STATUS_IGNORE);
588 MPI_Isend(v[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, &mpi_req);
589 MPI_Recv(buf[0], n * sizeof(rvec), MPI_BYTE, MSRANK(ms, b), 0, ms->mastersComm_, MPI_STATUS_IGNORE);
590 MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);
593 for (i = 0; i < n; i++)
595 copy_rvec(buf[i], v[i]);
601 static void exchange_state(const gmx_multisim_t* ms, int b, t_state* state)
603 /* When t_state changes, this code should be updated. */
605 ngtc = state->ngtc * state->nhchainlength;
606 nnhpres = state->nnhpres * state->nhchainlength;
607 exchange_rvecs(ms, b, state->box, DIM);
608 exchange_rvecs(ms, b, state->box_rel, DIM);
609 exchange_rvecs(ms, b, state->boxv, DIM);
610 exchange_reals(ms, b, &(state->veta), 1);
611 exchange_reals(ms, b, &(state->vol0), 1);
612 exchange_rvecs(ms, b, state->svir_prev, DIM);
613 exchange_rvecs(ms, b, state->fvir_prev, DIM);
614 exchange_rvecs(ms, b, state->pres_prev, DIM);
615 exchange_doubles(ms, b, state->nosehoover_xi.data(), ngtc);
616 exchange_doubles(ms, b, state->nosehoover_vxi.data(), ngtc);
617 exchange_doubles(ms, b, state->nhpres_xi.data(), nnhpres);
618 exchange_doubles(ms, b, state->nhpres_vxi.data(), nnhpres);
619 exchange_doubles(ms, b, state->therm_integral.data(), state->ngtc);
620 exchange_doubles(ms, b, &state->baros_integral, 1);
621 exchange_rvecs(ms, b, state->x.rvec_array(), state->natoms);
622 exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms);
625 static void copy_state_serial(const t_state* src, t_state* dest)
629 /* Currently the local state is always a pointer to the global
630 * in serial, so we should never end up here.
631 * TODO: Implement a (trivial) t_state copy once converted to C++.
633 GMX_RELEASE_ASSERT(false, "State copying is currently not implemented in replica exchange");
637 static void scale_velocities(gmx::ArrayRef<gmx::RVec> velocities, real fac)
639 for (auto& v : velocities)
645 static void print_transition_matrix(FILE* fplog, int n, int** nmoves, const int* nattempt)
650 ntot = nattempt[0] + nattempt[1];
651 fprintf(fplog, "\n");
652 fprintf(fplog, "Repl");
653 for (i = 0; i < n; i++)
655 fprintf(fplog, " "); /* put the title closer to the center */
657 fprintf(fplog, "Empirical Transition Matrix\n");
659 fprintf(fplog, "Repl");
660 for (i = 0; i < n; i++)
662 fprintf(fplog, "%8d", (i + 1));
664 fprintf(fplog, "\n");
666 for (i = 0; i < n; i++)
668 fprintf(fplog, "Repl");
669 for (j = 0; j < n; j++)
672 if (nmoves[i][j] > 0)
674 Tprint = nmoves[i][j] / (2.0 * ntot);
676 fprintf(fplog, "%8.4f", Tprint);
678 fprintf(fplog, "%3d\n", i);
682 static void print_ind(FILE* fplog, const char* leg, int n, int* ind, const gmx_bool* bEx)
686 fprintf(fplog, "Repl %2s %2d", leg, ind[0]);
687 for (i = 1; i < n; i++)
689 fprintf(fplog, " %c %2d", (bEx != nullptr && bEx[i]) ? 'x' : ' ', ind[i]);
691 fprintf(fplog, "\n");
694 static void print_allswitchind(FILE* fplog, int n, int* pind, int* allswaps, int* tmpswap)
698 for (i = 0; i < n; i++)
700 tmpswap[i] = allswaps[i];
702 for (i = 0; i < n; i++)
704 allswaps[i] = tmpswap[pind[i]];
707 fprintf(fplog, "\nAccepted Exchanges: ");
708 for (i = 0; i < n; i++)
710 fprintf(fplog, "%d ", pind[i]);
712 fprintf(fplog, "\n");
714 /* the "Order After Exchange" is the state label corresponding to the configuration that
715 started in state listed in order, i.e.
720 configuration starting in simulation 3 is now in simulation 0,
721 configuration starting in simulation 0 is now in simulation 1,
722 configuration starting in simulation 1 is now in simulation 2,
723 configuration starting in simulation 2 is now in simulation 3
725 fprintf(fplog, "Order After Exchange: ");
726 for (i = 0; i < n; i++)
728 fprintf(fplog, "%d ", allswaps[i]);
730 fprintf(fplog, "\n\n");
733 static void print_prob(FILE* fplog, const char* leg, int n, real* prob)
738 fprintf(fplog, "Repl %2s ", leg);
739 for (i = 1; i < n; i++)
743 sprintf(buf, "%4.2f", prob[i]);
744 fprintf(fplog, " %3s", buf[0] == '1' ? "1.0" : buf + 1);
751 fprintf(fplog, "\n");
754 static void print_count(FILE* fplog, const char* leg, int n, int* count)
758 fprintf(fplog, "Repl %2s ", leg);
759 for (i = 1; i < n; i++)
761 fprintf(fplog, " %4d", count[i]);
763 fprintf(fplog, "\n");
766 static real calc_delta(FILE* fplog, gmx_bool bPrint, struct gmx_repl_ex* re, int a, int b, int ap, int bp)
769 real ediff, dpV, delta = 0;
770 real* Epot = re->Epot;
773 real* beta = re->beta;
775 /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce
776 to the non permuted case */
780 case ReplicaExchangeType::Temperature:
782 * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
784 ediff = Epot[b] - Epot[a];
785 delta = -(beta[bp] - beta[ap]) * ediff;
787 case ReplicaExchangeType::Lambda:
788 /* two cases: when we are permuted, and not. */
790 ediff = E_new - E_old
791 = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
792 = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
793 = de[b][a] + de[a][b] */
796 ediff = E_new - E_old
797 = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
798 = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
799 = [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)]
800 = [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)]
801 = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]) */
802 /* but, in the current code implementation, we flip configurations, not indices . . .
803 So let's examine that.
804 = [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)]
805 = [H_b(x_ap) - H_a(x_ap)] + [H_a(x_bp) - H_b(x_pb)]
806 = (de[b][ap] - de[a][ap]) + (de[a][bp] - de[b][bp]
807 So, if we exchange b<=> bp and a<=> ap, we return to the same result.
808 So the simple solution is to flip the
809 position of perturbed and original indices in the tests.
812 ediff = (de[bp][a] - de[ap][a]) + (de[ap][b] - de[bp][b]);
813 delta = ediff * beta[a]; /* assume all same temperature in this case */
815 case ReplicaExchangeType::TemperatureLambda:
817 /* delta = reduced E_new - reduced E_old
818 = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
819 = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
820 = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
821 [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
822 = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
823 beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
824 = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
825 /* delta = beta[b]*de[b][a] + beta[a]*de[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
826 /* permuted (big breath!) */
827 /* delta = reduced E_new - reduced E_old
828 = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
829 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
830 = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
831 - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
832 - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
833 = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
834 [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
835 + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
836 = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
837 [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
838 + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b))
839 = ([beta_bp de[bp][a] - beta_ap de[ap][a]) + beta_ap de[ap][b] - beta_bp de[bp][b])
840 + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */
841 delta = beta[bp] * (de[bp][a] - de[bp][b]) + beta[ap] * (de[ap][b] - de[ap][a])
842 - (beta[bp] - beta[ap]) * (Epot[b] - Epot[a]);
844 default: gmx_incons("Unknown replica exchange quantity");
848 fprintf(fplog, "Repl %d <-> %d dE_term = %10.3e (kT)\n", a, b, delta);
852 /* revist the calculation for 5.0. Might be some improvements. */
853 dpV = (beta[ap] * re->pres[ap] - beta[bp] * re->pres[bp]) * (Vol[b] - Vol[a]) / gmx::c_presfac;
856 fprintf(fplog, " dpV = %10.3e d = %10.3e\n", dpV, delta + dpV);
863 static void test_for_replica_exchange(FILE* fplog,
864 const gmx_multisim_t* ms,
865 struct gmx_repl_ex* re,
866 const gmx_enerdata_t* enerd,
871 int m, i, j, a, b, ap, bp, i0, i1, tmp;
873 gmx_bool bPrint, bMultiEx;
874 gmx_bool* bEx = re->bEx;
875 real* prob = re->prob;
876 int* pind = re->destinations; /* permuted index */
877 gmx_bool bEpot = FALSE;
878 gmx_bool bDLambda = FALSE;
879 gmx_bool bVol = FALSE;
880 gmx::ThreeFry2x64<64> rng(re->seed, gmx::RandomDomain::ReplicaExchange);
881 gmx::UniformRealDistribution<real> uniformRealDist;
882 gmx::UniformIntDistribution<int> uniformNreplDist(0, re->nrepl - 1);
884 bMultiEx = (re->nex > 1); /* multiple exchanges at each state */
885 fprintf(fplog, "Replica exchange at step %" PRId64 " time %.5f\n", step, time);
889 for (i = 0; i < re->nrepl; i++)
894 re->Vol[re->repl] = vol;
896 if ((re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda))
898 for (i = 0; i < re->nrepl; i++)
903 re->Epot[re->repl] = enerd->term[F_EPOT];
904 /* temperatures of different states*/
905 for (i = 0; i < re->nrepl; i++)
907 re->beta[i] = 1.0 / (re->q[ReplicaExchangeType::Temperature][i] * gmx::c_boltz);
912 for (i = 0; i < re->nrepl; i++)
914 re->beta[i] = 1.0 / (re->temp * gmx::c_boltz); /* we have a single temperature */
917 if (re->type == ReplicaExchangeType::Lambda || re->type == ReplicaExchangeType::TemperatureLambda)
920 /* lambda differences. */
921 /* de[i][j] is the energy of the jth simulation in the ith Hamiltonian
922 minus the energy of the jth simulation in the jth Hamiltonian */
923 for (i = 0; i < re->nrepl; i++)
925 for (j = 0; j < re->nrepl; j++)
930 for (i = 0; i < re->nrepl; i++)
932 re->de[i][re->repl] =
933 enerd->foreignLambdaTerms.deltaH(re->q[ReplicaExchangeType::Lambda][i]);
937 /* now actually do the communication */
940 gmx_sum_sim(re->nrepl, re->Vol, ms);
944 gmx_sum_sim(re->nrepl, re->Epot, ms);
948 for (i = 0; i < re->nrepl; i++)
950 gmx_sum_sim(re->nrepl, re->de[i], ms);
954 /* make a duplicate set of indices for shuffling */
955 for (i = 0; i < re->nrepl; i++)
957 pind[i] = re->ind[i];
960 rng.restart(step, 0);
964 /* multiple random switch exchange */
968 for (i = 0; i < re->nex + nself; i++)
970 // For now this is superfluous, but just in case we ever add more
971 // calls in different branches it is safer to always reset the distribution.
972 uniformNreplDist.reset();
974 /* randomly select a pair */
975 /* in theory, could reduce this by identifying only which switches had a nonneglibible
976 probability of occurring (log p > -100) and only operate on those switches */
977 /* find out which state it is from, and what label that state currently has. Likely
978 more work that useful. */
979 i0 = uniformNreplDist(rng);
980 i1 = uniformNreplDist(rng);
984 continue; /* self-exchange, back up and do it again */
987 a = re->ind[i0]; /* what are the indices of these states? */
992 bPrint = FALSE; /* too noisy */
993 /* calculate the energy difference */
994 /* if the code changes to flip the STATES, rather than the configurations,
995 use the commented version of the code */
996 /* delta = calc_delta(fplog,bPrint,re,a,b,ap,bp); */
997 delta = calc_delta(fplog, bPrint, re, ap, bp, a, b);
999 /* we actually only use the first space in the prob and bEx array,
1000 since there are actually many switches between pairs. */
1010 if (delta > c_probabilityCutoff)
1016 prob[0] = exp(-delta);
1018 // roll a number to determine if accepted. For now it is superfluous to
1019 // reset, but just in case we ever add more calls in different branches
1020 // it is safer to always reset the distribution.
1021 uniformRealDist.reset();
1022 bEx[0] = uniformRealDist(rng) < prob[0];
1024 re->prob_sum[0] += prob[0];
1028 /* swap the states */
1030 pind[i0] = pind[i1];
1034 re->nattempt[0]++; /* keep track of total permutation trials here */
1035 print_allswitchind(fplog, re->nrepl, pind, re->allswaps, re->tmpswap);
1039 /* standard nearest neighbor replica exchange */
1041 m = (step / re->nst) % 2;
1042 for (i = 1; i < re->nrepl; i++)
1047 bPrint = (re->repl == a || re->repl == b);
1050 delta = calc_delta(fplog, bPrint, re, a, b, a, b);
1059 if (delta > c_probabilityCutoff)
1065 prob[i] = exp(-delta);
1067 // roll a number to determine if accepted. For now it is superfluous to
1068 // reset, but just in case we ever add more calls in different branches
1069 // it is safer to always reset the distribution.
1070 uniformRealDist.reset();
1071 bEx[i] = uniformRealDist(rng) < prob[i];
1073 re->prob_sum[i] += prob[i];
1077 /* swap these two */
1079 pind[i - 1] = pind[i];
1081 re->nexchange[i]++; /* statistics for back compatibility */
1090 /* print some statistics */
1091 print_ind(fplog, "ex", re->nrepl, re->ind, bEx);
1092 print_prob(fplog, "pr", re->nrepl, prob);
1093 fprintf(fplog, "\n");
1097 /* record which moves were made and accepted */
1098 for (i = 0; i < re->nrepl; i++)
1100 re->nmoves[re->ind[i]][pind[i]] += 1;
1101 re->nmoves[pind[i]][re->ind[i]] += 1;
1103 fflush(fplog); /* make sure we can see what the last exchange was */
1106 static void cyclic_decomposition(const int* destinations, int** cyclic, gmx_bool* incycle, const int nrepl, int* nswap)
1111 for (i = 0; i < nrepl; i++)
1115 for (i = 0; i < nrepl; i++) /* one cycle for each replica */
1126 for (j = 0; j < nrepl; j++) /* potentially all cycles are part, but we will break first */
1128 p = destinations[p]; /* start permuting */
1136 break; /* we've reached the original element, the cycle is complete, and we marked the end. */
1140 cyclic[i][c] = p; /* each permutation gives a new member of the cycle */
1146 *nswap = maxlen - 1;
1150 for (i = 0; i < nrepl; i++)
1152 fprintf(debug, "Cycle %d:", i);
1153 for (j = 0; j < nrepl; j++)
1155 if (cyclic[i][j] < 0)
1159 fprintf(debug, "%2d", cyclic[i][j]);
1161 fprintf(debug, "\n");
1167 static void compute_exchange_order(int** cyclic, int** order, const int nrepl, const int maxswap)
1171 for (j = 0; j < maxswap; j++)
1173 for (i = 0; i < nrepl; i++)
1175 if (cyclic[i][j + 1] >= 0)
1177 order[cyclic[i][j + 1]][j] = cyclic[i][j];
1178 order[cyclic[i][j]][j] = cyclic[i][j + 1];
1181 for (i = 0; i < nrepl; i++)
1183 if (order[i][j] < 0)
1185 order[i][j] = i; /* if it's not exchanging, it should stay this round*/
1192 fprintf(debug, "Replica Exchange Order\n");
1193 for (i = 0; i < nrepl; i++)
1195 fprintf(debug, "Replica %d:", i);
1196 for (j = 0; j < maxswap; j++)
1198 if (order[i][j] < 0)
1202 fprintf(debug, "%2d", order[i][j]);
1204 fprintf(debug, "\n");
1210 static void prepare_to_do_exchange(struct gmx_repl_ex* re, const int replica_id, int* maxswap, gmx_bool* bThisReplicaExchanged)
1213 /* Hold the cyclic decomposition of the (multiple) replica
1215 gmx_bool bAnyReplicaExchanged = FALSE;
1216 *bThisReplicaExchanged = FALSE;
1218 for (i = 0; i < re->nrepl; i++)
1220 if (re->destinations[i] != re->ind[i])
1222 /* only mark as exchanged if the index has been shuffled */
1223 bAnyReplicaExchanged = TRUE;
1227 if (bAnyReplicaExchanged)
1229 /* reinitialize the placeholder arrays */
1230 for (i = 0; i < re->nrepl; i++)
1232 for (j = 0; j < re->nrepl; j++)
1234 re->cyclic[i][j] = -1;
1235 re->order[i][j] = -1;
1239 /* Identify the cyclic decomposition of the permutation (very
1240 * fast if neighbor replica exchange). */
1241 cyclic_decomposition(re->destinations, re->cyclic, re->incycle, re->nrepl, maxswap);
1243 /* Now translate the decomposition into a replica exchange
1244 * order at each step. */
1245 compute_exchange_order(re->cyclic, re->order, re->nrepl, *maxswap);
1247 /* Did this replica do any exchange at any point? */
1248 for (j = 0; j < *maxswap; j++)
1250 if (replica_id != re->order[replica_id][j])
1252 *bThisReplicaExchanged = TRUE;
1259 gmx_bool replica_exchange(FILE* fplog,
1260 const t_commrec* cr,
1261 const gmx_multisim_t* ms,
1262 struct gmx_repl_ex* re,
1264 const gmx_enerdata_t* enerd,
1265 t_state* state_local,
1271 int exchange_partner;
1273 /* Number of rounds of exchanges needed to deal with any multiple
1275 /* Where each replica ends up after the exchange attempt(s). */
1276 /* The order in which multiple exchanges will occur. */
1277 gmx_bool bThisReplicaExchanged = FALSE;
1281 replica_id = re->repl;
1282 test_for_replica_exchange(fplog, ms, re, enerd, det(state_local->box), step, time);
1283 prepare_to_do_exchange(re, replica_id, &maxswap, &bThisReplicaExchanged);
1285 /* Do intra-simulation broadcast so all processors belonging to
1286 * each simulation know whether they need to participate in
1287 * collecting the state. Otherwise, they might as well get on with
1288 * the next thing to do. */
1289 if (haveDDAtomOrdering(*cr))
1292 MPI_Bcast(&bThisReplicaExchanged, sizeof(gmx_bool), MPI_BYTE, MASTERRANK(cr), cr->mpi_comm_mygroup);
1296 if (bThisReplicaExchanged)
1298 /* Exchange the states */
1299 /* Collect the global state on the master node */
1300 if (haveDDAtomOrdering(*cr))
1302 dd_collect_state(cr->dd, state_local, state);
1306 copy_state_serial(state_local, state);
1311 /* There will be only one swap cycle with standard replica
1312 * exchange, but there may be multiple swap cycles if we
1313 * allow multiple swaps. */
1315 for (j = 0; j < maxswap; j++)
1317 exchange_partner = re->order[replica_id][j];
1319 if (exchange_partner != replica_id)
1321 /* Exchange the global states between the master nodes */
1324 fprintf(debug, "Exchanging %d with %d\n", replica_id, exchange_partner);
1326 exchange_state(ms, exchange_partner, state);
1329 /* For temperature-type replica exchange, we need to scale
1330 * the velocities. */
1331 if (re->type == ReplicaExchangeType::Temperature || re->type == ReplicaExchangeType::TemperatureLambda)
1335 std::sqrt(re->q[ReplicaExchangeType::Temperature][replica_id]
1336 / re->q[ReplicaExchangeType::Temperature][re->destinations[replica_id]]));
1340 /* With domain decomposition the global state is distributed later */
1341 if (!haveDDAtomOrdering(*cr))
1343 /* Copy the global state to the local state data structure */
1344 copy_state_serial(state, state_local);
1348 return bThisReplicaExchanged;
1351 void print_replica_exchange_statistics(FILE* fplog, struct gmx_repl_ex* re)
1355 fprintf(fplog, "\nReplica exchange statistics\n");
1360 "Repl %d attempts, %d odd, %d even\n",
1361 re->nattempt[0] + re->nattempt[1],
1365 fprintf(fplog, "Repl average probabilities:\n");
1366 for (i = 1; i < re->nrepl; i++)
1368 if (re->nattempt[i % 2] == 0)
1374 re->prob[i] = re->prob_sum[i] / re->nattempt[i % 2];
1377 print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1378 print_prob(fplog, "", re->nrepl, re->prob);
1380 fprintf(fplog, "Repl number of exchanges:\n");
1381 print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1382 print_count(fplog, "", re->nrepl, re->nexchange);
1384 fprintf(fplog, "Repl average number of exchanges:\n");
1385 for (i = 1; i < re->nrepl; i++)
1387 if (re->nattempt[i % 2] == 0)
1393 re->prob[i] = (static_cast<real>(re->nexchange[i])) / re->nattempt[i % 2];
1396 print_ind(fplog, "", re->nrepl, re->ind, nullptr);
1397 print_prob(fplog, "", re->nrepl, re->prob);
1399 fprintf(fplog, "\n");
1401 /* print the transition matrix */
1402 print_transition_matrix(fplog, re->nrepl, re->nmoves, re->nattempt);