On the separability of large-amplitude motions in anharmonic frequency calculations

Abstract

Nuclear vibrational theories based upon the Watson Hamiltonian are ubiquitous in quantum chemistry, but are generally unable to model systems in which the wavefunction can delocalise over multiple energy minima, i.e. molecules that have low-energy torsion and inversion barriers. In a 2019 Chemical Reviews article, Puzzarini et al. note that a common workaround is to simply decouple these problematic modes from all other vibrations in the system during anharmonic frequency calculations. They also point out that this approximation can be “ill-suited”, but do not quantify the errors introduced. In this work, we present the first systematic investigation into how separating out or constraining torsion and inversion vibrations within potential energy surface (PES) expansions affects the accuracy of computed fundamental wavenumbers for the remaining vibrational modes, using a test set of 19 tetratomic molecules for which high quality analytic potential energy surfaces and fully-coupled anharmonic reference fundamental frequencies are available. We find that the most effective and efficient strategy is to remove the mode in question from the PES expansion entirely. This introduces errors of up to +10 cm⁻¹ in stretching fundamentals that would otherwise couple to the dropped mode, and ±5 cm⁻¹ in all other fundamentals. These errors are approximately commensurate with, but not necessarily additional to, errors due to the choice of electronic structure model used in constructing spectroscopically accurate PES.

---

1 Introduction

Harmonic normal mode analysis and scaled harmonic frequencies continue to enjoy great popularity in assisting the assignment of vibrational infrared and Raman spectra,¹ but can fail to accurately model anharmonically coupled systems, e.g. the OH stretching spectrum of carboxylic acid dimers.² Moving beyond the harmonic approximation, the complexity and computational cost associated with constructing anharmonic potential energy surfaces (PES) and solving the nuclear vibrational Schrödinger equation usually demand compromise, and the accuracy of anharmonic models is sensitive both to model ansatz and implementation details.

Arguably one of the most important ansatz choices is that of the coordinate system, because this impacts upon both the efficiency of PES construction schemes and the complexity of the nuclear vibrational Hamiltonian. As stressed by Sibert III et al., “exact implementations of the curvilinear and rectilinear descriptions are strictly equivalent. However, a key practical point is that they are not equivalent in low orders of approximation.”³ Usually, this is a choice between Scylla and Charybdis, because the kinetic energy operator is more concisely formulated in rectilinear coordinates, whereas molecular potential energy surface expansions converge faster in curvilinear coordinates.

Anharmonic Watson Hamiltonian models^4,5 are widely used but share common limitations when it comes to describing large-amplitude vibrational motions. These limitations arise from the aperiodic nature and shape of rectilinear normal mode coordinates. Watson Hamiltonian models are strictly only suitable for semi-rigid molecules that undergo relatively low-amplitude vibrations from their equilibrium structures, and can yield strongly divergent results when applied to molecules that undergo large-amplitude motions such as torsion or inversion vibrations.

The simplest and most rigorous approximation to circumvent this involves enforcing the required condition of semi-rigidity by ignoring any coupling between the large-amplitude vibrations and all other vibrations. In a recent review, Puzzarini et al. pointed out that this “represents a crude approximation, which can be ill-suited in many situations.”⁶ However, to the best of our knowledge, the resultant errors have not been systematically investigated or quantified before.

If the large-amplitude vibration itself is of interest, one needs to account for its curvilinear nature using more sophisticated models.⁷ In some special cases, the Watson Hamiltonian is still suitable if the potential energy surface (PES) and wavefunction can be expanded about a saddle point connecting two symmetric minima. This approach has been used to compute tunnelling splittings for the inversion modes of H₃O⁺ and NH₃,^8,9 but cannot be extended to other types of motions, such as torsions.⁹

A more general solution is to numerically map out the PES along one or two large-amplitude coordinates, constrained to remain orthogonal to the remainder of the PES. This idea forms the basis for a series of closely-related methods, including the Hougen–Bunker–Jones Hamiltonian,¹⁰ reaction path Hamiltonian,¹¹ reaction surface Hamiltonian¹² and internal coordinate path Hamiltonian^13,14 models. However, the process of mapping out the PES along these special coordinates is far from routine, and a specialised set of basis functions are required in which to expand the wavefunction. This complicates the evaluation of kinetic and potential energy integrals, decreasing computational efficiency relative to methods based purely on the Watson Hamiltonian.

To investigate the limits of applicability of “pure” Watson Hamiltonian models, we will quantify how omitting or approximating PES-mediated coupling between torsion and inversion vibrations affects fundamental transition wavenumbers of other vibrational modes. In principle, the best way to assess the accuracy of any computational model is to directly benchmark against experimental data. However, in the context of vibrational spectroscopy, problems with this approach arise from both theoretical and experimental perspectives.

From the experimental side, it is not always possible to obtain high-resolution gas phase vibrational spectra and, even if the data is available, unique and unambiguous assignments are not always possible, particularly where there are strong resonances between fundamentals and overtones or combination bands. From the theoretical side, the main problem is that there are multiple sources of error that all contribute to the overall accuracy of the final predictions, including choice of electronic structure model used in constructing the potential energy surface, form and/or completeness of potential energy surface representation, and choice of nuclear vibrational model or form of nuclear vibrational Hamiltonian. It is important to carefully disentangle contributions from each potential source of error to ensure results are robust; i.e. when predictions align with experiment, they are right for the right reasons, not simply due to fortuituous error cancellation.¹⁵

Therefore, we will perform a thorough and systematic benchmarking study on a series of 19 tetratomic molecules for which high quality semi-global analytical potential energy surfaces and reference fundamental wavenumbers obtained using alternative nuclear vibrational models are available.¹⁶

However, the acid test of any computational model is ultimately its ability to predict and/or interpret experimental data. Informed by the tetratomic benchmarking results, we will construct Watson Hamiltonian models in which the computational sources of error are carefully controlled, and assess their accuracy against high-resolution experimental data for a series of pentatomic and hexatomic molecules, including the textbook example methanol.

2 Theory

The first choice that must be made in any nuclear vibrational calculation is how to represent the potential energy surface. For harmonic frequency calculations this is straightforward; all required information is encoded within the molecular Hessian. However, for anharmonic frequency calculations, there are a number of options for constructing semi-global potential energy surfaces, including force field parameterisation,^17–20 numerical mapping on multi-modal grids,^8,21 and symmetry-invariant fitting or interpolation.^22,23

In this work, we choose to use the PyPES and PyVCI program packages^16,24 that utilise the Watson Hamiltonian⁴ in conjunction with sextic force field expansions of the potential energy. Although analytic force field expansions have gone out of fashion in variational calculations, in favour of m-mode grid-based representations of the PES,^8,21 we choose to use PyPES because it conveniently offers all necessary tools: PyPES easily facilitates the utilisation of spectroscopically accurate PES compiled from the literature and newly ab initio-generated force fields, allowing theory-theory and theory-experiment benchmarking to be carried out within the same methodological framework. The required PES and wavefunction manipulations to separate out large-amplitude vibrations are trivial to implement in the context of force field-expanded PES and wavefunctions constructed from a basis of harmonic oscillator product functions. Sextic force field expansions in rectilinear normal mode (RNM) coordinates are obtained via coordinate transformation from quartic force fields in curvilinear normal mode (CNM) coordinates, following the formalism established by Allen and Császár.¹⁹ This takes advantage of the rapid convergence of force field expansions in CNM coordinates, and provides an efficient and numerically stable procedure for approximating 5th and 6th order RNM force constants.¹⁶ Constraints that restrict how force fields describe large-amplitude modes – or indeed, omit them entirely – may be applied in either CNM or RNM space. Because the truncation schemes trialled in this work depend on how the force fields are constructed, we will briefly summarise our approach here.¹⁶

2.1 Force field construction

Curvilinear normal mode (CNM) coordinates, Q̃, are defined as linear combinations of (possibly redundant) internal coordinates, S, that reproduce rectilinear normal mode (RNM) coordinates, Q, to first order, i.e. whose first derivatives with respect to Cartesian displacements are the same:In this work, our choice of primitive internal coordinates for bond lengths, bond angles, dihedral angles, and out-of-plane angles arewhere r, θ and τ are the usual bond length, bond angle and dihedral coordinates,²⁵ sin(ω) is an out-of-plane umbrella coordinate,²⁵ and f_SPF(r) is the Simons–Parr–Finlan radial coordinate²⁶ that is defined relative to equilibrium and has appropriate asymptotic behaviour at both short- and long-range by design.

