Direct Computation of Influence Functional Coefficients from

Article pubs.acs.org/JCTC

Direct Computation of Influence Functional Coefficients from Numerical Correlation Functions Thomas C. Allen, Peter L. Walters, and Nancy Makri* Department of Chemistry, University of Illinois at Urbana−Champaign, 600 S. Goodwin Avenue, Urbana, Illinois 61801, United States ABSTRACT: Influence functional methods provide a powerful approach for simulating the dynamics of a system embedded in a harmonic bath, which may be parametrized to reflect a variety of environments and chemical processes. In this work, we develop a procedure for calculating the coefficients of the discretized influence functional using the classical approximation to the time correlation function, which is usually available through molecular dynamics simulations. This procedure circumvents the calculation of a spectral density via Fourier inversion of the correlation function. When the correlation function is available with high precision, we find that the direct procedure yields results just as accurate as those obtained using the spectral density expressions. However, when statistical noise is present, the direct procedure produces more accurate results. The direct procedure is efficient and easy to implement.

calculation,4 which destroy the Markovian character of the time-dependent Schrödinger equation in full space. Noting that the magnitude of these couplings should ultimately decay to zero when time points are well separated, our group has developed an iterative algorithm5 that employs a multi-time propagator (i.e., a tensor). This approach involves a physically motivated splitting of the evolution operator, valid over sizable time steps,6 along with relatively smooth energy-filtered system propagators 7 in a potential-optimized discrete variable representation of the influence functional,8 a collection of time-discretized system paths, and a set of two-time coefficients, ηk′k″, which determine the influence functional.9 The iterative quasi-adiabatic propagator path integral methodology (i-QuAPI) converges to the numerically exact dynamics once the memory tensor has incorporated enough time segments to span the bath-induced decorrelation time. The iQuAPI methodology is applicable to baths described by arbitrary spectral densities and most efficient when the bath consists mostly of high-frequency modes. Further acceleration of the i-QuAPI methodology is possible through path filtering techniques,10 propagator renormalization11 and the blip decomposition.12 The iterative path integral methodology has been extended to Fermionic baths13 and shown to be generalizable to arbitrary condensed-phase environments, provided that the influence functional is available (analytically or numerically).14 Quantum-semiclassical15 and quantumclassical16 formulations of the path integral exploit such iterative decompositions. In the original work,5b expressions for the discretized influence functional coefficients were obtained in terms of the

I. INTRODUCTION Owing to the difficulty of solving the Schrödinger equation for large systems, a number of approaches have been developed over the years to tackle problems involving the quantum dynamics of condensed-phase systems. A common initial step is to divide the model under study into a relatively compact set of “system” degrees of freedom as well as a more extended set of “bath” coordinates.1 A number of powerful analytical and computational tools are available for studying the dynamics of a system interacting with a bath of harmonic oscillators. For this reason, harmonic bath models have been formulated to describe the generic effects of molecular and condensedphase environments on the dynamics of barrier crossing, relaxation, tunneling, entanglement, or laser control. Moreover, the harmonic bath description of a system’s environment is intimately connected to the linear response approximation;2 so in many cases, these models are capable of accurately describing the quantitative physics underlying specific condensed-phase processes. Within this framework, one particularly appealing approach to evaluating the dynamics of system-bath Hamiltonians employs the path integral formalism.3 In this context, Feynman and Vernon have shown4 that the effects of a harmonic bath on the dynamics of the system may be integrated out and expressed in reduced form as an influence functional. This functional can be obtained via computational methods or parametrized by a handful of experimentally observable macroscopic properties of the bath, providing an extremely compact description of the changes induced in the quantum system’s dynamics by the presence of additional environmental degrees of freedom. However, as a consequence of this reduced description, couplings are introduced between all time points in the © 2016 American Chemical Society

Received: April 19, 2016 Published: July 21, 2016 4169

DOI: 10.1021/acs.jctc.6b00390 J. Chem. Theory Comput. 2016, 12, 4169−4177

Article

Journal of Chemical Theory and Computation

Here, s± (t) are the forward and backward paths of the system degrees of freedom, and α (t) is the correlation function of the bath,

bath spectral density;1c these expressions are particularly useful for theoretical work where various analytical models for the spectral density are commonly employed and the effect of different categories of densities on the system dynamics is of interest in its own right. Simulation of a particular chemical or biological process requires first extracting the spectral density from the environment’s correlation function. In this paper, we present an alternative derivation of the coefficients for the discretized influence functional directly from the correlation function of the environment, thus bypassing the numerical construction of a spectral density. This approach is particularly relevant to much current work, where the bath is frequently derived from systems characterized by experimental measurements or molecular dynamics simulations in which correlation functions are the primary objects of interest.17 In such situations, the method presented in this paper has several advantages. It is direct and easy to implement, avoiding subtle issues tied to boundary conditions and signal processing choices necessitated when actually performing a numerical Fourier inversion of the correlation function. Further, in common situations where statistical noise is present in the correlation function, the direct determination of the influence functional coefficients leads to increased accuracy. Section II reviews the discretized influence functional and presents the expression for direct evaluation of the coefficients in the time domain. Several applications on tunneling dynamics in model dissipative baths and baths that describe electron or proton transfer reactions in solution are presented in Section III. Finally, some concluding remarks are given in Section IV.

α (t ) =

Ĥ sb =

∑ j

α (t ) =

∫0

dt ′

∫0

t′

