Vibrational second-order perturbation theory based on curvilinear coordinates: Thermochemical applications

Mendolicchio, M.; Barone, V.

doi:10.1063/5.0252006

This work improves and extends a general and robust workflow for the computation of anharmonic vibrational frequencies to thermodynamic functions, paving the way toward the study of large flexible molecules. The key new feature is the extension of closed-form expressions for both zero-point vibrational energies and partition functions to second-order vibrational perturbation theory based on curvilinear internal coordinates. The use of curvilinear coordinates enables the reduction of couplings between different degrees of freedom, enriching the arsenal of existing vibrational approaches, and can lead to effective, low-dimensional linear-scaling models. The accuracy of the results obtained for some prototypical systems paves the way toward the systematic use of this new implementation in the study of molecules containing a few dozen atoms, as exemplified by the test cases of a molecular motor, a nucleoside, and two hormones.

I. INTRODUCTION

The challenge of performing accurate simulations of large systems has been a fundamental issue in theoretical and computational chemistry since its inception and is likely to remain a central concern in the future.¹ On the one hand, state-of-the-art quantum chemical (QC) methods are producing results with an accuracy rivaling that of experimental data for small semirigid systems.^2–4 For molecules not plagued by strong static correlation, the coupled cluster model, including single, double, and (perturbatively) triple excitations [CCSD(T)],⁵ is considered the golden standard of contemporary quantum chemistry, provided that core-valence correlation and complete basis set extrapolation are properly handled.¹ Unfortunately, the unfavorable scaling of this model with the size of the system significantly limits its application to large, flexible systems, which are of great fundamental and technological interest. This limitation is particularly pronounced in scenarios involving rugged potential energy surfaces (PESs), characterized by numerous low-energy minima separated by energy barriers of different heights. At the opposite end of the spectrum, while simplified QC models can accommodate much larger systems, they do so at the expense of substantial reductions in computational accuracy.

The latest developments of hardware and software are producing accurate low-scaling methods for the computation of electronic energies.^6–9 However, thermochemical studies also require, in addition to accurate electronic energies, partition functions that depend on precise molecular structures, zero-point vibrational energies (ZPVEs), and finite-temperature (FT) contributions. These contributions benefit only marginally from the developments mentioned above, which have not yet produced analytical energy derivatives. Therefore, methods rooted in density functional theory (DFT) offer at present the only viable route, with the latest double hybrid functionals yielding remarkably accurate molecular structures and harmonic force fields.^10,11 From another point of view, the standard harmonic oscillator (HO) model produces significant errors, mainly in the high-frequency regime (overestimating ZPVE contributions) and low-frequency regime (overestimating entropy contributions). The first issue can be effectively solved by the second-order vibrational perturbation theory (VPT2),^12–16 provided that the presence of resonances is properly taken into account.^17–20 The second problem requires the separation between different degrees of freedom, which can be performed safely only if the couplings between different modes are sufficiently small. Although Cartesian coordinates are usually employed in this connection, nonreactive potential energy surfaces can also be described in terms of a unique underlying bond pattern retained by all considered structures. These topological features allow for the automatic definition of the basic variables (bond lengths, valence, and dihedral angles),²¹ leading to a good description of both small- and large-amplitude internal motions (SAM and LAM, respectively), while minimizing intermode couplings at the same time.

This work is placed in the framework of a general computational strategy utilizing curvilinear internal coordinates,²² which we believe provides a robust, accurate, and effective solution to a variety of challenges in molecular modeling. The main steps of the strategy can be summarized as follows:

Black-box generation of a complete set of curvilinear internal coordinates, distinguishing between hard and soft degrees of freedom within the molecular system.²³
Exploration of soft degrees of freedom by means of constrained geometry optimizations. This is achieved through a fast semi-empirical method,²⁴ guided by suitable algorithms,^25–27 aimed at identifying low-lying energy minima.
Refinement of the initial results by full geometry optimizations with DFT methods.²⁸
Identification of inter-conversion paths between adjacent minima and removal of energy minima that can effectively relax into more stable structures under the relevant experimental conditions.²²
Final refinement of structures, electronic energies and, possibly, equilibrium values of physical–chemical properties by different variants of the new Pisa composite schemes (PCS).^29–31
Computation of ZPVEs and thermodynamic functions by the simple perturbation theory (SPT) model³² employing VPT2 data for SAMs and different effective one-dimensional models for individual LAMs.³³

Although these steps have recently been validated and incorporated into a general computational tool,^23,34 the calculation of ZPVEs and thermodynamic functions beyond the harmonic oscillator model requires further refinement. In fact, current implementations of VPT2 predominantly rely on rectilinear coordinates, which are not well-suited for accurate description of soft degrees of freedom such as torsions or inversions, resulting in significant couplings between SAMs and LAMs. The recent development of a generalized VPT2 engine that employs curvilinear internal coordinates^35–37 can address these challenges. The main aim of the present work is to extend this implementation to the computation of ZPVEs and thermodynamic functions. We will demonstrate that many couplings between different vibrational modes can be effectively neglected when using curvilinear internal coordinates, allowing for a clear-cut separation of LAMs from their small-amplitude counterparts, which can then be treated as effective one-dimensional modes.

In particular, the one-dimensional hindered rotor model proposed by Schlegel and co-workers,³⁸ when combined with the VPT2 framework for other vibrational modes, yields remarkably accurate vibrational entropies for both semi-rigid and flexible molecules.³⁹ At the same time, uncoupled non-periodic large-amplitude motions can be described accurately by one-dimensional variational approaches employing the same local information underlying the VPT2 model.⁴⁰

This paper is organized as follows: in Sec. II, after a short recap of VPT2 employing different kinds of coordinates, we generalize the well-known expressions of ZPVE and partition functions to curvilinear coordinates and introduce a general framework for the variational treatment of LAMs employing the same information needed for performing VPT2 computations. Then, after validating the proposed procedure with several small molecules for which accurate experimental results are available, we will analyze larger molecules, including a molecular motor, a nucleoside, and two hormones. The main conclusions and most promising perspectives are summarized in Sec. VI.

II. THEORETICAL BACKGROUND

In this section, the theoretical background underlying a general second-order perturbative treatment of vibro-rotational motions of isolated molecules encompassing both rectilinear (VPT2-R) and curvilinear (VPT2-C) coordinates will be summarized. Next, a resonance-free expression of ZPVE within the VPT2-C framework will be derived. The attention will then be focused on the development of a black-box algorithm for the approximate treatment of Fermi resonances (FRs) aimed at the calculation of thermochemical parameters without any threshold selection or variational treatment. Finally, the theoretical framework will be extended to the so-called reduced-dimensionality approach,^41,42 in which the number of modes investigated at the VPT2 level can be lower than that of normal modes. More in detail, the case of a single degree of freedom separated from the rest of modes will be explicitly considered, paving the route toward the extension of the VPT2-C framework to the treatment of first-order transition states (TSs).

A. Vibro-rotational problem in curvilinear coordinates

The vibro-rotational motions of isolated molecules have been treated by different methods,^43–65 ranging from perturbation theory^12–19 to variational and quasi-variational models,^66–73 not to speak of time-dependent approaches. In this work, we will employ the VPT2 model based on curvilinear coordinates (VPT2-C), whose main aspects are summarized in this section, with the aim of providing the theoretical background for the new developments presented in the present work.

As is well known, a molecule composed of N_a atoms has N normal modes of vibration (N = 3N_a − 6 for a nonlinear molecule and N = 3N_a − 5 for a linear one). Within the Born–Oppenheimer approximation,⁷⁴ the nuclear Hamiltonian can be expressed in the following way:

H = - \frac{ℏ^{2}}{2} \sum_{a = 1}^{N_{a}} \frac{\nabla_{a}^{2}}{m_{a}} + V,

(1)

where ℏ is the reduced Planck constant,

\nabla_{a}^{2}

is the Laplacian operator with respect to the Cartesian coordinates of the nucleus a in a laboratory-fixed frame (LF), m_a is the corresponding mass, and V is the potential energy. After placing the center of mass (COM) of the molecule at the origin of the axes, a body-fixed frame (BF) can be introduced, with the new set of Cartesian coordinates

x = {x_{1}, \dots, x_{N_{a}}}

being related to their LF counterparts by a set of N_r = 3 (2 for linear molecules) Euler angles

ϕ = {ϕ_{1}, \dots, ϕ_{N_{r}}}

⁠. The remaining N degrees of freedom represent molecular vibrations, and they can be expressed in terms of a non-redundant set of internal curvilinear coordinates s = {s₁, …, s_N}. Therefore, the general expression of the Hamiltonian can be derived from Eq. in terms of the generalized coordinates ξ = {X^COM, ϕ, s},⁷⁵^,

H = - \frac{ℏ^{2}}{2} \sum_{i = 1}^{3 N_{a}} \sum_{j = 1}^{3 N_{a}} {\tilde{G}}^{1 / 4} \frac{\partial}{\partial ξ_{i}} G^{- 1 / 2} G_{i j} \frac{\partial}{\partial ξ_{j}} {\tilde{G}}^{- 1 / 4} + V,

(2)

where X^COM is the vector that collects the Cartesian coordinates of the COM in the LF, while G_ij and

\tilde{G}

represent an element and the determinant of the well-known matrix G, respectively. The latter can be conveniently obtained as the inverse of the covariance metric tensor

\bar{G}

⁠, whose elements can be expressed within the t-vector formulation,^46,76

{\bar{G}}_{i j} = \sum_{a = 1}^{N_{a}} m_{a} t_{i, a} \cdot t_{j, a},

(3)

where t_i,a = ∂x_a/∂ξ_i. Given that the overall translations can be factored out exactly, the expression of the vibro-rotational Hamiltonian H^vr can be obtained by limiting the summations of Eq. to N′ + N_r, with N′ ≤ N being the number of active vibrations, which is different from the total number of normal modes when the calculations of reduced dimensionality are considered, while N′ = N for full calculations. Following previous work,^77,78 and in particular, the formulation proposed by Podolsky⁷⁹ and Lauvergnat,⁴⁷ the operator H^vr can be written as follows:

H^{v r} = - \frac{ℏ^{2}}{2} \sum_{i = 1}^{N^{'} + N_{r}} \sum_{j = 1}^{N^{'} + N_{r}} \frac{\partial}{\partial ξ_{i}} G_{i j} \frac{\partial}{\partial ξ_{j}} + V_{g} + V,

(4)

where the so-called extra-potential term V_g has been introduced,

V_{g} = \frac{ℏ^{2}}{32} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} [\frac{G_{i j}}{{\tilde{G}}^{2}} \frac{\partial \tilde{G}}{\partial s_{i}} \frac{\partial \tilde{G}}{\partial s_{j}} - 4 \frac{\partial \tilde{G}}{\partial s_{i}} (G_{i j} \frac{\partial \tilde{G}}{\partial s_{j}})] .

(5)

Then, the normal coordinates are obtained through Wilson’s GF method,⁸⁰

G^{v} F L = L Λ,

(6)

where G^v is the vibrational diagonal block of the G matrix, Λ is the diagonal matrix collecting the squared angular frequencies, and L is the matrix of eigenvectors. The dimensionless normal coordinates q = {q₁, …, q_N} can be obtained from their mass-weighted counterparts Q = {Q₁, …, Q_N},

q_{i} = γ_{i}^{1 / 2} Q_{i} = \sqrt{\frac{2 π c ω_{i}}{ℏ}} Q_{i},

(7)

where c is the light speed in vacuum, ω_i is the harmonic wave number (in cm⁻¹) of the normal mode i, and Q = L⁻¹s. The expansion of H^vr in terms of normal coordinates leads to the following expression:

H^{v r} = \sum_{f} \sum_{g} = H_{f g},

(8)

where f and g represent the degrees of the vibrational and rotational operators, respectively. The purely vibrational terms of Eq. up to the second order define the VPT2 Hamiltonian,

H = H^{(0)} + λ H^{(1)} + λ^{2} H^{(2)} .

(9)

In Eq. , λ is the perturbation parameter (set to 1 at the end of the calculation), while the perturbative terms contain kinetic

(g)

and potential

(f)

contributions,

\begin{gathered} H^{(0)} = \frac{1}{2} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} ω_{i} (q_{i}^{2} + p_{i}^{2}), \\ H^{(1)} = \frac{1}{6} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} \sum_{k = 1}^{N^{'}} (f_{ijk} q_{i} q_{j} q_{k} + 3 g_{i j, k}^{v} p_{i} q_{k} p_{j}), \\ H^{(2)} = \frac{1}{24} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} \sum_{k = 1}^{N^{'}} \sum_{l = 1}^{N^{'}} (f_{ijkl} q_{i} q_{j} q_{k} q_{l} + 6 g_{i j, k l}^{v} p_{i} q_{k} q_{l} p_{j}) + V_{g}, \end{gathered}

(10)

with p_i being the conjugate momentum of q_i (p_i = −i∂/∂q_i). The symbols g and V_g will be adopted hereafter to indicate the matrix G and the extra-potential term expressed in cm⁻¹, while

f_{i j k l \dots}

and

g_{i j, k l \dots}

represent the derivatives of potential energy and the matrix g^v with respect to the dimensionless normal coordinates,⁴⁵

f_{i j k \dots} = \frac{\partial^{n} V}{\partial q_{i} \partial q_{j} \partial q_{k} \partial q_{l} \dots}, g_{i j, k \dots} = \frac{\partial^{n} g_{i j}^{v}}{\partial q_{k} \partial q_{l} \dots} .

(11)

The application of second-order Rayleigh–Schrödinger (RSPT)⁸¹ or canonical Van Vleck (CVPT)^18,82,83 perturbation theory yields the same expression of the energy of a generic vibrational state R,

ε_{R} = χ_{0} + \sum_{i = 1}^{N^{'}} ω_{i} (v_{R, i} + \frac{1}{2}) + \sum_{i = 1}^{N^{'}} \sum_{j = i}^{N^{'}} χ_{i j} (v_{R, i} + \frac{1}{2}) (v_{R, j} + \frac{1}{2}),

(12)

where v_R,i is the quantum number of the normal mode i, and the χ₀ term together with the symmetric χ matrix have been introduced. Previous work was devoted to the determination of transition properties that require only an analytical expression of the χ matrix. In fact, the transition energy from R to S vibrational states is simply given by

\begin{align} ν_{S, R} & = ε_{S} - ε_{R} = \sum_{i = 1}^{N^{'}} ω_{i} Δ v_{S R, i} \\ + \sum_{i = 1}^{N^{'}} \sum_{j = i}^{N^{'}} χ_{i j} [Δ_{2} v_{S R, i j} + \frac{1}{2} (Δ v_{S R, i} + Δ v_{S R, j})], \end{align}

(13)

where Δ₂v_SR,ij = v_S,iv_S,j − v_R,iv_R,j and Δv_SR,i = v_S,i − v_R,i. Although the VPT2-C model^84–86 and its higher-order counterparts^87–89 have been analyzed by different authors, only recently, a fully automatized construction of curvilinear coordinates accompanied by compact general expression based on effective tensors has been derived, unifying VPT2-R and VPT2-C models.³⁵ As a result, any code implementing VPT2-R can be easily extended to its VPT2-C counterpart using the following expressions:

16 χ_{i i} = η_{iiii} - \frac{σ_{iii}^{2} / 2 + 9 ρ_{iii}^{2} / 2}{3 ω_{i}} - \frac{1}{2} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} [\frac{4 ρ_{jii}^{2}}{ω_{j}} + \frac{σ_{iij}^{2}}{2 ω_{i} + ω_{j}} - \frac{ρ_{iij}^{2}}{2 ω_{i} - ω_{j}}],

(14)

\begin{align} 4 χ_{i j} & = η_{iijj} - \frac{1}{2} [\frac{σ_{iij}^{2}}{2 ω_{i} + ω_{j}} + \frac{ρ_{iij}^{2}}{2 ω_{i} - ω_{j}}] - \frac{1}{2} [\frac{σ_{jji}^{2}}{2 ω_{j} + ω_{i}} + \frac{ρ_{jji}^{2}}{2 ω_{j} - ω_{i}}] \\ - \frac{ρ_{iii} ρ_{ijj}}{ω_{i}} - \frac{ρ_{jjj} ρ_{jii}}{ω_{j}} - \frac{1}{2} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} [\frac{σ_{ijk}^{2}}{ω_{i} + ω_{j} + ω_{k}} - \frac{ρ_{ijk}^{2}}{ω_{i} + ω_{j} - ω_{k}} \\ + \frac{ρ_{ikj}^{2}}{ω_{i} - ω_{j} + ω_{k}} - \frac{ρ_{jki}^{2}}{ω_{i} - ω_{j} - ω_{k}} + \frac{2 ρ_{kii} ρ_{kjj}}{ω_{k}}] + (C_{i j}^{χ}), \end{align}

(15)

where

η_{ijkl} = f_{ijkl} + g_{i j, k l}^{v} + g_{k l, i j}^{v},

(16)

σ_{ijk} = f_{ijk} - (g_{i j, k}^{v} + g_{i k, j}^{v} + g_{j k, i}^{v}),

(17)

ρ_{ijk} = f_{ijk} - (g_{i j, k}^{v} - g_{i k, j}^{v} - g_{j k, i}^{v}),

(18)

and

(C_{i j}^{χ})

is the Coriolis contribution to the χ_ij element, which is present only in the VPT2-R model and is given by

(C_{i j}^{χ}) = \frac{4 (ω_{i}^{2} + ω_{j}^{2})}{ω_{i} ω_{j}} \sum_{τ}^{x, y, z} B_{τ}^{e q} {ζ_{i j, τ}}^{2},

(19)

where

B_{τ}^{e q}

and ζ_ij,τ are the equilibrium rotational constants and Coriolis interaction term coupling modes i and j along the τ principal axis of inertia. In the VPT2-C model, the

(C_{i j}^{χ})

term is lacking, although it survives in the VPT2-R model, where, instead, the derivatives of the g matrix vanish (see the supplementary material for further details).

As is well-known, Eqs. and suffer from the presence of vanishing denominators when ω_i ≈ 2ω_j or ω_i ≈ ω_j + ω_k (referred to as FRs of type I and II, respectively). Most strategies aimed at determining FRs^90,91 are based on a two-step procedure,the first step being the estimation of the energetic proximity characterizing the interacting states,

| ω_{i} - (ω_{j} + ω_{k}) | \leq Δ ω^{1 - 2},

(20)

where Δω¹⁻² is a predefined threshold and j can be equal to k. In the second step, the “weight” of the potentially resonant term is calculated. In the present work, the test proposed by Martin et al.⁹⁰ has been adopted, in which the term is marked as resonant whenever the deviation of the perturbative contributions to the energies resulting from an ad hoc, variational model is beyond a defined threshold K¹⁻²,

\frac{ρ_{jki}^{4}}{64 {(1 + δ_{j k})}^{2} | ω_{i} - (ω_{j} + ω_{j}) |^{3}} \geq K^{1 - 2} .

(21)

Once the resonant terms have been identified, they can be simply removed, the resulting model being referred to as deperturbed VPT2 (DVPT2). In order to bypass the accuracy loss due to the elimination of interactions between interacting states, a more refined model, namely, the generalized VPT2 (GVPT2),^16,17 can be employed. Within this scheme, a variational matrix is built, whose diagonal entries are represented by the DVPT2 energies. The off-diagonal elements correspond to the pairs of interacting states linked by FRs, whose expression is provided by the corresponding matrix element of the canonical Van Vleck transformed Hamiltonian

\tilde{H}

⁠,

⟨ v_{R} + 1_{i} | \tilde{H} | v_{R} + 1_{j} + 1_{k} ⟩ = ρ_{jki} \sqrt{\frac{(v_{R, i} + 1) (v_{R, j} + 1) (v_{R, k} + 1 + δ_{j k})}{8 (1 + 3 δ_{j k})}} .

(22)

As already mentioned for the matrix χ, the VPT2-R results are recovered, neglecting the derivatives of the matrix g in Eqs. (21) and (22). Then, the eigenvalues arising from the diagonalization of the different blocks of the variational matrix provide the GVPT2 energies, while the eigenvectors collect the coefficients of the harmonic wave functions necessary to build the variational states.

Although the GVPT2 model produces remarkably accurate vibrational states for semi-rigid molecules, thermochemical (and kinetic) studies require lower accuracy but general procedures providing a continuous treatment of resonances for different stationary points and/or quantum chemical models. Therefore, we extended to curvilinear coordinates an approximate model with these characteristics^92,93 and further validated its accuracy.

B. Zero-point vibrational energy

As already mentioned, in previous work, the attention was focused on the calculation of vibrational transition energies^35,36 as well as vibrational and vibro-rotational averages of molecular properties.³⁷ However, a proper characterization of the ZPVE provides invaluable information for different applications, ranging from chemical kinetics and thermochemistry to conformational studies. For this reason, the first objective of the present study has been the development (and application) of an analytical expression of ZPVE within the VPT2-C framework. A first derivation of χ₀ was proposed by Isaacson,^85,86 also addressing the treatment of doubly degenerate vibrations, needed for the study of symmetric and linear tops. In the present context, a new derivation focusing solely on asymmetric tops has been performed, since non-Abelian symmetry groups can be effectively managed using the unified VPT2 model outlined in Ref. 94,

χ_{0} = χ_{0}^{V} + χ_{0}^{T} + χ_{0}^{\times},

(23)

with

χ_{0}^{V}

⁠,

χ_{0}^{T}

⁠, and

χ_{0}^{\times}

being the potential, kinetic, and cross term, respectively,

\begin{align} 64 χ_{0}^{V} & = \sum_{i = 1}^{N^{'}} f_{iiii} - \sum_{i = 1}^{N^{'}} \frac{7 {f_{iii}}^{2}}{9 ω_{i}} + \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \frac{3 {f_{iij}}^{2} ω_{j}}{4 ω_{i}^{2} - ω_{j}^{2}} \\ - \frac{8}{3} \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} \frac{{f_{ijk}}^{2} ω_{i} ω_{j} ω_{k}}{Δ_{ijk}}, \end{align}

(24)

\begin{align} 64 χ_{0}^{T} & = 10 \sum_{i = 1}^{N^{'}} g_{i i, i i} - \sum_{i = 1}^{N^{'}} \frac{3 {g_{i i, i}}^{2}}{ω_{i}} + 8 \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} g_{i j, i j} \\ + \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \frac{3 {g_{i i, j}}^{2} ω_{j} - 24 g_{i i, j} g_{i j, i} ω_{i} + 12 {g_{i j, i}}^{2} ω_{j}}{4 ω_{i}^{2} - ω_{j}^{2}} \\ - \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}} [\frac{8 g_{i j, k}^{2} ω_{i} ω_{j} ω_{k}}{Δ_{ijk}} \\ \frac{+ 4 g_{i j, k} g_{j k, i} ω_{j} (ω_{j}^{2} - ω_{i}^{2} - ω_{k}^{2})}{Δ_{ijk}} \\ \frac{+ 4 g_{i j, k} g_{i k, j} ω_{i} (ω_{i}^{2} - ω_{j}^{2} - ω_{k}^{2})}{Δ_{ijk}}], \end{align}

