Propagating multi-dimensional density operators using the multi-layer-ρ multi-configurational time-dependent Hartree method

Solving the Liouville–von-Neumann equation using a density operator provides a more complete picture of dynamical quantum phenomena than by using a wavepacket and solving the Schrödinger equation. As density operators are not restricted to the description of pure states, they can treat both thermalized and open systems. In practice, however, they are rarely used to study molecular systems as the computational resources required are even more prohibitive than those needed for wavepacket dynamics. In this paper, we demonstrate the potential utility of a scheme based on the powerful multi-layer multi-configurational time-dependent Hartree algorithm for propagating multi-dimensional density operators. Studies of two systems using this method are presented at a range of temperatures and including up to 13 degrees of freedom. The first case is single proton transfer in salicylaldimine, while the second is double proton transfer in porphycene. A comparison is also made with the approach of using stochastic wavepackets


I. INTRODUCTION
Quantum dynamics simulations, which solve the timedependent Schrödinger equation (TDSE), have become standard computational technique to understand the quantum behavior of molecular systems.In particular, they are used to great effect in describing the time-evolution of a system after photo-excitation and in helping to unravel the signal provided by time-resolved spectroscopic experiments.They also have utility in describing the fundamental reactivity of atoms and molecules and in other situations, such as proton and electron transfer where quantum effects dominate the dynamics.
Solving the TDSE yields the system wavefunction as an evolving superposition of states, known as a wavepacket.A wavepacket, however, is restricted to the description of pure states in which the superposition of states is coherent.While this is a complete description for closed-system dynamics at 0 K, what happens when the system is open (i.e., interacts with an environment) or is at finite temperature?Accurately modeling and predicting the dynamics of thermalized open systems are key to the understanding of almost all chemical and biological processes as no process occurring in nature is truly ever a closed system at 0 K.
Recent studies have shown that temperature effects can play an important role in the quantum dynamics in a range of systems, such as nanoparticles, 1 surface adsorption, 2 and liquid water. 3They are also important to correctly describe the electron transfer process in solar cells, which clearly do not have working temperatures of 0 K.There has even been evidence that dye sensitized solar cells, unlike first generation silicon based solar cells, actually increase in efficiency with temperature. 4,5Therefore, investigating these systems at working temperatures is critically important for the goal of creating more efficient, renewable energy sources.
The TDSE is no longer sufficient to simulate these types of systems, which must be described by an incoherent superposition of states.A density operator representation is, therefore, required for a full treatment, and the Liouville-von-Neumann (LvN) The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpequation must be solved in place of the TDSE.][8][9] The disadvantage, however, is that the computational effort required to solve the LvN scales quadratically with system size compared to the TDSE.Given that the latter is a hard computational problem with an underlying exponential scaling with system size, density operators have not seen much use in recent years as they are restricted to very small systems.
1][12] By providing a time-dependent tensor contraction scheme for the underlying primitive basis used in quantum dynamics simulations, it is able to treat much larger systems than conventional numerical solutions to the TDSE, which require huge grids.In the multi-layer (ML-MCTDH) form, [13][14][15] the TDSE for systems with over 1000 degrees of freedom has been treated. 168][19] Here, we revisit this direction of research using an ML-MCTDH expansion to show how multi-dimensional density matrices may be propagated.
A number of other methods have been developed, which allow for the simulation of multi-dimensional dynamics at finite temperatures.Some are based on the MCTDH method for solving the timedependent Schrödinger equation, [20][21][22][23] such as statistical MCTDH approaches 24,25 and thermofield MCTDH. 26Statistical approaches take an average over randomly sampled thermalized wavepackets to build up the density matrix, thus keeping the scaling as for solving the TDSE, although many runs may be required for a converged result.The Thermofield approach maps the density matrix dynamics onto a set of states, including an effective set of bath mode, turning the problem into propagating a wavefunction with a doubling of the number of vibrational modes.Other approaches are based on tensor trains 27 and the density matrix renormalization group (DMRG). 28,29hese methods use tensor contraction schemes to break the natural exponential scaling of solving the TDSE.They can be related to the ML-MCTDH method, which uses a more general tensor contraction scheme. 30roton transfer (PT) is a temperature sensitive process and is a common rate limiting step in enzymes and other biological systems. 31,32The light mass of the proton means tunneling will be important at temperatures below the reaction barrier, and so PT needs to be modeled using quantum dynamics simulations for a correct description.
In this paper, after an overview of the theory, two systems will be treated.The first is salicylaldimine, which is a rigid molecule with an asymmetric PT between a nitrogen and oxygen atom.It has been studied previously using wavepacket dynamics, 33 and the same simple model potentials will be used to show the scaling behavior of different density matrix approaches with up to 13 DOFs included in largest simulations.The temperature behavior as a function of system size will also be investigated.A second system based on porphycene, again using a Hamiltonian from previous wavepacket dynamics, 34 will then be investigated.This includes a double PT, and the behavior with the temperature including up to 10 degrees of freedom will be presented.

