Numerical Assessment for Accuracy and GPU Acceleration of TD-DMRG Time Evolution Schemes

TL;DR

This paper evaluates the accuracy of three TD-DMRG time evolution schemes and demonstrates significant GPU acceleration, enhancing the efficiency of quantum dynamics simulations for complex systems.

Contribution

It provides a comparative analysis of three TD-DMRG schemes and introduces GPU acceleration to significantly speed up tensor contractions in these methods.

Findings

TDVP-PS allows larger time steps than TDVP-MU.

TDVP-MU and TDVP-PS are more accurate than P&C-RK4.

GPU acceleration speeds up TDVP schemes by up to 73 times.

Abstract

Time dependent density matrix renormalization group (TD-DMRG) has become one of the cutting edge methods of quantum dynamics for complex systems. In this paper, we comparatively study the accuracy of three time evolution schemes in TD-DMRG, the global propagation and compression method with Runge-Kutta algorithm (P&C-RK), the time dependent variational principle based methods with matrix unfolding algorithm (TDVPMU) and with projector-splitting algorithm (TDVP-PS), by performing benchmarks on the exciton dynamics of Fenna-Matthews-Olson (FMO) complex. We show that TDVP-MU and TDVP-PS yield the same result when the time step size is converged and they are more accurate than P&C-RK4, while TDVP-PS tolerates a larger time step size than TDVP-MU. We further adopt the graphical processing units (GPU) to accelerate the heavy tensor contractions in TD-DMRG and it is able to speed up the TDVP-MU and TDVP-PS schemes by up to 73 times.

Tables

Table 1. Table 1: Summary of the computational complexity of P&C-RK4, TDVP-MU and TDVP-PS time evolution schemes.

Sub-step 111See the text in Section II for the definitions and labels of the sub-steps.	P&C-RK4	TDVP-MU	TDVP-PS
$\hat{H}|\Psi ⟩$	$M_{\text{S}}^2{M}_{\text{O}}^2{d}^2$	—	—
QR	$M_{\text{S}}^3{M}_{\text{O}}^{3}d$	—	$M_{\text{S}}^{3}d$
MatMul-QR	$M_{\text{S}}^3{M}_{\text{O}}^{3}d$	—	$M_{\text{S}}^{3}d$
SVD	$\min (M_{\text{S}}^3{M}_{\text{O}}^{2}d,M_{\text{S}}^3{M}_{\text{O}}{d}^2)$	$M_{\text{S}}^{3}d$	—
MatMul-SVD	$M_{\text{S}}^3{M}_{\text{O}}^{2}d$	$M_{\text{S}}^{3}d$	—
Get Env	—	$M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d$	$M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d$
Deriv	—	$M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d+M_{\text{S}}^3{d}^2$	$M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d$ (forward) $M_{\text{S}}^3{M}_{\text{O}}$ (backward)
Overall	$N{M}_{\text{S}}^3{M}_{\text{O}}^{3}d$	$N\left(M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d+M_{\text{S}}^3{d}^2\right)$ (VMF) $Nf\left(M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d+M_{\text{S}}^3{d}^2\right)$ (CMF)	$Nf\left(M_{\text{S}}^2{M}_{\text{O}}^2{d}^2+M_{\text{S}}^3{M}_{\text{O}}d\right)$

Table 2. Table 2: The mean cumulative deviation with various M S subscript 𝑀 S M_{\textrm{S}} and τ = 10 𝜏 10 \tau=10 by TDVP-PS at t 1 = 6560 subscript 𝑡 1 6560 t_{1}=6560 and t 2 = 41120 subscript 𝑡 2 41120 t_{2}=41120 . (Reference: M S = 256 , τ = 10 formulae-sequence subscript 𝑀 S 256 𝜏 10 M_{\textrm{S}}=256,\tau=10 )

	error ($t$)
$M_{\text{S}}$	$t_1=6560$	$t_2=41120$
32	$6\times {10}^{-7}$	$1\times {10}^{-3}$
64	$8\times {10}^{-9}$	$3\times {10}^{-4}$
128	$3\times {10}^{-10}$	$3\times {10}^{-5}$
192	$3\times {10}^{-11}$	$7\times {10}^{-6}$

Table 3. Table 3: The mean cumulative deviation with M S = 256 subscript 𝑀 S 256 M_{\textrm{S}}=256 and various τ 𝜏 \tau by TDVP-PS at t 1 = 6560 subscript 𝑡 1 6560 t_{1}=6560 and t 2 = 41120 subscript 𝑡 2 41120 t_{2}=41120 . (Reference: M S = 256 , τ = 10 formulae-sequence subscript 𝑀 S 256 𝜏 10 M_{\textrm{S}}=256,\tau=10 )

	error ($t$)
$\tau$	$t_1=6560$	$t_2=41120$
320	$4\times {10}^{-10}$	$8\times {10}^{-7}$
160	$7\times {10}^{-11}$	$7\times {10}^{-7}$
80	$9\times {10}^{-13}$	$8\times {10}^{-7}$
40	$7\times {10}^{-13}$	$6\times {10}^{-7}$

Equations

$$|\Psi ⟩=\sum _{\{a\},\{\sigma \}}{A}_{a_1}^{{\sigma }_1}{A}_{a_1{a}_2}^{{\sigma }_2}\cdots A_{a_{N-1}}^{{\sigma }_N}|{\sigma }_1{\sigma }_2\cdots {\sigma }_N⟩$$

$$|a_i[1:i]⟩$$

$$|a_j[j+1:N]⟩$$

$$S[1:i{]}_{a_i{a}_i^{\prime }}=⟨a_i[1:i]|a_i^{\prime }[1:i]⟩$$

$$S[j+1:N{]}_{a_j{a}_j^{\prime }}=⟨a_j[j+1:N]|a_j^{\prime }[j+1:N]⟩$$

$$|\Psi ⟩=\sum _{\{l\},\{r\}\{\sigma \}}$$

$$\sum _{{\sigma }_i,l_{i-1}}{L}_{l_{i-1}{l}_i^{\prime }}^{{\sigma }_i\ast }{L}_{l_{i-1}{l}_i}^{{\sigma }_i}$$

$$\sum _{{\sigma }_j,r_j}{R}_{r_{j-1}^{\prime }{r}_j}^{{\sigma }_j\ast }{R}_{r_{j-1}{r}_j}^{{\sigma }_j}$$

$$C_{l_{n-1}{r}_n}^{{\sigma }_n}\stackrel{\text{reshape}}{\to }{C}_{{\sigma }_n{l}_{n-1},r_n}\stackrel{\text{QR}}{=}\sum _{l_n}{L}_{{\sigma }_n{l}_{n-1},l_n}{D}_{l_n,r_n}\stackrel{\text{reshape}}{\to }\sum _{l_n}{L}_{l_{n-1}{l}_n}^{{\sigma }_n}{D}_{l_n{r}_n}$$

$$C_{l_{n-1}{r}_n}^{{\sigma }_n}\stackrel{\text{reshape}}{\to }{C}_{{\sigma }_n{l}_{n-1},r_n}\stackrel{\text{SVD}}{=}\sum _{s=1}^k{U}_{{\sigma }_n{l}_{n-1},s}{\Lambda }_{ss}{V}_{s,r_n}^{†}\stackrel{\text{reshape}}{\to }\sum _{s=1}^k{U}_{l_{n-1}s}^{{\sigma }_n}{\Lambda }_{ss}{V}_{s{r}_n}^{†}\approx \sum _{l_n=1}^{M_{\text{S}}<k}{L}_{l_{n-1}{l}_n}^{{\sigma }_n}{\Lambda }_{l_n{l}_n}{V}_{l_n{r}_n}^{†}$$

