Temperature REM

The Replica Exchange Method is based on multiple concurrent (parallel) canonical simulation that are allowed to occasionally exchange their configurations. For a system made of $N$ atoms, by ``configuration'' we mean a state defined by a $3N$ dimensional coordinate vector, independent of the momenta. Thus, in a replica exchange, only coordinates and not momenta are exchanged. In the standard implementation of the methodology, each replica, bearing a common interaction potential, is characterized by a given temperature and configurations between couple of replicas are tentatively exchanged at prescribed time intervals using a probabilistic criterion. The target temperature, i.e. the temperature corresponding to the thermodynamic state of interest, is usually the lowest among all replicas. In this manner, ``hot'' configurations from hot replicas, i.e. configurations where energy barrier are easily crossed, may be occasionally accepted at the target temperature. The canonical probability of a coordinate configuration $X$ for $m$-th replica is given by

$$P_m(X)=\frac{1}{Z_m}e^{-\beta_m V(X)}$$ (5.1)

where $m$ is the replica index, $\beta^{-1}=k_B T$, $V(x)$ is the potential of the system, and $Z_{m}=\int e^{-\beta_m V(X)}dX$ is the configurational partition function for $m$-th replica. Being the $M$ replicas independent, the probability distribution for a generic configuration of the $M$-fold extended system ${\bf X} = (X_1,..X_M)$ is

$$P_{\bf X}=\prod_m^MP_m(X_m)$$ (5.2)

As stated above, the global state ${\bf X}$ of the extended system may evolve in two ways: i) by evolving each replica independently (i.e. via MC or MD simulation protocols) and ii) by exchanging the configurations of two replicas. Regarding the second mechanism, we introduce the transition probability $W(X,\beta_m;X^\prime,\beta_n)$ for the exchange between the configuration $X$ of replica at $T_m$ and the configuration $X^\prime$ for the replica $T_n$. The probability for the inverse exchange is clearly given by $W(X^\prime,\beta_m;X,\beta_n)$. The detailed balance condition on the extended system for this kind of moves is given by

$$P_{\bf X}(...,X,\beta_m,X^\prime,\beta_n,...)W(X,\beta_m,X^\prime,\beta_n)=$$ (5.3)

$$P_{\bf X^{\prime}}(...,X^\prime,\beta_m;X,\beta_n,...)W(X^\prime,\beta_m;X,\beta_n)$$ (5.4)

which, using the expressions 5.2 and 5.1 for the global probability, is satisfied if the transition probability satisfies the equation

$$\frac{W(X,\beta_m,X^{\prime},\beta_n)}{W(X^{\prime},\beta_m;X,\beta_n)}=e^{-(\beta_m-\beta_n)(E(X^{\prime})-E(X))}.$$ (5.5)

The exchange of configurations of replicas obeying the detailed balance condition 5.5 can be as usual implemented by using the Metropolis algorithm

$$P_{\rm acc}=\min(1,e^{-\Delta})$$ (5.6)

with $\Delta=(\beta_m-\beta_n)(E(X^{\prime)}-E(X))$. Like in a standard MC technique, because of the detailed balance condition for the extended system, the sampling in the ${\bf X}$ multi-configuration space in REM evolves towards a global equilibrium defined by the multi-canonical probability distribution of the extended system Eq. 5.2.

In principle Eq. 5.6 refers to the probability of an exchange between any two replicas. In practice the exchanges are attempted between replicas that are contiguous in temperature. Let's see why. For any two replicas $m$, $n$, the total number of accepted exchanges between them is given by

$$N^{\rm acc}=N^{\rm acc}(\Delta E < 0 ) + N^{\rm acc}(\Delta E >0 )$$ (5.7)

where $\Delta E = (E(X^{\prime)}-E(X))$ and $N^{\rm acc}(\Delta E < 0 )$, $N^{\rm acc}(\Delta E >0 )$ are the number of accepted exchanges for which $\Delta E<0$ and $\Delta E >0$, respectively. When the extended system is at equilibrium, we clearly must have that

$$N^{\rm acc}(\Delta E < 0 )= N^{\rm acc}(\Delta E >0 )$$ (5.8)

Inserting the above equation into Eq. 5.7, we obtain

$$N^{acc} = 2 N^{acc}(\Delta E < 0)$$ (5.9)

Since, according to the prescription 5.6, the probability for accepting the move when $\Delta E<0$ is unitary, we may write that

$$\frac{N^{acc}}{N^{\rm tot}} = 2 P(\Delta E < 0)$$ (5.10)

where $N^{\rm tot}$ is the total number of attempted exchanges and $P(\Delta E < 0)$ is the cumulative probability that a $E(X^{\prime})< E(X)$.

