General Perturb-Then-Diagonalize Model for the Vibrational Frequencies and Intensities of Molecules Belonging to Abelian and Non-Abelian Symmetry Groups

In this paper, we show that the standard second-order vibrational perturbation theory (VPT2) for Abelian groups can be used also for non-Abelian groups without employing specific equations for two- or threefold degenerate vibrations but rather handling in the proper way all the degeneracy issues and deriving the peculiar spectroscopic signatures of non-Abelian groups (e.g., -doubling) by a posteriori transformations of the eigenfunctions. Comparison with the results of previous conventional implementations shows a perfect agreement for the vibrational energies of linear and symmetric tops, thus paving the route to the transparent extension of the equations already available for asymmetric tops to the energies of spherical tops and the infrared and Raman intensities of molecules belonging to non-Abelian symmetry groups. The whole procedure has been implemented in our general engine for vibro-rotational computations beyond the rigid rotor/harmonic oscillator model and has been validated on a number of test cases.


INTRODUCTION
The reliability of quantum chemical (QC) models to support experimental findings is related from one side to their accuracy and from the other side to their feasibility, robustness, and ease of use. 1,2 Concerning the accuracy, electronic structure computations of energies, geometries, and force fields are nowadays able to rival high-resolution spectroscopy for small systems and to help assignments and interpretation of all kinds of spectra for larger molecular systems, provided that nuclear motions and environmental effects are taken into proper account. 3−8 In the present contribution, we will be concerned with molecular vibrations, which are directly sampled by different conventional [infrared (IR), Raman] and chiral (VCD, ROA) spectroscopies, 2,9 but tune also the outcomes of other spectroscopies (e.g., distortion constants in microwave spectroscopy 4 or line shapes in electronic spectroscopies 10 ). Until quite recently, except for exceedingly small systems, the rigid rotor/harmonic oscillator model was nearly exclusively employed to describe molecular vibrations. However, the neglect of anharmonicity and ro-vibrational couplings can introduce significant errors, sometimes leading to even qualitatively wrong interpretations of experimental data.
Among the different approaches available to go beyond the rigid rotor/harmonic oscillator approximation, 11−35 those based on perturbation theory applied to the expansion of the nuclear Hamiltonian in the power series of products of vibrational and rotational operators (hereafter referred to as vibrational perturbation theory, VPT) are particularly appealing for their remarkable cost/performance ratio, at least for semi-rigid molecular systems. Moreover, some formulations of VPT, such as the Van Vleck contact transformation method, 36 fully justify a generalized model (GVPT2), 37,38 allowing to couple the advantages of perturbative (for weakly coupled modes) and variational (for strongly coupled modes) treatments in a wellsound and robust framework. Actually, the GVPT2 approach belongs to the class of perturb-then-diagonalize many-body models, which, although less widely used than diagonalize-thenperturb models, have in the present context some appealing advantages from both conceptual and implementative points of view. Implementations of VPT2 approaches in general-purpose QC software are now quite widespread, 39−50 although fully automatic and robust implementations of GVPT2 are less common. The situation is different for intensities of IR, Raman, VCD, and ROA spectra, which need to account for both mechanical and electrical/magnetic anharmonicity. To the best of our knowledge, the only general platform allowing GVPT2 computations for all the spectroscopic techniques mentioned above is the one implemented by some of the present authors in the Gaussian package. 51 Extension of perturbative/variational procedures to the treatment of flexible systems is underway along different avenues including approaches based on reaction paths or surfaces following the decoupling of one or two largeamplitude motions from a bath of small-amplitude modes and/ or the replacement of Cartesian coordinates by generalized internal coordinates. 52−55 Here, we tackle a different problem, related to the treatment of molecular systems with non-Abelian point group symmetries, namely, linear, symmetric, and spherical tops. As a matter of fact, a significant ensemble of molecular systems, ranging from small to large sizes and of interest in various research fields, belong to these classes, including, for instance, organic and organometallic compounds like coronene and ferrocene 30,56−58 or acetylene derivatives. 59−70 The presence of degenerate modes raises multiple practical and theoretical issues, which have been only marginally addressed until now. A possible workaround is to lift the degeneracies by reducing the symmetry of the molecular system to the closest Abelian group symmetry, 42 but this can lower the accuracy of the results. The rotational problem is actually simpler for non-Abelian groups because the rigid rotor approximation leads to analytical solutions, but a proper account of the degeneracy for the vibrational problem requires a careful check of the relative orientations of the degenerate modes to ensure that the anharmonic force fields and property surfaces are correctly built and an alternative derivation able to account for the couplings involving the degenerate modes. In a previous work, 71 we have presented a complete framework implementing different equations for the vibrational frequencies based on the symmetry for linear and symmetric tops and taking into proper account both intrinsic and accidental degeneracies, leading to additional terms in the Hamiltonian or to singularities in the perturbative expansion, respectively. However, the current situation is unsatisfactory at several levels. First, the formulation derived from the work of Pliva 72 cannot be extended straightforwardly to spherical tops, and to the best of our knowledge, no complete and correct derivation has been proposed. Second, the current implementation is particularly intricate due to the constant case switch depending on the symmetry and degeneracy of the modes in the calculation of the quantities of interest for the vibrational energies. Finally, it requires the use of complex algebra in the variational treatment, which leads to complex eigenvectors of the variational matrix, and thus the transition moments, even if the final intensities remain, of course, real.
Based on these premises, a unified treatment of Abelian and non-Abelian symmetries would be more efficient and general, provided that the results can be easily transformed to the more standard representation for compatibility and/or interpretative purposes. Furthermore, the simplicity of the new approach could allow a more straightforward implementation of new strategies based on VPT2 73 not only for linear and symmetric tops but also for spherical tops. Finally, a general and robust framework can be set for the calculation of both IR and Raman intensities for all point groups without any need of complex algebra. This has convinced us to follow a different route, that is, to employ the asymmetric-top formulation also for the other cases, handling in the proper way all the degeneracy issues and deriving the customary spectroscopic signatures of non-Abelian groups (e.g., -type doubling) by a posteriori transformations of the eigenvectors. As will be shown in the following, the results for frequencies are exactly the same as those delivered by our previous conventional implementation (including resonance contributions), but we are now able to compute also intensities for both IR and Raman spectra.
The paper is organized as follows: In the Theory section, we start with deriving energy expressions for linear and symmetric tops, including resonance-free expressions for the zero-point energy, and then move to a general treatment of resonances and to intensities for all vibrational spectroscopies. The Results and Discussion section is organized in the same way but also considers different levels of electronic structure theory, starting from Hartree−Fock (HF) and second-order Møller−Plesset perturbation theory (MP2), which do not involve any underlying noise connected to numerical integration, and then moving to methods rooted in the density functional theory (DFT) up to double hybrids. The main results, remaining challenges, and perspectives are shortly outlined in the concluding section.
2. THEORY 2.1. Framework. In order to set up the framework for our discussion, let us consider a system of N vibrational normal modes with a non-Abelian symmetry. These modes will generally be identified by the indexes i, j, k, l. Where relevant, and except if specified otherwise, different sets of indexes will be used to distinguish degenerate modes (s, t, u, followed by a subscript number for each mode) from the non-degenerate ones (m, n, o). As an example, modes s 1 and s 2 are two degenerate modes with the same harmonic wavenumber ω s , while m and n are two non-degenerate modes.
The formalism introduced by Plíva for symmetric and linear tops 72 relies on a special set of coordinates for degenerate modes based on a complex combination = ± ± q q iq s s s 1 2 (1) where q i represents the dimensionless normal coordinate associated to mode i. This form can actually be directly related to the definition of the harmonic vibrational wave function ψ. Thus, it is convenient to first recall its representation for a system with N′ non-degenerate modes and N″ sets of degenerate modes (N″ = 0 for asymmetric tops). In the canonical representation (C), ψ, associated to the vibrational state |v⟩, is given as a product of one-dimensional functions (2) where N s ″ is the degeneracy order associated to degenerate mode s, v i is the number of vibrational quanta corresponding to mode i, are the well-known one-dimensional harmonic oscillator wave functions. To make the discussion simpler, we will consider only molecules with at most doubly degenerate modes, that is, linear and symmetric tops. Extension to spherical tops is deferred to a dedicated section. Equation 2 can thus be written more explicitly as An alternative way to treat the degenerate coordinates is through the polar representation, 25 in which the wave function, namely, ψ v P (q), assumes the following form: is the harmonic wave function of the two-dimensional isotropic harmonic oscillator Hamiltonian expressed in terms of the polar coordinates ρ s and θ s arising from a pair of degenerate modes.
It is worth recalling that both representations, canonical and polar, have the same eigenvalues and hence vibrational energies but different eigenstates. From a theoretical perspective, the polar representation is preferable since it leads to an explicit quantization of the vibrational angular momentum stemming from each pair of degenerate vibrations. The formulation proposed by Plíva takes these properties into account through a transformation of the degenerate normal coordinates and the associated force field. 72,74 Having confirmed that the presence of degenerate modes does not preclude the development of VPT2 equations for symmetric and linear tops in the canonical representation (see Appendix A), the main task is to define a suitable transformation between the states obtained in each representation. At the harmonic level, it is possible to build a linear transformation between sets of degenerate states 38 where |ψ v P ⟩ and |ψ v C ⟩ are the column vectors containing the states with the same harmonic energy, which implies in practice that v m and v s = v s 1 + v s 2 are constant, and P v is a unitary matrix connecting the two sets of states. The complete set of rotation matrices required for the conversion of fundamentals, first overtones, and binary (1 + 1) combinations from the vibrational ground state is reported in Section S1 of the Supporting Information. Furthermore, following the recent extension of our computational framework to the inclusion of three-quanta states, 75 the full set of the corresponding rotations is reported for the sake of completeness.
Therefore, a unified framework can be set up, in which calculations are run in two steps: 1. anharmonic calculations in the canonical framework; 2. if degenerate modes are present, application of the necessary rotations to switch to the polar representation. 2.2. Vibrational Energies. Full derivations for the VPT2 energies of asymmetric 11,32,42 and symmetric or linear 71,72,74 tops have been reported elsewhere. Here, we will focus on the novel aspects using refs 48,71,76 as a basis. Based on Plíva's formalism, the vibrational energy for symmetric and linear tops can be written where ε 0 P is the zero-point vibrational energy (ZPVE) in the polar representation, d i is the degeneration of the ith normal mode, and χ P and g are respectively the anharmonic constant matrix in polar representation and the matrix containing the anharmonic contributions from the angular momenta, given in ref 71 and reported in Section S2 of the Supporting Information.
Let us consider a set of degenerate harmonic states in the canonical (|ψ v C ⟩) and polar (|ψ v P ⟩) representations. The associated blocks of the contact-transformed vibrational Hamiltonian are respectively H̃v ,v C and H̃v ,v P , where the subscript "v" indicates that only matrix elements between states with the same principal quantum numbers v are included in the block. These blocks will be simply referred to as diagonal blocks in the following. By construction, only the diagonal elements of the corresponding matrices contribute to the anharmonic energies. Using eq 5, the following identity can be writteñ Since the trace of a matrix is invariant under similarity transformation, we have Equation 7 can thus be systematically used to obtain the energies in the polar representation, starting from both energies and resonances (see later) evaluated within the canonical representation.
After demonstrating that the expression of the resonance-free ZPVE for asymmetric tops 77−79 can also be used for linear and symmetric tops (see Section S3 of the Supporting Information for details and notation), that is, Equation 7 can be applied to calculate the energies in the polar representation for one-and two-quanta states involving at least one degenerate mode. Combining the rotation matrices reported in Section S1 of the Supporting Information with the transformation of the Hamiltonian given in eq 7, it is straightforward to prove that the anharmonic fundamental energies do not vary under the change of representation, a property true for any excited state involving only one degenerate mode excited with a single quantum (v s = 1), for instance, , and so on. Following the procedure outlined for the fundamental states, the rotation-based framework yields the expressions for the polar energies of the first overtone associated to a degenerate mode s Ä as well as those of the binary (1+1) combination bands involving two degenerate modes, namely, s and t Contrary to fundamental bands, the energy terms in polar representation cannot be directly obtained from the canonical representation but require additional, off-diagonal terms, collectively called Darling−Dennison coupling terms, which will be illustrated later. Finally, three-quanta transitions are treated in the Section S1 of the Supporting Information.
2.3. Resonances and the Variational Correction. In the following, we will consider three models for the calculation of VPT2 energies. In the pure VPT2 approach (simply named VPT2), resonances are ignored, in the sense that all terms involved in the calculation of the quantities of interest are systematically included. As a result, this approach has a tendency to break down very quickly with the system size, often leading to unphysical results. An improvement is offered by the deperturbed VPT2 (DVPT2) scheme, where resonances are identified, usually through a multi-step procedure, 42,76,79,80 and the related terms are discarded. Strong discrepancies are prevented, but the selective removal of terms can result in an unbalanced account of the anharmonic contributions, with varying intensity. In the most refined version, the generalized VPT2 (GVPT2) scheme, 32,42,76,81 the resonant terms from DVPT2 are introduced back through an additional variational step, often together with other terms, broadly referred to as Darling−Dennison interaction terms. 82 For this reason, the latter should be systematically considered for spectroscopic applications like those here. The general construction and use of polyads will be reviewed together with GVPT2.
2.3.1. Fermi Resonances and DVPT2 Energies. At the VPT2 level, the energies in the polar representation can be systematically obtained through linear combinations of canonical energies and Darling−Dennison resonances. Let us first discuss the extension of the theoretical framework developed so far to the DVPT2 scheme, where each potentially resonant term in the χ matrix is analyzed by applying specific criteria, and those which are identified as resonant are discarded from the calculation. In the present work, the identification of Fermi resonances is done through a two-step procedure, which first considers the energetic proximity of the interacting states, namely, ω i ≈ 2ω j (type I) and ω i ≈ ω j + ω k (type II), and then the magnitude of the term, using Martin's test 80 to estimate the deviation of the term from the variational energy of a model, ad hoc system. Currently, DVPT2 can be employed for both canonical and polar representations, the only difference lying in the definition of the χ matrix and the construction of the g matrix in the second case. However, it can be demonstrated that the formalism described in Section 2.2 can be easily extended to DPVT2 calculations since eq 7 is valid for both resonant and nonresonant terms separatelỹ= † At this point, the transformation becomes similar to the VPT2 case presented above. An analogous procedure can be employed out to evaluate the χ P and g matrices.
In practice, subtle differences could be observed because of the numerical parameters and tests used to define the resonances. However, for an equivalent set of resonances, the results obtained through eqs 6 and 15 (χ P and g matrices being deprived of Fermi resonances) will converge.
2.3.2. Variational Correction in GVPT2. GVPT2 is built on top of DVPT2 by adding a final step to calculate the anharmonic energies as eigenvalues of a variational matrix, whose diagonal elements are the DVPT2 energies, discussed in the previous section, and the off-diagonal elements represent the corrective terms to Fermi resonances, complemented by Darling− Dennison interactions, evaluated over the basis of the canonical harmonic-oscillator wave functions.
On the premise that our reference will remain the canonical representation, which also fully fits the conditions of the application of the GVPT2 scheme at first, we will discuss the theory underlying the definition of the polar variational matrix and then how a full equivalence can be reached between the two representations.
The calculation of a specific resonant term can be carried out by applying a simple generalization of eq 8 to the off-diagonal block coupling states differing in terms of principal quanta (v ≠ v′), and it can be organized into two steps: Through the symmetry relations between the anharmonic force constants of symmetric and linear tops, it is possible to prove that the corrective terms due to the Fermi resonances are equivalent in the two representations. Hence, the focus in the following will be on the couplings between states, collectively referred to as Darling−Dennison resonances or interactions.
In this context, particular importance is given to the resonances between states for which the condition v = v′ holds, whose interaction generates the off-diagonal elements of the blocks H̃v ,v P . 2.3.3. -Type Doubling in the Polar Representation. -type doubling terms involve overtones of a given degenerate mode or combination bands of two degenerate modes. Amat derived a general rule 24 for the a priori identification of the non-vanishing off-diagonal terms s t s t s t s s t t (18) where the elements defined in eq 16 contribute to the energy only if the order of the principal symmetry axis (n in C n ) is a multiple of 4 and those in eq 18 only if n is even. The expressions of the terms U s ± , R st ± , and S st ± have been first derived by Grenier and Bresson 83,84 and then re-derived in ref 71. In the present derivation, there is no need of developing specific equations for computing the matrix elements defined in eqs 14−16 since they are off-diagonal elements of the diagonal blocksH P 2 ,2 s s andH P 1 1 ,1 1 s t s t (the full equations in terms of canonical quantities are reported in Section S5 of the Supporting Information). As a matter of fact, the calculation of -type terms can be performed concurrently with the conversion of the anharmonic energies through eq 7.
2.3.4. Diagonalization and GVPT2 Energies. In the preceding sections, it has been shown that the blocks of the polar variational matrix (H̃P) can always be expressed in terms of their canonical counterpart. In order to understand the effects of such a connection on the GVPT2 energies, we will first consider a H̃P matrix only containing the -type terms as off-diagonal elements. In this context, the variational problem simplifies to diagonalizing blocks of the typeH P 2 ,2 s s andH P 1 1 ,1 1 s t s t , whose eigenvalues are equal to those of their canonical counterparts. Therefore, the inclusion of -type doubling at the variational level implies the convergence of GVPT2 energies in the two representations.
This result can be easily generalized, given that an equivalent set of resonances is included in both representations, leading to the following identity: where P is a block-diagonal matrix composed of all rotation matrices required for converting the different diagonal blocks H̃v ,v C , and it is itself unitary.
In analogy with the treatment of -type doubling, the invariance of the eigenvalues of a matrix under unitary transformations can be exploited to state that canonical and polar GVPT2 energies converge to the same values. Let us remark that the -type terms are present even in the absence of accidental resonances. Consequently, their inclusion is mandatory in order to reach the convergence of the GVPT2 energies. From a practical point of view, eq 19 prevents any ambiguity connected to the representation choice since the set of GVPT2 energies is unique.
2.4. Transition Moments and Intensities. Starting from the available literature (e.g., the study of Tarrago and coworkers 85 on the transition dipole moments of C 3v -symmetry systems), a general computational framework to calculate both IR intensities and Raman activities of molecular systems with non-Abelian symmetries has been devised and implemented in our platform.
Let us start from the band intensities and the associated transition moments of linear and symmetric tops taking into account that thanks to the transformation shown in eq 5, it is possible to use the formulas obtained for asymmetric tops (and reported in Section S6 of the Supporting Information) to obtain their counterparts in the polar representation.
If we look at the fundamental bands (the full equation is reported in Sections S6.1.1 and S6.2.1 of the Supporting Information for degenerate modes), terms of the form Ä where P j collects the Cartesian components of the first derivative of the property P with respect to the jth dimensionless normal coordinate, will present a singularity whenever j = s 2 . For this reason, it is necessary to exclude those modes in the summation so that the degenerate modes are assumed to be resonant, and the "resonant" form (Section S6.2.1 in the Supporting Information) is used instead. For simplicity, only transitions from the ground state (noted 0) to a given final state f are considered, noted "0;f". The quantities of interest here for IR and Raman spectroscopies are the dipole strength and Raman activity, labeled in the canonical representation D 0;f C and S 0;f C , respectively, Here, the invariants {a 2 } 0;f C (isotropic) and {γ 2 } 0;f C (anisotropic) for the most general case of a complex tensor are defined as 86,87 α α 0; 0; (22) where τ and η run over the Cartesian axes, while ⟨μ⟩ 0;f C and ⟨α τη ⟩ 0;f C represent respectively the transition integrals of the electric dipole and a component of the polarizability tensor between the canonical states |ψ 0 C ⟩ and |ψ f C ⟩. Let us anticipate that a closure relation having the same form as eq 8 also holds for dipole strengths and Raman activities due to the unitarity of the rotation matrices.
The theoretical framework currently used for the calculation of transition moments can be straightforwardly extended to the polar representation (for more details, see Appendix B). Thus, once the transition dipole moments and polarizabilities are converted through eq 60, the calculation of both IR and Raman intensities in the polar representation is possible. It is worth mentioning that even though the transition moments evaluated in this way are generally complex, the corresponding intensities are always real.
2.4.1. Infrared Intensities and Raman Activities at the VPT2 Level. For readability, the initial-state label will be dropped, so D 0;f will be simply written D f . In analogy with energies, the dipole strengths and Raman activities for fundamental states are the same if they are degenerate so that their contribution to the anharmonic spectra is independent of the representation. This result simplifies considerably the whole conversion procedure since among states with up to two quanta, the only states potentially different with respect to the canonical representation are the sets ψ | ⟩ are respectively the vectors containing the Cartesian components of the transition dipole moments associated with the states ψ | ⟩ The Raman activities can be expressed in a compact notation through the introduction of the variables A r;s C , Γ r;s C , and S r;s (25) which can be interpreted respectively as the "off-diagonal" terms of {a 2 } r C , {γ 2 } r C , and S r C . As a consequence, the expressions of the Raman activities of the polar states of interest are A similar analysis applied to the transition moments of binary combination bands involving degenerate modes yields the dipole strengths of the polar states Finally, the corresponding Raman activities are From a comparison of eqs 25 and 26 with eqs 11a and 11b and eqs 28 and 30 with eqs 13a and 13b, it is possible to observe that the conversion of dipole strengths and Raman activities is ruled by expressions similar to those employed for the anharmonic energies.

