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)

Vibrational spectroscopy represents an active frontier for the identification and characterization of molecular species in the context of astrochemistry and astrobiology. As new missions will provide more data over broader ranges and at higher resolution, especially in the infrared region, which could be complemented with new spectrometers in the future, support from laboratory experiments and theory is crucial. In particular, computational spectroscopy is playing an increasing role in deepening our understanding of the origin and nature of the observed bands in extreme conditions characterizing the interstellar medium or some planetary atmospheres, not easily reproducible on Earth. In this connection, the best compromise between reliability, feasibility and ease of interpretation is still a matter of concern due to the interplay of several factors in determining the final spectral outcome, with larger molecular systems and non-covalent complexes further exacerbating the dichotomy between accuracy and computational cost. In this context, second-order vibrational perturbation theory (VPT2) together with density functional theory (DFT) has become particularly appealing. The well-known problem of the reliability of exchange-correlation functionals, coupled with the treatment of resonances in VPT2, represents a challenge for the determination of standardized or “black-box” protocols, despite successful examples in the literature. With the aim of getting a clear picture of the achievable accuracy and reliability of DFT-based VPT2 calculations, a multi-step study will be carried out here. Beyond the definition of the functional, the impact of the basis set and the influence of the resonance treatment in VPT2 will be analyzed. For a better understanding of the computational aspects and the results, a short summary of vibrational perturbation theory and the overall treatment of resonances for both energies and intensities will be given. The first part of the benchmark will focus on small molecules, for which very accurate experimental and theoretical data are available, to investigate electronic structure calculation methods. Beyond the reliability of energies, widely used for such systems, the issue of intensities will also be investigated in detail. The best performing electronic structure methods will then be used to treat larger molecular systems, with more complex topologies and resonance patterns.

1 Introduction

More than 200 molecules have been identified in the interstellar medium (ISM) during the last 50 years (McGuire, 2018), and a much larger number on comets or exoplanet atmospheres, thanks to their spectroscopic signatures. Among them, the so-called interstellar complex organic molecules (iCOMs), namely molecules containing more than five atoms with at least one carbon atom, are particularly significant because some of them are either biomolecular building blocks or their key precursors (Chyba and Sagan, 1992; Herbst and van Dishoeck, 2009; Hörst et al., 2012; Balucani, 2012; Saladino et al., 2012; Saladino et al., 2015).

While rotational spectroscopy plays a major role in characterizing molecules in the ISM (Puzzarini and Barone, 2020), this is not true on Earth and for exoplanet atmospheres, where the dominant regions are at shorter wavelengths. In this context, vibrational spectroscopies are often the methods of choice (Des Marais et al., 2002; Girish and Sony, 2008; Seager, 2014). Furthermore, chirality is one of the key features of life, at least as we know it. As a consequence, the study of prebiotic molecules cannot escape from this feature, and here chiroptical spectroscopies, (e.g., vibrational circular dichroism, VCD, and Raman optical activity, ROA) come into play. The intricacy of spectra and the difficulty of their interpretation due to the concomitant role of several factors has led to an increasing role of computational spectroscopy, which is able to match the extreme environments observable in space. Of course, a mandatory prerequisite to the use of computational approaches is their sufficient accuracy in order to prevent any risk of ambiguity in the interpretation. For the specific case of vibrational spectroscopies, this requirement necessarily includes the treatment of anharmonicity in nuclear motions, with the possible inclusion of the rotational sub-structure (Biczysko et al., 2018; Puzzarini et al., 2019). Beyond very sophisticated rovibronic models, still limited in practice to diatomic or triatomic systems in standard applications (Carter et al., 1990, 2000; Császár et al., 2012; Mitrushchenkov, 2012; Nauts and Lauvergnat, 2012; Yurchenko et al., 2016), variational methods in conjunction with electronic structure calculations at the coupled cluster (CC) or configuration interaction (CI) levels represent a gold standard (Biczysko et al., 2003; Mátyus et al., 2009; Papp et al., 2017; Erfort et al., 2020). Indeed, variational approaches can in principle be systematically improved to reach any target accuracy, by using higher-order force fields through an extensive sampling of the potential energy surface (PES). Besides the practical difficulties of attaining arbitrary precision, the convergence to the correct energy is generally slow, especially when the starting level is the harmonic approximation. In this case, an extensive description of the PES may be necessary, resulting in a high overall computational cost to obtain results in close agreement with the best experimental data.

As larger and more complex molecules gain greater astrochemical interest, for instance in connection with the issue of the exogenous origin of life, this cost can become prohibitive and less expensive alternatives are necessary. One possible strategy is to use a cheaper method to estimate a better starting point on which the variational correction is applied, like vibrational self-consistent field (VSCF). However, the efficiency of such an approach is strongly conditioned by the description of the nuclear vibrations (Bulik et al., 2017). Over the years, perturbative approaches, and in particular second-order vibrational perturbation theory (VPT2, Nielsen (1951)), have shown promising capabilities, offering significant improvements over the harmonic level at a more convenient cost than their variational counterparts (Puzzarini et al., 2019). From a computational perspective, the interest of these approaches is twofold. First, they require simpler descriptions of the PES, which can be achieved with a limited number of small displacements from the equilibrium along suitable coordinates. Second, few, analytic equations need to be derived and implemented to describe the vibrational energies and intensities, making their calculation very straightforward. As a result, VPT2 on top of electronic structures at a coupled-cluster level, in particular with single and double excitations, and a perturbative treatment of the triples (CCSD(T), Raghavachari et al. (1989)), has become a standard to get energy levels, notably in the mid-IR region. Several strategies have been proposed to lower the computer time required by the electronic structure calculations, which represent the bottleneck of accurate calculations for medium/large-size molecular systems, with hybrid schemes based on a separation of the harmonic problem and the anharmonic correction being particularly effective (Begue et al., 2005; Carbonniere et al., 2005; Puzzarini and Barone, 2008; Biczysko et al., 2018; Puzzarini et al., 2019). These approaches are rooted into the hypothesis that the largest source of errors in VPT2 results is related to the harmonic reference data, so that a good accuracy can be obtained by combining harmonic results produced by a very accurate electronic structure method, (e.g., CCSD(T)) with anharmonic contributions evaluated at a lower level, (e.g., DFT employing hybrid functionals) (Puzzarini et al., 2010; Barone et al., 2013a, 2015; Puzzarini et al., 2019). Of course, even harmonic contributions need to be computed by relatively cheap quantum mechanical methods, (e.g., DFT) for very large molecules, like polycyclic aromatic hydrocarbons or biomolecules (McGuire et al., 2021).

An important source of complications with DFT is the ever-growing range of exchange-correlation functionals, each with their own strengths and limitations. This has prompted an increasing activity of benchmarking in the literature to provide pointers for choosing the most suitable functionals for a given type of study (Goerigk et al., 2017; Mardirossian and Head-Gordon, 2017). Unfortunately, the quantities of interest have been predominantly structural parameters, energetics and thermodynamics, with vibrational spectroscopies limited to the harmonic level. Recently, the growing interest in VPT2 has led to more extensive studies with account for anharmonicity (Bloino et al., 2016; Barone et al., 2020; Mitra and Roy, 2020; Nejad et al., 2020). This paves the way to more robust computational protocols. However, two aspects remain insufficiently explored. First, intensities are still rarely considered, especially for Raman and chiral spectroscopies. As bioactivity in planetary atmospheres is being increasingly probed, the combination of multiple measurement techniques becomes appealing, which means having the capacity of predicting correctly multiple properties and their interaction, as it is the case in chiroptical spectroscopies. A further difficulty in assessing the quality of the DFT functionals and the effect of basis sets is that reference implementations in more sophisticated electronic structure calculation methods is often lacking, so direct comparison with experiment is the only viable path. On the other side, VPT2 can suffer from resonances, special conditions which lead to an incorrect account of the anharmonic effects. While several strategies have been proposed in the literature and automated schemes are being devised to overcome this issue (Martin et al., 1995; Kuhler et al., 1996; Barone, 2005; Bloino et al., 2012, 2015; Krasnoshchekov et al., 2014, 2020), the effect of these corrections on the final energies and intensities has still been scarcely investigated. With these considerations in mind, the study proposed here is divided in several steps, starting from small systems where VPT2 calculations are straightforward and accurate reference data are available. By increasing the size of the systems of interest and their structural complexity, more resonance patterns can appear, which can impact notably the results, independently of the quality of the underlying electronic structure calculation methods. The results of these studies will pave the way to the full simulation of vibrational spectra of larger systems.

The manuscript is organized as follows. First, key aspects of VPT2 for the calculation of energies and intensities, as well as the problem of resonances, are briefly recalled. Strategies to address the latter are presented and discussed. The computational protocol is then described, followed by an analysis of the quality of some representative functionals combined with a variety of basis sets. The impact of the VPT2 treatment is also assessed. Finally, larger molecular systems are used as case studies to investigate the overall quality of DFT and VPT2.

2 Theoretical Background

2.1 The VPT2 Framework

In this section, we highlight key aspects of VPT2, which can have important practical impacts on the results. This will also let us introduce the notation used in the manuscript. More details on the derivations of the equations can be found in References. (Nielsen, 1951; Plíva, 1990; Willets et al., 1990; Vázquez and Stanton, 2006; Bloino and Barone, 2012; Franke et al., 2021). For a similar reason, reference equations are recalled, but most of the concepts, especially regarding the problem of resonances, are introduced in a descriptive way, focusing on the physical concepts and implications.

First of all, it is noteworthy that a number of hypotheses are commonly made in the theoretical developments, namely that the Born-Oppenheimer approximation is valid for the separation of the electronic and nuclear wave functions, and the Eckart-Sayvetz conditions (Eckart, 1935) are met to separate translational and vibrational-rotational motions through specific position and orientation of the molecule, the Eckart frame. Among several alternatives, the most widely used starting point, chosen here as well, is the rovibrational Hamiltonian proposed by Watson and expressed in a basis of dimensionless normal coordinates $q$ and their conjugate momenta $p$ (Watson, 1968),

${\mathrm{\mathscr{H}}}_{rovib}=\frac{{\hslash }^2}{2hc}{\displaystyle \sum _{\tau ,\eta }{\mu }_{\tau \eta }}\left(J_{\tau }-{\pi }_{\tau }\right)\left(J_{\eta }-{\pi }_{\eta }\right)+\frac{1}{2}{\displaystyle \sum _{i=1}^N{\omega }_i{p}_i^2+V+U},(1)$

where h and c are respectively the Planck constant and the speed of light, $J$ and $\pi$, the rotational and vibrational angular momentum operators, and $\omega$ is the vector of harmonic wavenumbers. The summations run over the Cartesian coordinates of the Eckart frame, τ and η, and the N normal modes, i. U is a mass-dependent contribution, which vanishes for linear molecules (Bunker and Jensen, 1998),

$\mathit{U}=-\frac{{\hslash }^2}{\mathrm{8hc}}{\displaystyle \sum _{\tau }{\mu }_{\tau \tau }}.(2)$

V and μ are respectively the potential energy operator and the effective inverse molecular inertia tensor, for which analytic definitions are generally not available. Both can be expanded as Taylor series of $q$ about the equilibrium geometry (Császár, 2011),

${\mu }_{\tau \eta }={\mu }_{\tau \eta }^{\text{eq}}-{\displaystyle \sum _{i=1}^N{\mu }_{\tau \tau }^{\text{eq}}{\text{a}}_{i,\tau \eta }{\mu }_{\eta \eta }^{\text{eq}}{q}_i}+\frac{3}{4}{\displaystyle \sum _{\varsigma }{\displaystyle \sum _{i,j=1}^N{\mu }_{\tau \tau }^{\text{eq}}{\text{a}}_{i,\tau \varsigma }{\mu }_{\varsigma \varsigma }^{\text{eq}}{\text{a}}_{j,\varsigma \eta }{\mu }_{\eta \eta }^{\text{eq}}{q}_i{q}_j}}+\cdots ,(3)$

$V=\frac{1}{2}{\displaystyle \sum _{i=1}^N{\omega }_i{q}_i^2}+\frac{1}{6}{\displaystyle \sum _{i,j,k=1}^N{\text{f}}_{ijk}{q}_i{q}_j{q}_k}+\frac{1}{24}{\displaystyle \sum _{i,j,k,l=1}^N{\text{f}}_{ijkl}{q}_i{q}_j{q}_k{q}_l}+\cdots ,(4)$

with

${\text{a}}_{i,\tau \eta }=\frac{\partial {\text{I}}_{\tau \eta }^{\text{eq}}}{\partial q_i},{\text{f}}_{ijk}=\frac{{\partial }^{3}V}{\partial q_i\partial q_j\partial q_k},{\text{f}}_{ijkl}=\frac{{\partial }^{4}V}{\partial q_i\partial q_j\partial q_k\partial q_l}.$

${\text{f}}_{ijk}$ and ${\text{f}}_{ijkl}$ are commonly referred to as the cubic and quartic force constants, respectively, and ${\mathbf{I}}^{\text{eq}}$ is the equilibrium molecular inertia tensor.

In the present work, the systems under study are assumed to be in their rotational ground states, so the contributions of the rotational operator $J$ are ignored, and the development of the effective inverse molecular inertia tensor μ is truncated at the zeroth order. It is noteworthy that, even in these conditions, the coupling between vibrations and rotations does not vanish, due to the term,

$\frac{{\hslash }^2}{2hc}{\displaystyle \sum _{\tau ,\eta }{\mu }_{\tau \eta }^{\text{eq}}{\pi }_{\tau }{\pi }_{\eta }}={\displaystyle \sum _{\tau }{B}_{\tau }^{\text{eq}}}{\displaystyle \sum _{i,j,k,l=1}^N{\zeta }_{ij,\tau }{\zeta }_{kl,\tau }}\sqrt{\frac{{\omega }_j{\omega }_l}{{\omega }_i{\omega }_k}}{q}_i{p}_j{q}_k{p}_l,$

with ${\mathit{B}}^{\text{eq}}$ the equilibrium molecular rotational constants vector, and $\mathit{\zeta }$ the matrix of Coriolis couplings.

The vibrational Hamiltonian of interest here, truncated at the zeroth-order of Eq. 3 and second-order of Eq. 4, is thus,

${\mathrm{\mathscr{H}}}_{\text{vib}}=\frac{1}{2}{\displaystyle \sum _{i=1}^N{\omega }_i\left(p_i^2+q_i^2\right)}+\frac{1}{6}{\displaystyle \sum _{i,j,k=1}^N{\text{f}}_{ijk}{q}_i{q}_j{q}_k}+\frac{1}{24}{\displaystyle \sum _{i,j,k,l=1}^N{\text{f}}_{ijkl}{q}_i{q}_j{q}_k{q}_l}+{\displaystyle \sum _{\tau }{B}_{\tau }^{\text{eq}}}{\displaystyle \sum _{i,j,k,l=1}^N{\zeta }_{ij,\tau }{\zeta }_{kl,\tau }}\sqrt{\frac{{\omega }_j{\omega }_l}{{\omega }_i{\omega }_k}}{q}_i{p}_j{q}_k{p}_l+U.(5)$

The U term is usually accounted for through its zeroth-order expansion,

$U^{\left(0\right)}=-\text{Γ}{\displaystyle \sum _{\tau }\frac{B_{\tau }^{\text{eq}}}{4}}(6)$

where $\Gamma =0$ for linear molecules and $\Gamma =1$ otherwise. However, it does not play a role in transition energies, and will be ignored in the following.

Starting from the solution of the vibrational Schrödinger equation for the harmonic oscillator, labeled with the “(0)” superscript,

${\mathrm{\mathscr{H}}}_{\text{vib}}^{\left(0\right)}|{\psi }^{\left(0\right)}\rangle ={\epsilon }^{\left(0\right)}|{\psi }^{\left(0\right)}\rangle .(7)$