⎧ 1 F = exp⎨− ⎩ ℏ

j

δ(ω − ωj) (2.6)

⎫ (sk+′ − sk−′)(ηk ′ k ″sk+″ − ηk*′ k ″sk−″)⎬ ⎭ k ″= 0 k′

⎪

∑ ∑

⎪

k ′= 0

⎪

(2.7)

where the coefficients are given by the integrals ηk ′ k ″ =

(k′+ 1 )Δt 2 (k′− 1 )Δt 2

∫

dt ′

(k′′+ 1 )Δt 2 (k′′− 1 )Δt 2

∫

dt ″α(t ′ − t ″), k′ > k″ (2.8)

ηk 0 =

(k + 1 )Δt 2 (k − 1 )Δt 2

∫

dt ′

∫0

1/2Δt

dt ″α(t ′ − t ″), k ≠ 0, N (2.9)

(2.1) (k + 1 )Δt 2 1 (k − )Δt 2

N Δt

ηNk =

∫(N− 1 )Δt dt′ ∫ 2

dt ″α(t ′ − t ″), k ≠ 0, N (2.10)

(k + 1 )Δt

ηkk =

∫(k− 1 )2Δt

+ iΔt ∑ j

(2.2)

η00 = ηNN =

∫0

t′

dt ′

∫(k− 1 )Δt dt″α(t′ − t″) 2

c j2 2mjωj2 1/2Δt

dt ′

(2.11)

∫0

t′

dt ″α(t ′ − t ″)

c j2 i + Δt ∑ 2 2 j 2mj ωj

dt ″(s+(t ′) − s−(t ′))

dt ′ ∑

mjωj

N

⎪

⎞ (α(t ′ − t ″)s (t ″) − α*(t ′ − t ″)s (t ″))⎟ ⎠

∫0

j

c j2

Splitting the time t into path integral steps Δt = t/N gives rise to the discretized form of the influence functional, where the time integrals are replaced by sums. With the QuAPI discretization of the path integral,6 the influence functional has the form

−

t

∑

2

+

⎛ i × exp⎜⎜ − ⎝ ℏ

π 2

J(ω) ≡

the influence functional is given by Feynman and Vernon’s form,4 augmented by a phase associated with the counterterms: t

(2.5)

The latter can also be written in terms of the spectral density function,1b

Note that the “counterterms”1c are included with the bath Hamiltonian, ensuring that stable states (e.g., donor and acceptor in the case of charge transfer) are separated by the energy difference specified by the Hamiltonian of the bare system.6 If the bath is initially described by the thermal density matrix, e−βHb/Tr e−βHb, of the isolated bath Hamiltonian,

⎛ 1 F = exp⎜ − ⎝ ℏ

c j2 ⎡ ⎤ ⎛1 ⎞ ⎢⎣coth⎝⎜ ℏωjβ ⎠⎟cos ωjt − i sin ωjt ⎥⎦ 2mjωj 2

∑ j

pj2̂

⎞ ⎛ p2̂ 1 j Ĥ b = ∑ ⎜⎜ + mjωj2xĵ2⎟⎟ 2mj 2 j ⎝ ⎠

(2.4)

where f j = cjxj is the force on the system by the bath. Equation 2.4 is given by the expression

II. METHODOLOGY For a system described by the coordinate s, the Hamiltonian describing the harmonic bath and its bilinear interaction with the system in the quasi-adiabatic splitting of the propagator6 has the form 2 ⎛ cjs ̂ ⎞ 1 2⎜ ⎟ + mjωj ⎜xĵ − 2mj 2 mjωj2 ⎟⎠ ⎝

−βĤ 1 Tr(e bf ̂ (t )·f ̂ (0)) ̂ ℏ Tr(e−βHb)

N Δt

ηN 0 =

⎞ 2 2 ⎟ + − (s (t ′) − s (t ′) )⎟ 2mjωj2 ⎠

∫(N− 1 )Δt dt′ ∫0 2

c j2

(2.12)

1/2Δt

dt ″α(t ′ − t ″) (2.13)

Makri and Makarov obtained expressions for these coefficients in terms of frequency integrals involving the spectral density function.5a

(2.3) 4170


Article


functional coefficients requires the availability of J (ω) on a fine frequency grid, which necessitates knowledge of the correlation function over a rather long time interval. However, numerical errors in MD calculations tend to grow over time, making longtime data costly and less reliable. Thus, one must generally interpolate the sparse spectral density that results from the transformation. In addition, the global nature of the Fourier transform, eq 2.17, effectively spreads the noise found in the long-time tail of the correlation function over the entire frequency domain of the obtained spectral density. Thus, we consider here the alternative, direct evaluation of the influence functional coefficients from the correlation function. Equations 2.8−2.12 can be employed for this purpose. However, use of the real-valued classical correlation function in these equations would yield real-valued influence functional coefficients, which would not capture quantum decoherence contributions to the dynamics.25 Below, we assume that the classical correlation function for the effective harmonic bath on which the system’s environment is being mapped, eq 2.16, is available from Boltzmann-weighted molecular dynamics simulations. Taking the time derivative of eq 2.16, we obtain

Working in frequency space is especially attractive if the spectral density is available directly without necessitating a numerical transformation; this situation occurs primarily in model studies where general features of the bath and its impact on the system are of interest, rather than specific results for particular chemical processes. However, in simulations of a particular chemical or biological system, the available information often is the time correlation function. In those situations one would first obtain the spectral density through numerical Fourier transformation of the correlation function, ⎛1 ⎞ J(ω) = 2tanh⎜ ℏωβ ⎟ ⎝2 ⎠

