Non-equilibrium dynamics from RPMD and CMD

Welsch, Ralph; Song, Kai; Shi, Qiang; Althorpe, Stuart C.; Miller, Thomas F.

doi:10.1063/1.4967958

We investigate the calculation of approximate non-equilibrium quantum time correlation functions (TCFs) using two popular path-integral-based molecular dynamics methods, ring-polymer molecular dynamics (RPMD) and centroid molecular dynamics (CMD). It is shown that for the cases of a sudden vertical excitation and an initial momentum impulse, both RPMD and CMD yield non-equilibrium TCFs for linear operators that are exact for high temperatures, in the t = 0 limit, and for harmonic potentials; the subset of these conditions that are preserved for non-equilibrium TCFs of non-linear operators is also discussed. Furthermore, it is shown that for these non-equilibrium initial conditions, both methods retain the connection to Matsubara dynamics that has previously been established for equilibrium initial conditions. Comparison of non-equilibrium TCFs from RPMD and CMD to Matsubara dynamics at short times reveals the orders in time to which the methods agree. Specifically, for the position-autocorrelation function associated with sudden vertical excitation, RPMD and CMD agree with Matsubara dynamics up to $O (t^{4})$ and $O (t^{1})$ ⁠, respectively; for the position-autocorrelation function associated with an initial momentum impulse, RPMD and CMD agree with Matsubara dynamics up to $O (t^{5})$ and $O (t^{2})$ ⁠, respectively. Numerical tests using model potentials for a wide range of non-equilibrium initial conditions show that RPMD and CMD yield non-equilibrium TCFs with an accuracy that is comparable to that for equilibrium TCFs. RPMD is also used to investigate excited-state proton transfer in a system-bath model, and it is compared to numerically exact calculations performed using a recently developed version of the Liouville space hierarchical equation of motion approach; again, similar accuracy is observed for non-equilibrium and equilibrium initial conditions.

I. INTRODUCTION

Quantum time correlation functions (TCFs) play an important role in the description of chemical dynamics, which has led to the development of numerous approximate methods for their calculation.^1–23 Two widely used methods based on imaginary-time path integrals are centroid molecular dynamics (CMD)^24–28 and ring-polymer molecular dynamics (RPMD).^29–32 Although both methods have known artifacts, such as the spurious-mode effect in RPMD^33–35 and the curvature problem^35,36 in CMD, they have proven effective for a vast range of chemical applications including the calculation of thermal rate constants,^30,37–55 diffusion coefficients,^31,56–61 and vibrational spectra.^{33–36,62–65}

With only a few exceptions,^66–68 RPMD and CMD have been applied for the characterization of processes with thermal equilibrium initial conditions. The aim of this work is to systematically investigate whether RPMD and CMD can also be usefully employed in the context of non-equilibrium TCFs and expectation values. It is found that nearly all important properties of RPMD and CMD are preserved when calculating non-equilibrium TCFs, as are the relationships of the two methods to Matsubara dynamics.^69,70 Furthermore, the numerical performance of the two methods for the calculation of non-equilibrium quantum TCFs is tested for a range of model systems.

The paper is organized as follows. Sec. II describes the approaches to calculate non-equilibrium TCFs using RPMD and CMD and also investigates the performance of the approaches in important limiting cases analytically. Sec. III investigates the numerical performance of the methods for one-dimensional potentials and for non-equilibrium proton transfer in a system-bath model, and conclusions are presented in Sec. IV.

II. THEORY

A. RPMD and CMD for equilibrium TCFs

Both RPMD and CMD employ the machinery of classical molecular dynamics to approximate the equilibrium Kubo-transformed quantum TCF,

{\tilde{C}}_{A B} (t) = \frac{1}{β} \int_{0}^{β} d λ \times Tr [e^{- λ \hat{H}} \hat{A} e^{- (β - λ) \hat{H}} e^{i \hat{H} t / ℏ} \hat{B} e^{- i \hat{H} t / ℏ}],

(1)

where

\hat{H} = \frac{{\hat{p}}^{2}}{2 m} + V (\hat{q}) .

(2)

For simplicity, equations are presented in one dimension and are easily generalized.

In RPMD, ${\tilde{C}}_{A B} (t)$ is approximated using

{\tilde{C}}_{A B}^{(RP)} (t) = \frac{1}{{(2 π ℏ)}^{N}} \int d 𝐪_{0} \int d 𝐩_{0} e^{- β_{N} H_{N} (𝐩_{0}, 𝐪_{0})} \times A (𝐩_{0}, 𝐪_{0}) B (𝐩_{t}, 𝐪_{t}),

(3)

where $β_{N} = \frac{β}{N}$ ⁠,

H_{N} (𝐩, 𝐪) = \sum_{i = 1}^{N} \frac{p_{i}^{2}}{2 m} + U_{N} (𝐪)

(4)

is the ring-polymer Hamiltonian,

U_{N} (𝐪) = \sum_{i = 1}^{N} [\frac{m {(q_{i} - q_{i - 1})}^{2}}{2 {(β_{N} ℏ)}^{2}} + V (q_{i})],

(5)

and q₁ = q_N. The classical dynamics associated with the Hamiltonian in Eq. (4) determines the time evolution of the ring-polymer coordinates (p_t, q_t), such that

{\dot{p}}_{i} = \frac{- m}{{(β_{N} ℏ)}^{2}} (2 q_{i} - q_{i - 1} - q_{i + 1}) - \frac{\partial V (q_{i})}{\partial q_{i}},

(6)

{\dot{q}}_{i} = \frac{p_{i}}{m} .

(7)

In CMD, ${\tilde{C}}_{A B} (t)$ is approximated using

{\tilde{C}}_{A B}^{(CMD)} (t) = \frac{1}{2 π ℏ} \int d {\bar{Q}}_{0} \int d {\bar{P}}_{0} e^{- β ({\bar{P}}_{0}^{2} / 2 m + W ({\bar{Q}}_{0}))} \times A ({\bar{Q}}_{0}, {\bar{P}}_{0}) B ({\bar{Q}}_{t}, {\bar{P}}_{t}),

(8)

where

W (\bar{Q}) = - \frac{1}{β} \ln {{(\frac{m}{2 π β_{N} ℏ^{2}})}^{(N - 1) / 2} \times \int d 𝐪 e^{- β_{N} U_{N} (𝐪)} δ (\bar{q} (𝐪) - \bar{Q})}

(9)

is the centroid potential of mean force (PMF), and $\bar{q} (𝐪)$ is the ring-polymer position centroid. The centroid position and momenta evolve according to classical dynamics subject to the centroid PMF as

\dot{\bar{P}} = - \frac{\int d 𝐪 e^{- β_{N} U_{N} (𝐪)} \frac{\partial U_{N} (𝐪)}{\partial \bar{Q}} δ (\bar{q} (𝐪) - \bar{Q})}{\int d 𝐪 e^{- β_{N} U_{N} (𝐪)} δ (\bar{q} (𝐪) - \bar{Q})},

(10)

\dot{\bar{Q}} = \frac{\bar{P}}{m} .

(11)

Both RPMD and CMD preserve the quantum Boltzmann distribution and are exact for the description of ${\tilde{C}}_{A B} (t)$ in several important limits, including harmonic potentials, high temperature, and at t = 0 for autocorrelation functions of linear operators; RPMD is likewise exact at t = 0 for autocorrelation functions of non-linear operators. Several efforts to derive RPMD and CMD have been undertaken.^69–72 In particular, it has been found that both RPMD and CMD are related to Matsubara dynamics,^69,70 which approximates ${\tilde{C}}_{A B} (t)$ by a mixture of quantum statistics and classical dynamics. Matsubara dynamics is derived from first principles and has been shown to reproduce TCFs more accurately than RPMD, CMD, or the linearized semiclassical-initial value representation approach, although its slow numerical convergence makes it less broadly applicable than these other approximate methods.^69,70

B. Application of RPMD and CMD to non-equilibrium time-correlation functions

Here, we apply RPMD and CMD to non-equilibrium TCFs of the form

{\tilde{C}}_{A B} (t) = \frac{1}{β} \int_{0}^{β} d λ \times Tr [e^{- λ {\hat{H}}^{(0)}} \hat{A} e^{- (β - λ) {\hat{H}}^{(0)}} e^{i {\hat{H}}^{(1)} t / ℏ} \hat{B} e^{- i {\hat{H}}^{(1)} t / ℏ}] .

(12)

We will consider two different cases for which ${\hat{H}}^{(0)} \neq {\hat{H}}^{(1)}$ ⁠, associated with either a sudden vertical excitation or an initial momentum impulse. Although these two cases will be described separately, it is straightforward to apply them in combination.

1. Non-equilibrium initial conditions associated with a vertical excitation

We first consider the case of a sudden vertical excitation, for which the Hamiltonians differ only in the potential energy,

{\hat{H}}^{(0)} = \frac{{\hat{p}}^{2}}{2 m} + V^{(0)} (\hat{q}),

(13)

{\hat{H}}^{(1)} = \frac{{\hat{p}}^{2}}{2 m} + V^{(1)} (\hat{q}),

(14)

and we define $Δ V (\hat{q}) = V^{(1)} (\hat{q}) - V^{(0)} (\hat{q})$ ⁠. This case is directly related to the Condon approximation, as is shown in Appendix A.

The corresponding ring-polymer Hamiltonians are

H_{N}^{(j)} (𝐩, 𝐪) = \sum_{i = 1}^{N} \frac{p_{i}^{2}}{2 m} + U_{N}^{(j)} (𝐪),

(15)

where j = 0,1,

U_{N}^{(j)} (𝐪) = \sum_{i = 1}^{N} [\frac{m {(q_{i} - q_{i - 1})}^{2}}{2 {(β_{N} ℏ)}^{2}} + V^{(j)} (q_{i})],

(16)

and the RPMD approximation is