the vibrational wave function (ψ) and energy (ε) are expanded in perturbative series, and the Hamiltonian ${\mathrm{\mathscr{H}}}_{\text{vib}}^{\left(0\right)}$ replaced with the one in Eq. 5. Through a series of transformations, (e.g., Van Vleck. (1929); Papousek and Aliev. (1982)), a single formula is obtained to compute the VPT2 energy of any vibrational level. It should be noted that, by construction, the energies and wave functions are expressed in the initial normal-mode basis, so the vibrational states remain expressed as combinations of excited oscillators with associated quanta like it is done for harmonic levels, simply noted $|\mathit{v}\rangle$, where $\mathit{v}$ is the vector of vibrational quanta.

A similar approach can be applied for the intensities, where the transition moments of each property of interest, labeled generically as $\mathbf{P}$ for convenience here, are developed to generate analytic equations,

${\langle \mathbf{P}\rangle }_{I,F}=\frac{\langle {\psi }_I|{\mathbf{P}}^e|{\psi }_F\rangle }{\sqrt{\langle {\psi }_I|{\psi }_I\rangle \langle {\psi }_F|{\psi }_F\rangle }}.(8)$

The perturbed wave functions, ${\psi }_I$ and ${\psi }_F$, are not normalized, so the denominator cannot be simplified. Here, the full anharmonicity is considered, both of the wave function (mechanical) and of the property, so ${\mathbf{P}}^{\text{e}}$ itself is expanded as a Taylor series with respect to the normal coordinates. The derivation is quite tedious and depends on the initial (I) and final (F) states, but general equations can be obtained based on the nature of the initial and final states (Vázquez and Stanton, 2006; Bloino and Barone, 2012; Bloino, 2015).

2.2 Vibrational Energies

Let us consider first the case of vibrational energies. Following the resolution of the anharmonic vibrational Schrödinger equation at the VPT2 level, the vibrational energy of a state $|\mathit{v}\rangle$ is given by,

${\epsilon }_{\mathit{v}}={\epsilon }_0+{\displaystyle \sum _{i=1}^N{v}_i{\omega }_i}+{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=i}^N{\chi }_{ij}}}\left[v_i{v}_j+\frac{1}{2}\left(v_i+v_j\right)\right],(9)$

with ${\epsilon }_0$ the anharmonic zero-point vibrational energy.

In the present study, we are only interested in transition energies from the vibrational ground state, so,

${\nu }_{\mathit{v}}={\epsilon }_{\mathit{v}}-{\epsilon }_0={\displaystyle \sum _{i=1}^N{v}_i{\omega }_i}+{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=i}^N{\chi }_{ij}}}\left[v_i{v}_j+\frac{1}{2}\left(v_i+v_j\right)\right].(10)$

The anharmonic correction is contained in the $\chi$ matrix, whose elements are defined as,

$16{\chi }_{ii}={\text{f}}_{iiii}-\frac{5{\text{f}}_{iii}^2}{3{\omega }_i}-{\displaystyle \sum _{j=1\ne i}^N\frac{\left(8{\omega }_i^2-3{\omega }_j^2\right){\text{f}}_{iij}^2}{{\omega }_j\left(4{\omega }_i^2-{\omega }_j^2\right)},}(11)$

$4{\chi }_{ij}={\text{f}}_{iijj}-\frac{2{\omega }_i{\text{f}}_{iij}^2}{\left(4{\omega }_i^2-{\omega }_j^2\right)}-\frac{2{\omega }_j{\text{f}}_{ijj}^2}{\left(4{\omega }_j^2-{\omega }_i^2\right)}-\frac{{\text{f}}_{iii}{\text{f}}_{ijj}}{{\omega }_i}-\frac{{\text{f}}_{jjj}{\text{f}}_{iij}}{{\omega }_j}+{\displaystyle \sum _{k=1\ne i,j}^N\left[\frac{2{\omega }_k\left({\omega }_i^2+{\omega }_j^2-{\omega }_k^2\right){\text{f}}_{ijk}^2}{\left({\omega }_i+{\omega }_j+{\omega }_k\right)\left({\omega }_i-{\omega }_j-{\omega }_k\right)\left({\omega }_j-{\omega }_i-{\omega }_k\right)\left({\omega }_k-{\omega }_i-{\omega }_j\right)}-\frac{{\text{f}}_{iik}{\text{f}}_{jjk}}{{\omega }_k}\right]}+\frac{4\left({\omega }_i^2+{\omega }_j^2\right)}{{\omega }_i{\omega }_j}{\displaystyle \sum _{\tau =x,y,z}{B}_{\tau }^{\text{eq}}}{\left\{{\zeta }_{ij,\tau }\right\}}^2.(12)$

Equation 10 corresponds to the VPT2 energy. From a quick analysis of Eqs 11, 12, singularity conditions, collectively known as Fermi resonances (FRs), can be identified, where the energy can tend toward infinite values. For simplicity, let us consider a single term, ${\text{f}}_{iij}^2/8\left(2{\omega }_i-{\omega }_j\right)$, which arises from the development in partial fraction of $2{\omega }_i{\text{f}}_{iij}^2/4\left(4{\omega }_i^2-{\omega }_j^2\right)$, the other term being strictly non-resonant. Whenever $\left(2{\omega }_i-{\omega }_j\right)$ becomes negligible and ${\text{f}}_{iij}$ is not null, the contribution can in theory tend to infinity. Several strategies can be adopted to handle this problem. The most common schemes of potential interest here are reported in Table 1, and will be described below.

TABLE 1

TABLE 1. Description of the treatment of Fermi resonant terms in the calculation of energy, and Fermi and Darling–Dennison resonance in intensities for some VPT2 schemes. The condition indicates in which cases the alternative term is used in place of the potentially resonant term (reported for the original VPT2 scheme since no transformation is done). See text for details. ${\Sigma }_{ijk}={\omega }_i+{\omega }_j+{\omega }_k\text{and}{\Delta }_{ijk}={\omega }_i-{\omega }_j-{\omega }_k.$

In the deperturbed VPT2 (DVPT2) treatment, each potentially resonant term is screened, applying one or more criteria to establish if it is resonant and should be ignored, or not. The most common procedure involves two steps, first on the energetic proximity of the interacting states–here $|2_i\rangle$ and $|1_j\rangle$ –, that is the magnitude of the denominator, and then on the overall weight of the term. Several strategies have been proposed for the latter (Martin et al., 1995; Barone, 2005; Krasnoshchekov et al., 2020). Here, we will employ the test proposed by Martin and coworkers (Martin et al., 1995), which considers the deviation of the term from the variational energy of a model, ad hoc system. In practice, the overall procedure for the example case we have chosen can be summarized in two successive steps.

1. $\left|2{\omega }_i-{\omega }_j\right|\le \text{Δ}{\omega }^{1-2},$

2. ${\text{f}}_{iij}^4/256{\left|2{\omega }_i-{\omega }_j\right|}^3\ge K^{1-2}.$

As a result, the capacity of identifying correctly the resonances, and ultimately the quality of the DVPT2 results, will depend on the chosen thresholds, $\text{Δ}{\omega }^{1-2}$ and $K^{1-2}$. Conversely, no test is involved in degeneracy-corrected PT2 (DCPT2, Kuhler et al. (1996)). All potentially resonant terms are replaced by a non-resonant variant, as illustrated in Table 1 for the specific term discussed here. While satisfactory to treat resonance conditions, this model is ill-suited for well-separated, coupled states, where VPT2 is correct. HDCPT2 introduces a mixing of DCPT2 and the conventional VPT2 for each potentially resonant term (Bloino et al., 2012). Here again, no test is performed. Each term is computed twice, using the VPT2 and DCPT2 formulations, with the results combined by application of a merging function $\text{Λ}$ given in Table 1, where α and β are the mixing parameters, described in (Bloino et al., 2012).

Among these models, DVPT2 can lead to a truncated treatment since terms are selectively discarded. To recover them, a variational treatment is added, which reintroduces the terms as explicit couplings between the states with DVPT2 energies. In a simplified picture, a symmetric matrix is built, listing all fundamental $\left(|1_i\rangle \right)$, first overtone $\left(|2_i\rangle \right)$ and binary combination $\left(|1_i{1}_j\rangle \right)$ states as rows and columns (it is possible to add three-quanta states to reach NIR regions (Bloino, 2015), but this will not be discussed here). The DVPT2 energies are then put on the diagonal, while the coupling terms corresponding to the identified resonances, for instance here between $|1_j\rangle$ and $|2_i\rangle$, are added out of the diagonal. New energies are obtained by diagonalization of the resulting matrix. The full procedure will be referred to as generalized VPT2 (GVPT2, Barone (2005)). In practice, building the full matrix can be extremely expensive in terms of memory, since it scales as $N^4$ up to two-quanta transitions, and $N^6$ up to three quanta. A more practical and numerically stable solution is to build polyads, small ensembles of states connected through resonances, and diagonalize the corresponding matrices (Martin and Taylor, 1997). It should be noted that the diagonalization process introduces a mixing between the resonant states, which can depart significantly from their original definitions. If the transformation is small enough, the more intuitive oscillator-based notation used at the harmonic level and for standard VPT2 remains close to the true states. Otherwise, a different notation, similar to the one used in variational approaches, must be adopted.

To conclude on the vibrational energies, some couplings are not accounted for at the VPT2 level. In standard calculations, they are the interaction between first overtones and two-quanta binary combinations–generally noted 2–2 –, and between fundamentals (noted 1–1). As originally documented by Darling and Dennison (Darling and Dennison, 1940), they can give important contributions to vibrational energies. These terms, collectively known as Darling–Dennison resonances (DDRs) or interactions, can be added during the variational step in GVPT2 (Krasnoshchekov et al., 2014; Rosnik and Polik, 2014). The polyads are then incremented by including these new couplings between states, and the corresponding matrices completed with the new contributions. The introduction of these corrections can further transform the final states from the original VPT2 states. Because of the cost involved in computing and integrating these terms, and the limitations in the quality of the variational correction, only significant contributions are normally included. This can be done, like FRs in DVPT2, through a two-step procedure, first by considering only close states, and then by looking at the magnitude of the coupling to be included in the variational matrix (Bloino et al., 2016).

As a final remark, Fermi resonances-like conditions can be present in Darling–Dennison interaction terms. Because they do not correspond to energies, (i.e. the eigenvalues of a well-defined Hamiltonian), the strategies proposed above cannot be directly applied. The transformation inspired by HDCPT2 is used here, as described in (Bloino et al., 2016). Formally, DD resonances cannot be used in conjunction with DCPT2, and by extension HDCPT2, because the effective Hamiltonian corresponding to the generated vibrational energies is unknown. Since the coupling terms to be added do not overlap with Fermi resonances–but can be coupled to them through common states involved in both–it is possible to build a variational matrix with the HDCPT2 energies and only the off-diagonal terms related to DDRs included. The corresponding model will be called GHDCPT2, by analogy with GVPT2 for DVPT2/VPT2. Possible schemes available to compute the vibrational energies within VPT2 are summarized in Table 2.

TABLE 2

TABLE 2. Summary of the available VPT2 schemes to compute vibrational energies and intensities, and the strategies applied to deal with Fermi resonances and to include Darling–Dennison interaction terms. See text for details.

While relatively generic, the discussion deliberately omits the problem of true degeneracies present in linear-, symmetric- and spherical-top molecules. A proper treatment at the VPT2 level requires a special notation and an alternative derivation (Plíva, 1990; Piccardo et al., 2015), which has been explicitly carried out only for the energies of doubly degenerate representations. Because of these limitations and to avoid expanding too much this presentation, only asymmetric tops are considered in the following. The generalization of the VPT2 framework to any molecular symmetry through a unified framework is deferred to a forthcoming paper.

2.3 Transition Moments

Let us now consider the case of intensities, and the calculation of transition moments. Contrary to energies, multiple formulas are needed, depending on the type of transition, to fundamental, overtone or combination states, and the total number of quanta changing between the initial and final states (Willets et al., 1990; Vázquez and Stanton, 2006; Bloino and Barone, 2012; Bloino, 2015). Considering only transitions from the ground state, an interesting case is the transition moment associated to fundamental bands. The full equation can be found elsewhere (Bloino and Barone, 2012) and is reported in the Supplementary Material. Here, we consider only some characteristic terms.

$\text{A}:-\frac{s_1}{8}{\displaystyle \sum _{j,k=1}^N{\text{f}}_{ijk}}\left({\mathbf{P}}_{jk}+{\mathbf{P}}_{kj}\right)\left(\frac{1}{{\omega }_i+{\omega }_j+{\omega }_k}-\frac{\text{S}}{{\omega }_i-{\omega }_j-{\omega }_k}\right),$

$\text{B}:-\frac{s_0}{8}{\displaystyle \sum _{j,k=1}^N{\text{f}}_{ijkk}}{\mathbf{P}}_j\left[\frac{1}{{\omega }_i+{\omega }_j}-\frac{\text{S}\left(1-{\delta }_{ij}\right)}{{\omega }_i-{\omega }_j}\right],$

$C:+\frac{s_0}{16}{\displaystyle \sum _{j,k,l=1\ne i}^N{\text{f}}_{ikl}{\text{f}}_{jkl}}{\mathbf{P}}_j\left[\frac{1}{\left({\omega }_i+{\omega }_j\right)\left({\omega }_i+{\omega }_k+{\omega }_l\right)}-\frac{1}{\left({\omega }_i+{\omega }_j\right)\left({\omega }_i-{\omega }_k-{\omega }_l\right)}-\frac{\text{S}}{\left({\omega }_i-{\omega }_j\right)\left({\omega }_i+{\omega }_k+{\omega }_l\right)}+\frac{\text{S}}{\left({\omega }_i-{\omega }_j\right)\left({\omega }_i-{\omega }_k-{\omega }_l\right)}\right].$

In the above equations, S = 1 if the property is a function of the normal coordinates, like the electric dipole, (e.g., for IR) or the polarizability tensor (for Raman), and S = −1 if the property depend on their conjugate momenta, like the atomic axial tensor, (e.g., for VCD). The definition of $s_0$, $s_1$, $P_j$, $P_{jk}$ and $P_{kj}$ depends on the property of interest, which can be found in details in (Bloino and Barone, 2012; Bloino et al., 2015). i is the excited mode in the final state (transition from the ground state $|0\rangle$ to $|1_i\rangle$).

The first term (A) is similar to that found in the energy expression, with the potential presence of Fermi resonances. Strategies like those discussed in the previous section can in principle be employed. However, because the numerator $\left(s_0{\text{f}}_{ikl}{\text{f}}_{jkl}{\mathbf{P}}_j/16\right)$ and the quantity of interest (intensity instead of energy) are different here, the test proposed by Martin and coworkers is not fully adapted. Hence, a fundamental (ex: $1_i$) and an overtone (ex: $2_j$) or a combination band (ex: $1_j+1_k$) could potentially be deemed not in resonance for the energy, but give incorrect contributions to the intensity. Nonetheless, the terms are quite similar so that a sufficiently low threshold can prevent or significantly mitigate this risk. The potential caveat is for nearly uncoupled states, so the corresponding cubic force constant (${\text{f}}_{ijk}$ here) could be nearly null, but in this case, it is unlikely that the property derivatives (${\mathbf{P}}_{jk}$ and ${\mathbf{P}}_{kj}$) will be large enough to compensate. Another, more delicate problem, is that the strategy proposed in DCPT2, which relies on the variational energy of a model system, like Martin’s test, cannot be applied here. Indeed, the ad hoc transformations used to compute the vibrational energies prevent the analytic definition of a transformed Hamiltonian to be used in the derivation of new formulas. For these reasons, DCPT2, and by extension HDCPT2, are ill-suited for the calculation of transition moments. Within a DVPT2 scheme, a list of states in resonances is established from the energy calculations and used to flag resonant terms in intensity calculations as well. The variational correction used in GVPT2 is done here by projecting the “deperturbed” results onto the variational states produced by the variational treatment, by applying the eigenvector matrix.

