In “Replica-exchange-with-tunneling for fast exploration of protein landscapes” [F. Yaşar *et al.*, J. Chem. Phys. **143**, 224102 (2015)], a novel sampling algorithm called “Replica Exchange with Tunneling” was proposed. However, due to its violation of the detailed balance, the algorithm fails to sample from the correct canonical ensemble.

In the Metropolis’ algorithm,^{1} a new coordinate *x*_{k+1} is proposed from the previous coordinate *x _{k}* probabilistically, according to the distribution $p x k + 1 \u2223 x k $ that satisfies $p x k + 1 \u2223 x k =p x k \u2223 x k + 1 .$ At each step, the newly proposed coordinate

*x*

_{k+1}is accepted with probability

If the new coordinate is not accepted, then *x*_{k+1} ≔ *x _{k}*. On updating the coordinate repeatedly,

*x*is sampled from the distribution $\pi x k $, if the system is ergodic. The probability distribution corresponding to the canonical ensemble is often used, i.e., $\pi x = Z \u2212 1 exp \u2212 \beta E x $, where

_{k}*Z*is the partition function, $\beta = k B T \u2212 1 $,

*k*is the Boltzmann factor,

_{B}*T*is the temperature, and $E x $ is the potential energy or the total energy of the system.

The Hybrid Monte Carlo (HMC) method algorithm proposed by Duane *et al.*^{2} uses a deterministic simulation to generate proposals. Given a starting coordinate *x _{k}*, a velocity

*v*is generated according to $ P kin v k \u221dexp \u2212 \beta K v k $, where

_{k}*K*is the kinetic energy function. A new coordinate

*x*

_{k+1}is then generated by first performing

*τ*steps of a simulation from $ x 0 , v 0 = x k , v k $, namely, $ x 0 , v 0 \u21a6 \tau x \tau , v \tau $ and then taking $ x k + 1 = x \tau $ as the proposed coordinate. The newly proposed coordinate is accepted with probability

using the probability distribution depending on the total energy $ \pi tot x , v \u221dexp \u2212 \beta E x + K v $. Since the proposal is deterministic, the probability distribution of the proposal is written as

where $\delta \u22c5 $ is the Dirac delta function. It can be shown that $ p HMC x k + 1 , v k + 1 \u2223 x k , v k = p HMC x k , v k \u2223 x k + 1 , v k + 1 $, when three conditions hold: (a) the system is symmetric with respect to the sign inversion of velocities; (b) the simulation is time-reversible $ x \tau , \u2212 v \tau \u21a6 \tau x 0 , \u2212 v 0 $; (c) the phase space volume is preserved before and after the deterministic simulation $d x 0 d v 0 =d x \tau d v \tau $,

Equations (6) and (8) hold by condition (a); Eq. (7) holds by condition (b); and Eq. (5) holds by condition (c), respectively. When all three conditions are met, the distribution of *x _{k}* asymptotically converges to $\pi x k $.

^{2}An example of algorithms that satisfy all three conditions is the leapfrog algorithm, which was used in HMC. Both in HMC and in its many variants,

^{3}proposals are generated by carefully designed algorithms to ensure all three conditions.

In “Replica-exchange-with-tunneling for fast exploration of protein landscapes,”^{4} the algorithm RET is constructed as a combination of the HMC and the replica exchange molecular dynamics (REMD),^{5} realized by (d) a microcanonical simulation, (e) an exchange and scaling process, and (f) the other microcanonical simulation. Although the leapfrog algorithm preserves the phase-space volume for microcanonical simulation parts ((d) and (f)), the exchange and scaling part (e) does not preserve the phase space volume. Indeed, the velocity variables are scaled at the procedure (e) as in Eq. (3) of the paper, resulting in the (infinitesimal) phase space volume change of

where *N* is the number of particles in the system. The right hand side of Eq. (9) does not vanish, and thus it violates the phase-space volume preservation. It therefore does not satisfy the detailed balance condition, meaning that the resulting ensemble is not sampled from the canonical distribution. It is notable that in the original REMD,^{5} the phase-space volume does not change because scaling factors cancel each other.

To numerically illustrate the problem of the unpreserved phase space volume, I performed the RET simulation of a simple system of an ethanol molecule *in vacuo* modeled with the optimized potentials for liquid simulations—all-atom (OPLS-AA)^{6} force field, using the program in the supplementary material of Ref. 4. The simulations were run both in NVT (600 K) and in RET with temperatures set at 300 K, 340 K, 390 K, 440 K, 500 K, and 600 K. The velocity rescaling algorithm^{7} was used as a thermostat in all simulations. The time step was set to 0.5 fs, and the RET interval was set to 20 steps to demonstrate the problem discussed. The relaxation time in RET was set to 5 steps. The kinetic energies corresponding to simulations in 600 K were recorded every 10 fs. In RET, energies were recorded when microcanonical simulations start. The simulation was run to 100 ps for equilibrium and a further 5 ns for production, for each method. Kinetic energies of this system are expected to converge into the *χ*^{2}-distribution, independent of the methods used. Figure 1 presents the distributions of kinetic energies. The distributions clearly show that the RET exhibited a different kinetic energy distribution than NVT did. Further, means and standard errors of instantaneous temperatures derived from kinetic energies were 600.1 ± 0.3 K in NVT and 578.0 ± 0.2 K in the RET, respectively, indicating that the trajectory in the RET is not sampled from the canonical ensemble at the target temperature. The phase space volume is thus not negligible, even in practice. Additionally, in Fig. 4 of Ref. 4, the average fraction of the configurations in the RET systematically deviated from the REMD particularly in low-temperature regions, which may also be caused from the violation of detailed balance condition.

I thank Dr. Yasuhiro Matsunaga for discussions.