Stochastic Adaptive Single-Site Time-Dependent Variational Principle.

In recent years, the time-dependent variational principle (TDVP) method based on the matrix product state (MPS) wave function formulation has shown its great power in performing large-scale quantum dynamics simulations for realistic chemical systems with strong electron-vibration interactions. In this work, we propose a stochastic adaptive single-site TDVP (SA-1TDVP) scheme to evolve the bond-dimension adaptively, which can integrate the traditional advantages of both the high efficiency of the single-site TDVP (1TDVP) variant and the high accuracy of the two-site TDVP (2TDVP) variant. Based on the assumption that the level statistics of entanglement Hamiltonians, which originate from the reduced density matrices of the MPS method, follows a Poisson or Wigner distribution, as generically predicted by random-matrix theory, additional random singular values are generated to expand the bond-dimension automatically. Tests on simulating the vibrationally resolved quantum dynamics and absorption spectra in the pyrazine molecule and perylene bisimide (PBI) J-aggregate trimer as well as a spin-1/2 Heisenberg chain show that it can be automatic and as accurate as 2TDVP but reduce the computational time remarkably.

Based on a linear-chain matrix product state (MPS) representation for the many-body wave function with the advantages of high compression and local structure, the time-dependent variational principle (TDVP) approach[1,2] as a time evolution method for MPS (for other approaches, see refs (3−9)) has been shown to be a powerful tool for simulating the quantum dynamics and spectroscopy of large realistic chemical systems with electron–vibration (electron–phonon) interactions.[10−19]

In MPS-TDVP simulations, there are several types of errors, namely the projection error of projecting Ĥ|ψ⟩ onto the tangent space of the given MPS |ψ⟩, the truncation error in tensor singular value decomposition (SVD), and the Krylov and time step error for the evolution operation of on |ψ⟩.[3] The most widely used single-site TDVP (1TDVP) is symplectic and accordingly has zero truncation error, but the bond-dimension m of the time-evolving MPS via 1TDVP cannot increase to accommodate entanglement increasing with time. In many cases, numerous test calculations need to be performed for a given task to test the convergence with respect to the bond-dimension m and obtain a suitable value of m for 1TDVP simulation. However, it may be quite tedious and expensive, especially for realistic complicated systems. To overcome this disadvantage, one can resort to two-site TDVP (2TDVP). Compared to 1TDVP, 2TDVP can adaptively increase the bond-dimension m and has a smaller projection error but a nonzero truncation error. Generally, 2TDVP will be more accurate and robust but less efficient than 1TDVP. Therefore, the current TDVP calculations of strongly correlated large systems always face a dilemma: it is difficult to balance the accuracy and robustness of 2TDVP and the computational efficiency of 1TDVP, especially for simulations of vibronic problems, in which the size of local Hilbert spaces can be very large.

To break the above bottleneck, Yang and White[20] recently suggested a method to increase the bond-dimension of 1TDVP dynamically via temporarily creating an MPS representation of the current time-evolved state by expanding the MPS to represent both the current state and a sequence of Krylov vectors generated from it. Utilizing the benefit of the expanded manifold coming from this representation, the subsequent 1TDVP time step becomes more accurate and reliable. In addition, Dunnett and Chin[21] also proposed an autonomously adaptive variant of 1TDVP to capture the growing entanglement “on the fly”. To ensure that the local bond-dimension m can evolve, they utilized full QR factorization instead of the normally used thin or reduced QR factorization in updating the MPS local tensors. To measure the convergence behavior with increasing bond-dimension m, the norms of the effective Hamiltonian applied to their respective MPS site tensors are evaluated.