Cases B and C introduce a new form of resonances, between fundamentals. This can be related to the Darling–Dennison coupling in energy, which represent a correction to VPT2 shortcomings for the latter case, but can result in an incorrect account of the anharmonicity for intensities. A typical situation is the presence of one or more methyl groups, with the related C–H symmetric and antisymmetric stretching vibrations having very close energies. The treatment of DDRs is thus critical for intensity calculations, at least between fundamentals (1–1). Since Darling–Dennison couplings only appear in the variational treatment, they are absent in DVPT2, which makes it potentially unsuitable to deal with intensities without some improvement. To avoid confusion between the two formulations, the intensity-extended version, which performs exactly in the same way for the energies, will be labeled IDVPT2. Term C has the particularity of being sensitive to both Fermi and Darling–Dennison resonances. Robust strategies to identify DDRs in intensity calculations are still lacking in the literature, and may require the combination of multiple criteria. Indeed, the Darling–Dennison interaction terms being significantly different from the terms involved in the equation of transition moments, resonances between weakly coupled states could be ignored during the energy step. In practice, this problem usually arises for very close states. For this reason, a secondary test is added here on the magnitude of the DD term, by considering its magnitude weighted by the squared inverse of the energy difference between the potentially resonant states ($1/{\left({\omega }_i-{\omega }_j\right)}^2$ here). In case of resonance, a slightly different formulation must be used to properly discard all contributions, reported in Table 1 for the cases described here and in full in the Supplementary Material. The discarded terms are reintroduced through the projection described before onto the variational states. It should be noted that GVPT2 here refers by default to the treatment of both Fermi and Darling–Dennison resonances, so the non-resonant basis is IDVPT2. The transformations applied to each term (A, B, C) in the schemes described above are given in Table 1. Available schemes with a summary of their particularities are listed in Table 2.

To summarize this section, GVPT2 represents here the most complete and accurate method, provided FRs and DDRs are correctly identified. However, the resulting states may depart significantly from the original basis of harmonic states, making band assignment more complex. This problem does not exist for DVPT2 and its variant IDVPT2, but the truncated treatment can limit the overall quality of the results. If only energies are of interest, for instance in thermodynamics and reaction kinetics, HDCPT2/GHDCPT2 can be appealing for the absence of selection criteria.

As a final remark, it should be noted that the schemes presented here cannot overcome inherent limitations of VPT2, which derive from both the limited Taylor expansion and low-order perturbative formulation. Large amplitude motions (LAMs) like torsions and some out-of-plane vibrations can be poorly described by quartic force fields especially when employing Cartesian-based normal modes, and require extensive samplings of the PES, coupled with ad hoc Hamiltonians. Since the aim of this work is to assess the quality achievable by the combination of VPT2 and DFT, only semi-rigid systems devoid of LAMs are considered in the following. Mitigation strategies are discussed in conclusion.

3 Computational Details

Very tight criteria were used for geometry optimization, which means that the convergence was reached when the maximum forces and displacements were smaller than $5\times {10}^{-6}$ Hartrees/Bohr and $2\times {10}^{-5}$ Å, respectively. Then, analytic harmonic frequencies were computed to ensure that a true energy minimum was reached.

Anharmonic data were built by numerical differentiation of analytic second-derivatives of the potential energy and of first derivatives of the properties of interest, using a step of 0.01 $\sqrt{\text{amu}}\cdot \text{bohr}$ along the mass-weighted normal coordinates (Barone, 2005).

The cost of generating the anharmonic force field can be prohibitive for some electronic structure calculation methods, starting with double-hybrid DFT functionals. Conversely, if the anharmonic correction is small compared to the harmonic reference, the need for accuracy becomes less critical. As a consequence, a reasonable mitigation technique would be to compute the harmonic terms in the vibrational Harmiltonian (Eq. 5) and in the Taylor expansion of the properties of interest at the highest possible level of theory, and the anharmonic terms with a more cost-effective alternative. Studies in the literature have already pointed out the ability of such hybrid schemes to producing very accurate results (Begue et al., 2005; Carbonniere et al., 2005; Puzzarini and Barone, 2008; Biczysko et al., 2018; Puzzarini et al., 2019). While very convenient, these models involve some technical challenges related to the necessity of combining two levels of calculation, with their respective equilibrium structures and normal coordinates. In the best conditions, the latter are the same, making the definition of the hybrid scheme simple, since a single basis of normal modes need to be considered. In practice, subtle or even larger differences appear, in particular for vibrations, which are more sensitive than geometrical parameters. Two approaches can be chosen, depending on the reference geometries and displacement vectors for the numerical differentiation. The first one is to use only the normal coordinates from the higher level, computing the force constants and first derivatives of properties with respect to the Cartesian coordinates at each displaced geometry at the lower level. Otherwise, the anharmonic force field and property surfaces are built using only reference data from the lower level. The second technique has the advantage of making any combination simple to test, harmonic data from one level are combined with pure anharmonic data from another one. The problem is that the two sets of normal coordinates must be very close, and the geometries very similar, especially to combine properties. The consistency between the two methods can be assessed by computing the Duschinsky transformation to pass from one set of normal coordinates to the other (see Bloino et al. (2016) for details). This problem does not exist with the first method. However, if the two levels are significantly different, some approximation is still present if the second numerical differentiation is only diagonal, which is typically the case with standard schemes involving only two displacements ($+\delta Q$ and $-\delta Q$) along each coordinate. In the present case, the chosen levels in the hybrid schemes were close enough so that both techniques converged to the same results (overlap between the sets of normal modes above 90% based on the Duschinsky transformation, as described in Ref (Puzzarini et al., 2019)). For this reason, the second, more flexible approach was chosen.

For the VPT2 schemes, the following parameters were used. For Fermi resonances (DVPT2/GVPT2), states distant by less than 200 cm⁻¹ in energy, and with an associated value for Martin’s test greater than 1 cm⁻¹ were flagged as resonant. For Darling–Dennison resonances (IDVPT2/GVPT2/GHDCPT2), the following criteria were used:

• 1–1: Frequency difference lower than 100 cm⁻¹, minimum value of the coupling term: 10 cm⁻¹, and weighted coupling greater than 1 cm⁻¹.

• 2–2: Frequency difference lower than 100 cm⁻¹, and minimum value of the coupling term: 10 cm⁻¹.

For HDCPT2/GHDCPT2, the parameter for the $\text{Λ}$ function were the default ones in Gaussian, $\alpha =1.0$, $\beta =0.5\times {10}^6$.

Optimization and harmonic frequency calculations were done with Gaussian 16 (Frisch et al., 2016) and anharmonic calculations with a development version of Gaussian (Frisch et al., 2020).

4 Results and Discussions

4.1 Phase I: Electronic Structure Calculation Methods

4.1.1 DFT Functionals

To identify the best performing functional, the group of representative small molecules sketched in Figure 1 (top panel, G1) was selected. They were chosen so that the Fermi resonances were minimal or even absent to prevent any ambiguity in their identification. In these conditions, as noted in the theoretical details, GVPT2 is expected to perform the best, and was used as reference. Building upon previous analyses (Bloino et al., 2016; Barone et al., 2020), a shortlist of well-performing hybrid–B3LYP, B3PW91 (Becke, 1993), APF (Austin et al., 2012)–and double-hybrid–B2PLYP (Grimme, 2006) and revDSD-PBEP86-D3(BJ) (Santra et al., 2019)–exchange-correlation functionals was built. The ensemble of functionals/basis sets combinations and their respective abbreviations can be found in Supplementary Table S1 of the Supplementary Material. Let us first consider the influence of the functional, adopting the calendar-based basis set jun-cc-pVTZ as reference (Papajak et al., 2011), since it has shown very good performance in previous studies (Barone et al., 2020). It should also be noted that a basis set of at least triple-ζ quality is needed to reach converged results with double hybrid functionals. The impact of the empirical dispersion using Grimme’s DFT-D3 scheme (Grimme et al., 2010) with Becke-Johnson damping (Grimme et al., 2011) (noted D3(BJ)) was also evaluated. APF has its own empirical-dispersion scheme, which could not be used here due to technical issues. The mean (MAE) and maximum (|MAX|) absolute errors for the fundamental energies at both harmonic (ω) and anharmonic (GVPT2, ν) levels are reported in Figure 2. A first observation is that the empirical dispersion on such small systems is negligible. In order to make the effects of the dispersion more apparent, the differences between the errors obtained for each functional with and without empirical dispersion are plotted in Supplementary Figures S1–S3 in Supplementary Material. In this particular case, D3(BJ) doesn’t show particular advantages for these small molecules and can slightly worsen the results, so it was not considered in the following. The influence of empirical dispersion on the prediction of vibrational spectra for larger-size molecules and complexes will be deferred to a separate work. For hybrid functionals, we observe a decrease by two thirds of the maximum error between the harmonic and anharmonic fundamental energies. The gain is lower for B3LYP, which is slightly closer to experiment at the harmonic level and less scattered (|MAX| = 150 cm⁻¹). However, one should be careful in comparing harmonic computations with experimental data, since the latter refer actually to non-harmonic vibrations. Hence, an overestimation of the transition energies is expected. All functionals chosen here display similar trends at the anharmonic level, with an average error of about 13 cm⁻¹ and a maximum error ranging from 49 to 76 cm⁻¹. PW6B95 is performing slightly worse than the others, with the highest error, close to APF which, however, has an average unsigned error close to B3PW91 (13 vs 12 cm⁻¹), the best at this level. Considering now the double-hybrid functionals, a first observation is that they are not significantly different at the harmonic level, having even a higher average error than B3LYP and B3PW91. Adding the anharmonic correction leads to significantly better results with an average error within 10 cm⁻¹, and a maximum error below 20 cm⁻¹ for revDSD-PBEP86-D3(BJ) (25 for B2PLYP). This confirms the very good performance of double hybrid functionals already noted in the literature (Biczysko et al., 2010) for vibrational energies.

FIGURE 1

FIGURE 1. Structures of the molecules in this benchmark study. (A): Water, (B): H₂S, (C): formaldehyde, (D): thioformaldehyde, (E): ethylene, (F): 1,1-difluoroethylene, (G): fluoroethylene, (H): 1-fluoro-1-chloro-ethylene, (I): furan, (J): methyloxirane, (K): naphthalene, (L): 1-fluorobenzene, (M): Glycine.

FIGURE 2

FIGURE 2. Comparison of the mean (MAE) and maximum absolute (|MAX|) errors for different functionals with the jun-cc-pVTZ (JNTZ) basis set over the first set of molecules (G1): APF: APF, B3L: B3LYP, B3P: B3PW91, PW6: PW6B95, B2P: B2PLYP, RDS: revDSD-PBEP86-D3(BJ). ω (blue) corresponds to the harmonic approximation, ν (pink/red) to GVPT2. The values are indicated on top of the bars. Reference experimental data taken from: water (Tennyson et al., 2001), H₂S (Isoniemi et al., 1999), formaldehyde (Reisner et al., 1984), thioformaldehyde (Suzuki et al., 2007), ethylene (Nakanaga et al., 1979), 1,1-difluoroethylene (Krasnoshchekov et al., 2013), fluoroethylene (Tasinato et al., 2006), 1-fluoro-1-chloro-ethylene (Pietropolli Charmet et al., 2016).

4.1.2 Basis Sets

Let us now focus on the basis sets. Following the previous results, B3PW91 was chosen as model for the hybrid functionals, and revDSD-PBEP86-D3(BJ) for double hybrids. For B3PW91, three Pople-based basis sets were chosen, two double-ζ, 6–31+G(d) and SNSD, a basis set developed in our group with emphasis on the description of spectroscopic observables (SNSDT, 2021), and a large triple-ζ basis set, 6–311++G(3df,3pd). The SNSD basis set is obtained by adding to 6–31G(d,p) the diffuse functions of the aug-cc-pVDZ basis set and, for non-hydrogen atoms, a very tight s function. It is noteworthy that Cartesian, (e.g., 6d) functions were used in the original implementation and must be retained because, as is well known, this choice effectively increases also the space of s functions (details on the basis sets for the first 2 rows can be found in the Supplementary Material). Variations of the Dunning basis sets, of double-ζ (jul-cc-pVDZ) and triple-ζ (jun-cc-pVTZ and spaug-cc-pVTZ) quality were considered. Furthermore, the polarization consistent double- (pc-1) and triple-ζ (pc-2) basis sets proposed by Jensen (Jensen, 2001), as well as their augmented versions with diffuse functions, aug-pc-1 and aug-pc-2, were used. Some modifications were also applied to the latter two for this work. For aug-pc-1, the new version labeled naug-pc-1 had all d diffuse functions removed for all atoms present in the molecules of interest in this work. With the same aim of reducing the computational cost, all d and f diffuse functions were removed from the aug-pc-2 basis set, leading to the new naug-pc-2 set. For instance, in the case of 1-fluoro-1-chloroethylene, this reduces the number of basis functions from 114 to 94 (pc-1) and to 234 to 176 (pc-2). A decrease of about 31% (double-ζ) and 62% (triple-ζ) on the total computational cost needed to generate the anharmonic constants is observed on two octo-core processors machines (2.60 GHz). Provided that the overall accuracy is not very different, going from aug-pc to naug-pc basis sets represents an appealing strategy to tackle larger molecular systems. Finally, Ahlrichs and coworkers’ def2-TZVP basis set (Weigend and Ahlrichs, 2005) was also included. For revDSD-PBEP86-D3(BJ), only basis sets of at least triple-ζ quality were considered. In addition to spaug-cc-pVTZ, jun-cc-pVTZ and def2-TZVP already selected for B3PW91, Dunning’s aug-cc-pVTZ (Kendall et al., 1992) was used. Quadruple-ζ basis sets become quickly prohibitive when tackling larger molecular systems. For this reason, smaller and thus more affordable alternatives were explored, by combining s and p basis functions from cc-pVQZ with d and f basis functions from cc-pVTZ (Dunning, 1989). Recent investigations (Kruse et al., 2020) have shown that basis sets optimized for explicitly correlated methods do a very good job also in conventional computations. The main distinctive feature of these basis sets is to use an (n+1)ζ set for the sp space in conjunction with an nζ set for polarization functions. We have therefore tested also this family of basis sets. A version without diffuse functions was labeled M-cc-pV[T+Q]Z, and an extension with diffuse functions on all atoms except H was labeled M-aug-cc-pV[T+Q]Z. The errors on the first test suite with respect to the reference data are shown in Figure 3. Let us start with B3PW91. The harmonic results are relatively stable with all basis sets, with an average unsigned error ranging from 50 to 64 cm⁻¹, for a maximum absolute error of 141–218 cm⁻¹. With anharmonic contributions, a significant improvement in the mean unsigned error is observed for all basis sets, with values between 11 for SNSD and 16 cm⁻¹. If we look at the maximum error, however, the results are more disparate, some performing well, others exhibiting very high errors. Overall, SNSD gives the best average energies, but def2-TZVP has a lower maximum deviation, of 54 cm⁻¹. Basis sets of triple-ζ quality give very close results, with MAE of about 12–13 cm⁻¹ and a maximum error of about 60 cm⁻¹. The addition of diffuse functions tends to improve the results for double-ζ basis sets, as evidenced by comparison between pc-1 and aug-pc-1 results (|MAX| from 78 to 68 cm⁻¹), as well as 6–31+G(d) with SNSD (|MAX| from 70 to 59 cm⁻¹). Regarding the Pople basis sets, 6–31+G(d) seems to be clearly insufficient, in line with previous observations (Mitra and Roy, 2020). Since the main difference between 6–31+G(d) and SNSD basis sets is the presence in the latter of p functions on hydrogen and diffuse d functions on heavy atoms, our results confirm the importance of these contributions. For pc-1, the improvement on the maximum error with the addition of diffuse functions comes at the expense of a larger average error. It is noteworthy that the new naug-pc-1 basis set, despite being smaller than aug-pc-1, gives results substantially equal to the latter one. The maximum error of pc-1 and its related basis sets remains significantly higher than the others, which make them potentially ill-suited for spectroscopic applications. The improvement obtained with diffuse functions becomes less apparent at the triple-ζ level, with pc-2 vs aug-pc-2 being near equal in terms of performance. Here too, naug-pc-2 provides the same quality as aug-pc-2 for a significantly lower cost. The largest average error is found for jul-cc-pVDZ, with an average deviation of about 30% compared to the standard performance of the basis sets chosen here. The marginally higher maximum error hints at a more systematic issue, which could be investigated in the future on larger systems. The same procedure was applied to B3LYP for comparison purposes. The results can be found in the Supplementary Material (Supplementary Figure S6), with similar trends and conclusion.

