q-AQUA-pol Potential Energy Surface

• The model integrates CCSD(T) many-body corrections with the classical TTM-3F framework to achieve high accuracy in water cluster energetics and spectra.
• It employs permutationally invariant polynomials for 2-, 3-, and 4-body terms, yielding low RMS errors across extensive reference datasets.
• Its n-mode representation and heat-bath integral strategy enable efficient grid-based quadrature for quantum vibrational calculations and scalable simulations.

Searching arXiv for the specified papers and closely related q-AQUA/q-AQUA-pol work. arXiv search query: "q-AQUA-pol water potential energy surface" The q-AQUA-pol potential energy surface (PES) is a many-body water PES built from CCSD(T) reference data up to the 4-body level and used, in recent quantum vibrational work, in conjunction with an n-mode representation to enable grid-based quadrature, compressed integral storage, and selected-configuration-interaction treatments of water-cluster spectra (Yu et al., 2022, Tran et al., 16 Sep 2025). In the formulation used for cluster spectroscopy, q-AQUA-pol is expressed as a $\Delta$-correction to the classical polarizable TTM-3F model, while the underlying many-body representation resolves the total energy into 1-, 2-, 3-, and 4-body contributions, with higher-order terms treated as negligible at the CCSD(T) level (Yu et al., 2022).

1. Formal many-body structure

For an $N$-water system, the total potential energy in the q-AQUA framework is written as

$$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$

In the q-AQUA-pol form used for the exact mid-IR vibrational calculations, the PES is built as a $\Delta$-correction to TTM-3F: $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$ where $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$ are all atomic coordinates (Tran et al., 16 Sep 2025).

The TTM-3F reference contains intramolecular one-body monomer potentials, permanent electrostatics with point charges placed on each H and O, induced dipole interactions with Thole damping via site polarizabilities $\alpha_O$ and $\alpha_H$, and damped dispersion and exchange-repulsion at long range. The $\Delta n\mathrm B$ terms correct the TTM-3F $n$-body energies and are constructed so that each correction goes to zero when any subset of molecules is driven apart (Tran et al., 16 Sep 2025).

This decomposition is significant because it combines a physically motivated polarizable baseline with explicit many-body corrections fit to high-level electronic-structure data. A plausible implication is that the model inherits long-range physical structure from TTM-3F while using explicit $N$0 terms to recover the CCSD(T) many-body energetics required for cluster spectroscopy and condensed-phase simulation.

2. Parameterization and analytic components

The one-body monomer term is taken from the Partridge–Schwenke spectroscopic monomer PES. In the q-AQUA construction it is described as analytic, full-dimensional, and permutationally invariant with respect to the two H atoms (Yu et al., 2022).

The higher-body terms are represented by permutationally invariant polynomials (PIPs) and switching constructions that separate short- and long-range physics. The essential components are as follows.

Term	Representation	Reference data / fit quality
$N$1	Partridge–Schwenke monomer PES	Spectroscopic monomer PES
$N$2	7th-order PIP of 42 symmetry functions + $N$3 with switching	71 892 dimer geometries; typical RMS error 25 cm$N$4
$N$5	4th-order short-range PIP + 3rd-order longer-range PIP with switching	45 332 trimer configurations; RMS $N$6 cm$N$7 and $N$8 cm$N$9
$$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$0	Improved CCSD(T) PIP with 200 grouped PIPs plus selected 4th-order terms	3 692 tetramer configurations; RMS $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$1 cm$$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$2

For the 2-body term, the short-range part is fit by a 7th-order PIP in the six inter- and intra-monomer distances, while the long-range part is described by a high-level dipole-dipole interaction $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$3 built on an accurate CCSD(T) dipole-moment surface. A smooth switching function $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$4 blends the two in a window $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$5: $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$6

For the 3-body term, two complementary PIP fits are used: a short-range fit for $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$7 Å in Morse-type variables of OO, OH, and HH distances with 222111 symmetry, and a longer-range fit for $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$8 Å in inverse distances with the same symmetry; the two are joined by a second switching function $$V(1,\dots,N) =\sum_{i=1}^N V_{1\mathrm b}(i) +\sum_{i<j} V_{2\mathrm b}(i,j) +\sum_{i<j<k} V_{3\mathrm b}(i,j,k) +\sum_{i<j<k<l} V_{4\mathrm b}(i,j,k,l)\,.$$9 that ensures analytic continuity. The 4-body term is based on an improved version of a CCSD(T) PIP, fitted with a compact basis of 200 grouped PIPs plus a selected set of 4th-order terms, and multiplied by a final switch $\Delta$0 at large tetramer size (Yu et al., 2022).

