Response to “Comment on ‘Replica-exchange-with-tunneling for fast exploration of protein landscapes”’ [J. Chem. Phys. 145, 057101 (2016)] Free

In Replica Exchange with Tunneling (RET),¹ at temperature T₁, a configuration A is replaced by $Bˆ$⁠, and B by $Aˆ$ at temperature T₂, by a set of deterministic moves,

where A = (q_A, v_A) denotes a state characterized by coordinates q_A of all its atoms and the associated velocities v_A. The evolution of A (B) to A′ (B′) and B″ (A″) to $Bˆ$ (⁠$Aˆ$⁠) is by short microcanonical molecular dynamics. The crossing arrows mark the exchange step which involves a rescaling of velocities (and the associated kinetic energies),

such that

with $A″=(qA′,vA″)$ and $B″=(qB′,vB″)$⁠.

In his comment,² Dr. Sakuraba raises the interesting point that the exchange move does not conserve phase space volume (Eq. (9) of the comment). This is an important observation as it implies that the exchange move is asymmetric. This can be seen by integrating out the velocity degrees of freedom leading to an expression for the probability to find configurations with potential energy $Epot(qA′)$ in a microcanonical ensemble with total energy E₁,

with N the number of particles and $Ω(Epot(qA′))$ the density of states with potential energy E_pot. As the probability for finding a state with coordinates $qA′$ is given by $1/Ω(Epot(qA′))$⁠, and for finding a state with coordinates $qB′$ by $1/Ω(Epot(qB′))$⁠, it follows that in a microcanonical ensemble of energy E₁, the probability to be in a configuration with the coordinates $qA′$ and picking a configuration with coordinates $qB′$ is given by

A similar expression holds in the microcanonical ensemble of energy E₂ where the probability to be in a configuration with the coordinates $qB′$ and picking a configuration with coordinates $qA′$ is given by

Hence, while the exchange move is rejection-free (all states have equal weight in a microcanonical ensemble), the proposal selection probability is asymmetric and given by

which is Equation (9) of the comment.

In RET, the exchange move is combined with short microcanonical molecular dynamics trajectories in which the total energy stays constant but the distribution of kinetic and potential energy quickly changes according to the equipartition theorem. For large systems and sufficiently long runs this will lead to final states where the velocities $vˆB$ (⁠$vˆA$⁠) follow a Maxwell-Boltzmann distribution corresponding to temperature T₁ (T₂). However, Eq. (6) implies that the initial velocities $vB″$ (⁠$vA″$⁠) of these trajectories are not Maxwell-Boltzmann distributed. It follows that as the exchange move these microcanonical molecular dynamic simulations do not conserve phase space volume. The net effect of exchange move and microcanonical molecular dynamics trajectories is that at temperature T₁, starting from a configuration A with potential energy E_pot(q_A) and total energy E₁, a configuration $Bˆ$ is proposed with potential energy $Epot(qˆB)$ and the same total energy E₁; and at temperature T₂, a configuration B(E_pot(q_B), E₂) is compared with a configuration $Aˆ(Epot(qˆA),E2)$⁠. Hence, in a microcanonical picture, the RET move is accepted with a probability one, but the ratio of the proposal probabilities differs from one and—following the same argument that led from Eqs. (4) to (6)—is given by

The Metropolis-Hastings that guarantees convergence to the correct distribution evaluates the product of proposal and acceptance probabilities, leading to the following criterion for acceptance or rejection of the RET move:

This equation can be simplified using that in Eq. (3), and both functions on the right-hand side grow strongly with their arguments. Hence, for systems with large number of particles N, P_NVE(E_pot, E) is a sharply peaked function, and a saddle point expansion will lead to

with the inverse microcanonical temperature β_E = 1/k_BT_E = dlnΩ(E)/dE and $Eˆpot$ the most probable potential energy. Hence, for sufficiently large systems and long enough microcanonical segments, the RET move can be accepted with probability,

which is Eq. (4) of our paper (Ref. 1), and implies that in this limit phase space volume is conserved during the RET move.

These conditions for the validity of the RET move are not given in the example that Dr. Sakuraba presents in the comment. His system is a single ethanol molecule in vacuo. As this is a non-interacting system, the velocities at the end of the microcanonical segment are unchanged from the ones obtained by the rescaling of Eq. (1), and the deviations in Figure 1 of the comment are not surprising. While we had already in our paper outlined that microcanonical segment needs to be long enough for a redistribution of energies, i.e., leading to configurations with more favorable potential energies and velocities whose distribution corresponds to the target temperatures, we had not taken into consideration the case of a too small system. We show in Figure 1 the same kind of distribution for a system of 256 ethanol molecules simulated at the same temperature and with the same protocol, except that the number of microcanonical steps is raised to 1000 and 20 000 steps between RET moves. For this larger system and longer microcanonical segment, we find good agreement between RET and regular Replica Exchange Molecular Dynamics (REMD).

We are grateful to Dr. Sakuraba for pointing out that the replica exchange with tunneling approach is valid only for sufficiently large systems and long enough microcanonical segment length. A way to test that these conditions are realized is to monitor the distribution of kinetic energies as done by us in Figure 1 of our article¹ and (as a negative example for too small system) in the comment.² 

We acknowledge support from the National Science Foundation (Research Grant No. CHE-1266256) and by Hacettepe University Scientific Research Fund under Project No. FHD.2015.6939.