{\tilde{C}}_{A B}^{(RP)} (t) = {(\frac{1}{2 π ℏ})}^{N} \int d 𝐪 \int d 𝐩 \times e^{- β_{N} H_{N}^{(0)} (𝐩, 𝐪)} A (𝐩, 𝐪) B (𝐩_{t}, 𝐪_{t}),

(17)

where the time-evolution of p_t and q_t is governed by $H_{N}^{(1)}$ ⁠.

The CMD approximation for this case is

{\tilde{C}}_{A B}^{(CMD)} (t) = \frac{1}{2 π ℏ} \int d {\bar{Q}}_{0} \int d {\bar{P}}_{0} e^{- β ({\bar{P}}_{0}^{2} / 2 m + W^{(0)} ({\bar{Q}}_{0}))} \times A ({\bar{Q}}_{0}, {\bar{P}}_{0}) B ({\bar{Q}}_{t}, {\bar{P}}_{t}),

(18)

with the classical dynamics of the centroid coordinates governed by the centroid mean force as

F_{1} (\bar{Q}) = - \frac{\int d 𝐪 e^{- β_{N} U_{N}^{(1)} (𝐪)} \frac{\partial U_{N}^{(1)} (𝐪)}{\partial \bar{Q}} δ (\bar{q} (𝐪) - \bar{Q})}{\int d 𝐪 e^{- β_{N} U_{N}^{(1)} (𝐪)} δ (\bar{q} (𝐪) - \bar{Q})},

(19)

where $U_{N}^{(1)} (𝐪)$ is the ring-polymer potential associated with V⁽¹⁾(q), and $W^{(0)} ({\bar{Q}}_{0})$ is the PMF from Eq. (9) associated with V⁽⁰⁾(q). This non-equilibrium implementation of CMD is similar to the protocol given in Ref. 67 and corresponds to instantaneous thermalization of the non-centroid ring-polymer modes on the potential energy surface V⁽¹⁾(q).

The RPMD and CMD approximations for non-equilibrium TCFs given in Eqs. (17) and (18) preserve many of the formal properties as for the calculation of equilibrium TCFs. Both methods correctly reduce to the classical limit and capture the high-temperature regime exactly. Furthermore, in Appendix B, we use linearization of the difference between the forward and backward Feynman paths in real time to show that RPMD and CMD exactly reproduce non-equilibrium TCFs of linear operators in harmonic potentials, regardless of the initial potential V⁽⁰⁾(q).

Also unchanged in going from equilibrium to non-equilibrium initial conditions are the formal relationships of RPMD and CMD to Matsubara dynamics. To obtain RPMD from Matsubara dynamics, the Matsubara Liouvillian $L_{M}$ is written in terms of complex phase space as^69,70

L_{M} = L_{M}^{(RP)} + i L_{M}^{(I)},

(20)

where the RPMD Liouvillians $L_{M}^{(RP)}$ and $L_{M}^{(I)}$ are given in Ref. 70 and both $L_{M}^{(RP)}$ and $L_{M}^{(I)}$ conserve the ring-polymer Hamiltonian. The RPMD equations of motion are then obtained by discarding the imaginary component $L_{M}^{(I)}$ and analytically continuing $L_{M}^{(RP)}$ into real space.⁷³ This procedure can be done independently of the initial conditions, as the Matsubara Liouvillian $L_{M}$ only involves the Hamiltonian governing the time-evolution of the system (i.e., ${\hat{H}}^{(1)}$ for the present case). To obtain CMD from Matsubara dynamics, one invokes a centroid mean-field approximation of the exact force on the centroid^69,70

- \frac{\partial {\tilde{U}}_{M} (\tilde{𝐐})}{\partial \bar{q} (\tilde{𝐐)}} = F_{0} (\bar{q} (\tilde{𝐐})) + F_{fluct.} (\tilde{𝐐}),

(21)

where the Matsubara potential ${\tilde{U}}_{M} (\tilde{𝐐})$ is defined in Appendix C, $\bar{q} (\tilde{𝐐})$ is the centroid coordinate, and the fluctuating part of the force, $F_{fluct.} (\tilde{𝐐})$ ⁠, is neglected.⁷⁰ Again, this procedure is independent of the initial distribution. Notably, RPMD leaves the phase $θ_{M}$ appearing in Matsubara dynamics invariant during the time evolution for both equilibrium and non-equilibrium initial conditions, as the action of the ring-polymer Liouvillian $L^{(RP)}$ on the Boltzmann distribution,

L^{(RP)} e^{- β_{N} H_{N}^{(0)} (𝐩_{t}, 𝐪_{t})} = e^{- β_{N} H_{N}^{(0)} (𝐩_{t}, 𝐪_{t})} \times \sum_{j = 1}^{N} \frac{p_{t, j}}{m} \frac{\partial Δ V (𝐪_{t})}{\partial q_{t, j}},

(22)

does not involve the inter-bead harmonic spring-term, which is equivalent to preservation of the phase $θ_{M}$ ⁠.⁷⁰ Similarly, CMD leaves the phase $θ_{M}$ invariant during the time evolution for both equilibrium and non-equilibrium initial conditions, as the action of the CMD Liouvillian on the Boltzmann distribution,

L^{(CMD)} e^{- β ({\bar{P}}_{t}^{2} / 2 m + W^{(0)} ({\bar{Q}}_{t}))} = e^{- β ({\bar{P}}_{t}^{2} / 2 m + W^{(0)} ({\bar{Q}}_{t}))} - \frac{{\bar{P}}_{t}}{m} (\frac{\partial W^{(0)} ({\bar{Q}}_{t})}{\partial {\bar{Q}}_{t}} + F_{1} ({\bar{Q}}_{t})),

(23)

does not involve the inter-bead harmonic spring-term.⁷⁰ For equilibrium initial conditions, preservation of the Matsubara phase $θ_{M}$ in RPMD and CMD ensures that the quantum Boltzmann distribution is preserved in the dynamics; for non-equilibrium initial conditions, it ensures that for systems with dissipation, the system correctly relaxes to the equilibrium thermal distribution associated with ${\hat{H}}^{(1)}$ ⁠.

Finally, for non-equilibrium initial conditions associated with a sudden vertical excitation, we calculate the orders in time to which RPMD and CMD agree with Matsubara dynamics, which serves as a proxy for the exact quantum dynamics. It can be shown that for the case of equilibrium initial conditions, the orders to which RPMD and partially adiabatic CMD are consistent with Matsubara dynamics are the same orders to which both methods agree with the exact results.^71,74 We thus expect that Matsubara dynamics is a useful proxy for the exact quantum results in the current analysis of the short-time accuracy of RPMD and CMD.

In comparing to Matsubara dynamics, we first consider the case of RPMD. We expand the Matsubara dynamics time-evolution operator as

e^{- L_{M} t} = \sum_{j} \frac{{(- (L_{M}^{(RP)} + i L_{M}^{(I)}) t)}^{j}}{j!}

(24)

and obtain the highest power j for which ${(L_{M}^{(RP)})}^{j} B_{1} (\tilde{𝐐}) B_{2} (\tilde{𝐏})$ agrees with ${(L_{M})}^{j} B_{1} (\tilde{𝐐}) B_{2} (\tilde{𝐏})$ ⁠. To this end, we can exploit the property that $L_{M}^{(I)} B_{1} (\tilde{𝐐}) B_{2} (\tilde{𝐏}) = 0$ ⁠, as long as the function $L_{M}^{(I)}$ acts on can be written as a product of a pure, permutationally invariant function of $\tilde{𝐐}$ (such as the centroid, $B_{1} (\tilde{𝐐}) = \bar{q} (\tilde{𝐐})$ ⁠) times a similarly invariant function of $\tilde{𝐏}$ ⁠.⁶⁹ Thus, we need to identify the lowest order j for which ${(L_{M}^{(RP)})}^{j}$ generates a mixed function in $\tilde{𝐐}$ and $\tilde{𝐏}$ ⁠. Specifically, for the calculation of the non-equilibrium position-autocorrelation function (⁠ $B_{1} (\tilde{𝐐}) B_{2} (\tilde{𝐏}) = \bar{q} (\tilde{𝐐})$ ⁠), we find

L_{M}^{(RP)} \bar{q} (\tilde{𝐐}) = \frac{\bar{P}}{m},

(25)

L_{M}^{{(RP)}^{2}} \bar{q} (\tilde{𝐐}) = - \frac{1}{m} \frac{\partial {\tilde{U}}_{M} (\tilde{𝐐})}{\partial \bar{q} (\tilde{𝐐})},

(26)

L_{M}^{{(RP)}^{3}} \bar{q} (\tilde{𝐐}) = - \frac{1}{m^{2}} \sum_{n = - (M - 1) / 2}^{(M - 1) / 2} {\tilde{P}}_{n} \frac{\partial^{2} {\tilde{U}}_{M} (\tilde{𝐐})}{\partial \bar{q} (\tilde{𝐐}) \partial {\tilde{Q}}_{n}},

(27)

where M is the number of Matsubara modes defined in Appendix C. To obtain the order to which the RPMD non-equilibrium position-autocorrelation function agrees with Matsubara dynamics, we can use integration by parts. We find that acting with $L_{M}^{(RP)}$ on $\bar{q} (\tilde{𝐐}) e^{- β_{N} H_{N}^{(0)} (𝐩, 𝐪)}$ generates a mixed function of Q and P (see Eq. (22)). Therefore, for the non-equilibrium position-autocorrelation function, RPMD agrees with Matsubara dynamics up to t⁴, as compared with t⁶ for equilibrium initial conditions, which can be derived using the same scheme outlined above. The same technique can be used to generate the $t \to 0$ properties of other TCFs of Hermitian operators that are pure functions of $\hat{q}$ or $\hat{p}$ ⁠, as summarized in Table I.

TABLE I.

Orders in time to which RPMD and CMD reproduce Matsubara dynamics for non-equilibrium TCFs ${\tilde{C}}_{A B} (t)$ for a sudden vertical excitation (VE) or an initial momentum impulse (MI). For comparison, results for equilibrium (Eq) TCFs are given.

		CMD			RPMD