∫0

∞

Re α(t )cos ωt dt

(2.14)

Note that eq 2.14 requires only the real part of the correlation function. Molecular dynamics (MD) simulations yield the classical limit of this quantity, αcl(t ) = lim α(t )

(2.15)

ℏ→ 0

which in the case of a harmonic bath is given by αcl(t ) = lim α(t ) = ℏωjβ → 0

= (πβ ℏ)−1

∑ j

∞

∫−∞

c j2 ⎛ 1 ⎞−1 ⎜ ℏω β ⎟ cos ω t j j ⎠ 2mjωj ⎝ 2

J(ω) cos ωt dω ω2

2 d αcl(t ) = − dt ℏβ

(2.16)

∫0

2mjωj

sin(ωjt ) =

2 Im α(t ) ℏβ

Equation 2.18 gives the imaginary part of the correlation function in terms of its real part. Using this expression, it is now possible to rewrite the integrals which determine the discretized influence functional coefficients. We first consider the off-diagonal coefficients, eq 2.8, which can now be written in terms of the real-valued, classical bath correlation function alone:

∞

Re αcl(t )cos ωt dt

j

c j2

(2.18)

Inversion of this leads to the following expression for the spectral density in terms of the classical correlation function: J(ω) = ℏωβ

∑

(2.17)

If the medium of interest exhibits truly classical behavior (ℏωjβ ≪ 1 for strongly coupled degrees of freedom), one may use αcl in eq 2.17 to extract the spectral density. On the other hand, high frequency molecular vibrations often cause pronounced deviations from classical behavior. Force fields frequently employed in MD simulations can account for such effects to some extent. However, if the classical trajectories are integrated subject to forces obtained from ab initio potential energy surfaces, some correction of the classical correlation function may be necessary. Quantum effects in the correlation function of condensed-phase systems arise primarily from zero-point energy, and one can employ one of several available quantum correction factors18 to include such effects in an approximate fashion. Alternatively, quasiclassical methods, such as linearized semiclassical initial value or linearized path integral approximations,19 or forward−backward semiclassical dynamics20 (FBSD), may be used instead of MD. These methods, which employ classical trajectories sampled from a quantized phase space distribution (either the Wigner21 or the Husimi22 coherent state transform of the density operator) offer excellent alternatives to quantum dynamical calculations in such cases. The only challenge of such calculations is the determination of the quantized density distribution. Fortunately, the coherent state density is easily evaluated by imaginary time path integral techniques,23 and a variety of approximations are available for generating the Wigner phase space function.24 Regardless of the particular approximation used to obtain the bath correlation function, extraction of the spectral density from eqs 2.14 or 2.17 requires care. Most notably, these integrals require knowledge of the correlation function over times at least as long as the decay time of this function. Further, evaluation of the integrals required to compute the influence

Re ηk ′ k ″ =

(k′+ 1 )Δt 2 1 ′− (k )Δt 2

∫

dt ′

(k′′+ 1 )Δt 2 ′′ 1 (k − )Δt 2

∫

dt ″αcl(t ′ − t ″) (2.19)

Im ηk ′ k ″ =

1 ℏβ 2

(k′+ 1 )Δt 2 (k′− 1 )Δt 2

∫

dt ′

(k′′+ 1 )Δt 2 (k′′− 1 )Δt 2

∫

dt ″

d d (t ′ − t ″ )

αcl(t ′ − t ″)

(2.20)

The real components of the off-diagonal coefficients, eq 2.19, may be obtained by straightforward numerical quadrature. On the other hand, the imaginary components can be simplified. Equation 2.20 can be written as a single integral, Im ηk ′ k ″ =

1 ℏβ 2

(k′+ 1 )Δt

∫(k′− 1 )2Δt 2

⎡ ⎛ ⎛ 1 ⎞⎞ ⎢αcl ⎜t ′ − ⎝⎜k″Δt − Δt ⎠⎟⎟ ⎠ ⎣ ⎝ 2

⎛ ⎛ 1 ⎞⎞⎤ − αcl ⎜t ′ − ⎜k″Δt + Δt ⎟⎟⎥dt ′ ⎝ ⎝ 2 ⎠⎠⎦

(2.21)

and a change of variables leads to the expression Im ηk ′ k ″ =

1 ⎡ ℏβ ⎢ 2 ⎣

(k′− k′′+ 1)Δt

∫(k′−k )Δt ′′

(k′− k′′)Δt

− 4171

αcl(t )dt ⎤

∫(k′−k −1)Δt αcl(t )dt ⎥⎦ ′′

(2.22)


Article


trajectories must be available over long time intervals in order to obtain the spectral density on closely spaced points. Further, experience in our group indicates that one must exercise special care in the treatment of boundary points while performing the discrete cosine transform to avoid numerical error in the reorganization energy. On the other hand, the direct procedure which is based on integration of α(t), only requires short time data, at most up to the time where the correlation function itself becomes negligible, allowing the use of shorter molecular dynamics trajectories to map the given chemical environment on a harmonic bath. This is so because the influence functional coefficients are translationally invariant; i.e., ηk′k″ = ηk′+m k″+m, implying that only those with k″ = 0,1 need to be determined numerically. Further, the largest and most significant coefficients are those with small k′ (and k″), which are obtained from the early time values of the correlation function. The required values of k′ − k″ ≃ k′ are determined by the effective memory length ΔkmaxΔt, which tends to be shorter than the time span of the correlation function. Even if this is not the case, the noise in the tail of the correlation function will affect the accuracy of only those coefficients with k′ − k″ ≫ 1, which tend to be rather small and not as crucial in determining the value of the influence functional. Last, the time integration in the proposed scheme further smooths out the integrand, partially canceling numerical noise inherent in the MD data.