In this Letter, we propose a stochastic adaptive 1TDVP method to evolve the bond-dimensions automatically, which is computationally efficient without any additional tensor operations. This method is based on the statistical distribution law of tensor SVD singular values λ (i.e., the square root of eigenvalues of the left or right subsystem’s reduced density matrix) in the MPS decomposition. It is well-known that the λ spectrum shows a clear exponential decay as for the nth singular value in generic gapped one-dimensional systems, implying a linear relationship between the logarithmic singular values (s ≡ log(λ)) and .[22] However, quantitative deviations from this relationship have been observed for large m,[23,24] preventing the use of this property to estimate the necessary m value for a predetermined sufficiently small truncation threshold ε. There is, however, more information available: one can derive a so-called entanglement Hamiltonian[25−30] from the reduced density operator ρ of an MPS as Ĥen = −log ρ. Its eigenvalues (energy levels) are, up to a constant factor of no further importance, the s ≡ log(λ); their level spacings (or first-order differentials) are then given by Δs = s – s, and their second-order differentials are given by Δ2s = Δs – Δs. We assume here that, in general, entanglement Hamiltonians exhibit the same level spacing statistics as many-body Hamiltonians do. There, one finds that the level spacings are distributed according to a Poisson distribution as arises for randomly distributed eigenvalues in the case of integrable systems and according to a Wigner distribution in the case of nonintegrable systems, as a result of random-matrix theory.[31−33] This implies certain distributions (exponential and quasi-Gaussian) for the second-order differentials Δ2s = Δs – Δs (see the Supporting Information), whose parameters are fitted from the spectrum of the reduced density operators. Using these distributions, we can estimate the small singular values in TDVP. The tests on simulating the vibrationally resolved quantum dynamics and absorption spectra in the pyrazine molecule and perylene bisimide (PBI) J-aggregate trimer as well as spin-1/2 Heisenberg chain show that this can be automatized and be as accurate as 2TDVP but save a lot of computational time.

To analyze the singular value distribution behavior in real-time TDVP simulations, we take a 4-mode exciton–vibration model for the pyrazine molecule[34] as an exampleHere, the three terms describe the electronic and vibrational parts as well as their interactions, respectively. ↠and â represent the creation and annihilation operators of electronic state i, while ε denotes the onsite energy of state i (i = j) or the electronic coupling between state i and state j (i ≠ j). b̂† and b̂ are the raising and lowering operators for vibration mode K, while ω is the vibrational frequency of mode K. g is the coupling between vibrational mode K and the electronic Hamiltonian term. Due to the small number of the degrees of freedom in this model, we can perform an exact 2TDVP simulation (total time 4000 au, the time step is 20 au) without any SVD truncation errors to avoid discarding small singular values. The computational details can be found in the Supporting Information. Then, we extract the singular values {λ} on each bond at different time steps and show the logarithms of singular values s versus for two example SVDs with large bond-dimensions within short/long time ranges across multiple entanglement regimes in Figure a,b.

Figure 1

Singular value distribution versus log2n at (a) bond 4–5 at a time point of 300 au and (b) bond 6–7 at a time point of 2000 au. The red and green lines are fittings by with α = 1 for the largest 10 singular values and all singular values, respectively. Distributions of all the second-order differentials (c) at bond 4–5 during time range of 0–300 au and (d) at bond 6–7 during time range of 0–2000 au with exponential and quasi-Gaussian fittings.

There are obvious linear correlations for large s, i.e., a small n, but the linear relationship breaks down for a smaller s (i.e., large n). This agrees with previous findings in grand canonical ensembles and Ising models[23,24] and can be ascribed to the increasing degree of energetic degeneracy of higher excited states. However, interestingly, we find that s’s second-order differentials (Δ2s) fulfill an exponential (eq ) or quasi-Gaussian (eq ) distribution, as shown in Figure c,d.

This is related to random features in the vanishingly small singular values and the analogy between the occurrence of s’s first-order difference and inhomogeneous Poisson process or Wigner distribution (see the Supporting Information for detailed discussion). It should be also noted that the exponential/quasi-Gaussian distribution behavior of the second-order differentials of logarithms of singular values may be generic for other model Hamiltonians, such as the spin-1/2 Heisenberg model and Fermi–Hubbard model, as we show in Supplementary Figures S3 and S5.