FIGURE 3

FIGURE 3. Comparison of the mean (MAE) and maximum absolute (|MAX|) errors at the B3PW91 (B3P) and revDSD-PBEP86-D3(BJ) (RDS) levels with different basis sets over the first set of molecules (G1). G31: 6–31+G(d), SNSD: SNSD, G311: 6–311++G (3df, 3pd), PC1: PC-1, APC1: aug-PC-1, NPC1^a: naug-PC-1, JLDZ: jul-cc-pVDZ, DTP: def2-TZVP, PC2: PC-2, APC2: aug-PC-2, NPC2^b: naug-PC-2, SPTZ: spaug-cc-pVTZ, JNTZ: jun-cc-pVTZ, ATZ: aug-cc-pVTZ, MZ^c: M-cc-pV[T+Q]Z, MAZ^d: M-aug-cc-pV[T+Q]Z. ω (blue) corresponds to the harmonic approximation, ν (pink/red) to GVPT2. The values are indicated on top of the bars. See Figure 2 for details on the reference data. (a): basis set built from aug-PC-1 Jensen (2001), by removing the d diffuse functions. (b): basis set built from aug-PC-2 Jensen (2001), by removing the d and f diffuse functions. (c): basis set built by combining s and p basis functions from cc-pVQZ, and d and f from cc-pVTZ. (d): basis set built by combining s and p basis functions from cc-pVQZ, and d and f from cc-pVTZ, augmented with d diffuse functions on all atoms except for H.

Switching to revDSD-PBEP86-D3(BJ), all chosen basis sets exhibit very close behaviors at the harmonic level, with mean unsigned errors between 58 and 63 cm⁻¹, except for jul-cc-pVDZ, which performed slightly better here. At the GVPT2 level, these errors go down to 8–11 cm⁻¹ with maximum errors typically situated between 18 and 21 cm⁻¹, a third of what was found with B3PW91. The only exceptions are for the double-ζ basis sets, SNSD and jul-cc-pVDZ, with the largest error at 34 and 38 cm⁻¹, respectively, and aug-cc-pVTZ, which has the highest maximum absolute error (|MAX| = 72 cm⁻¹). This larger discrepancy is to be ascribed to the out-of-plane twisting motion of ethylene (mode 4, at about 1050 cm⁻¹ at the harmonic level). The lower accuracy is not related to any variational correction, but to the anharmonic force field itself. The new basis set, M-aug-cc-pV[T+Q]Z does not suffer from this problem, with results indistinguishable from jun-cc-pVTZ and spaug-cc-pVTZ. To confirm these trends, a parallel study was done with B2PLYP as reference, showing relatively close behaviors (Supplementary Figure S6 in the Supplementary Material). To conclude this part, jun-cc-pVTZ provides very good results with double hybrid functionals for the systems of interest here, while for hybrid functionals, double-ζ basis sets augmented with diffuse functions are sufficient.

4.1.3 Hybrid Schemes

As mentioned above, the cost of double hybrid functionals with the associated basis sets can become prohibitive for large molecules, so that more cost-effective approaches are needed, possibly based on the use of hybrid schemes, which allow the computation of the anharmonic contributions at a lower level of theory. The performance of such models is illustrated in Figure 4, where the errors on the fundamental energies of the first group of molecules with respect to experiment are compared with the best data obtained with a single electronic structure calculation method. More explicitly, B3LYP/SNSD, B3PW91/SNSD and B3PW91/jul-cc-pVDZ for hybrid, and B2PLYP/jun-cc-pVTZ and revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ for double hybrid functionals were chosen. With revDSD-PBEP86-D3(BJ) as the harmonic reference, the improvement is apparent, with a mean unsigned error of 9 cm⁻¹, on par with the full, more expensive double-hybrid calculations. revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ in conjunction with B3PW91/SNSD has also a maximum error almost on par with the full reference (23 vs. 20 cm⁻¹), while the other combinations performs a bit worse, at 47 cm⁻¹. This result is in contrast with what was observed with the pure hybrid functionals, where B3PW91/SNSD performed worse in terms of largest error. The most likely explanation is an error in the harmonic reference, well compensated by revDSD-PBEP86-D3(BJ) while the anharmonic force field is of very good quality. An excellent agreement is reached between the hybrid schemes and the full double hybrid anharmonic calculations, with the best performing combination being clearly revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD.

FIGURE 4

FIGURE 4. Comparison of the mean (MAE) and maximum absolute (|MAX|) errors of pure functionals and hybrid models over the first set of molecules (G1). ω (blue) corresponds to the harmonic approximation, ν (pink/red) to GVPT2. The values are indicated on top of the bars. The hybrid models are built by combining harmonic calculations with double hybrid functionals, with anharmonic force fields at a lower level, and noted “harmonic//anharmonic”. See Figures 2, 3 for details on the abbreviations and the reference data used.

4.2 Phase II: VPT2 Schemes

These preliminary tests have highlighted in some situations the sensitivity of the anharmonic correction. In the second phase of this study, the influence of the VPT2 schemes described in the theoretical details, first on the fundamental energies, and next on the intensities, will be analyzed. To be more representative of potential resonance patterns, larger molecular systems were chosen, starting from the derivatives of ethylene already included in G1 since they have simple resonance patterns. Another criterion for the selection of the molecules was the potential availability of reliable experimental data for the transition intensities in the literature. The final list is displayed in the middle panel (G2) of Figure 1. Based on their performance in the first part of this benchmark, three levels of theory were selected to proceed, B3PW91/SNSD, revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ and the hybrid scheme combining them, revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD. For the sake of readability, shortened labels will be used in the discussion, respectively B3PW91, REVDSD and HYBRID.

4.2.1 Ethylene and Halogenated Derivatives

Starting with the vibrational energies of ethylene (Table 3), the VPT2 fundamental energies are indeed fairly good for REVDSD and HYBRID, worsening with B3PW91. The largest error is to be ascribed to mode 9 (CH symmetric stretching), in resonance with the combination band 7 + 8. The energy gap between the resonant states with B3PW91 is half that of REVDSD and HYBRID (22.7 vs. 11.4 cm⁻¹) while the coupling term, the cubic force constant involving the three modes, remains unchanged, which leads to a larger contribution to the anharmonic $\chi$ matrix. Removing the resonant terms (DVPT2) reduces the difference between the two methods, with a mean unsigned error nearly equal (12–13 cm⁻¹) and the maximum absolute error remaining slightly higher for B3PW91. HYBRID and REVDSD give very close results. Regarding the variational correction (GVPT2), a further improvement in the maximum error, by about one third, and in the average error, from 12 to 10 cm⁻¹, is observed for REVDSD and HYBRID, while no changes are visible for B3PW91. Since the only significant Darling–Dennison terms appear in the CH stretching first overtone region, their inclusion does not lead to any notable changes (GVPT2 vs GVPT2^a, see Table for details). Regarding DCPT2 and its variants, GHDCPT2 is not expected to have any impact over HDCPT2 for the reasons mentioned for GVPT2, that is the lack of Darling–Dennison terms in the energy range of interest here. For B3PW91, the errors of DCPT2 and HDCPT2 are similar to DVPT2/GVPT2, while for REVDSD, GVPT2 performs a bit better. It is noteworthy that, for all methods, DCPT2 and HDCPT2 recover very well the vibrational energy of mode 9, the larger error being on mode 8 (C=C stretching) for all the DCPT2 variants. This mode is an interesting case of the problem of criteria for the identification of Fermi resonances. Indeed, for B3PW91, the transition energy to $|1_8\rangle$ remains unchanged between VPT2, DVPT2 and GVPT2, showing that at this level, it is not involved in any resonance. Conversely, for REVDSD and HYBRID, it decreases from VPT2 to DVPT2, getting closer to the experimental value. The variational correction further improves the result. This disparity of behaviors is caused by the identification or not of a Fermi resonance between $|1_8\rangle$ and $|2_1\rangle$. Investigating the output from the resonance test, the energy difference between these states at the REVDSD/HYBRID and B3PW91 levels are respectively 27.3 and 53.0 cm⁻¹, while the corresponding cubic force constant is roughly the same, around 53 cm⁻¹. Because of this, the value for Martin’s test is 1.5 cm⁻¹ for HYBRID (1.6 for REVDSD), but lower than the threshold for B3PW91, so the two states are not flagged in resonance at this level. By contrast, DCPT2 and HDCPT2 remain remarkably constant from one method to the other, the transition energy varying from 1657 to 1654 cm⁻¹ from B3PW91 to HYBRID, compared to a reduction from 1657 to 1642 cm⁻¹ for GVPT2, which explains the higher maximum error of the former schemes. As noted before (Bloino et al., 2012), this consistency makes DCPT2 and, especially, HDCPT2 valuable when studying the quality of the anharmonic correction for different electronic structure methods, since the risk of sudden changes between methods can be caused by different resonance patterns despite low variations in the force field. Nevertheless, the overall better quality of GVPT2 confirms its choice in the first phase. Another observation is that the behavior and performance of HYBRID very often matches REVDSD. For this reason, only the former will be mentioned in the following, except where they diverge. For the IR intensities, only the experimental bands, which could be unambiguously assigned to a single, fundamental transition have been retained (see Table 3). VPT2 applied to the hybrid scheme already gives qualitatively good results, further refined with GVPT2. In absence of Darling–Dennison couplings, IDVPT2 is not different from DVPT2, which explains their identical trends. In contrast, VPT2 with B3PW91 gives very large errors, related again to mode 9, which is closer in energy to the combination band 7 + 8 compared to REVDSD at the harmonic level. The result is a transition intensity more than 40 times larger than the experimental one. A qualitative agreement is reached by removing the resonances in DVPT2, and the intensity is further improved by application of the variational correction.

TABLE 3

TABLE 3. Harmonic and anharmonic vibrational energies (in cm⁻¹) and intensities (km/mol) of ethylene computed at the revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ, revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ and hybrid revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels of theory and with different VPT2 schemes, VPT2 (V), DVPT2 (DV), IDVPT2 (IDV), GVPT2 (GV), DCPT2 (DC), HDCPT2 (HDC), GHDCPT2 (GHDC). The maximum (MAX|) and mean (MAE) absolute errors, as well as the maximum relative error weighted by the reference value (|REL|, given in %) are reported with respect to experimental data from (Nakanaga et al., 1979).

Another interesting example is 1,1-difluoroethylene, whose results are collected in Table 4 (Supplementary Table S2 of the Supplementary Material for REVDSD). At the VPT2 level, this time, HYBRID gives a larger unsigned error compared to experiment than B3PW91, a behavior opposite to what was observed for ethylene. Upon investigation of the individual transition energies, similar patterns appear. The error is concentrated on one mode, 11, a CH symmetric stretching vibration. This mode is involved in a Fermi resonance with the 9 + 10 combination band. Between the two methods, HYBRID and B3PW91, the value of the cubic force constant between these modes is about the same, but the energy difference, ${\omega }_{11}-\left({\omega }_9+{\omega }_{10}\right)$ is halved for the former with respect to the latter (15.9 vs 29.4 cm⁻¹). This difference is mostly responsible for the changes observed at the VPT2 level. Again, a correct treatment of resonances lead to very close results between the three methods, with a maximum error of 30–31 cm⁻¹ and an average of 9–10 cm⁻¹. Compared to ethylene, the improvement with GVPT2 is even higher, and the accuracy reached by the hybrid scheme very good (|MAX| = 17 cm⁻¹, MAE = 7 cm⁻¹). These data highlight the quality reachable by GVPT2. At variance, DCPT2 and HDCPT2 show the worst performance after VPT2, with a maximum error of 53 cm⁻¹ for HYBRID. The cause here is again mode 11, for which the mitigation schemes used in DCPT2 and HDCPT2 give an underestimation of the anharmonic effect, with the transition energy to the associated fundamental state at 3110 cm⁻¹ with HYBRID to be compared to an experimental value of 3058 cm⁻¹. Some interesting observations can be made for mode 9 (CH₂ scissor) and 10 (C=C stretching). The former mode is indeed involved in a Fermi resonance at all levels, but its treatment leads to divergent results. For B3PW91, removal of the resonant terms increases the transition energy from 1341 (VPT2) to 1358 cm⁻¹ (DVPT2), making it very close to experiment (1359 cm⁻¹). Inclusion of the variational correction worsens the result, making it by chance very close to VPT2. For HYBRID, the trend is similar, but the starting VPT2 value is higher, at 1365 cm⁻¹. Exclusion of the resonance leads to a significant overestimation of the energy of the transition, to 1385 cm⁻¹, which is then lowered at 1366 cm⁻¹ at the GVPT2 level. DCPT2/HDCPT2 values tend to slightly overestimate the energy compared to GVPT2 by about 5 cm⁻¹. Mode 10 shows a pattern similar to that observed for mode 8 in ethylene, which is in part expected due to the similarities between the two vibrations. Here too, Fermi resonances are only detected at the HYBRID level, and ignored with B3PW91. However, the incidence of the schemes is lower than for ethylene, with variations smaller than 5 cm⁻¹ (13 cm⁻¹ before).

TABLE 4

TABLE 4. Harmonic and anharmonic vibrational energies (in cm⁻¹) and intensities (km/mol) of 1,1-difluoroethylene computed at the B3PW91/SNSD and hybrid revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels of theory and with different VPT2 schemes. The maximum (|MAX|) and mean (MAE) absolute errors, as well as the maximum relative error weighted by the reference value (|REL|, given in %) are reported with respect to experimental data from (Tasinato et al., 2006).

