Self-Consistent Filtering Scheme for Efficient Calculations of

Article pubs.acs.org/JCTC

Self-Consistent Filtering Scheme for Efficient Calculations of Observables via the Mixed Quantum-Classical Liouville Approach David Dell’Angelo and Gabriel Hanna* Department of Chemistry, University of Alberta, Edmonton, Alberta T6G 2G2, Canada S Supporting Information *

ABSTRACT: Over the past decade, several algorithms have been developed for calculating observables using mixed quantum-classical Liouville dynamics, which differ in how accurately they solve the quantum-classical Liouville equation (QCLE). One of these algorithms, known as sequential short-time propagation (SSTP), is a surface-hopping algorithm that solves the QCLE almost exactly, but obtaining converged values of observables requires very large ensembles of trajectories due to the rapidly growing statistical errors inherent to this algorithm. To reduce the ensemble sizes, two filtering schemes (viz., observable cutting and transition filtering) have been previously developed and effectively applied to both simple and complex models. However, these schemes are either ad hoc in nature or require significant trial and error for them to work as intended. In this study, we present a self-consistent scheme, which, in combination with a soundly motivated probability function used for the Monte Carlo sampling of the nonadiabatic transitions, avoids the ad hoc observable cutting and reduces the amount of trial and error required for the transition filtering to work. This scheme is tested on the spin-boson model, in the weak, intermediate, and strong coupling regimes. Our transition filtered results obtained using a newly proposed probability function, which favors the sampling of nonadiabatic transitions with small energy gaps, show a significant improvement in accuracy and efficiency for all coupling regimes over the results obtained using observable cutting and the original implementation of transition filtering and are comparable to those obtained using the combination of these two techniques. It is therefore expected that this novel scheme will substantially reduce ensemble sizes and simplify the computation of observables in more complex systems.

■

INTRODUCTION In the study of quantum processes, one is often primarily interested in a particular subset of the degrees of freedom (DOF) of a system. In such cases, the dynamics of this subset may be treated quantum mechanically, while the remaining degrees of freedom in its environment (or bath) may be treated classically to a good approximation. Over the past few decades, the field of theoretical chemistry has witnessed the development of a host of methods for simulating the nonadiabatic dynamics of such mixed quantum-classical systems, which essentially differ in the way they deal with quantum decoherence.1−8 One of the most accurate methods for accounting for quantum coherence effects in mixed quantumclassical simulations is Mixed Quantum-Classical Liouville (MQCL) dynamics.6,9 In contrast to the popular Fewest Switches Surface Hopping approach1 and its extensions,7,8,10−12 the MQCL approach is rigorously derived from the underlying quantum mechanical equations of motion and, therefore, properly incorporates quantum coherence effects by allowing evolution of the classical DOF on both single and coherently coupled pairs of potential energy surfaces (PESs). Owing to its rigorous foundation, however, this method is computationally expensive for many-particle systems and, hence, has been primarily restricted to minimal model systems. Nevertheless, © 2016 American Chemical Society

this approach holds great promise for unlocking new insights into a wide range of systems if the computational expenses can be reduced by making sensible approximations. Several algorithms have been developed over the years for performing MQCL dynamics, which differ in how accurately they solve the quantum-classical Liouville equation (QCLE).6,13,14 The Sequential Short-Time Propagation (SSTP)15,16 and Trotter-Based Quantum Classical (TBQC)17 algorithms involve generating an ensemble of trajectories in which the classical DOF hop between PESs dictated by the adiabatic state of the quantum subsystem. Although SSTP and TBQC solve the MQCL equation almost exactly, obtaining converged solutions is very challenging due to the rapidly growing statistical error inherent to these algorithms. On the other hand, the Poisson Bracket Mapping Equation (PBME)18−20 and the Forward−Backward Trajectory Solution (FBTS)21 algorithms involve generating an ensemble of trajectories in which both the quantum and classical DOF conveniently evolve according to classical-like equations of motion. Despite their simplicity and efficiency, these methods solve the QCLE more approximately and can yield inaccurate Received: October 15, 2015 Published: January 5, 2016 477

DOI: 10.1021/acs.jctc.5b00991 J. Chem. Theory Comput. 2016, 12, 477−485

Article

Journal of Chemical Theory and Computation

(super)operator, 3 , on the operator representing that observable6,9

results in situations where a mean-field picture of the dynamics is not valid.22−24 Thus, efforts toward making the SSTP and TBQC algorithms more efficient are worthwhile. Previously, two schemes, known as observable cutting16 and transition filtering,25−27 have been developed and used, both separately and together, to tame the temporal growth of the statistical weights in the calculation of observables via the SSTP algorithm. These weights lead to large fluctuations in the calculated observables at relatively short times and make it very difficult to obtain converged results at longer times. To attenuate the growth of these fluctuations, the observable cutting scheme imposes an upper bound or threshold on the magnitude of the weight. This cutting ensures that the magnitude of the weight never becomes too large, but a series of calculations is needed to determine the optimal choice of the threshold value. In contrast, transition filtering reduces the statistical fluctuations by disallowing those nonadiabatic transitions that would lead to large changes in the bath momenta from occurring. However, these schemes are either ad hoc in nature or require significant trial and error for them to work effectively. Therefore, a self-consistent scheme, which avoids the ad hoc observable cutting and reduces the amount of trial and error required for the transition filtering, is needed. In this study, we focus on developing such a scheme within the context of the SSTP algorithm for calculations of observables. Although TBQC has been shown to yield excellent results for stronger nonadiabatic couplings and longer times using fewer trajectories than SSTP,17 we choose to develop and test the new scheme on the spin-boson model28 within the context of the SSTP algorithm since it is easier to implement than TBQC. Nevertheless, this scheme is sufficiently general, such that it can also be applied to the TBQC algorithm. In developing this scheme, we will take into account the following considerations: (i) The scheme should provide a well-founded prescription for assigning values to any free parameters in the transition probability. (ii) The scheme, when combined with a suitable probability function used for the Monte Carlo (MC) sampling of the nonadiabatic transitions, should not rely on observable cutting to improve its performance. (iii) The scheme should rely solely on the exact form of the momentum jump approximation16,29 used in the evaluation of the nonadiabatic jump operator. (iv) The scheme should admit long-time simulations of observables for a wide range of subsystem-bath coupling strengths by reducing the number of trajectories required for convergence. The structure of the paper is as follows. Section II describes the MQCL approach for computing a general observable and highlights the elements of the SSTP algorithm that will be taken into consideration in the remainder of the paper. Section III introduces three new transition probabilities, an amended version of the original transition filtering implementation, and a self-consistent scheme for obtaining converged results. The model on which we test the new scheme and the simulation details are described in Section IV. The results of the testing are given in Section V. Finally, concluding remarks are made in Section VI.

