Quantum and classical molecular dynamics for H atom scattering from graphene.

This work presents systematic comparisons between classical molecular dynamics (cMD) and quantum dynamics (QD) simulations of 15-dimensional and 75-dimensional models in their description of H atom scattering from graphene. We use an experimentally validated full-dimensional neural network potential energy surface of a hydrogen atom interacting with a large cell of graphene containing 24 carbon atoms. For quantum dynamics simulations, we apply Monte Carlo canonical polyadic decomposition to transform the original potential energy surface (PES) into a sum of products form and use the multi-layer multi-configuration time-dependent Hartree method to simulate the quantum scattering of a hydrogen or deuterium atom with an initial kinetic energy of 1.96 or 0.96 eV and an incident angle of 0°, i.e., perpendicular to the graphene surface. The cMD and QD initial conditions have been carefully chosen in order to be as close as possible. Our results show little differences between cMD and QD simulations when the incident energy of the H atom is equal to 1.96 eV. However, a large difference in sticking probability is observed when the incident energy of the H atom is equal to 0.96 eV, indicating the predominance of quantum effects. To the best of our knowledge, our work provides the first benchmark of quantum against classical simulations for a system of this size with a realistic PES. Additionally, new projectors are implemented in the Heidelberg multi-configuration time-dependent Hartree package for the calculation of the atom scattering energy transfer distribution as a function of outgoing angles.

Optimal Mode Combination in the Multiconfiguration Time-Dependent Hartree Method through Multivariate Statistics: Factor Analysis and Hierarchical Clustering.

The multiconfiguration time-dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH) are powerful algorithms for the efficient computation of nuclear quantum dynamics in high-dimensional systems. By providing time-dependent variational orbitals and an optimal choice of layered effective degrees of freedom, one is able to reduce the computational cost to an amenable number of configurations. However, choices related to selecting properly the mode grouping and tensor tree are strongly system dependent and, thus far, subjectively based on intuition and/or experience. Therefore, herein we detail a new protocol based on multivariate statistics─more specifically, factor analysis and hierarchical clustering─for a reliable and convenient guiding in the optimal design of such complex "system-of-systems" tensor-network decompositions. The advantages of employing the new algorithm and its applicability are tested on water and two floppy protonated water clusters with large amplitude motions.

The coupling of the hydrated proton to its first solvation shell.

The Zundel ([Formula: see text]) and Eigen ([Formula: see text]) cations play an important role as intermediate structures for proton transfer processes in liquid water. In the gas phase they exhibit radically different infrared (IR) spectra. The question arises: is there a least common denominator structure that explains the IR spectra of both, the Zundel and Eigen cations, and hence of the solvated proton? Full dimensional quantum simulations of these protonated cations demonstrate that two dynamical water molecules and an excess proton constitute this fundamental subunit. Embedded in the static environment of the parent Eigen cation, this subunit reproduces the positions and broadenings of its main excess-proton bands. In isolation, its spectrum reverts to the well-known Zundel ion. Hence, the dynamics of this subunit polarized by an environment suffice to explain the spectral signatures and anharmonic couplings of the solvated proton in its first solvation shell.

Representation of Diabatic Potential Energy Matrices for Multiconfiguration Time-Dependent Hartree Treatments of High-Dimensional Nonadiabatic Photodissociation Dynamics.

Conventional quantum mechanical characterization of photodissociation dynamics is restricted by steep scaling laws with respect to the dimensionality of the system. In this work, we examine the applicability of the multi-configurational time-dependent Hartree (MCTDH) method in treating nonadiabatic photodissociation dynamics in two prototypical systems, taking advantage of its favorable scaling laws. To conform to the sum-of-product form, elements of the ab initio diabatic potential energy matrix (DPEM) are re-expressed using the recently proposed Monte Carlo canonical polyadic decomposition method, with enforcement of proper symmetry. The MCTDH absorption spectra and product branching ratios are shown to compare well with those calculated using conventional grid-based methods, demonstrating its promise for treating high-dimensional nonadiabatic photodissociation problems.

Comparison of the multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) method and the density matrix renormalization group (DMRG) for ground state properties of linear rotor chains.