Introduction of the Variational Correction.
It has been demonstrated that if the variational matrix H̃P in the polar representation can be expressed by a rotation of the canonical one H̃C (see eq 19), the anharmonic energies within the two representations converge to the same values. This result can be easily extended to vibrational intensities. Indeed, while the definition of variational states is dependent on the representation of the reference states, polar and canonical variational states are equivalent when projected onto the same basis. In this context, the canonical basis is chosen as a reference.
As a matter of fact, the variational states, and hence the intensities, do not depend on the representation. Let us remark that the equality of the variational states holds even when only -type doubling terms are included as off-diagonal elements of the variational matrix, that is, even in the absence of accidental resonances.
In summary, the anharmonic spectrum is completely independent of the representation, with the only difference being the harmonic-state basis. Thus, the computational protocol currently employed for asymmetric tops can be straightforwardly extended to the treatment of symmetric/linear tops without any loss of accuracy.
2.5. Extension to Spherical Tops. In this section, we will show how the framework devised for symmetric and linear tops can be extended to spherical tops. For the sake of concision, the rotation-based formulation will be applied to systems presenting at most threefold degenerate vibrations. Actually, the eigenstates of the three-dimensional isotropic harmonic oscillator are again complex linear combinations of the canonical wave functions. This type of combinations can be extended to higher degeneracy orders, in which case, the following derivation can be straightforwardly adapted to systems exhibiting such characteristics.
Back to our application, the starting point remains the canonical representation since the specificity of the spherical top lies in the definition of the rotation matrices necessary for the transformation of the wave functions. While the canonical wave function is still given by eq 2, the so-called spherical wave function has the following form, where non-degenerate and doubly and triply degenerate modes are gathered in different terms are the solutions of the three-dimensional isotropic harmonic oscillator Hamiltonian expressed with respect to the spherical coordinates stemming from a trio of degenerate modes, r s , γ s , and ϕ s , and the quantum numbers v s′ , k s′ , and m s′ can assume the following values: In analogy with symmetric and linear tops, the analysis of the spherical wave functions enables the definition of the rotation matrices Q v Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article where |ψ v S ⟩ contains the spherical states sharing the same principal angular number v and |ψ v C ⟩ is the corresponding canonical counterpart. Let us stress that when the states |ψ v S ⟩ only involve non-degenerate and doubly degenerate modes, the rotation matrix is P v .
Let us now analyze the effect of rotations in deeper detail, highlighting the analogies with symmetric and linear tops. The canonical expression of ZPVE, whose expression is reported in eq 9, can be equivalently expressed through the following formula: where ψ 0 C represents the harmonic ground-state wave function in the canonical representation. The main advantage of the rotation-based framework is that the specific properties of the rotor are collected in the definition of the wave function, more specifically through the use of rotation matrices to carry out the conversion procedure described above. As a result, the reference basis is always the canonical one, and the operators depending on the normal coordinates and their conjugate momenta are never subject to any modifications. Based on this, the expression of the spherical counterpart ε 0 S is where ψ 0 S is the harmonic ground-state wave function in the spherical representation. Since the ground-state wave function is independent of the representation (see Appendix B for more details), we can conclude that Hence, the resonance-free ZPVE in the spherical representation can be still evaluated through the customary expression, which can be then used for systems presenting both doubly and triply degenerate modes without any restriction.
In the previous sections, all different types of bands involving up to two-quanta excitations characterizing symmetric and linear tops have been derived and analyzed separately. In general terms, this separation is not necessary since only a few rotation matrices are actually sufficient to build all the other ones. With the aim of treating all the states up to two quanta, the largest value that the quantum numbers v m , v s , and v s′ reported in eq 34 can assume is 2 so that the matrices P 1 s , P 2 s , Q 1 s , and Q 2 s allow to express any spherical state in terms of canonical ones (the extension to three-quanta states would only require additional matrices, P 3 s and Q 3 s ).
In the same way as symmetric and linear tops, both transition energies and intensities for states only involving one excited quantum of a threefold degenerate vibration are still independent of the representation, and this is also true for the ZPVE. Concerning overtones and binary combination bands involving triply degenerate modes, a complete derivation of both energies and intensities in the spherical representation for states up to two quanta is reported in Section S7 of the Supporting Information.
Finally, it is worth mentioning that the equivalence of representations at the GVPT2 level is not affected by the presence of triply degenerate vibrations, with it being strictly related to the connection of the wave functions by unitary matrices. Therefore, the GVPT2 results remain independent of the representation.