The fundamental energies of the last two ethylene derivatives, 1-fluoroethylene and 1-fluoro-1-chloroethylene, are listed in Supplementary Tables S3, S4 of the Supplementary Material, respectively. For the former molecule, a very good agreement is observed with all schemes, with patterns similar to ethylene and 1,1-difluoroethylene. An interesting aspect here is the influence of Darling–Dennison couplings between the fundamentals of modes 10 and 11, both related to CH stretchings. The impact can be assessed by comparing GVPT2 with GVPT2^a, the latter having DDRs excluded, and GHDCPT2 with HDCPT2. For GVPT2, inclusion of these couplings increases the energy gap between the two states, lowering $|1_{10}\rangle$ and raising $|1_{11}\rangle$. For HYBRID, this improves the agreement of the former with experiment, but worsens the other one. For B3PW91, the theoretical energies match very closely their experimental counterparts, with a difference of 7 and 2 cm⁻¹ for $|1_{10}\rangle$ and $|1_{11}\rangle$, respectively. As a final comment on GVPT2, let us mention that $|1_{10}\rangle$ is also involved in a Fermi resonance with $|1_8{1}_9\rangle$, resulting in a polyad involving three states for the variational correction. If the coupling is not negligible, which is the case here—38 cm⁻¹ for the Fermi resonance, 16 cm⁻¹ for the Darling-Dennison coupling with the HYBRID scheme, and similar values for the other levels—non-negligible mixing between these states occur. As a result, the variational states do not overlap fully with the original DVPT2 states, and a “similarity rule” was applied, choosing the states with the largest overlap. While using an oscillator-based notation is more intuitive, when doing so, the description of the associated vibration for the assignment of the bands should be carefully checked. An effective vibration can be derived by computing the coefficients of the perturbed wave function (Puzzarini et al., 2019). The same splitting is observed for GHDCPT2, with transition energies being slightly better compared to the GVPT2 ones. Overall, the average error is low, with a maximum error related to modes not involved in resonances, mode 5 (H-C-H out-of-plane wagging) for REVDSD and 9 (C=C stretching) for B3PW91. It is noteworthy that for this system, the results with the hybrid scheme are slightly worse than their B3PW91 counterparts (MAE = 11 cm⁻¹ with GVPT2, compared to 7 cm⁻¹ for B3PW91). The cause is the superposition of the largest errors from REVDSD and B3PW91 in this case. For the intensities, a comparison with the experimental data for which the assignment was unambiguous shows a very good agreement. The transition intensities for $|1_1\rangle$ is perfectly reproduced, and the relative errors is of about 20% otherwise. While the B3PW91 energies were very good, the predicted intensities are not so close to experiment, and the hybrid scheme works as expected, bringing the values closer to experiment, with relative errors similar to those of REVDSD. Since modes 10 and 11 are involved in a Darling–Dennison resonance, some differences would be expected between DVPT2 and IDVPT2. However, the corresponding transition intensities are low and only marginal variations can be seen. Let us conclude this analysis of ethylene and its derivatives with 1-fluoro-1-chloroethylene. A number of observations made before can be confirmed here as well. The error with B3PW91 is again higher, and the hybrid scheme can significantly recover the quality of the anharmonic calculations done at a higher and more expensive level of theory. In analogy with 1-fluoroethylene, inclusion of Darling–Dennison interaction terms in the variational correction step improves the energies of the C-H stretching modes. A notable difference with previous cases is that DCPT2 and HDCPT2 perform better, especially for REVDSD and the hybrid scheme, where the maximum error is about one third lower than GVPT2. The largest GVPT2 error for REVDSD corresponds to mode 11 (CH symmetric stretching), where a correction pattern similar to mode 10 of 1,1-difluoroethylene is observed; VPT2 overestimates the energy, which is then lowered with DVPT2, but reverted back by application of the variational correction. Conversely, DCPT2/HDCPT2 provides a good estimate of the actual energy. It is noteworthy that, despite the larger maximum error, GVPT2 gives a low average error (7 cm⁻¹, the same as DCPT2/HDCPT2), confirming the overall quality of the scheme. The larger error for B3PW91 stems from another fundamental state, $|1_{10}\rangle$ (CH stretching), which is significantly overestimated by all schemes. This worse performance of B3PW91 is to be ascribed in part to the anharmonic force field since this state shows the largest error also for the hybrid model, decreasing from 37 to 23 cm⁻¹. The intensities are generally well reproduced, especially for the more intense bands. A notable exception is the intensity of $|1_{10}\rangle$ at the VPT2 level, which is almost null (2 km/mol for HYBRID, 9 for REVDSD), compared to an experimental estimate of 130 km/mol. This problem is less serious with B3PW91 and is likely related to the smaller energy gap between the resonant states $|1_{10}\rangle$ and $|2_6\rangle$ (30 cm⁻¹ for B3PW91, 12 cm⁻¹ for REVDSD/HYBRID). The effect is not so sharp for the energy, being partly hidden by the systematic overestimation of the band position. Removal of the resonance helps recovering the correct intensity in this case, while it causes opposite effects on energy for B3PW91 and HYBRID, lowering it for the former while increasing it for the latter. This difference of behaviors and importance in the error emphasizes the necessity for comprehensive resonance schemes, which are not only aimed at the energy, but also consider the potential impact on the intensity, as the latter can be far more sensitive to them. To conclude on the intensity of the resonant mode, once the variational correction is applied (GVPT2), the HYBRID intensity is significantly lowered, whereas this is not the case for the pure REVDSD. This behavior hints at some potential deficiencies of B3PW91/SNSD in describing higher-order derivatives of the electric dipole.

4.2.2 Furan and Methyloxirane

For furan, shown in Table 5 (Supplementary Table S5 of the Supplementary Material for REVDSD), the predicted energies are very good at all levels, with the hybrid scheme improving over B3PW91 to lower the average error at 4 cm⁻¹ with a maximum absolute error of 9 cm⁻¹, compared to 5 and 11 cm⁻¹, respectively. Because of the $C_{2v}$ symmetry of the molecule, only one Fermi resonance is present, involving the fundamental of mode 17 (in-plane ring deformation), with very mild effects due to the low coupling, the corresponding vibrational energy varying within 5 cm⁻¹ across the VPT2 schemes.

TABLE 5

TABLE 5. Harmonic and anharmonic vibrational energies (in cm⁻¹) and intensities (km/mol) of furan computed at the B3PW91/SNSD and hybrid revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels of theory and with different VPT2 schemes. The maximum (|MAX|) and mean (MAE) absolute errors, as well as the maximum relative error weighted by the reference value (|REL|, given in %) are reported with respect to experimental data from (Mellouki et al., 2001).

Methyloxirane has a comparable dimension, but $C_1$ symmetry. Because of this and the presence of a methyl group, which offers different pattern of CH stretchings, it is expected to be a more challenging case. For reasons that will be clearer later, the vibrations of the system are divided into two groups based on their energy, below 2000 cm⁻¹, which covers the so-called fingerprint region, and above this threshold, dominated by transitions related to the CH stretchings. The vibrational energies and IR intensities for the first 18 normal modes (out of 24) are reported in Table 6 (Supplementary Table S6 of the Supplementary Material for REVDSD). The error at the harmonic level is lower than for the previous examples, mainly because the higher-energy range, generally overestimated, is excluded. REVDSD tends to overestimate more systematically the transition energies compared to B3PW91, which is reflected by the higher average and maximum absolute errors (MAE = 27 cm⁻¹, |MAX| = 50 cm⁻¹, compared to 18 and 34). However, the behavior is not consistent and the distribution of errors differs between the two levels of electronic structure calculation methods. The largest error is thus found for the fundamental of mode 4 (CO asymmetric stretching) in B3PW91 and 17 (CH₃ asym. deformation) for REVDSD. These disparities have an important consequence on the resonances, including Darling–Dennison couplings, and can influence significantly their magnitude. As a matter of fact, only 1 Fermi resonance is detected in this list, but it involves the fundamental of mode 3 (CH₃ bending, with $|2_1\rangle$) for B3PW91, and mode 4 (with $|1_3{1}_2\rangle$) for HYBRID. For both, the effect is very small or even negligible, with a variation of 1 cm⁻¹ found between VPT2 and DVPT2 for B3PW91, compared to 7 cm⁻¹ for HYBRID. the consequence is that VPT2, DVPT2 and GVPT2 gives almost the same results. For HYBRID, a 1-1 Darling–Dennison coupling is also included, between the fundamental of modes 10 (CH₂ wagging) and 11 (CH bending). Comparison of the variational correction with (GVPT2) or without (GVPT2^a) this term shows a very minor effect, of 1 cm⁻¹, even if it is in the direction of a better agreement with experiment. Comparatively, DCPT2 and HDCPT2 perform a bit worse than DVPT2-based schemes, even if the average absolute error is about the same as the latter. The error stems from the fundamental of mode 18 (CH₂ scissor), not related to Fermi resonances with the standard test. However, the automatic transformation to non-resonant terms causes a divergence, that HDCPT2 fails to recover. Going to the higher-energy zone, in the CH-stretching region, an unequivocal equivalency between the variational states and VPT2 states is not possible anymore. Indeed, as shown in Supplementary Table S8 of the Supplementary Material, with HYBRID serving as reference, several fundamental states become interconnected inside the same resonance polyad. The diagonalization of the associated variational matrix leads to a strong mixing in some cases, shown in Supplementary Table S9 of the Supplementary Material, and associating a state to each fundamental becomes arbitrary. First of all, let us consider the VPT2 schemes for which this problem does no exist, with the fundamental energies and IR intensities reported in Supplementary Table S7 of the Supplementary Material. Because of the limited variational correction, GHDCPT2 does not suffer at such an extent of the difficulties highlighted for GVPT2, and was included as well. As expected, the error is distinctly higher, with a trend among the VPT2 schemes similar to what was observed in the lower-energy region. The largest VPT2 error in all cases is related to mode 20 (CH₂ symmetric stretching), and remains the largest source of errors even after removal of the Fermi resonances, but with a magnitude reduced by more than a half. The use of alternative schemes like DCPT2 and HDCPT2 does not improve the result, with an energy difference with respect to experiment of 54 cm⁻¹ for HYBRID, between DVPT2 (35 cm⁻¹) and VPT2 (82 cm⁻¹). Looking at the GVPT2 variational states, an alternative explanation arises, as the closest state in energy has a strong contribution from a combination band followed by an overtone, and is the most intense peak in the region (5 km/mol). For the hybrid scheme, an accidental degeneracy occurs at the VPT2 and DVPT2 levels between $|1_{22}\rangle$ (CH₃ asym. stretching) and $|1_{23}\rangle$ (CH stretch). The variational correction induces a splitting, with the most similar variational states having associated transition energies of 2995 and 3015 cm⁻¹, respectively, compared to 2995 and 3001 cm⁻¹ experimentally. To conclude this analysis, let us focus on the IR intensity, where a few phenomena occur. First, the VPT2 intensity of $|1_{19}\rangle$ (CH₃ sym. stretching) is almost null for REVDSD and HYBRID, and is recovered at the DVPT2 level, once Fermi resonances have been removed. Upon inspection of the FRs involving $|1_{19}\rangle$, the combination band $|1_{18}{1}_{17}\rangle$, within 10 cm⁻¹ of the former at the harmonic level, has a VPT2 energy of 2961 cm⁻¹ and an intensity of 12 km/mol. At the DVPT2 level, it becomes less than 1 km/mol, which means that the Fermi resonance causes a transfer of the intensity, with the fundamental one vanishing, while a connected combination or overtone becomes dominant in the surrounding. From a theoretical point of view, this is possible due to the relatively symmetric form of the potentially resonant terms between the transition moments of fundamentals and overtones/combinations, so a Fermi resonance will generally have an opposite effect on the corresponding moments (Bloino and Barone, 2012). This problem needs to be carefully considered when doing band assignment to avoid erroneous interpretations, and has been found to be critical when dealing with inorganic systems or large organic complexes (Fusè et al., 2019). The larger energy gap between $|1_{19}\rangle$ and $|1_{18}{1}_{17}\rangle$ at the B3PW91 level (22 cm⁻¹) is most likely the reason for which it does not happen at this level of theory. Another remarkable situation is the IR intensity for fundamental $|1_{23}\rangle$, which varies clearly between DVPT2 and IDVPT2 with all methods. The reason is that this state is involved in 1-1 Darling–Dennison resonances with two other fundamental states, $|1_{21}\rangle$ and $|1_{19}\rangle$ for B3PW91, $|1_{22}\rangle$ and $|1_{21}\rangle$ for REVDSD, and $|1_{22}\rangle$ and $|1_{19}\rangle$ for HYBRID, the variational terms for the latter being visible in Supplementary Table S9 of the Supplementary Material. All these are associated to CH stretching vibrations. The subtle changes in the resonance polyads illustrate the sensitivity of the variational treatment to the quality of the harmonic force field. The most striking example of the need to account for the 1-1 DDR in intensities is for $|1_{22}\rangle$ with REVDSD, which raises from 6 to 37 km/mol between the harmonic and VPT2 levels, and is unaffected by the removal of Fermi resonances. Once 1-1 DDRs are discarded, the intensity falls back to 3 km/mol, more in line with the harmonic values. Because of subtle variations in the anharmonic force field, such a steep reduction does not occur with the hybrid scheme when switching from DVPT2 to IDVPT2. While the transformation caused by the variational treatment can lead to redistributions of the energies and intensities, $|1_{22}\rangle$ is almost exclusively coupled with $|1_{23}\rangle$, of similar intensity, so not much difference is observed in the intensities of the final GVPT2 states. An important conclusion to draw from this example is that, while GVPT2 can provide accurate results, more complex resonance patterns, in the form of polyads, can occur with growing system size and structural complexity. The very concept of fundamentals, overtones and combination bands may prove to be unsuitable and source of misinterpretations.

TABLE 6

TABLE 6. Harmonic and anharmonic vibrational energies (in cm⁻¹) and intensities (km/mol) for the first 18 modes (by increasing energy) of methyloxirane computed at the B3PW91/SNSD and hybrid revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels of theory and with different VPT2 schemes. The maximum (|MAX|) and mean (MAE) absolute errors, as well as the maximum relative error weighted by the reference value (|REL|, given in %) are reported with respect to experimental data from (Polavarapu et al., 1985; Sunahori et al., 2010; Merten et al., 2013).

4.2.3 Naphthalene