$$\hat{O}=\sum _{\{w\},\{\sigma \},\{{\sigma }^{\prime }\}}{W}_{w_1}^{{\sigma }_1^{\prime },{\sigma }_1}{W}_{w_1{w}_2}^{{\sigma }_2^{\prime },{\sigma }_2}\cdots W_{w_{N-1}}^{{\sigma }_N^{\prime },{\sigma }_N}|{\sigma }_1^{\prime }{\sigma }_2^{\prime }\cdots {\sigma }_N^{\prime }⟩⟨{\sigma }_N{\sigma }_{N-1}\cdots {\sigma }_1|$$

$$\hat{O}|\Psi ⟩=\sum _{\{w,a\},\{\sigma \}}{A}_{\{w,a{\}}_1}^{\prime {\sigma }_1^{\prime }}{A}_{\{w,a{\}}_1\{w,a{\}}_2}^{\prime {\sigma }_2^{\prime }}\cdots A_{\{w,a{\}}_{N-1}}^{\prime {\sigma }_N^{\prime }}|{\sigma }_1^{\prime }{\sigma }_2^{\prime }\cdots {\sigma }_N^{\prime }⟩$$

$$A_{\{w,a{\}}_{i-1}\{w,a{\}}_i}^{\prime {\sigma }_i^{\prime }}=\sum _{{\sigma }_i}{W}_{w_{i-1}{w}_i}^{{\sigma }_i^{\prime },{\sigma }_i}{A}_{a_{i-1}{a}_i}^{{\sigma }_i}$$

$$|\Psi ⟩+|\Phi ⟩=\sum _{\{a,b\},\{\sigma \}}{A}_{\{a,b{\}}_1}^{\prime {\sigma }_1}{A}_{\{a,b{\}}_1\{a,b{\}}_2}^{\prime {\sigma }_2}\cdots A_{\{a,b{\}}_{N-1}}^{\prime {\sigma }_N}|{\sigma }_1{\sigma }_2\cdots {\sigma }_N⟩$$

$$A^{\prime {\sigma }_1}=\left[\begin{array}{cc}{A}^{{\sigma }_1}& B^{{\sigma }_1}\end{array}\right],A^{\prime {\sigma }_i}=\left[\begin{array}{cc}{A}^{{\sigma }_i}& 0\\ 0& B^{{\sigma }_i}\end{array}\right](i=2,3\cdots N-1),A^{\prime {\sigma }_N}=\left[\begin{array}{c}{A}^{{\sigma }_N}\\ B^{{\sigma }_N}\end{array}\right]$$

$$|k_1⟩$$

$$|k_2⟩$$

$$|k_3⟩$$

$$|k_4⟩$$

$$|\Psi (t+\tau )⟩$$

$$⟨\delta \Psi ⟩i\text{partialderivative}t-\hat{H}\Psi =0$$

$$\text{partialderivative}|\Psi ⟩t=-i\hat{P}\hat{H}|\Psi ⟩$$

$$\hat{P}=\sum _{i=1}^N\hat{P}[1:i-1]\otimes {\hat{I}}_i\otimes \hat{P}[i+1:N]-\sum _{i=1}^{N-1}\hat{P}[1:i]\otimes \hat{P}[i+1:N]$$

$$\hat{P}[1:i]$$

$$\hat{P}[i+1:N]$$

$${\hat{I}}_i$$

$$\hat{P}[1:0]$$

$$\hat{P}$$

$$\hat{Q}[1:i]$$

$$\displaystyle\cdot\big{(}S[1:i-1]^{-1}_{a^{\prime}_{i-1}a_{i-1}}\delta_{\sigma^{\prime}_{i}\sigma_{i}}-\sum_{a^{\prime}_{i},a_{i}}A^{\sigma^{\prime}_{i}}_{a^{\prime}_{i-1}a^{\prime}_{i}}S[1:i]^{-1}_{a^{\prime}_{i}a_{i}}A^{\sigma_{i}*}_{a_{i-1}a_{i}}\big{)}$$

$$\hat{Q}[i:N]$$

$$\displaystyle\cdot\big{(}S[i+1:N]^{-1}_{a^{\prime}_{i}a_{i}}\delta_{\sigma^{\prime}_{i}\sigma_{i}}-\sum_{a^{\prime}_{i-1},a_{i-1}}A^{\sigma^{\prime}_{i}}_{a^{\prime}_{i-1}a^{\prime}_{i}}S[i:N]^{-1}_{a^{\prime}_{i-1}a_{i-1}}A^{\sigma_{i}*}_{a_{i-1}a_{i}}\big{)}$$

$$i\text{partialderivative}{C}_{l_{N-1}^{\prime }}^{{\sigma }_N^{\prime }}t=\sum _{{\sigma }_N^{\prime \prime },l_{N-1}^{\prime \prime }}H[N{]}_{l_{N-1}^{\prime }{\sigma }_N^{\prime },l_{N-1}^{\prime \prime }{\sigma }_N^{\prime \prime }}{C}_{l_{N-1}^{\prime \prime }}^{{\sigma }_N^{\prime \prime }}$$

$$i\text{partialderivative}{L}_{l_{i-1}^{\prime }{l}_i^{\prime }}^{{\sigma }_i^{\prime }}t=$$

$$H[i{]}_{l_{i-1}^{\prime }{\sigma }_i^{\prime }{l}_i^{\prime },l_{i-1}{\sigma }_i{l}_i}=$$

$$h[1:i-1{]}_{\{l^{\prime },w,l{\}}_{i-1}}=$$

$$h[i+1:N{]}_{\{l^{\prime },w,l{\}}_i}=$$

$$h[i{]}_{\{l^{\prime },w,l{\}}_{i-1},\{l^{\prime },w,l{\}}_i}=$$

$$p[i{]}_{l_{i-1}^{\prime }{\sigma }_i^{\prime },l_{i-1}{\sigma }_i}=$$

$$\text{partialderivative}t(\sum _{l_{i-1},{\sigma }_i}{L}_{l_{i-1}{l}_i}^{{\sigma }_i\ast }{L}_{l_{i-1}{l}_i^{\prime }}^{{\sigma }_i})=0$$

Full text

††preprint: AIP/123-QED

I Introduction

Time dependent density matrix renormalization group (TD-DMRG) has emerged as a powerful tool to deal with many-body chemical and physical problems Ma, Luo, and Yao (2018); Paeckel et al. (2019). Although DMRG is initially designed to solve the ground state of one dimensional strongly correlated systems White (1992, 1993), its applications are successfully extended to dynamical properties both in the time and frequency domains, such as linear and nonlinear optical response of polyenes Shuai et al. (1998), polaron formation and diffusion Zhang, Jeckelmann, and White (1999); Yao (2017); Barford and Mannouch (2018), interconversion dynamics of pyrazine Greene and Batista (2017); Baiardi and Reiher (2019); Xie et al. (2019), ab initio electron dynamics in hydrogen chain Ronca et al. (2017), exciton dissociationYao et al. (2018), spectra of molecular aggregates Ren, Shuai, and Kin-Lic Chan (2018) and many other topics Chin et al. (2013); Schröder et al. (2019); Borrelli and Gelin (2017); Kloss and Lev (2019); Kloss, Reichman, and Tempelaar (2019); Frahm and Pfannkuche (2019). One of the key components in TD-DMRG is the time evolution scheme, which is essential to the numerical accuracy and efficiency. The available schemes could be roughly classified into three groups. The first group is based on globally approximating the formal propagator $e^{-iHt}$ or $e^{-iHt}\ket{\Psi}$, including time-evolving block decimation (TEBD) Vidal (2004); White and Feiguin (2004); Daley et al. (2004), WI,II method Zaletel et al. (2015), Runge-Kutta García-Ripoll (2006); Ren, Shuai, and Kin-Lic Chan (2018), Chebyshev expansion Halimeh, Kolley, and McCulloch (2015), Krylov subspace  García-Ripoll (2006); Wall and Carr (2012) methods and split operator method on the grid basis Greene and Batista (2017). The same feature shared in these schemes is that in each time step the wavefunction is firstly propagated as a whole globally and then compressed. The second group is more inspired by the original DMRG, which is formulated in the local renormalized space and the basis is adapted by the averaged reduced density matrix. The representatives are time step targeting method (TST)  Feiguin and White (2005) and some related variants Dutta and Ramasesha (2010); Ronca et al. (2017). The third group is based on the time dependent variational principle (TDVP) Dirac (1930). Depending on the different ways to derive the equations of motion (EOMs), this group includes original method with fixed gauge freedom  Haegeman et al. (2011) and the more recent projector-splitting method (PS) from a tangent space view  Haegeman et al. (2016). Among the above evolution schemes, all schemes are suited to long-range interactions except TEBD. In addition, the global evolution scheme is the most straightforward one when the modern framework of matrix product state / matrix product operator (MPS/MPO) is investigated, while the PS scheme seems to have become the most popular choice as it has been widely employed in the recent articles Borrelli and Gelin (2017); Schröder et al. (2019); Baiardi and Reiher (2019); Kloss and Lev (2019); Kloss, Reichman, and Tempelaar (2019); Kurashige (2018); Xie et al. (2019) and implemented in a number of TD-DMRG packages Hauschild and Pollmann (2018); Mendl (2018). Although many evolution schemes have been applied extensively to the simulation of quantum dynamics in both chemistry and physics, a pragmatic analysis of their relative accuracy and efficiency has not been reported yet.