(25)

\begin{align} 64 χ_{0}^{\times} & = \sum_{i = 1}^{N^{'}} \frac{2 f_{iii} g_{i i, i}}{3 ω_{i}} + \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \frac{f_{iij} (24 g_{i j, i} ω_{i} - 6 g_{i i, j} ω_{j})}{4 ω_{i}^{2} - ω_{j}^{2}} \\ + \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} \frac{8 ω_{k} (ω_{k}^{2} - ω_{i}^{2} - ω_{j}^{2}) f_{ijk} g_{i j, k}}{Δ_{ijk}}, \end{align}

(26)

where

Δ_{ijk} = ω_{i}^{4} + ω_{j}^{4} + ω_{k}^{4} - 2 ω_{i}^{2} ω_{j}^{2} - 2 ω_{i}^{2} ω_{k}^{2} - 2 ω_{j}^{2} ω_{k}^{2} .

(27)

Equations – are in full agreement with those reported by Isaacson, but an additional analysis based on the partial fraction decomposition allowed us to achieve once again a very compact expression valid for both the VPT2-C and VPT2-R variants,

\begin{array}{l} 64 χ_{0} & = \sum_{i = 1}^{N^{'}} a_{iiii} - \sum_{i = 1}^{N^{'}} \frac{5 {σ_{iii}}^{2} / 2 + 9 {ρ_{iii}}^{2} / 2}{9 ω_{i}} \\ + \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \frac{3}{2} [\frac{{ρ_{iij}}^{2}}{2 ω_{i} - ω_{j}} - \frac{{σ_{iij}}^{2}}{2 ω_{i} + ω_{j}}] \\ - \frac{1}{3} \sum_{i = 1}^{N^{'}} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} [\frac{{σ_{ijk}}^{2}}{ω_{i} + ω_{j} + ω_{k}} - \frac{ρ_{ikj}^{2}}{ω_{i} - ω_{j} + ω_{k}} \\ - \frac{{ρ_{ijk}}^{2}}{ω_{i} + ω_{j} - ω_{k}} + \frac{{ρ_{jki}}^{2}}{ω_{i} - ω_{j} - ω_{k}}] + V_{g} + (C_{0}^{χ}), \end{array}

where the

a_{iiii}

elements are given by

a_{iiii} = f_{iiii} + 2 g_{i i, i i} + 8 \sum_{j = 1}^{N^{'}} g_{i j, i j},

(28)

and

(C_{0}^{χ})

is the Coriolis contribution to χ₀, only present in the VPT2-C variant,

(C_{0}^{χ}) = - 16 \sum_{τ}^{x, y, z} B_{τ}^{e q} [1 + \sum_{i = 1}^{N^{'} - 1} \sum_{j = i + 1}^{N^{'}} {ζ_{i j, τ}}^{2}] .

(29)

The last contribution to χ₀ is due to the extra-potential term, whose expression in wave numbers,

\begin{align} V_{g} & = \frac{1}{32} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} [\frac{g_{i j}}{{\tilde{g}}^{2}} \frac{\partial \tilde{g}}{\partial q_{i}} \frac{\partial \tilde{g}}{\partial q_{j}} - 4 \frac{\partial}{\partial q_{i}} (\frac{g_{i j}}{\tilde{g}} \frac{\partial \tilde{g}}{\partial q_{j}})] \\ = \frac{1}{32} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} [\frac{5 g_{i j}}{{\tilde{g}}^{2}} \frac{\partial \tilde{g}}{\partial q_{i}} \frac{\partial \tilde{g}}{\partial q_{j}} - 4 \frac{g_{i j, i}}{\tilde{g}} \frac{\partial \tilde{g}}{\partial q_{j}} - 4 \frac{g_{i j}}{\tilde{g}} \frac{\partial^{2} \tilde{g}}{\partial q_{i} \partial q_{j}}] \end{align}

(30)

can be further rewritten by making use of the following notation:

{\tilde{g}}_{i j \dots} = \frac{1}{\tilde{g}} {(\frac{\partial^{n} \tilde{g}}{\partial q_{i} \partial q_{j} \dots})}_{e q}

(31)

to obtain a far more manageable equation,

V_{g} = \frac{1}{32} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} [5 g_{i j} {\tilde{g}}_{i} {\tilde{g}}_{j} - 4 g_{i j, i} {\tilde{g}}_{j} - 4 g_{i j} {\tilde{g}}_{i j}] .

(32)

In general, the extra-potential term is well-approximated by its value for the reference configuration,^84,85 which has the following expression in the case of stationary points:

V_{g}^{e q} = \frac{1}{32} \sum_{i = 1}^{N^{'}} [5 ω_{i} {\tilde{g}}_{i}^{2} - 4 ω_{i} {\tilde{g}}_{i i} - 4 \sum_{j = 1}^{N^{'}} g_{i j, i} {\tilde{g}}_{j}],

(33)

where the equivalence

g_{i j}^{v, e q} = ω_{i} δ_{i j}

has been exploited, so that only the diagonal second-order derivatives of

\tilde{g}

are needed. The terms

{\tilde{g}}_{i j \dots}

can be directly evaluated from the derivatives of the matrix g or by numerical differentiation of analytical determinants. For the sake of readability, a dedicated section has been included in the Appendix. The ZPVE is the energy associated with a null vector of quantum numbers (v_R,i = 0 with i = 1, …, N′), which is different from χ₀ and is denoted by the symbol ɛ₀,^19,83,95

ε_{0} = χ_{0} + \sum_{i}^{N^{'}} \frac{ω_{i}}{2} + \frac{1}{4} \sum_{i = 1}^{N^{'}} \sum_{j = i}^{N^{'}} χ_{i j} .

(34)

A combination of Eqs. , , , and leads to the first expression (to the best of our knowledge) of the resonance-free ZPVE obtained within the VPT2-C formulation. In fact, the following expression gathers both VPT2-C and VPT2-R models:

\begin{align} ε_{0} & = \sum_{i = 1}^{N^{'}} \frac{ω_{i}}{2} + \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} \frac{b_{iijj}}{32} - \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} \sum_{k = 1}^{N^{'}} [\frac{{σ_{ijk}}^{2}}{48 (ω_{i} + ω_{j} + ω_{k})} \\ + \frac{ρ_{kii} ρ_{kjj}}{32 ω_{k}}] + V_{g}^{e q} + (C_{0}), \end{align}

(35)

which is of course deprived of FRs and where an additional effective tensor has been introduced,

b_{ijkl} = f_{ijkl} + 2 g_{i j, k l} + 4 g_{i k, j l} .

(36)

The term (C₀) is the Coriolis contribution to ɛ₀, only present in the VPT2-R variant,

(C_{0}) = - \sum_{τ}^{x, y, z} B_{τ}^{e q} \sum_{i = 1}^{N^{'}} \sum_{j = 1}^{N^{'}} {ζ_{i j, τ}}^{2} \frac{{(ω_{i} - ω_{j})}^{2}}{ω_{i} ω_{j}} .

(37)

C. Soft degrees of freedom

As is well known, accurate descriptions of highly anharmonic degrees of freedom cannot be obtained by the VPT2-R model. This limitation stems from its reliance on a fourth-degree polynomial representation of the PES in terms of rectilinear coordinates. Previous studies have shown that a more effective approach involves leveraging the interplay between different methodologies (e.g., perturbative and variational techniques) to address each vibrational mode with the most suitable method.

From this standpoint, a curvilinear description of internal motions can significantly reduce the magnitude of anharmonic terms involving different degrees of freedom, a feat that is challenging when using rectilinear displacements. To establish this approach, two preliminary steps are essential. First, the VPT2-C formalism must be fully developed for obtaining an accurate description of SAMs. Second, it is crucial to examine how the removal of the “special” coordinates affects the other frequencies and molecular properties.

With the purpose of illustration, one can consider a system presenting a single degree of freedom F poorly described at the VPT2 level (N′ = N − 1), with q_F being the corresponding normal coordinate. Within a reduced-dimensionality scheme neglecting the overall contribution of F, the terms χ_FF and χ_iF (i ≠ F) vanish, while χ_ii and χ_ij (with i, j ≠ F) are deprived of the contributions defining the T matrix,

T_{i i} = - \frac{1}{2} [\frac{4 ρ_{F i i}^{2}}{ω_{F}} + \frac{σ_{i i F}^{2}}{2 ω_{i} + ω_{F}} - \frac{ρ_{i i F}^{2}}{2 ω_{i} - ω_{F}}],

(38)

T_{i j} = - \frac{1}{2} [\frac{σ_{i j F}^{2}}{ω_{i} + ω_{j} + ω_{F}} - \frac{ρ_{i j F}^{2}}{ω_{i} + ω_{j} - ω_{F}},

(39)

+ \frac{σ_{i F j}^{2}}{ω_{i} - ω_{j} + ω_{F}} - \frac{ρ_{j F i}^{2}}{ω_{i} - ω_{j} - ω_{F}} + \frac{2 ρ_{F i i} ρ_{F j j}}{ω_{F}}] .

(40)

Since the main problems faced by VPT2 in the treatment of this class of vibrations could be related not only to the coordinate choice but also to the fourth-order polynomial expansion of the potential energy, the overall elements of the χ matrix can be quite similar for the VPT2-C and VPT2-R variants, suggesting that the former approach does not provide any benefit. However, any shift of the offending contributions from inter- to intra-mode couplings (as expected when one goes from rectilinear to curvilinear coordinates) would permit the use of reduced-dimensionality variational treatments, which can be easily implemented and extended to large molecules. This strategy can be particularly useful in the presence of out-of-scale vibrational couplings that could seriously affect the vibrational properties of the LAM, as well as those of the remaining modes. Conversely, the information provided by those couplings can be preserved and accounted for within the perturbative calculation involving the SAMs, while only the diagonal terms linked to the LAM are removed and used in a subsequent variational framework. In the latter case, all the elements of the T matrix vanish, with the exception of χ_iF,

- \frac{ρ_{F F F} ρ_{F j j}}{ω_{F}} .

(41)

It is important to point out that when the diagonal cubic terms (⁠ $f_{F F F}$ and $g_{F F, F}^{v}$ ⁠) vanish (e.g., for symmetry reasons), this approach becomes equivalent to a full perturbative calculation for the SAMs subspace.

Once the coordinate of interest has been separated, the last step is represented by the variational solution of a one-dimensional problem, where the corresponding Hamiltonian matrix elements,

H_{v_{R, F}}^{v_{R, F} + n} = ⟨ v_{R, F} | H | v_{R, F} + n ⟩

(42)

can be easily obtained through the second-quantization formalism, leading to a set of compact expressions. In this context, the one-dimensional vibrational Hamiltonian H^F can be written in the following way:

H^{F} = T^{F} + V^{F},

(43)

where the kinetic and potential energies can be expressed as

\begin{align} T^{F} & = \frac{1}{2} ω_{F} p_{F}^{2} + \frac{1}{2} g_{F F, F} p_{F} q_{F} p_{F} + \frac{1}{4} g_{F F, F F} p_{F} q_{F}^{2} p_{F} \\ + \frac{1}{6} g_{F F, F F F} p_{F} q_{F}^{3} p_{F} + \frac{1}{24} g_{F F, F F F F} p_{F} q_{F}^{4} p_{F} + V_{g}^{F}, \end{align}

(44)

\begin{align} V^{F} & = \frac{1}{2} ω_{F} q_{F}^{2} + \frac{1}{6} f_{F F F} q_{F}^{3} + \frac{1}{24} f_{F F F F} q_{F}^{4} \\ + \frac{1}{120} f_{F F F F F} q_{F}^{5} + \frac{1}{720} f_{F F F F F F} q_{F}^{6} . \end{align}

(45)

The symbol

V_{g}^{F}

denotes the contribution to the extra-potential term arising only from the mode F,

V_{g}^{F} = \frac{1}{32} [5 ω_{F} {\tilde{g}}_{F}^{2} - 4 g_{F F, F} {\tilde{g}}_{F}^{2} - 4 ω_{F} {\tilde{g}}_{F F}] .

(46)

The expressions of the elements of interest, namely,

H_{v_{R, F}}^{v_{R, F} + n}

(with n = 1, …, 6) are

\begin{align} H_{v_{R, F}}^{v_{R, F}} & = \frac{1}{2} ω_{F} (S + 1) (v_{R, F} + \frac{1}{2}) + V_{g}^{F} \\ + \frac{f_{F F F F} + 2 g_{F F, F F}}{16} ({v_{R, F}}^{2} + v_{R, F} + \frac{1}{2}) \\ + \frac{g_{F F, F F}}{8} + \frac{g_{F F, F F F F}}{16} (v_{R, F} + \frac{1}{2}) \\ + \frac{f_{F F F F F F} + 6 g_{F F, F F F F}}{288} ({v_{R, F}}^{3} + \frac{3}{2} {v_{R, F}}^{2} + 2 v_{R, F} + \frac{3}{4}), \end{align}

(47)

\begin{align} H_{v_{R, F}}^{v_{R, F} + 1} & = [\frac{f_{F F F F F} + 4 g_{F F, F F F}}{48 \sqrt{2}} ({v_{R, F}}^{2} + 2 v_{R, F} + \frac{3}{2}) \\ + \frac{g_{F F, F F F}}{4 \sqrt{2}}] B_{R, F} (1) + [\frac{f_{F F F} + g_{F F, F}}{4 \sqrt{2}}] B_{R, F} {(1)}^{3}, \end{align}

(48)

\begin{align} H_{v_{R, F}}^{v_{R, F} + 2} & = [\frac{ω_{F}}{4} (S - 1) + \frac{f_{F F F F}}{24} (v_{R, F} + \frac{3}{2}) \\ + \frac{f_{F F F F F F} + 2 g_{F F, F F F F}}{384} ({v_{R, F}}^{2} + 3 v_{R, F} + 3) \\ + \frac{g_{F F, F F F F}}{16}] B_{R, F} (2), \end{align}

(49)

H_{v_{R, F}}^{v_{R, F} + 3} = [\frac{f_{F F F} - 3 g_{F F, F}}{12 \sqrt{2}} + \frac{f_{F F F F F} - 4 g_{F F, F F F}}{96 \sqrt{2}} (v_{R, F} + 2)] B_{R, F} (3),

(50)

\begin{align} H_{v_{R, F}}^{v_{R, F} + 4} & = [\frac{f_{F F F F F F} - 10 g_{F F, F F F F}}{960} (v_{R, F} + \frac{5}{2}) \\ + \frac{f_{F F F F} - 6 g_{F F, F F}}{96}] B_{R, F} (4), \end{align}

(51)

\begin{align} H_{v_{R, F}}^{v_{R, F} + 5} & = [\frac{f_{F F F F F} - 20 g_{F F, F F F}}{480 \sqrt{2}}] B_{R, F} (5), \end{align}

(52)

\begin{align} H_{v_{R, F}}^{v_{R, F} + 6} & = [\frac{f_{F F F F F F} - 30 g_{F F, F F F F}}{5760}] B_{R, F} (6) . \end{align}

(53)

Although this work deals with local minima of the PES, the above-mentioned equations are also applicable to transition states. In particular, the symbol S denotes the sign of the harmonic force constant (i.e., S = 1 for a minimum and S = −1 for a TS), while B_R,F(n) is a factor depending on the vibrational quantum numbers,⁴⁰

B_{R, F} (n) = \sqrt{\prod_{i = 1}^{n} (v_{R, F} + i)} .

(54)

D. Black-box treatment of Fermi resonances for thermochemistry

In the previous sections, the VPT2-C framework has been generalized to the calculation of the ZPVE, as well as to the treatment of a coordinate belonging to a different class, including the case of an imaginary frequency. However, the calculation of the energy levels still suffers from the presence of FRs, with this issue also affecting the calculation of the vibrational partition functions and, therefore, of thermochemical properties. Thrular and Isaacson proposed the so-called simple perturbation theory (SPT) in which all normal modes are assumed to be decoupled to build vibrational partition functions,³²

Q_{vib} (T) = \frac{e^{- h c ε_{0} / k_{B} T}}{\prod_{\begin{array}{c} i = 1 \\ (i \neq F) \end{array}}^{N} (1 - e^{- h c ε_{1_{i}} / k_{B} T})},

(55)

but the ZPVEs (ɛ₀) and fundamental energy levels

(ε_{1_{i}})

are evaluated at the VPT2 level, the latter requiring the detection and proper treatment of resonances. To this end, the matrix χ can be partitioned into three contributions,

χ = N + S + T,

(56)

where T is given in Eq. , N collects contributions of Eqs. and that are not affected by possible resonances;

16 N_{i i} = η_{iiii} - \frac{{σ_{iii}}^{2} / 2 + 9 {ρ_{iii}}^{2} / 2}{3 ω_{i}} - \frac{1}{2} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} [\frac{4 {ρ_{jii}}^{2}}{ω_{j}} + \frac{{σ_{iij}}^{2}}{2 ω_{i} + ω_{j}}],

(57)

\begin{align} 4 N_{i j} & = η_{iijj} - \frac{1}{2} [\frac{{σ_{iij}}^{2}}{2 ω_{i} + ω_{j}} + \frac{{σ_{jji}}^{2}}{2 ω_{j} + ω_{i}}] - \frac{ρ_{iii} ρ_{ijj}}{ω_{i}} - \frac{ρ_{jjj} ρ_{jii}}{ω_{j}} \\ - \frac{1}{2} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} [\frac{{σ_{ijk}}^{2}}{ω_{i} + ω_{j} + ω_{k}} + \frac{2 ρ_{kii} ρ_{kjj}}{ω_{k}}], \end{align}

(58)

and S includes potentially resonant contributions that deserve particular attention,

16 S_{i i} = \frac{1}{2} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N} \frac{{ρ_{iij}}^{2}}{2 ω_{i} - ω_{j}},

(59)

\begin{align} 4 S_{i j} & = - \frac{1}{2} [\frac{{ρ_{iij}}^{2}}{2 ω_{i} - ω_{j}} + \frac{{ρ_{jji}}^{2}}{2 ω_{j} - ω_{i}}] + \frac{1}{2} \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} [\frac{{ρ_{ijk}}^{2}}{ω_{i} + ω_{j} - ω_{k}} \\ - \frac{{ρ_{ikj}}^{2}}{ω_{i} - ω_{j} + ω_{k}} + \frac{{ρ_{jki}}^{2}}{ω_{i} - ω_{j} - ω_{k}}] . \end{align}