As the final example of the second phase, naphthalene is the standard prototype of polycyclic aromatic hydrocarbons (PAHs), which have been extensively studied in recent years (Bloino, 2015; Mackie et al., 2015, 2016). This is also a model system for various applications, thanks to its rigidity and high symmetry (Swofford et al., 1976; Dierksen and Grimme, 2004; Basire et al., 2008; Pirali et al., 2009, 2013). A practical consequence of the latter is the possibility to check straightforwardly potential numerical instability or noise in the generation of the anharmonic force field or property surface. The vibrational energies and IR intensities are reported in Table 7. Because of the larger size, a full anharmonic force field at the revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ level becomes very expensive, so only B3PW91 and the hybrid scheme were used here. Since some vibrations are both IR and Raman inactive, some experimental data are missing and were omitted in the evaluation of the maximum and average absolute errors. Because of the larger number of normal modes, critical conditions with very close resonant states are more likely to occur, albeit mitigated by the symmetry conditions in the present case. Like furan, a good agreement is generally observed, with a mean absolute error of 8–10 cm⁻¹. The main exception is VPT2, which exhibit large errors, resulting in an average twice as large as DVPT2. For both B3PW91 and HYBRID, the largest deviation is related to mode 41 (CH stretching), involved in two Fermi resonances. One of them is with a combination band very close in energy, within 3 cm⁻¹ ($|1_{40}{1}_{36}\rangle$ for HYBRID and $|1_{40}{1}_{35}\rangle$ for B3PW91). The consequence of such a small denominator can be a large contribution from the relative term to the anharmonic correction. A second case occurs, mostly for the hybrid scheme, for mode 43 (CH stretching), where the associated fundamental energy is underestimated by 48 cm⁻¹. Removal of resonances through the DVPT2 scheme or their transformation with DCPT2 leads to more accurate energies. Actually, the DCPT2/HDCPT2 vibrational energies for these modes match exactly the experiment while the GVPT2 ones are overestimated for $|1_{41}\rangle$ by about 20 cm⁻¹ (3030 for HYBRID and 3028 for B3PW91, compared to 3008 cm⁻¹), and for $|1_{43}\rangle$ by 9 (HYBRID) and 6 cm⁻¹ (B3PW91). Even with a proper treatment of resonances, the maximum error remains large, and has its origin in the 600–800 cm⁻¹ region. It should be noted beforehand that, while the vibrations are the same between REVDSD and B3PW91, their relative energies are not, and swaps can be observed between couples of modes (12–13, 26–27, 32–33, 35–36). However, since the overlap is very high, building an hybrid scheme is not a problem. In the first approach, where the numerical differentiation is done along the higher-level normal coordinates, the correct transformation is done automatically. For the latter approach, used here, a reassignment of the normal modes must be done. For consistency with the previous analyses and to facilitate some discussions on the errors, the normal modes remain here ordered by increasing energy, irrespective of their irreducible representation. Unless the symmetry is incompatible with experimental observations, for instance with the Raman and IR inactive $A_u$ vibrations, the same order is chosen for the experimental data. Considering first B3PW91, the largest error is for mode 13 (out-of-plane skeletal distortion, $B_{2g}$) symmetry with a GVPT2 value of 720 cm⁻¹, compared to 773 cm⁻¹. Since it is not involved in any resonance or Darling–Dennison coupling, as confirmed by the VPT2 value, the large anharmonic correction of 77 cm⁻¹ is due to the force field. A direct consequence is that the order of the energy levels is altered, with $|1_{13}\rangle$ swapping position with the fundamental of mode 12 (skeletal breath, $A_g$). Both states are Raman active, and have close Raman intensities, which means that they can in principle be interchanged in the band assignment. By modifying the energy order, a maximum error of 44 cm⁻¹ is obtained. The same occurs for the hybrid scheme. The out-of-plane skeleton distortion of $B_{2g}$ symmetry is here the 12th mode in order of increasing energy, and it is lowered by all VPT2 schemes to be below the fundamental of mode 11 (CH bend, of symmetry $B_{1g}$), which is also Raman active. Reordering the energy levels reduces the distance between the peak found experimentally at 764 cm⁻¹ and its proposed assignment from 93–48 cm⁻¹. The largest source of error becomes the previous transition, at 726 cm⁻¹, with a theoretical transition 55 cm⁻¹ lower. The other major source of error in the mid-IR region is on $|1_{14}\rangle$ (out-of-plane CH bending, $B_{2g}$), with a deviation at the GVPT2 level compared to experiment of 39 cm⁻¹ for HYBRID and 26 for B3PW91. If we consider only the 20 IR-active vibrations, of $B_{1u}$, $B_{2u}$, and $B_{3u}$ symmetries, the maximum absolute error drops to 26 cm⁻¹ for B3PW91 and 13 for HYBRID with the GVPT2 scheme, and the mean absolute error to 7 and 5 cm⁻¹, respectively. The DCPT2 and HDCPT2 models perform slightly worse, with maximum errors of 33 cm⁻¹ and 35 cm⁻¹, respectively. The mean error is also slightly higher, at 9 and 7 cm⁻¹. This shows the very good agreement achievable with the hybrid scheme, for a low computational cost, while improvements are still needed to describe correctly the out-of-plane motions. A detailed comparison between theory and experiment (taken from Cané et al. (2007)) cannot be performed for intensities, since the values are very different. A noteworthy aspect is the large impact Fermi resonances can have on the intensities. The most striking example is with B3PW91/SNSD, where the transition intensity to fundamental $|1_{31}\rangle$ (in-plane CH bend) increases from 7 (harmonic) to 1278 km/mol, a value one order of magnitude higher than any other transition. With the DVPT2 scheme, the intensity lowers to a more plausible value of 5 km/mol. This problem is not seen to such an extent within the hybrid model, but some excessive anharmonic correction can be noted for mode 45 (CH stretching), where the transition intensities to the corresponding fundamental is 65 km/mol at the VPT2 level, while it is otherwise negligible. Despite the relatively large number of CH stretching modes, very few 1-1 resonances with non negligible impact on the intensity can be found by comparing the DVPT2 and IDVPT2 results. The most plausible reason is the high symmetry of the molecule, which imposes strong restrictions on the couplings between the modes.

TABLE 7

TABLE 7. Harmonic and anharmonic vibrational energies (in cm⁻¹) and intensities (km/mol) of naphthalene computed at the B3PW91/SNSD and hybrid revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels of theory and with different VPT2 schemes. The maximum (|MAX|) and mean (MAE) absolute errors, as well as the maximum relative error weighted by the reference value (|REL|, given in %) are reported with respect to experimental data from (Cané et al., 2007).

4.3 Phase III: Spectral Band-Shapes

To conclude this study, we directly apply the most accurate GVPT2 level on top of the best-performing combination of electronic structure calculation methods to simulate full IR spectral band-shapes. The chosen molecules in this last step are furan, methyloxirane, and naphthalene, already extensively analyzed in the previous phase, as well as glycine (the four most stable conformers) and fluoro-benzene, depicted in Figure 1. Let us start with furan, shown in Figure 5. The earlier analysis showed a very good agreement with experimental data. The region above 2000 cm⁻¹ has been magnified by a factor of 10 to emphasize the low-intensity vibrational structure due to overtones and combination bands, and thus absent at the harmonic level. In the fingerprint region (below 2000 cm⁻¹), the performance of all methods with GVPT2 are confirmed in terms of band positions. The experimental spectrum shows a broader structure, with a splitting of the bands typical of a rotational sub-structure ignored here. The relative intensities of the peaks are qualitatively reproduced at all levels, with the heights of the peaks at about 750 cm⁻¹ slightly overestimated by theory. The region above 2000 cm⁻¹ exhibits a high density of peaks, dominated by the CH stretching zone above 3000 cm⁻¹, well predicted by calculations in terms of position and shape. The experimental band shows a clear shoulder on the high-energy side, which can be related to the presence of two transitions of similar intensity, $|1_{18}\rangle$ and $|1_{20}\rangle$. In general, the B3PW91 and revDSD-PBEP86-D3(BJ) results are very similar, and so is their combination. An interesting exception is observed for the peak at about 2350 cm⁻¹ related to the first overtone of mode 13, quite intense for B3PW91 and revDSD-PBEP86-D3(BJ) but lower with the hybrid model, matching more closely the experiment.

FIGURE 5

FIGURE 5. Infrared spectrum of furan at the B3PW91/SNSD (B3P/SNSD), revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ (RDS/JNTZ) and revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels using the GVPT2 scheme, compared to experiment taken from (Mellouki et al., 2001). Theoretical line-shapes were broadened by mean of Lorentzian distribution functions with half-widths at half-maximum of 4 cm⁻¹. The region beyond 2000 cm⁻¹ had been magnified to show the lower-intensity structure. To see more clearly the broader features of the experimental spectrum, a higher magnification was used compared to theory.

For methyloxirane as well, shown in Figure 6, the band-shape is well predicted by theory. Due to technical constraints related to the experimental setup, only the fingerprint region is available with a high level of confidence and used here as reference (Merten et al., 2013; Kreienborg et al., 2019). Some larger differences can be observed between the levels of theory, especially in the pattern around 1500 cm⁻¹. B3PW91 with the SNSD basis set shows a series of small peaks of close intensities after the more intense peak at about 1415 cm⁻¹$\left(|1_{15}\rangle \right)$, resulting in a broad pattern. Moreover, the peak experimentally observed at about 1515 cm⁻¹ is not clearly visible and its position seems slightly underestimated. Conversely, revDSD-PBEP86-D3(BJ) with the jun-cc-pVTZ basis set shows a more distinct pattern, with the doublet of peaks present at 1480 cm⁻¹ in the recorded spectrum discernible, and the higher-energy peak of this region is correctly reproduced. However, the intensity of the combination band $|1_8{1}_3\rangle$ is predicted to be strong, erroneously producing a new peak before the doublet. Increasing the basis set used in conjunction with B3PW91 alters remarkably the pattern there, with the doublet now perfectly visible, but the overall number of contributing transitions in the region drops, resulting in a poorer agreement in its higher-energy part. Using this level of theory within the standard hybrid scheme combining revDSD-PBEP86-D3(BJ) and B3PW91 does not alter much the results, thus confirming that the dominant elements in the anharmonic transition intensities are the harmonic components of the electric dipole, that is its first derivatives. While currently not of direct astrochemical interest, the vibrational circular dichroism (VCD) spectrum is also included (lower panel of Figure 6) to show the quality of the anharmonic approach achievable also for techniques with weaker signal, and the straightforward support offered by versatile formulations (Bloino et al., 2015). Because of the lack of implementations of atomic axial tensors at the second-order Møller–Plesset (MP2) level of theory, VCD is currently not available for double hybrid functionals. For this reason, only the electric dipole and harmonic potential energy surface are taken from revDSD-PBEP86-D3(BJ) in the hybrid model, the rest being filled with B3PW91 results. VCD is more sensitive than IR, which also explains the lower accuracy of theory compared to experiment. The sign pattern, that is the alternation of positive and negative bands, is generally well reproduced, with the major discrepancies found around 1500 cm⁻¹. It is noteworthy that the basis set can influence the intensity and even the sign of the bands as seen at 1500 cm⁻¹, where a second negative band appears at the B3PW91/jun-cc-pVTZ level. As a result, larger basis sets may be necessary to reach a convergence of the properties and support newly available spectroscopies in the future, even if smaller ones already provide accurate energies.

FIGURE 6

FIGURE 6. Infrared and vibrational circular dichroism spectra of methyloxirane at the B3PW91/SNSD (B3P/SNSD), B3PW91/jun-cc-pVTZ (B3P/JNTZ), revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ (RDS/JNTZ) and hybrid levels using the GVPT2 scheme, compared to experiment taken from (Merten et al., 2013). Theoretical line-shapes were broadened by mean of Lorentzian distribution functions with half-widths at half-maximum of 3 cm⁻¹. The spectra generated with the anharmonic constants at the B3PW91/jun-cc-pVTZ were inverted (intensity multiplied by −1) to make the difference with B3PW91/SNSD more visible.

Beside its high astronomical interest, another appeal of naphthalene from a theoretical perspective comes from its high symmetry, which gives well-defined bands, with intense overtones and combination bands, especially in the 1500–2000 cm⁻¹ region, as shown in Figure 7. It is thus easier to validate anharmonic models and highlight their improvement over the harmonic approximation, for which scaling factors are insufficient to reach even a qualitative agreement with experiment over the whole mid-IR range. Here, both B3PW91 and the hybrid scheme show very good agreement with the reference spectrum. The pattern of lower-intensity bands is well reproduced, and the band positions match very well their experimental counterparts, with only a minor underestimation of the peaks’ energy at about 1900 cm⁻¹, confirming the low average error found on IR-active fundamental transitions during the second phase of this study. The small bump in the red-wing of the CH-stretching band is consistent with the low-intensity combination bands at about 3000 cm⁻¹. It should be noted that a more thorough analysis of this region, supported by very accurate experimental data has been done in (Mackie et al., 2015). However, because of the density of peaks and the potential couplings between the CH stretching modes, such an extensive analysis would go beyond the purpose of the present study.

FIGURE 7

FIGURE 7. Infrared spectrum of naphthalene at the B3PW91/SNSD (B3P/SNSD), revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels using the GVPT2 scheme, compared to experiment taken from (Hudgins et al., 1994; NIST Mass Spec Data Center and Stein, 2015). The harmonic level was taken at the revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ level, which serves as basis for the hybrid model. Theoretical line-shapes were broadened by mean of Lorentzian distribution functions with half-widths at half-maximum of 4 cm⁻¹.

As an example of halogenated system, the IR spectrum of fluorobenzene has been simulated and is compared to experiment in Figure 8. The main experimental features are well reproduced by all the computational models. However, B3PW91/SNSD provides a general blue-shift of the band positions, which are, instead, closely reproduced at the revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ level. The latter model also predicts very well the band-shape in the CH stretching region. This quality is preserved by the hybrid scheme, which further improves the band position and corrects the band at 1500 cm⁻¹, split at the RDSD level. The low-intensity pattern of overtones and combination bands is also remarkably well described.

FIGURE 8

FIGURE 8. Infrared spectrum of fluorobenzene at the B3PW91/SNSD (B3P/SNSD), revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ (RDS/JNTZ) and revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ//B3PW91/SNSD levels using the GVPT2 scheme, compared to experiment taken from (NIST Mass Spec Data Center and Stein, 2015). Theoretical line-shapes were broadened by mean of Lorentzian distribution functions with half-widths at half-maximum of 4 cm⁻¹.

To conclude this analysis of the quality of vibrational spectra computed at the GVPT2 level with DFT-based electronic structure calculation methods, we consider the IR spectrum of glycine, which has been extensively studied by some of us (Biczysko et al., 2012; Barone et al., 2013a, 2013b, 2013c). The molecule has several stable conformers (Shu et al., 2020) and the first two modes, related to torsions around the C-N and C-C bonds, are large amplitude motions (LAMs), which can be poorly described by quartic force fields in rectilinear coordinates. VPT2 is known to perform badly in these conditions, giving very inaccurate energies for these vibrations, as well as for any mode coupled with them. The standard treatment is to identify the LAMs and exclude them entirely from the anharmonic treatment. However, this was not a problem in the present case for two reasons. First, the two torsional modes are outside the available experimental range, and can be ignored. Second, their coupling with the rest of the system is relatively weak and can be kept, so the full force field was used in the computations. The final spectrum is compared to experiment in Supplementary Figure S13 of the Supplementary Material. In the present study, the four most stable conformers were considered. The main bands are well reproduced, both in terms of position and intensity. revDSD-PBEP86-D3(BJ) and the hybrid scheme reproduce more accurately the band positions compared to B3PW91, the latter being systematically blue-shifted with respect to experiment, especially for the bands beyond 1000 cm⁻¹. Regarding the intensities, B3PW91 underestimate the intensities at about 800 cm⁻¹, while the other methods overestimate them. The largest discrepancies are found in the 1000–1250 cm⁻¹ zone, with revDSD-PBEP86-D3(BJ) and the hybrid model showing three peaks in place of the experimental doublet. Interestingly, the extra band at 1145 cm⁻¹ is very intense for the hybrid scheme compared to revDSD-PBEP86-D3(BJ) or B3PW91. While this feature is related to the combination band $|1_2{1}_{10}\rangle$ of the main conformer (with an abundance of 80% for revDSD-PBEP86-D3(BJ) and 78% at the B3PW91 level), the latter is directly connected through a Fermi resonance to $|1_{12}\rangle$, whose intensity is enhanced in the hybrid scheme compared to revDSD-PBEP86-D3(BJ) or B3PW91 and then redistributed through the variational correction. At variance, B3PW91 predicts only two bands, followed by a smaller one at about 1210 cm⁻¹, also present with revDSD-PBEP86-D3(BJ), but not their combination. The higher energy band-shape is well reproduced, including the bands between 1300 and 1700 cm⁻¹.

Conclusion

In this work, we have introduced different flavors of second-order vibrational perturbation theory, offering thus the possibility to easily choose the most convenient variant for any application in vibrational spectroscopy concerning both energies and intensities. Actually, to the best of our knowledge, this is the first completely general implementation of intensities in the framework of the double-perturbation theory. This general platform has been next employed to select the most effective combination of density functional and basis set for the study of molecular systems of medium-to-large dimensions. The main outcomes of our study can be summarized as follows:

1. Thanks to the availability of analytical second energy derivatives and first property derivatives, a very effective treatment of mechanical and electric/magnetic anharmonicities is possible in the framework of generalized second-order vibrational perturbation theory (GVPT2) for both hybrid and double-hybrid functionals.

2. Contrary to common claims, resonances can be effectively dealt for both energies and intensities in the framework of the GVPT2 model.

3. While the standard thresholds used for the recovery of resonances appear quite robust, some fine tuning of their value can (and should) be performed for specific cases especially for intensities.

4. The new GHDCPT2 variant provides a fully black-box tool for energies.

5. The most effective among the tested functionals are B3PW91 among the hybrid models and rev-DSD-PBEP86-D3(BJ) among the double-hybrid variants.

6. The quality of the results at the double-harmonic level is the most important factor for obtaining accurate frequencies and intensities. Here double-hybrid functionals are significantly more reliable than the hybrid ones, often approaching the accuracy of CCSD(T) computations with comparable basis sets.

7. Anharmonic contributions can be confidently computed by hybrid functionals in most cases.

8. Diffuse functions on non-hydrogen atoms always play a non-negligible role with at least single s, p, d sets being mandatory for intensities at all computational levels.

9. The smallest underlying basis set providing sufficiently converged results is of double-ζ quality for hybrid functionals and of triple-ζ quality for double-hybrid functionals.