Thus, the imaginary components of the influence functional coefficients can be obtained by evaluating two single integrals of the classical correlation function. The above procedure is straightforwardly extended to all off-diagonal coefficients, leading to the expressions Im ηk 0 =

1 ⎡ ℏβ ⎢ 2 ⎢⎣

(k + 1 )Δt

⎤

k Δt

∫(k− 1 )2Δt

αcl(t )dt −

2

∫(k−1)Δt αcl(t )dt ⎥⎥⎦ (2.23)

Im ηNk =

ℏβ ⎡ ⎢ 2 ⎣

(N − k + 1/2)Δt

∫(N−k)Δt

αcl(t )dt

(N − k − 1/2)Δt

−

Im ηN 0 =

∫(N−k−1)Δt

ℏβ ⎡ ⎢ 2 ⎣

⎤ αcl(t )dt ⎥ ⎦

N Δt

(2.24) (N − 1/2)Δt

∫(N−1/2)Δt αcl(t )dt − ∫(N−1)Δt

⎤ αcl(t )dt ⎥ ⎦ (2.25)

Next, we turn to the diagonal coefficients. In this case the imaginary components include the discretized counterterms, i.e., Im ηkk =

1 ℏβ 2

(k + 1 )Δt 2 (k − 1 )Δt 2

∫

t′

dt ′

∫(k− 1 )Δt dt″ d(t′ d− t″) 2

αcl(t ′ − t ″) + Δt ∑ j

III. NUMERICAL TESTS In this section, we examine the numerical performance of our method in comparison with the procedure based on integrating J(ω) for a variety of systems, including some based on realistic molecular models where the correlation functions were obtained from molecular dynamics simulations. First, consider a process for which the system-bath interaction is described by the Ohmic spectral density,26 π J(ω) = ξℏωe−ω / ωc (3.1) 2

c j2 2mjωj2

(2.26)

Proceeding as before, this can be rewritten as 1 ⎡ Im ηkk = ℏβ ⎢ 2 ⎣

∫0

Δt

c j2 ⎤ αcl(t )dt − αcl(0)Δt ⎥ + Δt ∑ ⎦ 2mjωj2 j (2.27)

Noting that −1

αcl(0) = (ℏβ)

∑ j

c j2 mjωj2

where ξ is a measure of the system-bath coupling strength and the characteristic frequency ωc corresponds to the maximum of J(ω). We assume that the spectral density for this bath is not directly available, forcing one to resort to classical molecular dynamics techniques to capture the correlation of the bath. In that case one would obtain the classical limit of the correlation function, which for this Ohmic spectral density has the Lorentzian form,

(2.28)

eq 2.27 takes the simple form Im ηkk =

1 ℏβ 2

∫0

Δt

αcl(t )dt

(2.29)

Similar cancellation occurs in the other diagonal coefficients. One finds Im η00 = Im ηNN =

1 ℏβ 2

∫0

αcl(t ) =

1/2Δt

αcl(t )dt

(2.30)

ξωc β(1 + ωc2t 2)

(3.2)

One could obtain the spectral density by evaluating eq 2.17 and then calculate the influence functional coefficients according to the expressions in ref 5b or compute the influence functional coefficients directly according to the prescription given in the previous section. Below, we present numerical results and compare the two procedures. Table 1 compares the values of the influence functional coefficients obtained via the direct procedure described in Section II to those obtained from the exact Ohmic spectral density using the expressions in ref 5b (evaluated on a fine frequency grid) for a temperature that corresponds to ℏωc β = 0.15. At this relatively high temperature, the bath’s behavior is practically classical, thus approximating the correlation function by eq 3.2 should not result in significant error. Indeed, we find

The real components of the diagonal coefficients are again obtained via a two-dimensional quadrature. The main advantage of using the derived relations is that the discrete Fourier transform required to obtain the spectral density from the correlation function is avoided. The original form of the influence functional coefficients involves integrals of the spectral density multiplied by trigonometric functions, which require a fine frequency grid. In addition, the zero frequency point must be treated with care, as the integrand of the influence functional coefficients contains removable singularities or (in the case of the diagonal coefficients) a pole. However, because of the reciprocal Fourier relationship between frequency spacing and time length, classical 4172


Article

Journal of Chemical Theory and Computation Table 1. Comparison of ηk′k″ Values Obtained from Integration of Spectral Density (ref 8) versus Present Method of Direct Calculation from αcl(t) for an Ohmic Spectral Densitya Re ηk′k″ η00 η10 η11 η20 η21 η30 η31 η20,1 η30,1 a

