A bivariational, stable, and convergent hierarchy for time-dependent coupled cluster with adaptive basis sets.

We propose a new formulation of time-dependent coupled cluster with adaptive basis functions and division of the one-particle space into active and secondary subspaces. The formalism is fully bivariational in the sense of a real-valued time-dependent bivariational principle and converges to the complete-active-space solution, a property that is obtained by the use of biorthogonal basis functions. A key and distinguishing feature of the theory is that the active bra and ket functions span the same space by construction. This ensures numerical stability and is achieved by employing a split unitary/non-unitary basis set transformation: the unitary part changes the active space itself, while the non-unitary part transforms the active basis. The formulation covers vibrational as well as electron dynamics. Detailed equations of motion are derived and implemented in the context of vibrational dynamics, and the numerical behavior is studied and compared to related methods.

Vibrational ADAPT-VQE: Critical points lead to problematic convergence.

Quantum chemistry is one of the most promising applications for which quantum computing is expected to have a significant impact. Despite considerable research in the field of electronic structure, calculating the vibrational properties of molecules on quantum computers remains a relatively unexplored field. In this work, we develop a vibrational Adaptive Derivative-Assembled Pseudo-Trotter Variational Quantum Eigensolver (vADAPT-VQE) formalism based on an infinite product representation (IPR) of anti-Hermitian excitation operators of the Full Vibrational Configuration Interaction (FVCI) wavefunction, which allows for preparing eigenstates of vibrational Hamiltonians on quantum computers. In order to establish the vADAPT-VQE algorithm using the IPR, we study the exactness of disentangled Unitary Vibrational Coupled Cluster (dUVCC) theory and show that dUVCC can formally represent the FVCI wavefunction in an infinite expansion. To investigate the performance of the vADAPT-VQE algorithm, we numerically study whether the vADAPT-VQE algorithm generates a sequence of operators that may represent the FVCI wavefunction. Our numerical results indicate frequent appearance of critical points in the wavefunction preparation using vADAPT-VQE. These results imply that one may encounter diminishing usefulness when preparing vibrational wavefunctions on quantum computers using vADAPT-VQE and that additional studies are required to find methods that can circumvent this behavior.

Efficient vibrationally correlated calculations using n-mode expansion-based kinetic energy operators.

Due to its efficiency and flexibility, the n-mode expansion is a frequently used tool for representing molecular potential energy surfaces in quantum chemical simulations. In this work, we investigate the performance of n-mode expansion-based models of kinetic energy operators in general polyspherical coordinate systems. In particular, we assess the operators with respect to accuracy in vibrationally correlated calculations and their effect on potential energy surface construction with the adaptive density guided approach. Our results show that the n-mode expansion-based operator variants are reliable and systematically improvable approximations of the full kinetic energy operator. Moreover, we introduce a workflow to generate the n-mode expanded kinetic energy operators on-the-fly within the adaptive density guided approach. This scheme can be applied in studies of species and coordinate systems, for which an analytical form of the kinetic energy operator is not available.

Quasi-direct Quantum Molecular Dynamics: The Time-Dependent Adaptive Density-Guided Approach for Potential Energy Surface Construction.

We present a new quasi-direct quantum molecular dynamics computational method which offers a compromise between quantum dynamics using a precomputed potential energy surface (PES) and fully direct quantum dynamics. This method is termed the time-dependent adaptive density-guided approach (TD-ADGA) and is a method for constructing a PES on the fly during a dynamics simulation. This is achieved by acquisition of new single-point (SP) calculations and refitting of the PES, depending on the need of the dynamics. The TD-ADGA is a further development of the adaptive density-guided approach (ADGA) for PES construction where the placement of SPs is guided by the density of the nuclear wave function. In TD-ADGA, the ADGA framework has been integrated into the time propagation of the time-dependent nuclear wave function and we use the reduced one-mode density of this wave function to guide when and where new SPs are placed. The PES is thus extended or updated if the wave function moves into new areas or if a certain area becomes more important. Here, we derive equations for the reduced one-mode density for the time-dependent Hartree (TDH) method and for multiconfiguration time-dependent Hartree (MCTDH) methods, but the TD-ADGA can be used with any time-dependent wave function method as long as a density is available. The TD-ADGA method has been investigated on molecular systems containing single- and double-minimum potentials and on single-mode and multi-mode systems. We explore different approaches to handle the fact that the TD-ADGA involves a PES that changes during the computation and show how results can be obtained that are in very good agreement with results obtained by using an accurate reference PES. Dynamics with TD-ADGA is essentially a black box procedure, where only the initialization of the system and how to compute SPs must be provided. The TD-ADGA thus makes it easier to carry out quantum molecular dynamics and the quasi-direct framework opens up the possibility to compute quantum dynamics accurately for larger molecular systems.

