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.

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.

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.

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.

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.

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.

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.

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.

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.

Toward Accurate Theoretical Vibrational Spectra: A Case Study for Maleimide.

We employ and combine a number of recent developments in vibrational structure methods to push their current size limitations toward molecules with tens of modes and showcase their availability for the maleimide molecule. In particular, we assess the use of different rectilinear vibrational coordinates, namely, normal coordinates, hybrid optimized and localized coordinates, and flexible adaptation of local coordinates of nuclei coordinates. These different coordinate parameterizations are employed in conjunction with the adaptive density-guided approach to generate potential energy surfaces (PESs). A screening procedure is furthermore introduced, which provides estimates of the importance of individual terms in the PES, resulting in significant reductions in the computational cost of the PES construction. We find that all three sets of coordinates provide approximately the same level of accuracy in vibrational structure calculations and report fundamental excitation energies with a mean absolute deviation of less than 12 cm-1 when compared to experimental data. We expect that similar accuracy in vibrational structure calculations can be achieved for molecules of larger size using the proposed procedures.

A Gaussian process regression adaptive density guided approach for potential energy surface construction.

We present a new iterative scheme for potential energy surface (PES) construction, which relies on both physical information and information obtained through statistical analysis. The adaptive density guided approach (ADGA) is combined with a machine learning technique, namely, the Gaussian process regression (GPR), in order to obtain the iterative GPR-ADGA for PES construction. The ADGA provides an average density of vibrational states as a physically motivated importance-weighting and an algorithm for choosing points for electronic structure computations employing this information. The GPR provides an approximation to the full PES given a set of data points, while the statistical variance associated with the GPR predictions is used to select the most important among the points suggested by the ADGA. The combination of these two methods, resulting in the GPR-ADGA, can thereby iteratively determine the PES. Our implementation, additionally, allows for incorporating derivative information in the GPR. The iterative process commences from an initial Hessian and does not require any presampling of configurations prior to the PES construction. We assess the performance on the basis of a test set of nine small molecules and fundamental frequencies computed at the full vibrational configuration interaction level. The GPR-ADGA, with appropriate settings, is shown to provide fundamental excitation frequencies of an root mean square deviation (RMSD) below 2 cm-1, when compared to those obtained based on a PES constructed with the standard ADGA. This can be achieved with substantial savings of 65%-90% in the number of single point calculations.

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.

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.

Adaptive density-guided approach to double incremental potential energy surface construction.

We present a combination of the recently developed double incremental expansion of potential energy surfaces with the well-established adaptive density-guided approach to grid construction. This unique methodology is based on the use of an incremental expansion for potential energy surfaces, known as n-mode expansion; an incremental many-body representation of the electronic energy; and an efficient vibrational density-guided approach to automated determination of grid dimensions and granularity. The reliability of the method is validated calculating potential energy surfaces and obtaining fundamental excitation energies for three moderate-size chain-like molecular systems. The use of our methodology leads to considerable computational savings for potential energy surface construction compared to standard approaches while maintaining a high level of accuracy in the resulting potential energy surfaces. Additional investigations indicate that our method can be applied to covalently bound and strongly interacting molecular systems, even though these cases are known to be very unfavorable for fragmentation schemes. We therefore conclude that the presented methodology is a robust and flexible approach to potential energy surface construction, which introduces considerable computational savings without compromising the accuracy of vibrational spectra calculations.

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.

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.

Vibrationally resolved emission spectra of luminescent conjugated oligothiophenes from anharmonic calculations.

We report on accurate and efficient calculations of vibrationally resolved emission spectra for oligothiophenes from anharmonic vibrational configuration interaction wave-function calculations in reduced vibrational spaces. These reduced spaces are chosen based on the independent mode displaced harmonic oscillator model. Good agreement with experiment is obtained for all-trans oligothiophenes with two to five rings also when employing only a few active modes. Vibrational modes incorporating inter-ring carbon-carbon stretches and a ring breathing mode are found to be the main players in the vibrational progression for the emission from the first excited electronic state for all investigated oligothiophene derivatives. The presented framework is here illustrated for oligothiophenes, but we have made no underlying system-dependent assumptions and believe it to become a valuable tool for the rational design of fluorescence biomarkers.

Machine learning for potential energy surfaces: An extensive database and assessment of methods.

On the basis of a new extensive database constructed for the purpose, we assess various Machine Learning (ML) algorithms to predict energies in the framework of potential energy surface (PES) construction and discuss black box character, robustness, and efficiency. The database for training ML algorithms in energy predictions based on the molecular structure contains SCF, RI-MP2, RI-MP2-F12, and CCSD(F12*)(T) data for around 10.5 × 106 configurations of 15 small molecules. The electronic energies as function of molecular structure are computed from both static and iteratively refined grids in the context of automized PES construction for anharmonic vibrational computations within the n-mode expansion. We explore the performance of a range of algorithms including Gaussian Process Regression (GPR), Kernel Ridge Regression, Support Vector Regression, and Neural Networks (NNs). We also explore methods related to GPR such as sparse Gaussian Process Regression, Gaussian process Markov Chains, and Sparse Gaussian Process Markov Chains. For NNs, we report some explorations of architecture, activation functions, and numerical settings. Different delta-learning strategies are considered, and the use of delta learning targeting CCSD(F12*)(T) predictions using, for example, RI-MP2 combined with machine learned CCSD(F12*)(T)-RI-MP2 differences is found to be an attractive option.

Approximate high mode coupling potentials using Gaussian process regression and adaptive density guided sampling.

We present a new efficient approach for potential energy surface construction. The algorithm employs the n-mode representation and combines an adaptive density guided approach with Gaussian process regression for constructing approximate higher-order mode potentials. In this scheme, the n-mode potential construction is conventionally done, whereas for higher orders the data collected in the preceding steps are used for training in Gaussian process regression to infer the energy for new single point computations and to construct the potential. We explore different delta-learning schemes which combine electronic structure methods on different levels of theory. Our benchmarks show that for approximate 2-mode potentials the errors can be adjusted to be in the order of 8 cm-1, while for approximate 3-mode and 4-mode potentials the errors fall below 1 cm-1. The observed errors are, therefore, smaller than contributions due to missing higher-order electron excitations or relativistic effects. Most importantly, the approximate potentials are always significantly better than those with neglected higher-order couplings.