Ali Al-Jaaidi, David Lauvergnat, Daniel Peláez; variational computation of anharmonic vibrational eigenstates using bound QFF in Orsay, France; CP-bQFF maintains accuracy with reduced cost

Ali Al-Jaaidi, David Lauvergnat, Daniel Peláez; variational computation of anharmonic vibrational eigenstates using bound QFF in Orsay, France; CP-bQFF maintains accuracy with reduced cost

Variational computation of anharmonic ground and excited vibrational eigenstates using bound Quartic Force Fields: Application with MCTDH and ElVibRot

Variational computation of anharmonic ground and excited vibrational eigenstates using bound Quartic Force Fields: Application with MCTDH and ElVibRot Ali Al-Jaaidi,† David Lauvergnat,∗,‡ and Daniel Pel´aez∗,† †Universit´e Paris-Saclay, CNRS, Institut des Sciences Mol´eculaires d’Orsay, 91405, Orsay, France. ‡Universit´e Paris-Saclay, CNRS, Institut de Chimie Physique, 91405, Orsay, France. E-mail: David.Lauvergnat@universite-paris-saclay.fr; Daniel.Pelaez-Ruiz@universite-paris-saclay.fr 1 arXiv:2609.18654v1 [physics.chem-ph] 16 Sep 2026

Abstract In this work we introduce the use of Quartic Force fields (QFF) potential expansions in the context of variational calculations. Such potentials are commonly employed in molecular Vibrational Second-Order Perturbation Theory (VPT2) studies, for which equations explicitly dependent on the QFF parameters exist. However, QFF are un- bound potentials for more or less large displacements from the reference point and, most of the time, this prevents their use in conjunction with variational wavepacket-based calculations. In this work, we propose a general correction to QFFs and introduce a fully automated numerical approach to avoid their unbound character. Our corrected potentials, bound QFF (bQFF), do not exhibit appreciable modification of the local topography around the region of interest for infrared spectroscopy. As a consequence of this, we can affirm that the vibrational eigenstructure (eigenvalues, eigenstates) re- mains essentially unaltered by our correction. To illustrate their numerical stability, we have interfaced our bQFF routines in combination with to two well-established quan- tum simulation software packages MCTDH and ElVibRot which feature variational approaches. More specifically, our bQFFs are separable and hence directly expressible as MCTDH operators. Furthermore, concerning the size of our bQFF expansion, we show that it is possible to tensor-decompose our bQFF in Canonical Polyadic form (CP-bQFF). We use the Monte Carlo Canonical Polyadic decomposition algorithm for this. CP-bQFF results are virtually identical to uncompressed bQFF, but the compu- tational efficiency is largely improved. Our approach paves the way for the automated variational study of anharmonic eigenstates in molecular systems within the reach of QFF-based potentials, using either time-dependent or time-independent schemes. Introduction Except for notable exceptions,1 the identification and assignment of chemical species from experimental vibrational spectra rely on the ability to accurately predict the molecular vibra- tional structure. This is more so in cases where recreating of the actual physical conditions is 2