Time-dependent coupled cluster with orthogonal adaptive basis functions: General formalism and application to the vibrational problem.

We derive equations of motion for bivariational wave functions with orthogonal adaptive basis sets and specialize the formalism to the coupled cluster Ansatz. The equations are related to the biorthogonal case in a transparent way, and similarities and differences are analyzed. We show that the amplitude equations are identical in the orthogonal and biorthogonal formalisms, while the linear equations that determine the basis set time evolution differ by symmetrization. Applying the orthogonal framework to the nuclear dynamics problem, we introduce and implement the orthogonal time-dependent modal vibrational coupled cluster (oTDMVCC) method and benchmark it against exact reference results for four triatomic molecules as well as a reduced-dimensional (5D) trans-bithiophene model. We confirm numerically that the biorthogonal TDMVCC hierarchy converges to the exact solution, while oTDMVCC does not. The differences between TDMVCC and oTDMVCC are found to be small for three of the five cases, but we also identify one case where the formal deficiency of the oTDMVCC approach results in clear and visible errors relative to the exact result. For the remaining example, oTDMVCC exhibits rather modest but visible errors.

Efficient time-dependent vibrational coupled cluster computations with time-dependent basis sets at the two-mode coupling level: Full and hybrid TDMVCC[2].

The computation of the nuclear quantum dynamics of molecules is challenging, requiring both accuracy and efficiency to be applicable to systems of interest. Recently, theories have been developed for employing time-dependent basis functions (denoted modals) with vibrational coupled cluster theory (TDMVCC). The TDMVCC method was introduced along with a pilot implementation, which illustrated good accuracy in benchmark computations. In this paper, we report an efficient implementation of TDMVCC, covering the case where the wave function and Hamiltonian contain up to two-mode couplings. After a careful regrouping of terms, the wave function can be propagated with a cubic computational scaling with respect to the number of degrees of freedom. We discuss the use of a restricted set of active one-mode basis functions for each mode, as well as two interesting limits: (i) the use of a full active basis where the variational modal determination amounts essentially to the variational determination of a time-dependent reference state for the cluster expansion; and (ii) the use of a single function as an active basis for some degrees of freedom. The latter case defines a hybrid TDMVCC/TDH (time-dependent Hartree) approach that can obtain even lower computational scaling. The resulting computational scaling for hybrid and full TDMVCC[2] is illustrated for polyaromatic hydrocarbons with up to 264 modes. Finally, computations on the internal vibrational redistribution of benzoic acid (39 modes) are used to show the faster convergence of TDMVCC/TDH hybrid computations towards TDMVCC compared to simple neglect of some degrees of freedom.

Optimizing the number of measurements for vibrational structure on quantum computers: coordinates and measurement schemes.

One of the primary challenges prohibiting demonstrations of practical quantum advantages for near-term devices amounts to excessive measurement overheads for estimating relevant physical quantities such as ground state energies. However, with major differences between the electronic and vibrational structures of molecules, the question of how the resource requirements of computing anharmonic, vibrational states can be reduced remains relatively unexplored compared to its electronic counterpart. Importantly, bosonic commutation relations, distinguishable Hilbert spaces and vibrational coordinates allow manipulations of the vibrational system that can be exploited to minimize resource requirements. In this work, we investigate the impact of different coordinate systems and measurement schemes on the number of measurements needed to estimate anharmonic, vibrational states for a variety of three-mode (six-mode) molecules. We demonstrate an average of 3-fold (1.5-fold), with up to 7-fold (2.5-fold), reduction in the number of measurements required by employing appropriate coordinate transformations, based on an automized construction of qubit Hamiltonians from a conventional vibrational structure program.