10. Empirical dispersion corrections (D3BJ) play a marginal role. However, their inclusion is always suggested because the situation could be different for larger molecules and, in any case, they do not increase the computational cost, nor degrade the results in significant ways.

11. The suggested strategy for reliable yet effective computations of vibrational spectra includes harmonic frequencies obtained at the revDSD-PBEP86-D3(BJ)/jun-cc-pVTZ level together with anharmonic contributions evaluated at the B3PW91/SNSD level.

In summary, DFT/VPT2 represents a very reliable and cost/effective model for the analysis of vibrational spectra of medium- to large-size molecules. Even if further developments are needed especially for flexible systems due to the presence of large amplitude motions, which are poorly described by quartic force fields, the robust implementation of such a platform in a general computer code paves the route toward widespread application also by non-specialists in connection with both assignment and interpretation of experimental vibrational spectra of molecular systems of broad astrochemical interest.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

Q.Y. and M.M. performed the computations and analyzed the results. J.B. supervised the work, analyzed the results and wrote part of the paper. V.B. devised the new basis sets introduced in this work, defined the general strategy and wrote part of the paper.

Funding

This work has been supported by the MIUR ‘PRIN 2017’ (Grant Number 2017A4XRCA) and by the Italian Space Agency (ASI; ‘Life in Space’ project, N. 2019-3-U.0).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The SMART@SNS Laboratory (http://smart.sns.it) is acknowledged for providing high-performance computing facilities.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2021.665232/full#supplementary-material

References

Austin, A., Petersson, G. A., Frisch, M. J., Dobek, F. J., Scalmani, G., and Throssell, K. (2012). A Density Functional with Spherical Atom Dispersion Terms. J. Chem. Theor. Comput. 8, 4989–5007. doi:10.1021/ct300778e

PubMed Abstract | CrossRef Full Text | Google Scholar

Balucani, N. (2012). Elementary Reactions of N Atoms with Hydrocarbons: First Steps towards the Formation of Prebiotic N-Containing Molecules in Planetary Atmospheres. Chem. Soc. Rev. 41, 5473–5483. doi:10.1039/c2cs35113g

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V. (2005). Anharmonic Vibrational Properties by a Fully Automated Second-Order Perturbative Approach. J. Chem. Phys. 122, 014108. doi:10.1063/1.1824881

CrossRef Full Text | Google Scholar

Barone, V., Biczysko, M., Bloino, J., Cimino, P., Penocchio, E., and Puzzarini, C. (2015). 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. Theor. Comput. 11, 4342–4363. doi:10.1021/acs.jctc.5b00580

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Biczysko, M., Bloino, J., and Puzzarini, C. (2013a). Accurate Structure, Thermodynamic and Spectroscopic Parameters from CC and CC/DFT Schemes: the Challenge of the Conformational Equilibrium in glycine. Phys. Chem. Chem. Phys. 15, 10094–10111. doi:10.1039/C3CP50439E

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Biczysko, M., Bloino, J., and Puzzarini, C. (2013b). Characterization of the Elusive Conformers of glycine from State-Of-The-Art Structural, Thermodynamic, and Spectroscopic Computations: Theory Complements Experiment. J. Chem. Theor. Comput. 9, 1533–1547. doi:10.1021/ct3010672

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Biczysko, M., Bloino, J., and Puzzarini, C. (2013c). Glycine Conformers: a Never-Ending Story? Phys. Chem. Chem. Phys. 15, 1358–1363. doi:10.1039/C2CP43884D

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Ceselin, G., Fusè, M., and Tasinato, N. (2020). Accuracy Meets Interpretability for Computational Spectroscopy by Means of Hybrid and Double-Hybrid Functionals. Front. Chem. 8, 859. doi:10.3389/fchem.2020.584203

PubMed Abstract | CrossRef Full Text | Google Scholar

Basire, M., Parneix, P., and Calvo, F. (2008). Quantum Anharmonic Densities of States Using the Wang-Landau Method. J. Chem. Phys. 129, 081101. doi:10.1063/1.2965905

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (1993). Density‐functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 98, 5648–5652. doi:10.1063/1.464913

CrossRef Full Text | Google Scholar

Begue, D., Carbonniere, P., and Pouchan, C. (2005). Calculations of Vibrational Energy Levels by Using a Hybrid Ab Initio and DFT Quartic Force Field: Application to Acetonitrile. J. Phys. Chem. A. 109, 4611–4616. doi:10.1021/jp0406114

PubMed Abstract | CrossRef Full Text | Google Scholar

Biczysko, M., Bloino, J., Carnimeo, I., Panek, P., and Barone, V. (2012). Fully Ab Initio IR Spectra for Complex Molecular Systems from Perturbative Vibrational Approaches: Glycine as a Test Case. J. Mol. Struct. 1009, 74–82. doi:10.1016/j.molstruc.2011.10.012

CrossRef Full Text | Google Scholar

Biczysko, M., Bloino, J., and Puzzarini, C. (2018). Computational Challenges in Astrochemistry. Wires Comput. Mol. Sci. 8, e1349. doi:10.1002/wcms.1349

CrossRef Full Text | Google Scholar

Biczysko, M., Panek, P., Scalmani, G., Bloino, J., and Barone, V. (2010). Harmonic and Anharmonic Vibrational Frequency Calculations with the Double-Hybrid B2PLYP Method: Analytic Second Derivatives and Benchmark Studies. J. Chem. Theor. Comput. 6, 2115–2125. doi:10.1021/ct100212p

PubMed Abstract | CrossRef Full Text | Google Scholar

Biczysko, M., Tarroni, R., and Carter, S. (2003). Variational Calculations of HBN Energy Levels in the X²Π and A²Σ⁺ States. J. Chem. Phys. 119, 4197–4203. doi:10.1063/1.1594174

CrossRef Full Text | Google Scholar

Bloino, J. (2015). A VPT2 Route to Near-Infrared Spectroscopy: The Role of Mechanical and Electrical Anharmonicity. J. Phys. Chem. A. 119, 5269–5287. doi:10.1021/jp509985u

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloino, J., Baiardi, A., and Biczysko, M. (2016). Aiming at an Accurate Prediction of Vibrational and Electronic Spectra for Medium‐to‐large Molecules: An Overview. Int. J. Quan. Chem. 116, 1543–1574. doi:10.1002/qua.25188

CrossRef Full Text | Google Scholar

Bloino, J., and Barone, V. (2012). A Second-Order Perturbation Theory Route to Vibrational Averages and Transition Properties of Molecules: General Formulation and Application to Infrared and Vibrational Circular Dichroism Spectroscopies. J. Chem. Phys. 136, 124108. doi:10.1063/1.3695210

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloino, J., Biczysko, M., and Barone, V. (2015). Anharmonic Effects on Vibrational Spectra Intensities: Infrared, Raman, Vibrational Circular Dichroism, and Raman Optical Activity. J. Phys. Chem. A. 119, 11862–11874. doi:10.1021/acs.jpca.5b10067

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloino, J., Biczysko, M., and Barone, V. (2012). General Perturbative Approach for Spectroscopy, Thermodynamics, and Kinetics: Methodological Background and Benchmark Studies. J. Chem. Theor. Comput. 8, 1015–1036. doi:10.1021/ct200814m

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulik, I. W., Frisch, M. J., and Vaccaro, P. H. (2017). Vibrational Self-Consistent Field Theory Using Optimized Curvilinear Coordinates. J. Chem. Phys. 147, 044110. doi:10.1063/1.4995440

PubMed Abstract | CrossRef Full Text | Google Scholar

Bunker, P. R., and Jensen, P. (1998). Molecular Symmetry and Spectroscopy. Ottawa, Canada: NRC Research Press.

Cané, E., Miani, A., and Trombetti, A. (2007). Anharmonic Force Fields of Naphthalene-H₈ and Naphthalene-D₈. J. Phys. Chem. A. 111, 8218–8222. doi:10.1021/jp071610p

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbonniere, P., Lucca, T., Pouchan, C., Rega, N., and Barone, V. (2005). Vibrational Computations beyond the Harmonic Approximation: Performances of the B3LYP Density Functional for Semirigid Molecules. J. Comput. Chem. 26, 384–388. doi:10.1002/jcc.20170

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter, S., Handy, N. C., Puzzarini, C., Tarroni, R., and Palmieri, P. (2000). A Variational Method for the Calculation of Spin-Rovibronic Energy Levels of Triatomic Molecules with Three Interacting Electronic States. Mol. Phys. 98, 1697–1712. doi:10.1080/00268970009483375

CrossRef Full Text | Google Scholar

Carter, S., Handy, N. C., Rosmus, P., and Chambaud, G. (1990). A Variational Method for the Calculation of Spin-Rovibronic Levels of Renner-Teller Triatomic Molecules. Mol. Phys. 71, 605–622. doi:10.1080/00268979000102001

CrossRef Full Text | Google Scholar

Chyba, C., and Sagan, C. (1992). Endogenous Production, Exogenous Delivery and Impact-Shock Synthesis of Organic Molecules: An Inventory for the Origins of Life. Nature 355, 125–132. doi:10.1038/355125a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Császár, A. G. (2011). Anharmonic Molecular Force Fields. Wires Comput. Mol. Sci. 2, 273–289. doi:10.1002/wcms.75

CrossRef Full Text | Google Scholar

Császár, A. G., Fábri, C., Szidarovszky, T., Mátyus, E., Furtenbacher, T., and Czakó, G. (2012). The Fourth Age of Quantum Chemistry: Molecules in Motion. Phys. Chem. Chem. Phys. 14, 1085–1106. doi:10.1039/C1CP21830A

PubMed Abstract | CrossRef Full Text | Google Scholar

Darling, B. T., and Dennison, D. M. (1940). The Water Vapor Molecule. Phys. Rev. 57, 128–139. doi:10.1103/PhysRev.57.128

CrossRef Full Text | Google Scholar

Des Marais, D. J., Harwit, M. O., Jucks, K. W., Kasting, J. F., Lin, D. N. C., Lunine, J. I., et al. (2002). Remote Sensing of Planetary Properties and Biosignatures on Extrasolar Terrestrial Planets. Astrobiology 2, 153–181. doi:10.1089/15311070260192246

PubMed Abstract | CrossRef Full Text | Google Scholar

Dierksen, M., and Grimme, S. (2004). The Vibronic Structure of Electronic Absorption Spectra of Large Molecules: A Time-dependent Density Functional Study on the Influence of "Exact" Hartree−Fock Exchange. J. Phys. Chem. A. 108, 10225–10237. doi:10.1021/jp047289h

CrossRef Full Text | Google Scholar

Dunning, T. H. (1989). Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 90, 1007–1023. doi:10.1063/1.456153

CrossRef Full Text | Google Scholar

Eckart, C. (1935). Some Studies Concerning Rotating Axes and Polyatomic Molecules. Phys. Rev. 47, 552–558. doi:10.1103/PhysRev.47.552

CrossRef Full Text | Google Scholar

Erfort, S., Tschöpe, M., and Rauhut, G. (2020). Toward a Fully Automated Calculation of Rovibrational Infrared Intensities for Semi-rigid Polyatomic Molecules. J. Chem. Phys. 152, 244104–244114. doi:10.1063/5.0011832

PubMed Abstract | CrossRef Full Text | Google Scholar

Franke, P. R., Stanton, J. F., and Douberly, G. E. (2021). 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. doi:10.1021/acs.jpca.0c09526

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2016). Gaussian 16 Revision A.03. Wallingford CT: Gaussian Inc.

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2020). Gaussian Development Version, Revision J.13. Wallingford CT: Gaussian, Inc.

Fusè, M., Mazzeo, G., Longhi, G., Abbate, S., Masi, M., Evidente, A., et al. (2019). Unbiased Determination of Absolute Configurations by Vis-À-Vis Comparison of Experimental and Simulated Spectra: The Challenging Case of Diplopyrone. J. Phys. Chem. B. 123, 9230–9237. doi:10.1021/acs.jpcb.9b08375

PubMed Abstract | CrossRef Full Text | Google Scholar

Girish, T. E., Sony, K. S., Vaidyan, V. K., and Jayakumar, V. S. (2008). Vibrational Spectroscopy and Search for Extraterrestrial Life. AIP Conf. Proc. 1075, 196–199. doi:10.1063/1.3046211

CrossRef Full Text | Google Scholar

Goerigk, L., Hansen, A., Bauer, C., Ehrlich, S., Najibi, A., and Grimme, S. (2017). A Look at the Density Functional Theory Zoo with the Advanced GMTKN55 Database for General Main Group Thermochemistry, Kinetics and Noncovalent Interactions. Phys. Chem. Chem. Phys. 19, 32184–32215. doi:10.1039/c7cp04913g

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Antony, J., Ehrlich, S., and Krieg, H. (2010). A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 132, 154104. doi:10.1063/1.3382344

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., Ehrlich, S., and Goerigk, L. (2011). Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 32, 1456–1465. doi:10.1002/jcc.21759

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S. (2006). Semiempirical Hybrid Density Functional with Perturbative Second-Order Correlation. J. Chem. Phys. 124, 034108. doi:10.1063/1.2148954

PubMed Abstract | CrossRef Full Text | Google Scholar

Herbst, E., and van Dishoeck, E. F. (2009). Complex Organic Interstellar Molecules. Annu. Rev. Astron. Astrophys. 47, 427–480. doi:10.1146/annurev-astro-082708-101654

CrossRef Full Text | Google Scholar

Hörst, S. M., Yelle, R. V., Buch, A., Carrasco, N., Cernogora, G., Dutuit, O., et al. (2012). Formation of Amino Acids and Nucleotide Bases in a Titan Atmosphere Simulation Experiment. Astrobiology 12, 809–817. doi:10.1089/ast.2011.0623

PubMed Abstract | CrossRef Full Text | Google Scholar

Hudgins, D. M., Sandford, S. A., and Allamandola, L. J. (1994). Infrared Spectroscopy of Polycyclic Aromatic Hydrocarbon Cations. 1. Matrix-Isolated Naphthalene and Perdeuterated Naphthalene. J. Phys. Chem. 98, 4243–4253. doi:10.1021/jp509985u

PubMed Abstract | CrossRef Full Text | Google Scholar

Isoniemi, E., Pettersson, M., Khriachtchev, L., Lundell, J., and Räsänen, M. (1999). Infrared Spectroscopy of H₂S and SH in Rare-Gas Matrixes. J. Phys. Chem. A. 103, 679–685. doi:10.1021/jp9838893

CrossRef Full Text | Google Scholar

Jensen, F. (2001). Polarization Consistent Basis Sets: Principles. J. Chem. Phys. 115, 9113–9125. doi:10.1063/1.1413524

CrossRef Full Text | Google Scholar

Kendall, R. A., Dunning, T. H., and Harrison, R. J. (1992). Electron Affinities of the First‐row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 96, 6796–6806. doi:10.1063/1.462569

CrossRef Full Text | Google Scholar

Krasnoshchekov, S. V., Craig, N. C., and Stepanov, N. F. (2013). Anharmonic Vibrational Analysis of the Gas-Phase Infrared Spectrum of 1,1-Difluoroethylene Using the Operator Van Vleck Canonical Perturbation Theory. J. Phys. Chem. A. 117, 3041–3056. doi:10.1021/jp311398z

PubMed Abstract | CrossRef Full Text | Google Scholar

Krasnoshchekov, S. V., Dobrolyubov, E. O., Syzgantseva, M. A., and Palvelev, R. V. (2020). Rigorous Vibrational Fermi Resonance Criterion Revealed: Two Different Approaches Yield the Same Result. Mol. Phys. 118, e1743887. doi:10.1080/00268976.2020.1743887

CrossRef Full Text | Google Scholar