II. THEORY
A. Solving the TDSE If a system can be described by a pure state with wavefunction Ψ(t), its time-evolution is given by the TDSE, where Ĥ is the Hamiltonian of the system.To numerically solve the TDSE, one approach is to expand the multi-dimensional wavefunction into a time-independent "grid" basis and use the Dirac-Frenkel variational principle, 35 ⟨∂Ψ|H − i Arguably, the most efficient of these methods to emerge is the MCTDH method. 20,21n the following, the MCTDH method is described along with its multi-layer (ML-MCTDH) extension, before moving on to the density matrix variants, ρMCTDH and ML-ρMCTDH.In the MCTDH formalism, there is a gauge constraint required to ensure that the time-dependent basis functions remain orthonormal. 21In the following, the standard gauge in which the constraint vanishes is chosen.

MCTDH
The MCTDH approach expands the wavefunction into a directproduct of some chosen number p sets of time-dependent basis functions, φ (κ) , known as single-particle functions (SPFs), where A j 1 ... j p (t) are the time-dependent expansion coefficients and κ is the mode number in the system with coordinate Q κ .These mode coordinates are typically a set of system coordinates Q κ = {q 1 , q 2 , . ..}.The SPFs themselves are further expanded into a linear combination of static primitive basis functions, with timeindependent coefficients.These static basis functions usually take the form of a discrete variable representation (DVR), 21,36 where χ α (Q κ ) are the time-independent DVR functions.Expanding the wavefunction in this way allows the results to be converged, from the fast, but inaccurate, limiting case of Time-Dependent Hartree (TDH) with n k = 1, up to the complete solution with n k = N k . 21The latter is a numerically exact solution of the TDSE within the limits of accuracy given by the basis set and integration schemes used.We now introduce the compact notation for the wavefunction expansion, The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpwhere ΦJ is an n-dimensional Hartree product configuration, AJ are the expansion coefficients, J is a multi-index subscript J = j 1 , . . ., j p , and Q = {Q 1 , . . ., Q p }.By applying the Dirac-Frenkel variational principle to the MCTDH wavefunction ansatz, two coupled equations of motion (EOM) can be derived.The first represents the evolving coefficients, where ΦJ is a p-dimensional Hartree product of the SPFs (a configuration) and AJ are the expansion coefficients.
The second EOM, derived for the SPFs, is where P κ is a projection operator defined as and ρ κ j,l is a density matrix defined as where Ψ κ l is a "single-hole" wavefunction, which excludes the DOF κ, (10) Therefore, by using this single-hole notation, summing over all l gives the full wavefunction, Thus, ⟨Ψ κ j |Ψ κ l ⟩ results in the integration over all the DOF except the κth.Here, ⟨ Ĥ⟩ is the mean-field operator acting on the SPFs,