(60)

As already mentioned, the most accurate treatment of these contributions is offered by the GVPT2 model, in which these terms are neglected in the perturbative treatment and reintroduced in a successive variational step.^17,20 However, this approach is not well suited for thermochemical applications, as the selection of resonant terms can be different for different methods and/or different stationary points. At the same time, the accuracy required for thermochemical applications is lower than that needed for high-resolution spectroscopy. In this framework, Kuhler et al.⁹² proposed a general approach based on the systematic replacement of all potentially resonant contributions with non-resonant ones. In particular, using an analogy with the simplified Hamiltonian,

H = [\begin{matrix} A - ε & k \\ k & A + ε \end{matrix}],

(61)

whose eigenvalues are

ε^{\pm} = A \pm ε \sqrt{1 + \frac{k^{2}}{ε^{2}}},

(62)

the authors showed that it is possible to rewrite the potentially resonant terms using the following transformation:

\frac{S k^{2}}{2 ε} \to S (\sqrt{k^{2} + ε^{2}} - ε),

(63)

where ɛ corresponds to half the absolute value of the frequency difference, which may cause the singularity, and k² is the absolute value of the constant term and S is the sign (i.e., 1 or −1). Application of Eq. to Eq. and leads to a version known as degeneracy corrected VPT2 (DCPT2) for matrix S,