very hard, such as interstellar medium,2 extreme environments,,3,4 as well as in the context of unstable or charged species, which require special experimental techniques which, in turn, may perturb the targeted signal.5–7 In this work by accurate prediction (or simulation) we solely refer to methodologies (and associated software packages) which enable a variational quantum description of nuclei, such as the Multiconfiguration Time-Dependent Hartree Method (MCTDH)8 or ElVibRot.9 In this respect, we also exclude the use of scaling factors (frequencies, intensities) except in the case of normalization or translation of the full spectrum. As the most prominent examples of these variational approaches, due to their relevance, one has the vibrational simulation of small water clusters,10–12 as well as their related charged counterparts6,13–15 Fully general vibrational simulations are very hard tasks typically involving several re- search groups and the use of highly specialized methodologies.16 Indeed, one needs to obtain a set of coordinates adapted to the energy regime and geometry of the system under consid- eration, the corresponding Kinetic Energy Operator (KEO)17,18 together with an accurate representation of the Potential Energy Surfaces (PES) and Dipole Operators (assuming a dipole mechanism). The latter quantities, in turn, involve accurate electronic structure cal- culations19 and cumbersome multidimensional fits, probably in the form of Deep Neural Networks.20–25 Arguably, the availability of a suitable PES is a major limiting factor. In this respect, it should be noted that some approaches might require a specific form (e.g., sum-of- products, SOP, as in some MCTDH implementations8,26). In such cases, the PES should be already of separable form,27,28 otherwise a further refit (i.e. tensor-decomposition) into a sep- arable sum-of-products is necessary.25,29–33 In contrast, some other methods are compatible with any type of potential expression (e.g. MCTDH using CDVR,34 ElVibRot,9TROVE35 or GENIUSH.36 Fortunately, such in-depth studies can be greatly simplified in the rather common case of semi-rigid molecules characterized by a single well. Indeed, for these, a local representation of the PES around a given minimum (assumed stable enough) suffices to determine its vibrational spectrum. Accurate local representations can be built using: 3

(i) many-body type expansion;37 (ii) adaptive methods guided by convergence of a molec- ular property (such as the Zero-Point Energy, ZPE);38,39 or, an even simpler strategy, (iii) Taylor-type expansions of the energy up to a given order (fourth, sixth) around a reference point. The latter, usually referred to as Force Fields (FF) in spite of the obvious ambigu- ity of the term, have been cleverly exploited in the context of, probably, the most popular family of fully quantum vibrational simulations, the Vibrational Second-order Perturbation Theory (VPT2).40–42 VPT2 is method of choice for the vast majority of users due to its relatively black box nature as epitomized by its successful application in the context astro- chemistry, in spite of its intrinsic inability to deal with resonances such as Coriolis, Fermi and Darling-Dennison.43 Note that it is possible to tackle such issues beyond a mere per- turbative approach.44 For the sake of completeness, it should be noted that to overcome the aforementioned limitations, a hierarchy of approaches, in the same spirit of quantum chem- istry VSCF, VCI, VCC, VMCSCF, VCASSCF,… exists.45 Finally, it should be highlighted that the interest in molecular vibrational eigenstates goes beyond the mere computation of fundamentals but also determination of IR spectra as well as initial conditions for Molecular Dynamics.46 In this work, we will focus on semi-rigid molecules whose PES is given by a so-called Quartic Force Fields (QFFs).41,47 Note that sextic force fields in the case of small molecules have been successfully employed.48–50 Furthermore, we will assume that an appropriate choice of the level of level of electronic structure theory has been chosen19,51 and that the necessary numerical precautions have been taken into account.52 Quartic Force Fields (QFF) are Taylor expansions of a PES up to fourth order around a given reference position, Q0, usually an equilibrium geometry, and is given by: V QFF(Q) = Vref + 1 2! f X i ΦiiQ2 i + 1 3! f X i≥j≥k ΦijkQiQjQk + 1 4! f X i≥j≥k≥l ΦijklQiQjQkQl (1) Here f represents the number of internal degrees of freedom (3N-6 or 3N-5), Qκ is 4

a (dimensionless, vide infra) displacement along normal mode κ, and Φii = ∂2V ∂Q2 i , Φijk = Pijk ∂3V ∂Qi∂Qj∂Qk , Φijkl = Pijkl ∂4V ∂Qi∂Qj∂Qk∂Ql are partial derivatives of the full potential, V (Q), along the indicated modes (i, j, k, . . .). Vref is the energy at the reference position, Vref = V (Q0) and P{µ} is the factor accounting for the total number of derivatives related by per- mutation of the {µ} set of indices. Note in passing that it is common to restrict to the this expression to the so-called semi-QFF, a truncation of the quartic terms considering only at most two-mode couplings.53 As evident from their definition, QFFs are unbound potentials, and this prevents their universal application in combination with variational methods, which are at the core of our interests. Indeed, while harmonic frequencies are strictly positive for minima, the cubic terms and/or the quartic ones for Φijkl < 0 and large displacements may become (largely) negative, i.e. unbound, hence exhibiting holes. To illustrate this, in Figure 1 we compare a Morse potential with its Taylor expansions up to 10th order. Figure 1: Comparison of Taylor expansions of increasing order of a Morse potential (stretch- ing mode) around its equilibrium geometry. As a consequence, any variational quantum simulation, using either time-dependent or time-independent methods, is doomed to failure since regions of configuration space asso- ciated to such unphysically distorted unbound QFF will be eventually populated. This is 5

epitomized in the context of vibrational spectroscopy.6,14,30 To illustrate this, we present below a trivial example: the relaxation of the ground state (GS) of the water molecule. Ob- viously, the expected result is a constant value, the Zero-Point Energy of water. In Figure 2, we compare the behavior of the relaxation energy with time for two different QFF potentials, namely: (i) an unbound one (red, lower curve at long times) and (ii) a bound one (green, upper curve). Note that our graph uses a double logarithmic scale. In both simulations, the same initial wave function, built as a product of the 1D anharmonic ground states was employed. In the upper curve (green), the bound PES shows the expected behavior. In contrast, for the case of the unbound potential, a sudden drop in the energy towards neg- ative (non-physical) values below the reference energy, Vref) is observed, thus signaling the presence of a hole. Both potentials are reported in the Supplementary Information (SI). Figure 2: Wavefunction relaxation54,55 on a bound and an unbound QFF water Potential. For a QFF (or any high-order polynomial expression), detecting holes by, for instance, finding polynomial roots as a signature of a hole, is extremely cumbersome. To avoid such explorations, it is common to resort to simple solutions such as applying potential cuts in terms of energy or coordinate value.56 Recently, Poirier and coworkers57,58 have introduced more systematic approaches to determine the actual locations of holes in full configuration space or in a subset of it (up to 30 DOFs). 6

In this work, we propose an automated system-independent procedure for the smooth (continuous, differentiable) correction of QFFs in semi-rigid systems, irrespective of their size. We denote the resulting PES as bound QFF (bQFF). To the best of our knowledge, no general correction method for a PES presenting holes have been presenting before. In a nutshell, our algorithm relies on analyzing the couplings among the various modes at different orders and applying separable corrections that preserve the description of the original vibrational levels up to combination bands while rendering the PES bound at large displacements. This effectively removes holes thus enabling its interface with variational methods, in general, as well as a symbolic expression as in the case of an MCTDH operator. The remaining part of the article is structured as follows: in the Methodology Section, we describe all algorithmic and numerical details of our approach, more specifically, the gen- eration of a bQFF, its efficient compression in Canonical Polyadic form, and all aspects con- cerning electronic structure calculations for our case studies. Section Results and Discussion presents two applications, a benchmark system which serves to illustrate straightforwardly our methodology and a full-dimensional study in 9D. Finally, in the Appendix we provide a detailed description of our quantum dynamical approaches and in the Supplementary Infor- mation we present extra data and graphs supporting our results.