We demonstrate the applicability of the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method to the problem of computing ground states of one-dimensional chains of linear rotors with dipolar interactions. Specifically, we successfully obtain energies, entanglement entropies, and orientational correlations that are in agreement with the Density Matrix Renormalization Group (DMRG), which has been previously used for this system. We find that the entropies calculated by ML-MCTDH for larger system sizes contain nonmonotonicity, as expected in the vicinity of a second-order quantum phase transition between ordered and disordered rotor states. We observe that this effect remains when all couplings besides nearest-neighbor are omitted from the Hamiltonian, which suggests that it is not sensitive to the rate of decay of the interactions. In contrast to DMRG, which is tailored to the one-dimensional case, ML-MCTDH (as implemented in the Heidelberg MCTDH package) requires more computational time and memory, although the requirements are still within reach of commodity hardware. The numerical convergence and computational demand of two practical implementations of ML-MCTDH and DMRG are presented in detail for various combinations of system parameters.

High-Dimensional Quantum Dynamics Study on Excitation-Specific Surface Scattering Including Lattice Effects of a Five-Atom Surface Cell.

In this work, high-dimensional (21D) quantum dynamics calculations on the mode-specific surface scattering of a carbon monoxide molecule on a copper(100) surface with lattice effects of a five-atom surface cell are performed through the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method. We employ a surface model in which five surface atoms near the impact site are treated as fully flexible quantum particles, while all other more distant atoms are kept at fixed locations. To efficiently perform the 21D ML-MCTDH wave packet propagation, the potential energy surface is transferred to a canonical polyadic decomposition form with the aid of a Monte Carlo-based method. Excitation-specific sticking probabilities of CO on Cu(100) are computed, and lattice effects caused by the flexible surface atoms are demonstrated by comparison with sticking probabilities computed for a rigid surface. The dependence of the sticking probability of the initial state of the system is studied, and it is found that the sticking probability is reduced when the surface atom on the impact site is initially vibrationally excited.

Importance of Appropriately Regularizing the ML-MCTDH Equations of Motion.

The multiconfiguration time-dependent Hartree (MCTDH) method and its generalization, the multilayer MCTDH (ML-MCTDH), result in equations of motion (EOMs) that are singular when there are virtual orbitals-the unoccupied single-particle functions-in the wave function expansion. For decades this singularity had been numerically removed by regularizing the reduced density matrix. In this Perspective we discuss our recent proposal to regularize the coefficient tensor instead, which has significant impact on both the efficiency and correctness of the EOMs in MCTDH and ML-MCTDH for challenging problems. We further demonstrate that when the system becomes large such that it is necessary to use ML-MCTDH with many layers, it is much more important to employ this new regularization scheme. We illustrate this point by studying a spin-boson model with a large bath that contains up to 100 000 modes. We show that even in the weak coupling regime the new regularization scheme is required to quickly rotate the virtual orbitals into the correct directions in Hilbert space. We argue that this situation can be common for applying a time-dependent tensor network approach to any large enough system.

Regularizing the MCTDH equations of motion through an optimal choice on-the-fly (i.e., spawning) of unoccupied single-particle functions.

The multi-configuration time-dependent Hartree method is a general algorithm to solve the time-dependent Schrödinger equation, in which the wavefunction is expanded in a direct product of self-adapting time-dependent Single-Particle Functions (SPFs) that are propagated in time according to the Dirac-Frenkel variational principle. In the current version of this approach, the size of the SPF basis is fixed at the outset so that singularities in the working equations resulting from unoccupied functions have to be removed by a regularization procedure. Here, an alternative protocol is presented, in which we gradually increase the number of unoccupied SPFs on-the-fly (i.e., spawning) and optimize their shape by variationally minimizing the error made by the finite size of the basis. An initial estimate for the respective new expansion coefficients is also computed, thus avoiding the need to regularize the equations of motion. The advantages of employing the new algorithm are tested and discussed in some illustrative examples.

Theoretical investigation of the H + HD → D + H2 chemical reaction for astrophysical applications: A state-to-state quasi-classical study.

