The article information
 Yunan Yan
 严运安
 LowStorage RungeKutta Method for Simulating TimeDependent Quantum Dynamics
 含时量子动力学模拟的低存储龙格库塔方法
 Chinese Journal of Chemical Physics, 2017, 30(3): 277286
 化学物理学报, 2017, 30(3): 277286
 http://dx.doi.org/10.1063/16740068/30/cjcp1703025

Article history
 Received on: March 3, 2017
 Accepted on: March 17, 2017
b. Beijing Computational Science Research Center, Beijing 100094, China
The physics laws of quantum dynamics, such as the Schrös wave equation [1]and the Heisenberg matrix equation [2]for isolated systems, take the form of linear differential equations with constant coefficients.In reality, a quantum system is inevitably coupled to an environment and thus becomes open. When the coupling is weak and the memory time of the environment is short compared to the intrinsic time scale of the system, a Markovian approximation ignoring the memory time leads to a good description [35].Otherwise, the hierarchical approach is a powerful tool to simulate the nonMarkovian dynamics by embedding the nonMarkovian system in a larger, possibly infinite Markovian system [69].Once again, the equations of motion for the open systems take a linear differential form with constant coefficients if external driving fields are not applied.Formally, the dynamics for all of these cases are dictated by the ordinary differential equation
$\begin{eqnarray}\dot{{\textbf{x}}} = \mathcal{M} {\textbf{x}}\end{eqnarray}$  (1) 
together with an initial condition
The numerical integration of such an equation is simple in concept from the mathematical perspective.For systems with small to moderate dimension, one can benefit from the extensive studies on the linear differential equations [1014].The accumulated results include the Lanczos iteration in Kryslov subspace [1517], the Chebyshev scheme [18, 19], matrix exponential method [2022], and the classical RungeKutta method [23].
In practice, however, the computations are complicated by the big size of the state vector in applying the quantum principles for complex systems.As the computer memory required to store the state vector grows exponentially with the size of the system, the length $N$ of the state vector can easily reach $10^9$ and even higher in various simulations.For large systems, the quantum dynamical simulation becomes expansive or even prohibited, which is called the curse of dimensionality [24].This curse with respect to the memory issue is enhanced by the numerical integration because the integrator itself needs additional memory to store the state vector.
For this reason, lowstorage methods draw a lot of attention.The split operator method, which is a symplectic integrator not requiring additional memory location for the wavefunction, is efficient for quantum mechanical systems without spacemomentum mixed terms in the Hamiltonian [25].For other systems, lowstorage RungeKutta (LSRK) methods, which use less memory locations than the classical RungeKutta schemes, are useful in reducing the memory use [2629].Actually, the LSRK scheme was proposed at the early time of computer simulations [3032], when the memory storage was small.Later, different groups made schemes using only two memory locations for any stage number (the number of derivative evaluations in each time step of the numerical integration) [27, 29, 33].Kennedy and coworkers systematically studied the lowstorage schemes and made the lowstorage idea popular [34, 35].
Regardless of how powerful the high performance computing becomes and how big the memory storage is available, we always ask more to simulate bigger systems to reveal new mechanism, to seek new phenomena and to design new materials.A partial, alternative remedy for this problem is to develop efficient algorithm with not only shorter computational time but also less memory demand so that we can simulate bigger systems with available computational resources.Thus, the lowstorage methods became increasingly popular during the past two decades although the available computer memory was dramatically improved.
The LSRK scheme assumes the same accuracy and robustness as the traditional explicit RungeKutta methods but only needs two memory locations for the state vector regardless of the order.By contrast, the classical explicit $s$stage RungeKutta method normally requires
Actually, we can develop more efficient LSRK schemes if we take advantage of the linearity of Eq.(1). Zingg and Chisholm [38] suggested a fourstage fourthorder LSRK method when adding an inhomogeneous timedependent term to Eq.(1). Indeed, more efficient lowstorage algorithm exists for linear differential equations with constant coefficients.
In this report, we show that
In this work, we show that
The RungeKutta methods are derivativefree schemes to find the numerical solution for the initial value problem of the ordinary differential equations
$\begin{eqnarray}\dot{{\textbf{x}}}= \textbf{f}(t,\textbf{x}) \end{eqnarray}$  (2) 
$\begin{eqnarray} \textbf{x}(t_0) =\textbf{x}_0\end{eqnarray}$  (3) 
for a given function
$\begin{eqnarray} {\textbf{k}}_1={\textbf{f}}(t_n,{\textbf{x}}_n)\end{eqnarray}$  (4) 
$\begin{eqnarray}{\textbf{k}}_j={\textbf{f}}\Big(t_n+c_j \Delta t,{\textbf{x}}_n+\Delta t\sum_{l=1}^{j1}a_{jl}\textbf{k}_l\Big),\hspace{0.3cm} j =2,\cdots,s \hspace{0.5cm}\end{eqnarray}$  (5) 
$\begin{eqnarray} \textbf{x}_{n+1}=\textbf{x}_n + \Delta t\sum_{j=1}^s b_j \textbf{k}_j\end{eqnarray}$  (6) 
where
$\begin{eqnarray}c_j = \sum_{l=1}^{j1} a_{j,l},\quad j=2,\cdots,s\end{eqnarray}$  (7) 
Therefore, only
$\begin{eqnarray} \label{eq:taylor}\textbf{x}_{n+1} = \textbf{x}_n + \sum_{j=1}^p \frac{\Delta t^j}{j!}\frac{\textrm{d}^j}{\textrm{d}t^j} \textbf{f}(t,\textbf{x})\Big_{\textbf{x}=\textbf{x}_n} + O(\Delta t^{p+1})\end{eqnarray}$  (8) 
If this condition holds, the RungeKutta scheme is called a
Eq.(3) clearly shows that the
$\begin{eqnarray} && \textrm{for} \hspace{0.5cm}j=1,\cdots,s,\nonumber\\ && \textrm{do the iteration} \nonumber \\ && \hspace{0.5cm} \Delta \textbf{y} = A[j] \Delta\textbf{y} + \Delta t\textbf{f}(t_n+c_j \Delta t,\textbf{y}) \end{eqnarray}$  (9) 
$\begin{eqnarray} \textbf{y} = \textbf{y} + B_j \Delta \textbf{y} \\ \textrm{end of the iteration} \nonumber\end{eqnarray}$  (10) 
Here the register
$\begin{eqnarray} a_{jk} &=& \left\{ \begin{array}{ll} A_{k+1}a_{j,k+1} + B_k,& \textrm{if~} k<j1\\ B_k,& \textrm{if~} k=j1 \\ 0,& \textrm{otherwise} \end{array} \right. \end{eqnarray}$  (11) 
$\begin{eqnarray} b_{j} &=& \left\{ \begin{array}{ll} A_{j+1} b_{j+1} + B_j,& \textrm{if~} j<s\\ B_j,& \textrm{if~} j=s \end{array} \right.\end{eqnarray}$  (12) 
For the autonomous linear equation Eq.(1), all secondorder partial derivatives vanish, namely,
$\begin{eqnarray} \frac{\partial^2 \textbf{f}}{\partial t^2} = \frac{\partial^2 \textbf{f}}{\partial t \partial \textbf{x}} = \frac{\partial^2 \textbf{f}}{\partial \textbf{x}\partial \textbf{x}}=0\end{eqnarray}$  (13) 
As a consequence, the order condition Eq.(8) is greatly simplified,
$\begin{eqnarray}\begin{split} &\sum_{j=1}^s b_j = 1 \\ &\sum_{j,l=1}^sb_j [\textbf{a}^{k1}]_{jl} c_l = \frac{1}{k!},\quad k=1,\cdots,p1\end{split}\end{eqnarray}$  (14) 
Therefore, there are
In numerical applications, another important issue is the stability of the solver, which requires that small errors in the initial state do not drive the RungeKutta approximation off from the real value. The stability analysis of the solver is complicated and casedependent for general systems.Thus a paradigm for the stability analysis is to use the onedimensional equation with a complex coefficient
$\begin{eqnarray} \dot{x} = λ x\end{eqnarray}$  (15) 
to obtain useful guidance. The coefficient
$\begin{eqnarray} G(z) = 1+\sum_{j=1}^s \gamma_j z^j\end{eqnarray}$  (16) 
with the coefficients
$\begin{eqnarray} \gamma_j = \left\{ \begin{array}{ll} \displaystyle{\sum_{k=1}^s b_k}, & \textrm{ if } j = 1 \\ \displaystyle{\sum_{k,l=1}^s} b_k [\textbf{a}^{j2}]_{kl} c_l,& \textrm{ if } j > 1 \end{array} \right.\end{eqnarray}$  (17) 
The iteration Eqs.(46) is robust against small perturbations if the characteristic polynomial
For schemes with
Following Ref.[37], we adopt a twostage process to find the LSRK method with optimized stability. For a given stage number
In principle, we have to optimize the stability with the exact area within the contour
$\begin{eqnarray} S[G] \approx \frac{x_{\textrm{max}}}{9} \sum_{j=1}^9 y_j\end{eqnarray}$  (18) 
where
Once
$\begin{eqnarray} & & \sum_{j=1}^s b_j c_j{^2} = \frac{1}{3}\end{eqnarray}$  (19) 
$\begin{eqnarray} \sum_{j=1}^s b_j c_j{^3} = \frac{1}{4} \end{eqnarray}$  (20) 
$\begin{eqnarray} \sum_{j=1}^s b_j c_j{^4} = \frac{1}{5}\end{eqnarray}$  (21) 
$\begin{eqnarray} \sum_{j,k=1}^s b_j c_j a_{jk} c_k = \frac{1}{8} \end{eqnarray}$  (22) 
$\begin{eqnarray}&& \sum_{j,k=1}^s b_j c_j{^2} a_{jk} c_k = \frac{1}{10}\end{eqnarray}$  (23) 
$\begin{eqnarray}\sum_{j,k=1}^s b_j c_j a_{jk} c_k{^2} = \frac{1}{15}\end{eqnarray}$  (24) 
$\begin{eqnarray} & \sum_{j,k=1}^s b_j c_j [\textbf{a}^2]_{jk} c_k = \frac{1}{30}\end{eqnarray}$  (25) 
$\begin{eqnarray} \sum_{j,k=1}^s b_j a_{jk} c_k{^2} = \frac{1}{12}\end{eqnarray}$  (26) 
$\begin{eqnarray} \sum_{j,k,l=1}^s b_j a_{jk} c_k a_{kl} c_l = \frac{1}{40}\end{eqnarray}$  (27) 
$\begin{eqnarray} &\sum_{j,k=1}^s b_j [\textbf{a}^2]_{jk} c_k{^2} = \frac{1}{60} \end{eqnarray}$  (28) 
$\begin{eqnarray} \sum_{j,k,l=1}^s b_j a_{jk} c_k a_{kl} c_l = \frac{1}{20}\end{eqnarray}$  (29) 
$\begin{eqnarray}\sum_{j,k=1}^s b_j a_{jk} c_k{^3} = \frac{1}{20}\end{eqnarray}$  (30) 
The nonlinear Eq.(14) and Eqs.(1930) are solved with the LevenbergMarquardt algorithm [42, 43]as implemented in the levmar package [44]with a modified long double precision.Multiple solutions are found in determining
FIG. 1 presents the stability contours normalized by the inverse of the stage number for all the methods listed in Tables Ⅰ and Ⅱ.As the number of stages coincides the order of the methods, the stability contours of LSRK44, LSRK66, LSRK88, LSRK1010 and LSRK1212 schemes are uniquely defined and do not depend on how the parameters
We check the efficiency of all the methods with a symmetric matrix
$\begin{eqnarray}λ_j \hspace{0.15cm}&=&\hspace{0.15cm} \frac{5 i (j1) }{N}  \frac{4(j1) }{N} \textrm{exp}\left[{1\frac{4(j1) }{N}}\right] \nonumber\\&&\quad 1\leq j \leq N\end{eqnarray}$  (31) 
where
For the simplicity of computation, we choose the discrete Fourier transform matrix
$\begin{eqnarray} \textbf{x}^{(\textrm{ex})}(t) = F λ(t) F^\dagger \textbf{x}(0) \end{eqnarray}$  (32) 
where
$\begin{eqnarray}\textrm{Err}(t) = \left\sum_{j=1}^N\Big[\textbf{x}_j^{(\textrm{ex})}(t)  \textbf{x}_j^{(\textrm{appr})}(t)\Big]^2\right^{1/2}\end{eqnarray}$  (33) 
In FIG. 2, we present the error of all the methods with the step size of
Not shown in FIG. 2, we also carry out the integration with the classical fourthorder RungeKutta method, the Fehlberg's 4(5) method [45], the Calvo 6(5) method [46], and the PrinceDormand 8(7) method [47].The 4thorder LSRK is identical to the classical fourthorder RungeKutta method in errors as well as the stability.Higherorder LSRK schemes are faster than the same order classical RK ones, but no more than a factor of two.However, a highorder LSRK (
Provided with a preselected error tolerance, higher order methods normally can take a bigger step size.For example, with the error tolerance of
Again, to have comparison, we consider the LSRK134(4) scheme which is designed for general differential equations.With a given step size, it systematically produces smaller errors than that of the LSRK44 method.However, provided the same accuracy (
We will apply the highorder lowstorage methods outlined in Section Ⅱ to a more realistic simulation, that is, the nonMarkovian exciton diffusion in a linear aggregate.The exciton dynamics plays a crucial role in organic photovoltaic cell [48], photosynthesis [49], organic lightemitting diodes [50]and many others [51, 52].In such systems, energy transfers to the reaction center along with the exciton diffusion.The natural photosynthesis systems exhibit surprisingly high efficiency in the exciton energy transfer [49].The theoretical understanding of the design features and the underline mechanism will greatly help the artificial solar energy application [53].We are still at the very beginning to understand the process and now expect that dissipation caused by the environment is the key for efficient energy transport [49, 53].
Here, we study the noiseassisted exciton diffusion in an artificial chain consisting of
As the time scale of the exciton dynamics is of the order of picosecond and longer, it is comparable to that of the environment. Thus, the dynamics is essentially nonMarkovian. Here, the nonMarkovian dynamics is investigated with the hierarchical approach, which introduces auxiliary density matrices to characterize the effect of the bath and is one of a few available methods to study nonMarkovianity [69, 6164].Usually, the number of auxiliary matrices are large and the requirement for the memory storage is high at low temperatures.The lowstorage scheme thus is ideal to integrate the hierarchical equations.
A. Hierarchical approach for dissipative exciton dynamicsIn order to focus on the exciton diffusion, the exciton is supposed to be stable and its dissociation as well as the decay is not taken into account.Moreover, we only investigate the dynamics in the oneexciton manifold, that is, only the diabatic electronic state corresponding to the
The Hamiltonian for the exciton diffusion under dissipation can be approximated with the systemplusbath model,
$\begin{eqnarray}\hat{H}_{\rm tot} = \hat{H}_{\textrm{ex}}+ \sum_m \hat{H}_m + \hat{H}_{\rm b}\end{eqnarray}$  (34) 
where
$\begin{eqnarray}\hat{H}_{\textrm{ex}} = \sum_{mn} J_{mn}m\rangle\langle n\end{eqnarray}$  (35) 
where
When the size of the monomer in the aggregate is big, each site may have a set of localized vibrations coupled to the exciton at this site only.Operator
$\begin{eqnarray}\hat{H}_m &=& \hat{h}_m \sum_\xi \hat{q}_{m,\xi}g_{m,\xi},\quad \hat{h}_m \equivm\rangle\langle m\end{eqnarray}$  (36) 
For linear dissipation, the effect of the vibrations on the exciton dynamics is fully captured by the bath response function,
$\begin{eqnarray} \alpha(t) = \frac{1}{\pi} \int^\infty_0 \textrm{d}\omega J(\omega) \left[\coth\left(\frac{\beta\hbar\omega}{2}\right) \cos(\omega t) i\sin(\omega t )\right]\end{eqnarray}$  (37) 
which in turn is determined by the spectral density function
$\begin{eqnarray}J(\omega) =\frac{\eta \omega{_0}^3\omegaλ} {[(\omega^2  \omega{_0}^2) ^2 + λ^2\omega^2]}\end{eqnarray}$  (38) 
where
For this spectral density the response function can be expanded into a sum of exponentials as follows [62]
$\begin{eqnarray}\alpha(t) &=& \frac{\eta\hbar\omega{_0}^3}{4\xi}\Big\{\textrm{e}^{\Omega_1 t} \Big[\coth\Big(\frac{i}{2}\hbar\beta\Omega_1\Big)  1\Big] \nonumber \\&&+\textrm{e}^{\Omega_2 t} \Big[\coth\Big(\frac{i}{2}\hbar\beta\Omega_2\Big)  1\Big] \Big\} \nonumber \\ &&  2 \textrm{e}\omega_0^3λ^2\hbar \beta \sum_{k=1}^\infty \frac{\nu_k\textrm{e}^{\nu_k t}}{(\nu_k^2 + \omega_0^2) ^2  λ^2\nu_k^2} \nonumber \\ &\equiv& \sum_{j=1}^{N_\textrm{e}} c_j \exp(\Omega_j t) + \alpha^{\textrm{(short)}}(t)\end{eqnarray}$  (39) 
where $\Omega_{1,2}=λ/2$±$i \xi$ with
With the exponential decomposition of the response functions, the hierarchical approach is efficient to propagate the reduced density matrix for the description of the nonMarkovian exciton dynamics.The hierarchical approach is a scheme to attack the nonperturbative nonMarkovian dynamics for the open quantum system by converting the memory of the bath into a series of auxiliary density matrices.Here we will follow the stochastic decoupling procedure of Shao and coworkers [9, 67]. Derivation based on fundamental Itô calculus can be also found somewhere else [68, 69]. The resulting equations are as follows
$\begin{array}{l} i\hbar \frac{{\partial {{\tilde \rho }_A}}}{{\partial t}} =  i\hbar (\sum\limits_{m = 1}^{{N_{{\rm{agg}}}}} {\sum\limits_{j = 1}^{{N_{\rm{e}}}} {{A_{m,j}}} } {\Omega _j} + {\cal R}){{\tilde \rho }_A}\\ \quad + \left[{{{\hat H}_{{\rm{ex}}}},{{\tilde \rho }_A}} \right] + \sum\limits_{m = 1}^{{N_{{\rm{agg}}}}} {\sum\limits_{j = 1}^{{N_{\rm{e}}}} {\left[{{{\hat h}_m},{{\tilde \rho }_{A_{m,j}^ + }}} \right]} } \\ \quad + \sum\limits_{m = 1}^{{N_{{\rm{agg}}}}} {\sum\limits_{j = 1}^{{N_{\rm{e}}}} {{A_{m,j}}} } {\rm{(}}{\alpha _{j; + }}\left[{{{\hat h}_m},{{\tilde \rho }_{A_{m,j}^  }}} \right] + \\ \quad {\alpha _{j;  }}\left\{ {{{\hat h}_m},{{\tilde \rho }_{A_{m,j}^  }}} \right\}{\rm{)}} \end{array}$  (40) 
$\begin{eqnarray}\alpha_{j;+}=\frac{c_j + c_{j^\prime}^\ast}{2},\quad \quad\alpha_{j;}=\frac{ i (c_j  c_{j^\prime}^\ast)}{2}\end{eqnarray}$  (41) 
with
In the above equation
Eq.(40) is an infinitedimensional differential equation with constant coefficients.In numerical simulations, it is an unattainable yet unnecessary task to explicitly treat an infinite system.In practice, one can safely truncate the hierarchy equation to the
As far as the present model is concerned, we will adopt the same parameters as those in Ref.[62].Specifically, in the excitonic part we will only consider the nearest neighbor coupling and use the calculated value
There are two key factors for the hierarchical approach to obtain converged results.One is the number of exponentials
With these settings, there are totally 1, 674, 992 auxiliary density matrices in the hierarchy and it takes about 5.6 GB memory to store all complex matrices.A LSRK method with any order will need 11.2 GB memory storage for the numerical propagation.By contrast, the classical fourthorder RungeKutta method requires 16.8 GB space; the fifthorder RungeKuttaFehlberg method needs 39.2 GB [45].
The numerical simulations are done with the Hyshe package [70], which employs a shared memory parallelization.Specifically, we use a 24core Intel(R) Xeon(R) E52650v3 CPU with a clock speed of
The exciton population at site one is shown in FIG. 3.The population follows a steady equilibration after an initial quick drop. Note that there exist three clear spikes around 0.24, 0.47, 0.69 ps, and minor ones around 0.92, 1.2, and 1.5 ps.The peak heights decrease as time elapses, but the peak separation is almost equal in time.One may ask why this quasiperiodic phenomenon emerges?
To answer this question, we scrutinize the exciton population at all sites and present the data in FIG. 4. The results reveal that the exciton transfers in the aggregate with a constant rate and is bounced by the two ends of the aggregate.Together with the dissipation induced by the vibrations, these two factors lead to an underdamped exciton motion.Based on this observation, the exciton hopping rate is about 123.4 ps
The wavelike motion of the exciton transfer can also be manifested by the coherence beating. FIG. 5plots the absolute value of the offdiagonal matrix elements
Previously, lowstorage RungeKutta methods were designed for integration of general nonlinear or timevariant systems, which is efficient in the use of the memory storage.To generate a lowstorage scheme for general equations, we have to take additional stages because there are more restrictions for the method.However, the dynamics of many quantum systems is governed by linear differential equations with constant coefficients.For these autonomous linear equations, there are only
Numerics shows that the 12stage 12thorder LSRK method could be an order of magnitude faster than the fourthorder RungeKutta method.
As an example, we have applied the suggested 12stage12thorder method for the hierarchical simulations of the exciton diffusion in a 15site linear aggregate, where it takes 5.6 GB memory to store the reduced and auxiliary density matrices. Thanks to the efficient lowstage scheme, we can propagate the nonMarkovian dynamics up to 2.4 ps within 6.5 h with a 24core Intel(R) Xeon(R) E52650v3 CPU.
Ⅴ. ACKNOWLEDGMENTSThis work is supported by the National Natural Science Foundation of China (No.21373064), the Program for Innovative Research Team of Guizhou Province (No.QKTD[2014]4021), and the Natural Scientific Foundation from Guizhou Provincial Department of Education (No.ZDXK[2014]18).All the calculations were performed at Guizhou Provincial HighPerformance Computing Center of Condensed Materials and Molecular Simulation in Guizhou Education University.
[1]  E. Schrödinger, Ann. Phys. 384 , 361 (1926). DOI:10.1002/(ISSN)15213889 
[2]  M. Born, W. Heisenberg, and P. Jordan, Zeit. Phys. 35 , 557 (1925). 
[3]  F. Bloch, Phys. Rev. 70 , 460 (1946). DOI:10.1103/PhysRev.70.460 
[4]  A. G. Redfield, IBM J. Res. Dev. 1 , 19 (1957). DOI:10.1147/rd.11.0019 
[5]  G. Lindblad, Commun. Math. Phys. 48 , 119 (1976). DOI:10.1007/BF01608499 
[6]  V. Shapiro, and V. Loginov, Physica A 91 , 563 (1978). DOI:10.1016/03784371(78)90198X 
[7]  Y. Tanimura, and R. Kubo, J. Phys. Soc. Jpn. 58 , 101 (1989). DOI:10.1143/JPSJ.58.101 
[8]  Y. Tanimura, Phys. Rev. A 41 , 6676 (1990). DOI:10.1103/PhysRevA.41.6676 
[9]  Y. A. Yan, F. Yang, Y. Liu, and J. Shao, Chem. Phys. Lett. 395 , 216 (2004). DOI:10.1016/j.cplett.2004.07.036 
[10]  L. Adrianova, Introduction to Linear Systems of Diffierential Equations, Translations of Mathematical Monographs, Vol. 146, Providence, Rhode Island: AMS, (1995). 
[11]  A. Greenbaum, Iterative Methods for Solving Linear Systems, Frontiers in Applied Mathematics, Vol. 17, Philadelphia: SIAM, (1997). 
[12]  E. Coddington and R. Carlson, Linear Ordinary Diffierential Equations. Philadelphia: SIAM, (1997). 
[13]  A. Forsyth, Part 3. Ordinary Linear Equations, Theory of Diffierential Equations Vol. 4, Cambridge: Cambridge University Press (1902). 
[14]  Y. Saad, Iterative Methods for Sparse Linear Systems 2nd Edn. , Philadelphia: SIAM, (2003). 
[15]  C. Lanczos, J. Res. Nat. Bur. Standards 45 , 255 (1950). DOI:10.6028/jres.045.026 
[16]  R. Haydock, Comput. Phys. Commun. 20 , 11 (1980). DOI:10.1016/00104655(80)901010 
[17]  V. Druskin, and L. Knizhnerman, J. Comput. Appl. Math. 50 , 255 (1994). DOI:10.1016/03770427(94)903050 
[18]  H. TalEzer, and R. Kosloff, J. Chem. Phys. 81 , 3967 (1984). DOI:10.1063/1.448136 
[19]  C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, J. Comput. Phys. 94 , 59 (1991). DOI:10.1016/00219991(91)90137A 
[20]  R. C. Ward, SIAM J. Numer. Anal. 14 , 600 (1977). DOI:10.1137/0714039 
[21]  M. Hochbruck, C. Lubich, and H. Selhofer, SIAM J. Sci. Comput. 19 , 1552 (1998). DOI:10.1137/S1064827595295337 
[22]  R. B. Sidje, ACM Trans. Math. Softw. 24 , 130 (1998). DOI:10.1145/285861.285868 
[23]  J. C. Butcher, Numerical Methods for Ordinary Diffierential Equations, 2nd Edn. Chichester: Wiley, (2008). 
[24]  R. E. Bellman, Dynamic Programming. Princeton: Princeton University Press (1957). 
[25]  M. Feit, J. Fleck, and A. Steiger, J. Comput. Phys. 47 , 412 (1982). DOI:10.1016/00219991(82)900912 
[26]  P. Houwen, Numer. Math. 20 , 149 (1972). DOI:10.1007/BF01404404 
[27]  P. J. Van Der Houwen, Construction of Integration Formulas for Initial Value Problems, Vol. 19, Amsterdam: NorthHolland, (1977). 
[28]  A. Wray, NASA Ames Research Center, California: Moffett Field, 202(1990). 
[29]  J. Williamson, J. Comput. Phys. 35 , 48 (1980). DOI:10.1016/00219991(80)900339 
[30]  S. Gill, Math. Proc. Cambridge Phil. Soc. 47 , 96 (1951). DOI:10.1017/S0305004100026414 
[31]  E. K. Blum, Math. Comput. 16 , 176 (1962). DOI:10.1090/S00255718196201456614 
[32]  D. J. Fyfe, Math. Comput. 20 , 392 (1966). DOI:10.1090/S00255718196602023170 
[33]  D. I. Ketcheson, SIAM J. Sci. Comput. 30 , 2113 (2008). DOI:10.1137/07070485X 
[34]  M. H. Carpenter and C. A. Kennedy, Fourthorder 2Nstorage RungeKutta schemes, Tech. Rep. NASATM109112, VA: National Aeronautics and Space Administration, NASA Langley Research Center, (1994). 
[35]  C. A. Kennedy, M. H. Carpenter, and R. Lewis, Appl. Numer. Math. 35 , 177 (2000). DOI:10.1016/S01689274(99)001415 
[36]  J. Butcher, BIT Numer. Math. 25 , 521 (1985). DOI:10.1007/BF01935372 
[37]  J. Niegemann, R. Diehl, and K. Busch, J. Comput. Phys. 231 , 364 (2012). DOI:10.1016/j.jcp.2011.09.003 
[38]  T. Becker, H. Kredel, and V. Weispfenning, Gröbner Bases: a Computational Approach to Commutative Algebra, Corrected Edn. , Graduate Texts in Mathematics, Vol. 141, London, UK: SpringerVerlag, (1993). 
[39]  D. Zingg, and T. Chisholm, Appl. Numer. Math. 31 , 227 (1999). DOI:10.1016/S01689274(98)001299 
[40]  J. A. Nelder, and R. Mead, Comput. J. 7 , 308 (1965). DOI:10.1093/comjnl/7.4.308 
[41]  S. Johnson, The Nlopt NonlinearOptimization Package. 
[42]  K. Levenberg, Quart. J. Appl. Math. 2 , 164 (1944). DOI:10.1090/qam/19440202 
[43]  D. W. Marquardt, J. Soc. Indust. Appl. Math. 11 , 431 (1963). DOI:10.1137/0111030 
[44]  M. Lourakis, Levmar: LevenbergMarquardt Nonlinear Least Squares Algorithms in C/C++, http://www.ics.forth.gr/~lourakis/levmar/July (2004). 
[45]  E. Fehlberg, Computing 6 , 61 (1970). DOI:10.1007/BF02241732 
[46]  M. Calvo, J. Montijano, and L. Randez, Comput. Math. Appl. 20 , 15 (1990). 
[47]  P. Prince, and J. Dormand, J. Comput. Appl. Math. 7 , 67 (1981). 
[48]  C. W. Tang, Appl. Phys. Lett. 48 , 183 (1986). DOI:10.1063/1.96937 
[49]  H. Lee, Y. C. Cheng, and G. R. Fleming, Science 316 , 1462 (2007). DOI:10.1126/science.1142188 
[50]  J. H. Burroughes, D. D. C. Bradley, A. R. Brown, R. N. Marks, K. Mackay, R. H. Friend, P. L. Burns, and A. B. Holmes, Nature 347 , 539 (1990). DOI:10.1038/347539a0 
[51]  G. Trinkunas, J. L. Herek, T. Polıvka, V. Sundström, and T. Pullerits, Phys. Rev. Lett. 86 , 4167 (2001). DOI:10.1103/PhysRevLett.86.4167 
[52]  V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd Edn. Weinheim: WILEYVCH, (2010). 
[53]  del Rey M., A. W. Chin, S. F. Huelga, and M. B. Plenio, J. Phys. Chem. Lett. 4 , 903 (2013). DOI:10.1021/jz400058a 
[54]  Z. Chen, V. Stepanenko, V. Dehm, P. Prins, L. Siebbeles, J. Seibt, P. Marquetand, V. Engel, and F. Wrthner, Chem. Eur. J. 13 , 436 (2007). DOI:10.1002/(ISSN)15213765 
[55]  X. Q. Li, X. Zhang, S. Ghosh, and F. Würthner, Chem. Eur. J. 14 , 8074 (2008). DOI:10.1002/chem.v14:27 
[56]  M. R. Wasielewski, Acc. Chem. Res. 42 , 1910 (2009). DOI:10.1021/ar9001735 
[57]  S. Wolter, J. Aizezers, F. Fennel, M. Seidel, F. Wurthner, O. Kühn, and S. Lochbrunner, New J. Phys. 14 , 105027 (2012). DOI:10.1088/13672630/14/10/105027 
[58]  H. Marciniak, X. Q. Li, F. Würthner, and S. Lochbrunner, J. Phys. Chem. A 115 , 648 (2011). 
[59]  D. Ambrosek, H. Marciniak, S. Lochbrunner, J. Tatchen, X. Li, F. Würthner, and O. Kühn, Phys. Chem. Chem. Phys. 13 , 17649 (2011). DOI:10.1039/c1cp21624d 
[60]  D. Ambrosek, A. Köhn, J. Schulze, and O. Kühn, J. Phys. Chem. A 116 , 11451 (2012). DOI:10.1021/jp3069706 
[61]  Q. Shi, L. Chen, G. Nan, R. X. Xu, and Y. Yan, J. Chem. Phys. 130 , 084105 (2009). DOI:10.1063/1.3077918 
[62]  Y. A. Yan, and O. Kühn, New J. Phys. 14 , 105004 (2012). DOI:10.1088/13672630/14/10/105004 
[63]  Y. A. Yan, and S. Cai, J. Chem. Phys. 141 , 054105 (2014). DOI:10.1063/1.4891798 
[64]  Y. Yan, J. Chem. Phys. 144 , 024305 (2016). DOI:10.1063/1.4939523 
[65]  O. Kühn, and V. Sundström, J. Chem. Phys. 107 , 4154 (1997). DOI:10.1063/1.474803 
[66]  H. Grabert, P. Schramm, and G. L. Ingold, Phys. Rep. 168 , 115 (1988). DOI:10.1016/03701573(88)900233 
[67]  J. Shao, J. Chem. Phys. 120 , 5053 (2004). DOI:10.1063/1.1647528 
[68]  M. Schröter, S. Ivanov, J. Schulze, S. Polyutov, Y. Yan, T. Pullerits, and O. Kühn, Phys. Rep. 567 , 1 (2015). DOI:10.1016/j.physrep.2014.12.001 
[69]  Y. Yan, and J. Shao, Front. Phys. 11 , 110309 (2016). DOI:10.1007/s1146701605709 
[70]  Y. A. Yan and Y. Zhou, Hybrid StochasticHierarchical Equations. http://nano.gznc.edu.cn/~yunan/hyshe.html (2012). 
b. 北京计算科学研究中心, 北京 100094