Krasnoshchekov, S. V., Isayeva, E. V., and Stepanov, N. F. (2014). 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. doi:10.1063/1.4903927

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreienborg, N. M., Bloino, J., Osowski, T., Pollok, C. H., and Merten, C. (2019). The Vibrational CD Spectra of Propylene Oxide in Liquid Xenon: a Proof-Of-Principle CryoVCD Study that Challenges Theory. Phys. Chem. Chem. Phys. 21, 6582–6587. doi:10.1039/C9CP00537D

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruse, H., Szabla, R., and Šponer, J. (2020). Surprisingly Broad Applicability of the Cc-pVnZ-F12 Basis Set for Ground and Excited States. J. Chem. Phys. 152, 214104. doi:10.1063/5.0006871

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhler, K. M., Truhlar, D. G., and Isaacson, A. D. (1996). General Method for Removing Resonance Singularities in Quantum Mechanical Perturbation Theory. J. Chem. Phys. 104, 4664–4671. doi:10.1063/1.471161

CrossRef Full Text | Google Scholar

Mackie, C. J., Candian, A., Huang, X., Maltseva, E., Petrignani, A., Oomens, J., et al. (2015). The Anharmonic Quartic Force Field Infrared Spectra of Three Polycyclic Aromatic Hydrocarbons: Naphthalene, Anthracene, and Tetracene. J. Chem. Phys. 143, 224314. doi:10.1063/1.4936779

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackie, C. J., Candian, A., Huang, X., Maltseva, E., Petrignani, A., Oomens, J., et al. (2016). The Anharmonic Quartic Force Field Infrared Spectra of Five Non-linear Polycyclic Aromatic Hydrocarbons: Benz[a]anthracene, Chrysene, Phenanthrene, Pyrene, and Triphenylene. J. Chem. Phys. 145, 084313. doi:10.1063/1.4961438

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardirossian, N., and Head-Gordon, M. (2017). Thirty Years of Density Functional Theory in Computational Chemistry: an Overview and Extensive Assessment of 200 Density Functionals. Mol. Phys. 115, 2315–2372. doi:10.1080/00268976.2017.1333644

CrossRef Full Text | Google Scholar

Martin, J. M. L., Lee, T. J., Taylor, P. R., and François, J. P. (1995). The Anharmonic Force Field of Ethylene, C₂H₄, by Means of Accurate Ab Initio Calculations. J. Chem. Phys. 103, 2589–2602. doi:10.1063/1.469681

CrossRef Full Text | Google Scholar

Martin, J. M. L., and Taylor, P. R. (1997). Accurate Ab Initio Quartic Force Field for trans-HNNH and Treatment of Resonance Polyads. Spectrochimica Acta A: Mol. Biomol. Spectrosc. 53, 1039–1050. doi:10.1016/S1386-1425(96)01869-0

CrossRef Full Text | Google Scholar

Mátyus, E., Czakó, G., and Császár, A. G. (2009). Toward Black-box-type Full- and Reduced-Dimensional Variational (Ro)vibrational Computations. J. Chem. Phys. 130, 134112. doi:10.1063/1.3076742

PubMed Abstract | CrossRef Full Text | Google Scholar

McGuire, B. A. (2018). 2018 Census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disk, and Exoplanetary Molecules. Astrophysical J. Suppl. Ser. 239, 17. doi:10.3847/1538-4365/aae5d2

CrossRef Full Text | Google Scholar

McGuire, B. A., Loomis, R. A., Burkhardt, A. M., Lee, K. L. K., Shingledecker, C. N., Charnley, S. B., et al. (2021). Detection of Two Interstellar Polycyclic Aromatic Hydrocarbons via Spectral Matched Filtering. Science 371, 1265–1269. doi:10.1126/science.abb7535

PubMed Abstract | CrossRef Full Text | Google Scholar

Mellouki, A., Liévin, J., and Herman, M. (2001). The Vibrational Spectrum of Pyrrole (C₄H₅N) and Furan (C₄H₄O) in the Gas Phase. Chem. Phys. 271, 239–266. doi:10.1016/S0301-0104(01)00447-5

CrossRef Full Text | Google Scholar

Merten, C., Bloino, J., Barone, V., and Xu, Y. (2013). Anharmonicity Effects in the Vibrational CD Spectra of Propylene Oxide. J. Phys. Chem. Lett. 4, 3424–3428. doi:10.1021/jz401854y

CrossRef Full Text | Google Scholar

Mitra, H., and Roy, T. K. (2020). Comprehensive Benchmark Results for the Accuracy of Basis Sets for Anharmonic Molecular Vibrations. J. Phys. Chem. A. 124, 9203–9221. doi:10.1021/acs.jpca.0c06634

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitrushchenkov, A. O. (2012). A New General Renner-Teller (Including Ɛ ≳ 1) Spectroscopic Formalism for Triatomic Molecules. J. Chem. Phys. 136, 024108. doi:10.1063/1.3672162

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakanaga, T., Kondo, S., and Saëki, S. (1979). Infrared Intensities and Coriolis Interactions in Ethylene. J. Chem. Phys. 70, 2471–2478. doi:10.1063/1.437709

CrossRef Full Text | Google Scholar

Nauts, A., and Lauvergnat, D. (2012). Quantum Dynamics of Floppy Molecular Systems with Elvibrot and Tnum. AIP Conf. Proc. 1504, 948–952. doi:10.1063/1.4771853

CrossRef Full Text | Google Scholar

Nejad, A., Meyer, E., and Suhm, M. A. (2020). Glycolic Acid as a Vibrational Anharmonicity Benchmark. J. Phys. Chem. Lett. 11, 5228–5233. doi:10.1021/acs.jpclett.0c01462

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, H. H. (1951). The Vibration-Rotation Energies of Molecules. Rev. Mod. Phys. 23, 90–136. doi:10.1103/RevModPhys.23.90

CrossRef Full Text | Google Scholar

NIST Mass Spec Data Center, d., and Stein, S. E. (2015). NIST Chemistry WebBook, NIST Standard Reference Database Number 69 (Gaithersburg MD, 20899: National Institute of Standards and Technology), Available at: http://webbook.nist.gov (Accessed April 12, 2021).

Google Scholar

Papajak, E., Zheng, J., Xu, X., Leverentz, H. R., and Truhlar, D. G. (2011). Perspectives on Basis Sets Beautiful: Seasonal Plantings of Diffuse Basis Functions. J. Chem. Theor. Comput. 7, 3027–3034. doi:10.1021/ct200106a

PubMed Abstract | CrossRef Full Text | Google Scholar

Papousek, D., and Aliev, M. R. (1982). Molecular Vibrational-Rotational Spectra. New York): Elsevier.

Papp, D., Szidarovszky, T., and Császár, A. G. (2017). 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. doi:10.1063/1.5000680

PubMed Abstract | CrossRef Full Text | Google Scholar

Piccardo, M., Bloino, J., and Barone, V. (2015). Generalized Vibrational Perturbation Theory for Rotovibrational Energies of Linear, Symmetric and Asymmetric Tops: Theory, Approximations, and Automated Approaches to Deal with Medium‐to‐large Molecular Systems. Int. J. Quan. Chem. 115, 948–982. doi:10.1002/qua.24931

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietropolli Charmet, A., Stoppa, P., Tasinato, N., Giorgianni, S., and Gambi, A. (2016). Study of the Vibrational Spectra and Absorption Cross Sections of 1-Chloro-1-Fluoroethene by a Joint Experimental and Ab Initio Approach. J. Phys. Chem. A. 120, 8369–8386. doi:10.1021/acs.jpca.6b07426

PubMed Abstract | CrossRef Full Text | Google Scholar

Pirali, O., Goubet, M., Huet, T. R., Georges, R., Soulard, P., Asselin, P., et al. (2013). The Far Infrared Spectrum of Naphthalene Characterized by High Resolution Synchrotron FTIR Spectroscopy and Anharmonic DFT Calculations. Phys. Chem. Chem. Phys. 15, 10141–10150. doi:10.1039/C3CP44305A

PubMed Abstract | CrossRef Full Text | Google Scholar

Pirali, O., Vervloet, M., Mulas, G., Malloci, G., and Joblin, C. (2009). High-resolution Infrared Absorption Spectroscopy of Thermally Excited Naphthalene. Measurements and Calculations of Anharmonic Parameters and Vibrational Interactions. Phys. Chem. Chem. Phys. 11, 3443–3454. doi:10.1039/B814037E

PubMed Abstract | CrossRef Full Text | Google Scholar

Plíva, J. (1990). Anharmonic Constants for Degenerate Modes of Symmetric Top Molecules. J. Mol. Spectrosc. 139, 278–285. doi:10.1016/0022-2852(90)90065-X

CrossRef Full Text | Google Scholar

Polavarapu, P. L., Hess, B. A., and Schaad, L. J. (1985). Vibrational Spectra of Epoxypropane. J. Chem. Phys. 82, 1705–1710. doi:10.1063/1.448403

CrossRef Full Text | Google Scholar

Puzzarini, C., and Barone, V. (2020). The Challenging Playground of Astrochemistry: an Integrated Rotational Spectroscopy—Quantum Chemistry Strategy. Phys. Chem. Chem. Phys. 22, 6507–6523. doi:10.1039/d0cp00561d

PubMed Abstract | CrossRef Full Text | Google Scholar

Puzzarini, C., and Barone, V. (2008). Toward Spectroscopic Accuracy for Organic Free Radicals: Molecular Structure, Vibrational Spectrum, and Magnetic Properties of F₂NO. J. Chem. Phys. 129, 084306. doi:10.1063/1.2969820

PubMed Abstract | CrossRef Full Text | Google Scholar

Puzzarini, C., Biczysko, M., and Barone, V. (2010). Accurate Harmonic/anharmonic Vibrational Frequencies for Open-Shell Systems: Performances of the B3LYP/N07D Model for Semirigid Free Radicals Benchmarked by CCSD(T) Computations. J. Chem. Theor. Comput. 6, 828–838. doi:10.1021/ct900594h

PubMed Abstract | CrossRef Full Text | Google Scholar

Puzzarini, C., Bloino, J., Tasinato, N., and Barone, V. (2019). Accuracy and Interpretability: The Devil and the Holy Grail. New Routes across Old Boundaries in Computational Spectroscopy. Chem. Rev. 119, 8131–8191. doi:10.1021/acs.chemrev.9b00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Raghavachari, K., Trucks, G. W., Pople, J. A., and Head-Gordon, M. (1989). A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 157, 479–483. doi:10.1016/S0009-2614(89)87395-6

CrossRef Full Text | Google Scholar

Reisner, D. E., Field, R. W., Kinsey, J. L., and Dai, H. L. (1984). Stimulated Emission Spectroscopy: A Complete Set of Vibrational Constants for X̃¹A₁ Formaldehyde. J. Chem. Phys. 80, 5968–5978. doi:10.1063/1.446677

CrossRef Full Text | Google Scholar

Rosnik, A. M., and Polik, W. F. (2014). 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. doi:10.1080/00268976.2013.808386

CrossRef Full Text | Google Scholar

Saladino, R., Botta, G., Pino, S., Costanzo, G., and Di Mauro, E. (2012). Genetics First or Metabolism First? The Formamide Clue. Chem. Soc. Rev. 41, 5526–5565. doi:10.1039/c2cs35066a

PubMed Abstract | CrossRef Full Text | Google Scholar

Saladino, R., Carota, E., Botta, G., Kapralov, M., Timoshenko, G. N., et al. (2015). Meteorite-catalyzed Syntheses of Nucleosides and of Other Prebiotic Compounds from Formamide under Proton Irradiation. Proc. Natl. Acad. Sci. USA 112, E2746–E2755. doi:10.1073/pnas.1422225112

PubMed Abstract | CrossRef Full Text | Google Scholar

Santra, G., Sylvetsky, N., and Martin, J. M. L. (2019). 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. doi:10.1021/acs.jpca.9b03157

PubMed Abstract | CrossRef Full Text | Google Scholar

Seager, S. (2014). The Future of Spectroscopic Life Detection on Exoplanets. Proc. Natl. Acad. Sci. 111, 12634–12640. doi:10.1073/pnas.1304213111

PubMed Abstract | CrossRef Full Text | Google Scholar

Shu, C., Jiang, Z., and Biczysko, M. (2020). Toward Accurate Prediction of Amino Acid Derivatives Structure and Energetics from DFT: glycine Conformers and Their Interconversions. J. Mol. Model. 26, 129. doi:10.1007/s00894-020-4342-7

PubMed Abstract | CrossRef Full Text | Google Scholar

SNSDT (2021). Double and Triple-ζ Basis Sets of SNS Family. Available at: https://smart.sns.it/?pag=download (Accessed April 12, 2021).

Google Scholar

Sunahori, F. X., Su, Z., Kang, C., and Xu, Y. (2010). Infrared Diode Laser Spectroscopic Investigation of Four C-H Stretching Vibrational Modes of Propylene Oxide. Chem. Phys. Lett. 494, 14–20. doi:10.1016/j.cplett.2010.05.072

CrossRef Full Text | Google Scholar

Suzuki, E., Yamazaki, M., and Shimizu, K. (2007). Infrared Spectra of Monomeric Thioformaldehyde in Ar, N₂ and Xe Matrices. Vibrational Spectrosc. 43, 269–273. doi:10.1016/j.vibspec.2006.02.007

CrossRef Full Text | Google Scholar

Swofford, R. L., Long, M. E., and Albrecht, A. C. (1976). C-H Vibrational States of Benzene, Naphthalene, and Anthracene in the Visible Region by Thermal Lensing Spectroscopy and the Local Mode Model. J. Chem. Phys. 65, 179–190. doi:10.1063/1.432815

CrossRef Full Text | Google Scholar

Tasinato, N., Stoppa, P., Pietropolli Charmet, A., Giorgianni, S., and Gambi, A. (2006). High-Resolution Infrared Study of Vinyl Fluoride in the 750−1050 cm⁻¹ Region: Rovibrational Analysis and Resonances Involving the ν₈, ν₁₀, and ν₁₁ Fundamentals. J. Phys. Chem. A. 110, 13412–13418. doi:10.1021/jp064013w

PubMed Abstract | CrossRef Full Text | Google Scholar

Tennyson, J., Zobov, N. F., Williamson, R., Polyansky, O. L., and Bernath, P. F. (2001). Experimental Energy Levels of the Water Molecule. J. Phys. Chem. Reference Data 30, 735–831. doi:10.1063/1.1364517

CrossRef Full Text | Google Scholar

Van Vleck, J. H. (1929). On σ-Type Doubling and Electron Spin in the Spectra of Diatomic Molecules. Phys. Rev. 33, 467–506. doi:10.1103/PhysRev.33.467

CrossRef Full Text | Google Scholar

Vázquez, J., and Stanton, J. F. (2006). Simple(r) Algebraic Equation for Transition Moments of Fundamental Transitions in Vibrational Second-Order Perturbation Theory. Mol. Phys. 104, 377–388. doi:10.1080/00268970500290367

CrossRef Full Text | Google Scholar

Watson, J. K. G. (1968). Simplification of the Molecular Vibration-Rotation Hamiltonian. Mol. Phys. 15, 479–490. doi:10.1080/00268976800101381

CrossRef Full Text | Google Scholar

Weigend, F., and Ahlrichs, R. (2005). Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 7, 3297–3305. doi:10.1039/B508541A

PubMed Abstract | CrossRef Full Text | Google Scholar

Willets, A., Handy, N. C., Green, W. H., and Jayatilaka, D. (1990). Anharmonic Corrections to Vibrational Transition Intensities. J. Phys. Chem. 94, 5608–5616. doi:10.1021/j100377a038

CrossRef Full Text | Google Scholar

Yurchenko, S. N., Lodi, L., Tennyson, J., and Stolyarov, A. V. (2016). Duo: A General Program for Calculating Spectra of Diatomic Molecules. Comput. Phys. Commun. 202, 262–275. doi:10.1016/j.cpc.2015.12.021

CrossRef Full Text | Google Scholar