$\hat{A}$	$\hat{B}$	VE	MI	Eq	VE	MI	Eq
$\hat{q}$	$\hat{q}$	t¹	t²	t³	t⁴	t⁵	t⁶
$\hat{p}$	$\hat{p}$	t⁰	t¹	t¹	t³	t⁴	t⁴
$A (\hat{q})$	$B (\hat{q})$	$\dots$	$\dots$	$\dots$	t²	t²	t²

		CMD			RPMD
$\hat{A}$	$\hat{B}$	VE	MI	Eq	VE	MI	Eq
$\hat{q}$	$\hat{q}$	t¹	t²	t³	t⁴	t⁵	t⁶
$\hat{p}$	$\hat{p}$	t⁰	t¹	t¹	t³	t⁴	t⁴
$A (\hat{q})$	$B (\hat{q})$	$\dots$	$\dots$	$\dots$	t²	t²	t²

To characterize the short-time behavior of CMD for the calculation of non-equilibrium TCFs, we apply the CMD Liouvillian to $B_{1} (\tilde{𝐐}) B_{2} (\tilde{𝐏}) = \bar{q} (\tilde{𝐐})$ once, which gives the Matsubara result to order t¹. The results of using integration by parts and acting on $\bar{q} (\tilde{𝐐}) e^{- β ({\bar{P}}_{0}^{2} / 2 m + W^{(0)} ({\bar{Q}}_{0}))}$ disagree with the Matsubara dynamics result. Hence, CMD reproduces the non-equilibrium position-autocorrelation function obtained from Matsubara dynamics up to t¹,⁷⁵ as compared with t³ for equilibrium initial conditions, which can be derived using the same scheme outlined above. The short-time behavior of other TCFs is summarized in Table I.

Finally, we emphasize that the formal short-time accuracy of RPMD and CMD is not necessarily indicative of the long-time accuracy of the methods for the calculation of a given TCF;⁷¹ for this reason, we supplement the formal analysis presented here with the numerical examples presented in Sec. III

2. Non-equilibrium initial conditions associated with a momentum impulse

We now consider the case for which the system is subjected to an initial momentum impulse, such that the Hamiltonians ${\hat{H}}^{(0)}$ and ${\hat{H}}^{(1)}$ in Eq. (12) differ only in terms of their kinetic energy,

{\hat{H}}^{(0)} = \frac{{(\hat{p} - Δ p)}^{2}}{2 m} + V (\hat{q}),

(28)

{\hat{H}}^{(1)} = \frac{{\hat{p}}^{2}}{2 m} + V (\hat{q}) .

(29)

Formally, this is a gauge transformation, such that the eigenvalues of ${\hat{H}}^{(0)}$ are unchanged with respect to the eigenvalues of ${\hat{H}}^{(1)}$ ⁠, and the eigenfunctions of the two Hamiltonians are related by a factor of $\exp (i Δ p q / ℏ)$ ⁠. Thus, the system described by ${\hat{H}}^{(0)}$ corresponds to the system described by ${\hat{H}}^{(1)}$ with an additional $Δ p$ momentum impulse (e.g., as a result of bond-cleavage). It should be noted that more general changes in the momentum function can destroy the simple ring-polymer form of the Hamiltonian ${\hat{H}}^{(1)}$ ⁠.^76,77

The ring-polymer Hamiltonian corresponding to Eqs. (28) and (29) can be derived using the usual Trotter splitting to give

H_{N}^{(0)} (𝐩, 𝐪) = \sum_{i = 1}^{N} \frac{{(p_{i} - Δ p)}^{2}}{2 m} + \sum_{i = 1}^{N} [\frac{m {(q_{i} - q_{i - 1})}^{2}}{2 {(β_{N} ℏ)}^{2}} + V (q_{i})],

(30)

H_{N}^{(1)} (𝐩, 𝐪) = \sum_{i = 1}^{N} \frac{p_{i}^{2}}{2 m} + \sum_{i = 1}^{N} [\frac{m {(q_{i} - q_{i - 1})}^{2}}{2 {(β_{N} ℏ)}^{2}} + V (q_{i})] .

(31)

Note that shifting the momentum operator in ${\hat{H}}^{(0)}$ leads to a corresponding shift in the momentum of the centroid in $H_{N}^{(0)} (𝐩, 𝐪)$ ⁠. The RPMD approximation for this case is given in Eq. (17), and the equations of motion are given in Eqs. (6) and (7).

The CMD approximation for this case is

{\tilde{C}}_{A B}^{(CMD)} (t) = \frac{1}{2 π ℏ} \int d {\bar{Q}}_{0} \int d {\bar{P}}_{0} e^{- β [{({\bar{P}}_{0} - Δ p)}^{2} / 2 m + W ({\bar{Q}}_{0})]} \times A ({\bar{Q}}_{0}, {\bar{P}}_{0}) B ({\bar{Q}}_{t}, {\bar{P}}_{t}) .

(32)

The centroid position and momenta evolve according to the equations of motion given in Eqs. (10) and (11).

Following the same steps outlined for the case of sudden vertical excitation, one can show that RPMD and CMD correctly reduce to the classical limit, reproduce non-equilibrium TCFs of linear operators in harmonic potentials, and retain their formal relation to Matsubara dynamics. The $t \to 0$ behavior of RPMD and CMD compared to Matsubara dynamics for an initial momentum impulse is given in Table I.

III. RESULTS AND DISCUSSION

A. 1D test systems

In this section, we provide numerical tests of RPMD and CMD for the calculation of non-equilibrium TCFs in one-dimensional potentials. Beginning with the case of sudden vertical excitation, we consider mildly anharmonic potentials of the form

V^{(0)} (q) = \frac{1}{2} {(q - Δ q)}^{2} + \frac{1}{10} {(q - Δ q)}^{3} + \frac{1}{100} {(q - Δ q)}^{4},

(33)

V^{(1)} (q) = \frac{1}{2} q^{2} + \frac{1}{10} q^{3} + \frac{1}{100} q^{4},

(34)

and strongly anharmonic quartic oscillators of the form

V^{(0)} (q) = \frac{1}{4} {(q - Δ q)}^{4},

(35)

V^{(1)} (q) = \frac{1}{4} q^{4} .

(36)

Note that reduced units (⁠ $ℏ = 1$ ⁠, m = 1, and k_B = 1) are used throughout this section. Non-equilibrium initial conditions are simulated by shifting the initial potential V⁽⁰⁾ by Δq with respect to the final potential V⁽¹⁾. The inverse temperature β is chosen as either 1 or 8, and converged path-integral results are obtained with n = 4 and n = 48 beads, respectively. For RPMD, a modified velocity-Verlet scheme³¹ is used to integrate the equations of motion using a time step of 0.05 a.u., and 5 × 10⁵ trajectories are used to converge each set of results. Initial distributions are obtained by running a long trajectory and periodically resampling the ring-polymer momenta from the Maxwell-Boltzmann distribution. Classical results are obtained by running 5 × 10⁵ trajectories, and CMD results are obtained by running 1 × 10⁶ trajectories.

Fig. 1 presents the non-equilibrium position-autocorrelation function ${\tilde{C}}_{q q} (t)$ for the mildly anharmonic potential at low temperatures and a shift of $Δ q = 0.3$ a.u. For comparison, the equilibrium (⁠ $Δ q = 0.0$ a.u.) position-autocorrelation functions are also provided, which have been previously studied.^27,29,78 While the classical mechanical results for the non-equilibrium TCFs do not match the exact quantum mechanical results at t = 0 and do not capture the correct oscillation frequency, both RPMD and CMD agree with the exact results. Overall, the accuracy for calculating the non-equilibrium TCFs is similar to the equilibrium case.

FIG. 1.

View large Download slide

Non-equilibrium (top) and equilibrium (bottom) position-autocorrelation function at low temperature (⁠ $β = 8$ ⁠) for the mildly anharmonic potential. Quantum results are shown in solid black, classical results as green long-dashed lines, RPMD results as blue short-dashed lines, and CMD results as orange dotted lines. Non-equilibrium initial conditions are introduced by sudden vertical excitation using $Δ q = 0.3$ a.u.

For the case of the quartic oscillator (Fig. 2), RPMD and CMD again perform similarly for equilibrium and non-equilibrium TCFs. As RPMD and CMD neglect the real-time coherences needed to fully treat the dynamics in this potential, they both diverge from the exact results after several oscillations. CMD captures more of the long-time behavior, which has been studied previously for equilibrium initial conditions.⁷⁹

FIG. 2.

View large Download slide

Same as Fig. 1, but for the strongly anharmonic quartic oscillator potential.

To quantify the results over a larger range of non-equilibrium initial conditions, we introduce a measure for the error of a given approximate method with respect to the exact quantum mechanical results,

E_{A B} (β Δ E) = \frac{\int_{0}^{T} d t {({\tilde{C}}_{A B}^{(exact)} (t) - {\tilde{C}}_{A B}^{(X)} (t))}^{2}}{\int_{0}^{T} d t {({\tilde{C}}_{A B}^{(exact)} (t))}^{2}},

(37)

where X indicates the approximate method. The term $β Δ E$ is a measure of how far the system is initially out of equilibrium; it is defined by the energy difference

Δ E = {⟨ {\hat{H}}^{(1)} ⟩}_{{\hat{H}}^{(0)}} - {⟨ {\hat{H}}^{(1)} ⟩}_{{\hat{H}}^{(1)}},

(38)

where ${⟨ \dots ⟩}_{{\hat{H}}^{(j)}}$ indicates equilibrium thermal averaging with the indicated Hamiltonian ${\hat{H}}^{(j)}$ ⁠. This energy difference is scaled by β to provide a unitless quantity. Note that $β Δ E = 0$ corresponds to the results for equilibrium initial conditions.