The fitting protocol uses PIP bases that enforce permutation invariance of monomers and H atoms. “Purification” removes basis functions that are exactly redundant under the molecular symmetry group, and “pruning” discards high-order monomials that contribute negligibly to the fit or lead to ill-conditioning. Coefficients are obtained by linear least squares, analytic gradients are generated efficiently via reverse-mode differentiation, and no empirical regularization is used; fit quality is judged by RMS errors on large test sets (Yu et al., 2022).

3. n-mode representation and quadrature realization

In the vibrational calculations, the q-AQUA-pol PES is not used directly on a full-dimensional grid. Instead, it is expanded in a normal-mode coordinate basis $\Delta$1 by an n-mode representation truncated at $\Delta$2: $\Delta$3

More formally, if $\Delta$4 with $\Delta$5, then

$\Delta$6

with each higher-order contribution defined by subtraction of all lower-order terms. For example,

$\Delta$7

This truncation avoids the full-dimensional grid while still capturing anharmonic coupling up to three modes simultaneously. The matrix elements required for the vibrational Hamiltonian,

$\Delta$8

factorize into products of one-, two-, or three-dimensional integrals because of the n-mode form (Tran et al., 16 Sep 2025).

One-dimensional integrals are evaluated by Gauss–Hermite quadrature,

$\Delta$9

and two-mode and three-mode terms are treated by corresponding product quadratures. Within selected CI, the two-mode integrals

$$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$0

are precomputed once, sorted in descending magnitude, and stored in compressed form. This “heat-bath” ordering allows the HCI selection step to test excitations in order of largest coupling rather than scanning all possible combinations; one-mode integrals are stored analogously (Tran et al., 16 Sep 2025).

4. Benchmark performance from clusters to liquid water

The original q-AQUA validation emphasizes both cluster benchmarks and condensed-phase observables. For water hexamer isomers, the dissociation-energy mean absolute error versus CCSD(T)/CBS benchmarks is under $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$1 kcal/mol, and harmonic-frequency MAEs for Prism, Cage, Book, and Ring are $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$2–$$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$3 cm$$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$4, with q-AQUA generally matching or slightly improving over MB-pol (Yu et al., 2022).

Unrestricted diffusion Monte Carlo (DMC) zero-point energies relative to the Prism equilibrium were reported for Prism, Cage, and Book1 as $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$5 cm$$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$6, $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$7 cm$$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$8, and $$V_{\mathrm{qAQUA}}(\mathbf R) = V_{\mathrm{TTM\text{--}3F}}(\mathbf R) +\sum_{i<j}\Delta V_{2\mathrm B}(\mathbf r_i,\mathbf r_j) +\sum_{i<j<k}\Delta V_{3\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k) +\sum_{i<j<k<\ell}\Delta V_{4\mathrm B}(\mathbf r_i,\mathbf r_j,\mathbf r_k,\mathbf r_\ell)\,,$$9 cm$\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$0, respectively. Removing the 4-body term shifts these by up to $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$1 cm$\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$2, which was used to confirm the physical importance of the 4-body interaction. For $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$3 isomers A3, A2d, and 9, q-AQUA gives binding energies of $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$4, $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$5, and $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$6 kcal/mol, within $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$7 kcal of the stated MP2/CBS benchmarks, whereas MB-pol underbinds by $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$8 kcal. Preliminary DMC on the 20-mer shows no “holes” and stable propagation (Yu et al., 2022).

Condensed-phase validation was carried out with classical MD, PIMD, and RPMD in i-PI for 256 waters under periodic boundary conditions. At 298 K, a model truncated at 1- and 2-body gives severe over-structuring and peak misplacement in the radial distribution function $\mathbf R=(\mathbf r_1,\dots,\mathbf r_N)$9, inclusion up to 3-body gives correct peak positions but amplitudes that are too large, and the full 4-body model with PIMD yields nearly quantitative agreement with the x-ray data of Skinner et al. For self-diffusion, the reported values at 298 K are $\alpha_O$0 Å$\alpha_O$1/ps for classical MD, $\alpha_O$2 Å$\alpha_O$3/ps for RPMD, and $\alpha_O$4 Å$\alpha_O$5/ps for experiment; nuclear quantum effects increase $\alpha_O$6 by roughly $\alpha_O$7 at 298 K. The orientational relaxation time $\alpha_O$8 at 298 K is $\alpha_O$9 ps classically and $\alpha_H$0 ps in RPMD, compared with an experimental $\alpha_H$1 ps. Turning the 4-body term off raises $\alpha_H$2 by $\alpha_H$3–$\alpha_H$4 (Yu et al., 2022).