Figure: Overlapping configurational energy distribution for two replicas. The shaded area is the acceptance probability for the configuration exchange. The overlap between the two distribution is a lower bound for the acceptance probability

Eq. 5.10 states that if the two normalized (configurational) energy distribution $P_m(E)$ of replica $m$ and $P_n(E)$ of replica $n$ are identical, then the probability for a successful exchange between the two replica is equal to the area of the overlap of the two distribution (i.e. the shaded area in Fig. 5.1). If $P_m(E)$ and $P_n(E)$ are not identical, we have in general that the overlap of the two distribution is a lower bound for the acceptance probability (the standard deviation $\delta E$ generally increases with the mean energy $\bar E$). Based on the above, and assuming that $M$, the total number of replicas, is even, one can then set up an exchange protocol periodically attempting $M/2$ simultaneous contiguous replica exchanges $m\leftrightarrow m+1$ with $m$ odd, or $M/2-1$ simultaneous contiguous replica exchanges $m\leftrightarrow m+1$ with $m$ even, accepting each of them with probability given by 5.6.

Given the above scheme, what is the optimal spacing in temperatures for enhanced sampling of the configuration space at the target temperature? First of all, the hottest temperature $T_M$, defining the full temperature range $\Delta T = T_M-T_1$ of the extended system, must be clearly selected such that $k_B T_M$ is of the order of the maximum height of the free energy barriers that must be overcome at the target temperature $T_1$. Concerning the temperature spacing, we have seen that acceptance probability for an exchange is larger, the larger is the overlap of the two energy distributions referring to the two contiguous replica, i.e. the closer are the temperatures. Of course, the closer are the temperatures and the larger is the number of replicas to be simulated, i.e. the heavier is the CPU cost of the simulation. For an optimal choice, we thus set

$$\bar {E}_{m+1} - \bar {E}_{m} = \sigma_{ E_m}$$ (5.11)

where $\bar {E}_{m}$ and $\sigma_{ E_m}$ are the mean energy and the standard deviation of energy distribution for the $m$-the replica. Assuming then that the system can be described by an ensemble of $N$ harmonic oscillators, we have that $\bar E_m =N kT_m$ and $\sigma_{E_m} = c N^{1/2} k^{1/2} T_{m}$.^5.1 Substituting these values in Eq. 5.11, we obtain the temperature spacing for optimal superposition:

$$\bar {T}_{m+1} - \bar {T}_{m} = \left ( \frac{c^2}{Nk_b} \right )^{1/2} T_m$$ (5.12)

In the parallel implementation of the temperature REM, in order to keep the communication overhead at the lowest possible level, we standardly exchange the temperatures and not the configurations. So the $m$-th slave process may explore the entire range of temperatures. When the $m$-th slave process periodically writes out the coordinates of the configuration (Typically in pdb or xyz format), one must also keep track of the current temperature (the program does this automatically) in order be able to reconstruct a posteriori the true $m$-th temperature configurational space of the $m$-th replica.

Figure: Typical REM Simulation with 8 replicas. Each process bear a particular color and the color follows the right scale, i.e. the replica index which is connected to the temperature. To reconstruct a trajectory at a given temperatures, one must combine the data form several processes.

In Fig. 5.2, we show a typical parallel REM simulation for a general system with 8 processes. In the x-axis we report the simulation time, in the left y-axis the process index and in the right y-axis the replica index which is bound to the actual temperature. Each color represents a process running in parallel with other processes with different colors. As it can be seen, on each process the temperature (i.e. the replica index) changes continuously.. So, for example the configurational sampling of the replica at the lowest temperature in the given time interval must be reconstructed combining the data for the slave processes 1,2,3,4,6 If the algorithm is working properly, (i.e. if the temperature spacing is chosen correctly and if there are no phase transition between $T_1$ and $T_M$), the temperature in each parallel process must perform a random walk in the temperature domain $[T_1,T_M]$.

Going back to equation 5.12, two important issues must be stressed: i) the temperature spacing for optimal overlap between contiguous replicas while keeping the total number of replicas not too high, is not uniform but grows with the replica temperature; ii) the temperature spacing between contiguous replicas must be decreased with increasing number of degrees of freedom. The latter is indeed a severe limitation of the standard REM technology, since, as the size of the system grows, a larger number of replicas must be employed for preserving a significant exchange acceptance probability. This is due to the inescapable fact that the energy fluctuations grow with $N^{1/2}$ while the energy grows with $N$. Moreover, in many important cases, one has to effectively samples reaction coordinates that are rather localized in the protein, like e.g. in the case of substrate-active site interactions. In the standard temperature REM, the extra heat in the hot replicas is clearly distributed among all the degrees of freedom of the system and therefore most of this heat is uselessly used for exchanging uninteresting configurations (e.g. solvent configurations).