\begin{gathered} S_{i i}^{D} = \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} O_{iij} [\sqrt{\frac{{ρ_{iij}}^{2}}{32} + \frac{Δ ω_{iij}^{2}}{4}} - \frac{Δ ω_{iij}}{2}], \\ \begin{aligned} S_{i j}^{D} & = - O_{iij} [\sqrt{\frac{{ρ_{iij}}^{2}}{32} + \frac{Δ ω_{iij}^{2}}{4}} - \frac{Δ ω_{iij}}{2}] \\ + \sum_{\begin{array}{c} k = 1 \\ (k \neq i, j) \end{array}}^{N^{'}} \{O_{ijk} [\sqrt{\frac{{ρ_{ijk}}^{2}}{16} + \frac{Δ ω_{ijk}^{2}}{4}} - \frac{Δ ω_{ijk}}{2}] \\ - O_{ikj} [\sqrt{\frac{{ρ_{ikj}}^{2}}{16} + \frac{Δ ω_{ikj}^{2}}{4}} - \frac{Δ ω_{ikj}}{2}] \\ - O_{jki} [\sqrt{\frac{{ρ_{jki}}^{2}}{16} + \frac{Δ ω_{jki}^{2}}{4}} - \frac{Δ ω_{jki}}{2}]\}, \end{aligned} \end{gathered}

(64)

where

O_{ijk} = s i g n (ω_{i} + ω_{j} - ω_{k}),

(65)

Δ ω_{ijk} = | ω_{i} + ω_{j} - ω_{k} |

(66)

has been introduced for the sake of readability. The similarity of the matrix χ for both formulations of VPT2 has allowed a straightforward extension of the DCPT2 scheme to the use of curvilinear coordinates replacing the anharmonic force constants with the elements of the tensors η, σ, and ρ. Within the DCPT2 scheme, the fundamental energies can be written as follows:

(67)

Although this model provides an expression of the χ matrix devoid of singularities, its accuracy is strongly dependent on the magnitude of k and ɛ. For example, when ɛ is small, the conventional VPT2 model becomes unreliable, and its replacement with DCPT2 is recommended. Otherwise, large values of ɛ imply that the VPT2 model yields correct results. A classification of all possible combinations of k and ɛ is reported in Table I.

TABLE I.

Applicability of the VPT2 and DCPT2 models in the function of the values of ɛ and k.

	VPT2		DCPT2
	k ≪ 1	k ≫ 1	k ≪ 1	k ≫ 1
ɛ ≪ 1	✗	✗	✓	✓
ɛ ≫ 1	✓	✓	✓	✗

	VPT2		DCPT2
	k ≪ 1	k ≫ 1	k ≪ 1	k ≫ 1
ɛ ≪ 1	✗	✗	✓	✓
ɛ ≫ 1	✓	✓	✓	✗

Based on the previous discussion, an improved approach can be based on a smooth transition from VPT2 to DCPT2, depending on the nature of the resonant term. This consideration led to the development of the hybrid degeneracy-corrected second-order perturbation theory (HDCPT2),⁹³ in which each quantity

f

is written in the following parametric form:

f^{H} = Λ f + (1 - Λ) f^{D},

(68)

where

f^{H}

⁠,

f

⁠, and

f^{D}

represent the term within HDCPT2, VPT2, and DCPT2 frameworks, respectively.

The parameter Λ depends not only on k and ɛ but also on two additional parameters α and β,

Λ = \frac{\tanh (α [\sqrt{k^{2} ε^{2}} - β]) + 1}{2},

(69)

with the former controlling the transition from VPT2 to DCPT2 and the latter its smoothness. The values of the parameters proposed in Ref. 93 were α = 1.0 and β = 5.0 × 10⁵.

Then, it is possible to apply the transformation introduced in Eq. to the elements of the matrix S,

\begin{gathered} S_{i i}^{H} = Λ S_{i i} + (1 - Λ) S_{i i}^{D} S_{i j}^{H} = Λ S_{i j} + (1 - Λ) S_{i j}^{D} . \end{gathered}

(70)

In analogy to Eq. , the fundamental frequencies within the HDCPT2 framework can be evaluated,

ε_{1_{i}} - ω_{i} = 2 (N_{i i} + S_{i i}^{H} + T_{i i}) + \frac{1}{2} \sum_{\begin{array}{c} j = 1 \\ (j \neq i) \end{array}}^{N^{'}} (N_{i j} + S_{i j}^{H} + T_{i j}),

(71)

and used in conjunction with Eq. for the calculation of the vibrational partition function deprived of singularities. The latter serves as a starting point to derive the thermochemical properties of interest, such as the vibrational contributions to the internal energy (U_v), entropy (S_v), and specific heat (c_v) accounting for anharmonic effects;

\begin{gathered} U_{v} = R \frac{ε_{0}}{k_{B}} + R \sum_{i = 1}^{N^{'}} \frac{h c ε_{1_{i}} / k_{B}}{e^{h c ε_{1_{i}} / k_{B} T} - 1}, \\ S_{v} = R \sum_{i = 1}^{N^{'}} [\frac{h c ε_{1_{i}} / k_{B}}{e^{h c ε_{1_{i}} / k_{B} T} - 1} - \ln (1 - e^{h c ε_{1_{i}} / k_{B} T})], \\ c_{v} = R \sum_{i = 1}^{N^{'}} {[\frac{h c ε_{1_{i}} / k_{B} T}{e^{h c ε_{1_{i}} / k_{B} T} - 1}]}^{2} . \end{gathered}

(72)

III. IMPLEMENTATION

The theoretical developments summarized in Sec. II have been implemented in a standalone program described in previous work,^35,37 including the DCPT2 and HDCPT2 extensions of the VPT2-C workflow, the calculation of ZPVE, and the machinery allowing for one-dimensional vibrational configuration interaction (VCI) calculations including up to sixth-order terms. The program has also been enhanced to include a general thermochemical analysis, which encompasses the calculation of partition functions at the harmonic and SPT levels, together with the derived thermodynamic functions. Such improvements lead, in our opinion, to a versatile tool for spectroscopic and thermochemical studies.

For the sake of clarity, let us summarize the organization of the whole workflow, which basically consists of the following steps:

Geometry optimization and calculation of the force constant matrix through a selected QC program.
Calculation of the harmonic vibrational frequencies and normal modes in curvilinear coordinates, through the standalone program, followed by the generation of the displaced geometries necessary for the finite differences and run of the corresponding single-point calculations by the QC package.
Extraction of the output data, construction of the anharmonic force constants and derivatives of the matrix g through the standalone program.
Anharmonic vibrational analysis, with particular reference to the calculation of transition energies, ZPVE, resonance diagnostics, and thermochemical analysis.

Although interfacing the code with any QC program is straightforward, in the present work, the Gaussian 16 package⁹⁶ has always been used.

IV. COMPUTATIONAL DETAILS

Different DFT methods have been used, in particular the hybrid B3LYP (hereafter B3) functional^97–99 and its double hybrid counterpart revDSD-PBEP86 (hereafter rDSD) functional.¹⁰⁰ The empirical dispersion corrections proposed by Grimme (D3 model with Becke–Johnson damping, hereafter, D3BJ^101,102) have been systematically included. The B3 model has been always employed in conjunction with the split valence polarized 6-31G* and 6-31+G* (from now on referred to as SVP) basis sets.^103–105 The rDSD functional has been employed in conjunction with the 3F12⁻ basis set,²⁹ which is obtained from its cc-pVTZ-F12 counterpart,¹⁰⁶ upon removal of the d function on first-row atoms and replacement of the two f functions on second- and third-row atoms by a single f function taken from the cc-pVTZ basis set.¹⁰⁷

V. RESULTS AND DISCUSSION

A panel of prototypical molecules is used to evaluate the accuracy of the proposed protocol against data issued from experiments and state-of-the-art computational methodologies. The focus will then shift to systems involving large-amplitude motions, taking formamide and its out-of-plane motion as a case study to highlight the interplay between perturbative and variational approaches. Finally, the anharmonic vibrational and thermochemical properties of larger molecules are analyzed, using possibly dual-level methods.

A. Validation

Semi-rigid molecules are the best choices for validation purposes, as in this case different coordinate choices should provide comparable results. In this context, five representative systems have been selected for a detailed study: formaldehyde, oxirane, imidazole, fluoroamine, and vinyl radical (see Fig. 1).

FIG. 1.

View large Download slide

Molecular structures of (a) formaldehyde, (b) oxirane, (c) vinyl radical, (d) fluoroamine, and (e) imidazole.

The ZPVE of these systems has been determined through both the VPT2-R and VPT2-C frameworks at the B3/SVP level of theory. The results collected in Table II offer a breakdown of the overall ZPVE into different contributions.

TABLE II.

Analysis and comparison of the anharmonic ZPVE (in cm⁻¹) of selected rigid systems between Cartesian and curvilinear coordinates at the B3/SVP level of theory.

		VPT2-R			VPT2-C
	Harm.a	Anh.b	W+Cc	Total	Anh.d	ExPot.e	Total
Formaldehyde	5 874	−80	−1	5 793	−77	−4	5 793
Oxirane	12 637	−172	1	12 466	−185	14	12 466
Imidazole	15 620	−189	1	15 432	−192	4	15 432
Fluoroamine	6 058	−101	−1	5 957	−93	−8	5 957
Vinyl radical	8 041	−123	3	7 921	−106	−14	7 921

		VPT2-R			VPT2-C
	Harm.a	Anh.b	W+Cc	Total	Anh.d	ExPot.e	Total
Formaldehyde	5 874	−80	−1	5 793	−77	−4	5 793
Oxirane	12 637	−172	1	12 466	−185	14	12 466
Imidazole	15 620	−189	1	15 432	−192	4	15 432
Fluoroamine	6 058	−101	−1	5 957	−93	−8	5 957
Vinyl radical	8 041	−123	3	7 921	−106	−14	7 921

^a

Harmonic ZPVE.

^b

Anharmonic correction arising from the expansion of the potential energy.

^c

Anharmonic correction arising from both the Coriolis effect and Watson’s terms.

^d

Anharmonic correction arising from the expansion of both the potential and kinetic energies.

^e

Extra-potential contribution evaluated at the reference geometry.

In particular, within the VPT2-R framework, the anharmonic correction has been partitioned into contributions arising solely from the anharmonic force constants and those arising from the additional Coriolis contribution. In the VPT2-C framework, the reported values correspond to the anharmonic contributions of the potential and kinetic energies, together with the extra-potential term evaluated at the equilibrium geometry.

The overall anharmonic correction is very close for different choices of coordinates, with the anharmonic contribution consistently dominating. Further analysis revealed that within the VPT2-C framework, the anharmonic correction issued from the potential energy typically exceeds its kinetic counterpart (see the supplementary material for further details).

The next step of the analysis concerns the determination of thermochemical properties. In this connection, the accuracy of the VPT2-C ZPVEs and HDCPT2-C energies that govern fundamental transitions will be analyzed, as both contributions enter the SPT approximation of the vibrational partition functions. Formaldehyde has been selected as a case study because of its well-known spectroscopic properties (above all the presence of a strong FR), which allow an unbiased assessment of the robustness and reliability of the computational workflow. The VPT2, DCPT2, and HDCPT2 models in conjunction with either rectilinear or curvilinear coordinates have been used for the calculation of the anharmonic ZPVE, fundamental wave numbers, absolute entropy (S), and specific heat at constant pressure (c_p). The results obtained at the rDSD/3F12⁻ level are compared to accurate computed values¹⁰⁸ and available experimental data¹⁰⁹ in Table III.

TABLE III.

Comparison of the experimental and theoretical fundamental wave numbers (in cm⁻¹), absolute entropy, and specific heat [in cal (mol K)⁻¹] of formaldehyde at the rDSD/3F12⁻ level of theory.

		Curvilinear coords.			Rectilinear coords.
Rectilinear state	Symm.	VPT2	DCPT2	HDCPT2	VPT2	DCPT2	HDCPT2	CCa	Expt.b
ZPVE		5785	5785	5785	5785	5785	5785	⋯	⋯
\|1₁⟩	A₁	2785	2786	2785	2785	2786	2785	2784	2783
\|1₂⟩		1755	1755	1755	1755	1755	1755	1735	1746
\|1₃⟩		1505	1506	1506	1505	1506	1506	1495	1500
\|1₄⟩	B₁	1180	1180	1180	1180	1180	1180	1163	1167
\|1₅⟩	B₂	2797	2816	2817	2797	2814	2813	2844	2843
\|1₆⟩		1248	1248	1248	1248	1248	1248	1241	1249
Sc		52.24	52.24	52.24	52.24	52.24	52.24	⋯	52.28d
$c_{p}$ e		8.45	8.45	8.45	8.45	8.45	8.45	⋯	8.46d

		Curvilinear coords.			Rectilinear coords.
Rectilinear state	Symm.	VPT2	DCPT2	HDCPT2	VPT2	DCPT2	HDCPT2	CCa	Expt.b
ZPVE		5785	5785	5785	5785	5785	5785	⋯	⋯
\|1₁⟩	A₁	2785	2786	2785	2785	2786	2785	2784	2783
\|1₂⟩		1755	1755	1755	1755	1755	1755	1735	1746
\|1₃⟩		1505	1506	1506	1505	1506	1506	1495	1500
\|1₄⟩	B₁	1180	1180	1180	1180	1180	1180	1163	1167
\|1₅⟩	B₂	2797	2816	2817	2797	2814	2813	2844	2843
\|1₆⟩		1248	1248	1248	1248	1248	1248	1241	1249
Sc		52.24	52.24	52.24	52.24	52.24	52.24	⋯	52.28d
$c_{p}$ e		8.45	8.45	8.45	8.45	8.45	8.45	⋯	8.46d

^a

GVPT2-R with CCSD(T)/aug-cc-pVTZ force field from Ref. 108.

^b

Reference 109.

^c

Molar entropy at 298.15 K and 1 atm.

^d

Reference 110.

^e

Specific heat at 298.15 K and 1 atm.

As expected, the value of the overall anharmonic ZPVE is very similar for the different VPT2 variants, since it does not involve any resonance. Equal wave numbers of the fundamental transitions are also obtained from standard VPT2-C and VPT2-R computations. Otherwise, a slight discrepancy of the DCPT2 and HDCPT2 results is observed for the |1₅⟩ state, which corresponds to the antisymmetric CH stretching. General trends are comparable in both theoretical frameworks, but there is no direct correspondence of rectilinear anharmonic force constants with elements of the σ and ρ tensors. At the same time, rectilinear and curvilinear k terms in Eq. (63) are not strictly equivalent, leading to different DCPT2 results. Of course, this discrepancy also extends to the HDCPT2 scheme, where the factor Λ introduces an additional variation.

As expected, the transition from the VPT2 model to the DCPT2 and HDCPT2 variants has a non-negligible impact only on the |1₅⟩ state. Even if the HDCPT2 result is closer to its experimental counterpart, the error (about 30 cm⁻¹) is still non-negligible. Formaldehyde highlights the limits of the HDCPT2 approximation. In fact, both the VPT2 and DCPT2 values of the fundamental under consideration are well below the experimental value of 2843 cm⁻¹, implying that even changing the value of the parameter Λ (with 0 ≤ Λ ≤ 1) would not solve the problem. Of course, a methodology exploiting a variational step for the inclusion of resonance effects (GVPT2)¹⁷ is more accurate, leading in this case to a value of 2852 cm⁻¹, which corresponds to a deviation of 9 cm⁻¹ with respect to the experimental data. The comparison of rDSD results with their CCSD(T) counterparts (see columns 5, 8 and 9 of Table III) confirm the remarkable accuracy of this double hybrid functional^10,11 that is applicable to much larger molecules. Note that the seemingly large discrepancy for the state |1₅⟩ (30 cm⁻¹) is reduced to 8 cm⁻¹ when the GVPT2 model is used in all computations.

It is worth mentioning that high-frequency modes have a limited impact on the vibrational partition function and consequently on the resulting thermochemical properties. This is well illustrated by the absolute entropy and specific heat of formaldehyde, whose computed values differ by 0.04 and 0.01 cal (mol K)⁻¹, respectively, from the experimental values reported by Gurvich and co-workers.¹¹⁰

The fundamental transition energies (which rule the SPT vibrational partition functions) have been analyzed for two other prototypical molecules, namely, formic acid and methanol [see Figs. 2(a) and 2(b)].

FIG. 2.

View large Download slide

Molecular structures of (a) formic acid, (b) methanol, and (c) acetaldehyde.

In order to put the new approach in a broader perspective, the results issued from the HDCPT2 approach employing curvilinear coordinates are compared not only with their HDCPT2-C and experimental counterparts but also with the values obtained by using the GVPT2-R method, which takes both FRs and Darling–Dennison resonances (DDRs) into full account.¹¹¹

Let us start by considering formic acid, whose fundamental wave numbers are collected in Table IV.

TABLE IV.

Comparison of the experimental and theoretical fundamental wave numbers (in cm⁻¹) of formic acid at the rDSD/3F12⁻ level of theory.

			Curvilinear coords.	Rectilinear coords.
Assignment	Symm.	Harm.	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		7443	7343	7343	7343	⋯
OH stretch.	A′	3757	3568	3568	3570	3570
CH stretch.		3093	2931	2930	2942	2943
C=O stretch.		1808	1773	1774	1770	1770
CH bend.		1412	1386	1385	1383	1387
OH bend.		1313	1253	1251	1230	1229
C–O stretch.		1135	1102	1102	1102	1105
OCO bend.		632	624	624	624	625
CH bend.	A″	1059	1037	1037	1033	1033
Torsion		676	652	652	652	638
MAE		65	7	7	2	⋯

			Curvilinear coords.	Rectilinear coords.
Assignment	Symm.	Harm.	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		7443	7343	7343	7343	⋯
OH stretch.	A′	3757	3568	3568	3570	3570
CH stretch.		3093	2931	2930	2942	2943
C=O stretch.		1808	1773	1774	1770	1770
CH bend.		1412	1386	1385	1383	1387
OH bend.		1313	1253	1251	1230	1229
C–O stretch.		1135	1102	1102	1102	1105
OCO bend.		632	624	624	624	625
CH bend.	A″	1059	1037	1037	1033	1033
Torsion		676	652	652	652	638
MAE		65	7	7	2	⋯

^a

Anharmonic calculation based on the GVPT2 scheme accounting for both FRs and DDRs effects.

^b

References 112 and 113.

All the methods tested strongly reduce the mean absolute deviation (MAE) of the harmonic approximation (from 65 to less than 10 cm⁻¹). This result underscores the non-negligible role of anharmonic contributions for a quantitative comparison with the experiment. While HDCPT2-C and HDCPT2-R results are very close, a non-negligible difference is observed with the GVPT2-R approach for CH stretching (11 cm⁻¹) and OH bending (23 cm⁻¹). These discrepancies can be traced back to the approximate treatment of FRs by using the HDCPT2 model, as a detailed analysis indicates that both vibrations are only marginally affected by DDRs. Although this aspect can be relevant for high-resolution spectroscopy, an average error below 10 cm⁻¹ on fundamental transition energies is more than acceptable for the computation of thermochemical properties.

A similar trend is shown by methanol, whose fundamental wave numbers are collected in Table V.

TABLE V.

Comparison of the experimental and theoretical fundamental wave numbers (in cm⁻¹) of methanol at the rDSD/3F12⁻ level of theory.

			Curvilinear coords.	Rectilinear coords.
Ass.	Symm.	Harm.	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		11 324	11 152	11 152	11 152	⋯
OH stretch.	A′	3865	3683	3683	3684	3681
CH₃ d-stretch.		3145	3007	3007	3016	3000
CH₃ s-stretch.		3022	2888	2888	2849	2844
CH₃ d-deform.		1527	1490	1490	1484	1477
CH₃ s-deform.		1491	1447	1445	1457	1455
OH bend.		1381	1332	1331	1325	1345
CH₃ rock.		1089	1071	1075	1074	1060
CO stretch.		1059	1033	1034	1032	1033
CH₃ d-stretch.	A″	3081	2953	2953	2974	2960
CH₃ d-deform.		1517	1484	1483	1476	1477
CH₃ rock.		1183	1156	1156	1156	1165
Torsion		291	232	232	232	200
MAE		67	9	9	7	⋯

			Curvilinear coords.	Rectilinear coords.
Ass.	Symm.	Harm.	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		11 324	11 152	11 152	11 152	⋯
OH stretch.	A′	3865	3683	3683	3684	3681
CH₃ d-stretch.		3145	3007	3007	3016	3000
CH₃ s-stretch.		3022	2888	2888	2849	2844
CH₃ d-deform.		1527	1490	1490	1484	1477
CH₃ s-deform.		1491	1447	1445	1457	1455
OH bend.		1381	1332	1331	1325	1345
CH₃ rock.		1089	1071	1075	1074	1060
CO stretch.		1059	1033	1034	1032	1033
CH₃ d-stretch.	A″	3081	2953	2953	2974	2960
CH₃ d-deform.		1517	1484	1483	1476	1477
CH₃ rock.		1183	1156	1156	1156	1165
Torsion		291	232	232	232	200
MAE		67	9	9	7	⋯

^a

Anharmonic calculation based on the GVPT2 scheme accounting for both FRs and DDRs effects.

^b

Reference 114.

Once again, the inclusion of anharmonic corrections results in a remarkable reduction of the error with respect to experiment. Although the MAEs of the different anharmonic calculations are comparable, the most significant deviations of the HDCPT2 results from their GVPT2-R counterparts are observed in the case of CH₃ symmetric stretching (39 cm⁻¹), CH₃ degenerate stretching (21 cm⁻¹), and CH₃ degenerate deformation (12 cm⁻¹). The discrepancies are always due to the different treatment of FRs, while the impact of DDRs is rather limited.

Then, the analysis has been extended to the calculation of entropies. In this case, we have also considered acetaldehyde [see Fig. 2(c)] that is characterized by a low-frequency rotation of the methyl group. In this context, the simulations have been performed using a curvilinear framework for the description of internal nuclear motions. In particular, the results obtained through harmonic, VPT2-C, and HDCPT2-C computations have been compared to their experimental counterparts to evaluate the effect of FRs on the thermochemical properties of interest (see Table VI).

TABLE VI.

Comparison between computed and experimental entropies [in cal (mol K)⁻¹] of selected systems at 298.15 and 1 atm. The rDSD/3F12⁻ level of calculation has been systematically adopted.

Molecule	Harm.a	VPT2b	HDCPT2c	Expt.d
Formic acid	59.33	59.41	59.48	59.48
Methanol	57.01	57.33	57.33	57.29
Acetaldehydee	62.66	62.67	62.67	63.06
		(63.10)	(63.10)

Molecule	Harm.a	VPT2b	HDCPT2c	Expt.d
Formic acid	59.33	59.41	59.48	59.48
Methanol	57.01	57.33	57.33	57.29
Acetaldehydee	62.66	62.67	62.67	63.06
		(63.10)	(63.10)

^a

Harmonic value of the thermochemical property, with the vibrational contribution being evaluated within the Wilson GF method.

^b

Curvilinear-based formulation of VPT2.

^c

Curvilinear-based formulation of HDCPT2.

^d

The experimental data have been reduced by 0.03 cal (mol K)⁻¹ to account for the pressure change from 1 bar = 0.1 MPa to 1 atm = 0.101 325 MPa.

^e

The values in parenthesis are obtained using the HR model for methyl torsion.

Among the investigated systems, formic acid is characterized by a difference between VPT2-C and HDCPT2-C schemes comparable with the overall anharmonic correction (0.07 vs 0.08 cal (mol K)⁻¹), with the HDCPT2-C result in excellent agreement with the experiment. Otherwise, the two perturbative variants provide the same anharmonic correction [0.32 cal (mol K)⁻¹)] in the case of methanol, leading once again to a final value very close to its experimental counterpart. Finally, as expected, the methyl torsion of acetaldehyde is poorly described by any perturbative treatment. In fact, the computed entropy of this system differs from the experimental value by about 0.39 cal (mol K)⁻¹. A viable strategy for addressing hindered rotations (HRs) leverages robust interpolation equations between harmonic oscillator and free-rotor models.^33,38 When this approach is used for hindered rotations, all other modes can be treated by the different VPT2 variants. The automatic procedure implemented in the Gaussian16 package detects without problems the torsional mode of acetaldehyde as a hindered rotation, also providing the corresponding fundamental transition energy and contribution to ZPVE. This integrated approach gives a standard entropy of acetaldehyde [63.10 cal (mol K)⁻¹] remarkably close to the experimental value [63.06 cal (mol K)⁻¹].

Now, a different problem is tackled, namely, the presence of soft degrees of freedom different from hindered rotations, which, as is well known, are poorly described by low-order perturbative methods. In such circumstances, an integrated strategy combining perturbative and variational methods can be effectively used, provided that negligible couplings are observed between the degrees of freedom treated at different levels. The approach is illustrated by the case of formamide (see Fig. 3), whose lowest frequency mode, namely, the out-of-plane motion of the NH₂ group, has a very large quartic force constant.

FIG. 3.

View large Download slide

Molecular structure of formamide.

To understand the extent of the problem, the wave numbers of the fundamental transitions have been evaluated with the variants VPT2-C or VPT2-R and compared to the corresponding experimental values. The different results are reported in Table VII.

TABLE VII.

Comparison of the experimental and theoretical fundamental wave numbers (in cm⁻¹) of formamide at the rDSD/3F12⁻ level of theory.

			Curvilinear coords.		Rectilinear coords.
Assignment	Symm.	Harm.	VPT2	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		9986	9936	9936	9936	9936	⋯
Asym. NH₂ stretch.	A′	3756	3554	3554	3554	3554	3564
Sym. NH₂ stretch.		3612	3436	3436	3435	3435	3440
CH stretch.		2991	2840	2840	2838	2866	2854
CO stretch.		1789	1756	1756	1756	1756	1754
In plane NH₂ scis.		1631	1680	1607	1606	1585	1579
In plane CH scis.		1424	1395	1393	1393	1396	1391
CN str.		1280	1487	1267	1267	1256	1258
i-NHO/NH₂ bend.		1059	1054	1054	1055	1056	1046
o-NHO/NH₂ bend.		568	567	567	569	569	566
Out of pl. CH def.	A″	1048	1021	1021	1023	1024	1033
Torsion		638	587	570	570	590	602
Out of pl. NH₂		175	576	576	576	568	289
			(342)c	(342)c	(608)c	(608)c
MAEd		64	36	11	7	7
MAEe		68	57	34	29	29
MAEf		(29)	(14)	(14)	(14)	(57)

			Curvilinear coords.		Rectilinear coords.
Assignment	Symm.	Harm.	VPT2	HDCPT2	HDCPT2	GVPT2a	Expt.b
ZPVE		9986	9936	9936	9936	9936	⋯
Asym. NH₂ stretch.	A′	3756	3554	3554	3554	3554	3564
Sym. NH₂ stretch.		3612	3436	3436	3435	3435	3440
CH stretch.		2991	2840	2840	2838	2866	2854
CO stretch.		1789	1756	1756	1756	1756	1754
In plane NH₂ scis.		1631	1680	1607	1606	1585	1579
In plane CH scis.		1424	1395	1393	1393	1396	1391
CN str.		1280	1487	1267	1267	1256	1258
i-NHO/NH₂ bend.		1059	1054	1054	1055	1056	1046
o-NHO/NH₂ bend.		568	567	567	569	569	566
Out of pl. CH def.	A″	1048	1021	1021	1023	1024	1033
Torsion		638	587	570	570	590	602
Out of pl. NH₂		175	576	576	576	568	289
			(342)c	(342)c	(608)c	(608)c
MAEd		64	36	11	7	7
MAEe		68	57	34	29	29
MAEf		(29)	(14)	(14)	(14)	(57)

^a

Anharmonic calculation based on the GVPT2 scheme accounting for both FRs and DDRs effects.

^b

Reference 115.

^c

Fundamental wave number determined through a one-dimensional VCI calculation along the normal mode associated with the out-of-plane motion of the NH₂ group.

^d

Mean absolute error of the theoretical wave numbers, with the fundamental transition of the out-of-plane motion of the NH₂ group being neglected.

^e

Mean absolute error of the theoretical wave numbers, with the fundamental transition of the out-of-plane motion of the NH₂ group being computed through perturbation theory.

^f

Mean absolute error of the theoretical wave numbers, with the fundamental transition of the out-of-plane motion of the NH₂ group being computed through a one-dimensional VCI.

First, the MAE was computed for the different sets of fundamentals without including the LAM. As expected, the corresponding values are in line with those obtained for formic acid and methanol, with the quite large MAE obtained at the harmonic level (68 cm⁻¹) being strongly reduced (to 11 and 7 cm⁻¹) when anharmonic corrections are taken into account. Inclusion of the NH₂ inversion in the error analysis leads to significant increases in MAEs related to the huge anharmonic correction associated with this LAM. In order to bypass this problem, a stepwise procedure has been adopted. First, the third- and fourth-order vibrational couplings between this vibration and the other ones have been grouped in different wave-number ranges. The corresponding histograms issued from rectilinear and curvilinear coordinates are reported in Fig. 4.

FIG. 4.

Graphical representation of the distribution of cubic (left panel) and quartic (right panel) couplings for the out-of-plane motion of the NH2 group in formamide. The results are presented as functions of the wave number range, based on calculations at the rDSD/3F12− level of theory.

View large Download slide

Graphical representation of the distribution of cubic (left panel) and quartic (right panel) couplings for the out-of-plane motion of the NH₂ group in formamide. The results are presented as functions of the wave number range, based on calculations at the rDSD/3F12⁻ level of theory.

As expected, only the rectilinear anharmonic force constants show values larger than 1000 cm⁻¹ and, in general, at wave numbers higher than those of their curvilinear counterparts. This confirms that the LAM is much less coupled to other modes when employing curvilinear internal coordinates. Next, a VPT2 calculation is carried out in the space of SAMs, retaining the couplings with the LAM, but removing the diagonal anharmonic force constants of this mode, as well as the corresponding derivatives of the g matrix. Note that, in general, this approximation also affects the transition frequencies of the SAMs, due to the neglect of the contribution given by Eq. (41), with F being the LAM. However, in the present case, this term vanishes since the diagonal cubic force constant, and the first-order g derivative along the LAM are null by symmetry, so that the fundamental absorption bands of the SAMs do not change when going from the full perturbative treatment to this reduced-dimensionality variant.

Finally, one-dimensional calculations involving the LAM have been carried out. In previous work, relaxed one-dimensional scans along selected coordinates were performed and the corresponding one-dimensional vibrational Hamiltonian were treated by the discrete variable representation (DVR).^62,116,117 Here, we employ a different approach, which does not require any additional computation with respect to those needed by the VPT2 approach. One-dimensional vibrational configuration interaction (VCI) computations have been performed with the same basis functions and anharmonic terms entering the VPT2 treatment in conjunction with force fields expressed in terms of either rectilinear or curvilinear coordinates (and g matrix in the latter case). The results of those computations (summarized in Table VII) show that the use of rectilinear coordinates leads to a fundamental transition energy (608 cm⁻¹) even worse than the corresponding perturbative value (568 cm⁻¹), with this behavior being related to the too large diagonal quartic force constant (more than 30 000 cm⁻¹). In fact, finite displacements along the rectilinear approximation of the out-of-plane distortion involve necessarily non-negligible contributions from NH stretchings, with the consequent large (spurious) energy increase. Otherwise, the out-of-plane motion is well-described by the corresponding normal mode expressed in terms of curvilinear coordinates, and the one-dimensional VCI treatment leads to a value of 342 cm⁻¹, lowering the error with respect to the experimental data from 287 cm⁻¹ to 53 cm⁻¹. For the sake of completeness, the convergence of the VCI procedure with respect to the size of the harmonic basis set has been analyzed. The results sketched in Fig. 5 evidence the robustness of the implemented protocol since convergence within 0.1 cm⁻¹ is achieved with a basis set of 19 Hermite polynomials, and convergence within 0.0001 cm⁻¹ is reached with 50 basis functions.

FIG. 5.

VCI wave number associated with the out-of-plane motion of the NH2 group of formamide at the rDSD/3F12− level of theory. The results are presented as functions of the basis set size.

View large Download slide

VCI wave number associated with the out-of-plane motion of the NH₂ group of formamide at the rDSD/3F12⁻ level of theory. The results are presented as functions of the basis set size.

In summary, the use of curvilinear internal coordinates paves the way toward the study of flexible molecules by means of a VPT2 treatment of SAMs in conjunction with effective hindered rotor and VCI one-dimensional treatments for periodic and non-periodic LAMs, respectively. Furthermore a recent benchmark study for 21 closed-shell and 10 open-shell small molecules for which accurate experimental data are available showed that a comparable procedure based on rectilinear coordinates and rDSD/3F12⁻ force fields achieved mean absolute errors of 0.05 kcal mol⁻¹ for ZPVEs and 0.05 cal (mol K)⁻¹ for entropies.³⁹ Therefore, the rDSD/3F12⁻ model will be our reference for larger systems.

B. Toward large molecules

The analysis presented in this section has been limited to relatively small molecular systems. When the focus shifts to large systems, an additional issue comes into play, namely, the problem of determining accurate vibrational properties while retaining an affordable computational cost. An effective solution is offered by the so-called dual-level method^1,20 in which the harmonic terms and the corresponding anharmonic corrections are evaluated at different levels of theory, since the accuracy of the former (generally much larger) contributions iscule usually more critical. Furthermore, the standard finite-difference approach employed for the construction of the quartic force field requires the computation of 2N′ + 1 analytical force constant matrices, whose cost increases steeply with molecular size, becoming rapidly prohibitive even for double hybrid functionals. Within the dual-level framework, the value of a generic vibrational property P is

P = P_{HL}^{H} + Δ P_{LL}^{A},

(73)

where P^H is the harmonic term, ΔP^A is the anharmonic correction, while HL and LL refer to generic high- and low-level methods, respectively. In order to clearly understand the impact of the harmonic contribution, it is convenient to rewrite Eq. in the following way:

P = P_{LL}^{H} + Δ P^{H} + Δ P_{LL}^{A},

(74)

where

Δ P^{H} = P_{HL}^{H} - P_{LL}^{H}

is the corrective factor that leads from the LL value of the harmonic term to its HL counterpart. In the present work, rDSD/3F12⁻ and B3/SVP (or alternatively, B3/6-31G*) have been used systematically as the HL and LL methods, respectively.

The reliability of this scheme has been tested on two systems involving either hindered rotations (pyruvic acid) or crowded moieties (norbornane) for which accurate computational results and/or experimental data are available.

The Tc conformer of pyruvic acid [see Fig. 6(a)] has been analyzed by different anharmonic models, including dual-level schemes. In particular, full calculations at the B3/SVP and rDSD/3F12⁻ levels have been performed. Furthermore, the anharmonic corrections at the B3/SVP level have been combined to harmonic frequencies at the rDSD/3F12⁻ level and with those at the CCSD(T)/cc-pVTZ level reported in Ref. 118. In addition, in this case, the largest errors characterize the purely harmonic calculations and a substantial improvement (fourfold reduction of the MAE) is produced by anharmonic calculation at the B3/SVP level. All other anharmonic calculations have comparable MAEs (5–6 cm⁻¹) corresponding to further twofold reduction. The most interesting aspect is that the accuracy of either the full rDSD/3F12⁻ calculation, or the corresponding one based on B3/SVP anharmonic corrections rivals that obtained through CCSD(T) harmonic frequencies. As already mentioned, the cost of this last method becomes rapidly prohibitive as molecular size increases, making dual-level approaches entirely based on DFT particularly appealing. Let us stress that this remark is further confirmed by the comparison of the anharmonic wave numbers at the rDSD/3F12⁻ level with those in which the anharmonic corrections are replaced by the corresponding B3/SVP values. From this perspective, only the harmonic frequencies require an enhanced level of theory, implying a significant reduction of the overall computational cost.

FIG. 6.

View large Download slide

Molecular structures of (a) Tc conformer of pyruvic acid and (b) norbornane.

Shaw et al.¹²¹ carried out a detailed investigation of norbornane [see Fig. 6(b)] through IR spectroscopy, reporting the experimental fundamentals within the spectral range from 600 to 1600 cm⁻¹. In this context, the computation of the harmonic wave numbers has been performed both at the B3/6-31G* and rDSD/3F12⁻ levels of theory, while the anharmonic corrections have been computed at the B3/6-31G* level in the framework of the HDCPT2 scheme. The results obtained using Eq. (74) are reported in Table VIII and compared to the available experimental data. The MAE obtained through harmonic calculations at the B3/6-31G* level (31 cm⁻¹) decreases to 7 cm⁻¹ after the refinement of the harmonic frequencies and the introduction of the anharmonic corrections, thus leading to remarkable agreement with the experimental data. It should be noted that ten different states exhibit a harmonic correction with absolute value $\geq 10$ cm⁻¹. As a matter of fact, a non-negligible reduction of the MAE is attributable to the correction of the harmonic terms, underscoring the substantial impact of this refinement.

TABLE VIII.

Comparison of the experimental and theoretical fundamental wave numbers of norbornane (in cm⁻¹) computed at the rDSD//B3 level of theory, through the HDCPT2 scheme.

State	Symm.	Harm.a	ΔHarm.b	ΔAnharm.c	Anharm.d	Expt.e
ZPVE		39 092	−98	−1321	37 673	⋯
\|1₁⟩	A₁	3113	2	−148	2967	⋯
\|1₂⟩		3108	−1	−149	2958	⋯
\|1₃⟩		3066	−4	−108	2954	⋯
\|1₄⟩		3057	−7	−135	2915	⋯
\|1₅⟩		1556	−19	−42	1495	1487
\|1₆⟩		1524	−20	−45	1459	1455
\|1₇⟩		1365	−11	−35	1319	1310
\|1₈⟩		1302	−6	−31	1265	1253
\|1₉⟩		1178	−6	−28	1144	1140
\|1₁₀⟩		1014	7	−23	998	989
\|1₁₁⟩		935	9	−19	925	916
\|1₁₂⟩		890	3	−19	874	864
\|1₁₃⟩		828	2	−16	814	818
\|1₁₄⟩		765	4	−15	754	758
\|1₁₅⟩		406	−3	−9	394	401
\|1₁₆⟩	A₂	3086	2	−149	2939	⋯
\|1₁₇⟩		3054	−4	−115	2935	⋯
\|1₁₈⟩		1513	−19	−42	1452	1451
\|1₁₉⟩		1345	−3	−37	1305	1294
\|1₂₀⟩		1321	−9	−33	1279	1290
\|1₂₁⟩		1253	−6	−33	1214	1212
\|1₂₂⟩		1151	−5	−20	1126	1111
\|1₂₃⟩		975	7	−19	963	956
\|1₂₄⟩		967	−4	−27	936	949
\|1₂₅⟩		548	−3	−11	534	546
\|1₂₆⟩		158	−13	−8	137	152
\|1₂₇⟩	B₁	3108	2	−149	2961	⋯
\|1₂₈⟩		3102	1	−147	2956	⋯
\|1₂₉⟩		3064	−2	−125	2937	⋯
\|1₃₀⟩		1532	−18	−58	1456	1466
\|1₃₁⟩		1368	−11	−33	1324	1315
\|1₃₂⟩		1256	−5	−39	1212	1211
\|1₃₃⟩		1141	−10	−16	1115	1103
\|1₃₄⟩		1050	1	−27	1024	1027
\|1₃₅⟩		972	5	−23	954	957
\|1₃₆⟩		905	9	−21	893	877
\|1₃₇⟩		811	−7	−19	785	786
\|1₃₈⟩		453	−6	2	449	435
\|1₃₉⟩	B₂	3103	−2	−137	2964	⋯
\|1₄₀⟩		3090	1	−144	2947	⋯
\|1₄₁⟩		3054	−4	−115	2935	⋯
\|1₄₂⟩		1527	−20	−51	1456	1462
\|1₄₃⟩		1365	−13	−34	1318	1320
\|1₄₄⟩		1306	−6	−32	1268	1272
\|1₄₅⟩		1288	−7	−34	1247	1244
\|1₄₆⟩		1199	−3	−31	1165	1177
\|1₄₇⟩		1104	−6	−22	1076	1081
\|1₄₈⟩		971	4	−17	958	956
\|1₄₉⟩		829	8	−18	819	814
\|1₅₀⟩		769	5	−19	755	760
\|1₅₁⟩		337	−5	−16	316	336
MAEf		31	⋯	⋯	7	⋯

State	Symm.	Harm.a	ΔHarm.b	ΔAnharm.c	Anharm.d	Expt.e
ZPVE		39 092	−98	−1321	37 673	⋯
\|1₁⟩	A₁	3113	2	−148	2967	⋯
\|1₂⟩		3108	−1	−149	2958	⋯
\|1₃⟩		3066	−4	−108	2954	⋯
\|1₄⟩		3057	−7	−135	2915	⋯
\|1₅⟩		1556	−19	−42	1495	1487
\|1₆⟩		1524	−20	−45	1459	1455
\|1₇⟩		1365	−11	−35	1319	1310
\|1₈⟩		1302	−6	−31	1265	1253
\|1₉⟩		1178	−6	−28	1144	1140
\|1₁₀⟩		1014	7	−23	998	989
\|1₁₁⟩		935	9	−19	925	916
\|1₁₂⟩		890	3	−19	874	864
\|1₁₃⟩		828	2	−16	814	818
\|1₁₄⟩		765	4	−15	754	758
\|1₁₅⟩		406	−3	−9	394	401
\|1₁₆⟩	A₂	3086	2	−149	2939	⋯
\|1₁₇⟩		3054	−4	−115	2935	⋯
\|1₁₈⟩		1513	−19	−42	1452	1451
\|1₁₉⟩		1345	−3	−37	1305	1294
\|1₂₀⟩		1321	−9	−33	1279	1290
\|1₂₁⟩		1253	−6	−33	1214	1212
\|1₂₂⟩		1151	−5	−20	1126	1111
\|1₂₃⟩		975	7	−19	963	956
\|1₂₄⟩		967	−4	−27	936	949
\|1₂₅⟩		548	−3	−11	534	546
\|1₂₆⟩		158	−13	−8	137	152
\|1₂₇⟩	B₁	3108	2	−149	2961	⋯
\|1₂₈⟩		3102	1	−147	2956	⋯
\|1₂₉⟩		3064	−2	−125	2937	⋯
\|1₃₀⟩		1532	−18	−58	1456	1466
\|1₃₁⟩		1368	−11	−33	1324	1315
\|1₃₂⟩		1256	−5	−39	1212	1211
\|1₃₃⟩		1141	−10	−16	1115	1103
\|1₃₄⟩		1050	1	−27	1024	1027
\|1₃₅⟩		972	5	−23	954	957
\|1₃₆⟩		905	9	−21	893	877
\|1₃₇⟩		811	−7	−19	785	786
\|1₃₈⟩		453	−6	2	449	435
\|1₃₉⟩	B₂	3103	−2	−137	2964	⋯
\|1₄₀⟩		3090	1	−144	2947	⋯
\|1₄₁⟩		3054	−4	−115	2935	⋯
\|1₄₂⟩		1527	−20	−51	1456	1462
\|1₄₃⟩		1365	−13	−34	1318	1320
\|1₄₄⟩		1306	−6	−32	1268	1272
\|1₄₅⟩		1288	−7	−34	1247	1244
\|1₄₆⟩		1199	−3	−31	1165	1177
\|1₄₇⟩		1104	−6	−22	1076	1081
\|1₄₈⟩		971	4	−17	958	956
\|1₄₉⟩		829	8	−18	819	814
\|1₅₀⟩		769	5	−19	755	760
\|1₅₁⟩		337	−5	−16	316	336
MAEf		31	⋯	⋯	7	⋯

^a

B3/6-31G* level of theory.

^b

Difference between the harmonic contributions at the rDSD/3F12⁻ level and those at the B3/6-31G* level.

^c

Anharmonic corrections at the B3/6-31G* level.

^d

Fully anharmonic results.

^e

Reference 121.

^f

Mean absolute error of the fundamental wave numbers with respect to their experimental counterparts, considering only those modes for which experimental values are available.

This kind of computation can be routinely performed for molecules containing up to about 25 atoms. For larger systems, anharmonic force fields can still be obtained at the B3 level, but harmonic computations at the rDSD level become prohibitively expensive in terms of both time and memory requirements. Therefore, one has to rely on B3 results for both harmonic and anharmonic contributions. The errors issued by this approximation have been tested by comparison of rDSD/3F12⁻ and B3/SVP harmonic zero point energies for norbornane together with the panel of molecules shown in Fig. 7.

FIG. 7.

View large Download slide

Molecular structures of (a) benzoquinone, (b) naphthoquinone, (c) (2CN-)naphthalene, (d) (9CN-)anthracene, (e) (9CN-)phenanthrene, and (f) (2CN-)pyrene. X = H, CN.

The results collected in Table X highlight once again not only the importance of an anharmonic description of vibrational motions but also the significant role played by an accurate description of the underlying harmonic force field. In particular, B3 harmonic frequencies are systematically larger than their rDSD counterparts and a scaling factor of 0.996 reduces the average error of ZPVE’s around 0.1 kcal mol⁻¹ for representative hydrocarbons and quinones. The OH stretching of pyruvic acid (whose frequency is underestimated by about 100 cm⁻¹ as shown in Table IX) is paradigmatic of pathological situations (including also CN stretchings in cyano compounds¹²² or some X–H stretching modes in the presence of hydrogen bridges¹²³) affected by much larger errors. Here, internal coordinates offer significant advantages as a limited number of energies (and/or gradients) along predefined directions is sufficient to determine the required corrections. Although the full development of this approach goes beyond the scope of the present paper, its potentialities are illustrated in Table X for the specific case of cyano derivatives of polyaromatic hydrocarbons (CN-PAHs), which are of increasing importance in the field of astrochemistry.^124,125 As a matter of fact, a constant correction of 0.3 kcal mol⁻¹ (1.3 times the nearly constant difference of 77 cm⁻¹ between the harmonic frequencies of CN stretchings computed at the B3 and rDSD level) for each CN group produces errors close to 0.1 kcal mol⁻¹ also for those molecules. Note that an analogous procedure (scaling by 0.996 for all the modes plus 1.3 times the difference of the OH stretching frequency) also performs a very good job for the ZPVE of pyruvic acid [the harmonic B3 value is increased from 46.48 to 46.68 kcal mol⁻¹ upon scaling, to be compared to the identical rDSD and CCSD(T) values of 46.69 kcal mol⁻¹]. Otherwise, the effect of these corrections is negligible on partition functions (hence, entropies and thermal contributions to enthalpies).

TABLE IX.

Comparison of the experimental and theoretical fundamental wave numbers (in cm⁻¹) of the Tc conformer of pyruvic acid.

		ωa		δνb			νc
Ass.	Symm.	B3d	rDSDe	CCf	B3d	B3d	rDSD//B3g	CC//B3h	rDSDe	exp.i
ZPVE		15 674	15 743	15 743	−191	15 483	15 552	15 552	15 553	⋯
OH str.	A′	3569	3665	3667	−205	3364	3460	3462	3465	3463
CH₃ a-str.		3176	3181	3174	−155	3021	3026	3019	3036	3025
CH₃ s-str.		3061	3062	3056	−117	2944	2945	2939	2948	2941
CO str.		1849	1835	1853	−38	1811	1797	1815	1801	1804
CO str.		1788	1768	1773	−28	1760	1740	1745	1733	1737
CH₃ a-bend.		1477	1468	1464	−34	1443	1434	1430	1433	1424
CC a-str.		1421	1423	1430	−35	1386	1388	1395	1393	1391
CH₃ s-bend.		1391	1394	1394	−34	1357	1360	1360	1344	1360
COH bend.		1248	1263	1271	−42	1206	1221	1229	1221	1211
CO str.		1163	1167	1166	−28	1135	1139	1138	1138	1133
CH₃ rock.		991	988	984	−15	976	973	969	972	970
CC s-str.		766	777	775	−17	749	760	758	761	761
CO bend.		606	611	608	−7	599	604	601	606	604
CO bend.		525	530	532	−1	524	529	531	531	⋯
CCO bend.		394	394	395	−6	388	388	389	393	389
CCC bend.		254	253	253	−1	253	252	252	267	258
CH₃ a-str.	A″	3118	3131	3126	−153	2965	2978	2973	2983	⋯
CH₃ a-bend.		1478	1474	1470	−40	1438	1434	1430	1431	⋯
CH₃ rock.		1040	1044	1040	−18	1022	1026	1022	1022	1030
CO rock.		741	741	741	−26	715	715	715	724	⋯
OH tors.		687	701	701	−26	661	675	675	681	668
CO rock.		395	396	394	−3	392	393	391	396	394
CH₃ tors.		119	126	124	21j	140	147	145	144k	134
CC tors.		90	93	94	5j	95	98	99	95k	90
MAEl		37	42	43	⋯	11	5	6	5	⋯

		ωa		δνb			νc
Ass.	Symm.	B3d	rDSDe	CCf	B3d	B3d	rDSD//B3g	CC//B3h	rDSDe	exp.i
ZPVE		15 674	15 743	15 743	−191	15 483	15 552	15 552	15 553	⋯
OH str.	A′	3569	3665	3667	−205	3364	3460	3462	3465	3463
CH₃ a-str.		3176	3181	3174	−155	3021	3026	3019	3036	3025
CH₃ s-str.		3061	3062	3056	−117	2944	2945	2939	2948	2941
CO str.		1849	1835	1853	−38	1811	1797	1815	1801	1804
CO str.		1788	1768	1773	−28	1760	1740	1745	1733	1737
CH₃ a-bend.		1477	1468	1464	−34	1443	1434	1430	1433	1424
CC a-str.		1421	1423	1430	−35	1386	1388	1395	1393	1391
CH₃ s-bend.		1391	1394	1394	−34	1357	1360	1360	1344	1360
COH bend.		1248	1263	1271	−42	1206	1221	1229	1221	1211
CO str.		1163	1167	1166	−28	1135	1139	1138	1138	1133
CH₃ rock.		991	988	984	−15	976	973	969	972	970
CC s-str.		766	777	775	−17	749	760	758	761	761
CO bend.		606	611	608	−7	599	604	601	606	604
CO bend.		525	530	532	−1	524	529	531	531	⋯
CCO bend.		394	394	395	−6	388	388	389	393	389
CCC bend.		254	253	253	−1	253	252	252	267	258
CH₃ a-str.	A″	3118	3131	3126	−153	2965	2978	2973	2983	⋯
CH₃ a-bend.		1478	1474	1470	−40	1438	1434	1430	1431	⋯
CH₃ rock.		1040	1044	1040	−18	1022	1026	1022	1022	1030
CO rock.		741	741	741	−26	715	715	715	724	⋯
OH tors.		687	701	701	−26	661	675	675	681	668
CO rock.		395	396	394	−3	392	393	391	396	394
CH₃ tors.		119	126	124	21j	140	147	145	144k	134
CC tors.		90	93	94	5j	95	98	99	95k	90
MAEl		37	42	43	⋯	11	5	6	5	⋯

^a

Harmonic result.

^b

Anharmonic correction.

^c

Full anharmonic result.

^d

B3/SVP level of theory.

^e

rDSD/3F12⁻ level of theory.

^f

CCSD(T)/cc-pVTZ level of theory from Ref. 118.

^g

Anharmonic fundamental wave numbers obtained by combining rDSD/3F12⁻ harmonic frequencies with B3/SVP anharmonic corrections.

^h

Anharmonic fundamental wave numbers obtained by combining CCSD(T)/cc-pVTZ harmonic frequencies with B3/SVP anharmonic corrections.

ⁱ

References 119 and 120.

^j

Correction computed through the HR model.

^k

Fundamental wave number computed through the HR model.

^l

Mean absolute error of the computed anharmonic wave numbers with respect to their experimental counterparts.

TABLE X.

Comparison of harmonic and anharmonic ZPVE (in kcal mol⁻¹) for medium-sized molecules.

Molecule	Harm.			ΔAnharm
Molecule	B3a	rDSDb	B3 scal.c	B3a
Norbornane	115.94	115.49	115.48	−3.92
Benzoquinone	53.41	53.26	53.20	−0.50
Naphtoquinone	83.05	83.31	82.98	−0.81
Naphtalene	92.67	92.27	92.30	−1.02
Anthracene	121.95	121.31	121.46	−1.33
Phenanthrene	122.24	121.68	121.75	−1.29
Pyrene	130.21	129.56	129.69	−1.14
2CN-naphtalene	91.97	91.23	91.30d	−0.89
9CN-anthracene	121.16	120.39	120.38d	−1.24
9CN-phenanthrene	121.39	120.61	120.60d	−1.21
2CN-pyrene	129.27	128.48	128.45d	−1.10

Molecule	Harm.			ΔAnharm
Molecule	B3a	rDSDb	B3 scal.c	B3a
Norbornane	115.94	115.49	115.48	−3.92
Benzoquinone	53.41	53.26	53.20	−0.50
Naphtoquinone	83.05	83.31	82.98	−0.81
Naphtalene	92.67	92.27	92.30	−1.02
Anthracene	121.95	121.31	121.46	−1.33
Phenanthrene	122.24	121.68	121.75	−1.29
Pyrene	130.21	129.56	129.69	−1.14
2CN-naphtalene	91.97	91.23	91.30d	−0.89
9CN-anthracene	121.16	120.39	120.38d	−1.24
9CN-phenanthrene	121.39	120.61	120.60d	−1.21
2CN-pyrene	129.27	128.48	128.45d	−1.10

^a

B3/SVP level of theory.

^b

rDSD/3F12⁻ level of theory.

^c

B3/SVP value scaled by a constant factor of 0.996.

^d

Reduction of 0.3 kcal mol⁻¹ for the presence of a cyano group.

The results obtained through this approach for molecules containing up to about 50 atoms are shown in Table XI for the two hormones, the molecular motor, and the nucleoside sketched in Fig. 8.

TABLE XI.

Comparison of the harmonic and anharmonic values of the ZPVE (in kcal mol⁻¹), enthalpy difference between standard conditions and 0 K (ΔH in kcal mol⁻¹), entropy [S in cal (mol K)⁻¹], and specific heat at constant pressure [c^p in cal (mol K)⁻¹] at 298 K and 1 atm, for different large-sized systems, evaluated within the VPT2-C framework through the HDCPT2 scheme at the B3/SVP level of theory.

	Property	Harm.a	ΔAnharm.b	Anharm.c
Androsterone	ZPVEd	291.76	−3.68	288.08
		(290.59)		(286.91)
	ΔH	12.66	0.09	12.75
	S	137.72	2.70	140.41
	c_p	82.72	1.00	83.73
β-estradiol	ZPVEd	244.25	−3.21	241.04
		(243.27)		(240.06)
	ΔH	11.53	−0.28	11.25
	S	130.79	−3.33	127.46
	c_p	74.74	−0.35	74.39
Mol. Motor	ZPVEd	242.78	−2.64	240.14
		(241.81)		(239.17)
	ΔH	12.91	0.35	13.26
	S	142.66	8.12	150.78
	c_p	85.09	1.41	86.50
Uridine	ZPVEd	142.62	−2.02	140.60
		(142.05)		(140.03)
	ΔH	10.01	−0.08	9.93
	S	123.19	−1.55	121.63
	c_p	59.94	−0.47	59.47

	Property	Harm.a	ΔAnharm.b	Anharm.c
Androsterone	ZPVEd	291.76	−3.68	288.08
		(290.59)		(286.91)
	ΔH	12.66	0.09	12.75
	S	137.72	2.70	140.41
	c_p	82.72	1.00	83.73
β-estradiol	ZPVEd	244.25	−3.21	241.04
		(243.27)		(240.06)
	ΔH	11.53	−0.28	11.25
	S	130.79	−3.33	127.46
	c_p	74.74	−0.35	74.39
Mol. Motor	ZPVEd	242.78	−2.64	240.14
		(241.81)		(239.17)
	ΔH	12.91	0.35	13.26
	S	142.66	8.12	150.78
	c_p	85.09	1.41	86.50
Uridine	ZPVEd	142.62	−2.02	140.60
		(142.05)		(140.03)
	ΔH	10.01	−0.08	9.93
	S	123.19	−1.55	121.63
	c_p	59.94	−0.47	59.47

^a

Harmonic value of the thermochemical property, with the vibrational contribution being evaluated within the Wilson GF method.

^b

Anharmonic correction to the thermochemical property.

^c

Curvilinear-based formulation of HDCPT2.

^d

In parenthesis: harmonic values scaled by 0.996.

FIG. 8.

View large Download slide

Molecular structure of (a) androsterone, (b) β-estradiol, (c) molecular motor, and (d) uridine.

The electronic energies of molecules of this size can be computed with an error within 1 kcal mol⁻¹ through last generation CCSD(T) methods using frozen natural orbitals (FNO) or pair natural orbitals (PNO) possibly paired with explicit correlation (F12) and other techniques to accelerate basis set convergence.^6,8,9 It is quite apparent that in this case, anharmonic contributions cannot be neglected and their inclusion at a reasonable cost, together with the scaling of harmonic contributions computed at the B3 level, paves the way toward systematic investigation of large molecules.

VI. CONCLUSION

This work extends and validates a recent second-order perturbative method based on curvilinear internal coordinates for the calculation of thermodynamic functions. In this way, inter-mode couplings are generally reduced with respect to rectilinear descriptions, and above all, the remaining contributions are shifted from potential to kinetic terms. This feature determines the great computational advantage of evaluating a coordinate-dependent kinetic energy operator (from analytical derivatives of internal curvilinear coordinates with respect to Cartesian ones) instead of high-order coupling potentials from expensive electronic structure calculations. Apart from that, the interpretation of mode couplings is simplified, and finally, the resulting deeper understanding of vibrational systems may yield better problem-adapted approximations.

The computational workflow starts with an accurate evaluation of anharmonic zero-point vibrational energies for semi-rigid molecular systems. The hybrid degeneracy-corrected second-order perturbation theory is then applied for the calculation of fundamental vibrational bands, enabling reliable calculations even in pathological cases and automatically handling soft degrees of freedom through accurate one-dimensional models. The simple perturbation theory is finally used to obtain reliable yet computationally inexpensive vibrational partition functions that provide all the needed vibrational contributions to thermodynamic properties. The accuracy of the proposed model has been validated by computing fundamental vibrational transitions for a panel of molecules of different sizes. The reliability of a global hybrid functional paired to a polarized split-valence basis set for anharmonic calculations is also highlighted. Additional improvements are achieved when B3 harmonic frequencies are replaced with their counterparts obtained by a double hybrid functional in conjunction with more extended basis sets. In general, this computational strategy enables a fully unsupervised evaluation of vibrational contributions to thermochemical properties with an overall accuracy close to 0.1 kcal mol⁻¹ in the case of enthalpies and 0.1 cal mol⁻¹ K⁻¹ in the case of entropies for molecules containing up to about 50 atoms.

The proposed approach can be applied alongside any quantum mechanical model for which analytical Hessians are available. This means that, in addition to isolated molecules, a black-box procedure for evaluating thermochemical properties beyond the harmonic approximation is also applicable to solution-phase systems described by the polarizable continuum model, whose analytical second derivatives are available.⁹⁶ Although further improvements are certainly required, particularly in relation to coupled large-amplitude motions, we think that the present version of the computational workflow can be confidently used for the study of medium-to-large molecular systems in both the gas phase and in solution by experiment-oriented researchers.

Let us finally remark that the methodology proposed in this work has been devised to provide an automatic, black-box approach for thermochemistry, retaining the intuitive system of harmonic oscillator quantum numbers, as well as a treatment of Fermi resonances deprived of parameters defined by the user. Although this framework is, in general, well-suited for the calculation of thermochemical data, it is not necessarily the most accurate approach for the simulation of vibrational spectra. Therefore,the perturbative treatment based on curvilinear coordinates will be shortly extended to its generalized version, explicitly accounting for both Fermi and Darling–Dennison resonances.

SUPPLEMENTARY MATERIAL

The supplementary material includes the definition of the anharmonic χ matrix within the different frameworks of VPT2.

ACKNOWLEDGMENTS

The authors acknowledge the funding from Gaussian, Inc.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

M. Mendolicchio: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Investigation (lead); Methodology (equal); Software (lead); Supervision (lead); Validation (lead); Writing – original draft (equal). V. Barone: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Validation (supporting); Writing – original draft (equal).

DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

APPENDIX: CALCULATION OF THE EXTRA-POTENTIAL TERM

Inspection of Eq. (30) shows that the extra-potential term (V_g) depends on first- and second-order derivatives of the determinant $\tilde{g}$ (or the ${\tilde{g}}_{i j \dots}$ coefficients) of the vibro-rotational matrix g. Moreover, the second-order off-diagonal derivatives of g are not required at a stationary point, as the matrix g is diagonal when the Eckart embedding is adopted $(g_{i j} = g_{i i} δ_{i j})$ ⁠. Therefore, it is possible to restrict the analysis to the coefficients ${\tilde{g}}_{i}$ and ${\tilde{g}}_{i i}$ ⁠.

One possible approach for calculating the desired terms relies entirely on the numerical differentiation of analytical determinants,

{\tilde{g}}_{i} = \frac{1}{{\tilde{g}}^{e q}} [\frac{\tilde{g} (+ δ q_{i}) - \tilde{g} (- δ q_{i})}{2 δ q_{i}}],

(A1)

{\tilde{g}}_{i i} = \frac{1}{{\tilde{g}}^{e q}} [\frac{\tilde{g} (+ δ q_{i}) + \tilde{g} (- δ q_{i}) - 2 \tilde{g} (q^{e q})}{δ q_{i}^{2}}],

(A2)

where δq_i is the displacement along the normal mode i. Alternatively, the coefficients

{\tilde{g}}_{i j \dots}

can be related to the derivatives of the matrix g. In particular,

{\tilde{g}}_{i} = \sum_{j = 1}^{N^{'} + N_{r}} \frac{g_{j j, i}}{g_{j j}},

(A3)

where N′ and N_r are the numbers of active vibrations and of rotational degrees of freedom, respectively. Similarly, the second-order terms can be obtained using the following expression:

{\tilde{g}}_{i i} = \sum_{j = 1}^{N^{'} + N_{r}} [\frac{g_{j j, i i}}{g_{j j}} + \sum_{k = 1}^{N^{'} + N_{r}} \frac{g_{j j, i} g_{k k, i} - g_{j k, i}^{2}}{g_{j j} g_{k k}}] .

(A4)

In summary, calculating V_g requires the complete set of first-order derivatives of the vibro-rotational matrix g, along with the diagonal derivatives of its vibrational and rotational blocks. As these terms are required for evaluating vibro-rotational spectroscopic parameters (see Ref. 37), V_g can be obtained without any additional computational effort.

REFERENCES

1.

C.

Puzzarini

,

J.

Bloino

,

N.

Tasinato

, and

V.

Barone

, “

Accuracy and interpretability: The devil and the holy grail. New routes across old boundaries in computational spectroscopy

,”

Chem. Rev.

119

,

8131

–

8191

(

2019

).

https://doi.org/10.1021/acs.chemrev.9b00007

Google Scholar

Crossref

PubMed

2.

J. M. L.

Martin

and

M. K.

Kesharwani

, “

Assessment of CCSD(T)-F12 approximations and basis sets for harmonic vibrational frequencies

,”

J. Chem. Theory Comput.

10

,

2085

–

2090

(

2014

).

https://doi.org/10.1021/ct500174q

Google Scholar

Crossref

3.

D.

Agbaglo

and

R. C.

Fortenberry

, “

The performance of explicitly correlated methods for the computation of anharmonic vibrational frequencies

,”

Int. J. Quantum Chem.

119

,

e25899

(

2019

).

https://doi.org/10.1002/qua.25899

Google Scholar

Crossref

4.

Z. N.

Heim

,

B. K.

Amberger

,

B. J.

Esselman

,

J. F.

Stanton

,

R. C.

Woods

, and

R. J.

McMahon

, “

Molecular structure determination: Equilibrium structure of pyrimidine (m-C₄H₄N₂) from rotational spectroscopy (⁠

(r_{e}^{SE})

⁠) and high-level ab initio calculation (r_e) agree within the uncertainty of experimental measurement

,”

J. Chem. Phys.

152

,

104303

(

2020

).

https://doi.org/10.1063/1.5144914

Google Scholar

Crossref

5.

K.

Raghavachari

,

G. W.

Trucks

,

J. A.

Pople

, and

M.

Head-Gordon

, “

A fifth-order perturbation comparison of electron correlation theories

,”

Chem. Phys. Lett.

157

,

479

–

483

(

1989

).

https://doi.org/10.1016/s0009-2614(89)87395-6

Google Scholar

Crossref

6.

Q.

Ma

and

H. J.

Werner

, “

Explicitly correlated local coupled-cluster methods using pair natural orbitals

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1371

(

2018

).

https://doi.org/10.1002/wcms.1371

Google Scholar

Crossref

7.

P. R.

Nagy

and

M.

Kállay

, “

Approaching the basis set limit of CCSD(T) energies for large molecules with local natural orbital coupled-cluster methods

,”

J. Chem. Theory Comput.

15

,

5275

–

5298

(

2019

).

https://doi.org/10.1021/acs.jctc.9b00511

Google Scholar

Crossref

8.

D. G.

Liakos

,

Y.

Guo

, and

F.

Neese

, “

Comprehensive benchmark results for the domain based local pair natural orbital coupled cluster method (DLPNO-CCSD(T)) for closed- and open-shell systems

,”

J. Phys. Chem. A

124

,

90

–

100

(

2020

).

https://doi.org/10.1021/acs.jpca.9b05734

Google Scholar

Crossref

9.

M.

Kállay

,

R. A.

Horváth

,

L.

Gyevi-Nagy

, and

P. R.

Nagy

, “

Basis set limit CCSD(T) energies for extended molecules via a reduced-cost explicitly correlated approach

,”

J. Chem. Theory Comput.

19

,

174

–

189

(

2023

).

https://doi.org/10.1021/acs.jctc.2c01031

Google Scholar

Crossref

10.

V.

Barone

,

G.

Ceselin

,

M.

Fusè

, and

N.

Tasinato

, “

Accuracy meets interpretability for computational spectroscopy by means of hybrid and double-hybrid functionals

,”

Front. Chem.

8

,

584203

(

2020

).

https://doi.org/10.3389/fchem.2020.584203

Google Scholar

Crossref

PubMed

11.

Q.

Yang

,

M.

Mendolicchio

,

V.

Barone

, and

J.

Bloino

, “

Accuracy and reliability in the simulation of vibrational spectra: A comprehensive benchmark of energies and intensities issuing from generalized vibrational perturbation theory to second order (GVPT2)

,”

Front. Astron. Space Sci.

8

,

665232

(

2021

).

https://doi.org/10.3389/fspas.2021.665232

Google Scholar

Crossref

12.

H. H.

Nielsen

, “

The vibration-rotation energies of molecules

,”

Rev. Mod. Phys.

23

,

90

–

136

(

1951

).

https://doi.org/10.1103/revmodphys.23.90

Google Scholar

Crossref

13.

I. M.

Mills

, “

Vibration-rotation structure in asymmetric- and symmetric-top molecules

,” in

Molecular Spectroscopy: Modern Research

, edited by

K. N.

Rao

and

C. W.

Mathews

(

Academic Press

,

New York

,

1972

), Chap. 3.2, pp.

115

–

140

.

Google Scholar

Crossref

14.

D. A.

Clabo

, Jr.,

W. D.

Allen

,

R. B.

Remington

,

Y.

Yamaguchi

, and

H. F.

Schaefer

III, “

A systematic study of molecular vibrational anharmonicity and vibration-rotation interaction by self-consistent-field higher-derivative methods. Asymmetric top molecules

,”

Chem. Phys.

123

,

187

–

239

(

1988

).

https://doi.org/10.1016/0301-0104(88)87271-9

Google Scholar

Crossref

15.

W. D.

Allen

,

Y.

Yamaguchi

,

A. G.

Császár

,

D. A.

Clabo

, Jr.,

R. B.

Remington

, and

H. F.

Schaefer

III, “

A systematic study of molecular vibrational anharmonicity and vibration-rotation interaction by self-consistent-fied higher-derivative methods. Linear polyatomic molecules

,”

Chem. Phys.

145

,

427

–

466

(

1990

).

https://doi.org/10.1016/0301-0104(90)87051-C

Google Scholar

Crossref

16.

V.

Barone

, “

Vibrational zero-point energies and thermodynamic functions beyond the harmonic approximation

,”

J. Chem. Phys.

120

,

3059

–

3065

(

2004

).

https://doi.org/10.1063/1.1637580

Google Scholar

Crossref

17.

V.

Barone

, “

Anharmonic vibrational properties by a fully automated second-order perturbative approach

,”

J. Chem. Phys.

122

,

014108

(

2005

).

https://doi.org/10.1063/1.1824881

Google Scholar

Crossref

18.

S. V.

Krasnoshchekov

,

E. V.

Isayeva

, and

N. F.

Stepanov

, “

Numerical-analytic implementation of the higher-order canonical Van Vleck perturbation theory for the interpretation of medium-sized molecule vibrational spectra

,”

J. Phys. Chem. A

116

,

3691

–

3709

(

2012

).

https://doi.org/10.1021/jp211400w

Google Scholar

Crossref

19.

A. M.

Rosnik

and

W. F.

Polik

, “

VPT2+K spectroscopic constants and matrix elements of the transformed vibrational Hamiltonian of a polyatomic molecule with resonances using Van Vleck perturbation theory

,”

Mol. Phys.

112

,

261

–

300

(

2014

).

https://doi.org/10.1080/00268976.2013.808386

Google Scholar

Crossref

20.

P. R.

Franke

,

J. F.

Stanton

, and

G. E.

Douberly

, “

How to VPT2: Accurate and intuitive simulations of CH stretching infrared spectra using VPT2+K with large effective Hamiltonian resonance treatments

,”

J. Phys. Chem. A

125

,

1301

–

1324

(

2021

).

https://doi.org/10.1021/acs.jpca.0c09526

Google Scholar

Crossref

21.

C.

Peng

,

P. Y.

Ayala

,

H. B.

Schlegel

, and

M. J.

Frisch

, “

Using redundant internal coordinates to optimize equilibrium geometries and transition states

,”

J. Comput. Chem.

17

,

49

–

56

(

1996

).

https://doi.org/10.1002/(sici)1096-987x(19960115)17:1<49::aid-jcc5>3.3.co;2-#

Google Scholar

Crossref

22.

V.

Barone

, “

Quantum chemistry meets high-resolution spectroscopy for characterizing the molecular bricks of life in the gas-phase

,”

Phys. Chem. Chem. Phys.

26

,

5802

–

5821

(

2024

).

https://doi.org/10.1039/d3cp05169b

Google Scholar

Crossref

PubMed

23.

F.

Lazzari

,

A.

Salvadori

,

G.

Mancini

, and

V.

Barone

, “

Molecular perception for visualization and computation: The proxima library

,”

J. Chem. Inf. Model.

60

,

2668

–

2672

(

2020

).

https://doi.org/10.1021/acs.jcim.0c00076

Google Scholar

Crossref

24.

C.

Bannwarth

,

S.

Ehlert

, and

S.

Grimme

, “

GFN2-xTB—An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions

,”

J. Chem. Theory Comput.

15

,

1652

–

1671

(

2019

).

https://doi.org/10.1021/acs.jctc.8b01176

Google Scholar

Crossref

25.

P.

Pracht

,

F.

Bohle

, and

S.

Grimme

, “

Automated exploration of the low-energy chemical space with fast quantum chemical methods

,”

Phys. Chem. Chem. Phys.

22

,

7169

–

7192

(

2020

).

https://doi.org/10.1039/c9cp06869d

Google Scholar

Crossref

PubMed

26.

G.

Mancini

,

M.

Fusè

,

F.

Lazzari

,

B.

Chandramouli

, and

V.

Barone

, “

Unsupervised search of low-lying conformers with spectroscopic accuracy: A two-step algorithm rooted into the island model evolutionary algorithm

,”

J. Chem. Phys.

153

,

124110

(

2020

).

https://doi.org/10.1063/5.0018314

Google Scholar

Crossref

27.

D.

Ferro-Costas

,

I.

Mosquera-Lois

, and

A.

Fernandez-Ramos

, “

TorsiFlex: An automatic generator of torsional conformers. Application to the twenty proteinogenic amino acids

,”

J. Cheminf.

13

,

100

(

2021

).

https://doi.org/10.1186/s13321-021-00578-0

Google Scholar

Crossref

28.

G.

Mancini

,

M.

Fusè

,

F.

Lazzari

, and

V.

Barone

, “

Fast exploration of potential energy surfaces with a joint venture of quantum chemistry, evolutionary algorithms and unsupervised learning

,”

Digital Discovery

1

,

790

–

805

(

2022

).

https://doi.org/10.1039/D2DD00070A

Google Scholar

Crossref

29.

V.

Barone

, “

Accuracy meets feasibility for the structures and rotational constants of the molecular bricks of life: A joint venture of DFT and wave-function methods

,”

J. Phys. Chem. Lett.

14

,

5883

–

5890

(

2023

).

https://doi.org/10.1021/acs.jpclett.3c01380

Google Scholar

Crossref

30.

V.

Barone

, “

PCS/Bonds and PCS0: Pick your molecule and get its accurate structure and ground state rotational constants at DFT cost

,”

J. Chem. Phys.

159

,

081102

(

2023

).

https://doi.org/10.1063/5.0167296

Google Scholar

Crossref

31.

S.

Di Grande

and

V.

Barone

, “

Toward accurate quantum chemical methods for molecules of increasing dimension: The new family of Pisa composite schemes

,”

J. Phys. Chem. A

128

,

4886

–

4900

(

2024

).

https://doi.org/10.1021/acs.jpca.4c01673

Google Scholar

Crossref

32.

D. G.

Truhlar

and

A. D.

Isaacson

, “

Simple perturbation theory estimates of equilibrium constants from force fields

,”

J. Chem. Phys.

94

,

357

–

359

(

1991

).

https://doi.org/10.1063/1.460350

Google Scholar

Crossref

33.

E.

Dzib

and

G.

Merino

, “

The hindered rotor theory: A review

,”

Wiley Interdiscip. Rev.:Comput. Mol. Sci.

12

,

e1583

(

2022

).

https://doi.org/10.1002/wcms.1583

Google Scholar

Crossref

34.

V.

Barone

, “

From perception to prediction and interpretation: Enlightening the gray zone of molecular bricks of life with the help of machine learning and quantum chemistry

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

15

,

e70000

(

2025

).

https://doi.org/10.1002/wcms.70000

Google Scholar

Crossref

35.

M.

Mendolicchio

,

J.

Bloino

, and

V.

Barone

, “

Perturb-then-diagonalize vibrational engine exploiting curvilinear internal coordinates

,”

J. Chem. Theory Comput.

18

,

7603

–

7619

(

2022

).

https://doi.org/10.1021/acs.jctc.2c00773

Google Scholar

Crossref

36.

M.

Mendolicchio

, “

Harnessing the power of curvilinear internal coordinates: From molecular structure prediction to vibrational spectroscopy

,”

Theor. Chem. Acc.

142

,

133

(

2023

).

https://doi.org/10.1007/s00214-023-03069-7

Google Scholar

Crossref

37.

M.

Mendolicchio

and

V.

Barone

, “

Accurate vibrational and ro-vibrational contributions to the properties of large molecules by a new engine employing curvilinear internal coordinates and vibrational perturbation theory to second order

,”

J. Chem. Theory Comput.

20

,

8378

–

8395

(

2024

).

https://doi.org/10.1021/acs.jctc.4c00857

Google Scholar

Crossref

38.

P. Y.

Ayala

and

H. B.

Schlegel

, “

Identification and treatment of internal rotation in normal mode vibrational analysis

,”

J. Chem. Phys.

108

,

2314

–

2325

(

1998

).

https://doi.org/10.1063/1.475616

Google Scholar

Crossref

39.

V.

Barone

,

J.

Lupi

,

Z.

Salta

, and

N.

Tasinato

, “

Development and validation of a parameter-free model chemistry for the computation of reliable reaction rates

,”

J. Chem. Theory Comput.

17

,

4913

–

4928

(

2021

).

https://doi.org/10.1021/acs.jctc.1c00406

Google Scholar

Crossref

40.

M. A.

Harthcock

and

J.

Laane

, “

Calculation of two-dimensional vibrational potential energy surfaces utilizing prediagonalized basis sets and Van Vleck perturbation methods

,”

J. Phys. Chem.

89

,

4231

–

4240

(

1985

).

https://doi.org/10.1021/j100266a017

Google Scholar

Crossref

41.

V.

Barone

,

M.

Biczysko

,

J.

Bloino

,

M.

Borkowska-Panek

,

I.

Carnimeo

, and

P.

Panek

, “

Toward anharmonic computations of vibrational spectra for large molecular systems

,”

Int. J. Quantum Chem.

112

,

2185

–

2200

(

2012

).

https://doi.org/10.1002/qua.23224

Google Scholar

Crossref

42.

M.

Fusè

,

G.

Mazzeo

,

G.

Longhi

,

S.

Abbate

,

Q.

Yang

, and

J.

Bloino

, “

Scaling-up VPT2: A feasible route to include anharmonic correction on large molecules

,”

Spectrochim. Acta, Part A

311

,

123969

(

2024

).

https://doi.org/10.1016/j.saa.2024.123969

Google Scholar

Crossref

43.

T.

Carrington

, Jr., “

Using collocation to study the vibrational dynamics of molecules

,”

Spectrochim. Acta, Part A

248

,

119158

(

2021

).

https://doi.org/10.1016/j.saa.2020.119158

Google Scholar

Crossref

44.

S.

Manzhos

,

X.

Wang

, and

T.

Carrington

, Jr., “

A multimode-like scheme for selecting the centers of Gaussian basis functions when computing vibrational spectra

,”

Chem. Phys.

509

,

139

–

144

(

2018

).

https://doi.org/10.1016/j.chemphys.2017.10.006

Google Scholar

Crossref

45.

A. G.

Császár

,

C.

Fábri

,

T.

Szidarovszky

,

E.

Mátyus

,

T.

Furtenbacher

, and

G.

Czakó

, “

The fourth age of quantum chemistry: Molecules in motion

,”

Phys. Chem. Chem. Phys.

14

,

1085

–

1106

(

2012

).

https://doi.org/10.1039/c1cp21830a

Google Scholar

Crossref

PubMed

46.

E.

Mátyus

,

G.

Czakó

, and

A. G.

Császár

, “

Toward black-box-type full- and reduced-dimensional variational (ro)vibrational computations

,”

J. Chem. Phys.

130

,

134112

(

2009

).

https://doi.org/10.1063/1.3076742

Google Scholar

Crossref

47.

D.

Lauvergnat

and

A.

Nauts

, “

Exact numerical computation of a kinetic energy operator in curvilinear coordinates

,”

J. Chem. Phys.

116

,

8560

–

8570

(

2002

).

https://doi.org/10.1063/1.1469019

Google Scholar

Crossref

48.

D.

Lauvergnat

,

E.

Baloïtcha

,

G.

Dive

, and

M.

Desouter-Lecomte

, “

Dynamics of complex molecular systems with numerical kinetic energy operators in generalized coordinates

,”

Chem. Phys.

326

,

500

–

508

(

2006

).

https://doi.org/10.1016/j.chemphys.2006.03.012

Google Scholar

Crossref

49.

Y.

Scribano

,

D. M.

Lauvergnat

, and

D. M.

Benoit

, “

Fast vibrational configuration interaction using generalized curvilinear coordinates and self-consistent basis

,”

J. Chem. Phys.

133

,

094103

(

2010

).

https://doi.org/10.1063/1.3476468

Google Scholar

Crossref

50.

J.

Tennyson

, “

Perspective: Accurate ro-vibrational calculations on small molecules

,”

J. Chem. Phys.

145

,

124112

(

2016

).

https://doi.org/10.1063/1.4962907

Google Scholar

Crossref

51.

B. T.

Sutcliffe

and

J.

Tennyson

, “

A general treatment of vibration-rotation coordinates for triatomic molecules

,”

Int. J. Quantum Chem.

39

,

183

–

196

(

1991

).

https://doi.org/10.1002/qua.560390208

Google Scholar

Crossref

52.

S. N.

Yurchenko

,

W.

Thiel

, and

P.

Jensen

, “

Theoretical ROVibrational Energies (TROVE): A robust numerical approach to the calculation of rovibrational energies for polyatomic molecules

,”

J. Mol. Spectrosc.

245

,

126

–

140

(

2007

).

https://doi.org/10.1016/j.jms.2007.07.009

Google Scholar

Crossref

53.

A. S.

Petit

and

A. B.

McCoy

, “

Diffusion Monte Carlo in internal coordinates

,”

J. Phys. Chem. A

117

,

7009

–

7018

(

2013

).

https://doi.org/10.1021/jp312710u

Google Scholar

Crossref

54.

I. W.

Bulik

,

M. J.

Frisch

, and

P. H.

Vaccaro

, “

Fixed-node, importance-sampling diffusion Monte Carlo for vibrational structure with accurate and compact trial states

,”

J. Chem. Theory Comput.

14

,

1554

–

1563

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00016

Google Scholar

Crossref

55.

K.

Yagi

,

M.

Keçeli

, and

S.

Hirata

, “

Optimized coordinates for anharmonic vibrational structure theories

,”

J. Chem. Phys.

137

,

204118

(

2012

).

https://doi.org/10.1063/1.4767776

Google Scholar

Crossref

56.

P. M.

Zimmerman

and

P.

Smereka

, “

Optimizing vibrational coordinates to modulate intermode coupling

,”

J. Chem. Theory Comput.

12

,

1883

–

1891

(

2016

).

https://doi.org/10.1021/acs.jctc.5b01168

Google Scholar

Crossref

57.

B.

Thomsen

,

K.

Yagi

, and

O.

Christiansen

, “

Optimized coordinates in vibrational coupled cluster calculations

,”

J. Chem. Phys.

140

,

154102

(

2014

).

https://doi.org/10.1063/1.4870775

Google Scholar

Crossref

58.

P. S. T.

Arnaud Leclerc

and

T.

Carrington

, “

Comparison of different eigensolvers for calculating vibrational spectra using low-rank, sum-of-product basis functions

,”

Mol. Phys.

115

,

1740

–

1749

(

2017

).

https://doi.org/10.1080/00268976.2016.1249980

Google Scholar

Crossref

59.

S. R.

White

, “

Density matrix formulation for quantum renormalization groups

,”

Phys. Rev. Lett.

69

,

2863

–

2866

(

1992

).

https://doi.org/10.1103/physrevlett.69.2863

Google Scholar

Crossref

PubMed

60.

U.

Schollwöck

, “

The density-matrix renormalization group in the age of matrix product states

,”

Ann. Phys.

326

,

96

–

192

(

2011

).

https://doi.org/10.1016/j.aop.2010.09.012

Google Scholar

Crossref

61.

A.

Baiardi

and

M.

Reiher

, “

The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges

,”

J. Chem. Phys.

152

,

040903

(

2020

).

https://doi.org/10.1063/1.5129672

Google Scholar

Crossref

62.

A.

Baiardi

,

C. J.

Stein

,

V.

Barone

, and

M.

Reiher

, “

Vibrational density matrix renormalization group

,”

J. Chem. Theory Comput.

13

,

3764

–

3777

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00329

Google Scholar

Crossref

63.

N.

Glaser

,

A.

Baiardi

, and

M.

Reiher

, “

Flexible DMRG-based framework for anharmonic vibrational calculations

,”

J. Chem. Theory Comput.

19

,

9329

–

9343

(

2023

).

https://doi.org/10.1021/acs.jctc.3c00902

Google Scholar

Crossref

64.

M.

Zou

,

Y.

Hassan

,

T. K.

Roy

,

A. B.

McCoy

, and

M. I.

Lester

, “

Infrared spectroscopy of the syn-methyl-substituted criegee intermediate: A combined experimental and theoretical study

,”

J. Chem. Phys.

160

,

204309

(

2024

).

https://doi.org/10.1063/5.0210122

Google Scholar

Crossref

65.

A. B.

McCoy

and

M. A.

Boyer

, “

Exploring expansions of the potential and dipole surfaces used for vibrational perturbation theory

,”

J. Phys. Chem. A

126

,

7242

–

7249

(

2022

).

https://doi.org/10.1021/acs.jpca.2c05792

Google Scholar

Crossref

66.

S.

Carter

,

A. R.

Sharma

,

J. M.

Bowman

,

P.

Rosmus

, and

R.

Tarroni

, “

Calculations of rovibrational energies and dipole transition intensities for polyatomic molecules using multimode

,”

J. Chem. Phys.

131

,

224106

(

2009

).

https://doi.org/10.1063/1.3266577

Google Scholar

Crossref

67.

G.

Rauhut

and

T.

Hrenar

, “

A combined variational and perturbational study on the vibrational spectrum of P₂F₄

,”

Chem. Phys.

346

,

160

–

166

(

2008

).

https://doi.org/10.1016/j.chemphys.2008.01.039

Google Scholar

Crossref

68.

O.

Christiansen

, “

Selected new developments in vibrational structure theory: Potential construction and vibrational wave function calculations

,”

Phys. Chem. Chem. Phys.

14

,

6672

–

6687

(

2012

).

https://doi.org/10.1039/c2cp40090a

Google Scholar

Crossref

PubMed

69.

D.

Papp

,

T.

Szidarovszky

, and

A. G.

Császár

, “

A general variational approach for computing rovibrational resonances of polyatomic molecules. application to the weakly bound H₂He⁺ and H₂·CO systems

,”

J. Chem. Phys.

147

,

094106

(

2017

).

https://doi.org/10.1063/1.5000680

Google Scholar

Crossref

70.

A.

Erba

,

J.

Maul

,

M.

Ferrabone

,

P.

Carbonnière

,

M.

Rérat

, and

R.

Dovesi

, “

Anharmonic vibrational states of solids from DFT calculations. Part I: Description of the potential energy surface

,”

J. Chem. Theory Comput.

15

,

3755

–

3765

(

2019

).

https://doi.org/10.1021/acs.jctc.9b00293

Google Scholar

Crossref

71.

A.

Erba

,

J.

Maul

,

M.

Ferrabone

,

R.

Dovesi

,

M.

Rérat

, and

P.

Carbonnière

, “

Anharmonic vibrational states of solids from DFT calculations. Part II: Implementation of the VSCF and VCI methods

,”

J. Chem. Theory Comput.

15

,

3766

–

3777

(

2019

).

https://doi.org/10.1021/acs.jctc.9b00294

Google Scholar

Crossref

72.

D.

Mitoli

,

J.

Maul

, and

A.

Erba

, “

Anharmonic terms of the potential energy surface: A group theoretical approach

,”

Cryst. Growth Des.

23

,

3671

–

3680

(

2023

).

https://doi.org/10.1021/acs.cgd.3c00104

Google Scholar

Crossref

73.

R. G.

Schireman

,

J.

Maul

,

A.

Erba

, and

M. T.

Ruggiero

, “

Anharmonic coupling of stretching vibrations in ice: A periodic VSCF and VCI description

,”

J. Chem. Theory Comput.

18

,

4428

–

4437

(

2022

).

https://doi.org/10.1021/acs.jctc.2c00217

Google Scholar

Crossref

74.

M.

Born

and

R.

Oppenheimer

, “

Zur quantentheorie der molekeln

,”

Ann. Phys.

389

,

457

–

484

(

1927

).

https://doi.org/10.1002/andp.19273892002

Google Scholar

Crossref

75.

E.

Mátyus

,

A. M.

Santa Daría

, and

G.

Avila

, “

Exact quantum dynamics developments for floppy molecular systems and complexes

,”

Chem. Commun.

59

,

366

–

381

(

2023

).

https://doi.org/10.1039/d2cc05123k

Google Scholar

Crossref

76.

S.

Yurchenko

,

Computational Spectroscopy of Polyatomic Molecules

(

CRC Press

,

2023

).

Google Scholar

Crossref

77.

R.

Meyer

and

H. H.

Gunthard

, “

General internal motion of molecules, classical and quantum-mechanical Hamiltonian

,”

J. Chem. Phys.

49

,

1510

–

1520

(

1968

).

https://doi.org/10.1063/1.1670272

Google Scholar

Crossref

78.

H. M.

Pickett

, “

Vibration—rotation interactions and the choice of rotating axes for polyatomic molecules

,”

J. Chem. Phys.

56

,

1715

–

1723

(

1972

).

https://doi.org/10.1063/1.1677430

Google Scholar

Crossref

79.

B.

Podolsky

, “

Quantum-mechanically correct form of Hamiltonian function for conservative systems

,”

Phys. Rev.

32

,

812

(

1928

).

https://doi.org/10.1103/physrev.32.812

Google Scholar

Crossref

80.

E. B.

Wilson

, Jr., “

Some mathematical methods for the study of molecular vibrations

,”

J. Chem. Phys.

9

,

76

–

84

(

1941

).

https://doi.org/10.1063/1.1750829

Google Scholar

Crossref

81.

S.

Califano

,

Vibrational States

(

John Wiley and Sons

,

1976

).

Google Scholar

82.

J. H.

Van Vleck

, “

On σ-type doubling and electron spin in the spectra of diatomic molecules

,”

Phys. Rev.

33

,

467

(

1929

).

https://doi.org/10.1103/physrev.33.467

Google Scholar

Crossref

83.

S. V.

Krasnoshchekov

,

E. V.

Isayeva

, and

N. F.

Stepanov

, “

Criteria for first- and second-order vibrational resonances and correct evaluation of the Darling-Dennison resonance coefficients using the canonical Van Vleck perturbation theory

,”

J. Chem. Phys.

141

,

234114

(

2014

).

https://doi.org/10.1063/1.4903927

Google Scholar

Crossref

84.

C. R.

Quade

, “

Internal coordinate formulation for the vibration–rotation energies of polyatomic molecules

,”

J. Chem. Phys.

64

,

2783

–

2795

(

1976

).

https://doi.org/10.1063/1.432577

Google Scholar

Crossref

85.

A. D.

Isaacson

, “

Including anharmonicity in the calculation of rate constants. 1. The HCN/HNC isomerization reaction

,”

J. Phys. Chem. A

110

,

379

–

388

(

2006

).

https://doi.org/10.1021/jp058113y

Google Scholar

Crossref

86.

A. D.

Isaacson

, “

Including anharmonicity in the calculation of rate constants. II. The OH + H₂ → H₂O + H reaction

,”

J. Chem. Phys.

128

,

134304

(

2008

).

https://doi.org/10.1063/1.2834934

Google Scholar

Crossref

87.

E. L.

Sibert

, “

Theoretical studies of vibrationally excited polyatomic molecules using canonical Van Vleck perturbation theory

,”

J. Chem. Phys.

88

,

4378

–

4390

(

1988

).

https://doi.org/10.1063/1.453797

Google Scholar

Crossref

88.

E. L.

Sibert

, “

Modeling vibrational anharmonicity in infrared spectra of high frequency vibrations of polyatomic molecules

,”

J. Chem. Phys.

150

,

054107

(

2019

).

https://doi.org/10.1063/1.5079626

Google Scholar

Crossref

89.

M. A.

Boyer

and

A. B.

McCoy

, “

A flexible approach to vibrational perturbation theory using sparse matrix methods

,”

J. Chem. Phys.

156

,

054107

(

2022

).

https://doi.org/10.1063/5.0080892

Google Scholar

Crossref

90.

J. M. L.

Martin

,

T. J.

Lee

,

P. R.

Taylor

, and

J.-P.

François

, “

The anharmonic force field of ethylene, C₂H₄, by means of accurate ab initio calculations

,”

J. Chem. Phys.

103

,

2589

–

2602

(

1995

).

https://doi.org/10.1063/1.469681

Google Scholar

Crossref

91.

S. V.

Krasnoshchekov

,

E. O.

Dobrolyubov

,

M. A.

Syzgantseva

, and

R. V.

Palvelev

, “

Rigorous vibrational fermi resonance criterion revealed: Two different approaches yield the same result

,”

Mol. Phys.

118

,

e1743887

(

2020

).

https://doi.org/10.1080/00268976.2020.1743887

Google Scholar

Crossref

92.

K. M.

Kuhler

,

D. G.

Truhlar

, and

A. D.

Isaacson

, “

General method for removing resonance singularities in quantum mechanical perturbation theory

,”

J. Chem. Phys.

104

,

4664

–

4671

(

1996

).

https://doi.org/10.1063/1.471161

Google Scholar

Crossref

93.

J.

Bloino

,

M.

Biczysko

, and

V.

Barone

, “

General perturbative approach for spectroscopy, thermodynamics, and kinetics: Methodological background and benchmark studies

,”

J. Chem. Theory Comput.

8

,

1015

–

1036

(

2012

).

https://doi.org/10.1021/ct200814m

Google Scholar

Crossref

94.

M.

Mendolicchio

,

J.

Bloino

, and

V.

Barone

, “

General perturb-then-diagonalize model for the vibrational frequencies and intensities of molecules belonging to abelian and non-abelian symmetry groups

,”

J. Chem. Theory Comput.

17

,

4332

–

4358

(

2021

).

https://doi.org/10.1021/acs.jctc.1c00240

Google Scholar

Crossref

95.

M. S.

Schuurman

,

W. D.

Allen

,

P.

von Ragué Schleyer

, and

H. F.

Schaefer

III, “

The highly anharmonic BH₅ potential energy surface characterized in the ab initio limit

,”

J. Chem. Phys.

122

,

104302

(

2005

).

https://doi.org/10.1063/1.1853377

Google Scholar

Crossref

96.

M. J.

Frisch

,

G. W.

Trucks

,

H. B.

Schlegel

,

G. E.

Scuseria

,

M. A.

Robb

,

J. R.

Cheeseman

,

G.

Scalmani

,

V.

Barone

,

G. A.

Petersson

,

H.

Nakatsuji

,

X.

Li

,

M.

Caricato

,

A. V.

Marenich

,

J.

Bloino

,

B. G.

Janesko

,

R.

Gomperts

,

B.

Mennucci

,

H. P.

Hratchian

,

J. V.

Ortiz

,

A. F.

Izmaylov

,

J. L.

Sonnenberg

,

D.

Williams-Young

,

F.

Ding

,

F.

Lipparini

,

F.

Egidi

,

J.

Goings

,

B.

Peng

,

A.

Petrone

,

T.

Henderson

,

D.

Ranasinghe

,

V. G.

Zakrzewski

,

J.

Gao

,

N.

Rega

,

G.

Zheng

,

W.

Liang

,

M.

Hada

,

M.

Ehara

,

K.

Toyota

,

R.

Fukuda

,

J.

Hasegawa

,

M.

Ishida

,

T.

Nakajima

,

Y.

Honda

,

O.

Kitao

,

H.

Nakai

,

T.

Vreven

,

K.

Throssell

,

J. A.

Montgomery

, Jr.,

J. E.

Peralta

,

F.

Ogliaro

,

M. J.

Bearpark

,

J. J.

Heyd

,

E. N.

Brothers

,

K. N.

Kudin

,

V. N.

Staroverov

,

T. A.

Keith

,

R.

Kobayashi

,

J.

Normand

,

K.

Raghavachari

,

A. P.

Rendell

,

J. C.

Burant

,

S. S.

Iyengar

,

J.

Tomasi

,

M.

Cossi

,

J. M.

Millam

,

M.

Klene

,

C.

Adamo

,

R.

Cammi

,

J. W.

Ochterski

,

R. L.

Martin

,

K.

Morokuma

,

O.

Farkas

,

J. B.

Foresman

, and

D. J.

Fox

, Gaussian 16, Revision C.01,

Gaussian, Inc.

,

Wallingford, CT

,

2016

.

97.

A. D.

Becke

, “

Density-functional exchange-energy approximation with correct asymptotic behavior

,”

Phys. Rev. A

38

,

3098

–

3100

(

1988

).

https://doi.org/10.1103/physreva.38.3098

Google Scholar

Crossref

98.

C.

Lee

,

W.

Yang

, and

R. G.

Parr

, “

Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density

,”

Phys. Rev. B

37

,

785

–

789

(

1988

).

https://doi.org/10.1103/physrevb.37.785

Google Scholar

Crossref

99.

A. D.

Becke

, “

Density-functional thermochemistry. III. The role of exact exchange

,”

J. Chem. Phys.

98

,

5648

–

5652

(

1993

).

https://doi.org/10.1063/1.464913

Google Scholar

Crossref

100.

G.

Santra

,

N.

Sylvetsky

, and

J. M. L.

Martin

, “

Minimally empirical double-hybrid functionals trained against the GMTKN55 database: revDSD-PBEP86-D4, revDOD-PBE-D4, and DOD-SCAN-D4

,”

J. Phys. Chem. A

123

,

5129

–

5143

(

2019

).

https://doi.org/10.1021/acs.jpca.9b03157

Google Scholar

Crossref

101.

S.

Grimme

, “

Density functional theory with London dispersion corrections

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

1

,

211

–

228

(

2011

).

https://doi.org/10.1002/wcms.30

Google Scholar

Crossref

102.

E.

Caldeweyher

,

C.

Bannwarth

, and

S.

Grimme

, “

Extension of the D3 dispersion coefficient model

,”

J. Chem. Phys.

147

,

034112

(

2017

).

https://doi.org/10.1063/1.4993215

Google Scholar

Crossref

103.

W. J.

Hehre

,

R.

Ditchfield

, and

J. A.

Pople

, “

Self—consistent molecular orbital methods. XII. Further extensions of Gaussian—type basis sets for use in molecular orbital studies of organic molecules

,”

J. Chem. Phys.

56

,

2257

–

2261

(

1972

).

https://doi.org/10.1063/1.1677527

Google Scholar

Crossref

104.

P. C.

Hariharan

and

J. A.

Pople

, “

The influence of polarization functions on molecular orbital hydrogenation energies

,”

Theor. Chim. Acta

28

,

213

–

222

(

1973

).

https://doi.org/10.1007/bf00533485

Google Scholar

Crossref

105.

T.

Clark

,

J.

Chandrasekhar

,

G. W.

Spitznagel

, and

P. V. R.

Schleyer

, “

Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li–F

,”

J. Comput. Chem.

4

,

294

–

301

(

1983

).

https://doi.org/10.1002/jcc.540040303

Google Scholar

Crossref

106.

K. A.

Peterson

,

T. B.

Adler

, and

H.-J.

Werner

, “

Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, B–Ne, and Al–Ar

,”

J. Chem. Phys.

128

,

084102

(

2008

).

https://doi.org/10.1063/1.2831537

Google Scholar

Crossref

107.

T. H.

Dunning

, “

Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen

,”

J. Chem. Phys.

90

,

1007

–

1023

(

1989

).

https://doi.org/10.1063/1.456153

Google Scholar

Crossref

108.

G.

Rauhut

,

G.

Knizia

, and

H. J.

Werner

, “

Accurate calculation of vibrational frequencies using explicitly correlated coupled-cluster theory

,”

J. Chem. Phys.

130

,

054105

(

2009

).

https://doi.org/10.1063/1.3070236

Google Scholar

Crossref

109.

T.

Nakanaga

,

S.

Kondo

, and

S.

Saëki

, “

Infrared band intensities of formaldehyde and formaldehyde-d₂

,”

J. Chem. Phys.

76

,

3860

–

3865

(

1982

).

https://doi.org/10.1063/1.443527

Google Scholar

Crossref

110.

L. V.

Gurvich

,

I. V.

Veyts

, and

C. B.

Alcock

,

Thermodynamic Properties of Individual Substances

(

Hemisphere Publishing Co.

,

New York

,

1989

).

Google Scholar

111.

B. T.

Darling

and

D. M.

Dennison

, “

The water vapor molecule

,”

Phys. Rev.

57

,

128

–

139

(

1940

).

https://doi.org/10.1103/physrev.57.128

Google Scholar

Crossref

112.

R. L.

Redington

, “

Vibrational spectra and normal coordinate analysis of isotopically labeled formic acid monomers

,”

J. Mol. Spectrosc.

65

,

171

–

189

(

1977

).

https://doi.org/10.1016/0022-2852(77)90186-2

Google Scholar

Crossref

113.

R. C.

Millikan

and

K. S.

Pitzer

, “

Infrared spectra and vibrational assignment of monomeric formic acid

,”

J. Chem. Phys.

27

,

1305

–

1308

(

1957

).

https://doi.org/10.1063/1.1743996

Google Scholar

Crossref

114.

T.

Shimanouchi

,

Tables of Molecular Vibrational Frequencies, Consolidated Volume I

(

National Institute of Standards and Technology

,

1972

).

Google Scholar

115.

D.

McNaughton

,

C. J.

Evans

,

S.

Lane

, and

C. J.

Nielsen

, “

The high-resolution FTIR far-infrared spectrum of formamide

,”

J. Mol. Spectrosc.

193

,

104

–

117

(

1999

).

https://doi.org/10.1006/jmsp.1998.7709

Google Scholar

Crossref

116.

V.

Barone

,

M.

Fusè

,

S. M. V.

Pinto

, and

N.

Tasinato

, “

A computational journey across nitroxide radicals: From structure to spectroscopic properties and beyond

,”

Molecules

26

,

7404

(

2021

).

https://doi.org/10.3390/molecules26237404

Google Scholar

Crossref

PubMed

117.

G.

Ceselin

,

Z.

Salta

,

J.

Bloino

,

N.

Tasinato

, and

V.

Barone

, “

Accurate quantum chemical spectroscopic characterization of glycolic acid: A route toward its astrophysical detection

,”

J. Phys. Chem. A

126

,

2373

–

2387

(

2022

).

https://doi.org/10.1021/acs.jpca.2c01419

Google Scholar

Crossref

118.

V.

Barone

,

M.

Biczysko

,

J.

Bloino

,

P.

Cimino

,

E.

Penocchio

, and

C.

Puzzarini

, “

CC/DFT route toward accurate structures and spectroscopic features for observed and elusive conformers of flexible molecules: Pyruvic acid as a case study

,”

J. Chem. Theory Comput.

11

,

4342

–

4363

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00580

Google Scholar

Crossref

119.

H.

Hollenstein

,

F.

Akermann

, and

H. H.

Günthard

, “

Vibrational analysis of pyruvic acid and D−, ¹³C− and ¹⁸O-labelled species: Matrix spectra, assignments, valence force field and normal coordinate analysis

,”

Spectrochim. Acta, Part A

34

,

1041

–

1063

(

1978

).

https://doi.org/10.1016/0584-8539(78)80127-5

Google Scholar

Crossref

120.

K. L.

Plath

,

K.

Takahashi

,

R. T.

Skodje

, and

V.

Vaida

, “

Fundamental and overtone vibrational spectra of gas-phase pyruvic acid

,”

J. Phys. Chem. A

113

,

7294

–

7303

(

2009

).

https://doi.org/10.1021/jp810687t

Google Scholar

Crossref

121.

R. A.

Shaw

,

C.

Castro

,

R.

Dutler

,

A.

Rauk

, and

H.

Wieser

, “

The vibrational spectra (100–1600 cm⁻¹) and scaled ab initio STO‐3G and 3‐21G harmonic force fields for norbornane, norbornene, and norbornadiene

,”

J. Chem. Phys.

89

,

716

–

731

(

1988

).

https://doi.org/10.1063/1.455195

Google Scholar

Crossref

122.

F.

Vazart

,

C.

Latouche

,

P.

Cimino

, and

V.

Barone

, “

Accurate infrared (IR) spectra for molecules containing the C≡N moiety by anharmonic computations with the double hybrid B2PLYP density functional

,”

J. Chem. Theory Comput.

11

,

4364

–

4369

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00638

Google Scholar

Crossref

123.

T.

Fornaro

,

D.

Burini

,

M.

Biczysko

, and

V.

Barone

, “

Hydrogen-bonding effects on infrared spectra from anharmonic computations: Uracil–water complexes and uracil dimers

,”

J. Phys. Chem. A

119

,

4224

–

4236

(

2015

).

https://doi.org/10.1021/acs.jpca.5b01561

Google Scholar

Crossref

124.

G.

Wenzel

,

I. R.

Cooke

,

P. B.

Changala

,

E. A.

Bergin

,

S.

Zhang

,

A. M.

Burkhardt

,

A. N.

Byrne

,

S. B.

Charnley

,

M. A.

Cordiner

,

M.

Duffy

,

Z. T. P.

Fried

,

H.

Gupta

,

M. S.

Holdren

,

A.

Lipnicky

,

R. A.

Loomis

,

H. T.

Shay

,

C. N.

Shingledecker

,

M. A.

Siebert

,

D. A.

Stewart

,

R. H. J.

Willis

,

C.

Xue

,

A. J.

Remijan

,

A. E.

Wendlandt

,

M. C.

McCarthy

, and

B. A.

McGuire

, “

Detection of interstellar 1-cyanopyrene: A four-ring polycyclic aromatic hydrocarbon

,”

Science

386

,

810

–

813

(

2024

).

https://doi.org/10.1126/science.adq6391

Google Scholar

Crossref

PubMed

125.

G.

Wenzel

,

T. H.

Speak

,

P. B.

Changala

,

R. H. J.

Willis

,

A. M.

Burkhardt

,

S.

Zhang

,

E. A.

Bergin

,

A. N.

Byrne

,

S. B.

Charnley

,

Z. T. P.

Fried

,

H.

Gupta

,

E.

Herbst

,

M. S.

Holdren

,

A.

Lipnicky

,

R. A.

Loomis

,

C. N.

Shingledecker

,

C.

Xue

,

A. J.

Remijan

,

A. E.

Wendlandt

,

M. C.

McCarthy

,

I. R.

Cooke

, and

B. A.

McGuire

, “

Detections of interstellar aromatic nitriles 2-cyanopyrene and 4-cyanopyrene in TMC-1

,”

Nat. Astron.

9

,

262

(

2024

).

https://doi.org/10.1038/s41550-024-02410-9

Google Scholar

Crossref

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/).

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License .

Vibrational second-order perturbation theory based on curvilinear coordinates: Thermochemical applications

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Vibro-rotational problem in curvilinear coordinates

B. Zero-point vibrational energy

C. Soft degrees of freedom

D. Black-box treatment of Fermi resonances for thermochemistry

III. IMPLEMENTATION

IV. COMPUTATIONAL DETAILS

V. RESULTS AND DISCUSSION

A. Validation

B. Toward large molecules

VI. CONCLUSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: CALCULATION OF THE EXTRA-POTENTIAL TERM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Vibrational second-order perturbation theory based on curvilinear coordinates: Thermochemical applications Open Access

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Vibro-rotational problem in curvilinear coordinates

B. Zero-point vibrational energy

C. Soft degrees of freedom

D. Black-box treatment of Fermi resonances for thermochemistry

III. IMPLEMENTATION

IV. COMPUTATIONAL DETAILS

V. RESULTS AND DISCUSSION

A. Validation

B. Toward large molecules

VI. CONCLUSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: CALCULATION OF THE EXTRA-POTENTIAL TERM

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Vibrational second-order perturbation theory based on curvilinear coordinates: Thermochemical applications