Now, we utilize the above exponential/quasi-Gaussian distribution behavior to implement the bond-dimension evolution in our stochastic adaptive 1TDVP (SA-1TDVP) method. First, we recall the MPS wave functionwhere σ represents the local d-dimensional basis on site i and N is the total number of sites. M(Mα) are rank-3 tensors, where α represents the MPS bond-dimension index linking site i and site i + 1. To optimize the MPS wave function with a finite value of m for the maximal bond-dimension, one needs to update M tensors successively. In the ith step of a 1TDVP left-sweep, only the tensor M on single site i is updated during the time evolution. In a standard 1TDVP, a reduced SVD is then performed for the (dm × m)-dimensional M[ (M)Here, the (dm × m)-dimensional matrix U has orthogonal columns, S contains the m singular values, and square (m × m)-dimensional V† has orthogonal rows. To further run the time evolution at site i + 1 and optimize Mσ in the right-direction sweeping process, M[ is then replaced by U, and SV† is contracted into M[. Because the number of nonzero singular values in S is always no more than m, the bond-dimension of M will not exceed m. In order to increase the bond-dimension to describe the growing entanglement adequately, we extend the reduced SVD in eq towhere the block Uex is constituted from Δm new column vectors. Gram–Schmidt orthogonalization is performed for the new randomly generated vectors in Uex with respect to the matrix U from reduced SVD, while the whole matrix is still column orthogonal but with a shape of (dm × (m + Δm)). In the language of MPS, the columns in Uex represent the newly added bases in (left) subspace. The subspace expansion is realized by random generation and the subsequent Gram–Schmidt orthogonalization. O is a block filled with 0. Note that if the sweep direction is reversed, one needs to extend the matrix V† instead.

To determine the column size Δm of Uex, we adopt a stochastic growth algorithm. The detailed SA-1TDVP algorithm for updating the bond-dimension of site i at an evolved time step then works as

Try 2TDVP for the initial T (∼20–50) time steps and save all the second-order difference values existing on each bond. Use the least-square method (using mean-square error as the cost function) to fit their distributions by exponential/quasi-Gaussian lines (eqs and 3, respectively) to get the parameter value for β. Here, we assume that β does not have a significant change during time evolution after the few initial steps, which is supported by Supplementary Figure S7. Then, the fitting will be performed only once, and the time cost of fitting is negligible compared to the TDVP calculation.

Perform a reduced SVD on (dm × m)-dimensional site tensor M[ and obtain nonzero singular values {λ} and matrices U, V†. Set N = m.

Generate a new random number Δ2s satisfying the above exponential/quasi-Gaussian distributions with the parameter β obtained from fitting and then calculate λ from λ, λ, and Δ2s.

If the new estimated singular value λ > λ or λ < ε, exit the stochastic iteration process and go to step 5; otherwise add λ to {λ} and set N ← N + 1; then, go back to step 3.

Set mnew = N and build mnew – m new columns of U (for left-sweep) or N rows of V† (for right-sweep) to update M[.

Our approach for growing m adaptively avoids the additional and computationally costly operation of applying the Hamiltonian operator to the MPS, whose complexity is O(m3wd + m2w2d),[3] where w represents the bond-dimension of the matrix product operator (MPO). Moreover, as shown in eqs and 6, the extension of the basis will not change the MPS itself, which is crucial for time evolution and is not fulfilled by those perturbative methods suggested for static single-site density-matrix renormalization-group (DMRG) calculations.[35,36] We now first test the SA-1TDVP method by using the 4-mode pyrazine model. We compare the performances of convergence and accuracy of conventional TDVP and our method in Figure . The autocorrelation function C(t) is defined asHere, |ψ(t)⟩ is the MPS at time t. Because the number of degrees of freedom in this model is small, we can perform an exact 2TDVP simulation with no truncation error, and we take this as a reference. It is shown that one can obtain very accurate results using SA-1TDVP with an SVD truncation cutoff ε of 10–8. The two different distribution formulas (exponential and quasi-Gaussian) give very similar performances. Figure b demonstrates that the error of time correlation functions for SA-1TDVP is of similar order of magnitude as the cutoff in the entire time evolution of this model. One may notice the slight increase of SA-1TDVP’s errors at long-time limit. This inefficiency may be caused by the fact the added random vectors in SA-1TDVP are nonoptimized and will be more likely orthogonal to the physically relevant renormalized states from 2TDVP when m is large. Figure c further shows that with the same cutoff value SA-1TDVP costs 90% less time than 2TDVP. At the same time, as a stochastic algorithm, SA-1TDVP is found to generate nearly the same bond-dimension as 2TDVP as shown in Figure d. In this particular case, the speedup of the two distributions is essentially the same, but this is a peculiarity. The efficiency gain for SA-1TDVP over 2TDVP was also found for the tests of the Heisenberg model (see Supplementary Figure S4), where the quasi-Gaussian distribution shows a slightly better performance than the exponential one, saving 84 and 75% computational time when compared with 2TDVP respectively.

Figure 2

Results for the 4-mode pyrazine model from SA-1TDVP (T = 50) and conventional 2TDVP methods. (a) The absolute correlation function from conventional 2TDVP and SA-1TDVP. (b) The error of correlation functions, measured by the absolute value of differences between them with those by exact 2TDVP. (c) The time cost of the SA-1TDVP and 2TDVP. (d) The increase of the max bond-dimension of SA-1TDVP and 2TDVP.

We now switch to the 24-mode pyrazine model[37] to show the performance of SA-1TDVP in larger systems. MPSs with different bond-dimensions in conventional 1TDVP are obtained from a 1 au preliminary 2TDVP calculation. For SA-1TDVP, all dynamics starts from an m = 10 MPS. We find that one can obtain nearly converged results using SA-1TDVP with a decreasing SVD truncation cutoff ε, as shown in Figure a. From Figure b, one can find that in the long-time range, the error of SA-1TDVP does not increase much because of the adaption of bond-dimensions in the latter, as illustrated in Figure d. From Figure c,d we find that the small change of cutoff may cause different behaviors. This is because a finite cutoff and relatively large bond-dimension result in a dense sequence so that one can find many singular values in a small interval. Other computational details can be found in the Supporting Information.

Figure 3

Results for the 24-mode pyrazine model from SA-1TDVP with different cutoffs. (a) The absolute autocorrelation function C(t) from SA-1TDVP. (b) The error of correlation versus reference SA-1TDVP (ε = 1 × 10–7). (c) The time cost of the SA-1TDVP. (d) The increase of the max bond-dimension of SA-1TDVP. All calculations use the exponential distribution here.

Finally, we test the calculation of the absorption spectra of 24-mode pyrazine and PBI trimer models[38,39] and compare our SA-1TDVP results with other accurate 2TDVP and multiconfiguration time-dependent Hartree (MCTDH)[40] calculations in Figure . The SA-1TDVP results are in good agreement with accurate methods. Our method reproduces the details successfully, such as the small oscillations in Figure a and the relative magnitudes of peaks in Figure b. This implies that SA-1TDVP provides an efficient and automatic computational tool for the full quantum dynamics simulation of realistic chemical problems.

Figure 4

Simulated absorption spectra for (a) the 24-mode pyrazine model and (b) PBI trimer model. Results of SA-1TDVP are compared with 2TDVP (ε = 10–8)[15] and MCTDH[38] reference results.

In conclusion, we observe an exponential/quasi-Gaussian distribution behavior for the second-order differentials of logarithmic singular values in DMRG and TDVP simulations for the first time, which can be ascribed to the Poisson or Wigner distribution of Δs. We propose an SA-1TDVP algorithm based on this distribution behavior to perform subspace expansion of MPS randomly and adaptively during the time evolution. We compare the accuracy and efficiency of our method in Heisenberg and exciton–vibration interaction models with conventional 2TDVP and 1TDVP as well as MCTDH. The tests on simulating the vibrationally resolved quantum dynamics and absorption spectra in the pyrazine molecule and PBI trimer show that it can be automatic and as accurate as 2TDVP but increase the computational efficiency significantly. The automatic adaption means that SA-1TDVP has the advantage of being able to avoid tedious tests of trying different bond-dimensions repeatedly in conventional 1TDVP. It is also worth to mention that this strategy can be easily incorporated in any time-dependent DMRG (tDMRG) code if there are SVD steps. For these reasons, SA-1TDVP provides an efficient, accurate, and user-friendly tool for quantum dynamics simulations of large strongly correlated systems.