High performance computing for DMRG algorithms has attracted much interest in recent years, including parallelization strategies with message passing interface (MPI) Chan (2004), open multi-processing (OpenMP) and hybrid MPI/OpenMPKurashige and Yanai (2009) on multiple central processing units (CPUs). Meanwhile, the application of graphical processing units (GPUs) in computational chemistry has also drawn much attention in the past decade, due to the tardy improvement of CPUs, and the more and more vibrant GPU software ecosystem, including electronic structure calculation Song and Martínez (2017); Kaliman and Krylov (2017); Shee et al. (2018); Liu et al. (2019), classical/ab initio molecular dynamics Penfold (2017); Lee et al. (2018), and open system quantum dynamics Kreisbeck et al. (2011); Tsuchimoto and Tanimura (2015). To the best of our knowledge, the only attempt to employ GPUs in DMRG was made by Nemes et al. in 2014 Nemes et al. (2014). They came up with a smart implementation exploiting both CPU and GPU which speeds up the Davidson diagonalization part in DMRG algorithm by 2 to 5 times. Despite their efforts, how GPUs can accelerate TD-DMRG algorithms and which evolution scheme is more suitable remain unclear.

To benchmark the performance of different TD-DMRG time evolution schemes, a proper model should be chosen. In this work, we focus on the vibronic coupling system represented by the Fenna-Matthews-Olson (FMO) complex from green sulfur bacteria. It is an extremely popular system for both experimental and theoretical studies of energy transfer in photosynthesis Engel et al. (2007); Cheng and Fleming (2009). The FMO model has become a “guinea pig” for comparing different computational methods Borrelli and Gelin (2017) and its exciton dynamics has been studied by many numerically exact methods, including quasi-adiabatic propagator path integral (QUAPI) Nalbach, Braun, and Thorwart (2011), hierarchical equations of motion (HEOM) Ishizaki and Fleming (2009); Shi et al. (2018), multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) Schulze and Kuhn (2015); Schulze et al. (2016), TD-DMRG Chin et al. (2013); Borrelli and Gelin (2017), etc.

In this paper, we select three time evolution schemes in TD-DMRG to simulate the exciton dynamics of 7-site FMO model, including global propagation and compression scheme with classical 4th order Runge-Kutta algorithm (P&C-RK4), TDVP with advanced matrix unfolding regularization scheme (TDVP-MU) and TDVP with projector-splitting algorithm (TDVP-PS). We firstly study the relative accuracy of the three schemes on the FMO model and then set out to explore the CPU-GPU heterogeneous computing to accelerate the TD-DMRG algorithms.

II Methodological Approaches

II.1 MPS Representation and TD-DMRG Algorithms

In the language of matrix product states (MPS) Schollwöck (2011), a quantum state $\ket{\Psi}$ under certain basis $\ket{\sigma_{1}\sigma_{2}\cdots\sigma_{N}}$, where $\ket{\sigma_{i}}$ is the basis for each degree of freedom (DOF) and is assumed to be orthonormal, can be represented as the product of a matrix chain, known as an MPS:

[TABLE]

$A^{\sigma_{i}}_{a_{i-1}a_{i}}$ are matrices in the chain connected by indices $a_{i}$. $\{\cdot\}$ in the summation represents the contraction of the respective connected indices, and $N$ is the total number of DOFs in the system. The dimension of $a_{i}$ is called (virtual) bond dimension denoted as $M_{\textrm{S}}$ or $|a_{i}|$, while the dimension of $\sigma_{i}$ is called physical bond dimension denoted as $d$. The graphical representation of an MPS is shown in Fig. 1(a). The many-body renormalized basis is defined as:

[TABLE]

where $[1:i]$ and $[j+1:N]$ represent that the renormalized bases are defined from the left and right side respectively. Hence, $A^{\sigma_{n}}_{a_{n-1}a_{n}}$ could be regarded as the coefficient matrix on the space spanned by $\ket{a_{n-1}[1:n-1]}\otimes\ket{\sigma_{n}}\otimes\ket{a_{n}[n+1:N]}$. However, in general, the renormalized basis is not necessary to be orthonormal, and the overlap matrix between them is defined as:

[TABLE]

The matrix product representation for a wavefunction is not unique in that inserting an identity matrix $I=GG^{-1}$ into the neighboring matrices $\cdots A^{\sigma_{i}}A^{\sigma_{i+1}}\cdots$ will obtain the same wavefunction but with different local matrices $\cdots A^{\prime\sigma_{i}}A^{\prime\sigma_{i+1}}\cdots=\cdots(A^{\sigma_{i}}G)(G^{-1}A^{\sigma_{i+1}})\cdots=\cdots A^{\sigma_{i}}A^{\sigma_{i+1}}\cdots$. Therefore, gauge conditions could be applied to eliminate the parameterization redundancy of an MPS. Among them, “mixed/left/right-canonical” gauge condition is usually adopted for convenience. A “mixed-canonical” MPS with gauge center at site $n$ is written as:

[TABLE]

where $L^{\sigma_{i}}_{l_{i-1}l_{i}}$ and $R^{\sigma_{j}}_{r_{j-1}r_{j}}$ satisfy:

[TABLE]

here $L^{\sigma_{i}}(R^{\sigma_{j}})^{*}$ represents the conjugate of $L^{\sigma_{i}}(R^{\sigma_{j}})$.

With Eq. 7 and  8, the renormalized basis fulfill the orthonormal relation $S[1:i]_{l_{i}l^{\prime}_{i}}=\delta_{l_{i}l^{\prime}_{i}}$ ($i=1,2,\cdots,n-1$) and $S[j+1:N]_{r_{j}r^{\prime}_{j}}=\delta_{r_{j}r^{\prime}_{j}}$ ($j=n,n+1,\cdots,N-1$). When the gauge center $n$ is at the right(left)-most site, the MPS is called left(right)-canonical state. With the ‘mixed/left/right-canonical” condition, the only redundancy left is that G is unitary. If we further decompose $C^{\sigma_{n}}$ by QR decomposition

[TABLE]