Rectilinear normal mode coordinates (Q) are defined as linear combinations (L⁻¹) of Cartesian displacements (X) that diagonalise the mass-weighted Hessian matrix:²⁷

Q = L⁻¹X. (3)

Curvilinear normal mode coordinates (Q̃) are defined analogously:

Q̃ = (BL)⁻¹S = RS, (4)

where B is the Wilson B matrix that contains derivatives of internal coordinates with respect to Cartesian displacements. Contracting the B and L matrices and inverting yields the R matrix that defines curvilinear normal mode coordinates as linear combinations of internals.

With curvilinear normal mode coordinates thus defined, the PES can be expanded as a Taylor series in Q̃ about equilibrium (unrestricted summation):

where the force constants, F̃, represent derivatives of the energy with respect to displacements along curvilinear normal modes. In general, these force constants may be stably computed to 4th order by numerical differentiation involving analytic ab initio Hessians. However, for molecules contained within the PyPES library, they may instead be obtained by analytical differentiation of implemented reference potentials.²⁸ Our procedures for computing the required force constants, either analytically or numerically, are summarised in Fig. 1.

Fig. 1 Two equivalent schematic representations of our process for deriving Taylor series expansions of potential energy surfaces, V, in curvilinear (CNM) and rectilinear normal mode (RNM) coordinates; QFF = quartic force field = 4th order Taylor series expansion; SFF = sextic force field = 6th order Taylor series expansion. In the bottom diagram, the bracketed superscripts indicate overall order of expansion, and order of expansion with respect to any given coordinate (mode representation), respectively. The symbol denotes a numerical differentiation process. Although curvilinear normal mode force constants may be computed directly via numerical differentiation, it is easier and more numerically robust to first construct full quartic force fields in rectilinear normal mode coordinates and then analytically transform into curvilinear normal mode coordinates.

Transforming the 4th order curvilinear normal mode force constants F̃ to 6th order rectilinear normal mode force constants F requires derivatives of Q̃ with respect to Q, to 5th order. These are computed via a two-step linear transformation from higher order analogues of the B matrix, using the L and R matrices defined above. Following a non-linear coordinate transformation procedure,^19,29 the desired sextic force field in rectilinear normal mode coordinates is obtained (unrestricted summation):

The 5th and 6th order force constants F_ijklm and F_ijklmn are not the true higher order derivatives of the potential energy, since they are missing contributions from F̃_ijklm and F̃_ijklmn. This approximation can result in errors of up to 10 cm⁻¹ in computed fundamental wavenumbers, but usually they are 1–2 cm⁻¹.¹⁶

2.2 Proposed force field truncation schemes

A natural consequence of the semi-local nature of the Watson Hamiltonian is that resultant nuclear vibrational models cannot generally capture tunnelling splittings. Nonetheless, full dimensional Watson Hamiltonian models can – in theory – provide reasonable approximations for higher-frequency modes, provided that tunnelling splittings are small. However, in practice, force field expansions often become unstable when applied to molecules with torsional or inversion modes of vibration.^16,20,30

Possible sources of instability include:

(1) Unphysical divergences or “holes” in the potential arise along the torsion/inversion coordinate, reflecting the inability of the truncated Taylor series to correctly model the shape of the PES along the vibrationally-accessible region.

(2) The inability of 6th order expansions in RNM coordinates to correctly reproduce the shape of 4th order expansions in CNM coordinates, resulting in unphysical couplings between torsion/inversion coordinates and others within the RNM force field.

(3) Both of the above.

We therefore investigate three possible truncation schemes:

(1) Ṽ_harm: constraining the problematic coordinate to remain harmonic in the CNM PES expansion, and then transforming to RNM coordinates.

(2) V_drop: dropping all terms involving the problematic coordinate in the RNM PES expansion, but retaining indirect coupling mediated through the curvilinear-to-rectilinear coordinate transformation process (see ref. 16).

(3) Ṽ_drop: dropping all terms involving the problematic coordinate during CNM PES construction prior to transformation into rectilinear normal mode coordinates. This ensures that the dropped mode remains frozen (harmonic and geometrically uncoupled) in the RNM PES expansion. It may become necessary to modify the redundant coordinate set S, which underpins the curvilinear normal mode coordinates Q̃ (see eqn (4)).

Transformation procedures used to obtain these truncated PES expansions are summarised in Table 1.

Table 1 Shorthand notation for potential energy surface expansions constructed in this work. The tilde superscript denotes that mode-dropping or harmonic constraints have been applied in curvilinear normal mode coordinates preceding transformation to obtain the final force field expansion in rectilinear normal mode coordinates. Values in superscripted parentheses are used to indicate both the overall expansion order and order with respect to specific coordinates. Computational scaling laws reflect the number of Hessian calculations that would be required to directly compute curvilinear normal mode force constants ab initio

Name	Definition	Scaling
V ref	V ^(4)(Q̃^(4)) → V^(6)(Q^(6))	[scr O, script letter O](Nvib^2)
Ṽ harm	V ^(4)(Q̃^(2)tor/inv,Q̃^(4)bend/str) → V^(6)(Q^(6))	[scr O, script letter O](Nbend/str^2)
V drop	V ^(4)(Q̃^(4)) → V^(6)(Q^(0)tor/inv,Q^(6)bend/str)	[scr O, script letter O](Nvib^2)
Ṽ drop	V ^(4)(Q̃^(0)tor/inv,Q̃^(4)bend/str) → V^(6)(Q^(0)tor/inv,Q^(6)bend/str)	[scr O, script letter O](Nbend/str^2)

---

2.3 Vibrational Hamiltonian and wavefunction

Anharmonic fundamental transition wavenumbers are computed with our freely available PyVCI program package, using vibrational configuration interaction (VCI) theory, as detailed in ref. 24. PyVCI utilises the J = 0 Watson Hamiltonian⁴ for non-linear molecules in rectilinear normal mode coordinates, including first-order Coriolis coupling between rotations and vibrations (in atomic units):where B^α_e is the equilibrium rotational constant about axis α and ζ_ij are zeta matrix elements as defined by Meal and Polo.^31,32

The VCI wavefunction is formed as a linear combination of harmonic oscillator product functions:

where each basis function is given by:in which ϕ_ni are single-mode harmonic oscillator functions, and n is a string of quantum numbers n₁,…,n_i,…,n_M that specify the vibrational state across all M modes. The maximum value that the sum of these quantum numbers can take is referred to as the excitation level, denoted VCI(n), where .

VCI fundamentals are identified as those singly-excited states with the largest overlap to the corresponding harmonic oscillator fundamentals, and fundamental transition wavenumbers are computed from the energy difference between these and the VCI ground state.

3 Methods

3.1 Benchmarking – tetratomics

Our benchmarking test set contains 19 tetratomic molecules (plus 5 deuterated isotopologues) that have either a torsion or inversion mode, for which spectroscopically accurate semi-global potential energy surfaces are available in the literature and implemented within the PyPES program package.²⁸ The electronic structure methods, basis sets and nuclear vibrational methods used to obtain the original PES and reference fundamentals are listed in Table 2, and molecular topology classes illustrated in Fig. 2. Trigonal pyramidal molecules (type A) possess symmetry-equivalent minima potentially accessible through inversion. Extended molecules (types C and D), on the other hand, each possess a torsional mode. Trigonal planar molecules (type B) are included as a control group, since they only have a single minimum energy structure accessible through vibrational motion.

Table 2 Key details of methods used to obtain computational reference data for all tetratomic molecules in our benchmarking data set together with (high-resolution) experimental references; c and t denote cis and trans isomers, respectively

		Expt.^a ref.	Electronic structure method	Nuclear vib. model	Comp. ref.