Gaussian process regression adaptive density-guided approach: Toward calculations of potential energy surfaces for larger molecules.

We present a new program implementation of the Gaussian process regression adaptive density-guided approach [Schmitz et al., J. Chem. Phys. 153, 064105 (2020)] for automatic and cost-efficient potential energy surface construction in the MidasCpp program. A number of technical and methodological improvements made allowed us to extend this approach toward calculations of larger molecular systems than those previously accessible and maintain the very high accuracy of constructed potential energy surfaces. On the methodological side, improvements were made by using a Δ-learning approach, predicting the difference against a fully harmonic potential, and employing a computationally more efficient hyperparameter optimization procedure. We demonstrate the performance of this method on a test set of molecules of growing size and show that up to 80% of single point calculations could be avoided, introducing a root mean square deviation in fundamental excitations of about 3 cm-1. A much higher accuracy with errors below 1 cm-1 could be achieved with tighter convergence thresholds still reducing the number of single point computations by up to 68%. We further support our findings with a detailed analysis of wall times measured while employing different electronic structure methods. Our results demonstrate that GPR-ADGA is an effective tool, which could be applied for cost-efficient calculations of potential energy surfaces suitable for highly accurate vibrational spectra simulations.

General exponential basis set parametrization: Application to time-dependent bivariational wave functions.

We present equations of motion (EOMs) for general time-dependent wave functions with exponentially parameterized biorthogonal basis sets. The equations are fully bivariational in the sense of the time-dependent bivariational principle and offer an alternative, constraint-free formulation of adaptive basis sets for bivariational wave functions. We simplify the highly non-linear basis set equations using Lie algebraic techniques and show that the computationally intensive parts of the theory are, in fact, identical to those that arise with linearly parameterized basis sets. Thus, our approach offers easy implementation on top of existing code in the context of both nuclear dynamics and time-dependent electronic structure. Computationally tractable working equations are provided for single and double exponential parametrizations of the basis set evolution. The EOMs are generally applicable for any value of the basis set parameters, unlike the approach of transforming the parameters to zero at each evaluation of the EOMs. We show that the basis set equations contain a well-defined set of singularities, which are identified and removed by a simple scheme. The exponential basis set equations are implemented in conjunction with the time-dependent modals vibrational coupled cluster (TDMVCC) method, and we investigate the propagation properties in terms of the average integrator step size. For the systems we test, the exponentially parameterized basis sets yield slightly larger step sizes compared to the linearly parameterized basis set.

Current and future machine learning approaches for modeling atmospheric cluster formation.

The formation of strongly bound atmospheric molecular clusters is the first step towards forming new aerosol particles. Recent advances in the application of machine learning models open an enormous opportunity for complementing expensive quantum chemical calculations with efficient machine learning predictions. In this Perspective, we present how data-driven approaches can be applied to accelerate cluster configurational sampling, thereby greatly increasing the number of chemically relevant systems that can be covered.

Bivariational time-dependent wave functions with biorthogonal adaptive basis sets: General formulation and regularization of equations of motion through polar decomposition.