$L^{\sigma_{n}}$ fulfills the relation in Eq. 7 and $D_{l_{n}r_{n}}$ is also a coefficient matrix defined between site $n$ and $n+1$ on space $\ket{l_{n}[1:n]}\otimes\ket{r_{n}[n+1:N]}$. Afterwards, $D_{l_{n}r_{n}}$ is combined with $R^{\sigma_{n+1}}_{r_{n}r_{n+1}}$ to obtain $C^{\sigma_{n+1}}_{l_{n}r_{n+1}}=\sum_{r_{n}}D_{l_{n}r_{n}}R^{\sigma_{n+1}}_{r_{n}r_{n+1}}$ and apparently the gauge center is moved one site to the right. The reverse process to move the gauge center to the left could be carried out in a similar way as Eq. 9 while QR is replaced with RQ decomposition. More generally, a canonical MPS with gauge center at $n$ as Eq. 6 could be prepared by preforming QR decomposition from site 1 to $n-1$ sequentially and RQ decomposition from site $N$ to $n+1$ sequentially on any MPS, which is called canonicalisation. In the following paper, QR and RQ are not distinguished for simplicity.

If the QR decomposition in Eq. 9 is replaced with singular value decomposition (SVD),

[TABLE]

$\Lambda$ is a real diagonal matrix with elements ordered from large to small $\Lambda_{11}\geq\Lambda_{22}>\cdots\geq\Lambda_{kk}\geq 0$ ( $k=\min{[|\sigma_{n}|\cdot|l_{n-1}|,|r_{n}|]}$) and $\sum_{s}\Lambda^{2}_{ss}=1$ if $\ket{\Psi}$ is normalized. $U^{\sigma_{n}}$ fulfills the relation in Eq. 7. If we retain all $\Lambda_{ss}$ ($M_{\textrm{S}}=k$) and move the gauge center to the right with $C^{\sigma_{n+1}}_{l_{n}r_{n+1}}=\sum_{l_{n}}\Lambda_{l_{n}l_{n}}V^{\dagger}_{l_{n}r_{n}}R^{\sigma_{n+1}}_{r_{n}r_{n+1}}$, the wavefunction is not changed and SVD plays the same role as QR in canonicalisation but with a little bit higher cost. However, if we only retain the largest $M_{\textrm{S}}$ terms ($M_{\textrm{S}}<k$) to obtain $\tilde{C}^{\sigma_{n+1}}_{l_{n}r_{n+1}}$, $\tilde{C}^{\sigma_{n+1}}_{l_{n}r_{n+1}}$ could be a good approximation to $C^{\sigma_{n+1}}_{l_{n}r_{n+1}}$ with a smaller bond dimension. With this algorithm, a left(right)-canonical MPS $\ket{\Psi}$ could be “compressed” to $\ket{\tilde{\Psi}}$ by successive approximate SVD decomposition from site $N$ to 1 (site 1 to $N$). At each local step of the compression, two truncation criteria are commonly used

(i) a fixed pre-defined $M_{\textrm{S}}$;

(ii) adaptive $M_{\textrm{S}}$ with all $\Lambda_{ss}$ larger than pre-defined $\zeta$ retained.

Apart from canonicalisation and compression, several other operations are usually met within the MPS algorithms, including $\hat{O}\ket{\Psi}$ and $\ket{\Psi}+\ket{\Phi}$. Similar to MPS, a common quantum operator could be expressed as a matrix product operator (MPO) Schollwöck (2011); Chan et al. (2016):

[TABLE]

The graphical representation of an MPO is shown in Fig. 1(b). With the local matrix representation of wavefunction and operator in Eq. 1 and Eq. 11, $\hat{O}\ket{\Psi}$ (MPO$\times$MPS) could be calculated by contracting $W^{\sigma^{\prime}_{i},\sigma_{i}}A^{\sigma_{i}}$ locally:

[TABLE]

Suppose the bond dimensions of the original MPS and the MPO are $M_{\textrm{S}}$ and $M_{\textrm{O}}$ respectively, then the new state $\ket{\hat{O}\Psi}$ has bond dimension $M_{\textrm{O}}M_{\textrm{S}}$. $\ket{\Psi}+\ket{\Phi}$ (MPS + MPS) is constructed by merging the local matrices $[A^{\sigma_{i}},B^{\sigma_{i}}]$ block-diagonally,

[TABLE]

Suppose the bond dimensions of the original two MPSs are $M_{\textrm{S}_{A}}$ and $M_{\textrm{S}_{B}}$ respectively, then the new state $\ket{\Psi+\Phi}$ has bond dimension $M_{\textrm{S}_{A}}+M_{\textrm{S}_{B}}$.

In the next subsections, we present the ideas and the algorithms of three different TD-DMRG time evolution schemes adopted in our benchmark calculations afterwards: P&C-RK4, TDVP-MU and TDVP-PS.

II.1.1 Schemes based on Global Propagation and Compression

The propagation and compression (P&C) scheme is a global time evolution scheme benefiting from the modern MPS/MPO representation of DMRG wavefunctions and operators. If $\ket{\Psi(t)}$ and the time derivative $\ket{\dot{\Psi}(t)}$ are explicitly known (according to Schrödinger equation, the time derivative is $-i\hat{H}\Psi$ in atomic units), any ordinary differential equation (ODE) integrator for the initial value problem (IVP) such as the classical 4th order Runge-Kutta algorithm (RK4) we used in our former work Ren, Shuai, and Kin-Lic Chan (2018), could be applied to obtain the MPS of the next time step $\ket{\Psi(t+\tau)}$. The equations of RK4 are:

[TABLE]

Each $\ket{k_{i}}$ is represented by an MPS and $\hat{H}(t_{i})$ is represented by an MPO. Therefore, only two types of operations including MPO$\times$MPS (Eq. 13) and MPS$+$MPS (Eq. 14) exist in Eq. 16. After each operation, the new MPS such as $\ket{k_{1}}=-i\hat{H}(t)|\Psi(t)\rangle$ will have a larger virtual bond dimension, which should be compressed using algorithm based on Eq. 10 for the further operations. It is worth mentioning that as the new MPS is not canonical, canonicalisation should be carried out before the actual compression.

The procedure of P&C-RK4 described above involves $\hat{H}\ket{\Psi}$ (“MPO $\times$ MPS”) scaling at $\order{M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}}$ for each site (labeled as “$\hat{H}\ket{\Psi}$”), “MPS + MPS” which does not require intensive computation and MPS compression. The compression consists of four kinds of operations (two for canonicalisation in Eq. 9 and two for compression in Eq. 10), namely QR decomposition of tensors on each site (labeled as “QR”), matrix multiplication to absorb the decomposed coefficients (labeled as “MatMul-QR”) and the SVD counterparts (labeled as “SVD” and “MatMul-SVD”). “QR” and “MatMul-QR” both scale at $\order{M_{\textrm{S}}^{3}M_{\textrm{O}}^{3}d}$ while their SVD counterparts scale at $\order{\min(M_{\textrm{S}}^{3}M_{\textrm{O}}^{2}d,M_{\textrm{S}}^{3}M_{\textrm{O}}d^{2})}$ for “SVD” and $\order{M_{\textrm{S}}^{3}M_{\textrm{O}}^{2}d}$ for “MatMul-SVD”, because a compression reduces bond dimension from $M_{\textrm{O}}M_{\textrm{S}}$ to $M_{\textrm{S}}$ during the sweep. Thus, the bottleneck of P&C-RK4 is the canonicalisation and the scaling of the method in a single time step is $\order{NM_{\textrm{S}}^{3}M_{\textrm{O}}^{3}d}$.

Compared to the local evolution scheme such as TST Feiguin and White (2005), P&C could exactly calculate the operation of high order Hamiltonian such as $\hat{H}^{n}\ket{\Psi}$ and it is natural to increase the bond dimension with entanglement during the time evolution. However, though P&C is in principle a general evolution scheme for TD-DMRG, it prefers Hamiltonian whose MPO is straightforward to construct and has small bond dimension $M_{\textrm{O}}$, such as Frenkel-Holstein type model whose $M_{\textrm{O}}$ is in a linear relationship to the number of electronic states Ren, Shuai, and Kin-Lic Chan (2018). With regard to the time step size, the error of P&C does not have a monotonic relationship with time step and there is an optimal time step depending on the system as discussed in the previous studies Ronca et al. (2017); Ren, Shuai, and Kin-Lic Chan (2018). A smaller time step will improve the accuracy of RK4 integrator but lead to more MPS compressions and then deteriorate the whole accuracy.

