- Email: [email protected]
Tensor propagator with weight-selected paths for quantum dissipative dynamics with long-memory kernels
2 February 1996
CHEMICAL
PHYSICS LE'r'rERS ELSEVIER
Chemical Physics Letters 249 (1996) 224-230
Tensor propagator with weight-selected paths for quantum dissipative dynamics with long-memory kernels Eunji Sim, Nancy Makri School of Chemical Sciences, University of Illinois, 505 S. Mathews Avenue, Urbana, IL 61801, USA Received 31 August 1995; in final form 22 November 1995
Abstract
The tensor multiplication scheme, which allows iterative evaluation of the path integral for quantum dissipative dynamics, is extended to processes involving long-memory kernels such as those occurring in low-frequency solvents. Rather than propagating the augmented reduced density tensor in full-dimensional space, we propose here including in the propagator tensor only those path segments which enter the path integral with appreciable weight. The appropriate path segments are selected according to their weights by importance sampling. Elimination of the vast majority of paths which make negligible contribution in the path integral results in dramatic reduction of the required storage.
1. Introduction
In recent years, much theoretical effort has been expended toward developing numerical methods for simulating the dynamics of quantum systems in dissipative environments. These are commonly modeled by harmonic baths coupled to the particle of interest and often provide a realistic description of condensed phase processes. Feynman's path integral formulation [1], in which harmonic degrees of freedom can be integrated out analytically, appears as the most promising approach to the dynamics of such systems. However, elimination of the harmonic bath from the path integral introduces an influence functional [2] which is nonlocal in time. As a consequence, the resulting path integral for the reduced density matrix cannot be mapped onto a differential equation and therefore cannot be evaluated by common matrix-vector multiplication methods which are widely used for wavefunction propagation in small molecules. With a finite time step, the path integral prescription amounts to evaluation of a multidimensional integral whose dimension increases with the total time. Global (multidimensional) evaluation of the path integral by Monte Carlo methods is not feasible for long real time due to the oscillatory nature of the integrand [3,4]. Our group has recently developed a successful alternative. Taking advantage of the finite memory length characteristic of environments with broad spectra, Makarov and Makri have shown [5,6] that the path integral can be broken into a series of successive low-dimensional operations. Although the evolution of the reduced density matrix is not Markovian in such cases, the finite length of nonlocal interactions in the influence functional implies that there exists a higher-dimensional object, a reduced density t e n s o r which obeys Markovian dynamics. The reduced density tensor has rank equal to the number Akma x of time steps necessary to span the memory length and is propagated iteratively by multiplication with an appropriate propagator tensor. 0009-26t4/96/$12.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0009-261 4 ( 9 5 ) 0 1 3 7 4 - 1
E. Sire,N. Makri/ ChemicalPhysicsLetters 249 (1996)224-230
225
The tensor multiplication scheme provides a viable means of calculating the quantum dynamics of dissipative systems if the memory length is not very long compared to the allowed time step. This is usually the case when the environment has a broad spectrum that extends to reasonably high frequencies. Following the dynamics in low-frequency solvents presents a greater challenge: the rank of the propagator tensor is (with a fixed time step) proportional to A kmax and therefore the required storage grows exponentially with the solvent memory time. In this Letter we propose a variant of the tensor multiplication scheme, which leads to very significant reduction of storage and thus enables iterafive calculation of the dynamics for processes occurring in low-frequency environments. Rather than including in the propagator tensor all conceivable path segments for a given value of Akmax (which, when combined through iteration give rise to all discretized paths connecting initial and final points) we propose taking into consideration only those path segments that contribute in the path integral with significant weight. These are selected from a Monte Carlo random walk employing a path of AkmaX time slices. While neglecting low-probability path segments does not affect the final result noticeably, the savings accomplished can be dramatic as the fraction of 'important' paths is usually extremely small. In Section 2 we describe the methodology that we propose. Section 3 illustrates this scheme with a numerical example involving the dynamics of a two-state system with parameters similar to those for primary charge separation in a bacterial photosynthetic reaction center. This process occurs in a low-frequency solvent with long-memory interactions [7]; yet, its dynamics can be obtained successfully with the present scheme. Finally, Section 4 presents some concluding remarks.
2. Tensor propagator with weight-selected paths The processes dealt with in this Letter are described by the following Hamiltonian: H=H°+~)
2mj P/
(1)
where H 0 is the Hamiitonian for the bare system, s is the system ('reaction') coordinate, and Qj are harmonic bath degrees of freedom linearly coupled to the system coordinate. The characteristics of the bath that are relevant to the dynamics of the system coordinate are captured in the spectral density function [8]
jfo,)=TEj
c j2
.a(o,-
mj(.oj
(2)
which for macroscopic environments is practically a continuous function. In order to obtain an accurate propagator we partition the Hamiltonian H into a reference part H 0 and a bath of harmonic oscillators whose equilibrium positions are adiabatically displaced along the 'adiabatic path'
Qj = cj s/mj o92, He. v = H - H0,
(3)
Employing a symmetric splitting of the time evolution operator based on the above adiabatic partitioning of the system-bath Hamiltonian leads to the following quasi-adiabatic propagator [9]: ( sk Q k [ e x p ( - i H A t / h ) Is,_ lQk-i ) -- ( sk [ e x p ( - i H o A t / h ) ISk_l ) (Qklexp[ - i H , . v (sk) At/2h]exp[ - i H~nv( sk_ ~) A t / 2 h ] IQk-~ ),
(4)
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 2 2 4 - 2 3 0
226
where A t - - t i N is the time step. The propagator for the adiabatically renormalized system Hamiltonian is calculated numerically in terms of the system eigenstates, r/max
( s k ~ e x p ( - i H o A t / h ) l s k - l ) = 12 exp( - i E . A t / h ) ( s k l ~ . ) ( q ~ . l s
k_ t)-
(5)
n=l
Here E. and qb are the eigenvalues and eigenfunctions of H 0 and the number nmax of states that must be included is determined primarily by the system eigenstates that contribute to the initial density matrix of the system. The propagators involving displaced harmonic oscillators are also known exactly; their parametric dependence on the reaction coordinate introduces Franck-Condon factors that account for nonadiabatic corrections to the exact system propagator of Eq. (5) [9]. If the interaction between system and bath is switched on at t = 0, i.e. if the initial density matrix has a product form, p(0) = ps(0) pb(0),
(6)
the evolution of the reduced density matrix of the system
~( s', s" ; t) - T r b [ ( s " l e x p ( - i H t / h ) p ( O ) exp( i H t / h ) ls' ) ]
(7)
is given by the quasi-adiabatic propagator path integral expression [10]
~(s',s",t) •
--
fdsfffds? ... f ds~ _,fdsofdsr ..fds;,_, -inoAt/h)ls~_ 1)... < s~lexp( - i n o A t / h ) l s f f ) . . . ×
× ( s"lexp(
X l ( s f f . s ~ ,..
. ,s~v_,,s",s . . o sI ....
(sffl p,(0) Is o )
- l, s ' ; A t ) , su_
(8)
where
I( s ~ , s ~ . . . .
, S N+_ | , s t t , S o
,S 1 ,
... s~v_t,s';At )
= Vrb{exp[ -- i n,.v ( s " ) A t / 2 h ] exp[ - i n¢.v ( s•_ 1) A t / h ] . . .
exp[ - i H,.~ ( sff ) A t / 2 h ]
× pb(O)exp[ine.v(So)at/Zt~] ...exp[inonv(S;,_~)at/h] exp[ino,~(s')at/2h]}
(9)
is the influence functional. Here the symbols {s~', s~ . . . . } and {so, s I . . . . } denote discretizations of the forward and backward paths with endpoints {s~, s"} and {so, s'}, respectively. In numerical calculations each integral in Eq. (8) is evaluated by an M-point quadrature. For this purpose these integrals are replaced by sums, bringing Eq. (8) into a form identical to that for a generalized spin-boson Hamiltonian (i.e. a discrete system defined in terms of a finite number of levels). We have found [11] the discrete variable representation (DVR) [12] that corresponds to H 0 and its time-dependent variants [13] to offer very economical quadratures for this purpose. If the bath is at temperature T = 1 / k s ~, the influence functional is given by the expression [2] 1=exp
-
I2 =
,
(10)
k'= 0
where S N÷ --- s" and s~ = s' and the coefficients ~kk' have been given in Ref. [10]. These coefficients constitute the discrete analog of the two-time kernel ~ ( t - t') in the path integral, where the function ~(t) is the bath response function. The tensor multiplication scheme exploits the observation that for condensed phase environments characterized by broad spectra the response function decays within a finite time interval. As a consequence, the
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
227
coefficients rlkk, rapidly approach zero as the difference Ik - k'l increases. The finite span of nonlocality implies that the inner sum in the influence functional can be truncated at some appropriate Ik - k'l = A kmax and that the multidimensional path integral can be broken into a series of lower-dimensional integrals. The evolution of the reduced density matrix at time t = NAt can then be expressed in terms of the rank-A kmax augmented reduced density tensor A through the projection /3(s~ ; N A t ) = a ' a k - ' ) ( s ~ , sff+, = " " = sN+ • a k _ , - , = o; Na,)10(s
).
(11)
The tensor A is propagated through a time increment A kmax A t according to the tensor multiplication [5,6] a(a'm~')( s±,+a,.,,, " • •, =
S~+aa,m.,-l;(k+Akmax) At)
1) f d s ~ . . , f d s ~ + a , . . _ , T (2a*"~')( Sk~.Sk~+l. . . . . Sk+2Akm.x_ *
×A(ak"~')(s~ . . . . . sk+ak.~ - l ; k A t ),
(12)
where each integral is in practice discretized on an M-dimensional DVR grid. Here T is a rank-2Akma x propagator tensor whose elements are given by the relation
T(2a'm~')(s~ , s~+ l . . . . . s~+2 a , ~ , - , ) k + Akmax - 1
~-"
N n=k
]o(S~)]l(Sn:t:,S~+l)12(SZ
'S~+2)'''[Akr~,(S~,S~+Akmax)
X K( s) , sff+ , )
(13)
(1
)
where the influence functional 'bonds' lak are defined as
lo( s~ ) = exp - ~ ( s~ - si )( rl, ks 2 - rl~ksi ) ,
(1+
I ~ , ( s ~ , s ~ + ~ ) = exp - ~ - ( s , + ~ k - s L ~ k ) ( n , + ~ k . k s ~ - n~'+~.,s~)
) ,
zx/c> 0.
(14)
In order to proceed, we point out that the above tensors can be 'unfolded', such that A (Akin) becomes a
vector of dimension (M2) ~ k ~ whose elements are the 2A/cm~x-dimensional arrays y~ that hold values of the coordinate configurations (s~-, sk+ + 1, .. ., S k + a k ~ - 1, Sk , Sk+ 1, . . . . S*-+ak~,,-l) and T (2ak~') becomes a matr/x of dimension (M2) 2akm~. Thus, Eq. (12) can be viewed as the matrix-vector multiplication
a( y; ( k + A/Cmax)At ) =
f d y' T( y', y) a( y' ; /cAt)
(15a)
or equivalently, in the DVR representation, (M2)2Akm~
a(y,;(k+Akma~)At)=
E
T(yj, y~)a(yj;kAt),
i = 1 . . . . . ( M 2 ) ak°~
(15b)
j=l
Each configuration Yi defines a path segment that spans a time length equal to A kmax A t. T h r o u g h iteration of Eq. (15) these path segments are combined in all possible ways to form, after N iterations, a Feynman path of length NAkmaxAt. Since in the original formulation of the tensor multiplication scheme all possible path segments Yi were included, the final result was completely equivalent to that of the full path integral with nonlocal interactions truncated at length equal to A kmax A t. One must recall, however, that not all paths contribute with the same weight in the path integral: The matrix elements of the (exact) system propagator have different magnitudes, and the presence of dissipation introduces significant damping such that the "important' paths constitute a very small fraction of the total number, allowing evaluation of the path integral by importance sampling for short time. Although Monte Carlo methods fail to
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
228
converge at long times due to phase cancellation, the fraction of path space corresponding to paths of negligible weight rapidly approaches unity as the total time increases. Taking advantage of this fact, we modify the tensor multiplication scheme by excluding from Eq. (15) all path segments whose weight in the path integral lies below a certain threshold. To select the important path segments, we perform a Monte Carlo random walk using as the sampling function the absolute value of the integrand in the path integral representation of Eq. (8) for N = Akmax - 1. Since the selected configurations must describe path segments near as well as away from the initial and total time, we do not constrain the endpoints in this calculation and omit the initial density matrix from the sampling function. At the end of this process repeated selections are discarded to yield R distinct path segments of length AkmaxAt which are used to construct the vector A. (The above selection procedure actually overestimates path weights, since it neglects inter-segment influence functional bonds which can only lower the weight of each full-length path.) Thus, the weight-selecting tensor multiplication scheme consists in the following iterative matrix-vector multiplication: (M2)2akm,
A(yi;(k+Akmax)At )=
~_,
T(yj, Y i ) A ( y j ; k A t ),
i = l . . . . . R.
(16)
j=l
Typically R >> (M 2)ak,~x for A kmax >~ 4, allowing propagation with far less computer storage compared to that required in full tensor multiplication. Even though this scheme employs importance sampling to select path segments, it differs from standard Monte Carlo path integral methods in a very fundamental aspect: As the selected path segments are combined in the multiplication process, using R path segments in the vector A will after m iterations produce a result equivalent to that obtained from the full path integral (with nonlocal interactions truncated at A kmax steps) with R N Monte Carlo points. Although the number of path segments that must be included for accurate results varies depending on the application, in most of our test calculations we were able to obtain converged results with R on the order of 105. After just 100 iterations (i.e. a calculation to total time equal to 100A kmax A t) our scheme will yield in such cases results equivalent to those of a Monte Carlo path integral calculation employing 10500 randomly selected points, a number clearly unapproachable by importance sampling. Finally, unlike the full tensor multiplication scheme, the modified scheme that uses weight-selected paths does not conserve the norm of the density matrix unless all non-negligible path segments have been included. We have found norm conservation to be an extremely useful criterion for convergence.
3. Numerical test
We demonstrate the ability of the weight-selecting tensor multiplication scheme to yield converged time evolution for processes involving long-range memory kernels by applying it to an asymmetric two-level system in contact with a sluggish dissipative bath. The energy levels, electronic couplings and dissipation parameter are similar to these for the initial step of electron transfer in bacterial photosynthesis [14-16] if the reduced accessory chlorophyll monomer is involved via a sequential mechanism. We chose this system as a test of the present scheme, as our group is presently attempting to perform long-time simulations in order to elucidate the mechanism of primary charge separation. The Hamiltonian has the form
H= with
V, 2
e2
+ . ~
c m - l , VI2 ~-~ 22
E l = O, ~2 =
--400
J(ta)
exp( - to/toe),
=
"rico
+ -~rnjto) ~ 2~j ^2+cj 0
-
Qj,
(17)
cm -~. The spectral density has the form (18)
1°!\
E. Sire, N. Makri / Chemical Physics Letters 249 (1996) 224-230
:~
0.75
~"
229
0.5
" x.
a7
, :::, 0.25
::~'..
0.0
L
, 7-
0.0
i
I
~
f--T~-
,
6.0
I
t
~
i
~-
12.0
18.0
t. ps
Fig. 1. Population of the two electronic states as a function of time obtained with the scheme presented in this paper for Akma x for various values of the number of path segments R. Chain-dotted line: R = 18 (threshold 0.01). Dotted line: R = 34 (threshold 0.001). Dashed line: R = 560 (threshold 10-4). Solid line: R = 2876 (threshold 10-5). Markers:R = 4261 (threshold 5 × 10-6).
with oJ¢ = 150 c m - 1 and the solvent reorganization energy is 2r/% = 500 c m - J. The initial state of the system is chosen as state 1. We have chosen path segments (spin configurations) according to the Monte Carlo procedure described in Section 2. Convergence was tested by lowering the threshold that determines which path segments are included in the dynamics. Fig. 1 shows the convergence of the scheme with respect to the number R of path segments included for A k m a x = 20. Shown there are the populations of the two states (i.e. the diagonal elements of the reduced density matrix) as a function of time. While the results are poor with only 34 path segments (corresponding to a threshold equal to 10-3), the calculation is essentially converged with 560 path segments (selected with threshold 10 -4) as can be seen from the comparison with results obtained with 2876 spin configurations (for a threshold 10-5). It is seen that the number of path segments required for convergence is indeed very small
:/:/
075
P".
0.5
0.25
7-7~F~7
0.0
,
0.0
,
7::::::: ~
i
r-~
~ ~ T-
,
,
6.0
12.0
18.0
t, ps Fig. 2. State populations as a function of time for different values of Akma x . Chain-dotted, dotted, dashed and solid lines and markers show results for Akma x = 5, 10, 15, 20 and 25, respectively.
230
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
compared to the total number o f available spin configurations which in the present case is 2 40 -- 1012. This rapid convergence is a consequence of the localization of the integrand in narrow parts of multidimensional space and is typical of strongly dissipative environments. In Fig. 2 we show the decay of the same populations as a function of time for different memory lengths corresponding to Akmax = 5, 10, 15, 20 and 25. The results shown for each value of Akmax are converged with respect to the number of path segments. Here the convergence is indeed slow because of the long range of solvent memory and calculation of the dynamics by full tensor multiplication would not have been feasible. For the longest memory length results the required number o f path segments was equal to 3606. Clearly, much longer arrays can easily be propagated on modem computers,
4. Concluding remarks In this Letter we have presented an efficient scheme for iterative evaluation of the real-time path integral in quantum dissipative systems. While the underlying idea is the same as in the earlier tensor multiplication scheme, the present method does not waste computer resources on paths that make negligible contribution to the dynamics. As a consequence, the required storage is reduced dramatically, making calculations feasible in cases where the solvent memory is very long. The numerical tests presented in Section 3 indicate that weight selection of paths can enable iterative calculation of the dynamics with full memory even if the latter spans many time steps. If the integrand is indeed localized, as is the case in strongly dissipative systems, the fraction of statistically significant paths is expected to decrease with increasing Akmax, resulting in slow increase of the required numerical effort with solvent memory length. We are currently using this scheme to investigate the dynamics of primary charge separation in bacterial photosynthetic reaction centers.
Acknowledgements This work has been supported by a National Science Foundation Young Investigator Award and by a David and Lucille Packard Foundation for Science and Engineering.
References [1] R.P. Feynman, Rev. Mod. Phys. 20 (1948) 367. [2] R.P. Feynman and F.L. Vernon Jr., Ann. Phys. 118 (1963). [3] J.D. Doll, D.L. Freeman and T.L. Beck, Advan. Chem. Phys. 78 (1990) 61. [4] N. Makri, Comput. Phys. Commun. 63 (1991) 389. [5] D.E. Makarov and N. Makri, Chem. Phys. Letters 221 (1994) 482. [6] N. Makri and D.E. Makarov, J. Chem. Phys. 102 (1995) 4600, 4611. [7] M. Marchi, J.N. Gehlen, D. Chandler and M. Newton, J. Amer. Chem. Soc. 115 (1993) 4178. [8] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P. Fisher, A. Garg and W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [9] N. Makri, Chem. Phys. Letters 193 (1992) 435. [10] N. Makri, J. Math. Phys. 36 (1995) 30. [11] M. Topaler and N. Makri, Chem. Phys. Letters 210 (1994) 448. [12] J.V. Lill, G.A. Parker and J.C. Light, Chem. Phys. Letters 89 (1982) 483; Z. Bacic and J.C. Light, Ann. Rev. Phys. Chem. 40 (1989) 469. [13] E. Sire and N. Makri, J. Chem. Phys. 102 (1995) 5616. [14] G.R. Fleming, J.L. Martin and J. Breton, Nature 333 (1988) 190. [15] C.-K. Chart, T.J. DiMagno, L.X.-Q. Chert, J.R. Norris and G.R. Fleming, Proc. Natl. Acad. Sci. US 88 (1991) 11202. [16] S. Schmidt, T. Arlt, P. Harem, H. Huber, T. Nagele, J. Wachtveitl, M. Meyer, H. Scheer and W. Zinth, Chem. Phys. Letters 223 (1994) 116.
CHEMICAL
PHYSICS LE'r'rERS ELSEVIER
Chemical Physics Letters 249 (1996) 224-230
Tensor propagator with weight-selected paths for quantum dissipative dynamics with long-memory kernels Eunji Sim, Nancy Makri School of Chemical Sciences, University of Illinois, 505 S. Mathews Avenue, Urbana, IL 61801, USA Received 31 August 1995; in final form 22 November 1995
Abstract
The tensor multiplication scheme, which allows iterative evaluation of the path integral for quantum dissipative dynamics, is extended to processes involving long-memory kernels such as those occurring in low-frequency solvents. Rather than propagating the augmented reduced density tensor in full-dimensional space, we propose here including in the propagator tensor only those path segments which enter the path integral with appreciable weight. The appropriate path segments are selected according to their weights by importance sampling. Elimination of the vast majority of paths which make negligible contribution in the path integral results in dramatic reduction of the required storage.
1. Introduction
In recent years, much theoretical effort has been expended toward developing numerical methods for simulating the dynamics of quantum systems in dissipative environments. These are commonly modeled by harmonic baths coupled to the particle of interest and often provide a realistic description of condensed phase processes. Feynman's path integral formulation [1], in which harmonic degrees of freedom can be integrated out analytically, appears as the most promising approach to the dynamics of such systems. However, elimination of the harmonic bath from the path integral introduces an influence functional [2] which is nonlocal in time. As a consequence, the resulting path integral for the reduced density matrix cannot be mapped onto a differential equation and therefore cannot be evaluated by common matrix-vector multiplication methods which are widely used for wavefunction propagation in small molecules. With a finite time step, the path integral prescription amounts to evaluation of a multidimensional integral whose dimension increases with the total time. Global (multidimensional) evaluation of the path integral by Monte Carlo methods is not feasible for long real time due to the oscillatory nature of the integrand [3,4]. Our group has recently developed a successful alternative. Taking advantage of the finite memory length characteristic of environments with broad spectra, Makarov and Makri have shown [5,6] that the path integral can be broken into a series of successive low-dimensional operations. Although the evolution of the reduced density matrix is not Markovian in such cases, the finite length of nonlocal interactions in the influence functional implies that there exists a higher-dimensional object, a reduced density t e n s o r which obeys Markovian dynamics. The reduced density tensor has rank equal to the number Akma x of time steps necessary to span the memory length and is propagated iteratively by multiplication with an appropriate propagator tensor. 0009-26t4/96/$12.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0009-261 4 ( 9 5 ) 0 1 3 7 4 - 1
E. Sire,N. Makri/ ChemicalPhysicsLetters 249 (1996)224-230
225
The tensor multiplication scheme provides a viable means of calculating the quantum dynamics of dissipative systems if the memory length is not very long compared to the allowed time step. This is usually the case when the environment has a broad spectrum that extends to reasonably high frequencies. Following the dynamics in low-frequency solvents presents a greater challenge: the rank of the propagator tensor is (with a fixed time step) proportional to A kmax and therefore the required storage grows exponentially with the solvent memory time. In this Letter we propose a variant of the tensor multiplication scheme, which leads to very significant reduction of storage and thus enables iterafive calculation of the dynamics for processes occurring in low-frequency environments. Rather than including in the propagator tensor all conceivable path segments for a given value of Akmax (which, when combined through iteration give rise to all discretized paths connecting initial and final points) we propose taking into consideration only those path segments that contribute in the path integral with significant weight. These are selected from a Monte Carlo random walk employing a path of AkmaX time slices. While neglecting low-probability path segments does not affect the final result noticeably, the savings accomplished can be dramatic as the fraction of 'important' paths is usually extremely small. In Section 2 we describe the methodology that we propose. Section 3 illustrates this scheme with a numerical example involving the dynamics of a two-state system with parameters similar to those for primary charge separation in a bacterial photosynthetic reaction center. This process occurs in a low-frequency solvent with long-memory interactions [7]; yet, its dynamics can be obtained successfully with the present scheme. Finally, Section 4 presents some concluding remarks.
2. Tensor propagator with weight-selected paths The processes dealt with in this Letter are described by the following Hamiltonian: H=H°+~)
2mj P/
(1)
where H 0 is the Hamiitonian for the bare system, s is the system ('reaction') coordinate, and Qj are harmonic bath degrees of freedom linearly coupled to the system coordinate. The characteristics of the bath that are relevant to the dynamics of the system coordinate are captured in the spectral density function [8]
jfo,)=TEj
c j2
.a(o,-
mj(.oj
(2)
which for macroscopic environments is practically a continuous function. In order to obtain an accurate propagator we partition the Hamiltonian H into a reference part H 0 and a bath of harmonic oscillators whose equilibrium positions are adiabatically displaced along the 'adiabatic path'
Qj = cj s/mj o92, He. v = H - H0,
(3)
Employing a symmetric splitting of the time evolution operator based on the above adiabatic partitioning of the system-bath Hamiltonian leads to the following quasi-adiabatic propagator [9]: ( sk Q k [ e x p ( - i H A t / h ) Is,_ lQk-i ) -- ( sk [ e x p ( - i H o A t / h ) ISk_l ) (Qklexp[ - i H , . v (sk) At/2h]exp[ - i H~nv( sk_ ~) A t / 2 h ] IQk-~ ),
(4)
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 2 2 4 - 2 3 0
226
where A t - - t i N is the time step. The propagator for the adiabatically renormalized system Hamiltonian is calculated numerically in terms of the system eigenstates, r/max
( s k ~ e x p ( - i H o A t / h ) l s k - l ) = 12 exp( - i E . A t / h ) ( s k l ~ . ) ( q ~ . l s
k_ t)-
(5)
n=l
Here E. and qb are the eigenvalues and eigenfunctions of H 0 and the number nmax of states that must be included is determined primarily by the system eigenstates that contribute to the initial density matrix of the system. The propagators involving displaced harmonic oscillators are also known exactly; their parametric dependence on the reaction coordinate introduces Franck-Condon factors that account for nonadiabatic corrections to the exact system propagator of Eq. (5) [9]. If the interaction between system and bath is switched on at t = 0, i.e. if the initial density matrix has a product form, p(0) = ps(0) pb(0),
(6)
the evolution of the reduced density matrix of the system
~( s', s" ; t) - T r b [ ( s " l e x p ( - i H t / h ) p ( O ) exp( i H t / h ) ls' ) ]
(7)
is given by the quasi-adiabatic propagator path integral expression [10]
~(s',s",t) •
--
fdsfffds? ... f ds~ _,fdsofdsr ..fds;,_, -inoAt/h)ls~_ 1)... < s~lexp( - i n o A t / h ) l s f f ) . . . ×
× ( s"lexp(
X l ( s f f . s ~ ,..
. ,s~v_,,s",s . . o sI ....
(sffl p,(0) Is o )
- l, s ' ; A t ) , su_
(8)
where
I( s ~ , s ~ . . . .
, S N+_ | , s t t , S o
,S 1 ,
... s~v_t,s';At )
= Vrb{exp[ -- i n,.v ( s " ) A t / 2 h ] exp[ - i n¢.v ( s•_ 1) A t / h ] . . .
exp[ - i H,.~ ( sff ) A t / 2 h ]
× pb(O)exp[ine.v(So)at/Zt~] ...exp[inonv(S;,_~)at/h] exp[ino,~(s')at/2h]}
(9)
is the influence functional. Here the symbols {s~', s~ . . . . } and {so, s I . . . . } denote discretizations of the forward and backward paths with endpoints {s~, s"} and {so, s'}, respectively. In numerical calculations each integral in Eq. (8) is evaluated by an M-point quadrature. For this purpose these integrals are replaced by sums, bringing Eq. (8) into a form identical to that for a generalized spin-boson Hamiltonian (i.e. a discrete system defined in terms of a finite number of levels). We have found [11] the discrete variable representation (DVR) [12] that corresponds to H 0 and its time-dependent variants [13] to offer very economical quadratures for this purpose. If the bath is at temperature T = 1 / k s ~, the influence functional is given by the expression [2] 1=exp
-
I2 =
,
(10)
k'= 0
where S N÷ --- s" and s~ = s' and the coefficients ~kk' have been given in Ref. [10]. These coefficients constitute the discrete analog of the two-time kernel ~ ( t - t') in the path integral, where the function ~(t) is the bath response function. The tensor multiplication scheme exploits the observation that for condensed phase environments characterized by broad spectra the response function decays within a finite time interval. As a consequence, the
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
227
coefficients rlkk, rapidly approach zero as the difference Ik - k'l increases. The finite span of nonlocality implies that the inner sum in the influence functional can be truncated at some appropriate Ik - k'l = A kmax and that the multidimensional path integral can be broken into a series of lower-dimensional integrals. The evolution of the reduced density matrix at time t = NAt can then be expressed in terms of the rank-A kmax augmented reduced density tensor A through the projection /3(s~ ; N A t ) = a ' a k - ' ) ( s ~ , sff+, = " " = sN+ • a k _ , - , = o; Na,)10(s
).
(11)
The tensor A is propagated through a time increment A kmax A t according to the tensor multiplication [5,6] a(a'm~')( s±,+a,.,,, " • •, =
S~+aa,m.,-l;(k+Akmax) At)
1) f d s ~ . . , f d s ~ + a , . . _ , T (2a*"~')( Sk~.Sk~+l. . . . . Sk+2Akm.x_ *
×A(ak"~')(s~ . . . . . sk+ak.~ - l ; k A t ),
(12)
where each integral is in practice discretized on an M-dimensional DVR grid. Here T is a rank-2Akma x propagator tensor whose elements are given by the relation
T(2a'm~')(s~ , s~+ l . . . . . s~+2 a , ~ , - , ) k + Akmax - 1
~-"
N n=k
]o(S~)]l(Sn:t:,S~+l)12(SZ
'S~+2)'''[Akr~,(S~,S~+Akmax)
X K( s) , sff+ , )
(13)
(1
)
where the influence functional 'bonds' lak are defined as
lo( s~ ) = exp - ~ ( s~ - si )( rl, ks 2 - rl~ksi ) ,
(1+
I ~ , ( s ~ , s ~ + ~ ) = exp - ~ - ( s , + ~ k - s L ~ k ) ( n , + ~ k . k s ~ - n~'+~.,s~)
) ,
zx/c> 0.
(14)
In order to proceed, we point out that the above tensors can be 'unfolded', such that A (Akin) becomes a
vector of dimension (M2) ~ k ~ whose elements are the 2A/cm~x-dimensional arrays y~ that hold values of the coordinate configurations (s~-, sk+ + 1, .. ., S k + a k ~ - 1, Sk , Sk+ 1, . . . . S*-+ak~,,-l) and T (2ak~') becomes a matr/x of dimension (M2) 2akm~. Thus, Eq. (12) can be viewed as the matrix-vector multiplication
a( y; ( k + A/Cmax)At ) =
f d y' T( y', y) a( y' ; /cAt)
(15a)
or equivalently, in the DVR representation, (M2)2Akm~
a(y,;(k+Akma~)At)=
E
T(yj, y~)a(yj;kAt),
i = 1 . . . . . ( M 2 ) ak°~
(15b)
j=l
Each configuration Yi defines a path segment that spans a time length equal to A kmax A t. T h r o u g h iteration of Eq. (15) these path segments are combined in all possible ways to form, after N iterations, a Feynman path of length NAkmaxAt. Since in the original formulation of the tensor multiplication scheme all possible path segments Yi were included, the final result was completely equivalent to that of the full path integral with nonlocal interactions truncated at length equal to A kmax A t. One must recall, however, that not all paths contribute with the same weight in the path integral: The matrix elements of the (exact) system propagator have different magnitudes, and the presence of dissipation introduces significant damping such that the "important' paths constitute a very small fraction of the total number, allowing evaluation of the path integral by importance sampling for short time. Although Monte Carlo methods fail to
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
228
converge at long times due to phase cancellation, the fraction of path space corresponding to paths of negligible weight rapidly approaches unity as the total time increases. Taking advantage of this fact, we modify the tensor multiplication scheme by excluding from Eq. (15) all path segments whose weight in the path integral lies below a certain threshold. To select the important path segments, we perform a Monte Carlo random walk using as the sampling function the absolute value of the integrand in the path integral representation of Eq. (8) for N = Akmax - 1. Since the selected configurations must describe path segments near as well as away from the initial and total time, we do not constrain the endpoints in this calculation and omit the initial density matrix from the sampling function. At the end of this process repeated selections are discarded to yield R distinct path segments of length AkmaxAt which are used to construct the vector A. (The above selection procedure actually overestimates path weights, since it neglects inter-segment influence functional bonds which can only lower the weight of each full-length path.) Thus, the weight-selecting tensor multiplication scheme consists in the following iterative matrix-vector multiplication: (M2)2akm,
A(yi;(k+Akmax)At )=
~_,
T(yj, Y i ) A ( y j ; k A t ),
i = l . . . . . R.
(16)
j=l
Typically R >> (M 2)ak,~x for A kmax >~ 4, allowing propagation with far less computer storage compared to that required in full tensor multiplication. Even though this scheme employs importance sampling to select path segments, it differs from standard Monte Carlo path integral methods in a very fundamental aspect: As the selected path segments are combined in the multiplication process, using R path segments in the vector A will after m iterations produce a result equivalent to that obtained from the full path integral (with nonlocal interactions truncated at A kmax steps) with R N Monte Carlo points. Although the number of path segments that must be included for accurate results varies depending on the application, in most of our test calculations we were able to obtain converged results with R on the order of 105. After just 100 iterations (i.e. a calculation to total time equal to 100A kmax A t) our scheme will yield in such cases results equivalent to those of a Monte Carlo path integral calculation employing 10500 randomly selected points, a number clearly unapproachable by importance sampling. Finally, unlike the full tensor multiplication scheme, the modified scheme that uses weight-selected paths does not conserve the norm of the density matrix unless all non-negligible path segments have been included. We have found norm conservation to be an extremely useful criterion for convergence.
3. Numerical test
We demonstrate the ability of the weight-selecting tensor multiplication scheme to yield converged time evolution for processes involving long-range memory kernels by applying it to an asymmetric two-level system in contact with a sluggish dissipative bath. The energy levels, electronic couplings and dissipation parameter are similar to these for the initial step of electron transfer in bacterial photosynthesis [14-16] if the reduced accessory chlorophyll monomer is involved via a sequential mechanism. We chose this system as a test of the present scheme, as our group is presently attempting to perform long-time simulations in order to elucidate the mechanism of primary charge separation. The Hamiltonian has the form
H= with
V, 2
e2
+ . ~
c m - l , VI2 ~-~ 22
E l = O, ~2 =
--400
J(ta)
exp( - to/toe),
=
"rico
+ -~rnjto) ~ 2~j ^2+cj 0
-
Qj,
(17)
cm -~. The spectral density has the form (18)
1°!\
E. Sire, N. Makri / Chemical Physics Letters 249 (1996) 224-230
:~
0.75
~"
229
0.5
" x.
a7
, :::, 0.25
::~'..
0.0
L
, 7-
0.0
i
I
~
f--T~-
,
6.0
I
t
~
i
~-
12.0
18.0
t. ps
Fig. 1. Population of the two electronic states as a function of time obtained with the scheme presented in this paper for Akma x for various values of the number of path segments R. Chain-dotted line: R = 18 (threshold 0.01). Dotted line: R = 34 (threshold 0.001). Dashed line: R = 560 (threshold 10-4). Solid line: R = 2876 (threshold 10-5). Markers:R = 4261 (threshold 5 × 10-6).
with oJ¢ = 150 c m - 1 and the solvent reorganization energy is 2r/% = 500 c m - J. The initial state of the system is chosen as state 1. We have chosen path segments (spin configurations) according to the Monte Carlo procedure described in Section 2. Convergence was tested by lowering the threshold that determines which path segments are included in the dynamics. Fig. 1 shows the convergence of the scheme with respect to the number R of path segments included for A k m a x = 20. Shown there are the populations of the two states (i.e. the diagonal elements of the reduced density matrix) as a function of time. While the results are poor with only 34 path segments (corresponding to a threshold equal to 10-3), the calculation is essentially converged with 560 path segments (selected with threshold 10 -4) as can be seen from the comparison with results obtained with 2876 spin configurations (for a threshold 10-5). It is seen that the number of path segments required for convergence is indeed very small
:/:/
075
P".
0.5
0.25
7-7~F~7
0.0
,
0.0
,
7::::::: ~
i
r-~
~ ~ T-
,
,
6.0
12.0
18.0
t, ps Fig. 2. State populations as a function of time for different values of Akma x . Chain-dotted, dotted, dashed and solid lines and markers show results for Akma x = 5, 10, 15, 20 and 25, respectively.
230
E. Sim, N. Makri / Chemical Physics Letters 249 (1996) 224-230
compared to the total number o f available spin configurations which in the present case is 2 40 -- 1012. This rapid convergence is a consequence of the localization of the integrand in narrow parts of multidimensional space and is typical of strongly dissipative environments. In Fig. 2 we show the decay of the same populations as a function of time for different memory lengths corresponding to Akmax = 5, 10, 15, 20 and 25. The results shown for each value of Akmax are converged with respect to the number of path segments. Here the convergence is indeed slow because of the long range of solvent memory and calculation of the dynamics by full tensor multiplication would not have been feasible. For the longest memory length results the required number o f path segments was equal to 3606. Clearly, much longer arrays can easily be propagated on modem computers,
4. Concluding remarks In this Letter we have presented an efficient scheme for iterative evaluation of the real-time path integral in quantum dissipative systems. While the underlying idea is the same as in the earlier tensor multiplication scheme, the present method does not waste computer resources on paths that make negligible contribution to the dynamics. As a consequence, the required storage is reduced dramatically, making calculations feasible in cases where the solvent memory is very long. The numerical tests presented in Section 3 indicate that weight selection of paths can enable iterative calculation of the dynamics with full memory even if the latter spans many time steps. If the integrand is indeed localized, as is the case in strongly dissipative systems, the fraction of statistically significant paths is expected to decrease with increasing Akmax, resulting in slow increase of the required numerical effort with solvent memory length. We are currently using this scheme to investigate the dynamics of primary charge separation in bacterial photosynthetic reaction centers.
Acknowledgements This work has been supported by a National Science Foundation Young Investigator Award and by a David and Lucille Packard Foundation for Science and Engineering.
References [1] R.P. Feynman, Rev. Mod. Phys. 20 (1948) 367. [2] R.P. Feynman and F.L. Vernon Jr., Ann. Phys. 118 (1963). [3] J.D. Doll, D.L. Freeman and T.L. Beck, Advan. Chem. Phys. 78 (1990) 61. [4] N. Makri, Comput. Phys. Commun. 63 (1991) 389. [5] D.E. Makarov and N. Makri, Chem. Phys. Letters 221 (1994) 482. [6] N. Makri and D.E. Makarov, J. Chem. Phys. 102 (1995) 4600, 4611. [7] M. Marchi, J.N. Gehlen, D. Chandler and M. Newton, J. Amer. Chem. Soc. 115 (1993) 4178. [8] A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P. Fisher, A. Garg and W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [9] N. Makri, Chem. Phys. Letters 193 (1992) 435. [10] N. Makri, J. Math. Phys. 36 (1995) 30. [11] M. Topaler and N. Makri, Chem. Phys. Letters 210 (1994) 448. [12] J.V. Lill, G.A. Parker and J.C. Light, Chem. Phys. Letters 89 (1982) 483; Z. Bacic and J.C. Light, Ann. Rev. Phys. Chem. 40 (1989) 469. [13] E. Sire and N. Makri, J. Chem. Phys. 102 (1995) 5616. [14] G.R. Fleming, J.L. Martin and J. Breton, Nature 333 (1988) 190. [15] C.-K. Chart, T.J. DiMagno, L.X.-Q. Chert, J.R. Norris and G.R. Fleming, Proc. Natl. Acad. Sci. US 88 (1991) 11202. [16] S. Schmidt, T. Arlt, P. Harem, H. Huber, T. Nagele, J. Wachtveitl, M. Meyer, H. Scheer and W. Zinth, Chem. Phys. Letters 223 (1994) 116.