Quantum dynamics simulation of nonadiabatic processes in complex systems is vital to the understanding of many photoinduced chemical reactions [1][2][3][4][5][6][7][8][9][10] as well as charge and energy transfer processes in solar energy conversion materials and photosynthesis. [11][12][13][14][15][16] While directly integrating the time-dependent Schrödinger equation (TDSE) offers the most accurate and detailed results, this approach is not practical for the simulation of complicated molecular systems due to the computational cost that grows exponentially with the number of degrees of freedom (DOF). The various computational approaches developed to circumvent the exponential scaling problems and enable the simulation of nonadiabatic dynamics in condensed phases can be roughly sorted into two categories: those that involve physical approximations and those that do not. Methods that do not involve physical approximation, henceforth referred to as the "numerically exact" methods, integrate the TDSE or the quantum Liouville equation (QLE) with costeffective numerical techniques that converge to the true answer with increased computational investment. 3,[17][18][19][20][21][22] For instance, techniques based on tensor networks 3,20,[23][24][25] or the Feynman-Vernon influence functional 21,22,[26][27][28][29][30][31] have been extensively utilized to study molecular systems. Despite their accuracy upon convergence, numerically exact methods are, in general, still quite expensive and are typically restricted to systems with a few tens to a few hundreds of DOF.
On the other hand, the approximate methods 2,[32][33][34][35] reduce the computational cost by invoking a classical-like trajectory description for part of the system, most typically the nuclear DOF, and are often able to treat systems with 10 3 -10 5 DOF. [36][37][38][39][40][41][42] These methods range from traditional mixed quantum-classical (MQC) methods, such as mean-field Ehrenfest and fewest-switches surface hopping (FSSH) methods, to a variety of nonadiabatic semiclassical dynamics. While approximate methods simulate the time-evolution of significantly more complicated systems, the accuracy for the description of quantum motion with classical trajectories remains an open question, and the answer might not only differ by the specific approximate methods but also depend on the characters of the simulated system, such as the shape of the potential energy surfaces (PESs), the type of electronic couplings, and the simulated temperature.
Recently, a family of nonadiabatic semiclassical dynamical methods with the mapping Hamiltonian has emerged, which provides a balanced performance in accuracy and computational efficiency. The Meyer-Miller Stock-Thoss (MMST) mapping Hamiltonian 32,34 provides a mapping between F discrete electronic states and a collection of F singly excited harmonic oscillators as a generalization of Schwinger's bosonization, 43,44 where the positions and momenta of these oscillators serve as Cartesian mapping variables for further trajectory-based semiclassical approximations. For example, the linearized semiclassical (LSC) dynamics can be derived by performing linearization approximation to the full semiclassical initial value representation (SC-IVR) as shown by Miller and co-workers, and the LSC method refrains from evaluating the computationally expensive stability matrix as required by SC-IVR. 35,45 An alternative way of arriving at LSC shown by Shi and Geva is to linearize the forward-backward action in the real-time pathintegral formulation 46 and thus can also be called the linearized path integral approach. [47][48][49] It is noted that the LSC expression for the time correlation function is identical to the expression obtained from the classical Wigner approach, featuring Wigner transformation of nuclear density used for the initial sampling of nuclear DOF. 46 The LSC strategy was demonstrated in many condensedphase systems to be effective in capturing the important quantum effects, for example, in the processes involving electronic transition and vibrational relaxation. 46,50 It can be shown that combining the quantum-classical Liouville equation (QCLE) with the mapping Hamiltonian generates the Poisson-bracket mapping equation (PBME) approximation, which is another version of the semiclassical mapping dynamics. 51 Recently, Saller et al. 42 and Gao et al. 52 have demonstrated that the LSC mapping dynamics can be combined with the resolution of identity operator to have better population dynamics in model systems.
The symmetrical quasiclassical (SQC) method proposed by Cotton and Miller 37 is also based on the MMST mapping Hamiltonian, but it utilizes a windowing approach for determining the electronic population. The idea of the SQC method originated from the Bohr-Sommerfeld quantization, where quantum states are associated with integer values of action variables. 53 In SQC dynamics, the positions and momenta mapping variables of the electronic DOF are transformed into the action-angle (a-a) variables, which evolve together with the nuclear DOF classically. The electronic population is obtained by filtering the action variables with a finite-width symmetrical window function, approximating the delta function. The square window function 53 serves as a coarse-grained filter for quantizing the populations of different electronic states, whose ensemble average over many trajectories provides continuous evolution of the electronic populations. The width of the window function is controlled by the zero-point energy (ZPE) parameter γ (the full quantum γ = 0.5), whose value can be fixed such as γ = ( √ 3 -1)2 ≈ 0.366 or adjusted on a per-trajectory basis such that the initial nuclear forces correspond to the potential energy surface (PES) of the initial electronic state. 54 Recently, the triangle window function (with optimal γ = 13) has been reported to give superior nonadiabatic dynamics than the original square window function for the weak coupling case. 38 SQC dynamics can also be applied to the spin-mapping Hamiltonian as proposed by Cotton and Miller, 55 where each electronic state is represented as a spin-1/2 DOF with spin-up and spin-down corresponding to the occupied and unoccupied states, and the creation and annihilation operators of electronic state are represented as the angular momentum ladder operators.
There are other mapping strategies for representing discrete electronic states, for instance, using a spin-S (such that 2S + 1 = F) to represent an F-level system was suggested by Meyer and Miller in 1979. 56 Recently, a general spin-mapping (SPM) approach has been proposed by Runeson and Richardson in 2019, which was formulated by using the Stratonovich-Weyl transform for a spin-1/2 for the two-level system 40 and a generalized F-level system that preserves the SU(F) symmetry. 41,57,58 Using the P, Q, and W representations of the Stratonovich-Weyl transform for the spin operators, the corresponding SPM approaches require specific ZPE parameters and hypersphere radius for the full-sphere sampling of the mapping variables and the W scheme is the recommended version of SPM. 41 As it turns out, the P, Q, and W schemes of SPM are equivalent to a special set of ZPE parameter choices in the extended classical mapping models (eCMMs) proposed by He and Liu in 2019, which are based on the Wigner-Weyl transform of the phase-space representation of quantum mechanics and a constraint sphere for the initial mapping phase-space variables, and the resulting ZPE parameters can be any values in γ ∈ (-1F, ∞). 59 It is noted that s may give accurate population dynamics for some systems when the ZPE parameter takes a negative value, such as γ = -0.2. The SPM/eCMM approaches are based on the initial sampling of the mapping variables on the full-sphere that reinforce the total population to be unity, but in some cases, the initial nuclear forces might cause unphysical nuclear dynamics due to the uniform sampling of all electronic states and thus their PESs. 14 To circumvent this issue, techniques ensuring the initial effective forces to be evaluated on the initial state's PES are developed like adjusting ZPE per trajectory by auxiliary commutator variables for eCMM 60 or focused mapping variable sampling for SPM. 41 However, the ZPE adjusting might not fully resolve the inverted potential issue for the mapping dynamics. 54 It is important to mark that the equations of motion for these spin-mapping approaches can be formulated to be the same in terms of Cartesian mapping variables as those for the LSC mapping dynamics defined in the MMST Hamiltonian, and the major differences between these mapping dynamics are the ZPE parameter, initial sampling, and observable evaluation. 33,34,61 In addition, the spin-mapping and phase-space-mapping strategies solve the problem of MMST mapping variables deviating from the singly excited oscillator (SEO) subspace, hence giving a better description of the identity operator and thus population dynamics. Although these
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp mapping relations would be exact for the electronic dynamics in the frozen-nuclei limit, 62 their accuracy has to be tested when applied to trajectory-based nonadiabatic dynamics with classical nuclear DOF since the coupled electronic-nuclear dynamics would be approximate. A detailed test on these different mapping dynamics as well as mixed quantum-classical dynamics is lacking, especially with more challenging models, ranging from spin-boson models to systems with anharmonic potentials and non-Condon electronic couplings.
To this end, we report a benchmark study for various approximate nonadiabatic mapping dynamics against the numerically exact quantum dynamics provided by tensor-train (TT) based methods. The tensor-train format for high-dimensional tensors provides a practical way for evolving the wavepacket of quantum systems with high-dimensional nuclear and electronic DOF, mitigating the exponential scaling problem, the so-called curse of dimensionality. In the TT format or the matrix product state (MPS), an Nth order wavefunction tensor can be approximated as a chain of N core tensors of order up to 3, and each core tensor corresponds to one physical dimension. A traditional split-operator Fourier transform (SOFT) method for wavepacket dynamics utilizes the Trotter expansion of the time-evolution operator and avoids implementing the nonlocal kinetic energy operator in the position space by Fourier transforming to the momentum space for kinetic energy propagation, but SOFT can handle very few DOF. Combining the TT format and SOFT propagation method, Greene and Batista proposed the tensor-train split-operator Fourier transform (TT-SOFT) method for direct nonadiabatic wavepacket propagation in highdimensional systems. 3 In TT-SOFT, the propagation involves the scaling and squaring that generate high-rank TT, and then, the rounding operation gives the low-rank approximation, but it may cause unstable norm conservation when used with limited rank. A more efficient way to propagate quantum dynamics is the TT implementation of the dynamical low-rank approximation, called the TT-KSL method 20,63 or the time-dependent variational principle (TDVP) with projector splitting, 25,[64][65][66] which stays in the fixed rank manifold by using projection onto the tangent space, giving rise to an optimal and efficient solution of propagation in the low-rank subspace. Recently, Lyu et al. have proposed tensor-train split-operator KSL (TT-SOKSL) 20 that takes advantage of the efficient propagation of TT-KSL on the low-rank manifold and also the local kinetic energy operator that only requires vector-vector multiplications in the momentum space as in SOFT. TT-SOKSL was shown to give better norm conservation than TT-SOFT in a 25-dimensional anharmonic and non-Condon model of retinal photoisomerization. 20 Solving the TDSE using the above-mentioned TT-based methods could only provide wavepacket dynamics at zero temperature, and thermofield dynamics (TFD) provides a way for solving the thermal Schrödinger equation that describes the thermal wavepacket defined in the double space. The combination of TT and TFD would allow for quantum dynamical treatment for multi-dimensional systems, and in particular, for harmonic systems, thermal Bogoliubov transformation enables an analytical expression for the thermal state in TT-TFD. 67 In addition, TT-TFD has been applied to compute finite-temperature quantum dynamics of spinboson models and the Fenna-Matthews-Olson (FMO) complex 68 as well as the kernels in the Nakajima-Zwanzig generalized quantum master equation. 69 In this work, we aim to test a broad range of models with qualitatively contrasting features (i.e., harmonic vs anharmonic PESs, Condon vs non-Condon electronic coupling, and zero temperature vs finite temperature) as well as quantitatively different sets of parameters (i.e., strong vs weak reorganization energies and biased vs unbiased electronic gap) to facilitate a comprehensive comparison among different application scenarios. We envision this work to provide practical recommendations for the choice of approximate mapping dynamics methods in simulating nonadiabatic dynamics in complex systems.
The rest of this paper is organized as follows. Section II gives an overview of all approximate and numerically exact methods employed in this work. Section III presents the benchmarked models. Section IV describes the simulation details. Section V reports and discusses the results. Section VI provides the concluding remarks.
We summarize the semiclassical and quasiclassical mapping dynamical methods that will be benchmarked. Consider a general multi-state system with F electronic states and N total nuclear DOF, whose total Hamiltonian in the diabatic basis is written as
where P = P1 , . . . , PN and R = R1 , . . . , RN denote the N massweighted nuclear coordinates and momenta; j (j = 1, . . . , F) represent the jth electronic state; V j j R ≡ V j R is the PES of the jth electronic state; V jk R is the electronic coupling between the jth and the kth electronic states (j ≠ k); and V( R) = ∑ F j,k V jk R jk is the electronic Hamiltonian. Within the Condon approximation, V jk R becomes constant V jk = Γ jk (j ≠ k) and V jk = V kj . The MMST mapping Hamiltonian is based on the mapping from the multi-state electronic Hilbert space to a collection of SEOs, 33,34
where only the jth oscillator is on its first excited state, while others are on the ground state. The electronic operators k j ≡ Mkj are mapped to â † k âj such that the commutator relation â j , â † k = δ kj is satisfied, where â † k and âj are the creation and annihilation operators of quantum oscillators. The quantum-mechanically exact mapping is given by
Here, we used that â †
, and [qj, pk ] = i hδ jk , where qj and pj are the coordinate and momentum operators of quantum oscillators. Then, taking a classical limit, we
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp obtain j j = 1 2 h (q 2 j + p 2 j )γ and k j = 1 2 h (q kip k )(q j + ip j ), where γ is the zero point energy (ZPE) parameter (when γ = 12 full quantum ZPE is included). Inserting the above relation into Eq. (1), we have the MMST mapping Hamiltonian,
Alternatively, it can be equivalently expressed in terms of the symmetrized version, 37,70 H(R, P, q, p) =
where
is the average potential energy. The equations of motion for the semiclassical mapping dynamics are the corresponding Hamilton's equations. The integrator for the mapping Hamiltonian is described in Ref. 71. It is noted that other integrators are also available. [72][73][74][75] In nonadiabatic dynamics, the electronic reduced density matrix (RDM) is a key property, which is defined as σ(t) = TrN[ρ(t)] = ∑ F j,k σ jk (t) jk, where ρ(t) is the overall density matrix and the diagonal σjj(t) and off-diagonal σ jk (t) (j ≠ k) elements correspond to the electronic population and the electronic coherence, respectively. TrN[⋅] and Tre[⋅] stand for the trace over the nuclear and electronic Hilbert spaces, respectively. In a typical nonadiabatic process, the initial state is assumed to be ρ(0) = ρN(0) ⊗ σ(0) = ∑ F j,k ρN(0)σ jk (0) jk, where the initial nuclear density is given by ρN(0) = Tre[ρ(0)] and the initial electronic state is a population, for example, σ(0) = mm. The RDM at time t is given by 52
where the quantum time correlation function has the form
The LSC approximation of the time correlation function C Â B(t) is given by 46,50,[76][77][78][79][80][81][82][83][84][85]
LSC mapping is based on the Wigner transform of the quantum operator defined as
There exist two mapping strategies, 52 including mapping No. 1,
W (q, p) =
and mapping No. 2,
W (q, p) =
where a phase-space function of Gaussian form is defined as
Here, G(q, p) can be used for initial sampling of the mapping variables (q 0 , p 0 ) and choosing mapping No. 1 or mapping No. 2 [superscript (1)(2)] leads to LSC1 and LSC2 approximations (equivalent to PBME 51 and LSC-IVR, 45 respectively) to the electronic RDM σ jk (t) in Eq. ( 6), 52
σmn(0) dR 0 dP 0 dq 0 dp 0
The initial sampling of the nuclear DOF is over the Wigner function
which is also used in the rest of mapping dynamics. For a better population dynamics, one casts the population in terms of the identity operator, 1, and the traceless operator, Qj = F Mjj -∑ F k Mkk , using the resolution of identity (RI) ∑ F j j j = 1 such that the population of state j is given by j j = 1 F 1 + Qj . 42,86 The population and coherence become
Using the RI tricks as in Eq. ( 14), the traceless operator can be written in terms of mapping No. 1 or No. 2 as follows:
W (q, p).
(
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp Using G(q, p) in observables Qj(t) or Mkj (t) for the initial sampling of the mapping variables, the three ways for treating Qm(0) and the identity operator 1 lead to the following methods: RI-LSC1 treats the identity operator as 1 and uses mapping No. 1 for Qm(0); RI-LSC2 treats the identity operator as 1 and uses mapping No. 2 for Qm(0); and RI-LSC3 treats the identity operator as G(q, p) given in Eq. ( 11) and uses mapping No. 2 for Qm(0). 14,71 RI-LSC1 ∼3 are also known as mLSC/ϕ 1 ϕ 1 , mLSC/ϕ 1 ϕ 2 , and mLSC/ϕ 2 ϕ 2 , respectively. 52,87-90
The SQC dynamics proposed by Cotton and Miller provides a windowing approach for formulating the initial sampling and the RDM. 37,38,53,70,91,92 The Cartesian mapping variables can be expressed equivalently in terms of action-angle (a-a) variables, (n, u), where n = (n 1 , . . . , nF) and u = (u 1 , . . . , uF) such that nj and uj are the action and angle variables of the jth harmonic oscillator: n j = 1 2 h q 2 j + p 2 jγ, u j =arctan (p j q j ). Taking the Wigner transform of electronic operators kj in terms of a-a variables,
In the SQC method, the delta functions in Eq. ( 16) are replaced with finite-width window functions, and the two popular choices are square windows 37 or triangle windows. 38 The SQC-square uses window function wa(nj) = h(γ -nja), where h(x) is the Heaviside step function and the reported optimal ZPE parameter value γ = 0.366 for the square window 37 is used,
The SQC-triangle uses the window function with the reported optimal value γ = 13, 38
[Mjj] (SQC-tri)
where wa(n j ) = (2γn j ) 2-F when (-γ < na < 1γ) and vanishes otherwise; w 0 (nj, n l ) = 1 when (n l < 2 -2γnj) and vanishes otherwise. It is typically recommended to use the triangle window within SQC. 38 The initial values of the action mapping variables n are sampled randomly within the initial window, while the angle mapping variables u are sampled uniformly in [0, 2π). During propagation, the a-a variables can be transformed back to the Cartesian mapping variable (q, p) via relations q j = 2(n j + γ) h cos (u j ) and p j = -2(n j + γ) h sin (u j ) such that the equations of motion are the same as LSC mapping dynamics. To ensure that the total population is one for all time, the RDM elements should be renormalized via σ jk (t) = σ raw jk (t)∑ F l σ raw ll (t), where σ raw jk (t) is the raw electronic reduced density matrix element at time t. The diabatic representation is chosen over the adiabatic representation for the mapping dynamics calculations in this work. 91,93
The eCMM method 39,59,60,75 can be formulated in terms of the mapping Hamiltonian in Eq. ( 4), but the initial conditions of the mapping variables are uniformly sampled over the full-sphere,
If the initial electronic RDM is σ(0) = mm, then the population of the nth state at time t is given by
where γ = 1-γ 1+Fγ and
The SPM method has Q, P, and W representations, 41 which are equivalent to eCMM with γ = 0, 1, ( √
F as in the Wrepresentation SPM, γ = γ and Q( γ, γ) = 1 such that the initial and final electronic population expressions become identical. It is noted that the W scheme of SPM is recommended in the literature, 41 and we keep the P, Q, W labels to mark the correspondence of the special γ values in eCMM.
The Q, P, W schemes of SPM could also use focused initial sampling for the mapping variables, which are denoted as Q-foc, P-foc, and W-foc, respectively. In the case of  = mm, B = k j, the focused distribution is given by
Then, AW(q, p) = 1, BW(q, p) = 1 2 h (q kip k )(q j + ip j )γδ kj , and the time correlation function is given by
where Z foc = ∫ dqdpρ foc (q, p). The MMST-foc scheme in our implementation corresponds to the case when γ = γ = 12 and Q( γ, γ) = 1. 40,41 The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
Ehrenfest dynamics treats the nuclear and electronic DOF in a mean-field manner. 94,95 The time-dependent wave function can be expanded as a linear combination of electronic states, i.e., Ψ(t) = ∑ j cj(t)ψ j with expansion coefficients {cj(t)j = 1, . . . , F}. The equation of motion for the electronic RDM, σ jk (t) = c j (t)c * k (t), is given by the quantum Liouville equation @ @t σ(t) =i h V(Rt), σ(t). The electronic propagation can be expressed as σ(t
Here, the time evolution operator Û(t + t, t) in the diabatic basis can be expressed in terms of adiabatic energies and transformation matrix
Moreover, the nuclear equations of motion are given by Ṙ 71 Although the equations of motion can also be formed the same way as in the above-mentioned mapping dynamics with γ = 0, Ehrenfest dynamics has no initial sampling for the electronic DOF since the initial electronic RDM elements σ jk (0) are known.
Within the FSSH method, the electronic expansion coefficient dynamics are dictated by the equation of motion derived from the time-dependent Schrödinger equation:
where the nonadiabatic coupling (NAC) vector is defined as d ab = ψ a ∇Rψ b , Ṙ is the nuclear velocity, and ψa @ @t ψ b = d ab ⋅ Ṙ. In surface hopping dynamics, only one electronic state can be the active state at a time and the nuclear dynamics are determined by the active state, ψ a , where the effective forces on the nuclear DOF from the active state's PES are given by
The active state can be determined by the hopping process from the current active state, ψ a , to another state, ψ b , with the transition probability,
where t serves as nuclear integration step length and electronic density matrix elements σ ba = c b c * a . To determine whether a hopping from state ψ a to ψ b occurs, a uniformly distribution random number ζ ∼ U(0, 1) is sampled and the transition is accepted if
When the adiabatic basis is used such as in Tully's models, the off-diagonal potential energy elements V ab vanish and thus the first term in g a→b vanishes, whereas in the diabatic basis that was employed in the other models here, the NAC vanishes, and thus, the second term in g a→b vanishes [Eq. ( 26)]. After a hopping, the nuclear momenta should be adjusted such that the total energy is conserved, which means that the kinetic energy should be added
, so the momenta are scaled as
However, if a hop is to an energetically higher state but the system has no sufficient kinetic energy, the transition will be aborted, also known as the frustrated hop.
The numerically exact methods that will be employed in this study are the tensor-train (TT) based wavepacket dynamics approaches. In this subsection, we first briefly explain the TT format that facilitates the storage of and operations on high-dimensional wavepacket and then introduce the TT-SOKSL and TT-TFD methods used for zero temperature and finite temperature real-time wavepacket propagations, respectively.
Consider a general time-dependent wavepacket with d degrees of freedom, Ψ(t), which can be expanded as follows:
where {j k } with k = 1, . . . , d is the basis set of size n k for the kth physical degree of freedom and X(t) is a d-dimensional array that stores all n 1 × ⋅ ⋅ ⋅ × n d time-dependent expansion coefficients. The size of X (time variable omitted henceforth) grows exponentially with d, causing the operation and storage of Ψ(t) with d > 10 intractable in the standard grid-based representation. The tensortrain methodology provides a numerically accurate way of representing the high-dimensional X with a train of three-dimensional arrays whose size scales linearly with d and therefore bypasses the so-called curse of dimensionality. The TT format of the complex
, such that any element of X, namely, X(j 1 , . . . , j d ), is approximated by the following tensor product:
which is equivalently expressed in the matrix product state form,
where the r k-1 × r k matrix X k (j k ) is the j k th slice of the r k-1 × n k × r k core tensor X k . According to Eq. ( 29), every element of X is obtained only with the three-dimensional tensors X 1 , . . . ,
then only dnr 2 elements are required to be stored explicitly. The set {r k } thus determines the storage with given n and d and is hence referred as TT-rank in analogy with the matrix rank. The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp 2. Tensor-train split-operator KSL (TT-SOKSL) method for zero-temperature real-time propagation
The TT-SOKSL method integrates the time-dependent Schrödinger equation,
where Ĥ = T + V is the total Hamiltonian and Ψ(t) is the TT wavepacket as introduced in Eqs. ( 28)-( 30) with {j k } in Eq. ( 28) being the standard position grid basis. Similar to the SOFT method, TT-SOKSL integrates Eq. ( 31) to the second-order accuracy with the Strang splitting method, which approximately evolves the wavepacket from Ψ(t k ) to Ψ(t k+1 ) for time step t = t k+1t k by solving the following sequence of differential equations:
(1) Integrate equation:
(2) Integrate equation:
(3) Integrate equation:
from t = t k to t = t k+1 , with initial condition Ψ 3 (t k ) = Ψ 2 (t k+1 ) to obtain Ψ 3 (t k+1 ), which gives the approximate solution Ψ(t k+1 ) = Ψ 3 (t k+1 ) + O(t 3 ).
TT-SOKSL integrates Eqs. ( 32)-( 34) sequentially with the TT-KSL integrator in combination with the SOFT approach. In Eq. ( 32), the wavepacket Ψ 1 (t) is in the position space representation, under which the position operator -i V2 h is a diagonal operator. Therefore, the integration of Eq. ( 32) is carried out with the TT-KSL integrator with an efficient TT diagonal matrix-vector multiplication or equivalently an elementwise TT vector-vector multiplication instead of the usual TT matrix-vector multiplication. The resulting Ψ 1 (t k+1 ) is then Fourier transformed to the momentum space so that the integration of Eq. ( 33) with kinetic operator -i T h is also carried out with an TT vector-vector multiplication. The resulting Ψ 2 (t k+1 ) is then inverse Fourier transformed back to the position space, and the integration of Eq. ( 34) is carried out in a manner similar to that of Eq. (32). To sum up, a complete TT-SOKSL integration step is performed with the following equation:
where KSL denotes one substep of rank-adaptive TT-KSL propagation.
The TT-TFD method solves the quantum Liouville equation,
by transforming the density operator ρ(t) from a matrix in the physical Hilbert space to a vector in the so-called "double space" and consecutively transforming Eq. ( 36) to a commutator-free differential equation that takes a form similar to Eq. ( 31), whose integration provides all dynamical information of the system at finite temperature. In this study, the density matrix is initiated at the initial thermal equilibrium distribution, and the Schrödinger-like propagation of its corresponding finite-temperature double space state vector is therefore seen as a finite-temperature wavepacket propagation.
To help writing the density matrix ρ(t) as a state vector, we define the double space unit vector I,
where i = ĩ indexes over the basis vectors of the system's Hilbert space H and ĩ is an exact copy of i, expanding a "fictional space" that is an exact copy of H. Therefore, I defined in the double space has dimensionality twice that of H. With I, we define the initial thermal wavepacket 0(β),
from which the thermal equilibrium density matrix ρ(0, β) = e -β Ĥ Tr {e -β Ĥ } ≡ e -β Ĥ Z(β) is fully recovered upon taking outer product with its dual and tracing out the fictional DOF,
Furthermore, we define the finite-temperature double-space wavepacket ψ(β, t) according to the initial value ψ(β, 0) = 0(β), and the equation of motion is given by dψ
where H = Ĥ ⊗ 1, with 1 being the fictional space identity operator. With this definition, it is readily verified that the time-evolved density matrix ρ(t) is recovered from ψ(β, t),
Therefore, ρ(t) as well as all corresponding dynamical observables are obtained by propagating ψ(β, t) according to Eq. (40). Furthermore, all physical dynamical observables are invariant upon redefining H in Eq. ( 40) as H = Ĥ ⊗ 1 -1 ⊗ H f , where H f is an arbitrary operator in the fictional space. This replacement provides mathematical flexibility necessary for a more compact equation of motion for model systems investigated in this study.
The preparation of the initial state ψ(β, 0) as in Eq. ( 38) requires imaginary time propagation, which could be computationally burdensome. However, when the Hamiltonian is harmonic, we can choose
This task is circumvented utilizing the thermal Bogoliubov transformation, which provides the matrix Ĝ that rotates the double-space vacuum state 0, 0 = 0 ⊗ 0 to 0(β),
where Ĝ = -i
θ j = arctanh(e -βω j 2 ), and âj, â † j (ãj, ã † j ) are the physical (fictional) bosonic creation and annihilation operators for the jth degree of freedom. Substituting Eqs. ( 43) and ( 44) into Eq. ( 40), we obtain
where
TT-TFD expresses 0, 0 as a rank-1 tensor train [Eq. (30)] and propagates Eq. ( 45) to obtain ψ θ (β, t). The propagation is carried out with the TT-KSL integrator. The expectation value of any physical observable  ∈ H is given by
where Āθ = e i Ĝ Āe -i Ĝ and Ā = Â ⊗ Ĩ.
The spin-boson model Hamiltonian for a two-level system with the donor state D and the acceptor state A can be expressed as follows:
Here, σx = DA + AD and σz = DD -AA are the Pauli matrices; Γ is the electronic coupling that is assumed to be constant in the Condon approximation; E = -2ε is the donor-to-acceptor reaction free energy; { R j , Pj, ω j } = { R j , Pj, ω j j = 1, . . . , N} are the mass-weighted coordinates, momenta, and frequencies associated with the N nuclear normal modes, respectively; and {cj} are the electronic-vibrational coupling coefficients. The mode frequencies and the electronic-vibrational coupling coefficients {ωj, cj} are determined by discretizing the spectral density defined as
The spectral density can be assumed to be in a closed-form expression with empirical parameters 6,97-100 or obtained from taking Fourier transform of the time correlation function of the energy gap fluctuations in all-atom molecular dynamics (MD). 9,101 In this work, we employ the Ohmic spectral density, which is given by
where ξ is the Kondo parameter and ωc is the cutoff frequency.
The displacement between the donor and acceptor equilibrium geometry along the jth mode is R eq j = 2c j ω 2 j ( j = 1, . . . , N). The reorganization energy between the two states is given by
In this work, we will examine nine spin-boson models, ranging from biased and unbiased energetic driving forces, small and large reorganization energies, as well as different electronic coupling strengths as summarized in Table I. Here, models No. 1-No. 4 and No. 9 are defined with the same electronic coupling Γ = 1, while and ξ are varied. Models No. 5-No. 8 are all biased cases with a downhill reaction free energy but different Γ and ξ parameters. It is noted that a larger ξ leads to a larger reorganization energy since Er ∼ 2hξω c . We categorize the spin-boson models in terms of adiabatic vs nonadiabatic regimes according to the ratio ωcΓ: the adiabatic regime (ωcΓ < 1) features large electronic coupling and slow nuclear motion, whereas the nonadiabatic regime (ωcΓ > 1) features small electronic coupling and fast nuclear motion. 102 As shown in Table S1, spin-boson model No. 1 is in the intermediate regime (ωcΓ = 1), models No. 2-6 are in the nonadiabatic regime, and models No. 7-9 are in the adiabatic regime. In the Marcus weak-coupling nonadiabatic limit, we further categorize the models into the normal region (-E < Er) and inverted region (-E > Er). Thus, models No. 2, No. 3, and No. 5 are in the inverted region, while models No. 4 and No. 6 are in the normal region. The nuclear tunneling effects are expected to be more significant in the inverted region than the normal region at the nonadiabatic limit and also expected when the barrier is thin at the adiabatic limit.
The initial electronic state is set to be in the donor state, whereas the initial nuclear state is in equilibrium with the donor-state PES for models No. 1-No. 4 and No. 9, in equilibrium with the unshifted bath Hamiltonian for models No. 5-No. 8,
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp
TABLE I. Spin-boson models' parameters, including energy offset , electronic coupling Γ, Kondo parameter ξ, and cutoff frequency ωc of Ohmic spectral density, maximum frequency ωmax, the number of nuclear modes N, and time step dt in propagation. Here, models No. 1 to No. 4 adopted from Ref. 103 and model No. 9 are given in the reduced unit, and models No. 5 to No. 8 adopted from Ref. 104 are given with energy unit in cm -1 . Model No. Γ Er ξ ω c ωmax N 1 1.0 1.0 0.198 0.1 1.0 5.0 60 2 1.0 1.0 0.397 0.1 2.0 10.0 60 3 1.0 1.0 1.488 0.1 7.5 36.0 60 4 0 1.0 0.992 0.2 2.5 12.0 60 9 1.0 1.0 0.080 0.08 0.5 5.0 60 5 50 20 39.73 20ωc 53 265 60 6 50 20 119.2 60ωc 53 265 60 7 50 100 79.46 40ωc 53 265 60 8 50 100 198.6 100ωc 53 265 60
In the nonadiabatic mapping methods, the initial nuclear positions and momenta, {R 0 , P 0 }, were sampled from the semiclassical Wigner function at finite temperature,
where β = 1kBT and Hi(R 0 , P 0 ) is the classical nuclear Hamiltonian of i, on which PES the nuclear initial conditions are sampled. The initial nuclear sampling at zero temperature corresponds to the limit as β → ∞,
Hi(R 0 , P 0 ).
We consider the linear vibronic coupling (LVC) models for describing electronic transitions through conical intersections in the fulvene molecule, the 2,6-bis(methylene) adamantyl (BMA) radical cation, and the 2-methylene-6isopropylidene adamantyl (MIA) radical cation, whose parameters were adopted from Ref. 105. The LVC Hamiltonian is given by
where the nuclear operators are
TABLE II. Linear vibronic coupling (LVC) model parameters, including the number of nuclear modes N, reorganization energy Er , reaction free energy E, initial electronic state, and initial nuclear state (energy unit in a.u.). Fulvene BMA MIA N 30 78 96 Er 0.0887 0.0297 0.0274 E 0.0989 -0.0004 -0.0250 Init. elec. state A D D Init. nucl. state VD VA VA
Here, { R j , Pj j = 1, . . . , N} are the mass-weighted coordinates and momenta, N is the number of nuclear DOFs, the two electronic states are D and A, {ωj} are the nuclear mode frequencies, {g j } are the electronic-vibrational couplings, and D/A are the minimum energies of PESs of the two electronic states. The important feature of the LVC model is the non-Condon diabatic coupling that depends on the current coordinates, and the linear diabatic coupling coefficients are {γ j }. The equilibrium position of the jth nuclear mode on the acceptor state A is R eq j = g j ω 2 j ( j = 1, . . . , N), and the reorganization energy is
The LVC model parameters of the three molecular systems and the initial electronic and nuclear states are summarized in Table II and the supplementary material with the mode-specific values of {ωj, g j , γ j }.
We consider a two-state 25-mode model Hamiltonian for the photoisomerization of the retinal chromophore in rhodopsin, which is a key process in vision. This retinal model adopted from Ref. 106 contains two primary modes, including the cis-trans reaction coordinate R1 = θ that describes the anharmonic C 11 === C 12
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp torsion and the linear coupling mode R2 = Rc that describes the polyene chain stretching, as well as 23 secondary harmonic bath modes { R j j = 3, . . . , 25}. The two electronic states 0, 1 correspond to S 0 and S 1 diabatic states that are linearly coupled. The retinal model Hamiltonian is given by
Here, I is the moment of inertia, W 0 and W 1 are PES slopes of the rotational degree of freedom for 0 and 1, respectively, E 1 serves as the vertical shift of 1, c is the frequency of the coupling mode, Rc, {ωjj = 3, . . . , 25} is the frequencies of the bath modes, κ is the coupling coefficient between Rc and the electronic DOF, cj are the coupling coefficients between the bath modes and the electronic DOF, and κ is the linear coupling coefficient in the off-diagonal diabatic coupling. The momenta of the primary modes are P1 = Pθ = -i @ @θ and P2 = Pc = -i @ @R c ( h = 1). The model parameters c = 1532 cm -1 and the following parameters are in eV: 1I = 4.84 × 10 -4 , E 1 = 2.48, W 0 = 3.6, W 1 = 1.09, κ = 0.1, and λ = 0.19. The positions and momenta operators will be dimensionless in this case. The bath mode frequencies and electronic-vibrational couplings {ωj, cjj = 3, . . . , 25} (in cm -1 ) are detailed in the supplementary material. For the photoisomerization reaction, the initial electronic state is the S 1 state and the initial nuclear condition is in equilibrium with the PES of the S 0 state. The initial state of the anharmonic primary mode θ corresponds to the normal distribution of ρ(θ) ∝ e -θ 2 (2σ 2 θ ) where σ θ = 0.152 282 752.
Tully's models refer to a set of two-state one-dimensional models described by the total Hamiltonian Ĥ = T + V and electronic Hamiltonian V = V00 00 + V11 11 + V01 01 + V10 10, which is originally for testing scattering dynamics using the FSSH method. 107 It can be employed to benchmark various models for their ability to describe nuclear quantum effects. 52
The potential and coupling of the SAC model are given by V00
where V 0 and V 1 represent PESs of 0 and 1, respectively, and V 01 and V 10 represent the diabatic couplings. Parameters A = 0.01, B = 1.6, C = 0.005, and D = 1.0. The initial state is on electronic state 0, and the nuclear initial position at x 0 = -9 and initial momentum p 0 = 5, 10, 20 is described by the Gaussian wavepacket,
which corresponds to the Wigner function for semiclassical initial sampling given by
with σx = 1 √ 2 and σp = h(2σx).
The potential and coupling of the SAC model are given by
where parameters A = 0.10, B = 0.28, C = 0.015, D = 0.06, and E 0 = 0.05. The initial state is the same as in the SAC model, except the initial momentum p 0 = 5, 15, 30.
The line cross parabola (LCP) model is a one-dimension twostate Hamiltonian 108 Ĥ = P 2 2 + V( R), where R and P are the nuclear position and momentum, respectively, and the potential energy V( R) is given by
Here, the electronic coupling is Γ = 0.5 and the parabola (ω = 1) and the linear potentials are
The initial electronic state is the population on the parabola, i.e., 11. The initial nuclear state is chosen to be a certain vibrational eigenstate of the parabola, such as n = 15 and n = 30. In semiclassical simulations, the initial nuclear phase-space points are sampled from the microcanonical distribution with the fixed energy of the corresponding eigenstate, as En = hωn + 1 2 for the nth eigenstate. With a random number ξ drawn from the uniform distribution in [0, 2π), the initial positions and momenta are sampled according to
ARTICLE pubs.aip.org/aip/jcp
In the numerically exact simulation using SOFT, the initial wave function corresponding to the nth eigenstate is
where α = ωh and Hn(R) is the nth-order Hermite polynomial.
The simulations of nonadiabatic dynamics using mixed quantum-classical, semiclassical, and quasiclassical approximate methods employed here include mean-field Ehrenfest, FSSH, five versions of LSC mapping dynamics (LSC1-2, RI-LSC1-3), SQC (with square and triangle window), eCMM, and SPM (P, Q, W, MMST schemes). Among these methods, LSC methods use Gaussian sampling of mapping variables [Eq. (11)] and eCMM/SPM use full-sphere sampling of mapping variables [Eq. ( 19)], while focused sampling of mapping variables can be implemented with SQC windowing [Eqs. ( 17) and ( 18)] and SPM approaches [Eq. ( 21)], and no sampling of mapping variables is needed in MQC methods of Ehrenfest and FSSH dynamics. It is noted that the Ehrenfest method can be viewed as a nonadiabatic mapping dynamics with ZPE parameter γ = 0 since they share the equations of motion. The initial nuclear conditions are obtained from sampling from the Wigner function as in Eq. ( 13) for all approximate methods here. The propagation of electronic DOFs, (q, p), was achieved by using the fourth-order explicit Runge-Kutta algorithm, and the nuclear DOF, (R, P), were propagated by the velocity Verlet integrator with time step t summarized in Table III. To improve the accuracy of integration, the electronic time step δt was chosen to be smaller than nuclear t such that 20δt = t. The nonadiabatic dynamics results of these mapping dynamics methods were obtained by averaging over 10 4 -10 6 trajectories, and the total averaged trajectories for different models are summarized in Table III. The convergence of the approximate results concerning the time step and number of trajectories has been checked. The ZPE parameter γ was chosen to be method specific: all five versions of LSC dynamics use γ = 12, SQC-square uses γ = 0.366, SQC-triangle uses γ = 13, SPM-W uses γ = √ F + 1 -1F, SPM-P uses γ = 1, SPM-Q uses γ = 0, MMST uses γ = 12, and eCMM uses γ = -0.2. An adiabatic FSSH algorithm is utilized to approximate dynamics of Tully's model. Nuclear propagation used the velocity-Verlet algorithm, mapping variables are propagated using the second-order Verlet algorithm, and nuclear variables share the same updating frequency as mapping variables, δt = t. Except for aforementioned simulation details, the adiabatic FSSH shares the same simulation parameters as its diabatic counterpart.
The computational costs for the approximate dynamical methods are summarized in Table S6, for example, the mapping dynamics of a spin-boson model averaged over 10 5 trajectories typically cost from 1.0 to 8.1 core-hours with Intel Xeon Gold 6132 CPU with 2.6 GHz base frequency, and other approximate dynamics nearly scales with the number of trajectories averaged. For the TT calculation, the computational cost heavily depends on the parameters chosen, including the number of basis, TT-rank, and time step, which in turn depends on the specific model. A typical range of computational cost for TT-TFD of spin-boson models is from 4.75 to 12 940 core-hours with dual Intel Xeon Gold 6132 CPU as shown in Table S7.
We note that the focus of the current work is the numerical performance of the approximate dynamical methods, rather than their theoretical rigor when they were proposed. For example, even if the Ehrenfest mean-field method can be rigorously derived, its accuracy cannot be automatically guaranteed. In addition, despite SQC windowing choice being ad hoc, its performance is decent in several cases tested here. The LSC series can be rigorously derived by linearizing the real-time path integral expression of TCF and is capable of capturing important but not all the nuclear quantum effects (NQEs). The accuracy of LSC mapping dynamics depends on using the resolution of identity tricks and the specific choices on the mapping strategies for the electronic operators. The Wigner initial nuclear sampling is naturally derived in LSC, but the original MQC methods were prescribed with classical sampling. To facilitate a direct comparison of methods, we utilized Wigner sampling in all approximate methods studied here. Still, it is important to point out
TABLE III. Simulation parameters of nonadiabatic semiclassical and quasiclassical dynamics of different models, including time step (t) of nuclear DOF and the number of trajectories averaged (N traj ) for MQC dynamics (Ehrenfest and FSSH), LSC mapping dynamics, and other mapping dynamics. Models t N traj (MQC) N traj (LSC) N traj (else) Spin-boson Nos. 1-4,9 0.001hΓ 10 4 10 6 10 5 Spin-boson Nos. 5-8 0.001 ps 10 4 10 6 10 5 Fulvene 1 × 10 -6 ps 10 4 10 6 10 5 BMA 5 × 10 -6 ps 10 4 10 6 10 5 MIA 5 × 10 -6 ps 10 4 10 6 10 5 SAC 1 × 10 -6 ps 10 4 10 6 10 5 DAC 1 × 10 -6 ps 10 4 10 6 10 5 LCP 0.001 a.u. 10 4 10 6 10 5 Retinal 1 × 10 -5 ps 10 4 10 6 10 5
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp that the accuracy of the MQC methods might be improved by using Wigner sampling that essentially captures some static NQE. In addition, in the eCMM method, any ZPE parameter γ ∈ (-1F, ∞) can be chosen, but the optimal value of γ may vary with the model.
For all tensor-train-based calculations, convergence tests were performed to ensure the population dynamics results converged to the quantum exact ones. The variables for the convergence test include the time step (t), TT-rank, and number of basis functions, as listed in Table IV. We performed the convergence tests by calculating population dynamics with these three parameters varied and concluded that the results are converged when the changes in population dynamics are no longer observed with added computational cost. We found that all three variables are highly system-dependent. The time step for the spin-boson models was set as t = 0.0036hΓ for models No. 1-4 and t = 10 a.u. for models No. 5-8. For the LVC models, a time step of t = 0.5 a.u. = 0.012 fs was used, and for the retinal model, the time step was set as t = 2.5 a.u. = 0.0301 fs. For the spin-boson models and the LVC models, the wavefunctions in the TT-KSL method and thermal wavepacket in the TT-TFD
TABLE IV. Tensor-train-based method and convergence parameters for different models. Model Method Temperature Time step t TT-rank Basis a No. of basis b Grid span c Spin-boson No. 1 TT-KSL 0 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 1 TT-TFD 0.2 Γk B 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 1 TT-TFD 2 ΓkB 0.05 hΓ 80 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 2 TT-KSL 0 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 2 TT-TFD 0.2 ΓkB 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 3 TT-KSL 0 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 3 TT-TFD 0.2 ΓkB 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 3 TT-TFD 2 ΓkB 0.012 5 hΓ 115 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 3 TT-TFD 10 ΓkB 0.006 25 hΓ 90 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 4 TT-KSL 0 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 4 TT-TFD 0.2 ΓkB 0.003 6 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 5 TT-KSL 0 10 a.u. 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 5 TT-TFD 300 K 10 a.u. 150 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 6 TT-KSL 0 10 a.u. 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 6 TT-TFD 300 K 10 a.u. 140 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 7 TT-KSL 0 10 a.u. 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 7 TT-TFD 300 K 10 a.u. 150 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 8 TT-KSL 0 10 a.u. 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 8 TT-TFD 300 K 10 a.u. 160 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 9 TT-KSL 0 0.1 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 9 TT-TFD 0.2 ΓkB 0.1 hΓ 20 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 9 TT-TFD 2 ΓkB 0.1 hΓ 40 osc 25 ⋅ ⋅ ⋅ Spin-boson No. 9 TT-TFD 10 ΓkB 0.1 hΓ 120 osc 25 ⋅ ⋅ ⋅ Fulvene TT-KSL 0 0.5 a.u. 40 osc 40 ⋅ ⋅ ⋅ Fulvene TT-TFD 300 K 0.5 a.u. 40 osc 40 ⋅ ⋅ ⋅ BMA TT-KSL 0 0.5 a.u. 40 osc 45 ⋅ ⋅ ⋅ BMA TT-TFD 300 K 0.5 a.u. 40 osc 45 ⋅ ⋅ ⋅ MIA TT-KSL 0 0.5 a.u. 45 osc 45 ⋅ ⋅ ⋅ MIA TT-TFD 300 K 0.5 a.u. 45 osc 45 ⋅ ⋅ ⋅ SAC P 0 = 5 SOFT 0 4 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) SAC P 0 = 10 SOFT 0 2 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) SAC P 0 = 20 SOFT 0 1 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) DAC P 0 = 5 SOFT 0 2 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) DAC P 0 = 15 SOFT 0 2 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) DAC P 0 = 30 SOFT 0 1 a.u. ⋅ ⋅ ⋅ pos 1024 [-15, 15) LCP SOFT 0 0.001 a.u. ⋅ ⋅ ⋅ pos 2048 [-30, 30) Retinal TT-SOKSL 0 2.5 a.u. 70 pos 256/32 (θR) d [-π, π)/[-5, 5)
a "osc" indicates the occupation number basis of a harmonic oscillator, and "pos" indicates the position grid basis. b Number of basis refers to the maximum occupation number when used with "osc" and the number of grid points when used with "pos." c The range of spatial basis in the atomic unit for SOFT simulation.
d For the retinal model, 256 grid points in the position basis are used for the rotational mode and 32 grid points are used for the rest harmonic modes.
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp method [Eq. ( 45)] are expanded by the harmonic oscillator eigenbasis. Due to the high proximity of the spin-boson model to the harmonic oscillator, ten basis functions are usually sufficient for the convergence of population dynamics. However, for the sake of consistency, we performed simulations on spin-boson models all with a basis size of 25. For the convergence of TT-rank on the spin-boson models, TT-TFD poses a slightly higher requirement than TT-KSL due to the added entanglement between the system and fictional components of the thermal wavepacket. The LVC models have more complexity than the spin-boson models by including electronic coupling terms that are linear to vibrational modes, among which the most demanding TT-TFD for MIA model required a TT-rank of 45 and 45 harmonic eigenbasis. Meanwhile, for the more complicated retinal model, a dense position grid was used to expand the wavefunction due to the highly anharmonic nature of the retinal primary reactive mode. The primary mode requires at least 200 spatial grid points to converge, and the TT-rank would converge at around 70. The TT-TFD results of spin-boson models No. 1-4 at β = 5 have been confirmed to agree with the numerically exact QuAPI results, 69,109 and the TT-KSL results of spin-boson models No. 5-8 at zero temperature have been confirmed to agree with numerically exact ML-MCTDH. 104 In the SOFT method for Tully's models, grid points of 256 in the position basis are usually sufficient and we used 1024 position grid points to ensure a great convergence. The time step in SOFT propagation varies for different kinetic energies, and higher initial momentum usually requires smaller t. For the SAC model, in the P 0 = 5 case, t = 4 a.u. is sufficient; in the P 0 = 15 case, t = 2 a.u. is proper; and in the P 0 = 30 case, t = 1 a.u. is required. For the DAC model, in the P 0 = 5 case, t = 4 a.u. is sufficient; in the P 0 = 15 case, t = 2 a.u. is proper; and in the P 0 = 30 case, t = 1 a.u. is required. Grid points in the position basis are placed in an equally spaced manner, and the first point is located at the left limit of the grid span. For the LCP model, time step t = 0.001 a.u., and 2048 grid points are utilized for convergence.
The spin-boson model is one of the most widely used models for understanding nonadiabatic dynamics in condensed-phase systems. A faithful prediction of its electronic dynamics is key for testing quantum dynamical methods that could be potentially useful for large and complex condensed phases. The spin-boson model features harmonic PES and Condon (constant) electronic coupling, which are decent approximations for condensed phase systems. Here, we consider eight spin-boson models that correspond to a range of thermodynamic energetic driving force, reorganization energy, and electronic coupling, and the key model parameters are summarized in Table I.
The population dynamics of spin-boson model No. 1 in terms of the population difference between the donor state and the acceptor state, i.e., σz(t) = σDD(t) -σAA(t), are reported in Fig. 1. Spinboson model No. 1 has a biased energy offset ( = 1.0) and a small reorganization energy (Er = 0.198) compared with electronic coupling Γ = 1.0), which is expected to have very oscillatory population dynamics. It is evident that the population oscillations are more significant at kBT = 0 (upper panels of Fig. 1) than at kBT = 0.2
No. 1 at zero temperature (upper panels) and at k B T = 2 (lower panels) obtained with various approximate mixed quantum-classical, semiclassical, and quasiclassical dynamics methods (solid lines) in comparison with numerically exact dynamics at zero temperature by TT-KSL and at finite-temperature by TT-TFD (dashed line).
(lower panels of Fig. 1), which can be attributed to the thermal damping due to the broader nuclear distribution at a higher temperature. First, MQC methods, including Ehrenfest and FSSH, seem able to capture only the general oscillatory trend but not the amplitude and phase quantitatively, compared with the numerically exact TT-SOKSL dynamics for zero temperature and TT-TFD dynamics
The Journal of Chemical Physics
ARTICLE pubs.aip.org/aip/jcp Here, LSC1 and LSC2 can capture the oscillation phase, but underestimate and overestimate the population oscillation, respectively. In contrast, the population oscillation is well reproduced when the RI trick for the identity operator is employed in RI-LSC1-3 methods, especially since the exact result can be accurately reproduced by RI-LSC2 and RI-LSC3 approaches. Third, the SQC dynamics give an excellent agreement with the exact population dynamics, where the SQC with triangle window function behaves better than the square window function by having a smaller oscillation amplitude [Figs. 1(c1) and 1(c2)].
Fourth, the eCMM/SPM approaches are seen to capture the oscillation phase of population accurately, and the oscillation amplitude shows an anticorrelation with the ZPE parameter: the amplitude grows to better match the exact result when decreasing γ values from 1 (P scheme), 0.5, 0.366 (W scheme), 0 (Q scheme), and -0.2 [Figs. 1(d1) and 1(d2)]. We are not expecting that different ZPE parameter choices, including P, Q, W schemes of SPM, would lead to the similar accuracy; instead, the best eCMM/SPM method for the tested spin-boson models is when a negative ZPE parameter is chosen γ = -0.2. In short, among all the approximate methods taken into consideration, we found that RI-LSC2, RI-LSC3, SPM-Q/eCMM with γ = 0, and eCMM with γ = 0-0. No. 3 displays a more damped oscillation in the population difference between donor and acceptor states than spin-boson model No. 1 does (Fig. 1) as expected since the reorganization energy depicts the strength of electronic-nuclear coupling. The zero-temperature (kBT = 0) and the finite-temperature (kBT = 0.2) dynamics obtained with the same method are very similar, which can be traced back to the fact that the Wigner initial sampling widths at the finite temperature [Eq. (53)] are close to those at zero temperature [Eq. ( 54)] for those high-frequency modes due to high ωc and ωmax values. The exact result of spin-boson model No. dynamics. Ehrenfest dynamics averages over quantum states, which disrupts the accurate representation of quantum superpositions and thus their coherence, while the original FSSH lacks explicit decoherence treatment. LSC methods provide better coherence dynamics than MQC methods, accurately reproducing the oscillation phase of both the real and imaginary parts of the coherence. Among these, LSC methods with RI treatment appear to perform particularly well. The SQC methods with triangle and square windows seem to behave identically and exhibit overdamped coherence, attributed to the inherent broadening effects of these window functions. The eCMM methods with γ = -0.2 and γ = 0 (SPM-Q) show the best agreement with the exact coherence dynamics among all approximate methods tested.
We show other spin-boson models' comparison between approximate and exact dynamical methods in the supplementary material, where the population dynamics of spin-boson models No. -S12. Generally, for all spin-boson models tested here, the traditional MQC methods, including Ehrenfest and FSSH, give the worst prediction to the population dynamics and even cannot capture the phase of the population oscillation. LSC1 and LSC2 tend to underestimate and overestimate the population difference, respectively, while LSC with RI trick gives a much better prediction, in which RI-LSC1 often has a smaller oscillation
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp amplitude than RI-LSC2 and RI-LSC3. SQC-triangle gives similar results as SQC-square, except for very early time square window might produce incorrect flat population dynamics due to the finite window width. SPM-Q and eCMM with negative ZPE parameters tend to give better agreement with the exact result for oscillatory population dynamics. The root-mean-square errors of the approximate methods from the exact results for all spin-boson models considered here are summarized in Table S5. For the coherence dynamics, similar performance can be observed that the LSC with RI treatment, eCMM with γ = 0 (SPM-Q), and γ = -0.2 are generally accurate, and the SQC methods would generate accurate coherence dynamics when they can capture the population dynamics. To sum up, RI-LSC2, RI-LSC3, SQC-triangle, eCMM with γ = 0 (SPM-Q), and γ = -0.2 are seen to provide excellent agreement with the exact population dynamics for the spin-boson models. Some model-specific comments are provided below. Spinboson model No. 2 shown in Fig. S2 has a biased energy offset and intermediate reorganization energy; thus, the population dynamics are expected to be in between models No. 1 and No. 3, where we observe intermediately oscillatory population decay and more damped oscillation at a higher temperature than at zero temperature, and the optimal dynamical method choices for model No. 2 are RI-LSC2,3 and SPM-Q/eCMM (γ = 0). Spin-boson model No. 4 shown in Fig. S4 has an unbiased energy offset and intermediate reorganization energy; thus, the population difference is expected to vanish in the long-time limit after a few damped oscillation periods, which is observed, and the best method for this model is eCMM (γ = -0.2) with the largest oscillation amplitude.
Spin-boson models No. 5-No. 8 have the same biased energy offset ( = 50 cm -1 ) but different electronic couplings and reorganization energies. Spin-boson models No. 5 and No. 6 shown in Figs. S5 and S6, respectively, have a small electronic coupling (Γ = 20 cm -1 ), model No. 5 has a small reorganization energy, and model No. 6 has a large reorganization energy. Both of the models display an overdamped population decay, especially a monotonic decay as seen in model No. 6. The population dynamics at zero temperature can be better captured by semiclassical dynamics methods than the population dynamics at finite temperature in both models No. 5 and No. 6, and the growing deviation between approximate dynamics and TT-TFD dynamics at 300 K at a long time seems to require double checking with another numerical exact method. In addition, the SQC-square method does not show population decay at the early time compared with the exact results, and in contrast, SQCtriangle does not suffer from this deviation. This phenomenon can be attributed to the fact that different window shapes and sizes affect the determination of the population of a particular state: the square windows with γ = 0.366 have a gap between two electronic states, whereas the triangle windows with γ = 13 contact at the vertices, and thus, SQC-square is seen to have a delay in population transfer, whereas SQC-triangle does not. Spin-boson models No. The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp limit with a large number of modes, the semiclassical and quasiclassical methods seem unable to capture this feature when the number of modes is limited. In Fig. 4, we show a much more friendly spin-boson model No. 9, where almost all approximate methods agree with the exact results well at temperatures kBT = 0 and kBT = 10, except for a slight out-of-phase deviation seen in MQC methods. An important feature of this model is the small reorganization energy (Er = 0.080) compared with the electronic coupling (Γ = 1.0) and the energy offset ( = 1.0). This model is also easy to get numerically exact results with the tensor-train approaches, and the time step can be relatively large (t = 0.1hΓ) compared with other spin-boson models. Similarly, in Fig. S9, it is evident that all approximate methods demonstrate excellent agreement with the exact results at kBT = 0.2 and kBT = 2. Model No. 9 is expected to alleviate the stress experienced by our readers resulting from the complexities observed in the more challenging spin-boson models No. 1-8.
The non-Condon effect arises when the electronic coupling is no longer a constant but dependent on nuclear configuration, and such an effect would be important for understanding electronic transition through conical intersections and long-range electron transfer. To this end, we report the benchmark test for fulvene, BMA, and MIA molecules that are described by LVC model Hamiltonians with electronic couplings linearly dependent on the coordinates. These gas-phase molecular systems would be particularly more challenging to simulate than condensed-phase systems, for the reason that the molecular vibrational modes, the electronic-nuclear couplings, and the off-diagonal electronic couplings may be more structured than the condensed-phase case. These LVC models have been benchmarked for the LSC and SQC dynamics at zero temperature, 89 but not at finite temperature, and also were not tested with SPM/eCMM methods before. FIG. 5. Electronic transition through conical intersections in fulvene (left), BMA (middle), and MIA (right) molecules at zero temperature (upper) and T = 300 K (lower) obtained with various approximate and numerically exact methods. The donor state population of linear vibronic coupling (LVC) models of the three molecules is plotted as a function of time with the initial state described in Table II. The numerically exact dynamics at zero temperature and finite temperature are computed with TT-KSL and TT-TFD, respectively (dashed), which is compared with mixed quantum-classical, semiclassical, and quasiclassical dynamics methods. In Fig. 5, we show the combined results for fulvene, BMA, and MIA at zero temperature as well as at the finite temperature of 300 K. First, the population decay in fulvene exhibits a similar signature stepwise pattern in both zero temperature and finite temperature cases, and most of the donor state population is transferred to the acceptor state in ultrafast time scale (0.03 ps) (left panels of Fig. 5). MQC methods of Ehrenfest and FSSH, as well as LSC2, SPM/eCMM with γ = 1 (P scheme), and γ = 0.5, tend to underestimate the population decay, whereas RI-LSC3, SPM/eCMM with γ = 0 (Q scheme), and γ = -0.2 tend to overestimate the population decay, leading to an unphysical negative population in a long time. Excellent agreement with the exact dynamics can be found with RI-LSC1, RI-LSC2, LSC1, SQC-square, SQC-triangle, and SPM/ eCMM with γ = 0.366 (W scheme), and the best agreement is seen with RI-LSC1 for fulvene. Second, in the BMA case, the best agreement with exact dynamics is found to be RI-LSC1 again, the excellent agreement can be produced with LSC1 and SPM/eCMM with γ = 0.366 (W scheme), and the rest methods display significant error. Third, in the MIA case, we observe a faster population decay at 300 K than at zero temperature and MIA with the largest number of modes shows a smoother decay profile in contrast to the abovementioned two cases. The best agreement with exact dynamics in MIA is achieved with SQC-triangle dynamics, while other methods tend to have large errors. To sum up, the electronic transition through conical intersections in the three gas-phase molecules can be captured accurately with SQC-triangle for MIA and RI-LSC1 for fulvene and BMA.
The 25-mode photoisomerization model of retinal from rhodopsin visual system consists two primary modes, including the anharmonic cis-trans torsion coordinate θ and the linear coupling harmonic mode Rc, and these primary modes are coupled to 23 bath modes. The electronic coupling is linearly dependent on Rc, which accounts for the non-Condon effect. With anharmonic potential and non-Condon coupling, the retinal photoisomerization model emerges as one of the most challenging types of model to test. Besides the above-mentioned semiclassical dynamics methods, we added a few dynamics with focused initial sampling for the mapping variables to this model for better nuclear dynamics. The retinal molecule is initially photoexcited to the excited state (S 1 ) and in the cis isomer; then, the electronic population gets transferred to the ground state (S 0 ) and the nuclear structure evolves to gain trans isomer fraction in the picosecond timescale. Figure 6 shows both electronic population dynamics (upper panels) and nuclear isomerization dynamics (lower panels) of photoexcited retinal at zero temperature where the numerically exact result was obtained with TT-SOKSL. Finite-temperature calculation of this model is difficult since the Wigner transformed distribution of the rotational degree of freedom is too wide even at 4 K to converge tensor networks numerically. Electronic population dynamics is reported as the ground state (S 0 ) population as a function of time, and the exact quantum mechanical dynamics shows features including that the initial rise in population starts at about 0.1 ps after a few small-amplitude Rabi oscillations, the first peak around ground state population of 0.5 is reached at 0.3 ps, then a slight decay is seen during t = 0.25-0.4 ps, and subsequently the ground state population is slowly accumulated FIG. 6. Photoisomerization dynamics of 25-mode retinal at zero temperature, including the ground electronic state (S 0 ) population dynamics (upper panels) and the trans isomer fraction dynamics (lower panels) after initially photoexcited to the S 1 excited state in the cis isomer, as obtained with numerically exact TT-SOKSL and various mixed quantum-classical, semiclassical, and quasiclassical dynamics methods.
to about 0.65 in 1 ps. The signature first peak at around 0.25-0.3 ps and the subsequent growth in ground-state population are seen in W-foc, Q-foc, MMST-foc, and Ehrenfest methods, of which W-foc and MMST-foc produce the most accurate electronic dynamics
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp compared with the exact result. SQC and LSC methods [Figs. 6(a1) and 6(b1)] fail to capture the first peak and underestimate the population growth, of which LSC1 and RI-LSC1 generate relatively better predictions since they tend to foster population transfer as also seen in the spin-boson models. It is evident that SPM/eCMM methods with sphere sampling of mapping variables systematically underestimate the population growth when compared with the corresponding SPM with focused sampling of mapping variables [Figs. 6(c1) and 6(d1)]
, and this can be traced back to the fact that the focused sampling would give more concentrated initial nuclear forces from the PES of the initially populated electronic state, which would drive the system away from the reactant region. The FSSH, LSC2, RI-LSC2/3, P-foc, and SPM/eCMM (γ = -0.2, 0) barely capture the qualitative trend of the electronic dynamics, where the worst method is eCMM
The isomerization reaction is characterized by the cis-trans dihedral angular nuclear dynamics as displayed in Figs. 6(a2)-6(d1). The exact result by TT-SOKSL shows that the first peak for the trans isomer fraction of 0.55 is reached at t = 0.2 ps, the second peak of about 0.45 trans fraction is reached at t = 0.3 ps, and then, the trans fraction decays gradually to about 0.3 at t = 1 ps. The best nuclear dynamics are reproduced by W-foc and MMST-foc methods as shown in Fig. 6(c2), where the two-peak signature is seen with a small time delay, which is much more accurate in nuclear dynamics compared with their sphere-sampling counterparts in Fig. 6(d2). The superb performance in generating nuclear dynamics with the W-foc and MMST-foc is because the initial nuclear forces are guaranteed to be determined by solely the PES of the initially populated electronic state when focused sampling of mapping variables is utilized. In contrast, when full sphere sampling of mapping variables is used, such as in SPM/eCMM methods, the initial nuclear forces result from averaging all electronic PESs uniformly, which include the unpopulated electronic state, thus making the nuclear dynamics unphysical, at least for the initial ultrafast timescale. Ehrenfest and Q-foc overestimate the peak height with a 0.5 ps delay, but the main trans fraction profile is reproduced. Finally, FSSH, LSC, and SPM/eCMM fail in producing the trans fraction profile, and the reason for LSC and SPM/eCMM is that uniform Gaussian sampling or full sphere sampling would lead to unphysical initial nuclear forces, thus hindering the cis to trans molecular structural change. Considering both electronic and nuclear dynamics, the best approximate methods for retinal with anharmonic potential and non-Condon coupling are W-foc and MMST-foc using focused sampling of mapping variables.
Tully's classic models were initially proposed for testing the FSSH method, 107 and here, we test single avoided crossing (SAC) and double avoided crossing (DAC) models, which mimic scattering problems in one dimension. It is noted that both SAC and DAC models have anharmonic potential and non-Condon coupling, but in a low-dimensional nuclear space, which is different from any condensed-phase environment. The DAC model is shown in Fig. 7, where we observe two different nonadiabatic population transfer dynamics with small and large initial momenta. When the initial momentum is small (p 0 = 5), a temporary population transfer to FIG. 7. Nonadiabatic dynamics σz(t) = σ 00 (t)σ 11 (t) of Tully's double avoid crossing (DAC) model at zero temperature with low initial momentum p 0 = 5 (upper panels) and high initial momentum p 0 = 30 (lower panels), where the exact result is obtained with the SOFT method and the approximate results are obtained with various mixed quantum-classical, semiclassical, and quasiclassical dynamics methods.
the other diabatic state occurs during 0.05-0.1 ps; then, the population gets transferred back to the original state. This behavior is well-captured by Ehrenfest, FSSH, RI-LSC2, LSC2, SQC, SPM-W, and eCMM (γ = 12), and relatively large errors are seen with LSC1 and SPM-P methods. When the initial momentum is large (p 0 = 30), a permanent population transfer is recognized and most
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp of the approximate methods can generate the correct dynamics, except for SPM-P, which gives relatively large errors.
The SAC model is shown in Fig. S19, where most of the approximate dynamical methods give correct semiquantitative predictions compared with exact quantum dynamics. The differences between various methods with small initial momentum (p 0 = 5) are more significant than those with large initial momentum (p 0 = 20). In the case of small initial momentum, accurate results are produced by RI-LSC1, LSC1, LSC2, SQC-square/triangle, SPM-W, and eCMM (γ = 12), and relatively large errors are observed in Ehrenfest, FSSH, RI-LSC2, SPM-P, SPM-Q, and eCMM (γ = -0.2). In the case of large initial momentum, there are no significant differences among different approximate methods, and they are all quite accurate but do not perfectly reproduce the exact result either. The SAC model with initial momentum p 0 = 10 and DAC model with initial momentum p 0 = 15 are shown in Figs. S20 and S21, respectively. It is noted that in the current implementation of FSSH, we used the wavefunction expansion coefficients for calculating population, while counting active states that could give better population have been reported in Ref. 110. To sum up, for these one-dimensional Tully's models, most approximate methods would work, and LSC2, SQC, and SPM-W typically provide accurate results. However, we notice that SPM-P is consistently bad in most cases we tested in this study, which might be because the large ZPE parameter may cause temporary flipped PES and thus unphysical forces when (q 2 j + p 2 j )(2 h) is too small compared with γ.
The LCP model is a challenging model for semiclassical or mixed quantum-classical methods. This one-dimensional model is far from any multidimensional condensed-phase system, where the semiclassical or mixed quantum-classical approaches tend to give reasonable predictions. Here, we used two vibrational eigenstates of the parabola, including a low-lying case (n = 15) and a high-lying case (n = 30) as the nuclear initial conditions.
In the low-lying case shown in Fig. S22, the numerically exact result shows a clear change of decay slope around 6 a.u. Although Ehrenfest and FSSH approaches exhibit less deviations from the exact results than the semiclassical mapping approaches, the Ehrenfest method does not capture the change of slope and FSSH gives an opposite slope trend. The LSC, SQC, and SPM methods show larger population deviations from the reference than the MQC methods. In the high-lying case shown in Fig. S23, the numerically exact results show a few oscillations. RI-LSC2 is seen to capture the first and half oscillations, better than the other methods. For the MQC approaches, FSSH gives better results than the Ehrenfest method. In the two tests, we observe that most of the MQC and semiclassical methods falter with this model, except for very early population decay. This failure may stem from inadequacies in the approximate methods for generating effective nuclear forces, which may not even exist in the multistate and extremely anharmonic PESs.
In this work, we systematically benchmarked various mixed quantum-classical, semiclassical, and quasiclassical dynamics methods in model systems against numerically exact tensortrain-based quantum dynamics methods, including TT-KSL for zero-temperature cases and TT-TFD for finite-temperature cases. Generally speaking, the choice of optimal approximate dynamical methods is a case-by-case problem, and for different natures of models, the answer varies dramatically. For the spin-boson model with harmonic potential and Condon coupling, we found that RI-LSC2/3, SQC-triangle, SPM-Q/eCMM (γ = 0), and eCMM (γ = -0.2) could give accurate population dynamics at zero and finite temperatures. For the LVC model of electronic transition through conical intersections parameterized by harmonic potential and non-Condon coupling, SQC-triangle worked the best for the MIA system with a larger number of modes and RI-LSC1 worked the best for fulvene and BMA systems with fewer modes. For the photoisomerization reaction in retinal that is described by anharmonic potential and non-Condon coupling, we found that W-foc and MMST-foc with focused sampling of mapping variables could produce accurate electronic and nuclear dynamics, which is attributed to the accurate initial nuclear forces from the focused sampling on the initial electronic state. For Tully's SAC and DAC scattering problem in one dimension, most approximate dynamical methods would work, such as LSC2, SQC-triangle/square, and SPM-W. It is challenging to capture the exact quantum dynamics in the LCP model using the approximate methods employed here. We note that SPM-P with γ = 1 consistently produces inaccurate results in most of the tested models and thus is not recommended to use.
There are several future directions based on the current work. The approximate method choice for benchmark could never be complete since many approximate methods emerged in recent years. For example, the partial linearized density matrix (PLDM) 111,112 and spin-PLDM 113,114 are based on treating the forward and backward propagation of electronic reduced density matrix separately and nuclear propagation with full linearization approximation, but it is beyond the current semiclassical and quasiclassical dynamics that share the equation of motion and thus beyond the scope of this work. In addition, quantum dynamics in a more sophisticated multi-state harmonic (MSH) model 14,71 and multi-state reaction coordinate (MRC) model 115 that are constructed from all-atom simulations could be evaluated with the tensor-train-based methods, which in turn pave the way for rational selection of approximate methods to all-atom nonadiabatic dynamics. The TCFs provided by the suitable approximate methods could be directly used for understanding transient electronic transition, such as photoinduced charge and energy transfer in the condensed phase. [116][117][118] There are also several ways to correct the initial sampling of mapping variables such that the initial nuclear forces are physical, including changing ZPE parameters statically 54 or dynamically. 113,119 When applying nonadiabatic dynamics, it is important to validate the approximate dynamical methods before applying them to more complicated systems where exact results are impossible to get. Benchmark studies have provided and will keep providing practical guidance in simulations of realistic nonadiabatic processes in complex condensed phases, such as organic photovoltaics, luminescence, photosynthesis, and photochemistry. 14,37,38,[40][41][42]45,51,52,59,[70][71][72]86,[89][90][91][92] In this case, a necessary and essential step is to test with model systems of a similar kind when the exact result can be known. The benchmark study would provide practical guidance in simulations of realistic nonadiabatic processes in complex condensed phases, such as organic photovoltaics, luminescence, photosynthesis, and photochemistry.
ARTICLE pubs.aip.org/aip/jcp SUPPLEMENTARY MATERIAL
See the supplementary material for the discretization of the spectral density, the LVC model parameters, the population and coherence dynamics as well as the root mean square error (RMSE) of spin-boson models No. 1-9 at different temperatures, ground state population dynamics of retinal with focused sampling methods, Tully's SAC and DAC models with different initial momenta, LCP model results, computational cost of approximate dynamics, as well as the previous works on relevant models. The numerically exact results of the tested models are available at https://github.com/xiangsunlab/benchmark-data.