ML-MCTDH
Although the MCTDH approach is able to treat larger systems compared to the standard approach, which expresses the wavefunction directly on the time-independent grid, it still has limitations with regard to memory requirements and computational efficiency. 22To improve its utility, we move to a multi-layer version, termed ML-MCTDH.The key to this method is to expand each SPF recursively in an "MCTDH-like" form to create "layers." 13,16In what follows, the nomenclature developed by Manthe 13 is used.
To form the top layer, the wavefunction is expanded as before in MCTDH into configurations, ϕ J , that are products of SPFs, where superscript 1 denotes the first layer, d 1 is the dimensionality of this layer (number of sets of SPFs), and q 1:κ = (q 1 . . .q d κ ) is the set of coordinates for the κth mode (set of SPFs).
In the next layer, each SPF is expanded into a further set of lower-dimensional SPFs, denoted by superscript 2, The superscript 2 : κ denotes functions (and coefficients) in the second layer coming from the κth mode from the top layer.These functions can, in turn, be expanded as Further layers take the same form as this with respect to the layer above.Note that the superscript keeps track of where the function comes from, running down the expansion tree from the top layer (κμ) in addition to the layer number [3 in Eq. ( 15)] and function number in that layer (e.g., k 1 ).Therefore, this can be simplified notationally into a single index so that (15) can be written as Note also that coefficients require a subscript, j, to denote which function in the layer above they relate to, or in general, φ m:κ j (q m:κ , t) Note that the relationship to the layer above on the left-hand side φ m:κ j is implicit-this is simply the jth function of the κth mode of the mth layer with dimension d.The history on the right-hand side, φ m+1:κλ k λ , only explicitly goes up one layer-this is the k λ th function of the mode κλ in layer m + 1, where κλ means the λth set of functions used to expand the κth mode of the layer above.The functions {φ κλ k } span the set of coordinates q m+1:κλ , which are a subset of q m:κ .Equation ( 17) also defines layer configurations, ϕ m:κ k .As a final note on the structure of the ML-MCTDH wavefunction, it should be mentioned that the lowest level of the expansion is usually a time-independent DVR, as with standard MCTDH, but it could also be a Gaussian wavepacket (GWP) basis. 37t is useful to define layer single-hole functions (SHFs).As for MCTDH, the top layer wavefunction can be written as The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp with the SHF ψ 1:κ j being the wavefunction excluding the SPFs for mode κ.Expanding φ 1:κ j using (17) gives where (19) and (20) define the relationship between SHFs on different layers.The generalization of ( 20) is The equations of motion are obtained as usual using the Dirac-Frenkel variational principle.First, the top layer coefficients δA ′ J is varied, and taking the common MCTDH gauge constraint ⟨ϕJ| φK⟩ = 0 to retain orthonormality of the SPFs, we obtain the usual MCTDH equations of motion, Varying the basis functions is equivalent to varying all the expansion coefficients on all the layers-these are the remaining "parameters" that define the wavefunction.These can be found by using the SHF of ( 22).
Analogously to MCTDH, it is useful to define the layer density matrices, the layer mean-field matrices, and finally the layer SPF projector, Using these, we can write the general form of the EOMs for the expansion coefficients for some arbitrary mth layer (including the lowest layer) as Note that ( 28) is identical to the usual MCTDH equations of motion for the SPFs except that a layer dependent basis representation is being used.In fact, standard MCTDH can be thought of the special case where the ML tree consists of just one layer.Equations ( 23) and ( 28) are the ML-MCTDH EOMs, which can be implemented using the recursive structure of Manthe. 13his recursive expansion can be represented as a tree (e.g., see Fig. 1 in the supplementary material), and each lower layer set of SPFs is represented by further "nodes."Each node is connected to nodes in the lower layer, representing the expansion into a lowerdimensional set of SPFs, and at the bottom layer by the standard primitive basis expansion of the SPFs.Expressing the nuclear wavefunction in this layered way is very flexible and can vary in how many layers are required. 14he ML-MCTDH approach is more complicated than standard MCTDH, and for smaller systems, it is much less efficient. 14However, as the system size increases, ML-MCTDH combines modes into smaller groups, now with a new layer of more manageable coefficients.Further reading on the specifics of the EOM for the ML-MCTDH method is widely available in the literature, for instance, Wang et al. 16,38,39 and Meyer et al. 14,[40][41][42] presented particularly thorough derivations.

B. Solving the LvN equation
If a system is not in a pure state, one can no longer express the system as a wavefunction and it can only be fully described by a density operator, 43 where Pi is the probability of being in the pure state ψ i .Representing the system using probabilities offers a way of describing an incoherent superposition of states.These are the outcome of the thermalization of a system interacting with an environment. 44nstead of solving the Schrödinger equation as one would for a pure state propagation, the time-evolution of the density matrix is given by the Liouville-von-Neumann equation, where L is the Liouvillian superoperator.
A variational principle analogous to the Dirac-Frenkel principle has been used to solve the LvN within the MCTDH formalism. 44he ρMCTDH method, therefore, operates in a Liouville space rather than in a Hilbert space and uses the Hilbert-Schmidt norm in place of the usual scalar product.For two matrices (operators), A and B, this norm is defined by the trace of the product Two formalisms for ρMCTDH have been developed, which are distinguished by the type of basic functions used.These are called types I and II, and their mathematical and numerical properties were examined by Raab and Meyer. 18,19

