Generalized Vibrational Perturbation Theory for Rotovibrational Energies of Linear, Symmetric and Asymmetric Tops: Theory, Approximations, and Automated Approaches to Deal with Medium-to-Large Molecular Systems.

Models going beyond the rigid-rotor and the harmonic oscillator levels are mandatory for providing accurate theoretical predictions for several spectroscopic properties. Different strategies have been devised for this purpose. Among them, the treatment by perturbation theory of the molecular Hamiltonian after its expansion in power series of products of vibrational and rotational operators, also referred to as vibrational perturbation theory (VPT), is particularly appealing for its computational efficiency to treat medium-to-large systems. Moreover, generalized (GVPT) strategies combining the use of perturbative and variational formalisms can be adopted to further improve the accuracy of the results, with the first approach used for weakly coupled terms, and the second one to handle tightly coupled ones. In this context, the GVPT formulation for asymmetric, symmetric, and linear tops is revisited and fully generalized to both minima and first-order saddle points of the molecular potential energy surface. The computational strategies and approximations that can be adopted in dealing with GVPT computations are pointed out, with a particular attention devoted to the treatment of symmetry and degeneracies. A number of tests and applications are discussed, to show the possibilities of the developments, as regards both the variety of treatable systems and eligible methods.

Introduction

Vibrational and rotational spectroscopies are among the most powerful tools for the study of chemical systems.1,2 The investigation of the rotational and rotovibrational spectra of polyatomic molecules has become of basic importance to determine accurate molecular geometries, as well as to get information on molecular force fields, rotovibrational interaction parameters and the relations between structure and chemical-physical properties. Nowadays, there is a constant interplay between molecular spectroscopy and computational chemistry. Indeed, computed data have become crucial for the interpretation of experimental results and, conversely, accurate spectroscopic measurements are used as benchmarks to validate theoretical approaches.1–6

The reliability of the theoretical models to support experimental findings is related to their accuracy. To this end, attention is usually concentrated on the choice of the method used to compute the electronic structure. However, the way in which nuclear motions are simulated is often basic, namely the harmonic approximation for vibrations and the rigid-rotor approximation for rotations. However, the neglect of anharmonicity and rotovibrational couplings can lead to significant errors and may result in incorrect interpretations of experimental data. To overcome such a limitation, various strategies have been devised.7–28 Among them, the approach based on perturbation theory applied to the expansion of the molecular Hamiltonian in power series of products of vibrational and rotational operators, also referred to as vibrational perturbation theory (VPT), is particularly appealing for its computational efficiency to treat medium-to-large semirigid systems.29–43 Moreover, some formulations of VPT, such as the Van Vleck contact transformation method, completely justify a generalized model (GVPT2),44,45 coupling the advantages of the perturbative development to deal with weakly coupled terms and those of the variational treatment to handle tight coupled ones. Implementation of VPT approaches in computational programs for chemistry has become common and black-box procedures have been devised to offer simple yet reliable ways of computing accurate rotovibrational spectra.3,31,46–54

Taking into account that the majority of chemical systems fall into the asymmetric top category and because of the simpler formulation, most developments in the last years have been focused on this case. As a result, a significant ensemble of molecular systems, ranging from small to large sizes, and of interest in various research fields, is excluded or approximately treated. Among others, we can mention organic and organometallic compounds as coronene and ferrocene38,55–57 or acetylene derivatives.58–69

The proper and effective introduction of symmetry leads to different developments for linear, symmetric, and spherical top systems with respect to the formulation of asymmetric tops. Though the rotational problem is simpler in the first three cases than in the last one, because the rigid rotor problem can be solved analytically, the theory of linear, symmetric, or spherical top molecules shows a number of complications due to the presence of degenerate vibrational modes, that makes analytical expressions for the vibrational interaction terms less simple.70,71

The aim of this work is to present a complete framework, able to handle asymmetric tops, as well as, linear and symmetric tops. Starting from the developments already presented in the literature,29,32,47,72–74 we review and generalize the formalism in order to completely support intrinsic and accidental degeneracies, where the first ones are generated by the molecular symmetry and lead to further terms in VPT developments, and the latter are not imposed by the symmetry of the Hamiltonian and lead to singularities in the perturbative formulation, for example, the well-known Fermi resonances.32,75,76 Particular attention is devoted to the latter singularities, presenting their treatment both within the rigorous variational-perturbative coupled GVPT approach, and within approximate methods. Moreover, a fully general formulation of the rotovibrational energies is presented to allow a unified treatment of both minima and first-order saddle points of the molecular potential energy surface (PES). Together with spectroscopic quantities, also thermodynamic functions and reaction rates are considered.

The general formulation can be used in two different ways. On an experimental level, once we have an effective Hamiltonian for a given vibrational state (or for a polyad of such states), we can attempt to determine the values of the spectroscopic constants by fitting them to the experimental frequencies of transitions between the rotation-vibration states.2,21 Such fitting means that we try to obtain the values of the spectroscopic constants that provide the best agreement with the experimental data. On the other hand, we can attempt to evaluate the spectroscopic quantities from a fully quantum mechanical (QM) approach.3,5 To do this, we need a molecular equilibrium geometry together with a set of second, third and semidiagonal fourth energy derivatives with respect to normal modes. The quantities entering VPT expressions can be computed by current electronic structure codes at different levels of sophistication. Hartree-Fock (HF), density functional theory (DFT), and second-order Møller-Plesset theory (MP2) models77–80 will be employed in this article but also other post-HF models (e.g., MCSCF, CCSD(T), etc.) could be used. In this frame, the expressions derived in the first sections can be used to reproduce and/or to predict the experimentally observed results. In the second part of this article, we will validate our implementation showing the feasibility and the limitations of the GVPT approach based on QM electronic computations in reproducing the experimental results.

Theory

Let us start by reminding that a symmetric top is defined by two properties; the equilibrium configuration of the nuclei has a symmetry axis of order 3 or higher and, if there is more than one axis satisfying the above condition, these axes are all coincident. If all the above conditions are present, the molecule has two equal moments of inertia. Otherwise, the molecule is either an asymmetric top (first condition not met, all moments of inertia are different) or a spherical top (second condition not satisfied, all moments of inertia are equal). Moreover, in a linear-top system all nuclei are aligned and the molecule has one vanishing moment of inertia and two non-null coincident ones.

Asymmetric tops have only nondegenerate harmonic vibrational frequencies, whereas linear and symmetric tops have both nondegenerate and doubly degenerate harmonic frequencies, and spherical tops can be affected by degenerations larger than two. The development presented in the following considers systems having at most doubly degenerate harmonic frequencies, letting aside the case of spherical tops.

As the general development of the theory relies on a significant number of equations, in order to make our presentation easier to follow, we have chosen to shift redundant formulas or the most cumbersome equations to specific appendices.

Molecular Hamiltonian and perturbation theory

Within the Born-Oppenheimer approximation,81,82 where the total Hamiltonian of a molecule can be separated into an electronic and a nuclear component, the Eckart-Sayvetz conditions are applied to minimize the coupling between the rotational and vibrational wavefunctions.76,83,84 The rotovibrational QM Hamiltonian for the nuclei in a given electronic state can be written,32,76,85

where is an element of the effective inverse molecular inertia tensor and and are, respectively, the components of the total and vibrational angular momentum operators along the molecule-fixed Cartesian axes τ or η.30,32,76,86 The explicit form of the latter is,

where is the matrix of the Coriolis coupling constants. Q and P are the mass-weighted vibrational normal coordinate and its conjugate momentum associated to the vibrational mode i, respectively, and the summations run on normal coordinates ( for linear systems). is the PES in which nuclei move and is a mass-dependent contribution, which vanishes for linear systems,76,86

In eq. (1) both and can be expanded as Taylor series of the mass-weighted normal coordinates about the equilibrium geometry,32,76

where is an element of the inverse of the equilibrium inertia moment of the molecule and.3,32,86, where ω is the classical frequency of vibrations, and and are respectively the third and fourth derivatives of the potential energy with respect to the normal modes, also referred to as the cubic and quartic force constants,32,33,76

After substitution of and in eq. (1) by their respective definitions in eqs. (2) and (5), the terms in can be written as,

where represents all the terms with a degree f in the vibrational operators (Q or P) and degree g in the rotational operators (). Hence, collect purely vibrational terms,

where,

are the terms of the expanded Hamiltonian corresponding to the zeroth-order development of written in term of the equilibrium molecular rotation constant. Note that all the constants in eqs. (8–10) are given by slightly non standard expressions based on mass-weighted vibrational normal coordinates, rather than on their reduced counterparts, since this allows a cleaner treatment when dealing with transition states (TS), rather than energy minima, avoiding complex force constants.36,40,87–89

and collect the Coriolis and rotovibrational terms, respectively. More complete expressions have been reported by Aliev and Watson (see Table 1 in Ref.86). Here, we reproduce only the lower-order terms,

Table 1

Non-zero off-diagonal variational elements involved in the first order vibrational (Fermi) resonances.

Type I Fermi resonances

Type II Fermi resonances

This way, can be treated perturbatively, taking as zeroth-order contribution the harmonic oscillator Hamiltonian,. The separation in perturbative orders of terms has been widely discussed in the literature, and different classification schemes have been proposed.3,6,29,30,32,86 A detailed assignment was proposed by Aliev and Watson (see Table 2 of Ref.86). It is noteworthy that the rigid-rotor term,, is usually treated as part of the perturbation to avoid rotational energy differences in the denominators of the perturbation development.

Table 2

Δ and terms involved in the DSPT2 treatment of diagonal elements

Type I Fermi resonances	Type II Fermi resonances

and the slash symbol (“”) between latin numbers is used as a separator between the possible force constants for which the relation stands.