We report a large set of state-to-state rate constants for the H + HD reactive collision, using Quasi-Classical Trajectory (QCT) simulations on the accurate H3 global potential energy surface of Mielke et al. [J. Chem. Phys. 116, 4142 (2002)]. High relative collision energies (up to ≈56 000 K) and high rovibrational levels of HD (up to ≈50 000 K), relevant to various non thermal equilibrium astrophysical media, are considered. We have validated the accuracy of our QCT calculations with a new efficient adaptation of the Multi Configuration Time Dependent Hartree (MCTDH) method to compute the reaction probability of a specific reactive channel. Our study has revealed that the high temperature regime favors the production of H2 in its highly rovibrationnally excited states, which can de-excite radiatively (cooling the gas) or collisionally (heating the gas). Those new state-to-state QCT reaction rate constants represent a significant improvement in our understanding of the possible mechanisms leading to the destruction of HD by its collision with a H atom.

On regularizing the ML-MCTDH equations of motion.

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.

On regularizing the MCTDH equations of motion.

The Multiconfiguration Time-Dependent Hartree (MCTDH) approach leads to equations of motion (EOM) which become singular when there are unoccupied so-called single-particle functions (SPFs). Starting from a Hartree product, all SPFs, except the first one, are unoccupied initially. To solve the MCTDH-EOMs numerically, one therefore has to remove the singularity by a regularization procedure. Usually the inverse of a density matrix is regularized. Here we argue and show that regularizing the coefficient tensor, which in turn regularizes the density matrix as well, leads to an improved performance of the EOMs. The initially unoccupied SPFs are rotated faster into their "correct direction" in Hilbert space and the final results are less sensitive to the choice of the value of the regularization parameter. For a particular example (a spin-boson system studied with a transformed Hamiltonian), we could even show that only with the new regularization scheme could one obtain correct results. Finally, in Appendix A, a new integration scheme for the MCTDH-EOMs developed by Lubich and co-workers is discussed. It is argued that this scheme does not solve the problem of the unoccupied natural orbitals because this scheme ignores the latter and does not propagate them at all.

Multidimensional Quantum Mechanical Modeling of Electron Transfer and Electronic Coherence in Plant Cryptochromes: The Role of Initial Bath Conditions.

A multidimensional quantum mechanical protocol is used to describe the photoinduced electron transfer and electronic coherence in plant cryptochromes without any semiempirical, e.g., experimentally obtained, parameters. Starting from a two-level spin-boson Hamiltonian we look at the effect that the initial photoinduced nuclear bath distribution has on an intermediate step of this biological electron transfer cascade for two idealized cases. The first assumes a slow equilibration of the nuclear bath with respect to the previous electron transfer step that leads to an ultrafast decay with little temperature dependence; while the second assumes a prior fast bath equilibration on the donor potential energy surface leading to a much slower decay, which contrarily displays a high temperature dependence and a better agreement with previous theoretical and experimental results. Beyond Marcus and semiclassical pictures these results unravel the strong impact that the presence or not of equilibrium initial conditions has on the electronic population and coherence dynamics at the quantum dynamics level in this and conceivably in other biological electron transfer cascades.

Transforming high-dimensional potential energy surfaces into sum-of-products form using Monte Carlo methods.

We propose a Monte Carlo method, "Monte Carlo Potfit," for transforming high-dimensional potential energy surfaces evaluated on discrete grid points into a sum-of-products form, more precisely into a Tucker form. To this end we use a variational ansatz in which we replace numerically exact integrals with Monte Carlo integrals. This largely reduces the numerical cost by avoiding the evaluation of the potential on all grid points and allows a treatment of surfaces up to 15-18 degrees of freedom. We furthermore show that the error made with this ansatz can be controlled and vanishes in certain limits. We present calculations on the potential of HFCO to demonstrate the features of the algorithm. To demonstrate the power of the method, we transformed a 15D potential of the protonated water dimer (Zundel cation) in a sum-of-products form and calculated the ground and lowest 26 vibrationally excited states of the Zundel cation with the multi-configuration time-dependent Hartree method.

Rotational effects on the dissociation dynamics of CHD3 on Pt(111).