We derive general bivariational equations of motion (EOMs) for time-dependent wave functions with biorthogonal time-dependent basis sets. The time-dependent basis functions are linearly parameterized and their fully variational time evolution is ensured by solving a set of so-called constraint equations, which we derive for arbitrary wave function expansions. The formalism allows division of the basis set into an active basis and a secondary basis, ensuring a flexible and compact wave function. We show how the EOMs specialize to a few common wave function forms, including coupled cluster and linearly expanded wave functions. It is demonstrated, for the first time, that the propagation of such wave functions is not unconditionally stable when a secondary basis is employed. The main signature of the instability is a strong increase in non-orthogonality, which eventually causes the calculation to fail; specifically, the biorthogonal active bra and ket bases tend toward spanning different spaces. Although formally allowed, this causes severe numerical issues. We identify the source of this problem by reparametrizing the time-dependent basis set through polar decomposition. Subsequent analysis allows us to remove the instability by setting appropriate matrix elements to zero. Although this solution is not fully variational, we find essentially no deviation in terms of autocorrelation functions relative to the variational formulation. We expect that the results presented here will be useful for the formal analysis of bivariational time-dependent wave functions for electronic and nuclear dynamics in general and for the practical implementation of time-dependent CC wave functions in particular.

Massive Assessment of the Binding Energies of Atmospheric Molecular Clusters.

Quantum chemical studies of the formation and growth of atmospheric molecular clusters are important for understanding aerosol particle formation. However, the search for the lowest free-energy cluster configuration is extremely time consuming. This makes high-level benchmark data sets extremely valuable in the quest for the global minimum as it allows the identification of cost-efficient computational methodologies, as well as the development of high-level machine learning (ML) models. Herein, we present a highly versatile quantum chemical data set comprising a total of 11 749 (acid)1-2(base)1-2 cluster configurations, containing up to 44 atoms. Utilizing the LUMI supercomputer, we calculated highly accurate PNO-CCSD(F12*)(T)/cc-pVDZ-F12 binding energies of the full set of cluster configurations leading to an unprecedented data set both in regard to sheer size and with respect to the level of theory. We employ the constructed benchmark set to assess the performance of various semiempirical and density functional theory methods. In particular, we find that the r2-SCAN-3c method shows excellent performance across the data set related to both accuracy and CPU time, making it a promising method to employ during cluster configurational sampling. Furthermore, applying the data sets, we construct ML models based on Δ-learning and provide recommendations for future application of ML in cluster configurational sampling.

Calculating vibrational excitation energies using tensor-decomposed vibrational coupled-cluster response theory.

The first implementation of tensor-decomposed vibrational coupled cluster (CP-VCC) response theory for calculating vibrational excitation energies is presented. The CP-VCC algorithm, which has previously been applied to solving the vibrational coupled cluster (VCC) ground-state equations without explicitly constructing any tensors of order three or higher, has been generalized to allow transformations with the Jacobian matrix necessary for computation of response excitation energies by iterative algorithms. A new eigenvalue solver for computing CP-VCC excitation energies is introduced, and the different numerical thresholds used for controlling the accuracy of the obtained eigenvalues are discussed. Numerical results are presented for calculations of the 20 lowest eigenvalues on a set of 10 four-atomic molecules, as well as for a number of polycyclic aromatic hydrocarbons (PAHs) of increasing size, up to PAH8 with 120 modes. It is shown that the errors introduced by the tensor decomposition can be controlled by the choice of numerical thresholds. Furthermore, all thresholds can be defined relative to the requested convergence threshold of the equation solver, which allows black-box calculations with minimal user input to be performed. Eigenstates of PAHs were efficiently computed without any explicitly constructed tensors, showing improvements in both memory and central processing unit time compared to the existing full-tensor versions.

A general implementation of time-dependent vibrational coupled-cluster theory.

The first general excitation level implementation of the time-dependent vibrational coupled cluster (TDVCC) method introduced in a recent publication [J. Chem. Phys. 151, 154116 (2019)] is presented. The general framework developed for time-independent vibrational coupled cluster (VCC) calculations has been extended to the time-dependent context. This results in an efficient implementation of TDVCC with general coupling levels in the cluster operator and Hamiltonian. Thus, the convergence of the TDVCC[k] hierarchy toward the complete-space limit can be studied for any sum-of-product Hamiltonian. Furthermore, a scheme for including selected higher-order excitations for a subset of modes is introduced and studied numerically. Three different definitions of the TDVCC autocorrelation function (ACF) are introduced and analyzed in both theory and numerical experiments. Example calculations are presented for an array of systems including imidazole, formyl fluoride, formaldehyde, and a reduced-dimensionality bithiophene model. The results show that the TDVCC[k] hierarchy converges systematically toward the full-TDVCC limit and that the implementation allows accurate quantum-dynamics simulations of large systems to be performed. Specifically, the intramolecular vibrational-energy redistribution of the 21-dimensional imidazole molecule is studied in terms of the decay of the ACF. Furthermore, the importance of product separability in the definition of the ACF is highlighted when studying non-interacting subsystems.