because the imaginary part is independent of temperature in the case of a harmonic bath. Next, we investigate the effects of noise in the correlation function, which arises from experimental or numerical sources. Figure 1 shows the Ohmic correlation function for the parameters used in Table 1, with added Gaussian noise whose amplitude increases as a Wiener process27 (i.e., as the square root of time,28 chosen to reflect the inherently increasing noise in typical MD correlation functions because of difficulties sampling large time separations). We assume this correlation function spans the time interval shown in the figure and that it is available at 5000 equally spaced time points. According to the Fourier theorem, the fine spacing in the time domain provides information in the frequency domain over a very large interval, while the frequency spacing is limited by the time length of available correlation function data. Thus, the useful part of the spectral density shown in Figure 1 includes a small number of data points. The sparsity of the spectral density introduces numerical error in the integrals required to obtain the influence functional coefficients, which contain oscillatory factors and (in the case of the real part of ηkk) pole behavior. Further, the Fourier transformation spreads the noise in the long-time tail of the correlation function over all frequencies. As shown in Tables 2 and 3, this procedure leads to non-negligible error in the influence functional coefficients, which becomes large as k′ − k″ increases. By contrast, the direct procedure presented in Section II produces very accurate results, with errors that are smaller by one or more orders of magnitude. There are two reasons for this impressive gain in accuracy. First, the direct procedure relies on values of the correlation function around the time (k′ − k″)Δt. For small values of k′ − k″, the correlation function is large and relatively free of statistical error in this range, allowing accurate evaluation of the integrals involved. Second, the time integration tends to smooth out the noise, which occurs on a time scale finer than the path integral time step. Because of this smoothing, the coefficients obtained via the direct procedure remain accurate even for large values of k′ − k″, which fall in the noisy region of the correlation function. Table 3 shows that the relative error of very small influence functional coefficients eventually grows and fluctuates in magnitude, but the absolute error is in the fifth decimal place and thus completely insignificant.

Im ηk′k″

from J(ω)

from αcl(t)

Δ (%)

from J(ω)

from αcl(t)

Δ (%)

0.02936 0.11492 0.11693 0.10580 0.22601 0.09258 0.20529 0.01713 0.00767

0.02925 0.11454 0.11651 0.10558 0.22532 0.09253 0.20494 0.01714 0.00767

0.37 0.33 0.36 0.21 0.31 0.05 0.17 0.06 0.00

0.04671 −0.00233 0.09263 −0.00463 −0.00596 −0.00560 −0.00990 −0.00067 −0.00020

0.04674 −0.00233 0.09267 −0.00463 −0.00596 −0.00560 −0.00990 −0.00067 −0.00020

0.06 0.00 0.04 0.00 0.00 0.00 0.00 0.00 0.00

Parameters: ℏωc β = 0.15, ξ = 1,Δt = 1.25 ℏβ

that the procedure described in this paper leads to values that match those obtained from the spectral density to three significant figures. We note that if the temperature is lowered, eq 3.2 will cease to provide a faithful approximation of the bath’s correlation function, and its use will introduce systematic error in the real part of the influence functional coefficients. This error is inherent in the classical approximation of the correlation function and not a flaw of the procedure presented in this paper. If this classical-limit correlation function were the only information available, the original procedure5a would first use it to obtain a spectral density through eq 2.14. This procedure would lead to an approximate spectral density, which would differ from eq 3.1 at low temperatures. The errors introduced through the use of a classical-limit correlation function can be minimized by adding some quantum effects in the trajectory calculation. This can be achieved either by using one of various quantum correction factors18 to arrive at an improved approximation to the correlation function or (more accurately) through the use of quasiclassical trajectory methods, which sample trajectories from a quantized (Wigner or coherent state transformed) phase space distribution. Note that the imaginary part of the influence functional coefficients is unaffected by errors in the classical approximation of the correlation function

Figure 1. Effects of Gaussian noise in the classical-limit correlation function (left panel) on the spectral density (right panel) for the Ohmic model with ℏωc β = 0.15, ξ = 1, Δt = 1.25 ℏβ. 4173


Article


Table 2. Comparison of the Real Part of the Influence Functional Coefficients, Re ηk′k″, Calculated from Sources with Added Gaussian Noise to the Exact Resultsa η00 η10 η11 η20 η21 η30 η31 η20,1 η30,1 a

from noise-free J(ω)

from noisy J(ω)

Δ (J, %)

from noisy αcl(t)

Δ (α, %)

0.02936 0.11492 0.11693 0.10580 0.22601 0.09258 0.20529 0.01713 0.00767

0.02817 0.11012 0.11218 0.10099 0.21641 0.08779 0.19566 0.00418 −0.01525

4.05 4.18 4.06 4.54 4.25 5.17 4.69 75.60 298.83

0.02926 0.11454 0.11652 0.10555 0.22531 0.09254 0.20488 0.01705 0.00776

0.34 0.33 0.35 0.23 0.31 0.04 0.20 0.47 0.13

Parameters are the same as those in Table 1.

Table 3. Comparison of the Imaginary Part of the Influence Functional Coefficients, Im ηk′k″, Calculated from Sources with Added Gaussian Noise to the Exact Resultsa η00 η10 η11 η20 η21 η30 η31 η20,1 η30,1 a

from noise-free J(ω)

from noisy J(ω)

Δ (J, %)

from noisy αcl(t)

Δ (α, %)

0.04671 −0.00233 0.09263 −0.00463 −0.00596 −0.00560 −0.00990 −0.00067 −0.00020

0.04047 −0.00235 0.08011 −0.00465 −0.00598 −0.00560 −0.00990 −0.00083 −0.00119

13.36 0.86 13.52 0.43 0.34 0.0 0.0 23.88 495.00