Dissociation of methane on metal surfaces is of high practical and fundamental interest. Therefore there is currently a big push aimed at determining the simplest dynamical model that allows the reaction dynamics to be described with quantitative accuracy using quantum dynamics. Using five-dimensional quantum dynamical and full-dimensional ab initio molecular dynamics calculations, we show that the CD3 umbrella axis of CHD3 must reorient before the molecule reaches the barrier for C-H cleavage to occur in reaction on Pt(111). This rules out the application of the rotationally sudden approximation, as explicitly shown through a comparison with calculations using this approximation. Further, we suggest that the observed umbrella swing should strongly affect the sensitivity of C-H cleavage to the initial alignment of the molecule relative to the surface as found experimentally for closely related systems. We find very large differences in reactivity for molecules pre-excited to different rotational states, particularly if these states are associated with different orientations of the C-H bond.

Isotope effects of ground and lowest lying vibrational states of H(3-x)D(x)O2(-) complexes.

Isotope effects of the H3O2(-) anion are investigated. For this, the 24 lowest excited vibrational states of the H(3-x)D(x)O2(-) complexes, with x = 0-3, are computed using two different Hamiltonians, namely, a 7D reduced-dimensionality one with a numerical representation of the kinetic energy operator (KEO) and a 9D full-dimensionality Hamiltonian with an exact analytic KEO. The computations are carried out with the multiconfiguration time-dependent Hartree method. The obtained results show that bridge and terminal H-D exchange cause a variation in energy with or without a rearrangement of states. A clear rearrangement of fundamental modes is observed in bridge H-D exchange of the H3O2(-) complex, where the frequency of bridge hydrogen stretching (z) is strongly lowered by substitution. The isotope effects show that rotation (ϕ), rocking (u1 + u2), wagging (u1 - u2), and O-O stretch (R) modes are sensitive to terminal H-D exchange, while the bridge-atom bending (x,y) and stretch (z) modes are sensitive to bridge H-D exchange. An influence coefficient, which measures the influence of an excitation of one mode on the various 1D reduced densities, is defined and analyzed in detail. It is shown that the D3O2(-) complex is more strongly correlated or coupled than the other isotopologues.

Expansion Hamiltonian model for a diatomic molecule adsorbed on a surface: Vibrational states of the CO/Cu(100) system including surface vibrations.

Molecular-surface studies are often done by assuming a corrugated, static (i.e., rigid) surface. To be able to investigate the effects that vibrations of surface atoms may have on spectra and cross sections, an expansion Hamiltonian model is proposed on the basis of the recently reported [R. Marquardt et al., J. Chem. Phys. 132, 074108 (2010)] SAP potential energy surface (PES), which was built for the CO/Cu(100) system with a rigid surface. In contrast to other molecule-surface coupling models, such as the modified surface oscillator model, the coupling between the adsorbed molecule and the surface atoms is already included in the present expansion SAP-PES model, in which a Taylor expansion around the equilibrium positions of the surface atoms is performed. To test the quality of the Taylor expansion, a direct model, that is avoiding the expansion, is also studied. The latter, however, requests that there is only one movable surface atom included. On the basis of the present expansion and direct models, the effects of a moving top copper atom (the one to which CO is bound) on the energy levels of a bound CO/Cu(100) system are studied. For this purpose, the multiconfiguration time-dependent Hartree calculations are carried out to obtain the vibrational fundamentals and overtones of the CO/Cu(100) system including a movable top copper atom. In order to interpret the results, a simple model consisting of two coupled harmonic oscillators is introduced. From these calculations, the vibrational levels of the CO/Cu(100) system as function of the frequency of the top copper atom are discussed.

Resonances of HCO Computed Using an Approach Based on the Multiconfiguration Time-Dependent Hartree Method.