Vibrationally resolved coupled-cluster x-ray absorption spectra from vibrational configuration interaction anharmonic calculations.

Vibrationally resolved near-edge x-ray absorption spectra at the K-edge for a number of small molecules have been computed from anharmonic vibrational configuration interaction calculations of the Franck-Condon factors. The potential energy surfaces for ground and core-excited states were obtained at the core-valence separated CC2, CCSD, CCSDR(3), and CC3 levels of theory, employing the adaptive density-guided approach scheme to select the single points at which to perform the energy calculations. We put forward an initial attempt to include pair-mode coupling terms to describe the potential of polyatomic molecules.

Time-dependent vibrational coupled cluster with variationally optimized time-dependent basis sets.

We develop time-dependent vibrational coupled cluster with time-dependent modals (TDMVCC), where an active set of one-mode basis functions (modals) is evolved in time alongside coupled-cluster wave-function parameters. A biorthogonal second quantization formulation of many-mode dynamics is introduced, allowing separate biorthogonal bases for the bra and ket states, thus ensuring complex analyticity. We employ the time-dependent bivariational principle to derive equations of motion for both the one-mode basis functions and the parameters describing the cluster (T) and linear de-excitation (L) operators. The choice of constraint (or gauge) operators for the modal time evolution is discussed. In the case of untruncated cluster expansion, the result is independent of this choice, but restricting the excitation space removes this invariance; equations for the variational determination of the constraint operators are derived for the latter case. We show that all single-excitation parts of T and L are redundant and can be left out in the case of variationally determined constraint-operator evolution. Based on a pilot implementation, test computations on Henon-Heiles model systems, the water molecule, and a reduced-dimensionality bi-thiophene model are presented, showing highly encouraging results for TDMVCC. It is demonstrated how TDMVCC in the limit of a complete cluster expansion becomes equivalent to multiconfiguration time-dependent Hartree for the same active-space size. Similarly, it is discussed how TDMVCC generally gives better and more stable results than its time-independent-modals counterpart, while equivalent results are obtained for complete expansions and full one-mode basis sets.

Extended vibrational coupled cluster: Stationary states and dynamics.

For the first time, equations are derived for computing stationary vibrational states with extended vibrational coupled cluster (EVCC) and for propagating nuclear wave packets using time-dependent EVCC (TDEVCC). Expressions for energies, properties, and auto-correlation functions are given. For TDEVCC, convergence toward the ground state for imaginary-time propagation is shown, as well as separability in the case of non-interacting subsystems. The analysis focuses substantially on the difference between bra and ket parameterizations for EVCC and TDEVCC compared to normal vibrational coupled cluster (VCC) and time-dependent VCC (TDVCC). A pilot implementation is presented within a new full-space framework that offers easy access to completely general, albeit not efficient, implementations of alternative VCC variants, such as EVCC. The new methods were tested on 35 three- and six-mode molecular systems. Both EVCC[k] and TDEVCC[k] showed good, hierarchical convergence toward the exact limit. This convergence was generally better than for normal VCC[k] and TDVCC[k] and better still than for (time-dependent) vibrational configuration interaction, though this should be balanced with the higher computational complexity of EVCC. The results highlight the importance of exponential parameterizations and separability in general, as seen, in particular, for the TDEVCC bra parameterization, which is in contrast to the partially linear one of TDVCC. With the results being rooted in the general structures of coupled cluster (CC) theory, they are expected to be relevant to other applications of both normal and extended CC theory as well.