0.04674 −0.00234 0.09268 −0.00464 −0.00598 −0.00558 −0.00989 −0.00067 −0.00029

0.06 0.43 0.05 0.22 0.34 0.36 0.10 0.00 45.00

Parameters are the same as those in Table 1.

Figure 2. Correlation function obtained from molecular dynamics simulations of a proton transfer complex in liquid methyl chloride at 247 K, and the corresponding spectral density.

correlation function is shown in Figure 2 for the initial 3 ps. This correlation function was transformed to yield the spectral density. As in the Ohmic model, most frequency points obtained this way were outside of the useful range of the spectral density. The interpolated spectral density, which was used in the evaluation of the influence functional coefficients, is also shown in Figure 2. We note that the characteristic frequencies of this model system are rather low, validating the use of a classical trajectory approximation to obtain the correlation function. The influence functional coefficients obtained with the two procedures are compared in Table 4. For consistency, we report the percent error with respect to the frequency domain result, although (based on the behaviors discussed in the model calculations) we believe that the values obtained with the direct procedure in terms of the correlation function are more accurate. Differences, and thus the error of

Next, we apply the procedure described in Section II to correlation functions obtained from actual molecular dynamics simulations of two reactions in solution. The first model we consider is based on a phenol-amine proton transfer model in liquid methyl chloride with the parametrization employed in ref 29b. Molecular dynamics simulations were performed at T = 247 K using the package LAMMPS. A cubic simulation box of side length 28 Å containing 255 solvent molecules with periodic boundary conditions was used. A Nosé−Hoover thermostat was applied and bonds were held fixed using the SHAKE algorithm.30 Classical trajectories were run with a time step of 1 fs and the correlation function corresponding to the reactant-product energy gap was obtained up to t = 40 ps. The noise was kept minimal by time-averaging the results over 1.5 ns and also with respect to 15 different initial conditions. The obtained classical 4174


Article

Journal of Chemical Theory and Computation Table 4. Comparison of results obtained for a model proton transfer system from the two methods at T = 247 K with Δt = 12.1 fs Im ηk′k″

Re ηk′k″ η00 η10 η11 η20 η21 η30 η31 η20,1 η30,1

Table 5. Comparison of results obtained for an electron transfer model based on ferrocene in benzene at T = 300 K with Δt = 19.3 fs Im ηk′k″

Re ηk′k″

J(ω)

α(t)

Δ (%)

J(ω)

α(t)

Δ (%)

0.14305 0.57051 0.57186 0.56386 1.13835 0.55215 1.12247 0.34136 0.17386

0.14269 0.56903 0.57040 0.56228 1.13530 0.55046 1.11930 0.32873 0.15983

0.25 0.26 0.26 0.28 0.27 0.31 0.28 3.70 8.07

0.72734 −0.00504 1.45299 −0.01157 −0.01339 −0.01765 −0.02625 −0.03796 −0.01068

0.72952 −0.00525 1.45730 −0.01195 −0.01389 −0.01816 −0.02708 −0.03854 −0.01002

0.30 4.17 0.30 3.28 3.73 2.90 3.16 1.53 6.18

η00 η10 η11 η20 η21 η30 η31

J(ω)

α(t)

Δ (%)

J(ω)

α(t)

Δ (%)

0.9654 3.4166 3.7073 3.4893 6.7326 3.5205 7.1364

0.9285 3.4877 3.6413 3.5061 6.9092 3.5074 7.0885

3.82 2.08 1.78 0.48 2.62 0.37 0.67

2.4013 −0.1906 4.6555 0.1534 −0.1613 −0.1792 0.1617

2.4220 −0.1905 4.6969 0.1534 −0.1613 −0.1792 0.1617

0.86 0.05 0.89 0.00 0.00 0.00 0.00

correlation function, i.e., the quantity obtained from molecular dynamics simulations. The first advantage of this direct approach is the simplification attained by bypassing the numerical evaluation of the cosine transform necessary to compute the spectral density, a procedure that requires special care of boundary points for accuracy. The second advantage of the direct evaluation in the time domain is increased accuracy if the correlation function contains numerical noise. Such noise, which often affects primarily the tail of the correlation function, is spread through the cosine transform over the entire frequency range of the spectral density. On the other hand, the most significant influence functional coefficients are based on information available in the earlier, less noisy part of the correlation function. In addition, the time integration involved in the proposed method tends to smooth out such noisy behavior, minimizing error. Even though the classical correlation function is purely real, the procedure we described does not neglect the imaginary part of the influence functional coefficients. This was possible by extracting the imaginary part from the available real part within the harmonic bath model. By retaining this imaginary part, the influence functional constructed this way captures the delicate quantum mechanical contributions to decoherence via spontaneous phonon emissions, which are critical to attain detailed balance. The derived expressions are easy to implement via numerical quadrature and allow immediate application to map realistic molecular and atomistic systems on simpler harmonic bath models.

the spectral density procedure, are seen to be quite small, owing to the high precision of the computed correlation function, which involved extensive averaging. Again, these differences become appreciable for large values of k′ − k″. Based on the conclusions from the model baths, the results of the direct time domain method should be more accurate. The second case we consider is a harmonic model constructed from atomistic data on the electron transfer reaction between ferrocene and ferrocenium in a benzene solvent at 300 K and pressure of 1 Atm. Molecular dynamics trajectories were obtained in the NPT ensemble using the package NAMD 2.931 with a time step of 2 fs for an initial simulation box of dimensions 42 × 42 × 47 Å, which included 459 benzene molecules with periodic boundary conditions. The CH bonds were held fixed. The correlation function was obtained up to 80 ps with results time-averaged over 1.2 ns and also with respect to 30 trajectory initial conditions. The correlation function, which in this case is the energy gap correlation function of the solvent, is shown in Figure 3, along with the spectral density obtained through the cosine transform with subsequent interpolation. The results for this model are summarized in Table 5. Once again, the high quality of the correlation function data led to close agreement in the influence functional coefficients obtained with the two procedures.