ρMCTDH(I)
Type I density operators are analogous to the MCTDH scheme for wavefunction expansion, as in Eq. ( 3), but instead, the density operator is expanded into single-particle density operators (SPDOs), . . .
The Journal of Chemical Physics As a density operator is Hermitian, the B coefficients must be real and the SPDOs are Hermitian.These properties are conserved throughout the propagation. 22As with the standard wavefunction MCTDH, this representation is not unique and constraints are necessary to ensure orthonormality.The following quantities in ρ-MCTDH are defined analogously to (8), (9), and (12).First, the new ρMCTDH density matrix is defined as The projector operator is and finally, the mean-field Liouvillian superoperator is where Π κ μ is a single-hole density matrix, similar to the MCTDH single-hole function.The double brackets in the above equations indicate the Hilbert-Schmidt norm.The operation takes place in Liouville space, rather than in the Hilbert space, taking the trace of the product of matrices rather than a vector product.
Given these definitions, the EOM for the coefficients and the SPDOs of the type I density matrix operator from Ref. 44 are where Ωτ is a Hartree product of SPDOs.

ρMCTDH(II)
The SPDOs of ρMCTDH(I) can further be expanded using the SPFs of the MCTDH method for wavefunctions σ , and SPFs {φ j }, leading to the ansatz This is known as the type II density operator.The EOM for the coefficients for the type II density matrices from Ref. 44 and the EOM for the SPFs for ρMCTDH(II) is where D is the single-particle reduced density matrix, where the aforementioned multi-index notation is used.
The extra effort of ρMCTDH(I) compared to wavepacket MCTDH is due to the extra dimension of the SPDOs compared to SPFs, and the bottleneck to this method lies in this high dimensionality.In contrast, as the type II density operator is expanded into SPFs, the bottleneck now lies in the expansion coefficients of these SPFs.ρMCTDH(I) will, therefore, be more efficient at high temperatures as the thermalization is automatically incorporated into the SPDOs rather than needing many configurations and coefficients as in ρMCTDH(II).
As a consequence of adhering to the Dirac-Frenkel variational principle, the MCTDH wavefunction conserves both total probability and energy.However, while this is still the case for ρMCTDH(II), for ρMCTDH(I), Tr(ρ 2 ) and Tr(ρ 2 H) are conserved rather than Tr(ρ) and Tr(ρH). 18The conservation of Tr(ρ) will be achieved in a converged calculation as an accurate description of the exact density matrix is obtained.

ML-ρMCTDH
A natural way to improve the scaling of ρMCTDH is to use the ML-MCTDH formalism for density matrices to enable larger calculations to be made.The ML-ρMCTDH method relates to ρMCTDH in the same way that ML-MCTDH relates to MCTDH.In essence, the SPDOs used to expand the full density matrix in the ρMCTDH equation [Eq.(32)] can be further expanded through the ML-algorithm.This new construction has analogous EOMs to the ML-MCTDH SPFs and coefficients, where the Hamiltonian operator is replaced with the Liouvillian operator.This novel technique can open the door to treating larger systems in a complete and accurate way, while also including temperature and solvent effects.
The derivation of the ML-ρMCTDH EOMS follows the same steps as ML-MCTDH.The ρMCTDH(I) density matrix is written in terms of SPDOs, (1) By using SPDOs on the top layer and each subsequent layer, the structure of the full ML expansion is thus the same as the standard ML formalism, but using a basis of density matrices rather than SPFs.Hence, using the Hilbert-Schmidt norm of ρMCTDH, the variational derivation above is valid with the following substitutions: The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpThus, it follows that the top layer equation [similar to (23)] is and the lower layer coefficients [as (28) for ML-MCTDH] where we have defined the following quantities for the layers: and the projector For type II, the density matrix is written in terms of SPFs rather than SPDOs, with configurations ϕ ′ K as for ML-MCTDH and a matrix of coefficients.Following the work of Raab et al. and extending to ML, the top layer coefficients of equations of motion come from (42) to give and the SPFs again come from (28) to give