COMPUTATIONAL DETAILS
The theoretical framework presented in the previous sections has been implemented in a development version of the Gaussian package. 88 Most of the available electronic structure computations allowing analytic computation of second energy and first property derivatives have been employed. These include HF and second-order Møller−Plesset (MP2) wave-function methods together with different flavors of DFT including representative hybrid (B3LYP 89−92 ) and double-hybrid (B2PLYP 93,94 ) exchange−correlation functionals, with the inclusion of empirical dispersion contributions by means of Grimme's D3 model with Becke−Johnson damping 95,96 (hereafter noted B3D3 and B2D3, respectively). As the core of this work concerns the development, functionals with a well-documented reliability on the molecular systems chosen here for illustration purposes were selected. The reliability of the B2PLYP functional in the calculation of both harmonic and anharmonic frequencies has been demonstrated in the literature. 97,98 The B3LYP functional has been used for the calculation of both harmonic and anharmonic frequencies only for linear systems, where B2PLYP results have also been shown. Concerning symmetric and spherical tops, the B3LYP functional has been only used for the calculation of the anharmonic corrections in the hybrid scheme. In this respect, it has been demonstrated that the quality of the harmonic frequencies is much more critical if compared with that of the corresponding anharmonic corrections. 2 The socalled calendar basis sets jun-cc-pVDZ and jun-cc-pVTZ 99 have been consistently employed (referred to in the following as JnDZ and JnTZ, respectively).
The anharmonic data required for the VPT2 calculation of frequencies and intensities have been obtained by finite differences of analytical force constants (full cubic and semidiagonal quartic force constants) and first-order derivatives of the properties (full second and semi-diagonal third derivatives), employing a default displacement of δQ i = 0.01 √amu·Å along each mass-weighted normal coordinate Q i . 42,100 From a practical point of view, the generation of the anharmonic force field and higher-order property derivatives is the most expensive step once the equilibrium structure has been found. Indeed, 2N frequency calculations are needed in addition to the one at the reference geometry, required to generate the displacement vectors. Being independent from one another, they can be run in parallel on separate machines, with the final constants built at the end of the process. Hence, in an optimal scenario, the computational cost can be reduced to roughly twice what is needed at the harmonic level. The VPT2 calculations themselves on systems of this size last a few minutes, independent of the scheme chosen.
For the treatment of resonances, the protocol detailed in ref 76 was used with the default parameters.
In addition to the basic formulation of VPT2, the DVPT2 and GVPT2 schemes are also used for the calculation of both anharmonic energies and intensities. The use of DVPT2 is necessary in the presence of Fermi resonances to avoid the unphysical results issuing from the standard VPT2 equations. In addition to recovering the discarded terms from DVPT2, the GVPT2 scheme allows a straightforward account of Darling− Dennison resonances, which are not explicitly considered within the purely perturbative approach. Besides improving the overall agreement with experimental energies, they can be critical to The so-called hybrid force field scheme has also been employed. 2,102−105 In this approach, harmonic and anharmonic contributions are treated at different levels of theory in view of their different contribution to the final VPT2 result. In particular, anharmonic contributions have been consistently computed with the JnDZ basis set and harmonic terms with JnTZ. In some cases, which will be explicitly mentioned in the discussion, the results at the CCSD(T) level in conjunction with extended basis sets, already available in the literature, were employed for the latter. It should be noted that the higher-level harmonic frequencies are not simply added to anharmonic corrections but are also employed in the conversion of the anharmonic force constants and property derivatives, in the construction of the χ matrix, the definition of the resonant terms, and for the intensity. The JnTZ basis set has been employed for the calculation of anharmonic contributions of aromatic systems since it is well known that out-of-plane vibrations of these molecules are particularly sensitive to the basis-set dimension. 106   Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article

