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

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 aﬀects 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 eﬀective and eﬃcient strategy is to remove the mode in question from the PES expansion entirely. This introduces errors of up to +10 cm (cid:2) 1 in stretching fundamentals that would otherwise couple to the dropped mode, and (cid:3) 5 cm (cid:2) 1 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.


Introduction
Harmonic normal mode analysis and scaled harmonic frequencies continue to enjoy great popularity in assisting the assignment of vibrational infrared and Raman spectra, 1 but can fail to accurately model anharmonically coupled systems, e.g. the OH stretching spectrum of carboxylic acid dimers. 2 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.'' 3 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 semirigidity 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.'' 6 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. 7 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 3 O + and NH 3 , 8,9 but cannot be extended to other types of motions, such as torsions. 9 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, 10 reaction path Hamiltonian, 11 reaction surface Hamiltonian 12 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. 15 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. 16 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.

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][18][19][20] numerical mapping on multi-modal grids, 8,21 and symmetryinvariant 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 4 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 gridbased 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. 19 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. 16 Constraints that restrict how force fields describe largeamplitude 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. 16

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 are where r, y and t are the usual bond length, bond angle and dihedral coordinates, 25 sin(o) is an out-of-plane umbrella coordinate, 25 and f SPF (r) is the Simons-Parr-Finlan radial coordinate 26 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 À1 ) of Cartesian displacements (X) that diagonalise the mass-weighted Hessian matrix: 27 Curvilinear normal mode coordinates (Q) are defined analogously: 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. 28 Our procedures for computing the required force constants, either analytically or numerically, are summarised in Fig. 1.
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 À1 in computed fundamental wavenumbers, but usually they are 1-2 cm À1 . 16

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 theoryprovide 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. 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.
(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.

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 4 for nonlinear molecules in rectilinear normal mode coordinates, including first-order Coriolis coupling between rotations and vibrations (in atomic units): where B a e is the equilibrium rotational constant about axis a and z 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 f n i are single-mode harmonic oscillator functions, and n is a string of quantum numbers n 1 ,. . .,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 P i n i n.
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.

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 semiglobal potential energy surfaces are available in the literature and implemented within the PyPES program package. 28 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.
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. 28 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.

Applications -pentatomics and hexatomics
For pentatomic and hexatomic molecules, CNM QFFs are obtained via exact, analytic coordinate transformation from Table 1 Shorthand notation for potential energy surface expansions constructed in this work. The tilde superscript denotes that modedropping 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 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][84][85][86][87][88][89] It provides an optimal balance between accuracy and computational cost, and typically reduces electronic structure model errors to o5 cm À1 . 9,[83][84][85][86][87][88][89] Frozen-core (fc) CCSD(T) (coupled cluster singles doubles perturbative triples) analytic Hessian calculations 90  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.

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.

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.'' 20 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): where n = n ref ,ñ harm , n drop orñ 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.
In agreement with previous work, 24 n ref values converge to within 1 cm À1 for most molecules, except those with largeamplitude torsion or inversion modes: NH 3 , DOOD and all nondeuterated 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 (|d cvge (n)| 4 5 cm À1 ) 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.
Constraining torsion or inversion modes to remain harmonic in curvilinear normal mode coordinates (n ref -ñ 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 o0.3 cm À1 .
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 d cvge (ñ harm ) would be similar to those of d cvge (n drop ) or d cvge (ñ 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 coordinatemediated 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 d cvge (ñ ref ) and d cvge (ñ 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.

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:    Table S1) and so these molecules are excluded from further analysis in this section. For NH 3 , although VCI(10) results diverge, VCI(9) results are monotonically converged to within 1 cm À1 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 reveals that error distributions are largely independent of molecule type, with all truncated models yielding maximum absolute errors o10 cm À1 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.
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 À1 . However, mode-dropped models are particularly accurate for bending modes, with randomly distributed errors always within AE5 cm À1 . 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.

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 (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.  (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. Table 4 Computed VCI(10) wavenumbers (cm À1 ) and percentage wavefunction contributions P (squared VCI coefficients) from respective VCI basis functions for the near-degenerate n 4 (COD bend) and 2n 6 (first overtone of torsion) states in c-DOCO, and near-degenerate n 2 (SiD stretch) and n 3 + n 5 (combination of SiO stretch and bend) states in c-DSiOD. VCI eigenstates are grouped according to their leading basis state wavefunction coefficients 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-ofplane 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.

and 4
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 À1 between n drop andñ 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 2 ), 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 posttransformation. Hence, the fundamentals computed from V drop and Ṽ drop agree to within 0.1 cm À1 .
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 blockdiagonal, 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.

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: where n = n ref ,ñ harm , n drop orñ drop . Mean absolute deviations for each molecular topology are illustrated in Fig. 6. For trigonal pyramidal and trigonal planar molecules (groups A and B), force field truncation errors generally increase with the extent of truncation: 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. Fig. 6 Mean absolute deviations h|d lit |i (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 hydrogenated molecules with rotatable bonds (groups C and D), n drop andñ drop are universally more accurate approximations to n lit than either n ref orñ harm . This is because droppedmode 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 n ref andñ 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][58][59][60][61] As such, all models exhibit large deviations from literature reference values. 62 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 n ref predictions tend to be about as accurate asñ harm , n drop andñ drop .

Internal validation
As in the benchmarking section above, we first assess force field stability by monitoring VCI convergence: d cvge (n) = n (9) À n (8) Consistent with previous computational results, 133 the fulldimensional 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 À1 . However, for all other molecules, divergences at higher VCI levels mean that reference fulldimensional 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 À1 with the exception of the symmetric NH 2 stretch n 2 in NH 2 OH, which undergoes resonance mixing with the first overtone of the NH 2 scissor 2n 3 and the second overtone of the NH 2 wag 3n 5 . Even in this case, fundamental wavenumbers converge to within 2 cm À1 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.
Aside from the n 2 mode of NH 2 OH, comparison of the V drop and Ṽ drop results in Fig. 8 reveals that these models otherwise agree to within 2 cm À1 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 2 wag of CH 2 NH and the NH 2 wag of NH 2 OH, with deviations of 1.6 cm À1 and 1.2 cm À1 , respectivelyvibrations for which the curvilinear normal coordinate definitions  Table S5 of the ESI. † change substantially upon dropping the torsional coordinate from the primitive internal basis.

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 singlereference 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 n 1 , which the dropped-mode models overestimate by B26 cm À1 . This is well outside the maximum +10 cm À1 overestimation error expected from the tetratomic benchmark. In this case, our full-dimensional data and reference computational results by Rauhut et al. 133 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 n 1 (V ref : 3280.7 cm À1 ) and the CN stretching overtone 2n 4 (V ref : 3260.0 cm À1 ), despite the fact that only a single band centre (3262.6 cm À1 ) is observed experimentally, with very weak rovibrational infrared transitions. 101 This raises the possibility that the 2n 4 overtone has been misassigned as the n 1 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 À1 for all modes except the O-H stretches. This even includes the NH 2 stretching vibration of hydroxylamine that is strongly resonance-coupled to overtones of its NH 2 scissor and NH 2 wag. Indeed, in most cases, errors to experiment are within AE5 cm À1 . O-H stretches are overestimated by B15 cm À1 , consistent with previous observations from benchmarking tetratomics.
In methanol (bottom panel), the O-H stretching error-toexperiment is unexpectedly small, at B2 cm À1 . 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 À1 . 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 À1 , cancelling to give the observed 2 cm À1 . The other obvious outlier in this data set is the COH bending mode n 6 that dropped-mode models overestimate by B17.5 cm À1 . 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 n 6 and torsional combitone states in methanol (n 7 + n 12 and n 8 + n 12 ). 119 It is therefore unsurprising that the missing resonance with torsional combitone states leads to larger deviations for the n 6 fundamental.

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 À1 in coupled stretching modes and absolute errors in other fundamentals of B5 cm À1 . 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 À1 on average but can be up to 15 cm À1 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.