The improved relaxation method with a complex absorbing potential (CAP) was used to compute resonance states of the formyl radical (HCO) using the Heidelberg multi-configuration time-dependent Hartree (MCTDH) program. To benchmark this approach, the same potential energy surface as was used in three other method development studies was used here. It was found that the MCTDH-based approach was able to accurately and efficiently compute 90 resonance states up to more than 1 eV above the dissociation limit. Extremely close agreement was obtained for energies and widths (lifetimes) calculated using MCTDH compared with those reported previously for three other CAP-based approaches that separately involved filter-diagonalization, a preconditioned complex-symmetric Lanczos algorithm, and a non-Hermitian real-arithmatic Lanczos method. The high accuracy achieved in this benchmark study supports the applicability of MCTDH to the study of resonances in larger systems in which increased dimensionality makes the efficiency of MCTDH advantageous.

A full-dimensional multilayer multiconfiguration time-dependent Hartree study on the ultraviolet absorption spectrum of formaldehyde oxide.

Employing the multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method in conjunction with the multistate multimode vibronic coupling Hamiltonian (MMVCH) model, we perform a full dimensional (9D) quantum dynamical study on the simplest Criegee intermediate, formaldehyde oxide, in five lower-lying singlet electronic states. The ultraviolet (UV) spectrum is then simulated by a Fourier transform of the auto-correlation function. The MMVCH model is built based on extensive MRCI(8e,8o)/aug-cc-pVTZ calculations. To ensure a fast convergence of the final calculations, a large number of ML-MCTDH test calculations is performed to find an appropriate multilayer separations (ML-trees) of the ML-MCTDH nuclear wave functions, and the dynamical calculations are carefully checked to ensure that the calculations are well converged. To compare the computational efficiency, standard MCTDH simulations using the same Hamiltonian are also performed. A comparison of the MCTDH and ML-MCTDH calculations shows that even for the present not-too-large system (9D here) the ML-MCTDH calculations can save a considerable amount of computational resources while producing identical spectra as the MCTDH calculations. Furthermore, the present theoretical B̃(1)A'←X̃(1)A' UV spectral band and the corresponding experimental measurements [J. M. Beames, F. Liu, L. Lu, and M. I. Lester, J. Am. Chem. Soc. 134, 20045-20048 (2012); L. Sheps, J. Phys. Chem. Lett. 4, 4201-4205 (2013); W.-L. Ting, Y.-H. Chen, W. Chao, M. C. Smith, and J. J.-M. Lin, Phys. Chem. Chem. Phys. 16, 10438-10443 (2014)] are discussed. To the best of our knowledge, this is the first theoretical UV spectrum simulated for this molecule including nuclear motion beyond an adiabatic harmonic approximation.

Rovibrational spectroscopy using a kinetic energy operator in Eckart frame and the multi-configuration time-dependent Hartree (MCTDH) approach.

For computational rovibrational spectroscopy the choice of the frame is critical for an approximate separation of overall rotation from internal motions. To minimize the coupling between internal coordinates and rotation, Eckart proposed a condition ["Some studies concerning rotating axes and polyatomic molecules," Phys. Rev. 47, 552-558 (1935)] and a frame that fulfills this condition is hence called an Eckart frame. A method is developed to introduce in a systematic way the Eckart frame for the expression of the kinetic energy operator (KEO) in the polyspherical approach. The computed energy levels of a water molecule are compared with those obtained using a KEO in the standard definition of the Body-fixed frame of the polyspherical approach. The KEO in the Eckart frame leads to a faster convergence especially for large J states and vibrationally excited states. To provide an example with more degrees of freedom, rotational states of the vibrational ground state of the trans nitrous acid (HONO) are also investigated.

Calculation of the vibrational excited states of malonaldehyde and their tunneling splittings with the multi-configuration time-dependent Hartree method.

We report energies and tunneling splittings of vibrational excited states of malonaldehyde which have been obtained using full dimensional quantum mechanical calculations. To this end we employed the multi configuration time-dependent Hartree method. The results have been obtained using a recently published potential energy surface [Y. Wang, B. J. Braams, J. M. Bowman, S. Carter, and D. P. Tew, J. Chem. Phys. 128, 224314 (2008)] which has been brought into a suitable form by a modified version of the n-mode representation which was used with two different arrangements of coordinates. The relevant terms of the expansion have been identified with a Metropolis algorithm and a diffusion Monte-Carlo technique, respectively.