∂ ̂ ̂ (X , t ) AW (X , t ) = i 3AW ∂t i ̂ ̂ (X , t )] ≡ [HW (X , t ), AW ℏ 1 ̂ (X , t ), AW ̂ (X , t )} − {AW ̂ (X , t ), HW ̂ (X , t )}) − ({HW 2

(1)

where the subscript W indicates that a partial Wigner transformation has been performed over the classical DOF, and X collectively refers to the positions R and momenta P of the classical DOF. In the above equation, Ĥ W is the system Hamiltonian and may be written as ̂ (t ) = Te + Ve + Tq̂ + Vq̂ + Vĉ HW

(2)

where Te and Ve are kinetic and potential energies of the classical environment, respectively; T̂ q and V̂ q are the kinetic and potential energy operators of the quantum subsystem, respectively; and V̂ c is the subsystem-environment coupling potential energy operator. The square brackets around the first term on the right-hand side of eq 1 denote the commutator, and the curly brackets around the last two terms denote the Poisson bracket. On the basis of adiabatic states, |α;R⟩, given by the solutions of the time-independent Schrödinger equation for a system with Hamiltonian Ĥ W (i.e., Ĥ W|α;R⟩ = Eα(R) |α;R⟩), the QCLE is ∂ αα ′ AW (X , t ) = ∂t

∑ [(iωαα′ + iLαα′)δαβδα′ β ′ − 1αα′,ββ ′]AWββ ′(X , t ) ββ ′

(3)

The first two terms in the square brackets on the right-hand side of this equation describe the evolution of the classical environment by the classical Liouville operator L accompanied by quantum mechanical phase oscillations with frequency ωαα′ = (Eα − Eα′)/ℏ. The classical Liouville operator propagates the classical subsystem via the Hellmann−Feynman forces, FαW = ∂V̂ (R )

−⟨α ; R | ∂cR |α ; R ⟩, on either the potential energy surface corresponding to one of the adiabatic states α with energy Eα or on the mean potential energy surface between states α and α′ with energy (Eα + E′α)/2, and is given by

iLαα ′ =

P ∂ 1 ∂ · + (FWα + FWα ′) · M ∂R 2 ∂P

(4)

where M is the vector of masses of the particles belonging to the environment. The third term, which is responsible for nonadiabatic transitions, is given by ⎛ ∂ ⎞ 1 P ·dαβ ⎜1 + Sαβ · ⎟δα ′ β ′ ⎝ ∂P ⎠ M 2 P * ⎛⎜ 1 * ∂ ⎟⎞ − ·dα ′ β ′ 1 + Sα ′ β ′· δαβ ⎝ ∂P ⎠ M 2

1αα ′ , ββ ′(t ) = −

(5)

∂

where dαβ = ⟨α ; R | ∂R |β ; R ⟩ is the nonadiabatic coupling

■

P

m a t r i x e l e m e n t , Sαβ = FWα δαβ − FWαβ( M ·dαβ)−1 =

MQCL APPROACH FOR COMPUTING OBSERVABLES The QCLE for the time evolution of an observable Â W(X, t) is given by operating with the quantum-classical Liouville

αβ (Eα − Eβ )dαβ( M ·dαβ)−1, and FW denotes the off-diagonal matrix elements of the force. The classical momentum derivatives are responsible for changes in the momenta of the

P

478


Article

Journal of Chemical Theory and Computation classical subsystem caused by energy exchange between the quantum and classical subsystems during nonadiabatic transitions of the quantum subsystem. From the derivation of the SSTP algorithm, the solution to 15,16 ′ eq 3 for Aαα W (X, t) is given by AWαα ′(X , t ) =

The multidimensional sums over the state indices and integrals over phase space variables are evaluated via MC sampling with initial conditions sampled according to ραW0′α0(X). The details of the SSTP algorithm for evaluating eq 10 are found in refs 15 and 16. However, we would just like to draw the reader’s attention to the transition probability function, Π0, originally defined in ref 15 and used for the MC sampling of the nonadiabatic transitions along an MQCL trajectory:

⎡ N ⎤ ⎢∏ (ei 3̂ Δt j) ⎥ αNαN′ αj − 1αj′− 1, αjαj′ AW (X ) ⎢ ⎥⎦ (α1α1′) ...(αN αN′ ) ⎣ j = 1

∑

(6)

Π0 =

where the time interval t is divided into N segments such that the jth segment has length Δtj = tj − tj−1. Propagation over a given time segment is performed by making the approximation (e

i 3̂ Δt j

(7)

where >αj−1αj′−1(t j − 1 , t j) = eiωαj−1αj′−1Δt j is the phase factor associated with that time segment. In simulating the quantum-classical dynamics of systems using this algorithm, the momentum-jump approximation is used to compute the action of 1αα ′ , ββ ′ on phase space functions. It is given by16,29

■

IMPROVING SAMPLING EFFICIENCY New Transition Probabilities. In this section, we propose new transition probabilities with the same general form as Π0 in eq 11

⎛ ⎞ ∂ ⎜1 + ΔE ⎟f (P) ≈ eΔEαβ ∂/∂(P̅·dαβ̂ )2f (P) αβ ⎜ ̂ )2 ⎟⎠ ∂(P ̅ ·dαβ ⎝

Πi =

= f (P + ΔP)

|Ai (X )|Δt 1 + |Ai (X )|Δt

̂ = d ̅ /|d ̅ |, the above, ΔEαβ = Eα − Eβ, P ̅ = P / M , and dαβ αβ αβ T M −1/2 , where T ̅ = dαβ / M . (N.B.: dαβ / M = dαβ where dαβ stands for the transpose, and M−1/2 is a diagonal matrix of the inverse square root masses.) The corresponding momentum shift is given by

̂ [sgn(P ̅ ·d ̂ ) (P ̅ ·d ̂ )2 + ΔE − (P ̅ ·d ̂ )] M dαβ αβ αβ αβ αβ

A1(X ) =

⎛ π⎞ P ·dαβ cos⎜ωαβ ⎟ ⎝ M 2⎠

A 2 (X ) =

ΔEαβ P ·dαβ M ΔEαβ

A3(X ) =

⎛ π ⎞ ΔEαβ P ·dαβ cos⎜ωαβ ⎟ ⎝ M 2 ⎠ ΔEαβ

(9)

(

αα ′

⎡

′) (α0α0′) ...(αN αN

N

∫ dX ⎢⎢∏ (ei3̂ Δt )α j

⎣ j=1

j − 1αj′− 1, αjαj′

(15)

A few comments about these expressions are in order. All of P them contain A 0(X ) = M ·dαβ as their base, with one or two extra factors. To arrive at A1(X), A0(X) is multiplied by π cos ωαβ 2 , which corresponds to the real part of eiωαβΔt appearing in eq 7. We chose Δt = π/2 since this corresponds to a point in the cosine function with the maximum slope, resulting in a function that is maximally sensitive to changes in ωαβ over the course of a trajectory. To arrive at A2(X), A0(X) is min multiplied by ΔEmin αβ /ΔEαβ, where ΔEαβ is taken to be the smallest energy gap between the α and β adiabatic surfaces and serves to establish a dimensionless energy scale. In practice, ΔEmin αβ may be determined from histograms of ΔEαβ, generated

∑ ∫ dXAWαα′(X , t )ρWα′α (X ) ∑

(14) min

αβ

the quantum transition to occur), then the argument of the square root is negative leading to imaginary momentum changes. If this occurs, the nonadiabatic transition is forbidden, and the trajectory is continued adiabatically. It should be noted that the total energy of the system is conserved along an MQCL trajectory when the momentum-jump approximation is invoked. Within the MQCL approach, the expectation value of an observable Â W(X, t) is given by

=

(13)

min

According to this equation, if ΔEαβ < 0 (i.e., an upward ̂ )2 < |ΔE | (i.e., there is transition from α → β) and (P ̅ ·dαβ αβ insufficient kinetic energy from the bath momenta along d ̂ for

⟨A(t )⟩ =

(12)

where A(X) should be chosen in such a way to efficiently sample the multidimensional sums and phase space integrals in eq 10 and thereby minimize the temporal growth of the statistical weights. Thus, one may pick forms of A(X) that are proportional to the terms which multiply the approximate ̂̅ 2 momentum jump operator eΔEαβ∂/∂(P̅·dαβ) (see eq 7 for these terms). To further minimize the weight growth, one may also introduce terms that favor the sampling of nonadiabatic transitions with smaller ΔEαβ’s. With this motivation in mind, we propose the following three new A(X)’s:

(8)

∂ 1 ̂ )2 was used. In where the fact that 2 Sαβ · ∂P = ΔEαβ ∂/∂(P ̅ ·dαβ

ΔP =

(11)

It should be noted that this definition is arbitrary but is chosen in such a way that minimizes the temporal growth of the statistical weights, ∏Nj=1(Πj)−1, associated with the sampling. This particular choice was motivated by the form of 1αα ′ , ββ ′ and the second line of eq 7. In the next section, we will propose other definitions of the transition probability, which are aimed at further minimizing the growth of the statistical weights.

)αj−1αj′−1, αjαj′ ≈ >αj−1αj′−1(t j − 1 , t j)eiLαj−1αj′−1Δt j × (δαj−1αjδ αj′−1αj′ + Δt 1αj−1αj′−1, αjαj′)

1

P ·d Δt M αβ P + M ·dαβ Δt

⎤ αN αN′ α0′α0 ⎦AW (X )ρW (X )

(10) 479

)


Article


corresponding to the maximum in the ΔEαβ histogram and will depend on the subsystem-bath coupling strength (since the width of the ΔEαβ histogram varies with the coupling strength (see Figure 1)).

by evolving the system adiabatically on its ground-state surface for a sufficiently long period of time. By placing ΔEαβ in the denominator of A2(X), one favors the sampling of nonadiabatic transitions with smaller energy gaps. This also has the effect of favoring the sampling of transitions which lead to larger Re(eiωαβΔt) terms since, in the small ωαβΔt limit, the cosine is maximized by smaller values of ωαβ. Finally, A3(X) simply combines the features introduced into A1(X) and A2(X). Transition Filtering. The original implementation of transition filtering was defined in terms of the following transition probability25−27 P ·d w(κϵ , ϵαβ ) M αβ P Δt M ·dαβ w(κϵ , ϵαβ )

Δt Π tf =

1+

(16)

where ⎧1, if |ϵαβ | ≤ κϵ w(κϵ , ϵαβ ) = ⎨ ⎪ ⎩ 0, otherwise

Figure 1. Histograms of the energy gaps between the ground and firstexcited state PESs for various subsystem-bath coupling regimes. From left to right: (a) weak, (b) intermediate, and (c) strong coupling.

⎪

(17)

Obtaining Converged Results. In this section, we present a self-consistent scheme for obtaining converged results using the SSTP algorithm with the transition probabilities and filtering presented in the previous subsections. This scheme is general and not meant to be coupled with observable cutting to achieve improved results. The steps of this scheme are as follows: 1. Evolve the system adiabatically on its ground-state surface for a sufficiently long period of time, so as to explore all of the configurations of relevance to the study. Calculate the various ΔEαβ’s at each time step of the evolution and generate histograms of the ΔEαβ’s. When using eq 14 or eq 15, ΔEmin αβ corresponds to the smallest value of ΔEαβ with a nonzero probability density in the histogram. 2. Choose initial trial values for the maximum number of nonadiabatic transitions allowed in a given trajectory, nmax, and the transition filtering thresholds, {καβ}. A reasonable starting point for nmax is two, while reasonable starting points for {καβ} are the values of ΔEαβ corresponding to the maxima in the ΔEαβ histograms. 3. Using the SSTP algorithm, calculate the time-dependent expectation value in eq 10 out to some desired time, using a sufficient number of trajectories for convergence. 4. As long as the computational time is reasonable, repeat step 3 with increased values of nmax and {καβ} until the result agrees with that of the previous calculation to within the statistical noise.

and ϵαβ =

⎞ ⎛ P2 P′2 + Eα(R ) − ⎜ + Eβ (R )⎟ 2M ⎠ ⎝ 2M

(18)

In the above equation, P′ = P + ΔP̃ , where ΔP̃ is an approximation to the momentum shift in eq 9, obtained by expanding the square root on the right-hand side of the ̂ )2 : equation to first order in ΔEαβ /(P ̅ ·dαβ Δ̃P =

ΔEαβ ̂ 2P ̅ ·dαβ

̂ M dαβ (19)

This particular form of w(κϵ,ϵαβ) filters out those nonadiabatic transitions with sufficiently large ΔEαβ’s and that would lead to excessively large changes in the bath momenta. More specifically, the transition probability in eq 16 limits the magnitude of the energy change associated with the virtual approximate momentum shift30 in eq 19 and the adiabatic energy gap, to be no greater than some threshold value κϵ. Whenever a nonadiabatic transition causes |ϵαβ| to be greater than or equal to κϵ, the transition probability becomes zero, and no transition can occur. As in observable cutting, κϵ is chosen via trial and error to reduce the statistical noise in the results. A more direct form of transition filtering, which avoids the use of the virtual approximate momentum shift and reduces the arbitrariness in choosing appropriate κ values, can be defined in terms of the following transition probability Π′i =

■

MODEL AND SIMULATION DETAILS The scheme proposed in this work will be tested on the spinboson model,28 for which numerically exact results exist.31−34 This model consists of a two-level system, with spin states {|↑⟩,|↓⟩} bilinearly coupled to a bath of M harmonic oscillators. The partially Wigner-transformed Hamiltonian for this model is given by

Δt |Ai (X )|w(καβ , ΔEαβ ) 1 + Δt |Ai (X )|w(καβ , ΔEαβ )

(20)

⎧ if (|ΔEαβ| ≤ καβ) ⎪1, w(καβ , ΔEαβ ) = ⎨ ⎪ ⎩ 0, otherwise

(21)

where

̂ = −Ωσx̂ + 1 HW 2

This form of the transition filtering strives to restrict nonadiabatic transitions to occur only in regions of an adiabatic PES where ΔEαβ is sufficiently small. In this case, a sensible initial choice of καβ can be the value of ΔEαβ corresponding to the maximum in the ΔEαβ histogram. In general, the optimal value of καβ will be greater than or equal to the value of ΔEαβ

M

∑ (Pj2 + ωj2R j2 − cjR jσẑ ) j=1

(22)

where σ̂x and σ̂z are the Pauli spin matrices, Rj and Pj are the mass-weighted coordinate and momentum, respectively, of the jth bath mode with frequency ωj, and cj is the subsystem-bath coupling coefficient. The bilinear coupling between the 480


Article


spin-boson model, which will help us to identify ΔEmin 12 ≡ ΔEmin and to choose initial trial values for κ12. These histograms are shown in Figure 1 for the weak, intermediate, and strong subsystem-bath coupling regimes. From the figures, we can see that ΔEmin = 0.66, 0.80, and 2.40 for the weak, intermediate, and strong coupling cases, respectively. We can also take these values as the initial trial values for κ12. In practice, the optimal value of κ12 will be κ12 = ΔEmin + δ, where δ depends on the subsystem-bath coupling strength. Figure S1 shows the effect of using ΔEmin values other than the actual ΔEmin value on the expectation value of σ̂z, for the system in the strong coupling regime. The results were generated using the Π2 transition probability and nmax = 6. As can be seen, very good agreement with the exact result37 is achieved for ΔEmin = 2.40 (i.e., the exact minimum energy gap from Figure 1), while the other values lead to deviations at some stage. Weak Coupling. To investigate the effect of the new transition probabilities as defined by eq 12 with i = 0...3, we plot the temporal growth of the (ensemble-averaged) statistical weights for several values of nmax up to t = 10. As seen in Figure 2, the new transition probabilities Π2 and Π3 yield slower

subsystem and bath is characterized by an Ohmic spectral density, such that cj = ξ(ω0 /ωc) ωj , where ξ is the Kondo parameter, ωc is a characteristic frequency of the spectral density, ωj = −ωc ln(1 − jω0/ωc), and ω0 = ωc[1−exp(−ωmax/ ωc)]/M (ωmax is the maximum frequency of the spectral density). It should be noted that since the subsystem-bath coupling is bilinear, the QCLE yields the exact quantum dynamics for this model. The analytical forms of the adiabatic states and nonadiabatic coupling matrix elements are known for this model and can be found in ref 35. For convenience, we use dimensionless variables and parameters with time scaled by ωc in all of our simulations. The initial density operator is assumed to be uncorrelated, with the subsystem in its |↑⟩ state and the bath in thermal equilibrium: ρŴ (X , 0) = ρs ̂ (0)ρe (X )

(23)

Thus, the subsystem density matrix is ⎛1 0⎞ ⎟ ρs (0) → ⎜ ⎝0 0⎠

(24)

and the Wigner distribution of the bath is given by M

ρe (X ) =

∏ j=1

36

⎡ ⎧ P2 ⎫⎤ ⎪ β⎪ 1 j ⎬⎥ ⎨ + ωj2R j2⎪ exp⎢ − ⎪ ⎢⎣ u″j ⎩ 2 2πu″j 2 ⎭⎥⎦ βωj

(25)

β ℏ ωj

where uj″ = uj coth uj, uj = 2 , and β is the reciprocal temperature. In this study, we consider three model parameter sets, corresponding to a wide range of subsystem-bath couplings (labeled as weak, intermediate, and strong) at different temperatures. The values of these parameters, given in Table 1, are the same as those used in ref 27. In Section V, we present

Figure 2. Effect of the various transition probabilities on the temporal growth of the ensemble-averaged MC weights for different numbers of allowed nonadiabatic jumps. These results were generated with the weak coupling parameter values and ΔEmin = 0.66. From left to right: (a) nmax = 2, (b) nmax = 4, and (c) nmax = 6 nonadiabatic jumps.

Table 1. Spin-Boson Model Parameter Sets, Corresponding to Weak, Intermediate, and Strong Subsystem-Bath Coupling Conditions at Different Temperatures weak intermediate strong

β

ξ

Ω

0.3 12.5 0.25

0.007 0.09 2.0

0.33 0.4 1.2

weight growths for all values of nmax than the original transition probability Π0 and Π1. Moreover, the distinction between the accumulated weights generated using Π2 and Π3 and those generated using Π0 and Π1 becomes more pronounced in the nmax = 4 and 6 cases. Next, we investigate the effect of transition filtering, as defined by eq 20 with i = 0...3, on the temporal growth of the (ensemble-averaged) statistical weights for several values of nmax. As seen in Figure 3, transition filtering suppresses the weight growth by orders of magnitude (cf. Figure 2). The new transition probabilities Π1′ , Π2′ , and Π3′ yield slower weight growth in the nmax = 2 case than Π0′ , but the situation changes in the larger nmax cases. In the nmax = 4 case, all four transition probabilities perform similarly, with Π1′ generating the slowest weight growth. However, in the nmax = 6 case, Π2′ clearly yields the slowest weight growth. Based on the results of Figures 2 and 3, we conclude that the Π2 transition probability, overall, yields the most stable weight growth. Therefore, in what follows, we focus on the comparisons between the Π0/Π0′ and Π2/Π2′ results. For completeness, we have included the results obtained using the other transition probabilities in the Supporting Information. The results for the expectation value of σ̂z obtained without transition filtering are shown in Figure 4(a). They were

results for the time-dependent population difference between the two spin states (i.e., the expectation value of σ̂z(t)), based on evaluating eq 10 via the SSTP algorithm. The integration of the QCLE, as outlined in Section II, is carried out with a dimensionless time step of Δt = 0.1. As mentioned in Section III, we set a bound nmax on the maximum number of nonadiabatic transitions allowed in any given trajectory. We found that 2 < nmax < 6 was ideal for obtaining converged results. For the simulation times considered herein, the number of trajectories required to obtain converged expectation values with transition filtering (i.e., using the transition probabilities defined by eq 20) is 1 × 106, while that used in the calculations without transition filtering (i.e., using the transition probabilities defined by eq 12) is 1 × 107.

■

RESULTS In keeping with the procedure outlined in Section III C, the first step involves generating ΔE12 histograms for the two-state 481


Article


Figure 5. Time-dependent population differences, ⟨σz (t)⟩, generated using the Π0 and Π2 transition probabilities and nmax = 6, for the system in the (a) weak (ΔEmin = 0.66), (b) intermediate (ΔEmin = 0.80), and (c) strong (ΔEmin = 2.40) coupling regimes. For comparison, the exact results from refs 33, 34, and 37 are shown, respectively.

Figure 3. Effect of the various transition probabilities on the temporal growth of the ensemble-averaged MC weights with different numbers of allowed nonadiabatic jumps and transition filtering. These results were generated with the weak coupling parameter values, ΔEmin = 0.66 and κ12 = 0.69. From left to right: (a) nmax = 2, (b) nmax = 4, and (c) nmax = 6 nonadiabatic jumps.

observable cutting. As can be seen, all of the results are quite comparable, with our result being most similar to the transition filtered result of ref 27 with observable cutting. Intermediate Coupling. We now move on to the results for the intermediate coupling regime. In general, as one increases the coupling strength, the number of lowerprobability nonadiabatic transitions increases (since the system will explore more solvent configurations with larger |ΔE12| values), leading to the onset of statistical noise at shorter times. In such cases, καβ plays an even more important role in obtaining the optimal result, when computational time is an issue. To investigate the effect of the new transition probabilities as defined by eq 12 with i = 0...3, we plot the temporal growth of the (ensemble-averaged) statistical weights for nmax = 2, 4, and 6. As seen in Figure S2, the new transition probabilities Π1, Π2, and Π3 yield slower weight growths for all values of nmax than the original transition probability Π0, especially in the nmax = 4 and 6 cases. Moreover, the difference between the accumulated weights generated using the new transition probabilities and those generated using Π0 becomes more pronounced in the nmax = 4 and 6 cases. However, of the three new transition probabilities, Π2 and Π3 yield the most stable weight growth in all cases. Next, we investigate the effect of transition filtering, as defined by eq 20 with i = 0...3, on the temporal growth of the (ensemble-averaged) statistical weights for nmax = 2, 4, and 6. As seen in Figure S3, transition filtering suppresses the weight growth by orders of magnitude (cf. Figure S2). The new transition probabilities Π1′ , Π2′ , and Π3′ yield slower weight growth in all cases than Π′0, with the differences between them decreasing with increasing nmax. In all cases, Π′2 generates the slowest weight growth. Based on the results of Figures S2 and S3, we conclude that the Π2 transition probability, overall, yields the most stable weight growth. Therefore, in what follows, we focus on the comparisons between the Π0/Π′0 and Π2/Π2′ results. For completeness, we have included the results obtained using the other transition probabilities in the Supporting Information. The results for the expectation value of σ̂z obtained without transition filtering are shown in Figure 4(b). They were calculated using the Π2 transition probability with ΔEmin = 0.80 and different numbers of allowed nonadiabatic transitions. For comparison, the numerically exact results are shown.34 Figure 5(b) presents the corresponding results generated using the Π0 and Π2 transition probabilities. As seen, the Π2 result is in

Figure 4. Time-dependent population difference, ⟨σz (t)⟩, generated using the Π2 transition probability and different numbers of allowed nonadiabatic jumps, for the system in the (a) weak (ΔEmin = 0.66), (b) intermediate (ΔEmin = 0.80), and (c) strong (ΔEmin = 2.40) coupling regimes. For comparison, the exact results from refs 33, 34, and 37 are shown, respectively.

calculated using the Π2 transition probability with ΔEmin = 0.66 and different numbers of allowed nonadiabatic transitions. For comparison, the numerically exact results from ref 33 are shown. As expected, increasing nmax from two to six yields better agreement with the exact results at longer times. Figure 5(a) presents the results obtained using the Π0 and Π2 transition probabilities. As seen, the results are almost identical up to t = 4, after which deviations begin to occur. We now move on to the results obtained with transition filtering. Figures 6(a) and 6(b) show the results for the expectation value of σ̂z. In Figure 6(a), we compare the results generated using Π′2, ΔEmin = 0.66, nmax = 2, and different transition filter values, starting from κ12 = 0.69 to κ12 = 0.80, in order to settle upon an appropriate value of κ12. Overall, we see that that small variations in καβ can lead to slower convergence at longer times. As can be seen, κ12 = 0.69 performs the best at longer times. We next check for convergence with respect to nmax. Figure 6(b) shows the results of increasing nmax from two to six, generated using Π′2 and κ12 = 0.69. Here, we see that the result changes very little as one increases nmax, indicating that the result is more or less converged with nmax = 2 (for the time period shown). Figure 7(a) presents the results obtained using the Π′0 and Π′2 transition probabilities with κ12 = 0.69. As seen, the results are very similar. In Figure 8(a), we compare our transition filtered result generated using Π2′ , κ12 = 0.69, and nmax = 2 to the transition filtered results of ref 27 with and without 482


Article


Figure 6. Effect of transition filtering with the Π2′ transition probability on the time-dependent population differences, ⟨σz (t)⟩, for the system in the weak [(a) and (b)], intermediate [(c) and (d)], and strong [(e) and (f)] coupling regimes. Results were generated using (a) ΔEmin = 0.66 and nmax = 2 for various values of the transition filter κ12; (b) ΔEmin = 0.66 and κ12 = 0.69 for various values of nmax; (c) ΔEmin = 0.80 and nmax = 2 for various values of the transition filter κ12; (d) ΔEmin = 0.80 and κ12 = 0.82 for various values of nmax; (e) ΔEmin = 2.40 and nmax = 2 for various values of the transition filter κ12; (f) ΔEmin = 2.40 and κ12 = 3.20 for various values of nmax.

transition filtered result generated with Π2′ , κ12 = 0.82, and nmax = 2 to the transition filtered results of ref 27 generated with and without observable cutting (and nmax = 2) and also to the TBQC result of ref 17 generated with observable cutting. As times goes on, the magnitude and phase of the oscillation in our transition filtered result and in the ref 27 result without observable cutting become increasingly different from those of the exact result. However, our result exhibits a pronounced decrease in the magnitude of the oscillations (reminiscent of the exact result), while the transition filtered result of ref 27 without observable cutting exhibits a very slight decrease in the magnitude. Combining the transition filtering with observable cutting improves the agreement with the exact result, but we believe that this improvement is fortuitous, since for such a long simulation time, a larger value of nmax is actually required for convergence. The TBQC result generated with observable cutting agrees best with the exact result, which suggests that implementing a form of transition filtering within the TBQC algorithm would lead to significant improvements. It should be noted that the ability of TBQC to more accurately simulate the dynamics for longer times than SSTP stems from the fact that TBQC is based on a Trotter factorization of the quantumclassical propagator. Strong Coupling. The strong coupling regime represents the most challenging case since the largest number of lower probability nonadiabatic transitions is expected within a given time compared to the other coupling cases. To investigate the effect of the new transition probabilities as defined by eq 12 with i = 0...3, we plot the temporal growth of the (ensembleaveraged) statistical weights for nmax = 2, 4, and 6. As seen in Figure S4, Π2 yields the slowest weight growth when nmax = 2, performing most similarly to Π3. However, the distinction between the rates of growth generated by the various transition probabilities becomes more pronounced with increasing nmax. When nmax = 4 and 6, the Π1, Π2, and Π3 transition probabilities yield substantially slower weight growths than Π0, with Π2 and Π3 yielding the slowest weight growth in the nmax = 6 case. Next, we investigate the effect of transition filtering, as defined by eq 20 with i = 0...3, on the temporal growth of the (ensemble-averaged) statistical weights for nmax = 2, 4, and 6.

Figure 7. Time-dependent population differences, ⟨σz (t)⟩, obtained using the Π0′ and Π2′ transition probabilities and nmax = 2, for the system in the (a) weak (ΔEmin = 0.66 and κ12 = 0.69), intermediate (ΔEmin = 0.80 and κ12 = 0.82), and (c) strong (ΔEmin = 2.40 and κ12 = 3.20) coupling regimes.

decent agreement with the exact result, while the Π0 result deviates quite substantially from the exact result starting at t ≈ 3.5. As in the weak coupling case, only the Π2 result is reliable at later times. We now move on to the results obtained with transition filtering. Figure 6(c) shows the results for the expectation value of σ̂z with ΔEmin = 0.80, Π2′ , nmax = 2, and different transition filter values, starting from κ12 = 0.82 to κ12 = 0.85, in order to settle upon an appropriate value of κ12. Based on these results, we see that κ12= 0.82 generates the best result in the intermediate coupling case. We next check for convergence with respect to nmax. Figure 6(d) shows the results of increasing nmax from two to six, with κ12 = 0.82. We see that the result changes slightly as one increases nmax, indicating that the result is more or less converged with nmax = 2 up to t ≈ 9.5. As expected, after t ≈ 9.5, the agreement with the exact result improves slightly with nmax = 6. Figure 7(b) presents the results obtained using the Π0′ and Π2′ transition probabilities with κ12= 0.82. As seen, the results are very similar only up to t ≈ 3. Beyond this time, the results differ from one other. Comparing to the exact results, we see that Π′2 is the most effective transition probability and that Π0′ leads to a substantial deviation at longer times. In Figure 8(b), we compare our 483


Article


We now move on to the results obtained with transition filtering. Figure 6(e) shows the strong coupling results for the expectation value of σ̂z generated using Π2′ , nmax = 2, and ΔEmin = 2.40 for different κ12 values ranging from 3.0 to 4.0. As can be seen, a value of κ12 = 3.20 yields the best agreement with the exact result. Figure 6(f) shows the results of increasing nmax from two to six, with κ12 = 3.20. We see that the best agreement with the exact result occurs when nmax = 2. When nmax increases, the results deviate more from the exact result (starting at t ≈ 0.6), in contrast to the unfiltered results where increasing nmax improves the agreement. Figure 7(c) presents the results obtained using Π0′ and Π2′ with κ12= 3.20, ΔEmin = 2.40, and nmax = 2. As can be seen, Π′2 performs the best, while Π′0 yields results which deviate from the exact result at very early times. In Figure 8(c), we compare our transition filtered result generated using Π2′ , κ12 = 3.20, ΔEmin = 2.40, and nmax = 2 to the transition filtered results of ref 27 generated with and without observable cutting (and nmax = 2) and also to the TBQC result of ref 17 generated with observable cutting. We first point out that our “benchmark” result, generated using the Π2 transition probability and nmax = 6, gives the best agreement with the exact result and is comparable to the result of ref 27 obtained using a combination of their implementation of transition filtering and observable cutting (see panel 4(c)). Our transition filtered result, generated using Π′2 and nmax = 2, and the result from ref 17, generated using observable cutting, agree very well with the exact result up to t ≈ 1.2 but deviate slightly for the remainder of the time. This should be contrasted with the transition filtered result of ref 27, which begins to deviate from the exact result much earlier (at t ≈ 0.5). Combining transition filtering with observable cutting does improve the agreement with the exact result at later times, as seen earlier, but this improvement is fortuitous due to the ad hoc nature of observable cutting.

■

Figure 8. Comparison of our results generated using Π2′ and nmax = 2 to the TBQC results of ref 17 obtained with observable cutting (o.c.), and to the results of ref 27 obtained using their implementation of transition filtering (t.f.) and a combination of t.f. and o.c., for the system in the (a) weak (ΔEmin = 0.66 and κ12 = 0.69), (b) intermediate (ΔEmin = 0.80 and κ12 = 0.82), and (c) strong (ΔEmin = 2.40 and κ12 = 3.20) coupling regimes.

CONCLUSIONS In this article, we considered three new probability functions for stochastically sampling nonadiabatic transitions in the numerical solution of the QCLE. These transition probabilities were chosen in an attempt to minimize the growth of the statistical weights and to be used within a self-consistent scheme that does not rely on ad hoc observable cutting to improve its performance. We tested these transition probabilities in the calculation of the time-dependent population difference, ⟨σz (t)⟩, in the spin-boson model in the weak, intermediate, and strong coupling regimes (for which numerically exact results are available). In particular, our results obtained using the Π′2 transition probability show a significant improvement in accuracy for all coupling regimes over the results obtained using the original transition probability with either observable cutting or the implementation of transition filtering from ref 27. However, the fact that the intermediate coupling TBQC result of ref 17 obtained with observable cutting agrees much better with the exact result suggests that further improvements over our SSTP results may be achieved by implementing a form of transition filtering within the TBQC algorithm. In terms of the improvement in efficiency, our transition filtered results were generated with an order of magnitude fewer trajectories than in the case of the results obtained without transition filtering. Most importantly, our results are comparable to those obtained using a combination of observable cutting and transition filtering from ref 27, which demonstrates that one can likely avoid the use of observable

As seen in Figure S5, transition filtering suppresses the weight growth by orders of magnitude (cf. Figure S4). For this coupling strength, all of the transition probabilities yield similar weight growths, with Π3′ yielding the slowest weight growth when nmax = 6. Based on the above results, the Π2 transition probability performs the best overall. Therefore, in what follows, we focus on the comparisons between the Π0/Π′0 and Π2/Π2′ results. For completeness, we have included the results obtained using the other transition probabilities in the Supporting Information. Figure 4(c) shows the results for the expectation value of σ̂z obtained without transition filtering. They were calculated using the Π2 transition probability with ΔEmin = 2.40 and different numbers of allowed nonadiabatic transitions. For comparison, the numerically exact results from ref 37 are shown. As in the weak and intermediate coupling cases, the agreement with the exact results improves upon increasing nmax from two to six. Figure 5(c) presents the corresponding results generated using the Π0 and Π2 transition probabilities. As seen, the Π2 results are in good agreement with the exact results, while the Π0 results deviate from the exact results starting at t ≈ 1.5. 484


Article

Journal of Chemical Theory and Computation cutting. Based on our results, we believe that the Π′2 transition probability, which was shown to most effectively sample nonadiabatic transitions with transition filtering, should be generally used for more efficient SSTP calculations of observables. Our findings hold promise for more efficient QCLE-based calculations of time-dependent observables in systems containing more than one quantum DOF in a complex, many-body, classical environment. In such cases, the solution of the QCLE may be complicated by the involvement of more than two adiabatic states in the dynamics. It is expected that the scheme outlined herein, in combination with the Π2′ transition probability, will substantially reduce the number of trajectories required to achieve accurate and stable results. Work along this direction is currently underway in our group.

■

(17) MacKernan, D.; Ciccotti, G.; Kapral, R. J. Phys. Chem. B 2008, 112, 424. (18) Kim, H.; Nassimi, A.; Kapral, R. J. Chem. Phys. 2008, 129, 084102. (19) Nassimi, A.; Bonella, S.; Kapral, R. J. Chem. Phys. 2010, 133, 134115. (20) Kelly, A.; van Zon, R.; Schofield, J. M.; Kapral, R. J. Chem. Phys. 2012, 136, 084101. (21) Hsieh, C.-Y.; Kapral, R. J. Chem. Phys. 2012, 137, 22A507. (22) Rekik, N.; Hsieh, C.-Y.; Freedman, H.; Hanna, G. J. Chem. Phys. 2013, 138, 144106. (23) Martinez, F.; Rekik, N.; Hanna, G. Chem. Phys. Lett. 2013, 573, 77−83. (24) Martinez, F.; Hanna, G. Mol. Simul. 2015, 41, 107−122. (25) Sergi, A.; Petruccione, F. Phys. Rev. E 2010, 81, 032101. (26) Uken, D. A.; Sergi, A.; Petruccione, F. Phys. Scr. 2011, T143, 014024. (27) Uken, D. A.; Sergi, A.; Petruccione, F. Phys. Rev. E 2013, 88, 033301. (28) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M.; Garg, A.; Zwerger, W. Rev. Mod. Phys. 1987, 59, 1. (29) Sergi, A.; MacKernan, D.; Ciccotti, G.; Kapral, R. Theor. Chem. Acc. 2003, 110, 49−58. (30) The approximate momentum shift is considered virtual, since it is not actually used during the course of the dynamics. It is only used for the purposes of the transition filtering. The exact momentum shift in eq 9 is used during the dynamics. (31) Makri, N.; Thompson, K. Chem. Phys. Lett. 1998, 291, 101. (32) Makri, N. J. Phys. Chem. B 1999, 103, 2823. (33) Thompson, K.; Makri, N. J. Chem. Phys. 1999, 110, 1343. (34) Makarov, D. E.; Makri, N. Chem. Phys. Lett. 1994, 221, 482. (35) MacKernan, D.; Ciccotti, G.; Kapral, R. J. Chem. Phys. 2002, 116, 2346−2353. (36) Kim, H.; Kapral, R. Chem. Phys. Lett. 2006, 423, 76. (37) Mak, C.; Chandler, D. Phys. Rev. A: At., Mol., Opt. Phys. 1991, 44, 2352.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b00991. Figure showing the effect of ΔEmin on the timedependent population difference. Figures showing the effect of the various transition probabilities on the temporal growth of the Monte Carlo weights and on the time-dependent population difference (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by grants from the University of Alberta and the Natural Sciences and Engineering Research Council of Canada (NSERC).

■

REFERENCES

(1) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (2) Webster, F.; Wang, E. T.; Rossky, P. J.; Friesner, R. A. J. Chem. Phys. 1994, 100, 4835−4847. (3) Schwartz, B. J.; Bittner, E. R.; Prezhdo, O. V.; Rossky, P. J. J. Chem. Phys. 1996, 104, 5942. (4) Jasper, A. W.; Stechmann, S. N.; Truhlar, D. G. J. Chem. Phys. 2002, 116, 5424−5431. (5) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. J. Chem. Phys. 2005, 123, 234106. (6) Kapral, R.; Ciccotti, G. J. Chem. Phys. 1999, 110, 8919. (7) Landry, B. R.; Subotnik, J. E. J. Chem. Phys. 2012, 137, 22A513. (8) Wang, L.; Trivedi, D.; Prezhdo, O. V. J. Chem. Theory Comput. 2014, 10, 3598. (9) Kapral, R. Annu. Rev. Phys. Chem. 2006, 57, 129. (10) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. J. Chem. Phys. 2012, 137, 22A545. (11) Wang, L.; Prezhdo, O. V. J. Phys. Chem. Lett. 2014, 5, 713. (12) Ouyang, W.; Dou, W.; Subotnik, J. E. J. Chem. Phys. 2015, 142, 084109. (13) Aleksandrov, I. V. Z. Naturforsch., A: Phys. Sci. 1981, 36, 902− 908. (14) Gerasimenko, V. I. Theor. Math. Phys. 1982, 50, 77−87. (15) MacKernan, D.; Kapral, R.; Ciccotti, G. J. Phys.: Condens. Matter 2002, 14, 9069−9076. (16) Hanna, G.; Kapral, R. J. Chem. Phys. 2005, 122, 244505. 485


Self-Consistent Filtering Scheme for Efficient Calculations of

Recommend Documents