5. Role in exact mid-IR quantum vibrational spectroscopy

The q-AQUA-pol PES is the central interaction model in selected-configuration-interaction calculations of the zero-temperature mid-infrared ($\alpha_H$5–$\alpha_H$6 cm$\alpha_H$7) vibrational spectra of the water monomer, dimer, trimer, and hexamer in its cage and prism geometries (Tran et al., 16 Sep 2025). In that setting, the PES is combined with the n-mode representation specifically to facilitate grid-based quadrature and integral storage.

Within selected configuration interaction, a new spectral strategy is introduced that is complementary to eigenstate enumeration. Instead of enumerating excited states directly, the spectrum is calculated with the response-vector method, and the resulting system of linear equations is solved using a basis of configurations that is optimally selected at each frequency of interest. The pre-sorted one- and two-mode integrals inherited from the q-AQUA-pol/n-mode representation are reused in this frequency-resolved HCI procedure (Tran et al., 16 Sep 2025).

The reported computational savings are explicit: the spectral HCI approach yields sizable savings in both storage, by a factor of $\alpha_H$8 in the trimer, and CPU time, with minutes per frequency point on a single node for $\alpha_H$9. The same study compares the resulting spectra to previous work and highlights limitations of the local monomer approximation. It also states that, to the best of its authors’ knowledge, the hexamer spectra are the most accurate ones reported to date (Tran et al., 16 Sep 2025).

6. Accuracy bounds, limitations, and computational significance

Several distinct accuracy statements are attached to q-AQUA-pol as used in vibrational applications. The underlying four-body-corrected fit reproduces benchmark CCSD(T) harmonic frequencies of the water hexamer to within $\Delta n\mathrm B$0 cm$\Delta n\mathrm B$1 on average, as cited to Howard and Tschumper. In the present vibrational calculations, after freezing low-frequency modes and truncating at three-mode coupling, absolute fundamental frequencies are converged to better than $\Delta n\mathrm B$2–$\Delta n\mathrm B$3 cm$\Delta n\mathrm B$4 compared to previous VCI on the HBB PES, with systematic behavior (Tran et al., 16 Sep 2025).

The main limitations stated for the vibrational implementation are associated with the n-mode truncation. The three-mode truncation can introduce spurious wells or unbound behavior in very floppy, low-frequency modes, especially torsions. To avoid this, all modes below $\Delta n\mathrm B$5 cm$\Delta n\mathrm B$6 were frozen in the reported work. As with all many-body expansions, $\Delta n\mathrm B$7 and higher effects are neglected, although for water they are described as very small beyond four bodies (Tran et al., 16 Sep 2025).

The computational advantages are equally specific. The n-mode expansion breaks a $\Delta n\mathrm B$8-dimensional integral into sums of at most three-dimensional integrals, making full grid quadrature tractable up to $\Delta n\mathrm B$9 modes. Using Gauss–Hermite quadrature and storing only one-, two-, and three-mode integrals minimizes memory and flops, while heat-bath ordering makes variational-space construction efficient. In the original q-AQUA implementation for periodic simulations, the total energy and gradient evaluation for 256 waters up to 4-body with periodic boundary conditions is $n$0 s per step on a single 2.4 GHz Intel Xeon core, and with 8-core OpenMP this drops by a factor of $n$1; the 4-body work is the most expensive, at $n$2 tetramers evaluated per step, and truncating at the 3-body level roughly halves the cost (Yu et al., 2022).

Taken together, the q-AQUA-pol many-body fit, the n-mode/Gauss–Hermite representation, and the heat-bath integral strategy make fully quantum, near-exact mid-IR spectra of water clusters up to the hexamer feasible and reliable at the $n$3–$n$4 cm$n$5 level, while the underlying q-AQUA model also reproduces cluster energetics, liquid structure, and diffusion with accuracy that remains anchored to CCSD(T) reference data (Tran et al., 16 Sep 2025, Yu et al., 2022).