In Fig. 3, the error of the approximate methods is shown for the mildly anharmonic potential and non-equilibrium initial conditions associated with vertical excitation (Fig. 3(a)) or a momentum impulse (Fig. 3(b)). Up to ten times the thermal energy of the system is added via the non-equilibrium initial conditions. First, the case of a sudden vertical excitation is discussed (Fig. 3(a)). In the high temperature case (⁠ $β = 1$ ⁠), all methods yield very small errors. For all non-equilibrium cases the error is lower than that for the equilibrium case, as the error decreases as a function of $β Δ E$ ⁠. In the low temperature regime, the results for the various approximate methods coalesce at high $β Δ E$ ⁠, but RPMD and CMD perform consistently better than the classical results over the range of non-equilibrium initial conditions. For all non-equilibrium initial conditions considered, the approximation error for RPMD and CMD is similar or lower than for equilibrium initial conditions. The same is true for an initial momentum impulse (see Fig. 3(b)).

FIG. 3.

View large Download slide

Approximation error for position-autocorrelation functions (see. Eq. (37)) for the mildly anharmonic potential as a function of the energy added to the system due to the non-equilibrium initial conditions (see Eq. (38)). (a) Initial conditions introduced by sudden vertical excitation (⁠ $Δ q \neq 0$ ⁠). (b) Initial conditions introduced by an initial momentum impulse.

Finally, an example of a non-linear TCF is considered. Fig. 4 displays the q² autocorrelation function for the mildly anharmonic potential at low temperature. Non-equilibrium initial conditions are introduced by a sudden vertical excitation and $Δ q = 0.3$ a.u. Again, RPMD performs similarly for the non-equilibrium and equilibrium cases, and it is correct in the short-time limit. For this case, CMD performs more similarly to classical dynamics, although various methods to alleviate this problem have been proposed.^25,80

FIG. 4.

View large Download slide

Same as Fig. 1, but for the q² autocorrelation function.

B. Excited-state proton transfer

To further test the applicability of RPMD for non-equilibrium reactions in the condensed phase, the dynamics of a double-well system coupled to a dissipative bath is studied. The potentials V⁽⁰⁾ and V⁽¹⁾ take the following form:

V^{(j)} (q, 𝐱) = V_{s}^{(j)} (q) + V_{b} (q, 𝐱),

(39)

where

