---
title: "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"
sdDatePublished: "2026-09-17T04:10:00Z"
source: "https://arxiv.org/pdf/2609.18654"
topics:
  - name: "science and technology"
    identifier: "medtop:13000000"
locations:
  - "Heidelberg"
  - "Île-de-France"
  - "France"
---


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.