III. COMPUTATIONAL DETAILS
The algorithms outlined above have all been implemented in the Quantics package 45,46 and, as shown below, are tested on a thermalized PT system.Initial thermalization of a density matrix can be achieved by propagating in imaginary time.This technique finds the lowest eigenvalue of a wavefunction and for wavepacket propagation is known as energy relaxation. 21,47e time-evolution of a wavepacket can be written as This equation can be converted to propagation in temperature by using the following limit to map from the time domain to the temperature domain it → 1 kT so that The same idea for a density matrix propagation results in a thermalized density matrix with the form This means that starting with a density matrix in which all states are equally populated (β = 0, i.e., T = ∞), propagation to temperature T yields the thermalized density.Systems of different sizes were propagated to examine the scaling properties of ρMCTDH(I), ρMCTDH(II), ML-ρMCTDH(I), and ML-ρMCTDH(II).Due to where the bottlenecks occur for the ρMCTDH(I) and ρMCTDH(II) methods, introduced in Sec.II B 2, ρMCTDH(I) much more readily lends itself to the ML formulation, as now the bottleneck lies at the bottom of the ML-tree, whereas the bottleneck with ML-ρMCTDH(II) lies at the top of the tree, where a large number of SPFs are needed.Therefore, the way that the ML-ρMCTDH(II) method is currently implemented means that it is a much less suitable method for larger systems.As a test, numerically exact propagations of the density matrix for the 2D system were also calculated to check that ρMCTDH was predicting the correct dynamics of the thermalized system.The dynamics of the proton transfer was extracted by evaluating the expectation value of a step operator placed at the proton transfer barrier along the reaction coordinate.This evaluates the wavepacket density that has crossed the barrier.
To provide a fair comparison of computational time for the different methods, all propagations used a standard Adams-Bashforth-Moulton (ABM) sixth-order adaptive time step predictor-corrector integration scheme.More efficient schemes have been developed for MCTDH wavepacket calculations, such as the constant mean-field (CMF) integrator, 21,48 which separates the different parts of the equations of motion and integrates them separately allowing longer time steps.However, no similar scheme is yet available for the ρMCTDH methods.

IV. TEST CASE I: PT IN SALICYLALDIMINE
Salicylaldimine exhibits a ground state intramolecular PT between an oxygen atom and a nitrogen atom, with the potential energy surface (PES) along the reaction coordinate with an asymmetric double well.The PES was computed in previous work as a fit to ab initio points at the Hartree-Fock theory level with a 6-31G * basis set. 33The PES function takes the form of a Taylor expansion of fourth order, and using mass-frequency scaled normal modes expanded around the transition state geometry, Q 0 .The thirteen modes that contribute the most to reaching the transition state from the minima were included.The energy barrier for conversion from the enol to the keto tautomer is 0.19 eV and for the reverse reaction is 0.058 eV. Figure 1 shows this potential surface along with the keto and enol tautomers.Details and parameters of the potential and all parameters can be found in the supplementary material of Ref. 33.
The focus of this test was the performance of the different methods with system size, and four systems of different sizes were studied.
From an analysis of the potential surface parameters, two modes have particular significance for the proton transfer.These are ν 1 , the proton transition mode, and, ν 36 , the in-plane perpendicular movement of the proton.This model is denoted as the 2D model.The next set of calculations used a larger 4D model.This added the modes ν 10 and ν 11 to the proton transfer modes ν 1 and ν 36 .A 7D model added modes ν 9 , ν 13 , and ν 23 .Finally, the full 13D model was studied.
Calculations used ρMCTDH on the 2D and 4D systems and ML-ρMCTDH on the 4D, 7D, and 13D models in both type I and type II formalisms where possible.Exact propagations of the density matrix were performed on the 2D system only.Dynamics simulations were initialized by relaxing to different temperatures.During this thermalization, the salicylaldimine was localized as the lower energy enol tautomer by using a harmonic potential where the minimum and frequency were chosen to match the position and curvature of the keto minimum of the double well.This potential had the form with q 0 = 1.790 75, ω 0 = 0.473 62 eV, and V 0 = −0.19246 eV.The harmonic oscillator potential energy surface can be seen in Fig. 1.Thus, a thermalized enol tautomer was created and subsequent dynamics could follow the PT process.

A. Scaling with system size
The first step is to demonstrate the scalings of the different methods by running propagations for 100 fs after relaxation to 1000 K.The zero point energy (ZPE) of the system along the reaction coordinate is 0.094 eV, which means that the energy barrier is effectively 1140 K, and thus, these simulations are just below the barrier.
The primitive basis for all calculations used a DVR.The ν 1 mode is represented as a sine DVR, with 61 grid points, and all the other modes use a harmonic oscillator basis on 21 grid points.

FIG. 1.
The ground state double well potential of the reaction coordinate for the proton transfer in salicylaldimine.The more stable tautomer is the enol form with the hydrogen on the oxygen.The left-hand plot is a cut through the PES along the PT mode, ν 1 .The red curve is a harmonic approximation to the enol minima.The right-hand plot shows the contours of the PES for this PT mode along with ν 36 , which is the motion of the proton perpendicular to ν 1 .

TABLE I.
A comparison of the timings and length of the vectors representing the density matrix, ρ, for the different methods for 100 fs propagations of different dimensional models of PT in salicylaldimine at 1000 K.