ICH: internal coordinate Hamiltonian, as implemented in TROVE^77 and RVIB4.^78 CVPTn: nth order Canonical Van Vleck perturbation theory, as implemented in VANVLK.^79,80 VPTn: nth order vibrational perturbation theory, as implemented in SPECTRO.^81 VCI: vibrational configuration interaction, as implemented in MULTIMODE.^82a High- and low-resolution gas phase experimental data.b With effective core potential (ECP) on metal atom, refer to original reference for details.c Near-complete atomic orbital basis, refer to original reference for details.
(A)
NH3	C 3v	33–36	CCSD(T)/CBS	ICH (TROVE)	37
PH3	C 3v	38, 39	CCSD(T)/CBS	ICH (TROVE)	40
SbH3	C 3v	41	CCSD(T)/aVQZ^b	ICH (TROVE)	42
BiH3	C 3v	43	CCSD(T)/aVQZ^b	ICH (TROVE)	42
SiH3^−	C 3v		CCSD(T)-R12/BS1^c	CVPT6	44
(B)
BF3	D 3h	45–47	CCSD(T)/VTZ	VPT2	48
CF3^+	D 3h		CCSD(T)/VTZ	VPT2	48
AlF3	D 3h		CCSD(T)/VTZ	CVPT4	49
SiF3^+	D 3h		CCSD(T)/VTZ	CVPT4	49
SO3	D 3h	50, 51	CCSD(T)/aV(Q+d)Z	VPT2	52
H2CO	C 2v	53, 54	CCSD(T)/aVQZ	ICH (TROVE)	55
H2SiO	C 2v		CCSD(T)/aV5Z	ICH (RVIB4)	56
(C)
HOOH	C 2	57–61	CCSD(T)-F12/aV7Z	ICH (RVIB4)	62
DOOD	C 2	63	CCSD(T)-F12/aV7Z	ICH (RVIB4)	62
HSOH	C 1	64–66	CCSD(T)/aV(Q+d)Z	ICH (TROVE)	67
(D)
t-HNNH	C 2h	68	CCSD(T)/VTZ	VPT2	69
c-HSiOH	C s		CCSD(T)/aV5Z	VPT2	70
t-HSiOH	C s		CCSD(T)/aV5Z	VPT2	70
c-DSiOD	C s		CCSD(T)/aV5Z	VPT2	70
t-DSiOD	C s		CCSD(T)/aV5Z	VPT2	70
c-HOCO	C s		CCSD(T)/CBS	VCI	71
t-HOCO	C s	72, 73	CCSD(T)/CBS	VCI	74, 75
c-DOCO	C s		CCSD(T)/CBS	VCI	71
t-DOCO	C s	73, 76	CCSD(T)/CBS	VCI	74, 75

---

Fig. 2 Topologies of molecules contained within our benchmark set.²⁸

Full dimensional sextic force fields in rectilinear normal mode coordinates (V_ref) and their truncated analogues (Ṽ_harm, V_drop and Ṽ_drop) are obtained via analytical coordinate transformation from curvilinear normal mode quartic force fields, following the procedures outlined in Section 2.2. The force constants that define these CNM QFFs are themselves obtained by analytical differentiation of the reference PES implemented within PyPES.²⁸

Anharmonic fundamental wavenumbers are computed using VCI at excitation levels (≡ maximum sum of vibrational quanta distributed amongst harmonic oscillator basis states) up to and including VCI(10). A complete set of computed fundamentals is reported in Table S1 of the ESI,† along with benchmark computational and experimental data compiled from the literature.

3.2 Applications – pentatomics and hexatomics

For pentatomic and hexatomic molecules, CNM QFFs are obtained via exact, analytic coordinate transformation from RNM QFFs, whose force constants are computed by numerical differentiation, with second derivative data supplied ab initio. In principle, differentiation could be carried out in curvilinear normal mode coordinates directly but in practice it is easier to make the required atomic displacements in rectilinear space, provided that the computational cost of generating the full QFF is not prohibitive. Full dimensional sextic force fields (V_ref) and mode-dropped analogues (V_drop and Ṽ_drop) are generated following exactly the same procedures as used in the benchmarking section.

Because the computed fundamentals will be compared directly to experiment, it is important to ensure that the computed RNM QFF is as accurate as possible, to avoid the choice of electronic structure model being the dominant source of error. For this reason, we employ a “hybrid” force field approach, in which the entire PES is computed at a high level of theory for which analytic Hessians are available (fc-CCSD(T)/aVTZ) and the harmonic force constants replaced with values obtained at an even higher level of theory (fc-CCSD(T)-F12a/VTZ-F12). This strategy has been widely used and extensively validated in the vibrational spectroscopy literature.^9,83–89 It provides an optimal balance between accuracy and computational cost, and typically reduces electronic structure model errors to <5 cm⁻¹.^9,83–89

Frozen-core (fc) CCSD(T) (coupled cluster singles doubles perturbative triples) analytic Hessian calculations⁹⁰ are carried out in an augmented triple-zeta basis^91,92 (fc-CCSD(T)/aVTZ) using the Cfour program package,^93,94 version 1. Higher-level fc-CCSD(T)-F12a⁹⁵/VTZ-F12⁹⁶ harmonic force constants are computed using Molpro^97,98 version 2018.1.

Anharmonic fundamental wavenumbers are computed at excitation levels up to and including VCI(9). A complete set of computed fundamentals is reported in Tables S4–S6 of the ESI,† along with benchmark computational and experimental data compiled from the literature.

4 Benchmarking – tetratomics

The overall aim of this section is to establish how excluding or constraining potentially problematic torsion and inversion modes within the Watson Hamiltonian affects computed fundamental transition frequencies of the remaining modes.

4.1 Force field stability

First, it is necessary to establish whether force field expansions can be stabilised by selective truncation along large-amplitude modes, either in curvilinear or rectilinear normal mode coordinates. As stressed by Császár and coworkers in previous works,^20,30 truncated force field expansions do not generally have the correct long-range behaviour for large-amplitude vibrational modes and thus should not be employed “for studies on highly excited (ro)vibrational states.”²⁰ However, by using carefully selected sets of internal coordinates and systematic PES construction procedures, it may be possible to ensure correct asymptotic behaviour in all regions of interest. Force field stability is relatively straightforward to assess, by monitoring convergence of computed anharmonic fundamental transition wavenumbers as a function of VCI excitation level.

Uncertainties in computed VCI(10) fundamental transition wavenumbers are computed as the difference between values obtained at VCI(9) and VCI(10):

δ_cvge(ν) = ν⁽¹⁰⁾ − ν⁽⁹⁾ (10)

where ν = ν_ref, [small nu, Greek, tilde]_harm, ν_drop or [small nu, Greek, tilde]_drop, and the superscript refers to the VCI excitation level at which each set of fundamentals was computed. VCI convergence uncertainties for all bending and stretching fundamentals are illustrated in Fig. 3. Torsion and inversion modes are excluded in all cases to ensure a fair comparison can be made between models.

Fig. 3 Signed convergence uncertainties δ_cvge in computed stretching and bending VCI fundamentals for all molecules in our benchmarking data set, excluding torsion and inversion fundamentals in all cases. Each box extends from the lower quartile to the upper quartile of the data, and whiskers extend out to twice the interquartile range. Outliers are marked with crosses (+) and each box is bisected by a line indicating the median.

In agreement with previous work,²⁴ν_ref values converge to within 1 cm⁻¹ for most molecules, except those with large-amplitude torsion or inversion modes: NH₃, DOOD and all non-deuterated molecules with large-amplitude torsions, i.e. HSiOH and HOCO, HSOH, and HOOH. The low-barrier, large-amplitude nature of these motions is evident spectroscopically, where tunnelling splittings arise due to wavefunction delocalisation over multiple minima. The largest errors occur in modes that can couple to these problematic modes by symmetry, where errors can propagate from poorly-described torsion/inversion modes into stretching and bending modes through mode coupling within the wavefunction.

The outliers in Fig. 3 (|δ_cvge(ν)| > 5 cm⁻¹) arise due to resonance mixing with excited inversion/torsion states that cause state assignments to alternate as a function of VCI excitation level. Two examples of this behaviour are shown in Table 3, and the remainder reported in Table S2 of the ESI.† The most extreme example is the OSH bending fundamental of HSOH, which couples strongly to the first torsional overtone via Fermi resonance. Similarly, a large convergence error occurs in the OH stretching fundamental of t-HSiOH, where state assignments vary due to resonance mixing to the seventh torsional overtone.

Table 3 From V_ref computed transition wavenumbers (cm⁻¹) that become near-degenerate with excited torsional states nν₆ as a function of VCI excitation level and percentage wavefunction contributions P (squared VCI coefficients) from respective VCI basis functions. States are automatically grouped according to their leading wavefunction coefficients. All other cases are reported in the ESI, Tables S2 and S3

V ref		VCI(8)	VCI(9)	VCI(10)
HSOH	Fundamental	1002.3	994.8	1016.3
	P(ν4)	87%	48%	52%
	P(2ν6)	8%	39%	37%
	Resonant	1040.4	1016.1	992.9
	P(ν4)	9%	40%	44%
	P(2ν6)	72%	24%	44%
t-HSiOH	Fundamental	3666.3	3651.0	3665.9
	P(ν1)	89	42%	57%
	P(3ν3 + ν4)	<1	<1%	2%
	P(8ν6)	<1	21%	<1%
	Resonant		3675.7	3661.2
	P(ν1)		25%	28%
	P(3ν3 +ν4)		20%	13%
	P(8ν6)		11%	<1%

