In a recent paper [H.-D. Meyer and H. Wang, J. Chem. Phys. 148, 124105 (2018)], we have examined the regularization of the equations of motion (EOMs) of the multiconfiguration time-dependent Hartree (MCTDH) approach. We could show that the standard regularization scheme used by almost all researchers in the field is not optimal. The improved regularization allows for larger values of the regularization parameter ϵ, is less sensitive to the actual choice of ϵ, and performs the rotation of initially unoccupied single-particle functions into the “correct” direction in Hilbert space much faster than the old scheme. The latter point increases both the accuracy and efficiency of time propagation for challenging problems. For simple problems, the new scheme requires some additional numerical work as compared with the old scheme, ranging from negligible to almost doubling the total numerical labor. For demanding problems, on the other hand, the additional numerical work of the new scheme is often overcompensated by less steps taken by the integrator. In the present paper, we generalize the new regularization scheme to the multi-layer (ML) extension of MCTDH. Although the principle idea of the new regularization scheme remains unaltered, it was not obvious how the new scheme should be implemented into ML-MCTDH. The ML-MCTDH EOMs are much more complicated than the MCTDH ones, and for optimal numerical performance it was necessary to derive a recursive algorithm for implementing the new regularization scheme.
I. INTRODUCTION
Over the last two decades, the multiconfiguration time-dependent Hartree (MCTDH) approach1–6 has been very successful. Several problems of molecular quantum dynamics can be solved in full dimensionality only with this method, for example, the computation of the photoexcitation spectrum of pyrazine,7 or other photoexcitation or photoionisation spectra of polyatomic molecules with vibronically coupled electronic states,8–10 the computation of the IR spectrum of the Zundel cation (),11–14 and studies on the tunneling splitting in malonaldehyde.15–18 Furthermore, difficult reactive scattering problems have been successfully tackled with MCTDH; see, e.g., the computation of reaction rates and cross sections for the H + CH4 → H2 + CH3 reaction.19–22 Moreover, studies on electron transfer reactions in condensed phases23–26 have been performed by employing high-dimensional models for the solid.
In 2003, Wang and Thoss introduced a multi-layer (ML) extension to the MCTDH algorithm,27 and in 2008 Manthe provided a recursive formulation28 of the rather intricate equations; the latter formulation simplifies the implementation of the cumbersome ML-equations. The ML-MCTDH method27–30 has been shown to be very successful when treating systems with many degrees of freedom (DOFs)—between 10 and 1000 DOFs, say—and when many of these DOFs are only weakly coupled. Examples of successful applications of ML-MCTDH are, e.g., the study of photo-dissociation of a molecule embedded in a complex,31 vibronic dynamics,32 reactive scattering,33,34 exciton vibrational dynamics,35–37 singlet fission,38 surface scattering,39 quantum transport,40–45 and charge transfer processes.46–59
Despite the great success of MCTDH and ML-MCTDH, there remains a problem with their equations of motion (EOMs). If there are unoccupied single particle functions (SPFs)—and starting from a Hartree product all SPFs, expect the first one in each degree of freedom, are unoccupied—then there are two problems. First, the initially unoccupied SPFs are undefined and they can be chosen at random. Second, the equations of motion (EOMs) are singular and need to be regularized to make them numerically solvable. The first problem was addressed by Manthe60,61 (see also the discussion of Lee and Fischer62 on bosonic systems), who showed how to define optimal initially unoccupied SPFs, where the optimal choice increases the accuracy of the MCTDH propagation.61 The second problem is usually solved by regularizing the density matrix, as the singularities of the density matrix cause the singularities of the EOM. However, we recently have shown63 that much better results are obtained when the A-coefficients (the expansion coefficients) are regularized, which in turn regularizes the density. The new regularization scheme quickly rotates the unoccupied SPFs toward their optimal form and thus lessens the need for starting with optimal SPFs. This could be demonstrated for MCTDH calculations on a series of examples. In the present paper, we extend our approach to ML-MCTDH. For ML-MCTDH calculations, the new regularization scheme is in general more important than for MCTDH ones, as in ML-MCTDH singularities are more frequent, because they may appear at every layer. As ML-MCTDH is a rather intricate formalism, this extension of the regularization is not simply straightforward, although the general idea of the new regularization scheme is identical to the one of our previous work.63
II. THEORY
A. Expansion of an ML-MCTDH wavefunction
We make use of the nomenclature introduced by Manthe,28 with slight modifications introduced in Ref. 29, and write an MCTDH wavefunction as
where are multi-dimensional mode-combined logical coordinates. There are p1 logical coordinates and denotes the number of SPFs of the κth mode. The upper index 1 indicates that the quantities are for the layer number 1, also called the top layer. In the MCTDH method, which is equivalent to a two-layer ML-MCTDH approach,64 the multi-dimensional SPFs are given by the second layer equation
where χ denotes a primitive basis function and is here the physical coordinate of the κ2th DOF of the node (2; κ1), which has branches.
Adding additional layers, one has to specify which node of the ML-tree is under consideration. For this, it is convenient to introduce the symbol
and similar for z + 1, etc. Here l denotes the layer and κ1, …, κl−1 denote the modes of the layers above which leads to the node under consideration. Hence, z precisely defines one particular node in the ML-tree.28,29
The ML expansion can now be written in full generality,
which implicitly defines the configurations and the composite index J, and where the coordinates are combined as
Equation (4) holds for the top layer if one defines q0 as the combination of all DOFs and sets , and it also holds for the bottom layer if one identifies with the primitive function , when z points to a bottom node. Hence Eq. (4) holds for all layers.
The EOMs for the propagation of the SPFs are formally the same for all layers
where is the projector onto the space spanned by the SPFs, is a density matrix, and is a matrix of mean-field operators acting on the functions. Note that this equation is identical to the MCTDH-EOM, except that all quantities carry in addition the node-index z.
The EOMs for the A-vectors follow from Eq. (6) by projecting the SPFs onto their corresponding basis,
There is no contribution from because, due to the projector, the SPFs are orthogonal to their time derivatives. However, the top A-vector obeys a different equation
B. Recursive formulation of the ML equations of motion
The density and the mean fields are given by
and
and a single-hole function (SHF), , is defined by erasing a SPF from the full wavefunction,
The top-layer SHF has a simple appearance
where is a composite index similar to J, but it misses the κ1th entry. Similarly, denotes a product of SPFs except the function , i.e., in general, . And is similar to but with the κ1th index replaced by m. (Compare with Ref. 29.) The SHFs of the deeper layers are more complicated and they obey the recursion28
where
The second layer SHF thus reads
It is convenient to abbreviate the product of A-vectors by B, starting with the top layer
and continuing recursively
where is a composite index, which can be extended recursively as .
We consider as a matrix and, due to the orthonormality of the SPFs, write the densities as
To derive a recursive generation of the mean-fields, we first introduce a sum-of-products (SOP) form of the Hamiltonian
where operates on the κth DOF of the bottom layer. The operator summands in Eq. (20) can be written as
where acts on the logical coordinate , while acts on the corresponding remaining space. As seen from Eq. (21b), can as well be written as a direct product of operators , which act on logical coordinates . In fact, both operators and can be broken down to a simple product of operators which act on one primitive coordinate only [see Eq. (20)].
Using the SOP form, we write the mean field operator as
where the mean-field matrix is given by
with
For an efficient implementation of the EOMs, the matrix elements of the operators with respect to the SPFs are needed because the SPFs of one layer (z) are expanded in the SPFs of the layer below (z + 1). We follow the discussion of Ref. 29 and define
Here the last line defines the recursion of the h-matrix elements, starting from the bottom layer
where lmax denotes the maximal layer (i.e., lmax = 2 for normal MCTDH) and proceeds bottom-up through the tree structure. In a practical implementation of the EOMs (6), one replaces with
because an SPF is expanded into the SPFs .
From Eq. (21b) follows the top-down recursion for ,
and inserting this into Eq. (24) gives the recursion for M,
where jν and are the indices which appear in and , respectively. With this result, we finally arrive at the desired recursion for the mean-field matrix [compare with Eqs. (17) and (23)],
where the symbol denotes a composite index similar to J, but with the κlth and νth index missing. And is to be interpreted as , but with the κlth and νth entry of J replaced with n and jν.
In an actual implementation, only the A-tensors of the SPFs are propagated. Using Eqs. (4), (6), (7), (22), and (27), the EOMs are given by
where the indices and are the ones which appear in the composite indices K and I, respectively.
C. Regularization of the ML-EOM
1. Formal expression
To regularize this EOM, one may proceed similarly as done in Ref. 63 and perform a singular value decomposition (SVD) of the matrix B,
where V is a unitary matrix, λ is a diagonal matrix with real non-negative singular values, and U is a matrix with orthonormal column vectors. Here
is the size of the index space {}.
which shows that the singular values are the square roots of the eigenvalues of ρ, i.e., of the so-called natural populations.
Proceeding as discussed in Ref. 63, we insert the SVD (33) into Eq. (32) (except for the very last matrix B), let several terms cancel each other, and obtain
To arrive at a singularity-free EOM, one has to regularize the singular values λ. For this, we may replace with [or use a variant of replacing with ], where ϵ denotes a small positive regularization parameter.63 The conventional approach, where only the density matrix is regularized, is equivalent63 to replacing in Eq. (36) with . The jth diagonal element of the latter product vanishes if the jth singular value vanishes. Adopting the new regularization scheme, however, the jth element of assumes in this case the large value ϵ−1/2. If there is an unoccupied SPF (more precisely an unoccupied natural orbital), then the old scheme will not move this SPF until it has obtained an occupation of the order of ϵ. The new scheme, on the other hand, will immediately start to rotate the unoccupied SPF with a high speed into its “optimal” direction in Hilbert space.63
However, the procedure (33) and (36) is not practical for two reasons. First, , the size of the index space of , may become very large, in particular, for the lower layers. An SVD of the matrix B in general will become very expensive, if not impossible. Second, as discussed in Sec. II B, the matrices B and M are not explicitly built in practice, but the densities and the mean-field matrices are built recursively starting from their simple top-layer forms.
2. Practical implementation
To derive a practical recursive evaluation of the SVD, we insert the diagonal form of the density, Eq. (35), into the recursion formula for the density, Eq. (19),
where in the first step we have introduced the definition
and in the second step made use of a singular value decomposition of , where is considered as a composite index,
Comparing the right- and left-hand side of Eq. (37) shows that
holds.
This is a quite remarkable result, as it shows how to compute the right-hand side SVD matrix V and the singular values λ of the matrix B without performing an SVD of B, but by performing an SVD of the much smaller matrix instead. The matrix has the dimension with . This number is usually not prohibitively large and is much smaller than the corresponding B index , because in an ML calculation the numbers of logical modes, pz, are usually small. (For a binary tree, pz = 2 holds throughout.) Finally we note that when building , Eq. (38) in practice, we employ the regularized singular values rather than the original ones. This improves the stability of the SVD of , Eq. (39).
What is still missing is a recursion for U. For this, we combine the SVD of B, Eq. (33), with the recursion of B, Eq. (17),
Comparing the first and the last line shows the recursion
where, as before, .
Of course, one does not want to generate the huge matrix explicitly; one rather computes the modified mean-field [compare with Eq. (23)]
recursively. Proceeding similarly as when deriving Eq. (30), and making use of Eqs. (17), (29), (42), and (43), one arrives at the desired recursion for the modified mean-field matrix
where, as before, the symbol denotes a composite index similar to J, but with the κlth and νth index missing. And tensor is to be interpreted as , but with the κlth and νth entry of J replaced with n and jν. And a similar notational convention is adopted for tensor u. Finally, the EOM for the coefficients [compare with Eq. (31)] reads
The working equations for the recursions are (25), (38)–(40), and (44). First, the h-matrix elements are built bottom-up (25) and then the top-level mean-field matrices are built. These are given by the usual MCTDH formulas and their generation requires a SVD of the top-layer coefficient tensor, . Then the top-down recursions for the mean-fields (44) are performed, which requires the building of and its singular-value decomposition (38)–(40). After that, the EOM for the coefficients, Eqs. (8) and (45), is to be solved.
III. NUMERICAL EXAMPLES
Similar as in the previous work,63 we test the new regularization scheme on the spin-boson model.65,66 This is a system-bath type Hamiltonian that has been widely used in the context of electron transfer theories. The model has two electronic states (|ϕ1⟩ and |ϕ2⟩) linearly coupled to a bath of harmonic oscillators. Using mass-weighted coordinates, the Hamiltonian reads
where σx and σz are Pauli matrices
The properties of the bath that influence the dynamics of the two-state subsystem are specified by the spectral density function65,66
In the test examples below, we use an Ohmic (linear) spectral density with an exponential cutoff
where α is the dimensionless Kondo parameter that characterizes the system-bath coupling strength and ωc is the cutoff frequency of the bath. The continuous bath spectral density of Eq. (49) can be discretized to the form of Eq. (48) via the relation24,27
in which the density of frequencies ρ(ω) is defined from the integral relation
In this work, ρ(ω) is chosen as
where ωN+1 = ∞ is removed from the simulation. The number of bath modes N in the discretization scheme is a convergence parameter and needs to be sufficiently large to represent the continuum within the time scale of interest. Unlike in our previous tests using MCTDH,63 ML-MCTDH is a much more powerful method that can easily converge N. The examples below employ a few hundred to a few thousand bath modes and treat physical regimes that have intermediate to very strong coupling strength. The number of layers ranges from 3 to 6 in the simulation. These examples represent challenging cases of realistic ML-MCTDH simulations.
In the examples below, we choose a simple initial state with the bath starting at its ground state. The dynamical property examined here is the time-dependent population (difference) of the two-level subsystem
We first consider an example with parameters D/Δ = 0, ωc/Δ = 25, and an intermediate coupling strength, α = 0.5. Figure 1 illustrates convergence of the result with respect to the regularization parameter ϵ for a bath of 250 modes. A three-layer ML-MCTDH was used in the simulation. The top layer has five single particle groups (logical coordinates) with six SPFs per group and a mode with two SPFs for the electronic DOF, resulting in a total of 65 × 2 = 15 552 configurations. For the second layer, there are four single particle groups per each top layer group and five SPFs are used per each second layer group. For the bottom layer, up to 13 modes are combined through use mode combination3 with the primitive basis functions (three per DOF) adiabatically contracted.23 Among the six SPFs for each top layer group and five SPFs for each second layer group, the first one represents the initial state. The remaining SPFs are randomly selected and orthonormalized against the first SPF. A variable stepsize predictor-corrector ordinary differential equation (ODE) solver was used in the time integration. Both the matrix diagonalization (of the reduced density matrix) in the standard scheme and the SVD in the new scheme are handled by the corresponding LAPACK (Linear Algebra Package) subroutines.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 250 modes is used in a three-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). The ϵ = 10−8 curve lies on top of the ϵ = 10−10 one.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 250 modes is used in a three-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). The ϵ = 10−8 curve lies on top of the ϵ = 10−10 one.
From Fig. 1, it looks like that for the standard regularization of the reduced density matrices, ϵ ≤ 10−8 is required to converge the result. A more quantitative measure of the convergence can be obtained by introducing a relative cumulative deviation (hereafter simply referred to as “deviation”), defined as
where τ is the time scale of the simulation, Pref(t) is the reference result obtained by using the smallest ϵ in the new regularization scheme, and Pmax and Pmin are the maximum and minimum values of Pref(t). The error tolerance of the predictor-corrector ODE solver was set such that there is ∼10−5 roundoff error in the deviation. Any deviation below, say 5 × 10−6, can be considered as allowable numerical noise. For the numerical convergence of P(t), we set the criterion that the deviation be less than ∼10−4, which is stringent enough as compared with other convergence parameters and is still sufficiently larger than the roundoff error.
In Table I, the convergence behavior using the new regularization is compared with the one using the standard regularization scheme. As discussed in our previous work,63 reducing ϵ may thus have two opposite effects: (1) a smaller ϵ results in a stiffer set of EOMs and thus will increase the computational cost; (2) by reducing ϵ, the correct SPF subspace is more quickly reached and the integrator will work at a true physical time scale thereafter and thus the overall effect is that the computational cost is reduced. Table I examines these two opposite effects by displaying the number of calls to the ML-MCTDH derivative subroutine in which time derivatives of the A-vectors of all layers are evaluated. The new scheme quickly rotates the unoccupied SPFs to the “correct” direction and is thus more effective, under which convergence is reached when ϵ ≤ 10−8. The standard scheme requires ϵ ≤ 10−9 to achieve the same convergence. Results with a larger ϵ in the new regularization scheme do not deviate as much as those obtained in the standard scheme. This shows that the new scheme is more robust.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 1. A three-layer ML-MCTDH is used to treat a bath of 250 modes.
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−4 | 18 238 | 2.4 × 10−2 | 17 744 | 6.1 × 10−3 |
| 10−6 | 17 895 | 1.2 × 10−2 | 20 840 | 2.8 × 10−3 |
| 10−7 | 21 411 | 4.1 × 10−3 | 21 314 | 1.8 × 10−4 |
| 10−8 | 17 265 | 4.3 × 10−4 | 16 541 | 1.4 × 10−4 |
| 10−9 | 15 730 | 1.6 × 10−4 | 16 198 | 7.5 × 10−5 |
| 10−10 | 15 798 | 9.0 × 10−5 | 15 758 | 2.9 × 10−5 |
| 10−11 | 12 588 | 3.7 × 10−5 | 12 555 | 1.2 × 10−5 |
| 10−12 | 10 483 | 1.5 × 10−5 | 10 282 | 4.9 × 10−6 |
| 10−13 | 10 885 | 8.4 × 10−6 | 10 259 | 2.7 × 10−6 |
| 10−14 | 19 356 | 4.2 × 10−6 | 10 502 | 1.8 × 10−6 |
| 10−15 | 44 565 | 2.1 × 10−6 | 10 782 | 1.2 × 10−6 |
| 10−16 | 191 629 | 2.7 × 10−6 | 12 583 | 9.3 × 10−7 |
| 10−18 | … | … | 13 305 | 2.4 × 10−6 |
| 10−20 | … | … | 17 841 | Reference |
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−4 | 18 238 | 2.4 × 10−2 | 17 744 | 6.1 × 10−3 |
| 10−6 | 17 895 | 1.2 × 10−2 | 20 840 | 2.8 × 10−3 |
| 10−7 | 21 411 | 4.1 × 10−3 | 21 314 | 1.8 × 10−4 |
| 10−8 | 17 265 | 4.3 × 10−4 | 16 541 | 1.4 × 10−4 |
| 10−9 | 15 730 | 1.6 × 10−4 | 16 198 | 7.5 × 10−5 |
| 10−10 | 15 798 | 9.0 × 10−5 | 15 758 | 2.9 × 10−5 |
| 10−11 | 12 588 | 3.7 × 10−5 | 12 555 | 1.2 × 10−5 |
| 10−12 | 10 483 | 1.5 × 10−5 | 10 282 | 4.9 × 10−6 |
| 10−13 | 10 885 | 8.4 × 10−6 | 10 259 | 2.7 × 10−6 |
| 10−14 | 19 356 | 4.2 × 10−6 | 10 502 | 1.8 × 10−6 |
| 10−15 | 44 565 | 2.1 × 10−6 | 10 782 | 1.2 × 10−6 |
| 10−16 | 191 629 | 2.7 × 10−6 | 12 583 | 9.3 × 10−7 |
| 10−18 | … | … | 13 305 | 2.4 × 10−6 |
| 10−20 | … | … | 17 841 | Reference |
Table I illustrates the similarity between the ML-MCTDH and MCTDH regularizations: a quick rotation (small ϵ) of the SPF subspaces is necessary in the beginning of the integration to increase both the accuracy and efficiency. Once ϵ is small enough, the number of calls to the derivative subroutine remains approximately the same (and in fact slightly goes down as for the MCTDH case) within a certain range. As ϵ becomes very small (ϵ = 10−16), the ML-MCTDH EOMs become very stiff under the standard regularization scheme and the computational cost increases by an order of magnitude. By contrast, under the new regularization scheme, the cost remains approximately the same over a large range of ϵ. This suggests robustness of the algorithm in terms of both accuracy and efficiency. Again, the difference lies on the fact that though ϵ formally represents the same parameter in the two schemes, it is used differently in the numerical implementation: in the standard scheme, ϵ−1 is the effective regularization parameter whereas in the new scheme it is ϵ−1/2. An even smaller ϵ = 10−24 may still be handled in the new regularization scheme.
We now consider a longer time dynamics than Fig. 1, which also requires a larger number of modes. Figure 2 illustrates convergence of P(t) for twice as long as in Fig. 1 and for a bath of 1000 modes. A four-layer ML-MCTDH was used in the simulation. The top layer still has five single particle groups, but eight SPFs are used per each group, resulting in a total of 85 × 2 = 65 536 configurations. For the other layers, there are up to four lower single particle groups per each upper layer group and five SPFs are used per each group. Mode combination and adiabatic basis contraction are again used for each bottom layer group, which holds up to 13 modes.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 1000 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−6 curve is completely off and that the ϵ = 10−10 and ϵ = 10−12 curves are on top of each other.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 1000 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−6 curve is completely off and that the ϵ = 10−10 and ϵ = 10−12 curves are on top of each other.
Since numerical singularities occur at all layers in ML-MCTDH, convergence requires smaller ϵ as the number of layers and/or SPFs increases. Figure 2 shows that under the standard regularization scheme, the simulation with ϵ = 10−6 predicts a completely erroneous dynamics—a coherent motion of P(t) near its initial value—that contradicts the converged result even qualitatively. The simulation with ϵ = 10−7 also predicts an incorrect initial delay. A more detailed analysis in Table II shows that a large ϵ also requires more integration steps due to its wrong dynamics.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 2. A four-layer ML-MCTDH is used to treat a bath of 1000 modes.
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−6 | 132 677 | 6.6 × 10−1 | 63 404 | 5.6 × 10−2 |
| 10−7 | 50 285 | 7.2 × 10−2 | 49 341 | 1.2 × 10−2 |
| 10−8 | 45 403 | 1.7 × 10−2 | 49 635 | 1.8 × 10−3 |
| 10−9 | 41 964 | 2.7 × 10−3 | 47 184 | 2.2 × 10−4 |
| 10−10 | 31 844 | 3.1 × 10−4 | 31 747 | 7.5 × 10−5 |
| 10−11 | 28 998 | 8.1 × 10−5 | 30 976 | 3.7 × 10−5 |
| 10−12 | 29 100 | 4.0 × 10−5 | 29 056 | 1.5 × 10−5 |
| 10−13 | 29 422 | 1.4 × 10−5 | 28 965 | 5.2 × 10−6 |
| 10−14 | 33 765 | 9.2 × 10−6 | 29 205 | 4.4 × 10−6 |
| 10−15 | 69 370 | 7.3 × 10−6 | 30 641 | 3.8 × 10−6 |
| 10−16 | 277 915 | 7.1 × 10−6 | 29 996 | 4.4 × 10−6 |
| 10−18 | … | … | 31 292 | 3.6 × 10−6 |
| 10−20 | … | … | 35 347 | 4.9 × 10−6 |
| 10−22 | … | … | 40 092 | Reference |
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−6 | 132 677 | 6.6 × 10−1 | 63 404 | 5.6 × 10−2 |
| 10−7 | 50 285 | 7.2 × 10−2 | 49 341 | 1.2 × 10−2 |
| 10−8 | 45 403 | 1.7 × 10−2 | 49 635 | 1.8 × 10−3 |
| 10−9 | 41 964 | 2.7 × 10−3 | 47 184 | 2.2 × 10−4 |
| 10−10 | 31 844 | 3.1 × 10−4 | 31 747 | 7.5 × 10−5 |
| 10−11 | 28 998 | 8.1 × 10−5 | 30 976 | 3.7 × 10−5 |
| 10−12 | 29 100 | 4.0 × 10−5 | 29 056 | 1.5 × 10−5 |
| 10−13 | 29 422 | 1.4 × 10−5 | 28 965 | 5.2 × 10−6 |
| 10−14 | 33 765 | 9.2 × 10−6 | 29 205 | 4.4 × 10−6 |
| 10−15 | 69 370 | 7.3 × 10−6 | 30 641 | 3.8 × 10−6 |
| 10−16 | 277 915 | 7.1 × 10−6 | 29 996 | 4.4 × 10−6 |
| 10−18 | … | … | 31 292 | 3.6 × 10−6 |
| 10−20 | … | … | 35 347 | 4.9 × 10−6 |
| 10−22 | … | … | 40 092 | Reference |
Table II shows that the standard regularization scheme requires ϵ ≤ 10−11 to achieve convergence, whereas the new scheme requires ϵ ≤ 10−10. Results with a larger ϵ in the new regularization scheme do not deviate as much as those obtained in the standard scheme. Furthermore, the computational cost under the new scheme remains almost the same spanning a large range of the regularization parameter, ϵ = 10−20–10−10. All this suggests again the stability of the new regularization scheme.
Examples above are typical applications of ML-MCTDH. Due to the fact that the new regularization scheme performs SVD on the matrix unfoldings of the A-vectors of all layers, instead of diagonalizing smaller density matrices in the standard scheme, it is expected that the new scheme is slightly more expensive if the numbers of calls to the derivative subroutines are the same under the two schemes. Indeed we find that the new scheme takes 25% more time than the standard scheme, e.g., for ϵ = 10−10 in the example above. However, for ϵ ≤ 10−15, the new scheme is clearly more efficient than the standard scheme. As shown below, only the new scheme is practical for investigating difficult regimes.
We now consider challenging applications. Figure 3 shows such an example for D/Δ = 0, α = 2, and ωc/Δ = 25. Dynamics of P(t) exhibit localization after a short transient relaxation. Convergence is shown with respect to the regularization parameter ϵ for a four-layer ML-MCTDH simulation on a model bath of 500 modes. The top layer has four single particle groups, with eight SPFs per each group. For the other layers, there are up to four lower single particle groups per each upper layer group and five SPFs are used per each group. Mode combination and adiabatic basis contraction are used for each bottom layer group. Due to the stronger coupling strength, four primitive basis functions are needed for each mode. Hence each bottom layer group contains only up to 8 modes.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 2, and ωc/Δ = 25. A bath of 500 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−8 curve is qualitatively wrong. The ϵ = 10−16 and ϵ = 10−17 curves are on top of each other.
Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 2, and ωc/Δ = 25. A bath of 500 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−8 curve is qualitatively wrong. The ϵ = 10−16 and ϵ = 10−17 curves are on top of each other.
As shown in Fig. 3 and Table III, convergence of P(t) is slow, and sometimes even oscillatory, versus ϵ. To achieve deviation less than 10−4 under the standard scheme, ϵ ≤ 10−16 is required. Such a small ϵ results in more than an order of magnitude more computational cost in the standard scheme. On the contrary, not only does the new scheme give almost the same convergence with ϵ ≤ 10−14, but it also requires much less computational effort. Even pushing to a very small ϵ = 10−22, the cost is only three times as expensive.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A four-layer ML-MCTDH is used to treat a bath of 500 modes.
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 7 166 | 3.4 × 10−1 | 7 168 | 3.7 × 10−2 |
| 10−9 | 5 713 | 1.3 × 10−2 | 6 181 | 3.1 × 10−2 |
| 10−10 | 5 658 | 4.0 × 10−2 | 5 943 | 1.3 × 10−2 |
| 10−11 | 5 673 | 1.9 × 10−2 | 6 205 | 4.6 × 10−3 |
| 10−12 | 5 644 | 7.0 × 10−3 | 6 274 | 1.5 × 10−3 |
| 10−13 | 6 370 | 2.4 × 10−3 | 6 065 | 4.9 × 10−4 |
| 10−14 | 12 252 | 7.8 × 10−4 | 5 541 | 1.6 × 10−4 |
| 10−15 | 43 570 | 2.7 × 10−4 | 7 560 | 4.5 × 10−5 |
| 10−16 | 217 340 | 1.0 × 10−4 | 6 865 | 1.8 × 10−5 |
| 10−17 | 1 153 089 | 5.0 × 10−5 | 8 200 | 2.4 × 10−6 |
| 10−18 | … | … | 9 867 | 3.6 × 10−6 |
| 10−19 | … | … | 10 943 | 8.3 × 10−6 |
| 10−20 | … | … | 13 284 | 4.8 × 10−6 |
| 10−21 | … | … | 16 019 | 3.2 × 10−6 |
| 10−22 | … | … | 19 501 | Reference |
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 7 166 | 3.4 × 10−1 | 7 168 | 3.7 × 10−2 |
| 10−9 | 5 713 | 1.3 × 10−2 | 6 181 | 3.1 × 10−2 |
| 10−10 | 5 658 | 4.0 × 10−2 | 5 943 | 1.3 × 10−2 |
| 10−11 | 5 673 | 1.9 × 10−2 | 6 205 | 4.6 × 10−3 |
| 10−12 | 5 644 | 7.0 × 10−3 | 6 274 | 1.5 × 10−3 |
| 10−13 | 6 370 | 2.4 × 10−3 | 6 065 | 4.9 × 10−4 |
| 10−14 | 12 252 | 7.8 × 10−4 | 5 541 | 1.6 × 10−4 |
| 10−15 | 43 570 | 2.7 × 10−4 | 7 560 | 4.5 × 10−5 |
| 10−16 | 217 340 | 1.0 × 10−4 | 6 865 | 1.8 × 10−5 |
| 10−17 | 1 153 089 | 5.0 × 10−5 | 8 200 | 2.4 × 10−6 |
| 10−18 | … | … | 9 867 | 3.6 × 10−6 |
| 10−19 | … | … | 10 943 | 8.3 × 10−6 |
| 10−20 | … | … | 13 284 | 4.8 × 10−6 |
| 10−21 | … | … | 16 019 | 3.2 × 10−6 |
| 10−22 | … | … | 19 501 | Reference |
ML-MCTDH simulations become more challenging for larger systems. Table IV analyzes the result with the same physical parameters as in Table III, but for a bath of 2000 modes. A five-layer ML-MCTDH is employed here. Under the standard scheme, convergence can be barely reached with an ϵ as small as 10−18. More severely, the computational effort increases here by more than two orders of magnitude. Things were quite different when the new regularization scheme was used. Convergence can be achieved within a broad range of ϵ = 10−24–10−16 and with very reasonable computational effort.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A five-layer ML-MCTDH is used to treat a bath of 2000 modes.
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 55 007 | 5.1 × 10−1 | 16 298 | 4.9 × 10−1 |
| 10−9 | 42 682 | 5.0 × 10−1 | 14 523 | 1.0 × 10−1 |
| 10−10 | 15 014 | 3.0 × 10−1 | 11 534 | 4.3 × 10−2 |
| 10−11 | 13 456 | 1.7 × 10−2 | 9 669 | 3.0 × 10−2 |
| 10−12 | 12 112 | 3.8 × 10−2 | 8 484 | 1.2 × 10−2 |
| 10−13 | 12 124 | 1.9 × 10−2 | 8 234 | 4.4 × 10−3 |
| 10−14 | 15 801 | 7.1 × 10−3 | 8 939 | 1.5 × 10−3 |
| 10−15 | 30 988 | 2.5 × 10−3 | 9 985 | 4.8 × 10−4 |
| 10−16 | 149 421 | 8.4 × 10−4 | 10 711 | 1.5 × 10−4 |
| 10−17 | 635 887 | 2.7 × 10−4 | 9 616 | 3.5 × 10−5 |
| 10−18 | 2 279 470 | 1.9 × 10−4 | 9 635 | 1.1 × 10−5 |
| 10−20 | … | … | 14 895 | 2.2 × 10−5 |
| 10−22 | … | … | 20 842 | 5.0 × 10−6 |
| 10−24 | … | … | 32 983 | Reference |
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 55 007 | 5.1 × 10−1 | 16 298 | 4.9 × 10−1 |
| 10−9 | 42 682 | 5.0 × 10−1 | 14 523 | 1.0 × 10−1 |
| 10−10 | 15 014 | 3.0 × 10−1 | 11 534 | 4.3 × 10−2 |
| 10−11 | 13 456 | 1.7 × 10−2 | 9 669 | 3.0 × 10−2 |
| 10−12 | 12 112 | 3.8 × 10−2 | 8 484 | 1.2 × 10−2 |
| 10−13 | 12 124 | 1.9 × 10−2 | 8 234 | 4.4 × 10−3 |
| 10−14 | 15 801 | 7.1 × 10−3 | 8 939 | 1.5 × 10−3 |
| 10−15 | 30 988 | 2.5 × 10−3 | 9 985 | 4.8 × 10−4 |
| 10−16 | 149 421 | 8.4 × 10−4 | 10 711 | 1.5 × 10−4 |
| 10−17 | 635 887 | 2.7 × 10−4 | 9 616 | 3.5 × 10−5 |
| 10−18 | 2 279 470 | 1.9 × 10−4 | 9 635 | 1.1 × 10−5 |
| 10−20 | … | … | 14 895 | 2.2 × 10−5 |
| 10−22 | … | … | 20 842 | 5.0 × 10−6 |
| 10−24 | … | … | 32 983 | Reference |
For an even larger bath of 5000 modes, Table V shows that the standard scheme is unable to achieve convergence in a practical sense (because of the two orders of magnitude more computational effort): for ϵ = 10−17, the simulation under the new regularization scheme takes less than 4 h (14 438 calls to the ML-MCTDH derivatives) and almost reaches convergence (deviation ∼10−4), whereas under the standard scheme, ϵ = 10−18 is required to achieve similar convergence and the simulation takes a month (3 046 997 calls to the ML-MCTDH derivatives). In such a situation, the new regularization scheme is the only viable choice.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A six-layer ML-MCTDH is used to treat a bath of 5000 modes.
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 77 628 | 5.1 × 10−1 | 80 769 | 5.0 × 10−1 |
| 10−9 | 88 948 | 5.1 × 10−1 | 28 753 | 4.8 × 10−1 |
| 10−10 | 89 720 | 4.8 × 10−1 | 25 155 | 7.7 × 10−2 |
| 10−11 | 23 541 | 2.0 × 10−1 | 20 498 | 4.4 × 10−2 |
| 10−12 | 23 581 | 3.7 × 10−2 | 14 441 | 2.8 × 10−2 |
| 10−13 | 16 476 | 3.5 × 10−2 | 11 120 | 1.1 × 10−2 |
| 10−14 | 17 973 | 1.5 × 10−2 | 9 946 | 3.8 × 10−3 |
| 10−15 | 38 166 | 5.4 × 10−3 | 8 907 | 1.3 × 10−3 |
| 10−16 | 176 864 | 1.8 × 10−3 | 10 738 | 4.2 × 10−4 |
| 10−17 | 664 745 | 6.2 × 10−4 | 14 438 | 1.2 × 10−4 |
| 10−18 | 3 046 997 | 1.3 × 10−4 | 13 250 | 3.0 × 10−5 |
| 10−20 | … | … | 13 698 | 1.5 × 10−5 |
| 10−22 | … | … | 19 282 | 7.8 × 10−6 |
| 10−24 | … | … | 41 213 | Reference |
| Standard regularization | New regularization | |||
|---|---|---|---|---|
| ϵ | Number of calls | Deviation | Number of calls | Deviation |
| 10−8 | 77 628 | 5.1 × 10−1 | 80 769 | 5.0 × 10−1 |
| 10−9 | 88 948 | 5.1 × 10−1 | 28 753 | 4.8 × 10−1 |
| 10−10 | 89 720 | 4.8 × 10−1 | 25 155 | 7.7 × 10−2 |
| 10−11 | 23 541 | 2.0 × 10−1 | 20 498 | 4.4 × 10−2 |
| 10−12 | 23 581 | 3.7 × 10−2 | 14 441 | 2.8 × 10−2 |
| 10−13 | 16 476 | 3.5 × 10−2 | 11 120 | 1.1 × 10−2 |
| 10−14 | 17 973 | 1.5 × 10−2 | 9 946 | 3.8 × 10−3 |
| 10−15 | 38 166 | 5.4 × 10−3 | 8 907 | 1.3 × 10−3 |
| 10−16 | 176 864 | 1.8 × 10−3 | 10 738 | 4.2 × 10−4 |
| 10−17 | 664 745 | 6.2 × 10−4 | 14 438 | 1.2 × 10−4 |
| 10−18 | 3 046 997 | 1.3 × 10−4 | 13 250 | 3.0 × 10−5 |
| 10−20 | … | … | 13 698 | 1.5 × 10−5 |
| 10−22 | … | … | 19 282 | 7.8 × 10−6 |
| 10−24 | … | … | 41 213 | Reference |
Finally we consider a variant of the spin-boson model, the polaron-transformed Hamiltonian of (46),
where is a translational operator for different electronic states
This is a pathological example and has been shown earlier by us63 that the standard regularization scheme completely fails in MCTDH for a bath of 40 modes and the new regularization scheme needs be used to obtain the correct result. The phase in the translational operator (55) requires a quick rotation of the single particle space to capture the correct dynamics. As shown in (55), this effect is more pronounced for a stronger coupling strength (larger cj’s) and/or a larger bath (more terms in the product and thus the exponential summation).
Here we consider a similar parameter set as in our previous MCTDH work,63 but with a larger bath of 512 modes. A four-layer ML-MCTDH was used in the simulation. The top layer contains five single particle groups, with eight SPFs used per each group, resulting in a total of 85 × 2 = 65 536 configurations. For other layers, there are 3–4 lower single particle groups per each upper layer group and eight SPFs are used per each group. Mode combination and adiabatic basis contract are again used for each bottom layer group, which holds up to 9 modes.
Under the standard scheme, P(t) remains at its initial value of 1 as time evolves even with the smallest possible ϵ that does not crash the ODE solver, indicating the inability of rotating the single particle space to the “correct” direction. On the contrary, the new scheme is capable of achieving convergence when ϵ ≤ 10−18. As shown in Table VI, the computational cost of the new regularization scheme is quite reasonable within a large range of ϵ. As the coupling strength and number of modes increase, this model will become so ill-conditioned that the new scheme will eventually fail. For practical simulations using ML-MCTDH, one needs to seek a better representation or reformulation of the problem when possible (e.g., in this case, go back to the standard spin-boson model). The fact that the new scheme can handle a much broader regime than the standard scheme provides more flexibility for realistic applications and reformulations of the problem.
Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ for the polaron-transformed spin-boson model. A four-layer ML-MCTDH is used to treat a bath of 512 modes. The physical parameter are D/Δ = 0, α = 2, and ωc/Δ = 2.5. Only the results under the new regularization scheme are shown since the standard scheme completely fails.
| ϵ | Number of calls | Deviation |
|---|---|---|
| 10−10 | 2747 | 6.7 × 10−1 |
| 10−12 | 2077 | 1.6 × 10−1 |
| 10−14 | 2498 | 1.6 × 10−2 |
| 10−16 | 3130 | 1.6 × 10−3 |
| 10−18 | 4386 | 1.6 × 10−4 |
| 10−20 | 6376 | 4.5 × 10−5 |
| 10−22 | 6144 | 2.9 × 10−5 |
| 10−24 | 6221 | Reference |
| ϵ | Number of calls | Deviation |
|---|---|---|
| 10−10 | 2747 | 6.7 × 10−1 |
| 10−12 | 2077 | 1.6 × 10−1 |
| 10−14 | 2498 | 1.6 × 10−2 |
| 10−16 | 3130 | 1.6 × 10−3 |
| 10−18 | 4386 | 1.6 × 10−4 |
| 10−20 | 6376 | 4.5 × 10−5 |
| 10−22 | 6144 | 2.9 × 10−5 |
| 10−24 | 6221 | Reference |
In passing, we note that Manthe has recently also discussed61 the polaron-transformed spin-boson model for a somewhat weaker coupling strength. He correctly solved the dynamical problem by using optimized60 initially unoccupied SPFs. Hence this approach may serve as an alternative to the use of the new regularization scheme. It is also possible, and presumably more effective, to combine Manthe’s choice of initial SPFs with our new regularization scheme.
IV. CONCLUDING REMARKS
In this work, our previously proposed new regularization of the MCTDH equations of motion is generalized to ML-MCTDH. In this scheme, the expansion coefficients of all layers, instead of the density matrices, in ML-MCTDH are regularized through their matrix unfoldings. This automatically regularizes the density matrices. Despite the principle idea remaining the same, the ML-MCTDH equations are much more complicated and this generalization is not obvious. For practical purposes, we have derived a recursive implementation of the new regularization scheme.
We have tested the new regularization scheme on several examples of the spin-boson model whose parameter set ranges from intermediate to very strong couplings and the size of the system ranges from a few hundred to 5000 degrees of freedom. These are very realistic applications of ML-MCTDH. Since singularities occur more frequently in ML-MCTDH, it is more important to adopt the new regularization scheme.
In general, it is found that the new scheme is much more robust with respect to changing the regularization parameter ϵ. For weak to intermediate coupling strength, it is still possible to achieve convergence in the standard scheme by reducing ϵ. However, this becomes prohibitively expensive in challenging applications where the coupling strength becomes strong. In those situations, the new regularization scheme becomes the only practical choice. Furthermore, it is often preferred not to check so many parameters in real applications of ML-MCTDH. The fact that the new regularization scheme is stable over a large range of the regularization parameter ϵ makes it possible to use a conservative default value and not worry about artifacts induced otherwise.
Because the coefficient tensors in ML-MCTDH are usually smaller than in MCTDH, the cost of SVDs in the new scheme is often only moderately more expensive than the density matrix diagonalization in the standard scheme for routine applications (for challenging applications, the latter is neither reliable nor practical.) Still, owing to the fact that orbital rotations are the most important at the beginning, it is possible to combine the new regularization scheme with the standard scheme such that the latter can be switched on when the natural populations become large enough.
In some previous work, the singularity for the short time (ML)-MCTDH propagation raises concerns and effort has been made to remove it. Our work shows that this singularity is needed to rapidly rotate the unoccupied SPFs into their optimal direction in Hilbert space. If this is not done properly, not only the accuracy but also the efficiency of the ML-MCTDH propagation will be affected. Many times, a more rapid orbital rotation at the beginning will actually reduce the overall computational cost (see, e.g., Tables II, IV, and V) because of the exploration of the correct regions of the Hilbert space. The new regularization offers an efficient way of achieving this goal.
ACKNOWLEDGMENTS
H.W. acknowledges the support from the National Science Foundation Grant No. CHE-1500285 and resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. H.-D.M. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG).