II.1.2 Schemes based on Time Dependent Variational Principle

The Rayleigh-Ritz variational principle is widely used in finding an approximate ground state in time independent Schrödinger equation. Similarly, time dependent variational principle (TDVP) also provides a strong tool to find an optimal time dependent wavefunction if the wavefunction ansatz and the initial state are known. The Dirac-Frenkel TDVP is Dirac (1930); Gatti et al. (2017)

[TABLE]

In a geometric fashion, TDVP could be understood as an orthogonal projection of $-i\hat{H}|\Psi\rangle$ onto the tangent space of $\ket{\Psi(t)}$ at the current time:

[TABLE]

where $\hat{P}$ is the projector constructed by the orthonormal vectors in the tangent space. For a general MPS in Eq.1, the tangent space projector could be defined as:

[TABLE]

Where

[TABLE]

The inversion of the overlap matrix $S^{-1}$ accounts for the non-orthogonality of the renormalized basis and the “-” terms are to eliminate the parameterization redundancy Haegeman et al. (2011); Wouters et al. (2013). The graphical representation of the projector is shown in Fig. 2. In the literature, there are two different time evolution schemes based on TDVP. They differ in choosing the specific gauge condition of the MPS and in solving Eq. 18, which will be discussed in detail in the following.

TDVP-MU

In the first TDVP evolution scheme, the gauge freedom of MPS is fixed. For convenience, the projector in Eq. 19 could be transformed to Eq. 24 by combining the neighbouring “+” term and “-” term together except one “+” term with $i=n$:

[TABLE]

where:

[TABLE]

This type of projector or the corresponding tangent space vectors were firstly proposed by Haegeman et al. Haegeman et al. (2011) in TD-DMRG and restated by Wouters et al. Wouters et al. (2013) to derive several post-DMRG methods for ground state and excited states.

Eq. 25 and  26 could be further simplified by adopting a specific gauge condition, and then some overlap matrices turn to identity. Assuming that the MPS is left-canonical with gauge center at site $N$, $S[1:i]$ is reduced to $I$ and it is most convenient to set $n=N$ in Eq. 24. Inserting the simplified projector into Eq. 18 yields:

[TABLE]

[TABLE]

where:

[TABLE]

Similar equations can be derived for right/mixed-canonical MPS. With Eq. 28 and Eq. 33, it is straightforward to prove $\sum_{l_{i-1},\sigma_{i}}L^{\sigma_{i}*}_{l_{i-1}l_{i}}\partialderivative{t}L^{\sigma_{i}}_{l_{i-1}l^{\prime}_{i}}=0$ and then:

[TABLE]

which ensures that the left-canonical condition preserves during the time evolution formally. Though in practical numerical calculation with a finite time step, the relation in Eq. 34 would not be rigorously fulfilled, this problem is not severe with a proper time step in our experience. Otherwise, the more general equations in Eq. 24 which have already considered the non-orthogonality of the left and right renormalized basis should be used. And it could be proved that $\partialderivative{t}S[1:i]=0\,(i=1,2,\cdots,n-1)$ and $\partialderivative{t}S[j:N]=0\,(j=n+1,n+2,\cdots,N)$.

Eq. 27 and Eq. 28 together form a set of coupled nonlinear equations that are very similar to the standard EOMs of (ML-)MCTDH Meyer, Manthe, and Cederbaum (1990); Beck et al. (2000); Wang and Thoss (2003). To integrate these equations, we borrow ideas called variable mean field (VMF) and constant mean field (CMF) from the MCTDH community Beck and Meyer (1997); Beck et al. (2000). VMF employs an all-purpose integrator to directly solve the coupled equations. While in CMF, it is assumed that $H[i]$ and $S[i+1:N]$ generally change much slower in time than the local matrices $C^{\sigma_{N}}$ and $L^{\sigma_{i}}$. As a result, during the integration of Eq. 27 and Eq. 28 one may hold the “mean field” $H[i]$ and $S[i+1:N]$ constant for $\tau$ and evolve only the local matrix with time step smaller than $\tau$. Hence, CMF is more efficient yet less accurate than VMF. In this work, we use a second order CMF with a midpoint scheme in which $L^{\sigma_{i}}(0)$ and $C^{\sigma_{N}}(0)$ are integrated to $L^{\sigma_{i}}(\tau)$ and $C^{\sigma_{N}}(\tau)$ with “mean field” constructed by $L^{\sigma_{i}}(\tau/2)$ and $C^{\sigma_{N}}(\tau/2)$. In VMF, adaptive Dormand-Prince’s 5/4 Runge-Kutta method (RK45) is adopted to integrate Eq. 27 and Eq. 28, while in CMF, since Eq. 27 is linear, it is integrated by Krylov subspace method instead.

Another aspect should be concerned is that the inversion of $S$ would be unstable numerically if some eigenvalues of $S$ are very small. This problem will be severe when the state is weakly correlated (such as a Hartree product state which is usually an initial state) and $M_{\textrm{S}}$ is much larger than what is required. To some extent, this instability problem makes this evolution scheme paradoxical in that large $M_{\textrm{S}}$ should in principle push the result to a numerically exact limit but in fact deteriorates it. The same problem also arises in (ML)-MCTDH, where in order to make the EOMs more “well-behaved”, $S$ is usually replaced with a regularized overlap matrix $\tilde{S}$ Beck et al. (2000):

[TABLE]

Here $\varepsilon$ is a small scalar commonly from $10^{-8}$ to $10^{-14}$. More recently, an improved regularization scheme based on the matrix unfolding (MU) of the coefficient matrix by SVD in (ML-)MCTDH is proposed by Meyer and Wang, which has been proved to make the time integration more accurate and robust Meyer and Wang (2018); Wang and Meyer (2018). The same idea is adopted here to integrate Eq. 28, giving the name of the scheme “TDVP-MU”. When calculating the overlap matrix $S[i+1:N]^{-1}_{a^{\prime}_{i}a_{i}}$, the gauge center is moved to the $(i+1)$th site and the matrix at this site is further decomposed by SVD:

[TABLE]

where $|r_{i}|$ equals $|l_{i}|$. Thus, the overlap matrix $S[i+1:N]$ and its inversion could be expressed as $U^{*}\Lambda^{2}U^{T}$ and $U^{*}\Lambda^{-2}U^{T}$ respectively. The Hamiltonian matrix in Eq. 29 is also reconstructed. For site from $i+1$ to $N$ the matrix $A^{\sigma_{i}}$ in Eq. 32 is replaced with matrix $R^{\sigma_{i}}$ in Eq. 36 and then Eq. 28 with matrix unfolding algorithm becomes:

[TABLE]

The expression inside “$[\cdots]$” is $S[i+1:N]^{-1}$. The key point of this new regularization scheme is that the underlined part could be contracted first, which is $\delta_{r^{\prime}_{i}r_{i}}$. Thus, only the singular matrix $\Lambda_{r^{\prime}_{i}r^{\prime}_{i}}$ instead of $\Lambda_{r^{\prime}_{i}r^{\prime}_{i}}^{2}$ should be regularized:

[TABLE]

The power $1/2$ here is for consistency with the original regularization scheme (Eq. 35) and $\Lambda_{r^{\prime\prime}_{i}r^{\prime\prime}_{i}}$ in Eq. 37 is untouched in order to be minimally invasive as stated in Ref. Meyer and Wang, 2018; Wang and Meyer, 2018. In this work we choose $\varepsilon=10^{-10}$ unless otherwise stated.