---

Constraining torsion or inversion modes to remain harmonic in curvilinear normal mode coordinates (ν_ref → [small nu, Greek, tilde]_harm in Fig. 3) decreases the frequency of outliers, but does not otherwise substantially alter the error distribution. The single outlier corresponds to the OH stretching fundamental of hydrogen peroxide, where a resonance mixing that does not directly involve the torsional mode causes state assignments to alternate as a function of VCI excitation level (see ESI,† Table S3).

From the trialled truncation schemes, the most robust and reliable way to ensure rapid and monotonic VCI convergence is to drop problematic modes entirely. This not only eliminates outliers that arise due to resonance mixing with excited torsion/inversion states, but also reduces all convergence uncertainties to <0.3 cm⁻¹.

Taken together, the results presented in Fig. 3 indicate that force field instabilities are primarily due to using rectilinear coordinates to describe intrinsically curvilinear motions, rather than the accuracy of the force field expansion per se.

If anharmonicity along torsion or inversion modes were the dominant cause of force field instability, then uncertainty profiles for δ_cvge([small nu, Greek, tilde]_harm) would be similar to those of δ_cvge(ν_drop) or δ_cvge([small nu, Greek, tilde]_drop), because in both cases the energy is effectively constrained to change harmonically along the problematic mode. However, Fig. 3 shows that these uncertainty profiles are quite different, indicating significant off-diagonal coordinate-mediated coupling in Ṽ_harm that is rigorously absent from V_drop and Ṽ_drop. Fundamentally, this arises from the fact that potential energy surface expansions converge slowly in rectilinear normal mode coordinates, particularly when applied to intrinsically curvilinear motions such as torsions and inversions. This limitation in principle affects all PES models based upon rectilinear normal mode coordinates, including numerical m-mode representations.

Conversely, the fact that the δ_cvge([small nu, Greek, tilde]_ref) and δ_cvge([small nu, Greek, tilde]_harm) error distributions are similar (excluding resonance-induced outliers) indicates that both models have equally well-behaved PES expansions along the large-amplitude coordinate. This suggests that our force field expansions have inherited appropriate long-range behaviour from the primitive set of internal coordinates from which our curvilinear normal mode quartic force fields are constructed.

Overall, dropping problematic modes is the most efficient and robust strategy for stabilising force field expansions, but also represents the most extreme approximation considered here. Thus, it is important to characterise how PES truncation schemes affect the accuracy of computed fundamentals.

4.2 Effect of PES truncation schemes

Errors in computed anharmonic fundamental transition wavenumbers due to force field truncation are computed relative to their fully-coupled all-mode counterparts:

δ_ref(ν) = ν⁽¹⁰⁾ − ν⁽¹⁰⁾_ref (11)

where ν = [small nu, Greek, tilde]_harm, ν_drop or [small nu, Greek, tilde]_drop, and the superscript refers to the sum of quanta used in the VCI calculation. ν_ref reference results cannot be converged for HOOH, DOOD, HSOH, HSiOH, HOCO (see ESI,† Table S1) and so these molecules are excluded from further analysis in this section. For NH₃, although VCI(10) results diverge, VCI(9) results are monotonically converged to within 1 cm⁻¹ so are used instead. Boxplot analysis of errors aggregated by molecular topology are illustrated in Fig. 4, while 5 illustrates error distributions sorted according to vibrational mode type.

Fig. 4 Signed errors δ_ref (see eqn (11)) in bending and stretching fundamentals due to applying harmonic constraints on or dropping the torsion/inversion are illustrated in boxplot format, aggregated according to molecular topology for each approximate force field model. Each box extends from the lower quartile to the upper quartile of the data, and whiskers extend out to twice the interquartile range. Outliers are marked with crosses (+) and each box is bisected by a line indicating the median.

Fig. 5 Signed errors δ_ref (see eqn (11)) in bending and stretching fundamentals due to applying harmonic constraints on or dropping the torsion/inversion are illustrated in boxplot format, aggregated according to type of vibration (bending, high frequency X–H and other stretches) for each approximate force field model. Each box extends from the lower quartile to the upper quartile of the data, and whiskers extend out to twice the interquartile range. Outliers are marked with crosses (+) and each box is bisected by a line indicating the median.

Fig. 4 reveals that error distributions are largely independent of molecule type, with all truncated models yielding maximum absolute errors <10 cm⁻¹ in computed fundamentals for all modes except the labelled outliers. In these cases, strong Fermi resonances cause state assignments to vary between V_ref, Ṽ_harm and V_drop/Ṽ_drop models, according to the VCI coefficient data shown in Table 4. However, situations like this can be readily detected by inspection of VCI coefficients, without recourse to reference calculations, and appropriate error bars applied. We suggest applying error bars that stretch between the predicted wavenumbers of the two near-degenerate VCI eigenstates. These outliers will not be included in any further analysis.

Table 4 Computed VCI(10) wavenumbers (cm⁻¹) and percentage wavefunction contributions P (squared VCI coefficients) from respective VCI basis functions for the near-degenerate ν₄ (COD bend) and 2ν₆ (first overtone of torsion) states in c-DOCO, and near-degenerate ν₂ (SiD stretch) and ν₃ + ν₅ (combination of SiO stretch and bend) states in c-DSiOD. VCI eigenstates are grouped according to their leading basis state wavefunction coefficients

VCI(10)		V ref	Ṽ harm	Ṽ drop
c-DOCO	Fundamental	969.9	944.1	960.1
	P(ν4)	62%	54%	96%
	P(2ν6)	31%	40%	—
	Resonant	933.7	982.3	—
	P(ν4)	54%	42%	—
	P(2ν6)	40%	53%	—
c-DSiOD	Fundamental	1375.0	1362.3	1374.8
	P(ν2)	50%	55%	55%
	P(ν3 + ν5)	41%	36%	36%
	Resonant	1359.0	1378.5	1358.9
	P(ν2)	40%	35%	35%
	P(ν3 + ν5)	50%	57%	57%

---

Although Fig. 4 suggests that all three truncated models are of approximately equal quality, there are subtle differences in their error distributions that are investigated in more detail by reanalysing the data according to vibrational mode type. Fig. 5 shows that the Ṽ_harm model tends to produce randomly distributed errors that are similar across all mode types. Mode-dropped models, on the other hand, systematically overestimate stretching wavenumbers, particularly in X–H stretching modes (X = B, C, N, O), which are generally overestimated by 5–10 cm⁻¹. However, mode-dropped models are particularly accurate for bending modes, with randomly distributed errors always within ±5 cm⁻¹. This suggests that neglect of geometric coupling between the dropped large-amplitude mode (e.g. XH torsion) and retained stretching mode (XH stretch) is one of the main accuracy-determining factors of torsion-dropped models.

4.3 Coordinate definitions & coordinate-mediated coupling

V _drop and Ṽ_drop are particularly attractive models because they are conceptually simple, fully stabilise PES expansions in the relevant energy regime, and/but avoid Fermi resonances between excited torsion/inversion modes and other fundamentals. However, differences and similarities between these truncation schemes are not immediately clear. Inspection of Fig. 3 and 4 even creates the impression that they are the same. Indeed, they often yield very similar results and, in many cases, completely identical results. Because Ṽ_drop has the potential to reduce computational cost if curvilinear normal mode force constants are computed ab initio, it is desirable to exploit this similarity where possible.

Identical results can only be obtained if the dropped mode is completely uncoupled to other modes geometrically, because the chances of complete energetic decoupling are practically zero. For planar molecules, the large-amplitude mode is the only out-of-plane mode, so is geometrically orthogonal to all other modes by symmetry. Therefore, it is trivial to exclude the inversion (group B) or torsional (group D) internal coordinate and drop the corresponding vibrational mode without affecting any of the remaining coordinate definitions or vibrational wavenumbers. In these cases, V_drop and Ṽ_drop are rigorously identical. For non-planar molecules, non-trivial mappings between internal valence coordinates and curvilinear normal mode coordinates may arise, encoded within the R matrix that defines curvilinear normal mode coordinates as linear combinations of internals (see eqn (4)). Example R matrices for all molecular topologies are provided in the ESI,† Section S2 and computed fundamentals in Table S1.

For example, in ammonia (C_3v), the symmetric stretching coordinate acquires substantial bending character after the inversion mode is dropped from CNM space. This is an artefact of redundancy in the internal coordinate basis when it is used to span a reduced CNM space, which leads to a deviation of 5.8 cm⁻¹ between ν_drop and [small nu, Greek, tilde]_drop for the symmetric stretching fundamental. In this case, V_drop appears to be the more reliable model, because this redundancy problem is avoided.