Various formulations of perturbation theory have been devised, such as the Rayleigh-Schrödinger method90,91 (RS), the Bloch projector formalism,92,93 or the Van Vleck contact transformation approach (CT).44,45 We recall here the main features of the CT method. Differences with the RS development, also commonly used in the literature, will be highlighted. The CT formalism is based on the transformation of the Schrödinger equation,3,30,32,86,94

where the original Hamiltonian and wavefunction are transformed as,

is an Hermitian operator so that is unitary. It is chosen to obtain an effective block-diagonal Hamiltonian in a given basis, in order to separate each vibrational level or block of degenerate or near-degenerate vibrational levels, with the property that the eigenvalues of these blocks are the same as for. The operator can be written as a product of successive contact transformations,

where is chosen in order to diagonalize up to the n-th order. Up to the second-order, eq. (16) for corresponds to,

where represents a commutator. Taking matrix elements in the basis of eigenfunctions of, let us first consider the terms to illustrate the choice of,

where the uppercase subscript represents states with different energies and the lowercase one differentiates degenerate states. This means that is the eigenvalue for all eigenstates of the zeroth-order Hamiltonian. For the case, which is also referred to as a diagonal matrix element of, the second term in the right-hand side of eq. (21) vanishes, that is,

which is identical to the result derived via RS first-order perturbation theory.90,91 For the off-diagonal elements with, the first-order interaction term will vanish if we choose satisfying the following equation,

In this case, will only contribute to the effective Hamiltonian for perturbation orders higher than the first one. If, the value of as defined in eq. (23) will be excessively large. In this case, and are said to be in resonance and is set to be null, so that,

The case of degenerate states, where, is treated in the same way as for states of near-equal energies, with the term set to be null, so we have,

It is noteworthy that this off-diagonal term can result in the lifting, also called doubling, of the zeroth-order energy degeneracy.

The same considerations apply for the choice of in eq. (20), with the difference that, now, we impose that the terms vanish and is the perturbation correction to that derives from the cancellation of the off-diagonal terms of. It can be shown that the general matrix element of is given by the expression,94

where the first summation, with the symbol, is only carried out over the nonresonant states. It is noteworthy that for the elements, be it b = a and, the above equation reduces to,

which is identical to the matrix element derived via RS second-order perturbation theory.90,91 Conversely, the derivation of the off-diagonal elements of with from the Rayleigh-Schrödinger development is less rigorous. For this reason, an alternative form with respect to eq. (26) has been often used for the treatment of the latter,95–97

where.

Vibrational energies for asymmetric, symmetric, and linear tops