The complexity analysis of the TDVP-MU scheme with VMF and CMF is as follows. Note that in TD-DMRG it is necessary to contract small matrices one by one instead of explicitly constructing the large tensors Chan (2004); Schollwöck (2011). VMF and CMF both require the calculation of environment matrices $h[1:i-1]$ and $h[i+1:N]$ scaling at $\order{M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d}$ (labeled as “Get Env”), SVD decomposition of the gauge center with the subsequent matrix multiplication scaling at $\order{M_{\textrm{S}}^{3}d}$ in Eq. 36 (labeled as “SVD” and “MatMul-SVD”), and the calculation of time derivative in Eq. 27 and Eq. 37 (labeled as “Deriv”). Evaluation of the time derivatives contains a number of matrix multiplications and the overall scaling is $\order{M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d+M_{\textrm{S}}^{3}d^{2}}$, where the first two terms and the last term correspond to the contraction of $H[i]$ and $(I-p[i])$ respectively. Thus the total scaling in VMF is $\order{N(M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d+M_{\textrm{S}}^{3}d^{2})}$, while in CMF the total scaling for a single step is $\order{Nf(M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d+M_{\textrm{S}}^{3}d^{2})}$ where $f$ is the number of steps required for the propagation of each local site within $\tau$. Although at first glance CMF is more time-consuming than VMF, much larger step size is possible for CMF which reduces the time cost spent on decomposing the matrix and constructing the environment, compensating the effect of the factor $f$. This time saving in CMF is important in MCTDH in that constructing the environment involves a contraction of the high order coefficient matrix  Beck and Meyer (1997); Beck et al. (2000); Gatti et al. (2017). However, in MPS context such advantage is not prominent because all matrices including the coefficient matrix $C^{\sigma_{N}}$ are of the same size. As a result, constructing the environment has a similar complexity as calculating the derivatives.

TDVP-PS

The second evolution scheme based on TDVP is called projector-splitting (PS). The idea of PS is that the tangent space projector in Eq. 19 is invariant under different gauge conditions. More specifically, after canonicalisation of a general MPS in Eq. 1 from site $N$ to $i+1$, $|r_{i}[i+1:N]\rangle$ becomes the right-hand orthonormal renormalized basis, which is related to $|a_{i}[i+1:N]\rangle$ by:

[TABLE]

The matrix $D$ is an upper triangular matrix in RQ decomposition after the canonicalisation described in Eq. 9. Therefore, the overlap matrix $S[i+1:N]$ equals $D^{*}D^{T}$ and the projector $\hat{P}[i+1:N]$ in Eq.21 defined for a general non-canonical MPS is transformed to:

[TABLE]

Similar result can be obtained for $\hat{P}[1:i]$:

[TABLE]

On the one hand, this definition of the tangent space projector does not contain any inversion operations of the overlap matrix, which seems to be a remarkable improvement over the first definition in Eq. 20 and Eq. 21. On the other hand, since the gauge is not fixed in different terms of this projector, the integration algorithm described in section TDVP-MU could not be directly applied. Lubich and Haegeman et al. proposed to use a symmetric second order Trotter decomposition to split the formal propagator into the individual terms Lubich, Oseledets, and Vandereycken (2015); Haegeman et al. (2016):

[TABLE]

Based on the propagator in Eq. 42, a single step of time evolution consists of a left-to-right sweep and a subsequent right-to-left sweep each with step size $\tau/2$. Taking left-to-right sweep as an example, the matrix at the gauge center $C_{l_{i-1}r_{i}}^{\sigma_{i}}$ is firstly evolved forward in time by applying the projector $\hat{P}[1:i-1]\otimes\hat{I}_{i}\otimes\hat{P}[i+1:N]$:

[TABLE]

where $H[i]$ and the ingredients $h[1:i-1]$, $h[i+1:N]$, $h[i]$ all have the same definitions as in Eq. 29 30, 31, 32 except that the $A^{\sigma_{i}}$ in Eq. 32 is replaced with $L^{\sigma_{i}}$ or $R^{\sigma_{i}}$ accordingly. Then, the evolved matrix $C_{l_{i-1}r_{i}}^{\sigma_{i}}$ is decomposed by QR as Eq. 9 to obtain the left-canonical matrix $L^{\sigma_{i}}_{l_{i-1}l_{i}}$ and the coefficient matrix $D_{l_{i}r_{i}}$. $D_{l_{i}r_{i}}$ is evolved backward in time by applying the projector $\hat{P}[1:i]\otimes\hat{P}[i+1:N]$:

[TABLE]

Afterwards, the gauge center is moved to site $i+1$ by contracting the evolved $D_{l_{i}r_{i}}$ and $R^{\sigma_{i+1}}_{r_{i}r_{i+1}}$ together to obtain $C_{l_{i}r_{i+1}}^{\sigma_{i+1}}=\sum_{r_{i}}D_{l_{i}r_{i}}R^{\sigma_{i+1}}_{r_{i}r_{i+1}}$. Following the procedure above, the sweep continues until all the individual projectors in Eq. 42 are applied. Inspired by the original two-site DMRG algorithm, it is also possible to formulate TDVP-PS into a two-site algorithm so that the bond dimension could grow up adaptively Haegeman et al. (2016). However, considering the two-site algorithm is much more expensive than the one-site algorithm both in the tensor contraction and QR decomposition, we use the one-site algorithm described above in this paper.

In principle, solving Eq. 43 and Eq. 44 could be accomplished by any ODE integrator such as the RK45 algorithm we use in TDVP-MU. However, since they are linear equations, the Krylov subspace method (Lanczos algorithm for Hermitian operator here) is preferred as it is unitary and is considered to be better than the explicit time-stepping integrators for matrix exponential operator Hochbruck and Lubich (1997); Haegeman et al. (2016). In our calculation, the dimension of the Krylov subspace is adaptive and the Lanczos iteration continues until $\ket{\Psi(t+\tau)}$ converges.

TDVP-PS has a similar computational complexity as TDVP-MU (CMF) except for the calculation of derivatives. With simpler EOMs in Eq. 43 and Eq. 44, the calculation of the derivative for forward evolution scales at $\order{M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d}$ (labeled as “Deriv-Forward”) and for backward evolution it scales at $\order{M_{\textrm{S}}^{3}M_{\textrm{O}}}$ (labeled as “Deriv-Backward”). So the total scaling for a single step is $\order{Nf(M_{\textrm{S}}^{2}M_{\textrm{O}}^{2}d^{2}+M_{\textrm{S}}^{3}M_{\textrm{O}}d)}$, where $f$ represents the number of Krylov subspace vectors. Besides calculating the derivatives, Krylov subspace method also includes other matrix operations of small size such as calculating the $\{\alpha\}$ and $\{\beta\}$ matrix elements with the Lanczos three-term recurrence relation Arnoldi (1951) and diagonalizing the tridiagonal Hamiltonian in the subspace (labeled as “Krylov”). Other labels for TDVP-PS is the same with that of TDVP-MU.

In summary, since the TDVP-MU scheme and the TDVP-PS scheme are both based on the time dependent variational principle, $\ket{\Psi(t)}$ should be the same if not considering the numerical error. Additionally, the two schemes both require to define a fixed $M_{\textrm{S}}$ a priori, and additional renormalized basis should be constructed smartly to complement the empty MPS space if the initial state is weakly correlated. Numerically, once $M_{\textrm{S}}$ is fixed, the error of these two TDVP-based schemes has a monotonic relationship with the time step size. The main difference lies in that TDVP-MU would introduce a minor artificial regularization, while TDVP-PS is inherently free of it. We also notice that TDVP-PS has already been implemented in MCTDH recently Lubich (2015); Kloss, Burghardt, and Lubich (2017); Bonfanti and Burghardt (2018).

II.2 CPU-GPU Heterogeneous Computing