For hydrogen peroxide (C₂), although the torsion can couple to three other modes by symmetry, there is a close correspondence between the character of the internal and CNM coordinates, i.e. the torsional mode is defined almost exclusively by displacement along the torsional coordinate. Therefore, excluding the torsional internal coordinate (to ensure numerical stability) and dropping the torsional mode in CNM coordinates prior to coordinate transformation yields almost exactly the same RNM force field as obtained by dropping the torsional mode post-transformation. Hence, the fundamentals computed from V_drop and Ṽ_drop agree to within 0.1 cm⁻¹.

Overall, it is safe to drop modes in CNM space if they are either rigorously orthogonal to other modes by symmetry or there is a (near) one-to-one mapping between the CNM mode to be dropped and a corresponding internal coordinate to be excluded from the internal coordinate basis. Mathematically, this is equivalent to stating that the R matrix must be exactly or very near block-diagonal, with the mode(s) to be dropped contained within its own block. We note that this approach trivially generalises to allow multiple mode-dropping in CNM coordinates, provided that appropriate sets of internal and CNM coordinates to be excluded can be identified separate as sub-blocks of the R matrix.

4.4 Comparison with computational literature

While the internal validation procedures described above have enabled us to assess truncation errors for molecules whose V_ref force fields are otherwise stable enough to compute converged fundamental wavefunctions and energies, we have not been able to directly assess accuracy of truncated models for molecules with large-amplitude torsional or inversion modes. However, it is exactly these problematic cases that are most interesting.

Fortunately, for many molecules in our data set, alternative reference values are available in the literature, obtained by using more complete and/or accurate representations of the potential energy surface and more appropriate and/or accurate methods for solving the nuclear vibrational Schrödinger equation. Cases in which literature values were obtained using second-order vibrational perturbation theory (VPT2) will be omitted from any further analysis.

Differences between computed VCI(10) fundamentals and literature reference values are defined as:

δ_lit(ν) = ν⁽¹⁰⁾ − ν_lit, (12)

where ν = ν_ref, [small nu, Greek, tilde]_harm, ν_drop or [small nu, Greek, tilde]_drop. Mean absolute deviations for each molecular topology are illustrated in Fig. 6.

Fig. 6 Mean absolute deviations 〈|δ_lit|〉 (see eqn (12)) between computed fundamentals and literature reference data, grouped by molecular topology for each different force field model used in this work. Molecules in groups C and D are partitioned into hydrogenated (-h) and deuterated (-d) isotopologues.

For trigonal pyramidal and trigonal planar molecules (groups A and B), force field truncation errors generally increase with the extent of truncation:

|ν_ref − ν_lit| < |[small nu, Greek, tilde]_harm − ν_lit| < |ν_drop − ν_lit| ≈ |[small nu, Greek, tilde]_drop − ν_lit|

For trigonal planar molecules (group B), this is primarily a consequence of neglecting energetic coupling between the dropped inversion mode and other vibrational modes, because the inversion coordinate is geometrically decoupled from the others. For trigonal pyramidal molecules (group A), additional errors arise due to mixing between inversion and symmetric stretching coordinates, with the largest errors for molecules with hydrogen substituents that undergo the largest amplitude vibrations.

For hydrogenated molecules with rotatable bonds (groups C and D), ν_drop and [small nu, Greek, tilde]_drop are universally more accurate approximations to ν_lit than either ν_ref or [small nu, Greek, tilde]_harm. This is because dropped-mode models rigorously eliminate instabilities in the PES both along the torsional mode and by removing coordinate-mediated coupling between the torsional mode and the remaining modes of interest. They have the added advantage of simplifying computed vibrational spectra by also eliminating artificial or mis-assigned resonances involving torsional modes, as found for ν_ref and [small nu, Greek, tilde]_harm.

Hydrogen peroxide represents a pathological case for Watson Hamiltonian models, with a particularly low-barrier torsional mode, energetic and geometric coupling between modes, and a number of (Fermi) resonances between vibrational states.^57–61 As such, all models exhibit large deviations from literature reference values.⁶² Nonetheless, even in this case, the V_drop and Ṽ_drop models are the “best of the worst”, which implies that ignoring energetic coupling between modes is a reasonable trade-off to make to avoid the pitfalls listed above. On the whole, Fig. 6 shows that ν_ref predictions tend to be about as accurate as [small nu, Greek, tilde]_harm, ν_drop and [small nu, Greek, tilde]_drop.

5 Applications – pentatomics & hexatomics

Having verified that dropped-mode models are accurate to within 10 cm⁻¹ for stretching and bending modes of tetratomics, it remains to determine whether these models are overall accurate and generalisable enough to reliably assign gas phase vibrational spectra for larger molecules. As a first test of extensibility, we choose pentatomic and hexatomic molecules that all have a torsion vibration and for which full sets of fundamental vibration wavenumbers have been obtained from high-resolution spectroscopy: methyleneimine CH₂NH,^99–103 hydroxylamine NH₂OH,¹⁰⁴ formaldoxime CH₂NOH,^105–108 and methanol CH₃OH^109–119 (see insets in Fig. 8). Computed fundamental wavenumbers and reference experimental values are reported in the ESI.†

Methanol is a useful test case, because it exhibits particularly strong coupling between its low-barrier, large-amplitude torsional mode and other vibrational degrees of freedom. This is reflected in substantial tunnelling splittings (up to 20 cm⁻¹, and inversion of sign) into states of A and E symmetry, as illustrated in Fig. 7.^{110–115,117–121} As a textbook example of such behaviour, methanol has been the subject of intensive experimental^109–119 and computational^13,120–126 investigations.

Fig. 7 Experimental^{110,112–115,117–119} and computed^120,121 tunneling splittings (Δ = ν(E) − ν(A) + Δ_g.s.) for all vibrational modes of methanol except the torsion. Raw data is listed in Table S5 of the ESI.†

Fig. 8 Error to experiment (δ_expt) for computed variational anharmonic fundamentals (this work and selected literature^120,121,133) of methyleneimine (CH₂NH),^99–103 hydroxylamine (NH₂OH),¹⁰⁴ formaldoxime (CH₂NOH)^105–108 and methanol (CH₃OH).^{110,113–115,117–119} Irreducible representations are encoded in the font, A′ normal and A′′ in italics. For methanol, computational and experimental data are averaged (E : A, 2 : 1). Raw data is listed in Tables S4–S6 of the ESI.†

At the other extreme, methyleneimine, is relatively rigid and rotationally restricted by the C=N double bond. This makes it particularly amenable to computational interrogation, and several perturbative^127–131 and variational^132–134 anharmonic frequency calculations are reported in the literature.

Hydroxylamine and formaldoxime are intermediate in flexibility, capable of rotation about N–O single bonds. For hydroxylamine, perturbative^135,136 and sub-space¹⁰⁴ anharmonic frequency calculations have been performed. To the best of our knowledge, there are no anharmonic frequency studies on formaldoxime in the literature. Experimentally, no cis–trans delocalisation effects or tunneling splittings have been observed for either molecule.

5.1 Internal validation

As in the benchmarking section above, we first assess force field stability by monitoring VCI convergence:

δ_cvge(ν) = ν⁽⁹⁾ − ν⁽⁸⁾ (13)

Consistent with previous computational results,¹³³ the full-dimensional SFF (V_ref) of methyleneimine is stable at all VCI levels up to and including VCI(9), and all fundamentals converge to within 0.1 cm⁻¹. However, for all other molecules, divergences at higher VCI levels mean that reference full-dimensional results cannot be obtained. Therefore, we turn to dropped-mode models V_drop and Ṽ_drop.

In all cases, fundamentals converge to within 0.1 cm⁻¹ with the exception of the symmetric NH₂ stretch ν₂ in NH₂OH, which undergoes resonance mixing with the first overtone of the NH₂ scissor 2ν₃ and the second overtone of the NH₂ wag 3ν₅. Even in this case, fundamental wavenumbers converge to within 2 cm⁻¹ for both V_drop and Ṽ_drop models.

More problematic is the fact that this resonance mixing leads to alternating state assignments between V_drop and Ṽ_drop models, as evidenced by the VCI mixing coefficients presented in Table 5. This is consistent with previous observations that resonance mixing is highly sensitive to even small changes in the underlying PES model, but that this can be readily detected by inspection of VCI coefficients.

Table 5 Computed VCI(9) wavenumbers (cm⁻¹) and percentage wavefunction contributions P (squared VCI coefficients) for ν₂ and resonant state in NH₂OH. VCI eigenstates are grouped according to their leading basis state wavefunction coefficients

