Improving the replica-exchange molecular-dynamics method for efficient sampling in the temperature space

Replica-exchange molecular dynamics (REMD) is a popular sampling method in the molecular simulation. By frequently exchanging the replicas at different temperatures, the molecule can jump out of the minima and sample efficiently in the conformational space. Although REMD has been shown to be practical in a lot of applications, it does have a critical limitation. All the replicas at all the temperatures must be simulated for a period between the replica-exchange steps. This may be problematic for the reaction with high free energy barriers. In that case, too many replicas are required in the simulation. To reduce the calculation quantity and improve its performance, in this paper we propose a modified REMD method. During the simulation, each replica at each temperature can stay in either the active or inactive state and only switch between the states at the exchange step. In the active state, the replica moves freely in the canonical ensemble by the normal molecular dynamics, and in the inactive state, the replica is frozen temporarily until the next exchange step. The number of the replicas in the active states (active replicas) depends on the number of CPUs in the computer. Using the additional inactive replicas, one can perform an REMD simulation in a wider temperature space. The practical applications show that the modified REMD method is reliable. With the same number of active replicas, this REMD method can produce a more reasonable free energy surface around the free energy minima than the standard REMD method.

Figure 1

(a) Canonical distribution function in the energy space of the quantum harmonic oscillator at $T=5.0$ K and $T=1000$ K. The data from the REM simulation (ten active replicas) and modified REM simulation (three active replicas) are shown by the symbols “×” and “○”, respectively. The lines in the figure represent the theoretical values [Eq. (5)]. (b) Average energies of the harmonic oscillator at all temperatures. (c) Heat capacities of the harmonic oscillator at all temperatures. Theoretical values in (b) and (c) are calculated by Eq. (7). Energies and heat capacities are in unit of $k_B$ (J) and (J/K) respectively.

Figure 2

Canonical probability distributions of ALA dipeptide in the energy space at 300 K (left), 624 K (middle), and 1300 K (right) in the three simulations. The lines, dash-dot lines, and dashed lines represent the data from the first, the second, and the third simulation, respectively.

Figure 3

Free energy surface of ALA dipeptide in the 2D collective variable space from the seven-replica REMD simulation (a), three-replica REMD simulation (b), and three-replica modified REMD simulation (c). From the information, the free energy profiles on the path connecting the minimum ${C7}_{\mathrm{eq}}$ and the minimum ${C7}_{\mathrm{ax}}$ [marked in (a)] are shown in (d) by line, dash-dot line, and dashed line, respectively.

Figure 4

(a) Convergence of the free energy surface in the first standard REMD simulation (line), the second standard REMD simulation (dash-dot line), and the third modified REMD simulation (dashed line). The data are updated every 1 ns by Eq. (8). (b) Average root-mean-square errors of the free energy profiles on the transition path in the three simulations (see location of the path in Fig. 3 and the profiles in Fig. 3]. The root-mean-square error of each simulation is calculated from five independent 200-ns trajectories.

Figure 5

Free energy surface of trpzip2 in the standard REMD simulation (a) and the modified REMD simulation (b). The two collective variables are backbone RMSD and average radius of gyration of the aromatic pairs [Eq. (9)].