In the TD-DMRG algorithms described above, the hotspots are usually the tensor contractions and the matrix decompositions. GPUs are well suitable for the former not the latter. Therefore, we adopt a CPU-GPU heterogeneous computing strategy. We store the matrices in the GPU memory instead of the host memory, and in most cases call cuBLAS to manipulate them. When doing matrix decomposition, we first transfer the matrix from the GPU memory to the host memory and then use a single CPU core to complete the decomposition, and in the end copy the matrix back. The overhead caused by data transfer is not negligible which we believe could be avoided by a smarter implementation because GPUs are able to run computation and data transfer at the same time.

For comparison, multi-core CPU calculation is also benchmarked. Here, we only consider the simplest strategy of CPU parallelism, which is breaking up the dense matrix computations into subblocks performed automatically by the standard linear algebra libraries such as Intel® Math Kernel Library in our case.

II.3 Computational Details of the FMO Model

The widely used Frenkel-Holstein Hamiltonian Holstein (1959); Schröter et al. (2015) to describe the 7-site FMO model is that:

[TABLE]

where $a_{m}^{\dagger}(a_{m})$ is the exciton creation (annihilation) operator on the $m$th site whose local excitation energy is $E_{m}$, $J_{mn}$ is the Coulomb interaction between the $m$th and $n$th site, $b_{m\lambda}^{\dagger}(b_{m\lambda})$ is the phonon creation (annihilation) operator of vibration mode $\lambda$ of $m$th site with vibration frequency $\omega_{m\lambda}$ and dimensionless electron-phonon coupling strength $g_{m\lambda}$. Fig. 3 shows a diagram of the model.

We obtain $E_{m}$ and $J_{mn}$ from previous literature Moix et al. (2011); Schulze et al. (2016) and each site has the same environment with $\omega_{m\lambda}$ and $g_{m\lambda}$ calculated by equally spaced discretization of the bath spectral density de Vega, Schollwöck, and Wolf (2015) from experiments Wendling et al. (2000); Schulze et al. (2016). For each site, 35 vibration modes are discretized, giving 252 DOFs in total (245 vibrational DOFs + 7 electronic DOFs). The dimension of basis for each mode varies with frequency and most modes have 8 or 4 phonon occupation levels. To minimize the entanglement, the electronic sites in the MPS chain are arranged in the order of [7, 5, 3, 1, 2, 4, 6]. All the local vibrational DOFs are arranged next to the local electronic site as our previous work Ren, Shuai, and Kin-Lic Chan (2018).

We simulate the exciton dynamics at zero temperature. At $t=0$, an excitation at electronic site 1 is prepared with all vibrations at their lowest energy level, which composes a Hartree product state. Therefore, for simulation with $M_{\textrm{S}}>1$, TDVP-based methods (TDVP-MU and TDVP-PS) should somehow find $M_{\textrm{S}}-1$ renormalized basis to complement the empty matrix elements in the initial MPS. To solve this, in all the simulations, we use P&C-RK4 to propagate the state with $M_{\textrm{max}}=256$ and time step size $\tau=80\,\textrm{a.u.}$ in the first 10 steps to $t=800\,\textrm{a.u.}$. Thereafter, the state is compressed to the pre-defined bond dimension $M_{\textrm{S}}$ and then evolves with different evolution schemes. Throughout this paper we use atomic unit (a.u.) as the unit of time unless otherwise stated and omit the unit of $t$ and $\tau$ for simplicity.

III Results and Discussions

III.1 Accuracy and Time Step Size

To quantitatively evaluate the relative accuracy in our calculations, the mean cumulative deviation of exciton populations at time $t$ is used:

[TABLE]

where $P^{(n)}$ represents the exciton population at the $n$th site of the 7-site FMO model. $P^{(n)}_{\textrm{ref}}$ is the reference result obtained from TD-DMRG calculation with an appropriate MPS bond dimension $M_{\textrm{S}}$ and evolution time step size $\tau$. The standard for choosing $M_{\textrm{S}}$ and $\tau$ is that the exciton populations have converged on both of them. Such convergence on $M_{\textrm{S}}$ for the TDVP-PS method is illustrated in Fig. 4 where all curves have already converged on $\tau$. The population obtained with $M_{\textrm{S}}=32$ has already been semi-quantitatively correct although the deviation is visible after $t>20000$. The curves for $M_{\textrm{S}}=64$ and $M_{\textrm{S}}=128$ almost overlap with each other. To be more discreet, we choose $M_{\textrm{S}}=256$ with $\tau=10$ as the reference for better accuracy. A more comprehensive and quantitative demonstration for the validity of our reference is included in Appendix A. Note that using reference calculated by TDVP-MU (VMF) with sufficiently small regularization parameter rather than TDVP-PS has little effect on results presented in this paper.

The mean cumulative deviation of the three schemes as a function of $t$ is plotted in Fig. 5 with $M_{\textrm{S}}=64$. For P&C-RK4 the time step size is 160 which strikes a balance between the RK4 integrator error and the MPS wavefunction compression error. The step size in TDVP-based schemes are chosen such that the error due to the integrator is negligible compared to the restricted $M_{\textrm{S}}$. For TDVP-MU (VMF) the step size controlled by the adaptive RK45 integrator (relative tolerance and absolute tolerance are $10^{-5}$ and $10^{-8}$ respectively) spans from 5 to 30. For TDVP-MU (CMF) and TDVP-PS we choose $\tau=2$ which is a relatively small value to reduce the error introduced by the CMF approximation or Trotter decomposition. We find that at long time scale the TDVP based methods exhibit almost the same error, as they share the same starting point during their derivation, while P&C-RK4 guided by a different philosophy shows larger error than TDVP based methods.

Although TDVP-MU and TDVP-PS are similar in long time limit, subtle differences between the two schemes when $t<15000$ are observed, which is illustrated in Fig. 6. TDVP-MU (CMF) is not shown in Fig 6 as it produces the same result with the VMF approach at small step size limit. Fig. 6 shows that TDVP-MU (VMF) is not as accurate as TDVP-PS at short time scale, which is probably caused by the artificial regularization in Eq. 38. Using TDVP-MU (VMF) with sufficiently small regularization parameter $\varepsilon$ as the reference will give the same tendency shown here. When the regularization parameter $\varepsilon$ is decreased from $10^{-7}$ to $10^{-16}$, which reduces the invasion to the original EOMs, the accuracy of TDVP-MU (VMF) is improved but a smaller time step size is required. However, it should be noted that the error caused by the regularization at the short time regime is several orders of magnitude smaller than the error at the long time regime controlled by the MPS approximation of the wavefunction with a restricted $M_{\textrm{S}}$ as illustrated in Fig 5. Thus in this case setting $\varepsilon$ to be $10^{-10}$ is more reasonable in terms of practical computation because this setup costs less simulation time than setting $\varepsilon$ to be $10^{-13}$ or $10^{-16}$ without significant precision deterioration. We also compare the new MU regularization scheme with the original regularization scheme and find that the new MU scheme is more accurate with the same $\varepsilon$ as the former studies in (ML-)MCTDH Meyer and Wang (2018); Wang and Meyer (2018).