VCI(9)		V drop	Ṽ drop
NH2OH	Fundamental	3291.2	3303.9
	P(ν2)	45%	50%
	P(3ν5)	21%	19%
	P(2ν3)	3%	9%
	Resonant	3309.5	3282.9
	P(ν2)	27%	21%
	P(3ν5)	37%	41%
	P(2ν3)	6%	<1%

---

Aside from the ν₂ mode of NH₂OH, comparison of the V_drop and Ṽ_drop results in Fig. 8 reveals that these models otherwise agree to within 2 cm⁻¹ for all other fundamentals, and are rigorously identical for modes that are orthogonal to the dropped torsional coordinate by symmetry. The largest errors occur in the CH₂ wag of CH₂NH and the NH₂ wag of NH₂OH, with deviations of 1.6 cm⁻¹ and 1.2 cm⁻¹, respectively – vibrations for which the curvilinear normal coordinate definitions change substantially upon dropping the torsional coordinate from the primitive internal basis.

5.2 External validation

Errors to experiment for all four molecules are illustrated in Fig. 8, and corresponding numerical data provided in Tables S4–S6 of the ESI.† Where alternative high quality computational predictions are available in the literature,^120,121,133 these are also compared to experiment in Fig. 8. In order to compare the single-reference PyVCI results for methanol with experimental data, the low-resolution band centre is computed by averaging E and A fundamental transitions two to one.

For methyleneimine (top panel), all models agree well with experiment and each other, except for the NH stretch fundamental ν₁, which the dropped-mode models overestimate by ∼26 cm⁻¹. This is well outside the maximum +10 cm⁻¹ overestimation error expected from the tetratomic benchmark. In this case, our full-dimensional data and reference computational results by Rauhut et al.¹³³ suggest that this disagreement does not stem entirely from the torsion-dropping procedure, because these more rigorous models also significantly overestimate the predicted NH stretching fundamental, albeit to a lesser extent. Inspection of the VCI coefficients reported in Table S5 of the ESI,† reveals a Fermi resonance between ν₁ (V_ref: 3280.7 cm⁻¹) and the CN stretching overtone 2ν₄ (V_ref: 3260.0 cm⁻¹), despite the fact that only a single band centre (3262.6 cm⁻¹) is observed experimentally, with very weak rovibrational infrared transitions.¹⁰¹ This raises the possibility that the 2ν₄ overtone has been misassigned as the ν₁ fundamental. We encourage experimental reinvestigation of this spectral region, perhaps using Raman spectroscopy, because NH stretches tend to be very Raman-active.

For hydroxylamine and formaldoxime (middle two panels), dropped-mode models comfortably reproduce experimental fundamentals to within 10 cm⁻¹ for all modes except the O–H stretches. This even includes the NH₂ stretching vibration of hydroxylamine that is strongly resonance-coupled to overtones of its NH₂ scissor and NH₂ wag. Indeed, in most cases, errors to experiment are within ±5 cm⁻¹. O–H stretches are overestimated by ∼15 cm⁻¹, consistent with previous observations from benchmarking tetratomics.

In methanol (bottom panel), the O–H stretching error-to-experiment is unexpectedly small, at ∼2 cm⁻¹. However, this is almost certainly due to fortuituous error cancellation. Previous benchmark computational studies using full-dimensional electronic structure models for methanol with PES of approximately the same quality underestimate the O–H stretching fundamental by 5–10 cm⁻¹.^120,121 If this represents the error due to choice of electronic structure model used in constructing the PES, it suggests that mode-dropping incur an error of 7–12 cm⁻¹, cancelling to give the observed 2 cm⁻¹. The other obvious outlier in this data set is the COH bending mode ν₆ that dropped-mode models overestimate by ∼17.5 cm⁻¹. The fact that this fundamental is well modelled in previous full-dimensional studies^120,121 suggests that this error is due to mode-dropping. This is supported by the spectroscopic analysis of Lees et al., who report strong resonance mixing between ν₆ and torsional combitone states in methanol (ν₇ + ν₁₂ and ν₈ + ν₁₂).¹¹⁹ It is therefore unsurprising that the missing resonance with torsional combitone states leads to larger deviations for the ν₆ fundamental.

6 Conclusions

Three different approaches have been trialled to stabilise rectilinear normal mode force field expansions of the potential energy surface for molecules with large-amplitude torsion or inversion modes which can be performed using our freely available, open-source PyPES and PyVCI program packages.

Constraining the problematic mode to remain harmonic and uncoupled in curvilinear normal mode coordinates prior to transformation is not completely effective in stabilising the resultant sextic force field in rectilinear normal mode coordinates. This shows that problems with torsion and inversion modes are not just due to the shape of the potential energy surface along these coordinates, but are also an inherent problem that arises when curvilinear coordinates are expanded in a finite-order rectilinear normal mode basis. We note that this finding applies to all schemes in which potential energy surfaces are expanded or mapped out in rectilinear coordinates.

Removing large-amplitude torsions completely stabilises PES expansions in all cases, introducing errors of up to +15 cm⁻¹ in coupled stretching modes and absolute errors in other fundamentals of ∼5 cm⁻¹. Dropping problematic modes following coordinate transformation (V_drop) is always numerically robust, but comes with no computational gain during PES construction. Although dropping modes prior to coordinate transformation (Ṽ_drop) has the potential to reduce the computational cost of computing the requisite force constants ab initio, care must be taken to exclude matching coordinates from the internal coordinate basis or modify the basis to ensure that the process remains numerically stable. More generally, the R matrices that define curvilinear normal mode coordinates as linear combinations of internal coordinates provide a useful diagnostic of geometric separability of vibrational modes. We anticipate that this could be further utilised to break the “curse of dimensionality” associated with simulating IR spectra more generally by decomposing the full dimensional nuclear vibrational problem into a series of weakly coupled lower dimensional sub-problems.

The rigorous benchmarking procedures carried out in this work confirm the general applicability of the mode-dropping approach. Introduced errors in the remaining fundamentals are approximately commensurate with, but not necessarily additional to, errors associated with using CCSD(T) theory with at least a triple-zeta basis to construct “spectroscopically accurate” PES. Overall, the combined electronic and nuclear vibrational model error when comparing predicted fundamentals directly to experiment is around 5 cm⁻¹ on average but can be up to 15 cm⁻¹ for X–H stretches that couple geometrically to the dropped mode or slightly higher for modes that couple energetically to the dropped mode.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. Nejad thanks the German Academic Scholarship Foundation for funding a six month research stay in Christchurch New Zealand, the Fonds der Chemischen Industrie (FCI) for a generous scholarship, and Martin A. Suhm, Ricardo A. Mata, Benjamin Schröder for many valuable discussions. D. L. Crittenden thanks the DFG and BENCh RTG 2455 for funding a research visit to Germany under the Mercator Fellowship scheme. Funding by the Deutsche Forschungsgemeinschaft (DFG, project numbers 389479699/GRK2455 and 405832858) is gratefully acknowledged.

Notes and references