IV. CONCLUDING REMARKS In this work, we have described a procedure for calculating the discretized influence functional coefficients directly from information available through the classical limit of the

Figure 3. Correlation function obtained from molecular dynamics simulations of the ferrocene-ferrocenium charge transfer in benzene at T = 300 K, and the corresponding spectral density. 4175


Article


■

time dissipative charge transport dynamics. Mol. Phys. 2012, 110, 1967−1975. (11) Makri, N. Path integral renormalization for quantum dissipative dynamics with multiple timescales. Mol. Phys. 2012, 110, 1001−1007. (12) Makri, N. Blip decomposition of the path integral: Exponential acceleration of real-time calculations for quantum dissipative systems. J. Chem. Phys. 2014, 141, 134117. (13) (a) Segal, D.; Millis, A. J.; Reichman, D. R. Numerically exact path integral simulation of nonequilibrium quantum transport and dissipation. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 205323. (b) Segal, D.; Millis, A. J.; Reichman, D. R. Nonequilibrium transport in quantum impurity models: exact path integral simulations. Phys. Chem. Chem. Phys. 2011, 13, 14378−14386. (c) Simine, L.; Segal, D. Path integral simulations with fermionic and bosonic reservoirs: Transport and dissipation in molecular electronic junctions. J. Chem. Phys. 2013, 138, 214111. (14) Makri, N. Iterative evaluation of the path integral for a system coupled to an anharmonic bath. J. Chem. Phys. 1999, 111, 6164−6167. (15) Makri, N.; Thompson, K. Semiclassical influence functionals for quantum systems in anharmonic environments. Chem. Phys. Lett. 1998, 291, 101−109. (16) (a) Lambert, R.; Makri, N. Quantum-classical path integral: Classical memory and weak quantum nonlocality. J. Chem. Phys. 2012, 137, 22A552. (b) Lambert, R.; Makri, N. Quantum-classical path integral: Numerical formulation. J. Chem. Phys. 2012, 137, 22A553. (c) Makri, N. Quantum-classical path integral: A rigorous approach to condensed phase dynamics. Int. J. Quantum Chem. 2015, 115, 1209− 1214. (17) Frenkel, D.; Smit, B. Understanding Molecular Simulation from Algorithms to Applications; Elsevier: San Diego, 2002. (18) (a) Schofield, P. Phys. Rev. Lett. 1960, 4, 239. (b) Egelstaff, P. A. Adv. Phys. 1962, 11, 203. (c) Berne, B. J.; Harp, G. D. On the calculation of time correlation functions. Adv. Chem. Phys. 1970, 17, 63. (d) Oxtoby, D. W. Vibrational population relaxation in liquids. Adv. Chem. Phys. 1981, 47, 487−519. (e) Bader, J. S.; Berne, B. J. Quantum and classical relaxation rates from classical simulations. J. Chem. Phys. 1994, 100, 8359. (f) Egorov, S. A.; Everitt, K. F.; Skinner, J. L. Quantum dynamics and vibrational relaxation. J. Phys. Chem. A 1999, 103, 9494−9499. (g) Lawrence, C. P.; Nakayama, A.; Makri, N.; Skinner, J. L. Quantum dynamics of simple fluids. J. Chem. Phys. 2004, 120, 6621−6624. (19) (a) Wang, H.; Sun, X.; Miller, W. H. Semiclassical approximations for the calculation of thermal rate constants for chemical reactions in complex molecular systems. J. Chem. Phys. 1998, 108, 9726−9736. (b) Miller, W. H. The semiclassical initial value representation: a potentially practical way for adding quantum effects to classical molecular dynamics simulations. J. Phys. Chem. A 2001, 105, 2942−2955. (c) Poulsen, J. A.; Nyman, G.; Rossky, P. J. Practical evaluation of condensed phase quantum correlation functions: A Feynman–Kleinert variational linearized path integral method. J. Chem. Phys. 2003, 119 (23), 12179−12193. (20) (a) Shao, J.; Makri, N. Forward-backward semiclassical dynamics without prefactors. J. Phys. Chem. A 1999, 103, 7753− 7756. (b) Shao, J.; Makri, N. Forward-backward semiclassical dynamics with linear scaling. J. Phys. Chem. A 1999, 103, 9479−9486. (c) Makri, N.; Nakayama, A.; Wright, N. Forward-backward semiclassical simulation of dynamical properties in liquids. J. Theor. Comput. Chem. 2004, 3, 391−417. (21) Wigner, E. J. Calculation of the Rate of Elementary Association Reactions. J. Chem. Phys. 1937, 5, 720. (22) Husimi, K. P. Phys. Math. Soc. Jpn. 1940, 22, 264. (23) (a) Makri, N. Monte Carlo evaluation of forward-backward semiclassical correlation functions with a quantized coherent state density. J. Phys. Chem. B 2002, 106, 8390−8398. (b) Wright, N. J.; Makri, N. Forward-backward semiclassical dynamics for condensed phase time correlation functions. J. Chem. Phys. 2003, 119, 1634− 1642. (24) (a) Shi, Q.; Geva, E. Semiclassical theory of vibrational energy relaxation in the condensed phase. J. Phys. Chem. A 2003, 107, 9059−

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This material is based upon work supported by the National Science Foundation under Award CHE 13-62826. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana−Champaign and its National Center for Supercomputing Applications. T.C.A. is a Blue Waters Graduate Fellow.