RESULTS AND DISCUSSION
The framework previously discussed has been validated through a series of applications to linear and symmetric tops, illustrative of representative non-Abelian symmetry groups. 4.1. Linear Molecules. First, a set of three-and four-atom molecules, including hydrogen cyanide (HCN), hydrogen isocyanide (HNC), and acetylene (C 2 H 2 ), has been considered. The absence of Fermi and 1−1 Darling−Dennison resonances (between fundamental states) implies that the anharmonic fundamentals do not vary going from VPT2 to DVPT2 or GVPT2 schemes. The fundamental harmonic and anharmonic frequencies obtained with the hybrid force-field model described above are reported in Table 1.
As can be seen from Table 1, the computed fundamental energies are in good agreement with the experimental counterparts, with the largest discrepancy concerning the degenerate mode of HNC. At the B3D3 and B2D3 levels of theory, such an error is most likely due to the underestimation of the corresponding harmonic frequency, which is lower of the experimental value in both cases.
As an example of comparison between the polar and canonical representations, the vibrational frequencies of CO 2 for all states up to two quanta in the polar representation have been calculated and analyzed. The set of wavenumbers in the polar representation is reported in Table 2 and compared with reference experimental values.
For comparison purposes, the VPT2, DVPT2, and GVPT2 wavenumbers in the canonical representation are reported in Table 3.
As expected, the states involving at most one vibrational quantum in the degenerate bending mode have the same frequency irrespective of the chosen representation. Conversely, the frequencies of the first overtones related to the bending mode change between the representations, even though the energies converge to the same values when the GVPT2 model is applied.
By comparing Tables 2 and 3, it is straightforward to verify that the sum of the energies of the states ψ | ⟩ ± P 2 , 2 . Such an outcome is in full agreement with eq 8. Due to the strong Fermi resonance between the overtone of the bending mode and the fundamental of the symmetric stretching (2ω 1 ≈ ω 2 ), the CO 2 molecule has been widely used as a prototype in the study of resonances. Let us underline that only the state ψ | ⟩ , can be conveniently expressed in terms of canonical interaction terms where the expressions of the matrix elements associated to the Fermi resonances in the canonical representation 76 have been used in conjunction with the symmetry rules and the definition of the force constants f mss (σ) (σ = I, III, IV) reported in The results obtained at the B3D3/JnTZ level of theory have been employed in the calculation of both IR and Raman spectra, Table 3. Computed Anharmonic VPT2, DVPT2, and GVPT2 Wavenumbers (in cm −1 ) of CO 2 in the Canonical Representation a The canonical vibrational states are indicated as |v i v j ⟩, and the subscripts "a" and "b" distinguish degenerate modes. b Basis set: JnTZ.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article with the harmonic and anharmonic results obtained within the VPT2, DVPT2, and GVPT2 schemes being compared with the experimental data in Figure 1. As expected, the best agreement between theoretical and experimental spectra is reached within the GVPT2 scheme. This is particularly evident in the Raman spectrum, in which the inclusion of the Fermi resonance discussed above at the variational level leads to an excellent reproduction of the relative intensities characterizing the Fermi diad present in the experimental spectrum. As a last example of linear molecules, dicyanoacetylene (C 4 N 2 , see Figure 2) is considered.
Dicyanoacetylene has been detected on the Titan moon of Saturn by IR spectroscopy, and it is used as a prototype for similar astrochemical systems, such as cyanopolyynes. 120 The IR and Raman spectra of this system have been the object of several experimental works, 121−128 while a theoretical analysis has been recently presented by Dargelos and Pouchan. 120 With the aim of showing the application of our computational framework to longer chain systems, the VPT2, DVPT2, and GVPT2 fundamentals at different levels of theory have been calculated and compared with their experimental counterparts in Table 4.
At both B3D3 and B2D3 levels, two Fermi resonances of the first type coupling the states |1 4 ⟩ and |2 5 ⟩, and |1 3 ⟩ and |2 7 ,0 7 ⟩ are   The best agreement between theoretical and experimental fundamentals is reached at the B2D3 level, despite an out-ofscale discrepancy detected for the degenerate states |1 6 ,±1 6 ⟩ regardless of the electronic level of theory. The corresponding normal mode is depicted in Figure 3.
At the MP2 level, such a difference can be ascribed to an underestimation of the harmonic frequency (which is lower than the experimental value), whereas density functional calculations show an excessive anharmonic correction (229 and 158 cm −1 for B3D3 and B2D3, respectively). Interestingly, the MP2 harmonic energies of the |1 6 ,±1 6 ⟩ fundamentals are very close to those of the experiment, a trend observed for all modes below 1000 cm −1 .
4.2. Symmetric Tops. Shifting to symmetric-top systems, we analyze a set of six molecules, namely, bromotrifluoromethane (CF 3 Br), mono-and tri-deuterated methane (CH 3 D and CHD 3 , respectively), the cyclopentadienyl anion (C 5 H 5 − ), benzene (C 6 H 6 ), and pentaborane (B 5 H 9 ), which are characterized by a principal axis of order ranging from 3 to 6 and belong to the molecular point groups C nv and D nh . As already pointed out, anharmonic corrections at the MP2 level are free from the issues of numerical integration, which is not always true for methods rooted into DFT. In this respect, MP2 is more suitable for validation purposes of the new rotation-based framework, significantly reducing the possibility of errors in the symmetry relations proposed by Amat and Henry, and characterizing the anharmonic force field of symmetric and linear tops. As a matter of fact, while MP2 corrections will be used for all symmetric tops studied in the following, the use of DFT will be limited to the C 3v systems.
Let us start from CF 3 Br (see Figure 4), whose anharmonic IR spectrum has been simulated at the MP2/JnDZ level of theory with harmonic frequencies corrected by CCSD(T)/aug-cc-pVTZ-PP 129 calculations. In the former case, no Fermi resonances are present, while in the latter case, a single Fermi resonance of the second type, involving the non-degenerate fundamental |1 1 ⟩ and combination |1 2 1 3 ⟩, has been detected. In both cases, the only Darling−Dennison resonances detected correspond to the -type doubling of type R (which is the only one present when the principal axis order is 3).
The set of anharmonic fundamental frequencies, together with the energy of state |1 2 1 3 ⟩ (the only two-quanta state involved in a Fermi resonance within the CC//MP2 scheme), is compared with the experimental data in Table 5 where, as expected, the improvement of the results due to the use of coupled-cluster (CC) harmonic frequencies leads to excellent agreement.
The hybrid CC//MP2 results have also been used in the calculation of the anharmonic IR spectrum. Furthermore, in order to show the importance of anharmonic effects in the reproduction of both the position and intensity of the bands of the spectrum, the theoretical spectra have been evaluated by a stepwise inclusion of the anharmonic corrections. More specifically, the IR spectrum has been first evaluated at the purely harmonic level (HH), followed by the inclusion of the anharmonic corrections to the energies (AH) and finally to both energies and intensities (AA). A full comparison of theoretical and experimental spectra is reported in Figure 5.
As can be seen from the top panel of Figure 5, a remarkable improvement in the position of the bands is obtained by correcting the transition energies. Conversely, the bottom panel is characterized by a total absence of theoretical peaks unless the anharmonic contributions to the intensities are included. Indeed, in the spectral window corresponding to the bottom panel, only overtones and combination bands have been detected experimentally, while only fundamental bands have    Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article non-vanishing intensities within the harmonic oscillator model. The anharmonic corrections to the intensities allow not only to bypass this limit but also to produce a spectral profile in full agreement with the experimental one.
With the aim of showing the application of the new framework to the calculation of isotopomers belonging to different point groups, the CH 3 D and CHD 3 molecules (symmetric tops) will be compared to the CH 2 D 2 isotopomer (asymmetric top), whereas the spherical top isotopomers (CH 4 and CD 4 ) will be analyzed in a later section. In this context, four hybrid schemes have been used for computing the wavenumbers of the systems under consideration, namely, the MP2//MP2, B2D3//B3D3, CC//MP2, and CC//B3D3 models in conjunction with the JnDZ (B3D3 and MP2), JnTZ (B2D3), and cc-pVQZ (CC) basis sets. 134 A full comparison of all sets of theoretical data with the corresponding experimental counterparts is reported in Table 6.
Concerning CH 2 D 2 , two Fermi resonances, namely, ω 2 ≈ 2ω 7 and ω 8 ≈ ω 4 + ω 9 , have been found at all levels of calculation, with an additional one (ω 1 ≈ 2ω 3 ) detected at the B2D3//B3D3 level. Conversely, 1−1 and 2−2 Darling−Dennison resonances have not been identified. The resonance analysis carried out for CH 3 D shows the presence of a Fermi resonance of type I (ω 1 ≈ 2ω 5 ) and a 1−1 Darling−Dennison resonance (ω 1 ≈ ω 4 ) at all the computational levels employed here. A second Fermi resonance of type I (ω 2 ≈ 2ω 6 ) is also present at the MP2//MP2 level. Furthermore, a series of 2−2 Darling−Dennison interactions have been detected (including -type terms). Finally, CHD 3 is consistently characterized by two Fermi resonances (ω 2 ≈ 2ω 3 and ω 4 ≈ ω 3 + ω 5 ), with the addition of ω 4 ≈ 2ω 6      Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Lee and co-workers, 134 the state |1 1 ⟩ of CH 3 D is strongly coupled with |2 7 ,0 7 ⟩ (but not with |2 7 ,±2 7 ⟩ since the coupling element vanishes based on eq 38) so that at the GVPT2 level, the states with energies (see Table 6) 2931, 2995, 2978, and 2983 cm −1 for MP2//MP2, B2D3//B3D3, CC//MP2, and CC// B3D3, respectively, and those with energies 3017, 2911, 2890, and 2899 cm −1 are basically equal mixtures of |1 1 ⟩ and |2 7 ,0 7 ⟩. As expected, the results closer to those of the experiment are those based on CC harmonic frequencies, although the B2D3// B3D3 scheme leads to quite satisfactory results. Besides, the results presented here are in good agreement with those reported in ref 134, as confirmed by a comparison of the χ P and g matrices, for which only the diagonal elements are reported in Table 7 for readability.
The next case studies are two planar aromatic systems, namely, the cyclopentadyenil anion and benzene (see Figures 6  and 7), with the former being noted Cp − in the following.
As already anticipated, the presence of out-of-plane vibrations makes the use of a triple-ζ basis set mandatory for the calculation of the anharmonic force field so that the MP2/JnTZ level of theory has been employed for both systems. The assignment of the normal vibrations of the cyclopentadyenil anion has been recently revisited by Bencze and co-workers 137 on the basis of previous studies. 138 As outlined in ref 137, the cyclopentadienyl anion cannot exist without a counter cation under normal Table 7. Comparison of Extrapolated and Computed Anharmonic VPT2 and DVPT2 χ P and g Diagonal Elements (in cm −1 ) of CH 2 D 2 , CH 3 D, and CHD 3   Table 8. The results reported in Table 8 show that the VPT2 results are in good agreement with the experiment and that inclusion of the variational correction improves the agreement. The largest discrepancies concern the |1 12 ⟩ and |1 13 ,±1 13 ⟩ states. There are, of course, systematic shifts related to the difference between the structure of the cyclopentadienyl anion in the experimental complexes and that of the free anion, with the lowest frequencies being, as usual, the most sensitive to environmental effects.
The next molecule studied is benzene (see Figure 7), a D 6h symmetric-top system, which has been extensively studied by both IR and Raman spectroscopies. 30,140−143 The anharmonic calculations have been carried out at the MP2/JnTZ level or coupling the MP2/JnTZ anharmonic contributions to harmonic frequencies evaluated at the CCSD-(T)/ANO4321′ level. 144 The VPT2, DVPT2, and GVPT2 fundamentals of benzene computed at both levels of theory are compared with experimental data in Table 9.
Concerning the GVPT2 scheme, Fermi resonances of type II (ω i ≈ ω j + ω k ) as well as 2−2 DD resonances, were found, with the latter including -type doubling of types R and S. On the other hand, 1−1 DD resonances have not been identified. The use of the hybrid scheme leads to a remarkable improvement of the fundamentals, and this is particularly true for the state |1 14 ⟩, with the corresponding MP2/JnTZ harmonic frequency being clearly overestimated. On top of this, the inclusion of the variational correction further improves the agreement with the experimental data for all the transitions involved in resonances. The IR spectrum evaluated by means of the CC//MP2 hybrid scheme is compared with that of the experiment in Figure 8.
All spectra show a very strong band between 650 and 700 cm −1 corresponding to the fundamental transition to |1 12 ⟩, associated to the out-of-plane bending vibration sketched in Figure 9.
Since such a state is not involved in any resonance, its position is independent of the adopted VPT2 scheme.   Finally, the simulation of the IR spectrum of pentaborane (see Figure 10), a C 4v system, has been carried out.
Among the systems considered in the present study, pentaborane is the only one showing all types of -type doubling within the GVPT2 scheme since the order of its principal axis is a multiple of 4. The simulation has been performed by employing the same hybrid scheme as for the other systems, namely, by coupling MP2/JnDZ anharmonic contributions to MP2/JnTZ harmonic frequencies. A full comparison of the theoretical and experimental fundamentals is reported in Table 10.
At the DVPT2 level, several Fermi resonances (most of type II) have been detected, which have been successively included at the variational level within the GVPT2 scheme, together with the proper -type doubling terms and the other identified 2−2 Darling−Dennison resonances. Most of the frequencies are qualitatively correct and consistent with those reported in a Table 9. Comparison of Experimental and Computed Anharmonic Fundamental VPT2, DVPT2, and GVPT2 Wavenumbers (in    Figure 11, where the theoretical data obtained by means of the VPT2, DVPT2, and GVPT2 schemes are compared. While DVPT2 and GVPT2 spectra show a common pattern, the atypical behavior of the VPT2 spectrum is due to the strong Fermi resonances between the fundamentals |1 22 ,±1 22 ⟩ and the combination bands |1 18 1 27 ,±1 27 ⟩, which lead to a huge value for the transition dipole moment of the |1 22 ,±1 22 ⟩ states so that the corresponding band is the only one clearly visible in the theoretical spectrum. This problem is fixed within the DVPT2 scheme through the elimination of the resonant term, and the result is further refined at the GVPT2 level by the successive variational treatment. 4.3. Paving the Route to Spherical Tops. With the aim of showing the extension of our computational framework to spherical tops, we now analyze a series of systems of both tetrahedral and octahedral symmetries, including, for instance, tetraphosphorus (P 4 ), methane (CH 4 ) and its fully deuterated isotopologue (CD 4 ), and sulfur hexafluoride (SF 6 ), whose structures are sketched in Figure 12.
In the case of linear and symmetric tops, our computational protocol performs by default a full check of all the symmetry relations present between the anharmonic force constants. 71 Such a procedure has not yet been implemented for tetrahedral XY 4 - 136,151 and octahedral XY 6 -like 152 systems so that some slight (and generally negligible) discrepancies can be detected between anharmonic force constants which should be in principle identical. With the aim of limiting this issue as much as possible, the pruned (99,590) grid employed before for DFT calculations (which are in principle the most sensitive to this problem) will be replaced with a larger one (175,974 and 250,974 for first-row atoms and atoms in the second and later rows, respectively).
4.3.1. Tetrahedral Molecules. The fundamental frequencies of the P 4 molecule obtained through the MP2//MP2 and B2D3//B3D3 hybrid schemes are compared with their experimental counterparts in Table 11. The P 4 system does not show any Fermi or 1−1 Darling− Dennison resonance so that the values of the fundamentals do not vary going from VPT2 to DVPT2 or GVPT2 schemes. As can be seen from Table 11, the frequencies of the states |1 1 ⟩ and |1 2 ,±1 2 ⟩ are significantly improved within the B2D3//B3D3 scheme, reaching values within 1 cm −1 from the experimental counterparts. Concerning the triply degenerate states |1 1 ,1 1 ,±1 or 0 1 ⟩, the best estimate is reached by the MP2// MP2 scheme, even though in both cases, the value of the frequency is not as accurate as the previous ones. The origin of such a discrepancy can be traced back to the experimental conditions at which the gas-phase spectrum has been recorded, 153 Table III of ref 154 for more  details).
A strategy often used to deal with systems presenting degeneracies is that of modifying slightly one or more masses/ coordinates in order to lower the symmetry 134,140 and then   Table 10. Comparison of Experimental and Computed (VPT2, DVPT2, and GVPT2) Anharmonic Fundamental Wavenumbers (in cm −1 ) of Pentaborane employ a theoretical model able to treat only lower-or nondegenerate modes. This procedure has been applied to P 4 in order to perform a consistency check of our new implementation. More specifically, the T d symmetry of this system has been gradually reduced without any geometry modification by slightly increasing (by 0.1%) the masses of a single and then a couple of phosphorous atoms, thus obtaining a symmetric (C 3v point group) or an asymmetric (C 2v point group) top, respectively. A comparison of the VPT2 fundamentals of the lower-symmetry systems with those of the fully symmetric molecule is reported in Figure 13.
As expected, triply degenerate states |1 3 ,1 3 ,±1 3 or 0 3 ⟩ are no more present, being replaced with a non-degenerate state of A 1 symmetry and a couple of doubly degenerate states belonging to the irreducible representation E in the C 3v system and with three non-degenerate states with irreducible representations A 1 , B 1 , and B 2 in the C 2v structure. This is clearly visible since the blue bars (C 3v ) corresponding to the states |1 3 ,1 3 ,±1 3 ⟩ are equal to each other but different from that of |1 3 ,1 3 ,1 0 ⟩, while the Anharmonic calculations performed with the JnDZ basis set based on a set of harmonic frequencies evaluated with the JnTZ basis set.
corresponding orange bars (C 2v ) are all different. The couple of states |1 2 ,±1 2 ⟩ belonging to the irreducible representation E are still present in the C 3v system, while they are replaced with two non-degenerate states with representations A 1 and A 2 in the C 2v system. Again, the blue bars of the states |1 2 ,±1 2 ⟩ are equal, while the corresponding orange ones are different, following the symmetry breaking. As a further test, CH 4 and CD 4 are considered. In addition to the MP2//MP2 and B2D3//B3D3 schemes employed for P 4 , two hybrid calculations have been performed at the MP2/JnDZ and B3D3/JnDZ levels in conjunction with a set of harmonic frequencies at the CCSD(T)/cc-pVQZ level. 134 The VPT2 wavenumbers are compared with the experimental data in Table  12.
Both systems do not show any Fermi or 1−1 Darling− Dennison resonance, except for CD 4 within the CC//MP2 hybrid scheme, where a Fermi resonance of type I (ω 1 ≈ 2ω 4 ) has been detected, leading the vibrational frequency of |1 1 ⟩ shifts from 2064 to 2099 cm −1 between the DVPT2 and GVPT2 levels. The presence of 2−2 Darling−Dennison resonances has been detected in the GVPT2 model, and the resonances between states possessing the same principal quanta have been properly included. Concerning CH 4 , the best agreement between theory and experiments is reached with the CC// B3D3 hybrid scheme (MAE = 7 cm −1 ), even if all sets of VPT2 theoretical wavenumbers are close to those of the experiment, with the only exception being the MP2//MP2 scheme (MAE = 28 cm −1 ), for which larger discrepancies are observed.
The results for CD 4 present a trend similar to that of CH 4 , with the MP2//MP2 scheme showing again the worst agreement with the experiment, although it is able to match exactly the experimental value of the fundamental band Figure 13. Errors between anharmonic fundamentals of P 4 and the corresponding counterparts of the C 3v and C 2v symmetry-broken geometries obtained through the MP2//MP2 hybrid scheme. associated to the state |1 1 ⟩. The hybrid schemes based on CC harmonic frequencies (particularly CC//B3D3) are quite close to the experimental data, even though in this case, the best agreement is reached by the B2D3//B3D3 scheme (MAE = 6 cm −1 ).
The anharmonic frequencies and intensities at the CC// B3D3 level have been employed in the simulation of the IR spectrum, which is compared with the harmonic and experimental ones in Figure 14.
As expected, the inclusion of anharmonic effects improves the agreement between theory and experiments, especially for the band around 3000 cm −1 , corresponding to the excitation of the triply degenerate fundamental |1 3 ,1 3 ,±1 3 or 0 3 ⟩. 4.3.2. Octahedral molecules. Finally, sulfur hexafluoride has been chosen as a test case for the O h point group. The anharmonic calculations have been performed at the MP2/ JnDZ and B3D3/JnDZ levels of theory, in conjunction with harmonic frequencies at the MP2/JnTZ and B2D3/JnTZ levels, respectively. A comparison between the VPT2 fundamental wavenumbers and their experimental counterparts is reported in Table 13.
This system does not show any Fermi or 1−1 Darling− Dennison resonance at the levels of theory employed here, while 2−2 Darling−Dennison resonances have been identified and included at the GVPT2 level. The best agreement with the experimental fundamentals is obtained at the MP2//MP2 level (MAE = 11 cm −1 ), while the B2D3//B3D3 scheme is in this case characterized by a significant underestimation of all harmonic wavenumbers. For this reason, the anharmonic results obtained at the MP2//MP2 level have been employed in the simulation of the IR spectrum. Both theoretical harmonic and GVPT2 spectra are compared with the experimental one in Figure 15.
Similar to CH 4 , the inclusion of anharmonic contributions leads to a spectral profile closer to the experimental one, especially for the position of the band around 1000 cm −1 , corresponding to the triply degenerate fundamentals |1 3 ,1 3 ,±1 3 or 0 3 ⟩.

CONCLUSIONS
In this work, we have shown how the canonical representation used for the development of VPT2 equations of asymmetric tops can be extended to linear and symmetric tops, followed by a series of a posteriori transformations, to give results identical to those obtained with the polar representation, thus offering the possibility of an ease of choice of the most convenient form for any application in vibro-rotational and vibrational spectroscopies. Such a strategy offers a number of advantages with respect to previous, ad hoc procedures. The first aspect concerns the ease of implementation since the new approach does not require any heavy modification of the codes already supporting VPT2 for asymmetric tops. The second aspect is the simplicity  of the extension to spherical tops. Once the transformation matrix between the representations is known, it is possible to derive the necessary equations for any quantity of interest, which can be coded in small, specialized routines. However, the most important advantage is the availability of general equations for the intensities of all vibrational spectroscopies without the need of resorting to complex numbers. Actually, to the best of our knowledge, this is the first completely general implementation of intensities in the framework of the double-perturbation theory.
The results show that we dispose now of a general and robust implementation of GVPT2 for Abelian and non-Abelian groups allowing the effective treatment of medium-to large-sized molecules for all electronic structure methods for which analytical Hessians and first derivatives of properties are available. Hybrid methods in which harmonic and anharmonic contributions are treated at different levels can further extend the range of applications of the general platform. Studies in condensed phases can also be performed by means of mixed discrete-continuum models in which the solute and, possibly, some molecules of its cybotactic region are embedded in a polarizable continuum mimicking bulk solvent effects. Also in this case, the availability of analytical Hessians and first derivatives of properties allows an effective GVPT2 treatment.
Of course, all the intrinsic problems of a low-order perturbative treatment based on Cartesian normal modes are still present, especially concerning large-amplitude motions. Besides, the harmonic-oscillator wave functions do not always provide a suitable basis for the representation of vibrations, regardless of both the method used for the inclusion of anharmonic effects and the level of electronic theory employed. However, semi-rigid molecules can be routinely analyzed with remarkable results, largely sufficient for interpretation and assignment tasks. Extension to flexible systems can be pursued by coupling reduced-dimensionality treatments of largeamplitude motions to GVPT2 for small-amplitude motions. In this connection, use of curvilinear in place of rectilinear coordinates is an appealing option. While work in this and related connections is underway in our laboratory, we think that already the present implementation offers a number of interesting perspectives for the study of molecular systems of current scientific and technological interest.

A.1 Degeneracies in the Perturbative Treatment
One of the obvious issues in the VPT2 expansion based on the canonical representation is the risk of singularity arising from the degeneracy of the vibrational modes and the related harmonic states. Let us consider the calculation of the energy ε R of a vibrational state |ψ R ⟩ (the superscript "C" has been dropped for clarity). Following the Rayleigh−Schrodinger perturbation theory, the following relations are obtained: (1) where the vectorial form |v R ⟩ = |v R,1 ...v R,i ...v R,N ⟩ was preferred to represent the harmonic state |ψ R (0) ⟩. The problem lies in the second term of the right-hand side of eq 39c since if |v R ⟩ and |v S ⟩ are degenerate, the denominator would be null. For the sake of simplicity but without loss of generality, let us consider that both states involve the same set of degenerate modes, {s 1 , s 2 }. In order to have ε R (0) = ε S (0) , the total number of quanta v s = v s 1 + v s 2 must be kept constant, and the number of quanta involving non-degenerate modes cannot vary. Mathematically, this translates in an even number of bosonic if |v R ⟩ and |v S ⟩ are degenerate. In other words, the representation matrix of ̂( 1) over the degenerate state basis is null so that the canonical representation can be safely applied to treat systems with degenerate modes.

B.1 Calculation of Properties in the Polar Representation
In the following, we will describe the theoretical framework underlying the calculation of molecular properties (including transition ones) in the polar representation. Let us point out that even if the derivation reported below focuses on symmetric and linear tops, it still holds for systems presenting threefold degenerate vibrations. Let us consider the operator Θ̂describing a molecular property depending on normal coordinates or their conjugate momenta. Following previous works, 34 (0) ; (1) ; (2) ; (1) (1) (0) (0) (1) Similar expressions stand for the polar representation. Since we already know that polar and canonical harmonic wave functions are linked through the customary expression we will focus on the firstand second-order corrections to the wave functions.
The diagonal elements of the first-order Hamiltonian over the harmonic states are null as well as the off-diagonal ones when two degenerate states are considered. As a matter of fact, we can define the matrix H̅ C as so that eq 44 can be rewritten as By defining |ψ C(0) ⟩ and |ψ C(1) ⟩, the vectors collecting respectively the harmonic states and the corresponding firstorder corrections in the canonical representation, eq 46 can be recast in the matrix form and a similar expression holds in the polar representation The diagonal blocks of H̅ P and H̅ C are null, and the energy difference ε i (0) − ε j (0) reported above is constant within an offdiagonal block. As a matter of fact, we obtain the following identity: By substituting eqs 43 and 49 into eq 48, we obtain that Equation 51 can be rewritten as follows:  (53) Due to the presence of H SR C (2) in eq 53, the terms for which the harmonic states |ψ R C(0) ⟩ and |ψ S C(0) ⟩ are degenerate must be excluded from the summation. From this, an analysis analogous to the one reported for the first-order correction leads to the following identity: (0) ; (1) ; (2) (55) where Θ P(0) , Θ P(1) , and Θ P(2) are given by expressions similar to eq 42 that, recast in the matrix form, have the following expressions in the polar representation: where the subscript N indicates a compact matrix notation for the normalization term. By combining eqs 43, 50 and 54 with eq 56, the following identity is obtained: where ⟨Θ⟩ P and ⟨Θ⟩ C are respectively the representation matrices of the operator Θ̂over the perturbed polar and canonical states. Consequently, the conversion of molecular properties at the anharmonic level is ruled by expressions similar to those used for the transformation of the contact-transformed Hamiltonian. B.2.1.1 Transition Properties from the Ground State. Finally, let us consider the ground state as the initial state. It is worth mentioning that the matrix P only mixes polar and canonical states that are degenerate at the harmonic level so that the perturbed wave function associated to the ground state is independent of the representation ψ ψ | ⟩ = | ⟩ P C 0 0 As a result, the row vector ⟨Θ⟩ 0;S P , containing the values of the transition property (e.g., a component of the dipole moment or polarizability) from the ground state to a set of states sharing the set of principal quantum numbers v S , is related to its canonical counterpart ⟨Θ⟩ 0;S C as follows: 0; (59) or equivalently, where in eq 60, ⟨Θ⟩ 0;S P and ⟨Θ⟩ 0;S C are column vectors.
■ ASSOCIATED CONTENT