1. Over the last decade, the average citation of the paper by Scott and Radom on harmonically scaled harmonic wavenumbers [ A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502–16513 ] was 250 per year, according to a Web of Science search.
2. C. Emmeluth, M. A. Suhm and D. Luckhaus, J. Chem. Phys., 2003, 118, 2242–2255.
3. E. L. Sibert III, J. T. Hynes and W. P. Reinhardt, J. Phys. Chem., 1983, 87, 2032–2037.
4. J. K. G. Watson, Mol. Phys., 1968, 15, 479–490.
5. J. K. G. Watson, Mol. Phys., 1970, 19, 465–487.
6. C. Puzzarini, J. Bloino, N. Tasinato and V. Barone, Chem. Rev., 2019, 119, 8131–8191.
7. A. G. Császár, C. Fábri, T. Szidarovszky, E. Mátyus, T. Furtenbacher and G. Czáko, Phys. Chem. Chem. Phys., 2012, 14, 1085–1106.
8. J. M. Bowman, T. Carrington and H.-D. Meyer, Mol. Phys., 2008, 106, 2145–2182.
9. M. Neff and G. Rauhut, Spectrochim. Acta, Part A, 2014, 119, 100–106.
10. J. T. Hougen, P. R. Bunker and J. W. C. Johns, J. Mol. Spectrosc., 1970, 34, 136–172.
11. W. H. Miller, N. C. Handy and J. E. Adams, J. Chem. Phys., 1980, 72, 99–112.
12. T. Carrington Jr. and W. H. Miller, J. Chem. Phys., 1984, 81, 3942–3950.
13. D. P. Tew, N. C. Handy, S. Carter, S. Irle and J. Bowman, Mol. Phys., 2003, 101, 3513–3525.
14. D. P. Tew, N. C. Handy and S. Carter, J. Chem. Phys., 2006, 125, 084313.
15. R. A. Mata and M. A. Suhm, Angew. Chem., Int. Ed., 2017, 56, 11011–11018.
16. M. Sibaev and D. L. Crittenden, J. Chem. Phys., 2016, 144, 214107.
17. K. Machida, J. Chem. Phys., 1966, 44, 4186–4194.
18. A. R. Hoy, I. M. Mills and G. Strey, Mol. Phys., 1972, 24, 1265–1290.
19. W. D. Allen and A. G. Császár, J. Chem. Phys., 1993, 98, 2983–3015.
20. A. G. Császár, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 273–289.
21. O. Christiansen, Phys. Chem. Chem. Phys., 2007, 9, 2942–2953.
22. J. M. Bowman, G. Czakó and B. Fu, Phys. Chem. Chem. Phys., 2011, 13, 8094–8111.
23. J. Li, K. Song and J. Behler, Phys. Chem. Chem. Phys., 2019, 21, 9672–9682.
24. M. Sibaev and D. L. Crittenden, Comput. Phys. Commun., 2016, 203, 290–297.
25. R. E. Tuzun, D. W. Noid and B. G. Sumpter, Macromol. Theory Simul., 1996, 5, 771–788.
26. G. Simons, R. G. Parr and J. M. Finlan, J. Chem. Phys., 1973, 59, 3229–3234.
27. E. B. Wilson Jr., J. C. Decius and P. C. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, Dover Publications Inc., New York, 1980.
28. M. Sibaev and D. L. Crittenden, J. Comput. Chem., 2015, 36, 2200–2207.
29. W. D. Allen, A. G. Császár, V. Szalay and I. M. Mills, Mol. Phys., 1996, 89, 1213–1221.
30. G. Czakó, T. Furtenbacher, A. G. Császár and V. Szalay, Mol. Phys., 2004, 102, 2411–2423.
31. J. H. Meal and S. R. Polo, J. Chem. Phys., 1956, 24, 1119–1125.
32. J. H. Meal and S. R. Polo, J. Chem. Phys., 1956, 24, 1126–1133.
33. Š. Urban, V. Špirko, D. Papoušek, J. Kauppinen, S. P. Belov, L. I. Gershtein and A. F. Krupnov, J. Mol. Spectrosc., 1981, 88, 274–292.
34. C. Cottaz, I. Kleiner, G. Tarrago, L. R. Brown, J. S. Margolis, R. L. Poynter, H. M. Pickett, T. Fouchet, P. Drossart and E. Lellouch, J. Mol. Spectrosc., 2000, 203, 285–309.
35. G. Guelachvili, A. H. Abdullah, N. Tu, K. N. Rao, Š. Urban and D. Papoušek, J. Mol. Spectrosc., 1989, 133, 345–364.
36. I. Kleiner, L. R. Brown, G. Tarrago, Q.-L. Kou, N. Picqué, G. Guelachvili, V. Dana and J.-Y. Mandin, J. Mol. Spectrosc., 1999, 193, 46–71.
37. S. N. Yurchenko, J. Zheng, H. Lin, P. Jensen and W. Thiel, J. Chem. Phys., 2005, 123, 134308.
38. O. N. Ulenikov, E. S. Bekhtereva, V. A. Kozinskaia, J.-J. Zheng, S.-G. He, S.-M. Hu, Q.-S. Zhu, C. Leroy and L. Pluchart, J. Mol. Spectrosc., 2002, 215, 295–308.
39. A. Ainetschian, U. Häring, G. Spiegl and W. A. Kreiner, J. Mol. Spectrosc., 1997, 181, 99–107.
40. R. I. Ovsyannikov, W. Thiel, S. N. Yurchenko, M. Carvajal and P. Jensen, J. Chem. Phys., 2008, 129, 044309.
41. L. Fusina and G. Di Lonardo, J. Mol. Spectrosc., 2002, 216, 493–500.
42. S. N. Yurchenko, W. Thiel and P. Jensen, J. Mol. Spectrosc., 2006, 240, 174–187.
43. W. Jerzembeck, H. Bürger, J. Breidung and W. Thiel, J. Mol. Spectrosc., 2004, 226, 32–44.
44. K. Aarset, A. G. Császár, E. L. Sibert III, W. D. Allen, H. F. Schaefer III, W. Klopper and J. Noga, J. Chem. Phys., 2000, 112, 4053–4063.
45. S. Yamamoto, R. Kuwabara, M. Takami and K. Kuchitsu, J. Mol. Spectrosc., 1986, 115, 333–352.
46. S. Yamamoto, K. Kuchitsu, T. Nakanaga, H. Takeo, C. Matsumura and M. Takami, J. Chem. Phys., 1986, 84, 6027–6033.
47. S. G. W. Ginn, D. Johansen and J. Overend, J. Mol. Spectrosc., 1970, 36, 448–463.
48. Y. Pak and R. C. Woods, J. Chem. Phys., 1997, 106, 6424–6429.
49. Y. Pak, E. L. Sibert III and R. C. Woods, J. Chem. Phys., 1997, 107, 1717–1724.
50. J. Ortigoso, R. Escribano and A. G. Maki, J. Mol. Spectrosc., 1989, 138, 602–613.
51. N. F. Henfrey and B. A. Thrush, Chem. Phys. Lett., 1983, 102, 135–138.
52. J. M. L. Martin, Spectrochim. Acta, Part A, 1999, 55, 709–718.
53. F. K. Tchana, A. Perrin and N. Lacome, J. Mol. Spectrosc., 2007, 245, 141–144.
54. A. Perrin, A. Valentin and L. Daumont, J. Mol. Struct., 2006, 780, 28–44.
55. A. Yachmenev, S. N. Yurchenko, P. Jensen and W. Thiel, J. Chem. Phys., 2011, 134, 244307.
56. J. Koput, S. Carter and N. C. Handy, Chem. Phys. Lett., 1999, 301, 1–9.
57. W. B. Olson, R. H. Hunt, B. W. Young, A. G. Maki and J. W. Brault, J. Mol. Spectrosc., 1988, 127, 12–34.
58. A. Perrin, A. Valentin, J.-M. Flaud, C. Camypeyret, L. Schriver, A. Schriver and P. Arcas, J. Mol. Spectrosc., 1995, 171, 358–373.
59. C. Camy-Peyret, J.-M. Flaud, J. W. C. Johns and M. Noel, J. Mol. Spectrosc., 1992, 155, 84–104.
60. W. B. Cook, R. H. Hunt, W. N. Shelton and F. A. Flaherty, J. Mol. Spectrosc., 1995, 171, 91–112.
61. J.-M. Flaud, C. Camy-Peyret, J. W. C. Johns and B. Carli, J. Chem. Phys., 1989, 91, 1504–1510.
62. P. Małyszek and J. Koput, J. Comput. Chem., 2013, 34, 337–345.
63. J.-M. Flaud, J. W. C. Johns, Z. Lu, G. Winnewisser and H. Klein, Can. J. Phys., 2001, 79, 367–374.
64. H. Beckers, S. Esser, T. Metzroth, M. Behnke, H. Willner, J. Gauss and J. Hahn, Chem. – Eur. J., 2006, 12, 832–844.
65. O. Baum, HSOH: an elusive species with many different traits, PhD thesis, Cuvillier, Göttingen, 1st edn, 2008.
66. O. Baum, T. F. Giesen and S. Schlemmer, J. Mol. Spectrosc., 2008, 247, 25–29.
67. S. N. Yurchenko, A. Yachmenev, W. Thiel, O. Baum, T. F. Giesen, V. V. Melnikov and P. Jensen, J. Mol. Spectrosc., 2009, 257, 57–65.
68. F. Hegelund, H. Burger and O. Polanz, J. Mol. Spectrosc., 1994, 167, 1–10.
69. J. M. L. Martin and P. R. Taylor, Spectrochim. Acta, Part A, 1997, 53, 1039–1050.
70. J. M. L. Martin, J. Phys. Chem. A, 1998, 102, 1394–1404.
71. R. C. Fortenberry, X. Huang, J. S. Francisco, T. D. Crawford and T. J. Lee, J. Chem. Phys., 2011, 135, 214303.
72. J. T. Petty and C. B. Moore, J. Mol. Spectrosc., 1993, 161, 149–156.
73. T. J. Sears, W. M. Fawzy and P. M. Johnson, J. Chem. Phys., 1992, 97, 3996–4007.
74. R. C. Fortenberry, X. Huang, J. S. Francisco, T. D. Crawford and T. J. Lee, J. Chem. Phys., 2011, 135, 134301.
75. X. Huang, R. C. Fortenberry, Y. Wang, J. S. Francisco, T. D. Crawford, J. M. Bowman and T. J. Lee, J. Phys. Chem. A, 2013, 117, 6932–6939.
76. J. T. Petty and C. B. Moore, J. Chem. Phys., 1993, 99, 47–55.
77. S. N. Yurchenko, W. Thiel and P. Jensen, J. Mol. Spectrosc., 2007, 245, 126–140.
78. S. M. Colwell, S. Carter and N. C. Handy, Mol. Phys., 2003, 101, 523–544.
79. E. L. Sibert III, Comput. Phys. Commun., 1988, 51, 149–160.
80. E. L. Sibert III, J. Chem. Phys., 1988, 88, 4378–4390.
81. J. F. Gaw, A. Willetts, W. H. Green and N. C. Handy, Advances in molecular vibrations and collision dynamics, Greenwich, CT, 1991, p. 169.
82. J. M. Bowman, S. Carter and X. Huang, Int. Rev. Phys. Chem., 2003, 22, 533–549.
83. A. D. Boese and J. M. L. Martin, J. Phys. Chem. A, 2004, 108, 3085–3096.
84. N. Gohaud, D. Begue and C. Pouchan, Int. J. Quantum Chem., 2005, 104, 773–781.
85. K. Pflüger, M. Paulus, S. Jagiella, T. Burkert and G. Rauhut, Theor. Chem. Acc., 2005, 114, 327–332.
86. K. Yagi, S. Hirata and K. Hirao, Theor. Chem. Acc., 2007, 118, 681–691.
87. M. Sparta, M. B. Hansen, E. Matito, D. Toffoli and O. Christiansen, J. Chem. Theory Comput., 2010, 6, 3162–3175.
88. C. Puzzarini, M. Biczysko and V. Barone, J. Chem. Theory Comput., 2011, 7, 3702–3710.
89. S. V. Krasnoshchekov, R. S. Schutski, N. C. Craig, M. Sibaev and D. L. Crittenden, J. Chem. Phys., 2018, 148, 084102.
90. J. Gauss and J. F. Stanton, Chem. Phys. Lett., 1997, 276, 70–77.
91. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007–1023.
92. R. A. Kendall, T. H. Dunning Jr. and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
93. J. F. Stanton, J. Gauss, L. Cheng, M. E. Harding, D. A. Matthews and P. G. Szalay, CFOUR, Coupled-Cluster techniques for Computational Chemistry, a quantum-chemical program package, For the current version, see https://www.cfour.de.
94. D. A. Matthews, L. Cheng, M. E. Harding, F. Lipparini, S. Stopkowicz, T.-C. Jagau, P. G. Szalay, J. Gauss and J. F. Stanton, J. Chem. Phys., 2020, 152, 214108.
95. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
96. K. A. Peterson, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2008, 128, 084102.
97. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
98. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al., MOLPRO, version 2018.1, a package of ab initio programs, 2018, see https://www.molpro.net.
99. G. Duxbury, H. Kato and M. L. Le Lerre, Faraday Discuss. Chem. Soc., 1981, 71, 97–110.
100. G. Duxbury and M. L. Le Lerre, J. Mol. Spectrosc., 1982, 92, 326–348.
101. L. Halonen and G. Duxbury, Chem. Phys. Lett., 1985, 118, 246–251.
102. L. Halonen and G. Duxbury, J. Chem. Phys., 1985, 83, 2091–2096.
103. L. Halonen and G. Duxbury, J. Chem. Phys., 1985, 83, 2078–2090.
104. D. Luckhaus, J. Chem. Phys., 1997, 106, 8409–8426.
105. G. Duxbury, R. M. Percival, D. Devoy and M. R. M. Mahmoud, J. Mol. Spectrosc., 1988, 132, 380–392.
106. G. Duxbury, J. Mol. Spectrosc., 1988, 132, 393–406.
107. R. A. Bannai and G. Duxbury, J. Opt. Soc. Am. B, 1994, 11, 170–175.
108. R. A. Bannai, G. Duxbury, G. Ritchie and S. Klee, J. Mol. Spectrosc., 1996, 178, 84–92.
109. I. Kleiner, G. T. Fraser, J. T. Hougen and A. S. Pine, J. Mol. Spectrosc., 1991, 147, 155–172.
110. R. H. Hunt, W. N. Shelton, W. B. Cook, O. N. Bignall, J. W. Mirick and F. A. Flaherty, J. Mol. Spectrosc., 1991, 149, 252–256.
111. L. Xu and J. Hougen, J. Mol. Spectrosc., 1995, 173, 540–551.
112. G. Moruzzi, B. P. Winnewisser, M. Winnewisser, I. Mukhopadhyay and F. Strumia, Microwave, infrared, and laser transitions of methanol: atlas of assigned lines from 0 to 1258 cm⁻¹, CRC Press, Boca Raton, London, New York, 1995.
113. L.-H. Xu, X. Wang, T. J. Cronin, D. S. Perry, G. T. Fraser and A. S. Pine, J. Mol. Spectrosc., 1997, 185, 158–172.
114. X. Wang and D. S. Perry, J. Chem. Phys., 1998, 109, 10795–10805.
115. R. H. Hunt, W. N. Shelton, F. A. Flaherty and W. B. Cook, J. Mol. Spectrosc., 1998, 192, 277–293.
116. R. M. Lees and L.-H. Xu, Phys. Rev. Lett., 2000, 84, 3815–3818.
117. R. M. Lees, M. Mollabashi, L.-H. Xu, M. Lock and B. P. Winnewisser, Phys. Rev. A, 2002, 65, 042511.
118. M. A. Temsamani, L.-H. Xu and R. M. Lees, J. Mol. Spectrosc., 2003, 218, 220–234.
119. R. M. Lees, L.-H. Xu, J. W. C. Johns, Z.-F. Lu, B. P. Winnewisser, M. Lock and R. L. Sams, J. Mol. Spectrosc., 2004, 228, 528–543.
120. E. L. Sibert III and J. Castillo-Chará, J. Chem. Phys., 2005, 122, 194306.
121. J. M. Bowman, X. Huang, N. C. Handy and S. Carter, J. Phys. Chem. A, 2007, 111, 7317–7321.
122. A. Miani, V. Hänninen, M. Horn and L. Halonen, Mol. Phys., 2000, 98, 1737–1748.
123. V. Hänninen and L. Halonen, Mol. Phys., 2003, 101, 2907–2916.
124. B. Fehrensen, D. Luckhaus, M. Quack, M. Willeke and T. R. Rizzo, J. Chem. Phys., 2003, 119, 5534–5544.
125. J. Castillo-Chará and E. L. Sibert III, J. Chem. Phys., 2003, 119, 11671–11681.
126. Y. Scribano, D. M. Lauvergnat and D. M. Benoit, J. Chem. Phys., 2010, 133, 094103.
127. K. Zaki, M. Gélizé-Duvignau and C. Pouchan, J. Chim. Phys., 1997, 94, 37–53.
128. G. de Oliveira, J. M. L. Martin, I. K. C. Silwal and J. F. Liebman, J. Comput. Chem., 2001, 22, 1297–1305.
129. P. Carbonniere and V. Barone, Chem. Phys. Lett., 2004, 399, 226–229.
130. P. Carbonniere and V. Barone, Chem. Phys. Lett., 2004, 392, 365–371.
131. L. Margulès, J. Demaison, P. Sreeja and J.-C. Guillemin, J. Mol. Spectrosc., 2006, 238, 234–240.
132. C. Pouchan and K. Zaki, J. Chem. Phys., 1997, 107, 342–345.
133. G. Rauhut, G. Knizia and H.-J. Werner, J. Chem. Phys., 2009, 130, 054105.
134. S. Heislbetz and G. Rauhut, J. Chem. Phys., 2010, 132, 124102.
135. A. T. Kowal, Spectrochim. Acta, Part A, 2002, 58, 1055–1067.
136. M. Kumarasiri, C. Swalina and S. Hammes-Schiffer, J. Phys. Chem. B, 2007, 111, 4653–4658.

---

Footnote

† Electronic supplementary information (ESI) available: Computational details, computed wavenumbers for all investigated molecules, R matrix definitions, and computational and available gas phase experimental literature data. See DOI: 10.1039/d0cp03515g

---