Timings
Length of ρ vector Convergence was achieved with respect to the property of interest: the PT dynamics.Sufficient convergence was taken to be the point where the dynamics no longer significantly changed by including more basis functions, i.e., adding additional SPFs did not change the PT fraction plots shown below.This is a non-trivial problem and required large numbers of basis functions for good results.The basis sets used for the ρMCTDH calculations and the ML-trees are given in the supplementary material.The groupings of the modes were based on the coupling strengths and grouping those of similar strength together.The cpu-times on a Xeon(R) CPU E5-2620 v3 compute node for 2D, 4D, 7D, and 13D calculations are shown in Table I.The first thing to note is the speed-up of the 2D calculation of the ρMCTDH calculations compared to the numerically exact calculation.The ρMCTDH(II) propagation is particularly efficient here, requiring only 35 s in comparison to 6 h taken by the full solution.However, due to the number of coefficients required by this method, ρMCTDH(II) was not able to treat the 7D system, which was still achievable using ρMCTDH(I).
The differing effort of the methods can be represented by the length of the vector required to represent the density matrices as for the same Hamiltonian, other parts of the computational effort (e.g., calculating integrals and mean-fields) will be either the same or dominated by the of the density matrix.For MCTDH, the length of the vector required to store a wavepacket, and hence the effort, can be written as where p is the number of particles, N is a representative value for the length of a one-dimensional primitive grid, n is a representative value for the number of SPFs for a particle, and d is the dimensionality of a particle.The first term is due to storing the expansion coefficients, and the second term is due to the SPFs.For ML-MCTDH, the length is a sum of terms of the form of Eq. (54) for each node of the multi-layer tree, with appropriate dimensions, noting that for the top and intermediate layers, N represents the number of SPFs in the layer below.It is thus highly dependent on the tree.For ρMCTDH(I), the effort is changed by squaring the primitive grid size due to the use of SPDOs, effort(pMCTDH(I)) ∼ n p + pnN 2d .( For ML-ρMCTDH(I), the form for the vector length for nodes in the top and intermediate layers is that of the MCTDH wavepacket, Eq. (54).Thus, the effort is only increased on the bottom layer by an effective doubling of degrees of freedom.Finally, for ρMCTDH(II), the effort is like that of MCTDH for the SPFs, but the size of the expansion vector is squared, For ML-ρMCTDH(II), the top layer has this size, and all other nodes have the form of Eq. (54).Thus, naturally type I should avoid combining degrees of freedom together on the lowest layer, while type II must keep the number of coefficients in the top layer low.It can also be seen from Table I that using the ML-ρMCTDH method does not automatically increase efficiency, when compared to the ρMCTDH method.Due to the additional layers of coupled calculations required to propagate the wavefunction in the multilayer form, for the smaller systems, it becomes more expensive to solve the wavefunction in this way.This is not specific to the ML-ρMCTDH method, but the ML-MCTDH method, in general, and the way the ML-tree is set up can greatly affect the efficiency of the ML calculations for both density matrices and pure wavefunctions.However, for the 7D system, the ML-ρMCTDH calculations are not only feasible, but much faster than ρMCTDH(I).The unfavorable scaling of ML-ρMCTDH(II) with system size meant that the 13D system could not be treated by this method with the tree used here due to the large number of top layer coefficients.A different tree may be more successful.This was the largest system investigated and was only accessible using ML-ρMCTDH(I).
An alternative method for simulating dynamics at finite temperatures is to use the thermal wavepacket approach. 25This stochastic approach takes the dynamics averaged over many randomly generated wavepackets.These wavefunctions are initially set to an infinite temperature, i.e., equal populations of all configurations but with random phases.This set of wavefunctions is then relaxed collectively to the correct temperature, and then, each of the resulting wavepackets is propagated simultaneously in time.Thermalizing the MCTDH wavefunction was achieved by setting the length of the relaxation, t final , to where kB is the Boltzmann constant and T is the desired temperature.A random number generator sets the relative phases of the configurations in the initial MCTDH wavefunction. 40Only after propagating all the generated wavefunctions and averaging out the results can the final dynamics of the system at that temperature be determined.To assess how accurately this thermal wavepacket method predicts dynamics when compared to the ρMCTDH approaches described above, the dynamics of the 4D system (defined by the PT fraction at 1000 K) was calculated and is compared in Fig. 2.
It can be seen that convergence was not reached in Fig. 2 using the thermal wavepacket method, despite 100 propagations FIG. 2. A comparison of the fraction of density that undergoes proton transfer for a 4D model of salicylaldimine at 1000 K for the thermal wavepacket method and the ρMCTDH method.The convergence of the thermal WP method is shown using 30, 50, 70, and 100 samples.
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp(and thus, the same number of relaxations) being carried out.Thus, although each ML-MCTDH propagation is quicker than the analogous ML-ρMCTDH calculation, taking only 50 cpu minutes, many more would be required for complete convergence.Indeed, this number greatly increases with temperature and system size.In comparison, while the density matrix approach is initially much more time consuming, only one relaxation and one propagation are necessary for the equivalent simulation, regardless of system size and temperature.Therefore, as temperature and system size increase, the ML-ρMCTDH method potentially becomes the more efficient, and more accurate, approach for simulating thermalized dynamics.

B. Temperature dependent dynamics
To get insights into the PT dynamics of salicylaldimine, the fraction of PT as a function of time was calculated over a range of temperatures by evaluating the expectation value of a step function placed at the transition barrier (Q = 0).This is shown for the 2D, 4D, and 13D systems in Fig. 3.In each case, the most efficient method, as discussed in Sec.IV A, was used, i.e., for 2D ρMCTDH(II), 4D ρMCTDH(I), and 13D ML-ρMCTDH(I).
In all cases, at 0 K, the only movement of the system occurs through proton tunneling through the barrier after being released from the harmonic well, allowing a small part of the wavepacket to cross from the global enol-minimum into the less stable ketominimum.This density oscillates with a period of ∼40 fs.For the 2D system, increasing the temperature from 0 → 500 K → 1000 K does not change the dynamics of the system significantly.There is, however, a notable increase in motion between the wells as the temperature is increased from 1000 K → 2000 K, which allows enough energy to cross the barrier.In the latter case, the oscillations remain, but the peak fraction of transfer increases from 20% at 1000 K to 25% at 2000 K and further to 30% at 3000 K.
The 4D system shows the effect of the increased zero point energy as there is an increase of PT at 1000 K compared to the 2D model.At higher temperatures still, the oscillations disappear and the density spreads across both wells.The PT dynamics of the 13D system is similar, except for a larger gap between the 1000 and 2000 K transfer fraction.

V. TEST CASE II: DOUBLE PT IN PORPHYCENE
As a further demonstration of the abilities of the ML-ρMCTDH method, a double proton transfer system was studied.Abdel-Latif and Kühn 34 provided a two-dimensional model Hamiltonian for the double-proton transfer in porphycene using symmetric and antisymmetric proton transfer coordinates and a potential form with four minima.The Hamiltonian is The parameters, taken from Ref. 34, are listed in the supplementary material.The potential is shown in Fig. 4(c).This Hamiltonian was extended to include the remaining 106 vibrations of porphycene in the form of a bath of harmonic oscillators coupled to the system modes, All modes are mass-frequency scaled.They are, therefore, defined as of Chemical Physics 0][51] This structure has D 2h symmetry, and the "symmetric" xs and "anti-symmetric" xa double-proton vibrations have imaginary frequencies with b 2u and b 3g symmetry, respectively.These vibrations are shown in Fig. 4.
Taking symmetry into account, as well as the coupling between the two system modes and the bath, the system-bath Hamiltonian is defined as where i are the 16 b 2u vibrations, j are the 17 b 3g vibrations, and k are the 17 a 1g vibrations.The final model Hamiltonian thus includes 55 modes.The parameters for the system-bath coupling were obtained by calculating the minimum energy structures and relating them to the potential at the high symmetry point.The values of the coupling parameters, along with more details of how they were obtained, are listed in supplementary material.
Next, the porphycene PT model Hamiltonian was used in a series of calculations using the ρMCTDH methods.Again, this consisted of different dimensional systems at a range of temperatures to demonstrate the scaling as well as the changes in physical behavior.In all cases, an initial energy relaxation calculation was made to thermalize the system to the desired temperature, followed by a propagation of 500 fs.The number of basis functions was chosen such that all natural populations were below 0.001 after 250 fs.That is, the initial dynamics is well represented.
As a benchmark for the system dynamics, initially, 2D calculations, including the xs and xa PT modes, were run at 0, 500, 1061, 1500, and 2000 K.For comparison, the barrier height is 1061 K.The initial relaxation localized the density in one well centered at (xs, xa) = (−0.607,0.0) by including a step function in the Hamiltonian during the relaxation to create a high wall at xs = 0.0.This wall was then removed for the propagation, and the expectation value of the step function as a function of time was used to determine the amount of proton transfer taking place.
The proton transfer at varying temperatures is shown in Fig. 5.At 0 K (purple line), there is a slow transfer, which is almost complete by 300 fs before returning to the initial well.This effect is likely to be due to tunneling.It can also be seen that as the temperature increases, the transfer is faster, and the amount transferred decreases.Snapshots of the density at different times are shown in Fig. 6.The density is taken from the trace of the time-evolving density matrix.At 0 K, the initially localized density is seen to move toward and then cross the barrier.At 1061 K, the initial density is again localized in the one well and vibrationally excited in the xa mode, as seen by the structure.In contrast to the 0 K case, this hot density flows quickly around the maximum of the potential at (xs, xa) = 0, 0 to undergo the double proton transfer.At 20 fs, it is seen to be occupying all four wells, after which the density oscillates back and forth between the two minima at (xs, xa = ±0.6,0).A 4D system was then studied, including the two bath modes with the strongest coupling ν 13 and ν 18 .The amount of proton transfer for this system as a function of time at different temperatures is shown in Fig. 5(b).The bath modes have a significant effect in slowing down the dynamics, and population transfer due to tunneling at 0 K is still increasing at 500 fs.Even at the barrier height temperature, the transfer is increasing in a similar way to the tunneling.Only when well above the barrier at 2000 K, does the transfer now take place in less than 100 fs, reaching 50% transfer and staying there, indicating that the density is spread across the minima equally.Finally, a 6D system adding the modes ν 17 and ν57 and a 10D system adding ν 12 , ν 20 , ν 8 , and ν 25 were examined to show the performance of the ML-ρMCTDH methods.The cpu-times for calculations on the four systems at 0 K using the different algorithms are listed in Table II.All four methods could be applied to the 4D system, and the efficiency of the ML-ρMCTDH algorithm is clear.This difference in effort is even more pronounced for the 6D case, for which regular ρMCTDH is not feasible.For the 10D system, only the ML-ρMCTDH(I) was feasible in the present implementation, and surprisingly fast, taking under 4 h for a 500 fs propagation.