MR-MCTDH[n]: Flexible Configuration Spaces and Nonadiabatic Dynamics within the MCTDH[n] Framework.

Solving the time-dependent Schrödinger equation (TDSE) for large molecular systems is a complicated task due to the inherent exponential scaling of the problem. One of the most successful and versatile methods for obtaining numerically converged solutions for small to medium-sized systems is multiconfiguration time-dependent Hartree (MCTDH). In a recent publication [J. Chem. Phys. 2020, 152, 084101] we introduced a hierarchy of approximations to the MCTDH method which mitigate the exponential scaling by truncating the configuration space based on a maximum excitation level w.r.t. a selected reference configuration. The MCTDH[n] methods are able to treat large systems, but the single-reference Ansatz is not optimal in cases where one (or a few) degrees of freedom are special. Examples could be double-well systems, intramolecular vibrational-energy redistribution (IVR) calculations, or nonadiabatic dynamics. In this work we introduce a multireference (MR) extension to the MCTDH[n] methods where selected higher-order excitations for the special degrees of freedom can be introduced in a simple but flexible way. The resulting MR-MCTDH[n] methods allow for, for example, treating nonadiabatic dynamics within the single-set formalism with the wave packets on each electronic surface described using the same level of approximation. Example calculations are performed on formyl fluoride (IVR), salicylaldimine (double well), and pyrazine (nonadiabatic dynamics). The results show that fast convergence is achieved by extending the configuration space in the special modes that govern the quantum dynamics.

Vibrational Coupled Cluster Computations in Polyspherical Coordinates with the Exact Analytical Kinetic Energy Operator.

We present the first use of curvilinear vibrational coordinates, specifically polyspherical coordinates, in combination with vibrational coupled cluster theory. The polyspherical coordinates are used in the context of both the adaptive density-guided approach to potential energy surface construction and in the subsequent vibrational coupled cluster calculations of anharmonic vibrational states. Results obtained based on the polyspherical coordinate parametrization are compared to results obtained with the use of rectilinear vibrational coordinates, namely, normal coordinates and hybrid optimized and localized coordinates for the formaldehyde molecule. This comparison is carried out with the full vibrational configuration interaction model, using the respective fully coupled potential energy surfaces and untruncated kinetic energy operators. The polyspherical coordinates are shown to facilitate an acceleration of convergence for truncated methods when compared to the use of normal coordinates. We furthermore report on calculations on the hydrogen peroxide molecule in the polyspherical coordinate parametrization. The polyspherical vibrational coordinates are shown to perform very well, even for truncated methods, especially when considering the difficulty that rectilinear vibrational coordinates can exhibit in treating complicated internal molecular motion.

Systematic and variational truncation of the configuration space in the multiconfiguration time-dependent Hartree method: The MCTDH[n] hierarchy.

The multiconfiguration time-dependent Hartree (MCTDH) method is a powerful method for solving the time-dependent Schrödinger equation in quantum molecular dynamics. It is, however, hampered by the so-called curse of dimensionality which results in exponential scaling with respect to the number of degrees of freedom in the system and, thus, limits its applicability to small- and medium-sized molecules. To avoid this scaling, we derive equations of motion for a series of truncated MCTDH methods using a many-mode second-quantization formulation where the configuration space is restricted based on mode-combination levels as also done in the vibrational configuration interaction and vibrational coupled cluster methods for solving the time-independent Schrödinger equation. The full MCTDH wave function is invariant with respect to the choice of constraint (or gauge) operators, but restricting the configuration space removes this invariance. We, thus, analyze the remaining redundancies and derive equations for variationally optimizing the non-redundant matrix elements of the constraint operators. As an alternative, we also present a constraint that keeps the density matrices block diagonal during the propagation and the two choices are compared. Example calculations are performed on formyl fluoride and a series of high-dimensional Henon-Heiles potentials. The results show that the MCTDH[n] methods can be applied to large systems and that an optimal choice of constraint operators is key to obtaining the correct physical behavior of the wave function.