Despite the fact that TDVP based methods will obtain the same results in the small step size $\tau$ limit and small regularization parameter $\varepsilon$ limit, they behave differently to the adjustment of $\tau$ because their EOMs and integration algorithms differ from each other. Since TDVP-MU (CMF) and TDVP-PS are both second order methods, we investigate how large the time step $\tau$ can be used in these two schemes with $M_{\textrm{S}}=64$. The upper panel of Fig. 7 illustrates the error at $t_{1}=6560$ which is at the early stage of the evolution. The error of TDVP-MU (CMF) shows a linear relationship with $\tau$ under logarithmic scale and its slope of linear fitting is 2.03, which is in accordance with the second order approximation of TDVP-MU (CMF) . It also indicates that for CMF the integrator error is dominant in this time scale. To the contrary, the error of TDVP-PS is insensitive to $\tau$ within the range of $\tau$ shown in Fig. 7, which means the error due to second order Trotter decomposition is smaller than the error due to MPS ansatz with $M_{\textrm{S}}=64$. In the lower panel, the time to measure the error is at $t_{2}=41120$ and we observe generally the same trend as in the upper panel, except that when $\tau\leq 20$ the error of TDVP-MU (CMF) and TDVP-PS are close to each other, indicating that the time step has already converged at this time point and the error is controlled by the MPS ansatz. However, as the upper panel, TDVP-PS allows a time step at least 16 times larger than TDVP-MU (CMF) for converged result. Therefore, though both TDVP-PS and TDVP-MU (CMF) implemented here are second order methods, the prefactor of the error term of TDVP-PS is much smaller than that of TDVP-MU (CMF).

III.2 GPU Acceleration

According to the analysis of computational complexity of the intensive steps in the three time evolution schemes in Section  II (a summary is shown in Table 1), we explore how these steps can be accelerated by multi-core CPU and particularly CPU-GPU heterogeneous computing. Although the time cost of each algorithm should vary with different implementations and there certainly are rooms for optimizations in our codes, we believe the data presented in this section reflects the correct tendency. Our benchmark platforms are Intel® Xeon® CPU E5-2680 v4 @ 2.40GHz for CPU-only calculations and Intel® Xeon® Gold 5115 CPU @ 2.40GHz with NVIDIA® Tesla® V100-PCIE-32GB for CPU-GPU heterogeneous calculations. Only one GPU is used in the benchmark.

For the three schemes discussed in this paper with $M_{\textrm{S}}=128$, the wall time cost of a single evolution step ($\tau=160$) as well as the composed sub-steps are presented in Fig. 8. For TDVP-MU (VMF), there are several steps within $\tau$. We firstly focus on the total time cost of each scheme in Fig. 8 from which we can learn that for all schemes the parallelization implemented in the standard linear algebra library on CPU is not optimal for TD-DMRG. Using 4 cores can merely double the computational speed, and further improvement is not significant up to using all 28 cores of the CPU. A more inspiring fact is that the CPU-GPU heterogeneous computing is able to boost the speed of time evolution up to 40 times for TDVP based algorithms, although for P&C-RK4 the effect of the GPU is comparable with 4 cores of CPU. We note that when doing calculations on GPU, the GPU usage is not always 100%, so more drastic speedup is expected with even larger bond dimension. Indeed, when $M_{\textrm{S}}=256$ the time costs of TDVP-MU (CMF) on single core CPU and the heterogeneous CPU-GPU are 4572 s and 63 s respectively, indicating a 73-fold acceleration.

A closer inspection of the sub-steps shown in Fig. 8 reveals more details about the acceleration. For P&C-RK4 , the bottlenecks of the algorithm are “QR” and “MatMul-QR” as discussed in Section II.1.1. Although “MatMul-QR” can be efficiently accelerated by multi-core CPU and CPU-GPU heterogeneous computing, the same is not true for “QR” and “SVD”, limiting the overall acceleration efficiency. When the calculation runs on the GPU, the matrix multiplication virtually costs no time while “QR” and “SVD” become the bottlenecks. When the GPU is incoperated in TDVP-MU (VMF) , the SVD decomposition for regularization becomes prominent even though it only scales at $\order{NM_{\textrm{S}}^{3}d}$ in a single step. Thus, it is suggested to return back to the original equation in Eq. 28 for higher efficiency once the regularization does not take effect anymore, owing to the fact that calculating the inverse of $S$ directly only costs $\order{NM_{\textrm{S}}^{3}}$. Even better acceleration efficiency can be achieved in TDVP-MU (CMF) because the SVD decomposition is carried out less frequently. The “Get Env” part also takes less time in CMF than VMF due to the same reason. The time costs of bottleneck sub-steps of TDVP-PS are similar with that of TDVP-MU (CMF) except that the integrator seems to take a large fraction of time especially when running on the GPU. This is considered to be a demonstration of the GPU latency which makes the total time cost of lots of operations on small matrices larger than the total time cost of a few operations on large matrices. When $M_{\textrm{S}}$ is increased to 256 the relative time cost of the “Krylov” part is diminished. The exact abnormality with GPU acceleration is also manifested to a lesser extent in TDVP-MU (CMF) where the time-consuming “other” part is primarily composed of operations in the Runge-Kutta integrator. Hence, if a smaller bond dimension is employed, devoting more resource on the computation will not gain much benefit. In fact, for TDVP based schemes the absolute wall time with GPU are roughly the same between $M_{\textrm{S}}=128$ and $M_{\textrm{S}}=32$. For example, the time costs for a single step of evolution in TDVP-PS with $M_{\textrm{S}}=128$ and $M_{\textrm{S}}=32$ are 10.8 s and 7.9 s respectively. This phenomenon implies that when $M_{\textrm{S}}\leq 128$ the latency of GPU calculation is of the same order with the actual computational time cost. Therefore, with the CPU-GPU heterogeneous computation, TD-DMRG simulation with $M_{\textrm{S}}>100$ for a large-scale system would become a routine task in the future.

IV Conclusion

To summarize, in this paper, we carry out numerical benchmarks on three different TD-DMRG time evolution schemes, which are P&C-RK4, TDVP-MU and TDVP-PS, in terms of accuracy based on the vibronic coupling model represented by the FMO complex. After defining the tangent space projector for a general non-canonical MPS, the EOMs of TDVP-MU and TDVP-PS are re-derived from the same starting point but with different gauge conditions. The numerical results demonstrate that TDVP-MU and TDVP-PS could indeed obtain similar accuracy with a converged time step size $\tau$ , while P&C-RK4 is not as accurate as them. However, surprisingly, the converged time step size of TDVP-PS is at least 16 times larger than that of TDVP-MU (CMF), although both of them implemented in this work are second order methods. Regarding computational efficiency, we firstly analyze the complexity of each evolution scheme. To further accelerate the intensive tensor computations in TD-DMRG, we adopt a CPU-GPU heterogeneous computing strategy, in which CPU and GPU are respectively responsible for the tensor decomposition and contraction. We find that for the TDVP-based schemes where the tensor contraction is the main bottleneck, this heterogeneous computing approach is able to speed up these two schemes by up to 73 times ($M_{\textrm{S}}=256$). Finally, we conclude that with the accurate and efficient TDVP-based evolution schemes (especially TDVP-PS) and CPU-GPU heterogeneous algorithm/hardware, TD-DMRG has been a promising method for the large-scale quantum dynamics simulation of real chemical and physical problems.

Acknowledgements.

This work is supported by the National Natural Science Foundation of China through the project “Science CEnter for Luminescence from Molecular Aggregates (SCELMA)” Grant Number 21788102, as well as by the Ministry of Science and Technology of China through the National Key R&D Plan Grant Number 2017YFA0204501. J.R. is also supported by the Shuimu Tsinghua Scholar Program. The authors are indebted to Prof. Garnet Chan for the stimulating discussions. The authors also gratefully thank Mr. Hongde Yu for support on crystal structure of the FMO complex.

Appendix A The validity of the reference state

Here we present the data for the numerical convergence of the reference obtained by TDVP-PS with $M_{\textrm{S}}=256$ and $\tau=10$. The mean cumulative deviation with various $M_{\textrm{S}}$ and $\tau=10$ at $t_{1}=6560$ and $t_{2}=41120$ is shown in Table 2, while the same data for $M_{\textrm{S}}=256$ and various $\tau$ is shown in Table 3.

From Table 2 and Table 3 we can see that the reference is accurate to a level of approximately $10^{-11}$ for short time scale ($t_{1}=6560$) and $10^{-6}$ for long time scale ($t_{2}=41120$), which is considered to be accurate enough for the error scales discussed in the main text.