VI. CONCLUSIONS
The results presented demonstrate not only that the ρMCTDH method can successfully and efficiently model the dynamics of a proton transfer in a ground state system at finite temperatures but also further that the ML-ρMCTDH formulation can treat larger density matrices than previously possible.This ability to predict accurate quantum dynamics on systems above 0 K, combined with solvent and/or environmental effects, has many possible applications.For small systems (2D), the standard ρMCTDH method is able to match the exact result from 0 to 3000 K.For larger systems, it was further shown that the new multilayer formulation was able to treat up to a 13D system over this range of temperatures, particularly with ML-ρMCTDH(I).Both single and double PT systems were investigated, and it was seen that at temperatures above the barrier, the dynamics changes from being dominated by tunneling through the barrier to direct transfer.
It is clear that future work to further improve the efficiency of the ML-ρMCTDH propagations would make this an even more powerful tool to tackle many interesting and important problems.One factor to consider is the choice of integrator.So far, a simple predictor-corrector integration scheme was used for the propagations.A more tailored integrator, such as the constant mean-field integrator, often used for the wavepacket MCTDH method, would lead to further savings in time.The new Hamiltonian for double PT coupled to a bath will also provide a suitable model for future tests of system-bath dynamics.Another factor to consider is the structure of the ML tree used for the ML-ρMCTDH calculations.The ML-MCTDH method, and therefore the ML-ρMCTDH methods, enormously depend on the how the ML trees are set up, and investigating further into the best structure would greatly increase the time efficiency of the ML-ρMCTDH calculations.

SUPPLEMENTARY MATERIAL
See the supplementary material for details on the basis sets and multi-layer trees used in the calculations and the parameters of the model Hamiltonian for porphycene and how they were obtained.

FIG. 3 .
FIG. 3. The fraction of density that undergoes proton transfer for a (a) 2D model, (b) 4D model, and (c) 13D model of salicylaldimine at a range of temperatures from 0 → 3000 K.

FIG. 4 .
FIG. 4. The (a) symmetric, xs, and (b) anti-symmetric, xa, proton transfer vibrational modes of porphycene, along with (c) the potential surface, calculated at the B3LYP/6-31+G * * level of theory.The minima for the cis and trans isomers are marked on the potential.

FIG. 5 .FIG. 6 .
FIG. 5. Transfer of density between PT minima in the porphycene model at different temperatures, starting in one minima.(a) 2D model and (b) 4D model.

TABLE II .
CPU time in seconds, and the length of the vectors representing the density matrices, ρ vectors, for various dimensional porphycene model systems propagated for 500 fs at 0 K using different algorithms.All calculations were run in serial on a Xeon(R) CPU E5-2620 v3 compute node.