■

REFERENCES

(1) (a) Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction path Hamiltonian for polyatomic molecules. J. Chem. Phys. 1980, 72, 99− 112. (b) Caldeira, A. O.; Leggett, A. J. Path integral approach to quantum Brownian motion. Phys. A 1983, 121, 587−616. (c) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, M. Dynamics of the dissipative two-state system. Rev. Mod. Phys. 1987, 59, 1−85. (2) Makri, N. The linear response approximation and its lowest order corrections: an influence functional approach. J. Phys. Chem. B 1999, 103, 2823−2829. (3) (a) Feynman, R. P. Space-time approach to non-relativistic quantum mechanics. Rev. Mod. Phys. 1948, 20, 367−387. (b) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGrawHill: New York, 1965. (4) Feynman, R. P.; Vernon, F. L. The theory of a general quantum system interacting with a linear dissipative system. Ann. Phys. 1963, 24, 118−173. (5) (a) Makri, N.; Makarov, D. E. Tensor multiplication for iterative quantum time evolution of reduced density matrices. I. Theory. J. Chem. Phys. 1995, 102, 4600−4610. (b) Makri, N.; Makarov, D. E. Tensor multiplication for iterative quantum time evolution of reduced density matrices. II. Numerical methodology. J. Chem. Phys. 1995, 102, 4611−4618. (c) Shao, J.; Makri, N. Iterative path integral calculation of quantum correlation functions for dissipative systems. Chem. Phys. 2001, 268, 1−10. (d) Shao, J.; Makri, N. Iterative path integral formulation of equilibrium correlation functions for quantum dissipative systems. J. Chem. Phys. 2002, 116, 507−514. (6) Makri, N. Improved Feynman propagators on a grid and nonadiabatic corrections within the path integral framework. Chem. Phys. Lett. 1992, 193, 435−444. (7) (a) Makri, N. Effective non-oscillatory propagator for Feynman path integration in real time. Chem. Phys. Lett. 1989, 159, 489−498. (b) Makri, N. On smooth Feynman propagators for real time path integrals. J. Phys. Chem. 1993, 97, 2417−2424. (8) Topaler, M.; Makri, N. System-specific discrete variable representations for path integral calculations with quasi-adiabatic propagators. Chem. Phys. Lett. 1993, 210, 448. (9) (a) Makri, N. Numerical path integral techniques for long-time quantum dynamics of dissipative systems. J. Math. Phys. 1995, 36, 2430−2456. (b) Makri, N. Quantum dissipative systems: a numerically exact methodology. J. Phys. Chem. A 1998, 102, 4414−4427. (10) (a) Sim, E.; Makri, N. Tensor propagator with weight-selected paths for quantum dissipative dynamics with long-memory kernels. Chem. Phys. Lett. 1996, 249, 224−230. (b) Sim, E.; Makri, N. Filtered propagator functional for iterative dynamics of quantum dissipative systems. Comput. Phys. Commun. 1997, 99, 335−354. (c) Sim, E. Quantum dynamics for a system coupled to slow baths: on-the-fly filtered propagator method. J. Chem. Phys. 2001, 115, 4450−4456. (d) Lambert, R.; Makri, N. Memory path propagator matrix for long4176


Article

Journal of Chemical Theory and Computation 9069. (b) Frantsuzov, P. A.; Neumaier, A.; Mandelshtam, V. A. Gaussian resolutions for equilibrium density matrices. Chem. Phys. Lett. 2003, 381, 117. (c) Frantsuzov, P. A.; Mandelshtam, V. A. Gaussian resolutions of equilibrium density matrices. J. Chem. Phys. 2004, 121, 9247. (d) Bose, A.; Makri, N. Evaluation of the Wigner distribution via classical adiabatic switching. J. Chem. Phys. 2015, 143, 114114. (25) Makri, N. Exploiting classical decoherence in dissipative quantum dynamics: Memory, phonon emission, and the blip sum. Chem. Phys. Lett. 2014, 593, 93−103. (26) Weiss, U. Quantum Dissipative Systems; World Scientific: Singapore, 1993. (27) Gardiner, C. Stochastic Methods: A Handbook for the Natural and Social Sciences; Springer-Verlag: Berlin, 2009. (28) Lee, S.; Karplus, M. Browinian dynamics simulations: Statistical error of correlation functions. J. Chem. Phys. 1984, 81, 6106−6118. (29) (a) Azzouz, H.; Borgis, D. A quantum molecular-dynamics study of proton transfer reactions along asymmetrical H bonds in solution. J. Chem. Phys. 1993, 98, 7361−7374. (b) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: molecular dynamics with quantum transitions. J. Chem. Phys. 1994, 101, 4657−4667. (30) Andersen, H. C. Rattle: a ’velocity ’ version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys. 1983, 52, 24−34. (31) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J. C.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802.

4177


Direct Computation of Influence Functional Coefficients from

Recommend Documents