V_{s}^{(0)} (q) = {\begin{matrix} - \frac{1}{2} m ω_{b}^{2} q^{2} + \frac{m^{2} ω_{b}^{4}}{16 V_{0}^{‡}} q^{4} & if q < 0 \\ \infty & otherwise \end{matrix},

(40)

V_{s}^{(1)} (q) = - \frac{1}{2} m ω_{b}^{2} q^{2} + \frac{m^{2} ω_{b}^{4}}{16 V_{1}^{‡}} q^{4},

(41)

V_{b} (q, 𝐱) = \sum_{j = 1}^{N} \frac{1}{2} ω_{j}^{2} {(x_{j} - \frac{c_{j} q}{m ω_{j}^{2}})}^{2},

(42)

m = 1836 a.u., the barrier frequency is $ω_{b} = 500 {cm}^{- 1}$ ⁠, and the barrier height in the excited potential $V_{1}^{‡}$ is chosen such that the zero-point corrected barrier height remains approximately 7 k_BT. As the locations of the potential minima depend on the barrier height, $x_{min} = \pm \sqrt{4 V_{j}^{‡} / m ω_{b}^{2}}$ ⁠, different values of $V_{0}^{‡}$ and $V_{1}^{‡}$ create non-equilibrium initial conditions during the sudden vertical excitations from V⁽⁰⁾ to V⁽¹⁾. As shown in Fig. 5, the potential energy surfaces differ both in shape and in the position of the reactant minima; see Table II for the employed parameters. For comparison, calculations with equilibrium initial conditions (i.e., $V_{0}^{‡} = V_{1}^{‡}$ ⁠) are also performed. The system potentials are coupled to a dissipative bath with a Debye spectral density,

J_{D} (ω) = \frac{η γ ω}{ω^{2} + γ^{2}},

(43)

a low system-bath coupling of $η = 0.2 m ω_{b}$ ⁠, and a bath cutoff frequency of $γ = 500$ cm⁻¹ taken from the DW1 model.⁸¹ This system-bath coupling value is chosen to be low enough to exhibit substantial non-equilibrium effects, but high enough to limit the number of trajectory recrossings.⁸¹ The bath is discretized as⁸²

ω_{j} = γ \tan (\frac{π}{2} \frac{j}{N + 1}),

(44)

c_{j} = ω_{j} \sqrt{\frac{m η γ}{N + 1}} .

(45)

FIG. 5.

View large Download slide

System potentials $V_{s}^{(0)} (q)$ and $V_{s}^{(1)} (q)$ employed in the system-bath model for $V_{0}^{‡} = 3000$ cm⁻¹ and $V_{1}^{‡} = 1500$ cm⁻¹ (see Eqs. (40) and (41)). For better visibility the potentials are shifted vertically, and the arrow illustrates sudden vertical excitation from the ground-state reactant.

TABLE II.

Parameters employed in the excited-state proton transfer simulations. Energies reported in units of cm⁻¹.

T (K)	$β Δ E$	$V_{0}^{‡}$	$V_{1}^{‡}$
230	0	1500	1500
230	6	2600	1500
230	9	3000	1500
77	0	700	700
77	6	1150	700
77	9	1300	700

The non-equilibrium time-dependent side expectation value is calculated as

⟨ h ⟩ (t) = \frac{1}{Z} Tr [e^{- β {\hat{H}}^{(0)}} e^{- i {\hat{H}}^{(1)} t / ℏ} h (\hat{q}) e^{i {\hat{H}}^{(1)} t / ℏ}],

(46)

where h is a Heaviside function,

h (q) = {\begin{matrix} 1 & if q < 0 \\ 0 & else \end{matrix}

(47)

The transfer time τ for the proton transfer reaction is obtained from a linear fit of the side expectation value between t₁ = 200 fs and t₂ = 600 fs. Calculations are performed at two different temperatures (77 K and 230 K), such that the system is either above or below the cross-over temperature of $T_{c}$ = 115 K. N = 50 bath modes are used to converge the bath discretization. The RPMD trajectories are propagated using a step size of $Δ t = 0.2$ fs and either n = 16 (T = 230 K) or n = 64 (T = 77 K) ring-polymer beads. For each reported result, 10⁶ RPMD trajectories are performed, with initial conditions sampled via Monte Carlo.

For comparison, numerically exact quantum mechanical results are obtained using a newly developed version of the Liouville space hierarchical equation of motion (HEOM) method.^83–85 The HEOM approach is both non-perturbative and capable of describing non-Markovian effects of the bath.^86–88 The present simulations employ a tolerance filter for the auxiliary density operators⁸⁹ and the [R − 1/R] Padé decomposition scheme.⁹⁰ It was found that the calculations converge at R = 4 for the T = 230 K case and R = 5 for the T = 77 K case.

Fig. 6 presents the time-dependent non-equilibrium side expectation value for the proton-transfer reaction at T = 230 K $> T_{c}$ ⁠; the corresponding reaction transfer times are given in Table III. For the equilibrium case (Fig. 6(a), $β Δ E = 0$ ⁠), good agreement between the RPMD and the exact HEOM results is found. The transfer time obtained from classical dynamics differs by a factor of three. These findings are consistent with previous studies of thermal rate constants of double-well systems using RPMD.³⁰ Reasonable agreement between RPMD and exact results is also found for two cases with non-equilibrium initial conditions (Figs. 6(b) and 6(c), $β Δ E > 0$ ⁠). The amount of initial population transfer, as well as the subsequent transfer times, compare well for both cases. The classical transfer times differ from the exact results by approximately a factor of two, while the RPMD results agree to within 20%.

FIG. 6.

View large Download slide

Non-equilibrium side expectation value for a double well system coupled to harmonic bath at T = 230 K. The system is initialized in one of the two wells and reaction to the other well is monitored. Quantum results are shown in solid black, classical results as green short-dashed lines, and RPMD results as blue long-dashed lines. The numbers indicate the transfer time τ in ps. (a) Equilibrium initial conditions, (b) non-equilibrium initial conditions adding 6 k_BT of energy to the system, (c) non-equilibrium initial conditions adding 9 k_BT of energy to the system.

TABLE III.

Transfer times τ (in ps) and corresponding standard errors at 230 K obtained from linear fits to the non-equilibrium side expectation value to $f (t) = A (1 - \frac{t}{τ})$ between 200 fs and 600 fs. The number in parenthesis indicates the error in the last reported digit. Standard errors are obtained from fitting to non-equilibrium side expectation values obtained from ten independent subsets of trajectories.

$β Δ E$	Classical	RPMD	HEOM
0	8.2(1) $\times 10^{2}$	2.89(4) $\times 10^{2}$	2.82 $\times 10^{2}$
6	2.33(1) $\times 10^{2}$	1.31(3) $\times 10^{2}$	1.25 $\times 10^{2}$
9	2.35(2) $\times 10^{2}$	1.7(1) $\times 10^{2}$	1.39 $\times 10^{2}$

$β Δ E$	Classical	RPMD	HEOM
0	8.2(1) $\times 10^{2}$	2.89(4) $\times 10^{2}$	2.82 $\times 10^{2}$
6	2.33(1) $\times 10^{2}$	1.31(3) $\times 10^{2}$	1.25 $\times 10^{2}$
9	2.35(2) $\times 10^{2}$	1.7(1) $\times 10^{2}$	1.39 $\times 10^{2}$

Fig. 7 and Table IV present the result for the deep tunneling regime, T = 77 K $< T_{c}$ ⁠. For the equilibrium initial conditions (Fig. 7(a)), the classical transfer time deviates by over three orders of magnitude from the exact results. However, the RPMD transfer time remains within a factor of four of the exact results. The degree to which RPMD overestimates the transfer time in this equilibrium deep-tunneling case is consistent with earlier analysis.⁴⁴ For the two non-equilibrium cases shown in Figs. 7(b) and 7(c), the RPMD results again compare well with the exact results. The initial population transfer is accurately reproduced, and the RPMD transfer times agree with the exact results to within a factor of four. As seen for the equilibrium case, the classical transfer times deviate from the exact results by over two orders of magnitude.

FIG. 7.

View large Download slide

Same as Fig. 6, but for T = 77 K.

TABLE IV.

Same as Table III but for 77 K.

$β Δ E$	Classical	RPMD	HEOM
0	2.5(1) $\times 10^{4}$	3.23(7) $\times 10^{1}$	8.3 $\times 10^{0}$
6	2.19(6) $\times 10^{3}$	3.3(1) $\times 10^{1}$	9.9 $\times 10^{0}$
9	1.50(4) $\times 10^{3}$	4.2(3) $\times 10^{1}$	1.1 $\times 10^{1}$

$β Δ E$	Classical	RPMD	HEOM
0	2.5(1) $\times 10^{4}$	3.23(7) $\times 10^{1}$	8.3 $\times 10^{0}$
6	2.19(6) $\times 10^{3}$	3.3(1) $\times 10^{1}$	9.9 $\times 10^{0}$
9	1.50(4) $\times 10^{3}$	4.2(3) $\times 10^{1}$	1.1 $\times 10^{1}$

IV. CONCLUDING REMARKS

In this paper, we demonstrate that both RPMD and CMD can be used to approximate non-equilibrium TCFs and non-equilibrium time-dependent expectation values associated with a sudden vertical excitation or an initial momentum impulse. Both methods are exact at high-temperatures and in the classical limit. For the calculation of non-equilibrium TCFs of linear operators, both methods are exact for t = 0 and for harmonic potentials; RPMD is also exact for the calculation of TCFs of general non-linear operators for t = 0. For both methods, the connection to Matsubara dynamics found for equilibrium initial conditions^69,70 is preserved for non-equilibrium initial conditions. Furthermore, the orders in time to which RPMD and CMD agree with Matsubara dynamics are determined. Specifically, for the position-autocorrelation function associated with sudden vertical excitation, RPMD and CMD agree with Matsubara dynamics up to $O (t^{4})$ and $O (t^{1})$ ⁠, respectively; for the position-autocorrelation function associated with an initial momentum impulse, RPMD and CMD agree with Matsubara dynamics up to $O (t^{5})$ and $O (t^{2})$ ⁠, respectively. Similarly, the short-time comparison of RPMD and CMD to Matsubara dynamics is derived for more general TCFs as presented in Table I. Extensive numerical tests employing one-dimensional models show that RPMD and CMD give similar accuracy for calculating equilibrium and non-equilibrium correlation functions. Furthermore, the applicability of RPMD to non-equilibrium excited-state proton transfer processes is assessed using a system-bath model. Both above and below the cross-over temperature for deep tunneling of the transferring proton, good agreement is found between RPMD and numerically exact quantum dynamical results, even for cases in which the corresponding classical results are in error by over two orders of magnitude. The accuracy of RPMD for non-equilibrium initial conditions is found to be similar to that for equilibrium initial conditions. The results provided here indicate that the path-integral based methods allow for the approximate quantum dynamical study of photo-excited reactions in complex systems.⁹¹ In future work, it will be worth determining whether non-adiabatic extensions of RPMD^92–96 are similarly successful for the calculation of non-equilibrium TCFs.

ACKNOWLEDGMENTS

R.W. acknowledges financial support from the Deutsche Forschungsgemeinschaft under Grant No. We 5762/1-1. S.C.A. acknowledges support from the Leverhulme Trust, the Chinese Academy of Sciences visiting professorship for senior international scientists (Grant No. 2013T2J0022) program, and the hospitality of the Miller group at Caltech. Q.S. acknowledges support from the National Natural Science Foundation of China (Grant No. 21290194). T.F.M. acknowledges support from the National Science Foundation (NSF) CAREER Award under Grant No. CHE-1057112 and from the Office of Naval Research (ONR) under Grant No. N00014-16-1-2761.

APPENDIX A: CONNECTION TO CONDON APPROXIMATION

Let $| Ψ_{j} ⟩$ be the jth vibrational eigenstate of the Hamiltonian of the ground electronic state ${\hat{H}}^{(0)} = \hat{T} + {\hat{V}}^{(0)}$ ⁠. The state-resolved expectation value of an arbitrary operator $\hat{A}$ is

{⟨ A ⟩}^{(j)} = ⟨ Ψ_{j} | \hat{A} | Ψ_{j} ⟩,

(A1)

and the equilibrium thermal expectation value of operator $\hat{A}$ is given by

⟨ A ⟩ = \frac{1}{Z} \sum_{j} e^{- β E_{j}} ⟨ Ψ_{j} | \hat{A} | Ψ_{j} ⟩ = \frac{1}{Z} \sum_{j} e^{- β E_{j}} {⟨ A ⟩}^{(j)} = \frac{1}{Z} Tr (e^{- β {\hat{H}}^{(0)}} \hat{A}) .

(A2)

Invoking the Condon approximation (i.e., vertical excitation), the time-dependent expectation value of operator $\hat{A}$ following the excitation of the initial thermal distribution to V⁽¹⁾ is

⟨ A ⟩ (t) = \frac{1}{Z} \sum_{j} e^{- β E_{j}} ⟨ Ψ_{j} | e^{i {\hat{H}}^{(1)} t / ℏ} \hat{A} e^{- i {\hat{H}}^{(1)} t / ℏ} | Ψ_{j} ⟩ = \frac{1}{Z} Tr (e^{- β {\hat{H}}^{(0)}} e^{i {\hat{H}}^{(1)} t / ℏ} \hat{A} e^{- i {\hat{H}}^{(1)} t / ℏ}),

(A3)

with ${\hat{H}}^{(1)}$ being the Hamiltonian of the excited electronic state, ${\hat{H}}^{(1)} = \hat{T} + {\hat{V}}^{(1)}$ ⁠. The relation given in Eq. (A3) corresponds to the case of the non-equilibrium correlation function ${\tilde{C}}_{1 A} (t)$ from Eqs. (12)–(14).

Similarly, we can take

{\tilde{C}}_{A B} (0) = \frac{1}{β Z} \int_{0}^{β} d λ \sum_{i, j} ⟨ Ψ_{j} | e^{- λ {\hat{H}}^{(0)}} \hat{A} e^{- (β - λ) {\hat{H}}^{(0)}} | Ψ_{i} ⟩ \times ⟨ Ψ_{i} | \hat{B} | Ψ_{j} ⟩ = \frac{1}{β Z} \int_{0}^{β} d λ Tr (e^{- λ {\hat{H}}^{(0)}} \hat{A} e^{- (β - λ) {\hat{H}}^{(0)}} \hat{B})

(A4)

and invoke the Condon approximation to give the non-equilibrium TCF

{\tilde{C}}_{A B} (t) = \frac{1}{β Z} \int_{0}^{β} d λ \sum_{i, j} ⟨ Ψ_{j} | e^{- λ {\hat{H}}^{(0)}} \hat{A} e^{- (β - λ) {\hat{H}}^{(0)}} | Ψ_{i} ⟩ \times ⟨ e^{- i {\hat{H}}^{(1)} t / ℏ} Ψ_{i} | \hat{B} | e^{- i {\hat{H}}^{(1)} t / ℏ} Ψ_{j} ⟩ = \frac{1}{β Z} \int_{0}^{β} d λ Tr (e^{- λ {\hat{H}}^{(0)}} \hat{A} e^{- (β - λ) {\hat{H}}^{(0)}} e^{i {\hat{H}}^{(1)} t / ℏ} \hat{B} e^{- i {\hat{H}}^{(1)} t / ℏ}) .

(A5)

The relation given in Eq. (A5) corresponds to the case of the non-equilibrium Kubo-transformed correlation function ${\tilde{C}}_{A B} (t)$ from Eqs. (12)–(14). There is currently no general transformation between Kubo-transformed non-equilibrium TCFs and standard non-equilibrium TCFs.⁸⁰ However, this poses no problem for the calculation of time-dependent non-equilibrium expectation values, as is often of interest for the study of non-equilibrium chemical processes.

APPENDIX B: HARMONIC LIMIT

First the case of a vertical excitation is considered with H⁽⁰⁾ = T + V⁽⁰⁾ and H⁽¹⁾ = T + V⁽¹⁾, such that V⁽⁰⁾ is arbitrary and $V^{(1)} = \frac{1}{2} k^{2} q^{2}$ ⁠. The Kubo-transformed position-autocorrelation function can be written as

{\tilde{C}}_{q q} (t) = \frac{1}{Z} lim_{N \to \infty} \int d T N Q T N Q 𝐪 \int d 𝚫 \frac{1}{N} \sum_{k = 0}^{N - 1} q_{k} \times \prod_{j = 0}^{N - 1} ⟨ q_{j - 1} - \frac{Δ_{j - 1}}{2} | e^{- β_{N} {\hat{H}}^{(0)}} | q_{j} + \frac{Δ_{j}}{2} ⟩ \times ⟨ q_{j} + \frac{Δ_{j}}{2} | e^{i {\hat{H}}^{(1)} t / ℏ} \hat{q} e^{- i {\hat{H}}^{(1)} t / ℏ} | q_{j} - \frac{Δ_{j}}{2} ⟩ .

(B1)

Using linearization of the difference between the forward and backward Feynman paths in real time, which is exact for a harmonic potential V⁽¹⁾,^2,97 one obtains

{\tilde{C}}_{q q} (t) = \frac{1}{Z} lim_{N \to \infty} \frac{1}{{(2 π ℏ)}^{N}} \int d 𝐪 \int d 𝐩 \int d 𝚫 \times \frac{1}{N^{2}} \sum_{k = 0}^{N - 1} q_{k} \sum_{l = 0}^{N - 1} (q_{l} \cos (k t) + \frac{p_{l}}{m k} \sin (k t)) \times \prod_{j = 0}^{N - 1} e^{- i p_{j} Δ_{j} / ℏ} ⟨ q_{j - 1} - \frac{Δ_{j - 1}}{2} | e^{- β_{N} {\hat{H}}^{(0)}} | q_{j} + \frac{Δ_{j}}{2} ⟩ .

(B2)

Following Hele and Althorpe,⁹⁸ we transform to normal modes and do the integration over the N $-$ 1 delta functions in $D_{1} \dots D_{N - 1}$ ⁠, obtaining

{\tilde{C}}_{q q} (t) = \frac{1}{Z} lim_{N \to \infty} \frac{1}{{(2 π ℏ)}^{N}} \int d 𝐪 \int d P_{0} \frac{1}{N^{2}} \sum_{k = 0}^{N - 1} q_{k} \times \sum_{l = 0}^{N - 1} q_{l} \cos (k t) + \frac{\sqrt{N} P_{0}}{m k} \sin (k t) \times \sqrt{\frac{2 π β_{N} ℏ^{2}}{m}} e^{\frac{- β_{N}}{2 m} P_{0}^{2}} \prod_{j = 0}^{N - 1} ⟨ q_{j - 1} | e^{- β_{N} {\hat{H}}^{(0)}} | q_{j} ⟩ .

(B3)

In the limit of $N \to \infty$ ⁠,

⟨ q_{j - 1} | e^{- β_{N} {\hat{H}}^{(0)}} | q_{j} ⟩ = \sqrt{\frac{m}{2 π β ℏ_{2}^{N}}} \times e^{- β_{N} (\frac{m}{2 β_{N}^{2} ℏ} {(q_{j} - q_{j - 1})}^{2} + \frac{1}{2} (V^{(0)} (q_{j}) + V^{(0)} (q_{j - 1})))} .

(B4)

Furthermore, each term $\sqrt{\frac{m}{2 π β ℏ_{2}^{N}}}$ can be rewritten as ${(2 π ℏ)}^{- 1}$ $\int d p_{j} e^{\frac{- β_{N}}{2 m} P_{j}^{2}}$ ⁠; inserting these N - 1 terms into Eq. (B3) and back-transforming from normal modes yields

{\tilde{C}}_{q q} (t) = lim_{N \to \infty} \frac{1}{Z_{N}} \frac{1}{{(2 π ℏ)}^{N}} \int d 𝐪 \int d 𝐩 \frac{1}{N^{2}} \sum_{k = 0}^{N - 1} q_{k} \times \sum_{l = 0}^{N - 1} (q_{l} \cos (k t) + \frac{p_{l}}{m k} \sin (k t)) \times \prod_{j = 0}^{N - 1} e^{- β_{N} (\frac{p_{j}^{2}}{2 m} + \frac{m}{2 β_{N}^{2} ℏ^{2}} {(q_{j} - q_{j - 1})}^{2} + V^{(0)} (q_{j}))} .

(B5)

This exactly matches the RPMD result presented in Eq. (17) of the main text.

Similarly, for the case of a momentum impulse with $H^{(0)} = \frac{1}{2 m σ_{p}} {(\hat{p} - Δ p)}^{2} + \frac{1}{2} m k^{2} {\hat{q}}^{2}$ and $H^{(1)} = \frac{1}{2 m} p^{2} + \frac{1}{2} m k^{2} {\hat{q}}^{2}$ ⁠, we find for the non-equilibrium momentum-autocorrelation function

{\tilde{C}}_{p p} (t) = \frac{1}{Z} lim_{N \to \infty} \frac{1}{{(2 π ℏ)}^{N}} \int d 𝐪 \int d 𝐩 \frac{1}{N^{2}} \sum_{k = 0}^{N - 1} p_{k} \times \sum_{l = 0}^{N - 1} (- k m q_{l} \sin (k t) + \frac{p_{l}}{m} \cos (k t)) \times \prod_{j = 0}^{N - 1} e^{- β_{N} (\frac{{(p_{j} - Δ p)}^{2}}{2 m σ_{p}} + \frac{m σ_{p}}{2 β_{N}^{2} ℏ^{2}} {(q_{j} - q_{j - 1})}^{2} + \frac{1}{2} m k^{2} q_{j}^{2}) .}

(B6)

This exactly matches the RPMD approach presented in Sec. II B 2 of the main text.

Similar steps can be taken to show that CMD is exact for harmonic potentials, as the centroid and non-centroid modes decouple for harmonic potentials.²⁴

APPENDIX C: REVIEW OF MATSUBARA DYNAMICS

We briefly review the relevant aspects of Matsubara dynamics for the analysis of RPMD and CMD in the current paper. Full details are available in Ref. 69.

Matsubara dynamics approximates the quantum Kubo-transformed time-correlation function of Eq. (1) by

C_{A B}^{Mats} (t) = lim_{M \to \infty} C_{A B}^{(M)} (t),

(C1)

where

C_{A B}^{(M)} (t) = \frac{α_{M}}{2 π ℏ} \int d \tilde{𝐏} \int d \tilde{𝐐} A (\tilde{𝐐}) e^{- β ({\tilde{H}}_{M} (\tilde{𝐏}, \tilde{𝐐}) - i θ_{M} (\tilde{𝐏}, \tilde{𝐐}))} \times e^{L_{M} t} B (\tilde{𝐐})

(C2)

and $α_{M} = ℏ^{(1 - M)} {((M - 1) / 2)!}^{2}$ ⁠. The position coordinates $\tilde{𝐐} \equiv {{\tilde{Q}}_{n}}$ ⁠, with $n = - (M - 1) / 2, \dots, (M - 1) / 2$ ⁠, are the M Matsubara modes defined as⁶⁹

{\tilde{Q}}_{n} = lim_{N \to \infty} \frac{1}{\sqrt{N}} \sum_{l = 1}^{N} T_{l n} q_{l}, n = 0, \pm 1, \dots \pm (M - 1) / 2,

(C3)

where M is odd⁶⁹ and satisfies $M ≪ N$ ⁠; $𝐪 \equiv {q_{l}}, l = 1, \dots, N$ are a set of discrete path-integral coordinates distributed at equally spaced intervals $β ℏ / N$ of imaginary time, and

T_{l n} = {\begin{matrix} N^{- 1 / 2} & n = 0 \\ \sqrt{2 / N} \sin (2 π l n / N) & n = 1, \dots, (M - 1) / 2 \\ \sqrt{2 / N} \cos (2 π l n / N) & n = - 1, \dots, - (M - 1) / 2 \end{matrix} .

(C4)

The Matsubara phase $θ_{M}$ is given by

θ_{M} (\tilde{𝐏}, \tilde{𝐐}) = \sum_{n = - (M - 1) / 2}^{(M - 1) / 2} {\tilde{P}}_{n} {\tilde{ω}}_{n} {\tilde{Q}}_{- n} .

The momentum coordinates $\tilde{𝐏}$ are similarly defined in terms of p. ${\tilde{Q}}_{0} = \bar{q} (\tilde{𝐐})$ and ${\tilde{P}}_{0} = \bar{p} (\tilde{𝐏})$ are the position and momentum centroid coordinates. The associated Matsubara frequencies are ${\tilde{ω}}_{n} = 2 n π / β ℏ$ ⁠.

The functions $A (\tilde{𝐐})$ and $B (\tilde{𝐐})$ in Eq. (C2) are obtained by making the substitutions

q_{l} = \sqrt{N} \sum_{n = - (M - 1) / 2}^{(M - 1) / 2} T_{l n} {\tilde{Q}}_{n}

(C5)

into the functions

A (𝐪) = \frac{1}{N} \sum_{l = 1}^{N} A (q_{l}) and B (𝐪) = \frac{1}{N} \sum_{l = 1}^{N} B (q_{l})

(C6)

The Matsubara potential ${\tilde{U}}_{M} (\tilde{𝐐})$ is obtained similarly by substituting for q_l into the ring-polymer potential

U_{N} (𝐪) = \frac{1}{N} \sum_{l = 1}^{N} V (q_{l}) .

(C7)

The propagator $e^{- L_{M} t}$ contains the Matsubara Liouvillian

L_{M} = \sum_{n = - (M - 1) / 2}^{(M - 1) / 2} \frac{{\tilde{P}}_{n}}{m} \frac{\partial}{\partial {\tilde{Q}}_{n}} - \frac{\partial {\tilde{U}}_{M} (\tilde{𝐐})}{\partial {\tilde{Q}}_{n}} \frac{\partial}{\partial {\tilde{P}}_{n}},

(C8)

The formulas presented above result from just one approximation: decoupling the Matsubara modes from the non-Matsubara modes in the exact quantum Liouvillian (which causes all Liouvillian terms of $O (ℏ^{2})$ to vanish).⁶⁹ Within this assumption, the dynamics conserves the Hamiltonian H, and also the phase $θ_{M} (\tilde{𝐏}, \tilde{𝐐})$ ⁠, and hence the quantum Boltzmann distribution.

One can similarly obtain Matsubara dynamics for non-equilibrium initial conditions for the two cases discussed in the main text. In the case of a sudden vertical excitation, the Matsubara Liouvillian reads

L_{M}^{(1)} = \sum_{n = - (M - 1) / 2}^{(M - 1) / 2} \frac{{\tilde{P}}_{n}}{m} \frac{\partial}{\partial {\tilde{Q}}_{n}} - \frac{\partial {\tilde{U}}_{M}^{(1)} (\tilde{𝐐})}{\partial {\tilde{Q}}_{n}} \frac{\partial}{\partial {\tilde{P}}_{n}},

(C9)

where ${\tilde{U}}_{M}^{(1)}$ is obtained from the excited-state potential

U_{N}^{(1)} (𝐪) = \frac{1}{N} \sum_{l = 1}^{N} V^{(1)} (q_{l}) .

(C10)

In the case of an initial momentum impulse, TCFs are obtained as in Eq. (C2) with

{\tilde{H}}_{M}^{(0)} (\tilde{𝐏}, \tilde{𝐐}) = \frac{{({\tilde{P}}_{0} + Δ p)}^{2}}{2 m} + \frac{1}{2 m} \sum_{\begin{matrix} n = - (M - 1) / 2 \\ n \neq 0 \end{matrix}}^{(M - 1) / 2} {\tilde{P}}_{n}^{2} + U_{M} (\tilde{𝐐}) .

(C11)

For the non-equilibrium initial conditions, the phase $θ_{M} (\tilde{𝐏}, \tilde{𝐐})$ is still conserved.

REFERENCES

1.

E. J.

Heller

,

J. Chem. Phys.

65

,

1289

(

1976

).

https://doi.org/10.1063/1.433238

Google Scholar

Crossref

2.

M.

Hillery

,

R.

O’Connell

,

M.

Scully

, and

E.

Wigner

,

Phys. Rep.

106

,

121

(

1984

).

https://doi.org/10.1016/0370-1573(84)90160-1

Google Scholar

Crossref

3.

H.

Wang

,

X.

Sun

, and

W. H.

Miller

,

J. Chem. Phys.

108

,

9726

(

1998

).

https://doi.org/10.1063/1.476447

Google Scholar

Crossref

4.

W. H.

Miller

,

J. Phys. Chem. A

105

,

2942

(

2001

).

https://doi.org/10.1021/jp003712k

Google Scholar

Crossref

5.

T.

Yamamoto

,

H.

Wang

, and

W. H.

Miller

,

J. Chem. Phys.

116

,

7335

(

2002

).

https://doi.org/10.1063/1.1464539

Google Scholar

Crossref

6.

J. A.

Poulsen

,

G.

Nyman

, and

P. J.

Rossky

,

J. Chem. Phys.

119

,

12179

(

2003

).

https://doi.org/10.1063/1.1626631

Google Scholar

Crossref

7.

Q.

Shi

and

E.

Geva

,

J. Chem. Phys.

118

,

8173

(

2003

).

https://doi.org/10.1063/1.1564814

Google Scholar

Crossref

8.

Q.

Shi

and

E.

Geva

,

J. Phys. Chem. A

107

,

9059

(

2003

).

https://doi.org/10.1021/jp030497+

Google Scholar

Crossref

9.

S.

Bonella

and

D. F.

Coker

,

J. Chem. Phys.

122

,

194102

(

2005

).

https://doi.org/10.1063/1.1896948

Google Scholar

Crossref

PubMed

10.

H.

Torii

,

J. Phys. Chem. A

110

,

9469

(

2006

).

https://doi.org/10.1021/jp062033s

Google Scholar

Crossref

PubMed

11.

R.

Kapral

,

Ann. Rev. Phys. Chem.

57

,

129

(

2006

).

https://doi.org/10.1146/annurev.physchem.57.032905.104702

Google Scholar

Crossref

12.

J. L.

Skinner

,

B. M.

Auer

, and

Y.-S.

Lin

,

Adv. Chem. Phys.

142

,

59

(

2009

).

https://doi.org/10.1002/9780470475935.ch2

Google Scholar

Crossref

13.

J.

Liu

,

W. H.

Miller

,

F.

Paesani

,

W.

Zhang

, and

D. A.

Case

,

J. Chem. Phys.

131

,

164509

(

2009

).

https://doi.org/10.1063/1.3254372

Google Scholar

Crossref

PubMed

14.

J.

Liu

and

W. H.

Miller

,

J. Chem. Phys.

131

,

074113

(

2009

).

https://doi.org/10.1063/1.3202438

Google Scholar

Crossref

PubMed

15.

Y.

Wang

and

J. M.

Bowman

,

Chem. Phys. Lett.

491

,

1

(

2010

).

https://doi.org/10.1016/j.cplett.2010.03.025

Google Scholar

Crossref

16.

E.

Kamarchik

and

J. M.

Bowman

,

J. Phys. Chem. A

114

,

12945

(

2010

).

https://doi.org/10.1021/jp108255k

Google Scholar

Crossref

PubMed

17.

E.

Kamarchik

,

Y.

Wang

, and

J. M.

Bowman

,

J. Chem. Phys.

134

,

114311

(

2011

).

https://doi.org/10.1063/1.3567186

Google Scholar

Crossref

PubMed

18.

Y.

Wang

and

J. M.

Bowman

,

J. Chem. Phys.

134

,

154510

(

2011

).

https://doi.org/10.1063/1.3579995

Google Scholar

Crossref

PubMed

19.

Y.

Wang

and

J. M.

Bowman

,

J. Chem. Phys.

136

,

144113

(

2012

).

https://doi.org/10.1063/1.3700165

Google Scholar

Crossref

PubMed

20.

P.

Huo

,

T. F.

Miller

, and

D. F.

Coker

,

J. Chem. Phys.

139

,

151103

(

2013

).

https://doi.org/10.1063/1.4826163

Google Scholar

Crossref

PubMed

21.

J.

Beutier

,

D.

Borgis

,

R.

Vuilleumier

, and

S.

Bonella

,

J. Chem. Phys.

141

,

084102

(

2014

).

https://doi.org/10.1063/1.4892597

Google Scholar

Crossref

PubMed

22.

K. K. G.

Smith

,

J. A.

Poulsen

,

G.

Nyman

, and

P. J.

Rossky

,

J. Chem. Phys.

142

,

244112

(

2015

).

https://doi.org/10.1063/1.4922887

Google Scholar

Crossref

PubMed

23.

J.

Liu

,

Int. J. Quantum Chem.

115

,

657

(

2015

).

https://doi.org/10.1002/qua.24872

Google Scholar

Crossref

24.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

99

,

10070

(

1993

).

https://doi.org/10.1063/1.465512

Google Scholar

Crossref

25.

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

100

,

5093

(

1994

);

https://doi.org/10.1063/1.467175

Crossref

Google Scholar

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

100

,

5106

(

1994

);

Crossref

Google Scholar

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

101

,

6157

(

1994

);

Crossref

Google Scholar

J.

Cao

and

G. A.

Voth

,

J. Chem. Phys.

101

,

6168

(

1994

).

Crossref

Google Scholar

26.

S.

Jang

and

G. A.

Voth

,

J. Chem. Phys.

111

,

2357

(

1999

).

https://doi.org/10.1063/1.479514

Google Scholar

Crossref

27.

S.

Jang

and

G. A.

Voth

,

J. Chem. Phys.

111

,

2371

(

1999

).

https://doi.org/10.1063/1.479515

Google Scholar

Crossref

28.

G. A.

Voth

,

Adv. Chem. Phys.

93

,

135

(

2007

).

https://doi.org/10.1002/9780470141526.ch4

Google Scholar

Crossref

29.

I. R.

Craig

and

D. E.

Manolopoulos

,

J. Chem. Phys.

121

,

3368

(

2004

).

https://doi.org/10.1063/1.1777575

Google Scholar

Crossref

PubMed

30.

I. R.

Craig

and

D. E.

Manolopoulos

,

J. Chem. Phys.

122

,

084106

(

2005

).

https://doi.org/10.1063/1.1850093

Google Scholar

Crossref

31.

T. F.

Miller

and

D. E.

Manolopoulos

,

J. Chem. Phys.

122

,

184503

(

2005

).

https://doi.org/10.1063/1.1893956

Google Scholar

Crossref

PubMed

32.

S.

Habershon

,

D. E.

Manolopoulos

,

T. E.

Markland

, and

T. F.

Miller

,

Annu. Rev. Phys. Chem.

64

,

387

(

2013

).

https://doi.org/10.1146/annurev-physchem-040412-110122

Google Scholar

Crossref

PubMed

33.

M.

Shiga

and

A.

Nakayama

,

Chem. Phys. Lett.

451

,

175

(

2008

).

https://doi.org/10.1016/j.cplett.2007.11.091

Google Scholar

Crossref

34.

S.

Habershon

,

G. S.

Fanourgakis

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

129

,

074501

(

2008

).

https://doi.org/10.1063/1.2968555

Google Scholar

Crossref

PubMed

35.

A.

Witt

,

S. D.

Ivanov

,

M.

Shiga

,

H.

Forbert

, and

D.

Marx

,

J. Chem. Phys.

130

,

194510

(

2009

).

https://doi.org/10.1063/1.3125009

Google Scholar

Crossref

PubMed

36.

S. D.

Ivanov

,

A.

Witt

,

M.

Shiga

, and

D.

Marx

,

J. Chem. Phys.

132

,

031101

(

2010

).

https://doi.org/10.1063/1.3290958

Google Scholar

Crossref

PubMed

37.

S.

Jang

and

G. A.

Voth

,

J. Chem. Phys.

112

,

8747

(

2000

).

https://doi.org/10.1063/1.481490

Google Scholar

Crossref

38.

E.

Geva

,

Q.

Shi

, and

G. A.

Voth

,

J. Chem. Phys.

115

,

9209

(

2001

).

https://doi.org/10.1063/1.1412870

Google Scholar

Crossref

39.

I.

Navrotskaya

,

Q.

Shi

, and

E.

Geva

,

Isr. J. Chem.

42

,

225

(

2002

).

https://doi.org/10.1560/46PR-4LGF-EHNY-NDHL

Google Scholar

Crossref

40.

Q.

Shi

and

E.

Geva

,

J. Chem. Phys.

116

,

3223

(

2002

).

https://doi.org/10.1063/1.1445120

Google Scholar

Crossref

41.

E.

Geva

,

S.

Jang

, and

G. A.

Voth

, in

Encyclopedia of Materials Modeling: Vol. I, Fundamental Models and Methods

, edited by

S.

Yip

(

Springer

,

Berlin

,

2005

).

Google Scholar

42.

M. W.

Dzierlenga

,

D.

Antoniou

, and

S. D.

Schwartz

,

J. Phys. Chem. Lett.

6

,

1177

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b00346

Google Scholar

Crossref

PubMed

43.

R.

Collepardo-Guevara

,

Y. V.

Suleimanov

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

130

,

174713

(

2009

).

https://doi.org/10.1063/1.3127145

Google Scholar

Crossref

PubMed

44.

J. O.

Richardson

and

S. C.

Althorpe

,

J. Chem. Phys.

131

,

214106

(

2009

).

https://doi.org/10.1063/1.3267318

Google Scholar

Crossref

PubMed

45.

N.

Boekelheide

,

R.

Salomón-Ferrer

, and

T. F.

Miller

,

Proc. Natl. Acad. Sci. U. S. A.

108

,

16159

(

2011

).

https://doi.org/10.1073/pnas.1106397108

Google Scholar

Crossref

PubMed

46.

A. R.

Menzeleev

,

N.

Ananth

, and

T. F.

Miller

,

J. Chem. Phys.

135

,

074106

(

2011

).

https://doi.org/10.1063/1.3624766

Google Scholar

Crossref

PubMed

47.

Y. V.

Suleimanov

,

R.

Collepardo-Guevara

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

134

,

044131

(

2011

).

https://doi.org/10.1063/1.3533275

Google Scholar

Crossref

PubMed

48.

R.

Perez de Tudela

,

F. J.

Aoiz

,

Y. V.

Suleimanov

, and

D. E.

Manolopoulos

,

J. Phys. Chem. Lett.

3

,

493

(

2012

).

https://doi.org/10.1021/jz201702q

Google Scholar

Crossref

PubMed

49.

J. S.

Kretchmer

and

T. F.

Miller

,

J. Chem. Phys.

138

,

134109

(

2013

).

https://doi.org/10.1063/1.4797462

Google Scholar

Crossref

PubMed

50.

Y.

Li

,

Y. V.

Suleimanov

,

M.

Yang

,

W. H.

Green

, and

H.

Guo

,

J. Phys. Chem. Lett.

4

,

48

(

2013

).

https://doi.org/10.1021/jz3019513

Google Scholar

Crossref

PubMed

51.

Y. V.

Suleimanov

,

W. J.

Kong

,

H.

Guo

, and

W. H.

Green

,

J. Chem. Phys.

141

,

244103

(

2014

).

https://doi.org/10.1063/1.4904080

Google Scholar

Crossref

PubMed

52.

Q.

Meng

,

J.

Chen

, and

D. H.

Zhang

,

J. Phys. Chem.

143

,

101102

(

2015

).

https://doi.org/10.1063/1.4930860

Google Scholar

Crossref

53.

K. M.

Hickson

,

J.-C.

Loison

,

H.

Guo

, and

Y. V.

Suleimanov

,

J. Phys. Chem. Lett.

6

,

4194

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b02060

Google Scholar

Crossref

PubMed

54.

J. S.

Kretchmer

and

T. F.

Miller

III,

Inorg. Chem.

55

,

1022

(

2015

).

https://doi.org/10.1021/acs.inorgchem.5b01821

Google Scholar

Crossref

PubMed

55.

Y. V.

Suleimanov

and

J.

Espinosa-Garcia

,

J. Chem. Phys. B

120

,

1418

(

2016

).

https://doi.org/10.1021/acs.jpcb.5b02103

Google Scholar

Crossref

56.

D.

Kang

,

H.

Sun

,

J.

Dai

,

W.

Chen

,

Z.

Zhao

,

Y.

Hou

,

J.

Zeng

, and

J.

Yuan

,

Sci. Rep.

4

,

5484

(

2014

).

https://doi.org/10.1038/srep05484

Google Scholar

Crossref

PubMed

57.

E.

Guarini

,

M.

Neumann

,

U.

Bafile

,

M.

Celli

,

D.

Colognesi

,

E.

Farhi

, and

Y.

Calzavara

,

Phys. Rev. B

92

,

104303

(

2015

).

https://doi.org/10.1103/PhysRevB.92.104303

Google Scholar

Crossref

58.

R.

Biswas

,

Y.-L. S.

Tse

,

A.

Tokmakoff

, and

G. A.

Voth

,

J. Phys. Chem. B

120

,

1793

(

2016

).

https://doi.org/10.1021/acs.jpcb.5b09466

Google Scholar

Crossref

PubMed

59.

T. F.

Miller

and

D. E.

Manolopoulos

,

J. Chem. Phys.

123

,

154504

(

2005

).

https://doi.org/10.1063/1.2074967

Google Scholar

Crossref

PubMed

60.

T. E.

Markland

,

S.

Habershon

, and

D. E.

Manolopoulos

,

J. Chem. Phys.

128

,

194506

(

2008

).

https://doi.org/10.1063/1.2925792

Google Scholar

Crossref

PubMed

61.

Y. V.

Suleimanov

,

J. Phys. Chem. C

116

,

11141

(

2012

).

https://doi.org/10.1021/jp302453z

Google Scholar

Crossref

62.

F.

Paesani

,

S. S.

Xantheas

, and

G. A.

Voth

,

J. Phys. Chem. B

113

,

13118

(

2009

).

https://doi.org/10.1021/jp907648y

Google Scholar

Crossref

PubMed

63.

F.

Paesani

and

G. A.

Voth

,

J. Phys. Chem.

132

,

014105

(

2010

).

https://doi.org/10.1063/1.3291212

Google Scholar

Crossref

64.

G. R.

Medders

and

F.

Paesani

,

J. Chem. Theory Comput.

11

,

1145

(

2015

).

https://doi.org/10.1021/ct501131j

Google Scholar

Crossref

PubMed

65.

M.

Rossi

,

H.

Liu

,

F.

Paesani

,

J.

Bowman

, and

M.

Ceriotti

,

J. Chem. Phys.

141

,

181101

(

2014

).

https://doi.org/10.1063/1.4901214

Google Scholar

Crossref

PubMed

66.

A. R.

Menzeleev

and

T. F.

Miller

,

J. Chem. Phys.

132

,

034106

(

2010

).

https://doi.org/10.1063/1.3292576

Google Scholar

Crossref

PubMed

67.

S.

Jang

,

Y.

Pak

, and

G. A.

Voth

,

J. Phys. Chem. A

103

,

10289

(

1999

).

https://doi.org/10.1021/jp9925177

Google Scholar

Crossref

68.

S.

Jang

,

J. Chem. Phys.

124

,

064107

(

2006

).

https://doi.org/10.1063/1.2162887

Google Scholar

Crossref

69.

T. J. H.

Hele

,

M. J.

Willatt

,

A.

Muolo

, and

S. C.

Althorpe

,

J. Chem. Phys.

142

,

134103

(

2015

).

https://doi.org/10.1063/1.4916311

Google Scholar

Crossref

PubMed

70.

T. J. H.

Hele

,

M. J.

Willatt

,

A.

Muolo

, and

S. C.

Althorpe

,

J. Chem. Phys.

142

,

191101

(

2015

).

https://doi.org/10.1063/1.4921234

Google Scholar

Crossref

PubMed

71.

S.

Jang

,

A. V.

Sinitskiy

, and

G. A.

Voth

,

J. Chem. Phys.

140

,

154103

(

2014

).

https://doi.org/10.1063/1.4870717

Google Scholar

Crossref

72.

S.

Jang

, preprint arXiv:1308.3805 (

2013

).

73.

T. J.

Hele

,

Mol. Phys.

114

,

1461

(

2016

).

https://doi.org/10.1080/00268976.2015.1136003

Google Scholar

Crossref

74.

B. J.

Braams

and

D. E.

Manolopoulos

,

J. Chem. Phys.

125

,

124105

(

2006

).

https://doi.org/10.1063/1.2357599

Google Scholar

Crossref

PubMed

75.

Note that the choice of centroid mean force F₁ in the CMD approximation for the case of a sudden vertical excitation gives lower-order agreement with Matsubara dynamics at short times than if we had chosen instead to use a centroid mean force as

F_{01} (\bar{Q}) = - \frac{\int d 𝐪 e^{- β_{N} U_{N}^{(0)} (𝐪)} \frac{\partial U_{N}^{(1)} (𝐪)}{\partial \bar{Q}} δ (\bar{q} (𝐪) - \bar{Q})}{\int d 𝐪 e^{- β_{N} U_{N}^{(0)} (𝐪)} δ (\bar{q} (𝐪) - \bar{Q})}

to generate the force on the centroid. Using F₀₁(⁠

\bar{𝑄}

⁠), the agreement with Matsubara dynamics can be shown to be correct to order t² for the position- and t¹ for the momentum-autocorrelation functions.

76.

D.

Chandler

, in

Liquids, Freezing and Glass Transition, PT 1, Les Houches Summer School Session

, edited by

J. P.

Hansen

and

D.

Levesque

, and

J.

Zinn-Justin

(

Elsevier Science

,

Amsterdam

,

1991

), Vol. 51, pp.

193

–

285

.

77.

We note that one can also scale the initial momentum. This corresponds to the Hamiltonians ${\hat{H}}^{(0)} = \frac{{\hat{p}}^{2}}{2 m σ_{p}} + V (\hat{q})$ and ${\hat{H}}^{(1)} = \frac{{\hat{p}}^{2}}{2 m} + V (\hat{q})$ and leads to ring-polymer Hamiltonians $H_{N}^{(0)} = \sum_{j = 1}^{N} \frac{p_{j}^{2}}{2 m σ_{p}} + \frac{m σ_{p}}{2 β_{N}^{2} {\bar{h}}^{2}} \sum_{j = 1}^{N} {(q_{j} - q_{j - 1})}^{2} + \sum_{j = 1}^{N} V (q_{j})$ and $H_{N}^{(1)} = \sum_{j = 1}^{N} \frac{p_{j}^{2}}{2 m σ_{p}} + \frac{m σ_{p}^{2}}{2 β_{N}^{2} {\bar{h}}^{2}} \sum_{j = 1}^{N} {(q_{j} - q_{j - 1})}^{2} + \sum_{j = 1}^{N} V (q_{j})$ ⁠. Again, this case can be easily combined with the two cases presented in the main text and exhibits similar formal properties.