A pure vibrational Hamiltonian is obtained by correcting with and, followed by the transformation step described before.32,86 An additional term is usually included to account for the zeroth-order expansion of [(see eqs. (3) and (4)],30,32,76,86

where Γ = 1 for asymmetric and symmetric top systems, and Γ = 0 for linear systems. It should be noted that, due to its small contribution, this term is generally neglected.

If no resonance occurs, the first-order effect of does not contribute to the energy of any vibrational state, since both diagonal [eq. (22)], and off-diagonal [eq. (25)], terms vanish. Hence, the perturbative corrections to the energy up to the second order are all due to, with the largest contribution related to the diagonal elements. Nielsen first derived the solution for the latter,29 which was subsequently refined with more general formulas.29,30 Later, Plíva fixed omissions for symmetric tops with a principal axis of order higher than three,72 mainly due to missing force constants. His formulas were in turn corrected by Willetts and Handy.73 Following those works, we present here a new derivation, taking advantage of the framework built previously for asymmetric tops,47,53 done with an ad hoc tool, based on a symbolic algebra program.98

By applying specific rules to orient the degenerate normal modes,70,71 simple symmetry relations can be established between sets of related cubic and quartic force constants, as well as Coriolis constants. A first detailed classification was done by Henry and Amat in Refs.60, 99, for the first, and Refs.70–71 for the latter. For the force constants, at variance with eqs. (9) and (10), restricted sums were used in the potential energy expansions. Remembering that the commutator of the two normal coordinates associated to the same harmonic frequency is null, the nonvanishing cubic and quartic force constants with at least one degenerate normal mode for the case of unrestricted summations have been reordered and reported in TablesA1–A9 of Appendix A. The notation adopted in those Tables is similar to the one used by Plíva.72 Moreover, assuming hereafter the highest-order axis of symmetry to be along the z axis in the molecule-fixed reference frame, the symmetry relations affecting the Coriolis terms are given in Appendix A.

Table A1

Symmetry groups labels. I and II are non-abelian and abelian, respectively

I	Ia:	CNv, DN, DNh (any N); DNd (N odd);
	Ib:	( even);
II	IIa:	CN, CNh (any N); (N odd);
	IIb:	SN ( even).

Table A9

Non-vanishing quartic energy derivatives K with respect to and their symmetry relations.60

Symmetry	w.r.t.		Group	N
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			II
			II
			II
			II
			II
			II
			II
			II
			II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

From here on, the subscripts will be used to indicate generic vibrational modes, degenerate or not, while will be reserved to nondegenerate modes and to degenerate ones. When needed, a second subscript γ,, ι, which takes the values 1 or 2, is used to distinguish the two different normal coordinates associated to the same two-fold degenerate harmonic frequencies. For TSs, the transition vector (i.e., the normal mode with the nondegenerate imaginary frequency) is labeled by the subscript F. In this framework, the vibrational second-order perturbation theory leads to the following expression for the energies,

with,

δ is the Kronecker's delta, and are respectively the principal and angular vibrational quantum numbers, and d is the degeneracy of mode i. In the above expression, all- and-independent terms are collected in E0, a term which can be written in a form devoid of resonances,

with,

and (see Appendix A),

The elements of the anharmonic matrices and are given by,

with and (see Appendix A),

In the formulation adopted here, it is easy to see from eqs. (36) to (42) that the matrix elements χ, with, are imaginary. They are excluded from the vibrational energy, which contains only real terms, and enter, together with the imaginary frequency ω, in the expression providing tunneling and non classical reflection contributions to reaction rates.53

It is noteworthy that, at variance with eq. (30), the anharmonic contribution to the vibrational energy is usually expressed in the literature as the sum of and χ0 (or G0) terms. In the specific case of symmetric and linear tops, the χ0 term was omitted by Plíva, Willetts and Handy in their respective works.72,73 It was included in the derivation proposed by Truhlar and coworkers39 but it was based on a less general treatment than the one proposed by Plíva, which led to discrepancies with respect to the formulas given by Willetts and Handy and obtained in the present work. To the best of our knowledge, this is the first time that all terms needed to compute the vibrational energy as given in eq. (30) for symmetric, asymmetric and linear tops are gathered in a single work.

From eq. (30), it is possible to calculate the energy of any vibrational state. The energy of the vibrational ground state, that is the zero-point vibrational energy (ZPVE), is. It is straightforward to determine transition energies governing vibrational spectra (i.e., at constant n) with the relation,

Explicit expressions for the energies of fundamentals, first overtones and combination bands are given in the Appendix B.

Finally, the tunnel probability P, of interest in chemical rate constants computations, can be evaluated using the microcanonical ensemble with the semiclassical TS theory of Miller and coworkers.100,101 They used the definitions,

to invert the relation, where,

and obtain the generalized barrier penetration integral in terms of the n and l quantum numbers of the activated system, with, and the total energy E,

where,

In this framework, the semiclassical tunneling probability P for a one-dimensional barrier is given by,

Vibrational l-type doubling and l-type resonance

If no resonances occur, vibrational energies of nondegenerate states can be determined directly from eq. (30). On the other hand, for degenerate zeroth-order states, as seen above, the interaction terms cannot be canceled out with and must be treated variationally. The presence of those off-diagonal elements in the variational matrix will result in a further lifting of the degeneracy of the vibrational energies, initiated with the application of the second-order correction. This splitting is called l-type doubling or l-type resonance, depending if the diagonal energies involved have equal or different values, respectively. Using symmetry considerations, Amat derived a general rule to identify a priori the possible non-null off-diagonal matrix elements.32,102 It depends on the N-fold principal symmetry axis and the difference of quanta in the principal () and angular () vibrational quantum numbers between the states involved in the interaction term. The ensemble of non-zero l-type off-diagonal terms is obtained from the following relations,

where, as usual, only the modes undergoing a change in their quantum numbers between the two states involved in the matrix elements are shown. The off-diagonal elements given in eq. (53) are non-null if N is a multiple of 4, those given in eq. (54) for any symmetric top molecule and the elements of eq. (55) if N is even.

The first expressions of U, R and S for the various point groups have been given by Grenier-Besson.103,104 The formulas have been re-derived here, with the notation introduced in this work, and validated with respect to those obtained by Grenier-Besson. They are gathered in the Appendix C.

Vibrational first-order resonances

It has been shown that if two states are in resonance it is not possible to make the corresponding off-diagonal term vanish. A resonance can connect two or several vibrational levels and, moreover, multiple resonances can connect a network of levels. The sub-matrices where the resonances are involved are called polyads.86,97

As has only diagonal elements, its off-diagonal terms are all null. The presence of off-diagonal first-order terms due to is related to the so-called Fermi resonances. The latter are characterized by a strong interaction between two states that differ by one quantum in one mode and two quanta in either one (type I) or two different (type II) modes.32,33,75 Due to the creation of one vibrational quantum and the annihilation of two others, or conversely, these singularities are also called vibrational 1-2 resonances.97 They can appear when in eq. (23) is excessively large or in eq. (21), condition which can occur in two cases: (type I) or (type II).

Different methods have been developed to overcome the problem of Fermi resonances. One possible route is to solve the Dyson equation with the frequency-dependent self-energy.54 In this way, one need not to classify the different types of resonances or lose size-consistency, but to perform a root search of a nonlinear, recursive equation. The most common approach, called deperturbed VPT2 (DVPT2), consists in simply removing from the perturbative treatment the resonant terms after their identification. The explicit expressions of the potentially resonant terms in eqs. (36–42) are given in Appendix D. However, this treatment is incomplete due to the neglect of the resonant terms. An improvement can be obtained by treating variationally the levels involved in the resonance, reintroducing the removed terms as off-diagonal interaction elements. This method has been called generalized VPT2 (GVPT2)29,30,32,47 or, more recently, CVPT2+K94 or CVPT2+WK.105 The list of possible off-diagonal first-order interaction terms generalized to linear, symmetric and asymmetric tops is given in Table 1.

Although those methods have been widely discussed in the literature, less attention has been devoted to the identification of a general strategy to determine when an interaction term has to be considered in resonance. Indeed, all the methods presented above rely directly on the identification of the resonant terms. The definition of a singularity giving rise to unphysical contributions is far from straightforward, and different schemes have been proposed. The simplest approach is to check the magnitude of the denominator (i.e., and) with respect to a fixed threshold. If the value is below this limit, the term is considered resonant. Such a scheme does not account for the magnitude of the numerator, which makes difficult the definition of a reliable threshold adapted to a wide range of molecular systems. A more robust solution to this problem has been suggested by Martin and coworkers.106 Considering two resonant states and, we can write down the interaction between the two states as a variational matrix,

where and is the complex conjugate of ρ. If ρ tends to zero, the eigenvalues of the matrix in eq. (56) can be written as the following Taylor series,

where must be non-null. Up to the second-order, coincides with the vibrational energies or corrected with a second-order perturbation term, which arises from the interaction between and (here the case),107

where is precisely the possible resonant term in the VPT2 equations, that is, one of the terms in the summation in the right-hand side of eq. (27). Based on those considerations, the importance of the higher-order perturbative terms can be estimated from the fourth-order expansion term in eq. (57),

where for type I Fermi resonances and for type II Fermi resonances. Consequently, a threshold on the term can be a good marker to evaluate the importance of higher order effects and then if the second-order term has to be treated as resonant. Moreover, this term accounts not only for the energy difference but also for the magnitude of ρ. In a slightly different formulation, the threshold used to evaluate the presence of first-order resonances is calculated taking into account all high-order expansion terms, obtained subtracting the first two expansion terms from the square root of eq. (57),105

A general approach can be derived from the development presented above, which is to apply to all potentially resonant terms in the VPT2 formulas the transformation described previously,

An interesting feature of this approach is that there is no need for an identification of the resonant terms, which can be inconsistent whenever one has to consider a series of force fields for a given system, or a series of geometries along a reaction path. Indeed, variations in the set of resonant terms can make difficult any comparison of the VPT2 results between two or more simulations. This scheme is similar to the second-order degeneracy-corrected perturbation theory (DCPT2) introduced by Kuhler and coworkers,108 which will be discussed afterwards. The interest is to prevent the appearance of singularities in the calculation of anharmonic contributions using a simplified variational approach, since the right-hand side of eq. (61) cannot diverge if Δ becomes small. Far from resonance, the substitution still accounts for the interaction between the vibrational states and. It is noteworthy that, at variance with what has been done in Refs.53 and108, this time we apply the transformation of eq. (61) directly on all possibly resonant terms in the effective Hamiltonian, that is all terms in the summation in the right-hand side of eq. (27) which have frequencies differences (i.e., or) in the denominator. For this reason, we will refer to this approach as degeneracy-smeared vibrational perturbation theory (DSPT2). After the complete development of eq. (27), the possibly resonant terms can be grouped in sets of 2 or 4 components sharing the same Δ. For the two terms with the same Δ the substitution given in eq. (61) leads to,

with and. Since and are opposite, the last term of the transformation disappears.

As an example, let us consider the terms involving,

The substitution given in eq. (62) can be carried out with the following definitions,

The transformation to be applied in the case of 4 terms having the same Δ is straightforwardly derived,

with. As before, the previous transformation can be further simplified as the term is null.

All potentially resonant terms and the definition required to apply the transformation given above are gathered in Table 2. The extension of the DSPT2 treatment to the off-diagonal elements requires further discussion. Let us consider one of the terms in the summation in the right-hand side of eq. (27) with. This contribution can be related to the eigenvalues of the following matrix,

where and, with associated eigenvalues,

and. This matrix differs slightly from the one obtained with the proper variational description, which has the form,

Nevertheless, the matrix given in eq. (65) is more convenient for the mathematical derivation of the possible resonant terms, on which the previous substitution is applied,

where accounts for the signs of both and.

To illustrate this point, let us consider the resonant term with in,

We then apply the relation given in eq. (68) after the proper identification of the terms involved in the transformation,

The other identification sets to be used in the transformations of the possibly resonant terms in U, R and S are gathered in Table 3. An alternative way to treat resonances was proposed by Kuhler and coworkers in 1995 and slightly modified by some of us. The difference with the DSPT2 development lies in the terms on which the substitution given in eq. (61) is applied. Indeed, in degeneracy-corrected PT2 (DCPT2), the elements of the matrix are derived first and the possibly resonant terms are identified within the elements of χ [eqs. 36–40] and transformed. Further details can be found in Refs.53 and108. For degenerate modes, not treated in those previous works, we use the same transformation as for nondegenerate modes. To illustrate this point, let us consider the last term in the right-hand side of eq. (39), developed in partial fractions,

Table 3

Δ, s, and terms involved in the DSPT2 treatment of l-doubling off-diagonal elements

U l-type doubling

R l-type doubling

S l-type doubling

and the slash symbol (“”) between latin numbers is used as a separator between the possible force constants for which the relation stands.

By setting and, we obtain the following transformation,

The new matrix obtained by replacing possibly resonant terms in nonresonant ones is then used in the calculations of the vibrational energies.

However, both DSPT2 and DCPT2 transformations can give poor results far from resonance when both numerator and denominator become large. Indeed, when ρ is large, the equivalence of eq. (57) is not true and, while the VPT2 term can be still valid due to a large Δ, the DSPT2 and DCPT2 transformations are incorrect. To cope with this shortcoming, an hybrid scheme called hybrid DCPT2-VPT2 (HDCPT2) has been proposed by some of us. In this method, a switch function, Λ, is used to mix the results from the original VPT2 and the DCPT2 approaches for all possibly resonant terms in as follows,53

where represents the value of a possibly resonant term calculated with the original VPT2 formulation [left-hand side term in eq. (71)], and its counterpart calculated by mean of DCPT2 [right-hand side term in eq. (71)]. Λ is defined as,

where β controls the transition threshold between DCPT2 and VPT2, and α the “smoothness” of the transition. The same scheme applies for the hybrid DSPT2-VPT2 (HDSPT2),

where is the true VPT2 term [e.g., in eqs. (62) or (64)] and is its DSPT2 counterpart (i.e.,).

Vibrational second-order resonances

In analogy with first-order resonances, when two zeroth-order states involved in the contact transformation given by are close to each other, the off-diagonal elements cannot be canceled out and have to be treated variationally. Many types of resonances lead to off-diagonal second-order energy corrections. According to the classification of the total change of quanta, there are 1-1, 2-2 and 1-3 second-order resonances. For asymmetric tops, a detailed description of all these off-diagonal terms has been recently given by Rosnik and Polik.94 The total number of non-zero second-order off-diagonal elements becomes very large when doubly degenerate normal modes are also taken into account, because of the large number of combinations of nondegenerate/doubly degenerate normal modes that can be obtained when considering all states involved in the matrix elements. In this work, we have generalized the expression for the 2-2 vibrational second-order resonances to support also doubly degenerate states, in the specific case of the annihilation of two quanta in one mode and the creation of two quanta in another one. Known also as Darling-Dennison resonances,109 the non-zero off-diagonal elements for this situation are given by,

The definition of the κ terms is reported in Appendix E. The second-order off-diagonal elements are then used within the GVPT2 approach in the variational treatment of the polyads.

Therefore, each polyad contains the deperturbed vibrational energies of the resonances interacting states as diagonal elements, the first- and second-order resonances off-diagonal elements, as well as the possibly l-doublings and l-resonances, also off-diagonal terms. Note that, up to the second order, we will never have, with both and non-null, because the couples of states interacting within first-order resonances are always different from the couples interacting by second-order resonances.

Vibrational partition function for thermodinamics and kinetics

The partition function of a system is the sum of the Boltzmann factors of the energy levels each weighted by its degeneracy,110

where, k and T are the Boltzmann constant and the temperature, respectively, and the summation is on all possibly states σ. We treat here the vibrational molecular partition function, for which and are the energies and degeneracies of vibrational levels. Starting from eq. (80) and focusing on at most doubly degenerate vibrational modes, the harmonic vibrational partition function is obtained by,

where is the harmonic formulation of the vibrational energy and is the degeneracy due to the degenerate mode s. Developing the previous expression,

where is the harmonic ZPVE and we have used the relations and when.

Unfortunately, an analytical development of is not available beyond the harmonic level. Several routes have been proposed to deal with this situation.111–113 Here, we employ the approximated method proposed by Truhlar and Isaacson, called simple perturbation theory (SPT), in which the formal expression of the harmonic partition function is retained, but the ZPVE and ω terms are replaced with their anharmonic counterparts,35,108,114

E0 is the anharmonic ZPVE given in eq. (32), and, defined in eq. (B1), is reduced to ν below for the sake of readability. This approximation leads to analytical expressions for the vibrational contributions to the internal energy U, entropy S, and constant volume specific heat c,91,110

where R is the Boltzmann universal gas constant.

Rotational Hamiltonian

The terms ( in the effective Hamiltonian are the pure rotational and centrifugal contributions to the energy, which describe the rotational energy levels for the zero-point vibrational state. Their complete treatment has been widely discussed in the literature3,6,32,115–120 and we will recall here only some key aspects.

The quartic centrifugal terms form the simplest second-order contribution to. Their expression results from the second-order effect of,

where the tensor was originally introduced by Wilson.121,122 The sextic centrifugal distortion constants are from the term. The perturbation terms required for its calculation are (harmonic), (anharmonic) and (Coriolis), where the last two Coriolis contributions should be considered even if they have a degree in J greater than six because they can be reduced to sixth degree terms by the rotational commutators (i.e., (see Refs.32, 115 for further details). With the assignment and, all contributions reported above appear in the fourth order perturbative development. The final expression for was obtained by Chung and Parker123,124 and collected by Aliev and Watson (see Table 3 of Ref.120).

The vibrational contact transformation then leads to the rotational Hamiltonian,

where now both and contain terms that can be reduced by the use of rotational commutation relations. Taking as an example the explicit form of given in eq. (87), there are terms that can be reduced to,

where

As a consequence of this reduction, is corrected by a small contribution from the quartic terms,

and are obtained by cyclic permutation of the indices.

Further contact transformations with purely rotational operators, thus diagonal in the vibrational quantum numbers, are required in order to achieve a complete reduction of. In the completely reduced Hamiltonian, combinations of quartic and sextic distortion parameters are strictly related to the eigenvalues of, and then to physical observables. Different results can be obtained depending on the arbitrary choice applied to fix the reduction's parameters. The general form of the reduced Hamiltonian of an arbitrary molecule has been given by Watson.118,119,125 With the choice called by Watson asymmetric top (A) reduction, the matrix representation of in the symmetric top basis has the same form as that of a rigid asymmetric top,

where and are the total angular momentum and the ladder operators, respectively126 and represents an anticommutator. Δ and δ refer to the quartic distortion constants,, and to the sextic ones. The latter coefficients are given in Refs.86, 115. The disadvantage of the asymmetric top reduction is that it fails for both genuine and accidental symmetric tops. For the latter, the symmetric top (S) reduction suggested by Winnewisser and Van Eijck can be used,127,128

where the expression for the quartic (D and d) and sextic (H and h) distortion constants are presented in Ref.115. General expressions for sextic distortion constants have been recently revised in Ref.129.

For linear molecules, the angular momentum J is null. In this case, Watson has shown that the molecular Hamiltonian in eq. (1) becomes,130

for linear molecules is then given by,30,76,86,115

in which B is the equilibrium rotational constant and the explicit formulation of the quartic (D) and sextic (H) centrifugal distortion constants are given in Refs.32, 115. is already in a fully reduced form. The rotational energies for linear tops are obtained by replacing with and then by their eigenvalues,

where J is the total angular momentum quantum number and l the total vibrational angular momentum.

Vibrational dependence of the rotational Hamiltonian

The operators contain the terms describing the dependence of the rotational and centrifugal constants on the vibrational quantum numbers. The vibrational dependence of the rotational constants in the quartic approximation is described by,

where now indicates a specific vibrational state. The vibrational correction derives from the diagonal matrix elements of, specifically by the second-order corrections, considering and. For asymmetric tops, the constants are given by,30,32,86

Using the symmetry relations for and given in Refs.70–71 and accounting for the doubly degenerate normal modes, the α coefficients for linear and symmetric tops are,32

with. The first contribution in eqs. (98–102) is a corrective term related to the moment of inertia, the second one is due to the Coriolis interactions, and the last is an anharmonic correction. It is noteworthy that the Coriolis coupling term may be affected by resonances. In analogy with vibrational first-order resonance, the strategy that is adopted when a resonance occurs is to expand the Coriolis term and neglect the resonant part, as shown in Appendix D. By contrast, the summed Coriolis coupling term is not affected by resonances, as it is possible to write,

Taking as an example the resonance, we have (d = 1 and d = 1),

Similar simplifications can be applied for [note that the factor 1/2, which multiplies the Coriolis terms in eq. (102), is simplified by d = 2] and resonances. Taking these considerations into account, it easy to see that eq. (97) for the vibrational ground state is devoid of resonances, that is,.

Computational Details

The theoretical approach presented in the previous section has been included in a development version of the Gaussian package.131 The implementation can be used with any QM procedure for which analytical second derivatives are available, among which HF,77 DFT,78 and MP279 will be explicitly considered in the following. Examples of applications with each model will be given in the next section. Within DFT, the standard B3LYP functional132–134 has been used in conjunction with the SNSD basis set,135 that has been validated for vibrational studies.136–139 The double-hybrid functional B2PLYP140 and MP2 have been used in conjunction with the Dunning correlation-consistent valence aug-cc-pVTZ (AVTZ) and aug-cc-pVQZ (AVQZ) basis sets.141,142 For ferrocene, an organometallic compound taken as an example of medium-size systems, the B3LYP functional has been used in conjunction with the SNSD basis set for H and C atoms and the double-ζ ECP basis set of Hay and Wadt augmented with polarization functions (p type with exponent) (aug-LANL2DZ) for Fe, with the LANL2DZ pseudo potential to describe core electrons.143 The hybrid B3PW91 functional133 has been also employed in conjunction with the m6-31G basis set, based on 6-31G and improved for first-row transition metals.144 For triphenylamine, the B3LYP functional has been coupled with the valence double-ζ polarized basis set 6-31G*.145–148 Frequency calculations have been systematically carried out at the equilibrium geometry obtained at the same level of theory, using respectively tight () and very-tight (on force: Hartree/Bohr, estimated displacement: Bohr) convergence criteria for the self-consistent field and geometry optimization steps, respectively. For all DFT computations, an ultra-fine grid (199 radial points, 590 angular points) was used for the numerical integration of the two-electron integrals and their derivatives. The third and semidiagonal fourth derivatives of the PES have been obtained by numerical differentiation of the analytical second derivatives along the mass-weighted normal coordinates, with the default step, as,47,149

It should be noted that the calculation of the cubic and quartic force constants is the most demanding step in terms of computational cost. It can be sped up by using a reduced-dimensionality scheme where the numerical differentiations are done along a subset of normal coordinates corresponding to the modes to be treated anharmonically. In this case, the averaging done for and is applied over the number of elements actually calculated (1, 2 or 3 for and 1 or 2 for). Note that, if finite differentiation is performed along mode i, but not along modes j and k, the force constants and can not be evaluated. The anharmonic corrections for fundamental and combination bands of will still be given by eq. (B1) and eqs. (B2) and (B3), respectively, where χ and g terms are unchanged, whereas χ terms differ from the fully-dimensionality ones for the absence of the elements [see eqs. (36–42)],

More details on those schemes are available in Refs.150–151, while an example of application will be given in the next Section.

A hybrid CCSD(T)/DFT approach has also been used to carry out VPT2 calculations,137,152–154 where the harmonic frequencies are evaluated at the CCSD(T) level and the anharmonic correction at the DFT level. This scheme is based on the observation that most of the discrepancy with experimental results is due to the harmonic frequencies, which can be corrected by employing a higher level of theory. The CCSD(T) harmonic frequencies are inserted in eq. (30) in place of the DFT ones. In order to get reliable results, the equilibrium geometries and the normal coordinates at the CCSD(T) and DFT levels must be consistent. This is automatically checked by our procedure when applying the hybrid scheme.

To overcome the problem of 1-2 resonances in VPT2 calculations, the computational strategies presented in the previous section have been employed. For the DVPT2 and GVPT2 approaches, a term is identified as resonant if the absolute frequency difference in the denominator, Δ, is smaller than 200 cm−1 and Ξ in eq. (59) is larger than 1 cm−1. The default parameters previously used for HDCPT2 ( with ρ and Δ in cm−1) have been used to compute Λ for both HDCPT2 and HDSPT2, see Ref.53. Vibrational second-order 2-2 resonances are identified by two criteria: the absolute frequency difference between the two resonant states must be smaller than 10 cm−1, and the off-diagonal term greater than 20 cm−1. For Coriolis resonances, the terms in eqs. (95–98) with an absolute frequency difference lower than 20 cm−1 are discarded.

Results and Discussion

Full DFT and hybrid methods for the vibrational energies of small- to medium-sized linear systems

A set of linear molecules, that is, HCN, HNC, OCS, HCP, CO2, C2H2 and C4H2, have been selected to test the performance of full DFT and hybrid CCSD(T)/DFT methods to calculate the anharmonic corrections to the vibrational frequencies. On these molecules, all the schemes presented in the previous section to treat first-order resonances have been employed, and the results for the l-doubling interaction terms have been directly compared with the experimental data when present in the literature.

The VPT2 anharmonic corrections for the linear systems HCN, HNC, OCS, and HCP, shown in Table 4, were calculated at the MP2, B3LYP and B2PLYP levels of theory, in conjunction with AVTZ and AVQZ, as well as SNSD for B3LYP, basis sets. In the Table, the best theoretical results, computed at the CCSD(T) level, and experimental data are also reported for comparison purposes. For those systems, which are not affected by resonances, the anharmonic corrections calculated with the different methods are very close to one another. The main discrepancies with experimental results are found to be related to the harmonic part. More precisely, the corrections to the nondegenerate frequencies are very close to the observed values, while the corrections to the low-degenerate wavenumber show a greater sensitivity to the electronic methods and the size of the basis set. For HCN, OCS, and HCP, B3LYP/SNSD gives very good result, while, for HNC, the large anharmonic correction for the degenerate wavenumber is due to its underestimation of the quartic force constants.

Table 4

Comparison of computed and experimental harmonic ω and anharmonic fundamental VPT2 wavenumbers ν for the linear molecules HCN, HNC, OCS, HCP (in cm−1)

		MP2		B3LYP			B2PLYP		CCSD(T)	Expt.
		AVTZ	AVQZ	SNSD	AVTZ	AVQZ	AVTZ	AVQZ
HCN[a]
ω1	∏	718	721	747	759	758	745	745	729	727
ω2	Σ	2022	2034	2196	2200	2201	2125	2129	2125	2129
ω3		3467	3466	3449	3444	3440	3460	3456	3435	3442
		715	718	729	745	744	733	733	717	714[e]
		1987	1999	2169	2173	2175	2094	2098	2096	2097[e]
		3334	3339	3317	3312	3312	3327	3328	3309	3312[e]
		−3	−3	−18	−14	−13	−12	−12	−12	−13
		−35	−35	−27	−26	−26	−31	−30	−29	−32
		−133	−127	−132	−132	−128	−133	−128	−126	−130
HNC[b]
ω1	∏	485	488	477	468	467	467	467	471	490
ω2	Σ	2016	2027	2097	2103	2104	2059	2063	2044	2067
ω3		3818	3824	3801	3799	3801	3815	3818	3837	3842
		505	497	355	463	463	469	470	474	477
		1983	1993	2063	2069	2070	2023	2027	2008	2029
		3656	3661	3631	3634	3635	3650	3652	3666	3653
		+20	+9	−122	−5	−4	+2	+3	+3	−13
		−33	−34	−34	−34	−34	−36	−36	−36	−36
		−162	−163	−170	−165	−165	−165	−165	−171	−189
OCS[c]
ω1	∏	506	524	518	527	527	523	523	524	524
ω2	Σ	888	893	865	874	876	872	875	872	876
ω3		2124	2092	2116	2108	2110	2079	2083	2095	2093
		502	520	514	523	524	519	520	520	521
		869	876	849	858	860	855	859	855	863
		2097	2064	2084	2078	2080	2048	2052	2064	2060
		−4	−4	−4	−4	−3	−4	−3	−4	−3
		−19	−17	−16	−16	−16	−16	−16	−17	−13
		−27	−28	−32	−31	−30	−31	−31	−31	−33
HCP[d]
ω1	∏	677	689	697	712	720	699	707	689	688
ω2	Σ	1245	1255	1322	1338	1342	1291	1297	1299	1298
ω3		3355	3360	3345	3349	3348	3359	3359	3345	3346
		678	680	682	700	704	689	693	675	675
		1226	1236	1304	1319	1323	1272	1278	1281	1278
		3231	3233	3216	3219	3219	3231	3231	3213	3217
		+1	−9	−15	−13	−16	−9	−14	−14	−13
		−19	−19	−18	−19	−18	−19	−19	−18	−20
		−124	−128	−129	−130	−129	−128	−129	−132	−129

Δ represents the anharmonic correction.

Reference values were taken from:

CCSD(T)/AVTZ and experimental values from Ref.155.

CCSD(T)/ANO1 and experimental values from Ref.156.

CCSD(T)/CVQZ and experimental values from Ref.157.

CCSD(T)/CV5Z and experimental values from Ref.158.

experimental values from Ref.68.

CO2 represents an interesting test to validate the DCPT2 and DSPT2 schemes in presence of resonances. It has been one of the first molecules used in infrared and Raman measurements and has served as a prototype for the study of resonances. Vibrational wavenumbers for fundamental, overtones and combination bands obtained at the B2PLYP/AVQZ level and with the hybrid scheme, where the CCSD(T)-F12a/AVTZ harmonic frequencies taken from Ref.155 are used in conjunction with the B2PLYP/AVQZ force field, are shown in Table 5. The states are grouped based on the polyads. The well-known type I Fermi resonance that affects this system is due to, with normal modes 1 and 2 of (and (symmetry, respectively. The lowest energy states that are affected are collected in the following four polyads: with with and with, and with. Note that the states are not involved in the latter polyad since their interaction with is symmetry forbidden. From a numerical point of view, this is due to the fact that only is non-null for linear systems (see Tables 1 and 17). The discrepancies of the GVPT2 frequencies at the B2PLYP/AVQZ level with respect to the experimental results are mostly due to the underestimation of the ω2 harmonic frequency (1344 cm−1 vs. 1351 cm−1), as confirmed by the improvements obtained with the GVPT2 hybrid scheme, which leads to satisfactory agreements (the discrepancies never exceed 5 cm−1 and are on average 1-2 cm−1).

Table 5

Comparison of experimental and computed harmonic ω and anharmonic ν wavenumbers for CO2 (in cm−1)

	B2PLYP[a]				HYBRID[b]				Best theo.		Expt.
State	ω				ω				ω	ν	ν
	669	664	642	664	673	668	646	668	673[c]	668[c]	668[c,d,e,f]
	1344	1275	1197	1285	1351	1284	1202	1293	1351[c]	1285[c,g]	1285[c,d,e,f,g]
	1337	1382	1374	1374	1346	1390	1381	1381		1388[g]	1388[d,e,f,g]
	1337	1330	1286	1330	1346	1338	1293	1338		1336[d]	1336[d,e,f]
	2006	1918	1752	1934	2018	1931	1759	1947		1933[d]	1934[d]
	2013	2070	2061	2055	2024	2082	2072	2066		2077[d]	2077[d,f]
	2390	2345	2345	2345	2391	2347	2347	2347	2391[c]	2349[c,g]	2349[c,d,e,f,g]
	2688	2526	2220	2581	2702	2543	2227	2597		2548[g]	2548[g]
	2674	2656	2660	2679	2691	2671	2673	2694		2671[g]	2671[g]
	2681	2791	2742	2714	2696	2797	2757	2729		2797[g]	2797[g]
	4780	4666	4666	4666	4782	4670	4670	4670		4673[g]	4673[g,f]
	3059	2997	2975	2997	3064	3003	2980	3003			3004[g,f]
	3734	3600	3517	3605	3742	3610	3524	3615		[d,e,g]	3613[d,e,f,g]
	3727	3706	3701	3701	3737	3715	3710	3710		[d,e,g]	3714[d,e,f,g]

The vibrational states are grouped by polyads.

AVQZ basis set.

harmonic CCSD(T)-F12a/AVTZ, from Ref.155, and anharmonic B2PLYP/AVQZ force fields.

Refs.: [c]155, [d]160, [e]161, [f]162, [g]159.

Table A2

Non-vanishing cubic energy derivatives with respect to and their symmetry relations.99

Symmetry	w.r.t.	Qm		Group	N
		A1	any cs	I, II
		B1		I, II
		B2		I, II
		A1	cs = ct	I, II
		A2	cs = ct	I, II
		B1		I, II
		B2		I, II
				I, II
				II
				I, II
				II
				I, II
				II
				I, II
				II
				I, II
				II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table 16.

DSPT2 and DCPT2 treatments of resonances deserve some considerations. DSPT2 results coincide with their GVPT2 counterparts for all the states that are not affected by resonances. Conversely, DCPT2 provides values equal to GVPT2 ones just for the states that do not contain excitations on degenerate normal modes 1 and 2 (i.e., and), while the energies for the states and, which should also be unaffected by the resonance, are underestimated. In the perturbative treatment, these states do not involve resonant terms because those present in the elements of are exactly erased by those in when the summations in eq. (30) are performed. DSPT2 reproduces correctly this behavior, while the DCPT2 results are slightly different due to a noncomplete cancellation of the transformed resonant terms.

DSPT2 reproduces well the energies of the states involved in 2-dimensional polyads, while the results are not satisfactory for energies involved in larger dimensionality polyads. This is due to the approximation of treating the interactions terms by simplified two-state interacting matrices, then losing in DSPT2 the simultaneous interactions between more than two states. Despite this, DSPT2 can be used to estimate the energies for the fundamental states, since the latter are usually involved in at most 2-dimensional polyads.

Shifting to longer chain linear systems, the results for acetylene and diacetylene are shown in Figure 1 and Table 7, respectively. Acetylene is a well-known system, for which fundamentals, first overtones, combination bands, and l-doublings have been largely studied in the literature. The results for the vibrational frequencies calculated at the MP2 and B2PLYP levels, with the AVTZ basis set, and B3LYP, with the SNSD basis set, are graphically reported in Figure 1, together with the results obtained with the hybrid CCSD(T)/B2PLYP scheme. For each wavenumber value, the series of five marks corresponds, from left to right side, to VPT2, DCPT2, HDCPT2 and DSPT2, and HDSPT2 results. In line with our previous comments, the deviations from experimental values are mainly due to the harmonic part. This error is strongly reduced with hybrid schemes, which yield very good results. The perturbative correction reproduces well the partial lifting of the zeroth-order degeneracy, as can be observed for and, as well as, for and. Moreover, the inclusion of l-doubling is necessary to lift the degeneracy between and and to obtain accurate energies for the combination energies involving degenerate normal modes. For all electronic methods, no first-order resonances are found with Martin's test. Therefore, the purely perturbative VPT2 approach gives good results, slightly improved with the DSPT2 and DCPT2 methods. This is due to the approximate inclusion of higher-order perturbative terms in the treatment of the possibly resonant terms.

Figure 1

Deviations of harmonic ω and anharmonic ν wavenumbers from experimental values (the origin of the y axis) for acetylene (in cm−1). Experimental values are reported in the x axis at the bottom and the corresponding assignment at the top. The series of four values for each anharmonic frequency stands for, from left to right, VPT2, DCPT2, HDCPT2, DSPT2, and HDSPT2 treatments for possibly resonant terms. Computational methods: MP2 and B2PLYP with AVTZ basis set and B3LYP with SNSD. CCSD(T)/A'CVQZ harmonic and anharmonic frequencies from Table 5 of Ref.163. In the hybrid method, the harmonic frequencies are from CCSD(T)/A'CVQZ and the anharmonic force-field from B2PLYP/AVTZ calculations. Experimental values are taken from Ref.163 for fundamental frequencies, and from Ref.59 for overtones and combination bands. MAE stands for mean absolute error.

Table 7

Experimental and computed harmonic ω and anharmonic ν fundamental wavenumbers for diacetylene (in cm−1)

State	Symm.	B3LYP[a]		B2PLYP[b]		HYBRID[b]		Expt.
		ω	ν	ω	ν	ω	ν	ν
	Σg	3466	3343	3477	3352	3463	3338	3332
		2278	2238	2234	2189	2243	2197	2189
		915	901	908	890	894	872	872
	Σu	3467	3344	3478	3353	3465	3339	3334
		2111	2078	2064	2028	2064	2027	2022
	∏g	659	647	645	638	636	627	626
		529	522	507	512	485	491	483
	∏u	665	654	651	640	640	628	628
		237	237	231	232	221	222	220

The vibrational states are indicated as. DFT calculations were done in conjunction of the AVTZ basis set. Within the hybrid scheme, the harmonic wavenumbers, obtained at the AE-CCSD(T)/cc-pCVQZ level, were taken from Ref.61, and the anharmonic force-field calculated in this work at the B2PLYP/AVTZ level. The experimental values were taken from Refs.61 and63.

VPT2 values, no Fermi resonances identified with Martin's test.

GVPT2 values, one weakly interaction between and states.

High-resolution infrared and Raman spectra of C2H2 reported in the literature show the presence of fairly weak couplings between vibrational levels of the same symmetry due to second-order resonances.59,164–167 The 2-2 resonances between the two degenerate normal modes of acetylene were first reported by Huet and coworkers for 12C2D2.165 In their work, the off-diagonal interaction energies between and, and between and, which involve respectively and, are expressed with the and terms (see Table 2 in Ref.165). It has been found that those resonances are particularly relevant for the isotopomers of acetylene, whose two bending vibrations are very close in energy. Furthermore, the need to account for these interactions appears crucial in the study of the highly excited trans-bend levels in 12C2H2, observed by Field and coworkers using the stimulated emission pumping technique.168 Our results obtained at the MP2/AVTZ, B3LYP/SNSD and B2PLYP/AVTZ levels show a very good agreement with those of Huet et al. (see Table 6). Another case of interacting states, between and, was first considered by Mills.167 In Mills' formalism, the interacting energy is reported in Table 1 of Ref.167 as. This coupling ought to be considered in all symmetric isotopes of C2H2, in particular for 13C2H2. Even in this case, the agreement between our computational results and the observed values is remarkable.

Table 6

VPT2 second-order 2-2 interactions (Darling-Dennison) for 12C2H2 and 12C2D2 (in cm−1)

	MP2[a]	B3LYP[b]	B2PLYP[a]	Expt.
12C2H2
	102.0	98.8	100.7	—
	−53.2	−50.0	−52.4	−49.0[c]
				−52.4[d]
				−51.5[e]
12C2D2
	52.3	21.5	25.6	—
	−6.2	−14.1	−8.2	−8.0[f]
	1.5	0.1	1.0	0.4[g]
	269.9	302.8	287.7	—
	−25.7	−22.8	−25.0	−23.9[c]

Basis sets: [a] AVTZ, [b] SNSD. [c] term in Ref.167, [d] term in Ref.166, [e] term in Ref.59, [f] term in Ref.165, [g] term in Ref.165.

Diacetylene has been extensively studied from both experimental and theoretical points of view, because of its prevalence in hydrocarbon combustion and pyrolysis and is known to be present in the interstellar medium and in the atmospheres of several planets and moons of our solar system.61,63,169 The fundamental frequencies for diacetylene have been calculated at the B3LYP/AVTZ and B2PLYP/AVTZ levels, and with the hybrid scheme, where the harmonic frequencies obtained at the AE-CCSD(T)/cc-pCVQZ level61 are coupled with the B2PLYP/AVTZ force-field. The results are reported in Table 7. For this system, Martin's test reveals a weak interaction due to for B2PLYP and hybrid calculations, which is not found for B3LYP computations. The B2PLYP result for (890 cm−1), calculated with the GVPT2 approach, is in better agreement with the experimental data (872 cm−1) than the B3LYP result (901 cm−1), where the interaction term between the and states is treated at the perturbative level (VPT2). As expected, the hybrid values show a very good agreement with the observed ones.

From medium to large symmetric top systems

The wavenumbers calculated at the B2PLYP/AVTZ level for the fundamental, first overtones and combination bands for cyclopropane, which is an oblate symmetric top belonging to the symmetry point group, are reported in Table 8. Also in this case, the states are ordered by polyads. Martin's test identifies for this system three weak Fermi resonances, related to the interaction between and and and, and one tight Fermi resonance, involving and states. GVPT2, DCPT2, and DSPT2 results are reported in Table 8, together with the VPT2 values. HDCPT2 and HDSPT2 give results equal to DCPT2 and DSPT2, respectively, and are, therefore, omitted. The agreement with the experimental values is good for most of the energies, and both DCPT2 and DSPT2 show good results for the states not affected by resonances, as well as, for the states involving resonant interaction terms. Some discrepancies are found for, for which all methods slightly overestimate the experimental value, and, that is underestimated by the theoretical results with respect to the experimental one. VPT2 reproduces well the energy of, but slightly overestimates. GVPT2, which treats variationally the interaction between the latter two states, overestimates the energy of, whereas that of is in agreement with the experimental value. For this case, DCPT2 and DSPT2 reproduce well the energy of, while for the overestimation is similar to that of VPT2. At variance, the results are very good for the combination states involving the excitations of the normal modes labeled as 10 and 14.

Table 8

Harmonic ω, anharmonic ν wavenumbers for cyclopropane (in cm−1)

State	Symm.	ω
		3163	3042	3040	3041	3046	3027
		3061	2993	2983	2982	2954	—
		1531	1502	1497	1498	1515	1499
		1487	1471	1478	1475	1459	1461
		1218	1191	1191	1191	1191	1189
		1162	1129	1129	1129	1129	1127
		1095	1072	1072	1072	1072	1067
		3254	3108	3108	3108	3108	3102
		863	860	860	860	860	854
		3154	3006	3005	3005	3016	3019
		3016	2909	2924	2926	2907	—
		1486	1422	1441	1441	1446	1440
		1487	1515	1505	1495	1491	1480
		1056	1030	1030	1030	1030	1028
		887	854	854	854	854	868
		3233	3087	3087	3087	3087	3082
		1219	1194	1194	1194	1194	1191
		744	742	744	742	742	738
		1775	1714	1714	1714	1714	1727
		1775	1690	1690	1690	1690	1734
		2150	2097	2097	2097	2097	2090
		1838	1814	1817	1814	1814	1805
		1799	1772	1775	1779	1772	1766
		1799	1771	1774	1771	1771	1767
		1799	1772	1775	1773	1772	—

Computed values at the B2PLYP/AVTZ level.

The vibrational states are indicated as. Observed values were taken from Ref.170.

Note that the l-doubling between and has not been taken into account in the experimental values.

As shown above, the hybrid method allows to reduce the computational costs leading to satisfactory results. Table 9 shows the fundamental frequencies for benzene obtained with the hybrid model. Benzene is an oblate symmetric top ( symmetry), which has been widely studied in the literature by both Raman and infrared spectroscopy.38,171–174 In the hybrid computation, the harmonic frequencies have been calculated at the CCSD(T)/ANO4321′ level,175 and the anharmonic force field at the B3LYP/SNSD level. In Table 9, the fundamental frequencies at the B3LYP/SNSD level are also reported. B3LYP/SNSD calculations show a qualitatively good agreement with the experimental values for the majority of the frequencies. Martin's test identifies two weak type II Fermi resonances, the first affecting and states, the second and, and a slightly stronger one, involving and. The latter resonance leads to wrong VPT2 results for (3143 cm−1), that shows a discrepancy of about 100 cm−1 with respect to the observed value (3057 cm−1). At variance, the coupling between and is small, and the VPT2 result for (1604 cm−1) is closer to the observed value (1601 cm−1) than the GVPT2 one (1588 cm−1). The result of the DSPT2 and DCPT2 treatments (1599 cm−1) is also very good. Some discrepancies are present also for; VPT2, DSPT2 and DCPT2 (3070 cm−1) overestimate the reference value (3047 cm−1), while the opposite is true for GVPT2 (3029 cm−1). This frequency is close to the experimental value in the hybrid models, showing once again that the error is mainly due to the unsatisfactory treatment of the harmonic part. In the hybrid method, the two vibrational states and are still affected by resonance, showing similar results to those obtained by full DFT calculations. On the other hand, Martin's test does not identify the resonance affecting the and states in the hybrid case, because of the differences between the CCSD(T) and DFT harmonic frequencies. Moreover, two new weak couplings are identified, the first involving and, the second and. Consequently, is not variationally treated and shows coincident values for all methods (1598 cm−1), that is in good agreement with the observed one (1601 cm−1), while and are very satisfactory in all VPT2, DSPT2, DCPT2 and GVPT2 approaches. These considerations show that a good description of the harmonic frequencies is also important to identify correctly the resonant terms affecting the system. In this case as well, HDCPT2 and HDSPT2 treatment of resonances have been omitted from Table 9 since they are equivalent to DCPT2 and DSPT2.

Table 9

Computed harmonic ω and experimental and calculated anharmonic fundamental wavenumbers ν for benzene (in cm−1)

		B3LYP/SNSD					HYBRID					Expt.
State	Symm.	ω					ω					ν
		1011	997	997	997	997	1003	989	989	989	989	993
		3195	3054	3055	3055	3054	*3210	3069	3070	3070	3073	3074
		1375	1349	1349	1349	1349	*1380	1348	1351	1351	1350	(1350)
		717	692	692	692	692	709	684	684	684	684	(707)
		1015	980	980	980	980	1009	974	974	974	974	(990)
		864	842	842	842	842	865	843	843	843	843	847
		616	612	612	612	612	611	607	607	607	607	608
		1193	1179	1179	1179	1179	1194	1179	1179	1179	1179	1178
		*1635	1604	1599	1599	1588	1637	1598	1598	1598	1598	1601
		3169	3008	3008	3008	3008	3183	3023	3022	3022	3023	3057
		688	673	673	673	673	687	673	673	673	673	674
		1013	1009	1009	1009	1009	1020	1016	1016	1016	1016	(1010)
		*3159	3143	3069	3069	2996	*3173	3105	3076	3076	3009	(3057)
		1169	1158	1158	1158	1158	1163	1152	1152	1152	1152	1150
		1349	1323	1323	1323	1323	1326	1304	1302	1302	1304	1309
		1056	1038	1038	1038	1038	1056	1038	1038	1038	1038	1038
		1509	1479	1479	1479	1479	1509	1479	1479	1479	1479	1484
		*3185	3073	3067	3069	3029	*3200	3083	3080	3081	3040	3047
		411	402	402	402	402	406	397	397	397	397	398
		987	968	968	968	968	985	966	966	966	966	976
MAE			13	9	9	12		10	9	9	9

The vibrational states are indicated as. In the hybrid method, the harmonic frequencies are calculated at the CCSD(T)/ANO4321′ level, from Table 1 of Ref.175, and the anharmonic force field at the B3LYP/SNSD one.

The experimental values are from Ref.38.

The values in parentheses have not been observed directly but have been deduced from combination bands.

The frequencies treated as resonant (DVPT2/GVPT2) are indicated with a *.

MAE stands for Mean Absolute Error.

Following Amat's rule, and l-doublings are found to be non-null for benzene. The B3LYP/SNSD results for the R and S constants are shown in Table 10, together with the values calculated at the B3LYP/TZ2P level, taken as benchmark from Ref.172. Note that in Ref.172, R and S are reported as and. In both sets of results the resonances are treated at the DVPT2 level. The agreement between the two series of data is remarkable.

Table 10

R and S l-type doublings for C6H6, (in cm−1)

Const.	Modes		This work	Lit.	Const.	Modes		This work	Lit.
S	7	6	0.10	—	S	18	8	0.60	0.64
S	8	6	0.19	0.27	S	18	9	−1.50	−1.56
R	8	7	0.18	—	S	18	10	−10.09	−10.26
S	9	6	−0.11	—	R	18	16	0.36	0.36
R	9	7	0.32	0.26[a]	R	18	17	−0.26	−0.39
R	9	8	−0.72	−0.78	S	19	6	1.51	1.92
S	10	6	0.61	0.71	R	19	7	0.03	−0.09
R	10	7	0.04	—	R	19	8	−0.60	−0.39
R	10	8	0.68	0.70	R	19	9	−0.21	−0.20
R	10	9	−1.66	−1.84	R	19	10	0.03	—
R	16	6	−0.01	—	S	19	16	0.69	0.68
S	16	7	−0.13	—	S	19	17	−0.33	−0.33
S	16	8	−0.53	−0.49	S	19	18	0.03	0.04
S	16	9	0.62	0.64	S	20	6	−0.71	−0.35
S	16	10	0.33	0.34	R	20	7	0.05	—
R	17	6	−0.09	—	R	20	8	−0.24	−0.21
S	17	7	0.05	—	R	20	9	−0.28	−0.30
S	17	8	−0.71	−0.68	R	20	10	0.75	0.85
S	17	9	0.13	0.10	S	20	16	0.19	0.23
S	17	10	−0.44	−0.47	S	20	17	0.04	—
R	17	16	−0.40	−0.46	S	20	18	0.77	0.80
R	18	6	0.59	0.77	R	20	19	−0.50	−1.21
S	18	7	−0.02	—

Calculations at the B3LYP/SNSD level, with resonant terms treated within the DVPT2 approach. The reference values are calculated at the B3LYP/TZ2P level, from Table 6 of Ref.172.

Note that in the reference the values are reported as and.

indicates that the value corresponds with the one reported between parentheses in Ref.172.

Moving to larger systems, the importance of taking into account the anharmonicity appears clearly in Tables11 and 12. In the first Table, both the harmonic and anharmonic computational results for triphenylamine are compared with the observed frequencies. Triphenylamine has a D3 three-bladed propeller structure, with a planar central NCCC moiety (see Fig. 2), and has found applications in different fields, including for instance photoconductors and semiconductors.177–180 With 96 vibrational normal modes, the determination of the complete anharmonic force field for this system is computationally very expensive even at the DFT level. However, within the reduced-dimensionality approach, it is possible to calculate the anharmonic corrections for a small selection of vibrational energies of interest. If the harmonic energy of the latter are well separated from the energies of the vibrations ignored in the anharmonic treatment, the cubic and quartic forces involving normal modes of both sets can be assumed to be negligible. In Table 11, the anharmonic corrections have been applied to fundamental vibrational states having harmonic wavenumbers larger than 3000 cm−1 which correspond to the CH stretchings region. The calculation has been done at the B3LYP/6-31G* level, and the resonances have been treated with the DSPT2 method. In Table 11, the empirical fundamental frequencies, obtained scaling the B3LYP/AVTZ harmonic frequencies by a factor of 0.986 (see Ref.181), are also reported, together with the experimental results, measured by FTIR spectroscopy of triphenylamine monomers isolated in an argon matrix.181 The inclusion of anharmonic effects leads to a significantly better agreement between the theoretical and experimental results with respect to the scaled values.

Table 11

Fundamental vibrational wavenumbers for triphenylamine (in cm−1)

	B3LYP/6-31G*		Scaled	Expt.
Symm.	ω	ν[a]	ν[b]	ν[c]
E	3182	3029	3127	3016
				3043
A2	3190	3072	3135	3067
E	3190	3074	3135
E	3205	3070	3150
E	3217	3069	3159
A2	3214	3096	3158	3096
E	3214	3097	3157	3107

Anharmonic correction computed within the reduced dimensionality approach (see text), applying the DSPT2 method for resonances.

Harmonic values at the B3LYP/AVTZ level and scaled with a factor equal to 0.986, from Ref.181.

Observed values from Ref.181.

Table 12

Computed harmonic ω, GVPT2 anharmonic ν, and experimental wavenumbers for staggered and eclipsed ferrocene (in cm−1)

	B3LYP[a]		B3PW91[b]			B3LYP[a]		B3PW91[b]		Expt. [c]
Symm.	ω	ν	ω	ν	Symm.	ω	ν	ω	ν	ν
	453	440	488	475		448	437	482	470	480
E1′	466	457	501	492		436	427	470	460	496
	828	815	830	829		827	813	830	820	816
E1′	845	837	857	840		844	841	855	839	840
E1′	1022	1000	1025	1006		1021	1002	1026	1006	1012
	1130	1112	1142	1125		1131	1113	1142	1126	1112
E1′	1449	1415	1451	1418		1450	1417	1451	1419	1416
E1′	3239	3106	3245	3116		3238	3107	3245	3115	3106
	3250	3116	3256	3126		3249	3118	3256	3126	—

SNSD/aug-LANL2DZ basis set.

SNSD/m6-31G basis set.

Observed values from Ref.176.

Figure 2

Medium-sized symmetric top systems of interest.

As a last example, we report the results for ferrocene, an organometallic compound of great interest in biotechnologies and nanotechnologies, with important applications of its derivatives in catalysis, molecular electronics, polymer chemistry, nonlinear optical, and solar engineering.182–187 Its geometry has been studied by several theoretical methods and shows a sandwich structure with the metal situated between two parallel cyclopentadienyl rings. A small energy barrier separates the staggered and eclipsed rotational orientation of the two rings (see Fig. 2), with an energy difference of 0.9 kcal mol−1 from gas phase electron diffraction measurements.188–190 In gas phases calculations, the eclipsed conformer is a global minimum, whereas the staggered conformer is a saddle point with an imaginary frequency. In a recent study, a quite good agreement was obtained between the harmonic vibrational frequencies of ferrocene calculated at the B3LYP/m6-31(d) level and the observed values.190 A noticeable improvement in the theoretical results is obtained by taking into account the anharmonicity. From B3LYP calculations, with the hybrid SNSD/aug-LANL2DZ basis set as discussed in the computational details section, the anharmonic fundamental wavenumbers show a quantitative agreement with the experimental ones, especially for the range above 800 cm−1, where vibrations involving C and H atoms are excited. The lowest wavenumbers (480 and 496 cm−1) are due to the excitations of vibrational modes involving the metal. The latter are better described by the B3PW91 functional, coupled with the SNSD/m6-31G basis set. It is noteworthy that B3PLYP and B3PW91 anharmonic corrections are not significantly different, showing that the discrepancies between the observed and B3LYP values are due again to deficitary description of the harmonic vibrations associated to Fe.

Rotovibrational interaction terms

The importance of including the vibrational corrections to the rotational constants to achieve both accurate rotational energies and accurate geometrical parameters has been widely illustrated in the literature.191–195 The vibrational corrections α, the equilibrium B and ground vibrational state rotational constants B0 for the linear systems HCP, OCS, and C2H2 obtained at different computational levels are reported in Table 13, together with the equilibrium quartic distortion constants. Like for vibrational energies, the discrepancies with the reference values are mainly associated to B, while the differences show a lower sensitivity to the change of the computational level. On the other hand, the centrifugal distortion constants have a slightly larger variability. Accurate values for the latter are obtained from calculations involving accurate geometrical parameters and equilibrium rotational constants. The rotational constants for the symmetric top C3H6 at the B3LYP/SNSD level are shown in Table 14. Those results are compared with experimental and theoretical data, the latter obtained at the highly reliable CCSD(T) level. For this system, and are affected by a Coriolis resonance, due to (see Table 8) and the two associated states, that is, and, are not prevented by symmetry to interact. On the other hand, and are not affected by resonance, since vanishes for symmetry reasons. On the other hand, the total rotovibrational corrections to the rotational constants are not affected by resonances. The B3LYP/SNSD calculation shows good results also for equilibrium quartic distortion constants.

Table 13

Vibrational corrections α, rotational constants B and B0 and quartic and sextic distortion constants for HCP, OCS, and C2H2 (in cm−1)

	MP2	B3LYP[a]	B2PLYP	Best theo.	Expt.
HCP
	−0.00046[b]	−0.00058	−0.00047[b]	−0.00045[c]	−0.00045[d]
	0.00409	0.00346	0.00388	0.00362	0.00362[e]
	0.00313	0.00307	0.00313	0.00322	0.00318[d]
Be	0.65822	0.66174	0.66702	0.66931	—
B0	0.65506	0.65905	0.66400	0.66634	0.66633[c]
	−0.00316	−0.00269	−0.00302	−0.00297	—
	0.70884	0.64986	0.69155	0.70545[c]	0.70420[c,f]
OCS
	−0.00037[b]	−0.00035	−0.00036[a]	−0.00035[g]	−0.00034[g]
	0.00056	0.00068	0.00064	0.00066	0.00067
	0.00125	0.00121	0.00125	0.00123	0.00125
Be	0.20219	0.20030	0.20247	—	—
B0	0.20166	0.19971	0.20188	—	—
	−0.00053	−0.00059	−0.00059	—	—
	0.04063	0.04164	0.04223	0.04203	0.04270
C2H2
	−0.00137[h]	−0.00145	−0.00135[h]	−0.00141[i]	−0.00135[i,j]
	−0.00201	−0.00221	−0.00218	−0.00220	−0.00223
	0.00653	0.00556	0.00609	0.00584	0.00588
	0.00579	0.00575	0.00586	0.00601	0.00618
	0.00693	0.00672	0.00697	0.00686	0.00690
Be	1.16883	1.17463	1.18369	1.18245	—
B0	1.16259	1.16928	1.17775	1.17670	1.17665
	−0.00624	−0.00535	−0.00594	−0.00575	—
	1.58695	1.46786	1.56394	1.5902	1.627[f,k]
	0.89529	1.11004	1.08214	1.2631	1.6[f,k]

Basis sets: [a] SNSD; [b] AVQZ; [h] AVTZ.

Refs.: [c]196; [d]198; [e]199; [g]157; [i]163; [j]59; [k]197.

[f] Ground state observed values.

Table 14

Rotational constants and quartic distortion constants for C3H6 (in cm−1)

	B3LYP[a]	Best theo.	Expt.
Be	0.67034	0.67807[b]	—
B0	0.66342	0.67104	0.67024[b,c,d]
	0.00692	0.00702
Ce	0.41890	0.42414[b]	—
C0	0.41405	0.41914	0.41770[b,c]
			0.41881[d]
	0.00485	0.00500
	0.95346	0.93288[e]	0.96668[d,e,f,g]
	−1.23376	−1.18929	−1.24924
	0.47968	0.45564	0.48619

[a] SNSD basis set.

Refs.: [b]200; [c]201; [d]170; [e]202; [f]203.

[g] Ground state observed values.

Thermodynamics

If the fundamental, overtone and combination energies have to be handled with care because of resonances, it has been shown in the theoretical section that the ZPVE is not affected. Both harmonic and anharmonic ZPVEs of linear (HCN, CO2, C2H2), and symmetric top molecules (PH3, ClCH3, FCH3) are shown in Table 15. On overall, the mean anharmonic correction with respect to the harmonic ZPVE is about 0.4% for CO2, 1% for HCN, 1.2% for C2H2, and 1.4% for the symmetric top systems. It is noteworthy that for all these molecules the magnitudes of the anharmonic corrections are little affected by the choice of the computational method and the basis set, at least in the present cases. From the ZPVE and the anharmonic fundamental energies, a comparison with the experimental thermodynamic data can be achieved by the SPT model.40,53 The calculated and experimental absolute entropies at 298.15 K and 1 atm are also reported in Table 15. Under those thermodynamic conditions, the absolute entropies calculated with all methods available to treat the resonances lead to very close results. Compared to accurate experimental values,the inclusion of anharmonic corrections in the calculated thermodynamic values improves the accuracy of the results by about J mol−1 K−1.

Table 15

Comparison of computed and experimental harmonic (H) and anharmonic (A) ZPVE (in KJ mol−1) and absolute entropies at 298.15 K and 1 atm (in J mol−1 K−1), for linear and symmetric top molecules

	MP2[a]		B3LYP[b]		B2PLYP[a]		Best theo.	Expt.
	H	A	H	A	H	A
HCN
ZPVE	41.42	41.05	42.46	42.02	42.31	41.90	[c]	41.61[c]
Δ		−0.37		−0.44		−0.41
S	201.79	201.82	201.99	202.20	201.39	201.50	—	201.83[d,g]
CO2
ZPVE	30.18	30.08	30.49	30.36	30.24	30.12	[e]	—
Δ		−0.10		−0.13		−0.12
S	213.88	213.95	213.74	213.78	213.72	213.79	—	213.69[d,g]
C2H2
ZPVE	69.65	68.93	70.73	69.58	70.59	69.85	[e]	—
Δ		−0.72		−1.15		−0.74
S	200.82	200.92	200.08	201.14	200.06	200.34	—	200.85[d,g]
PH3
ZPVE	64.39	63.59	62.21	61.37	63.49	62.67	[f]	—
Δ		−0.80		−0.84		−0.82
S	209.90	209.98	210.25	210.33	209.97	210.05	—	210.13[d,g]
ClCH3
ZPVE	100.75	99.36	98.73	97.34	99.79	98.43	—	—
Δ		−1.39		−1.39		−1.36
S	233.77	233.92	234.42	234.58	234.06	234.22	—	234.26[d,g]
FCH3
ZPVE	104.80	103.37	102.61	101.18	103.72	102.30	—	—
Δ		−1.43		−1.43		−1.42
S	222.52	222.62	222.75	222.85	222.59	222.69	—	222.73[d,g]

Δ's are the anharmonic corrections.

Basis sets: [a] AVTZ, [b] SNSD.

Refs.: [c]156; [d]204; [e]155; [f]55.

[g] The tabulated values have been lowered by 0.11 J mol−1 K−1, to pass from the original 1 MPa values to 1 MPa (see “reference part” in [204]).

Conclusion

The VPT for rotovibrational energies and thermodynamic functions for asymmetric, symmetric and linear top systems has been revised and fully generalized to allow for the treatment of both minima and first-order saddle points of the PES. A particular attention has been devoted to the treatments of off-diagonal elements of the Hamiltonian and the perturbative equations in the presence of resonances. Previous strategies for dealing with first-order resonances (i.e., GVPT2, DCPT2, and HDCPT2) have been generalized and a new treatment (i.e., DSPT2 and its hybrid counterpart HDSPT2), has been presented and validated. A versatile implementation has been included in the Gaussian package.

Several case studies ranging from triatomic to large molecular systems have been explicitly treated by different QM approaches to fully validate the computational tool. The results show that the perturbative developments are very effective and reasonably accurate, and can be applied easily to DFT and DFT/CCSD(T) hybrid levels in conjunction with medium sized basis sets, and with reduced-dimensionality schemes. The latter approximations are particular appealing when dealing with medium- to large-molecules, allowing the inclusion of anharmonicity also in the cases otherwise unpractical due to prohibitive computational cost.

Table A3

Non-vanishing quartic energy derivatives K, K and K with respect to and their symmetry relations.60

Symmetry	w.r.t.	Qm	Qn		Group	N
		A1 or A2		any cs	I, II
		B1 or B2		any cs	I, II
		A1 or A2		cs = ct	I, II
		B1 or B2		cs = ct	I, II
		A1 or A2		cs = ct	II
		B1 or B2		cs = ct	II
		A1	A1	any cs	I, II
		A2	A2
		B1	B1	any cs	I, II
		B2	B2
		A1	B1		I, II
		A2	B2
		A1	B2		I, II
		A2	B1
		A1	A1	cs = ct	I, II
		A2	A2
		B1	B1	cs = ct	I, II
		B2	B2
		A1	B1		I, II
		A2	B2
		A1	A2	cs = ct	I, II
		A2	A1
		B1	B2	cs = ct	I, II
		B2	B1
		A1	B2		I, II
		A2	B1

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

Table A4

Non-vanishing quartic energy derivatives K, K and K with respect to and their symmetry relations.60

Symmetry	w.r.t.	Qm		Group	N
		A1		I, II
		B1		I, II
		A2		I, II
		B2		I, II
		A1		I, II
		B1		I, II
		A1		I, II
		B1		I, II
		A2		I, II
		B2		I, II
		A2		I, II
		B2		I, II
		A1		I, II
		B1		I, II
		A1		I, II
		B1		I, II
		A2		I, II
		B2		I, II
		A2		I, II
		B2		I, II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

Table A5

Non-vanishing quartic energy derivatives K and K with respect to and their symmetry relations60

Symmetry	w.r.t.		Group	N
			I, II
			I, II
			II
			I, II
			I, II
			I, II
			I, II
			II
			II
			II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

Table A6

Non-vanishing quartic energy derivatives K with respect to and their symmetry relations60

Symmetry	w.r.t.		Group	N
			I, II
			I, II
			I, II
			I, II
			I, II
			II
			II
			II
			II
			II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

Table A7

Non-vanishing quartic energy derivatives K with respect to and their symmetry relations60

Symmetry	w.r.t.		Group	N
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II
			I, II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.

Table A8

Non-vanishing quartic energy derivatives K with respect to and their symmetry relations60

Symmetry	w.r.t.		Group	N
			II
			II
			II
			II
			II
			II
			II
			II
			II
			II

c is the subscript labelling the degenerate representation of mode i, for example c = 1 for E or E1, c = 2 for E2, etc. p is a non zero integer number and N indicates the order of the principal symmetry axis. For I and II Group classification see Table A1.