# Accuracy Meets Interpretability for Computational Spectroscopy by Means of Hybrid and Double-Hybrid Functionals

- SMART Laboratory, Scuola Normale Superiore di Pisa, Pisa, Italy

Accuracy and interpretability are often seen as the devil and holy grail in computational spectroscopy and their reconciliation remains a primary research goal. In the last few decades, density functional theory has revolutionized the situation, paving the way to reliable yet effective models for medium size molecules, which could also be profitably used by non-specialists. In this contribution we will compare the results of some widely used hybrid and double hybrid functionals with the aim of defining the most suitable recipe for all the spectroscopic parameters of interest in rotational and vibrational spectroscopy, going beyond the rigid rotor/harmonic oscillator model. We will show that last-generation hybrid and double hybrid functionals in conjunction with partially augmented double- and triple-zeta basis sets can offer, in the framework of second order vibrational perturbation theory, a general, robust, and user-friendly tool with unprecedented accuracy for medium-size semi-rigid molecules.

## 1. Introduction

Spectroscopic techniques are unique tools to non-invasively probe the properties of complex molecular systems in a variety of environments and conditions. In fact, the increasing sophistication of well-established techniques like nuclear magnetic and electron paramagnetic resonance (NMR and EPR), microwave (MW), infrared (IR), Raman, visible (Vis), ultra-violet (UV), or fluorescence and the parallel blooming of new ones, e.g., vibrational, electronic, and magnetic circular dichroism (VCD, ECD, and MCD), Raman optical activity (ROA), circularly polarized luminescence (CPL), multi-photon and time-resolved methods have a huge impact in several fields of science and technology (He et al., 2007; Barone, 2012; Berova et al., 2012; Tasinato et al., 2015; Sugiki et al., 2017; Lane, 2018). In addition to being widely used to infer information about molecular structure and dynamics in both gas and condensed phases (Sugiki et al., 2017; Puzzarini and Barone, 2018), spectroscopy allows for the unequivocal identification of chemical species in hostile environments, e.g., the interstellar space (Baiano et al., 2020), or in samples of unknown composition (He et al., 2007; Lindon et al., 2017; Lane, 2018) and plays a pivotal role in the study of photochemical mechanisms in biological systems and in the development of new technological devices, including photovoltaic cells, optoelectronic devices, eco-sustainable solutions, UV-resistant materials, dyes, and fluorescent probes (Berova et al., 2012; Drummen, 2012; Lindon et al., 2017; Lane, 2018). Unfortunately, interpretation of experimental data is often a difficult task: the observed spectroscopic behavior results from the subtle interplay of stereo-electronic, dynamic and environmental effects, whose specific roles are difficult to disentangle. Furthermore, although conveying additional information, spectral congestion makes quantitative interpretation of experimental data even more difficult. To disentangle these complex signatures and disclose the underlying molecular properties, detailed molecular simulations are crucial (Barone, 2012). The last decade has witnessed an increasing interaction between experiment, theory, and simulation in the field of molecular spectroscopy and in all related applications. These have revealed the need of computational tools and theoretical methods not only to interpret the spectra, but also to design new experiments that would be impossible or very expensive to perform in a blind way. The widespread use of computational techniques in many areas of science and by an ever increasing number of non-experienced users (e.g., experimental chemists) has prompted the development, in our group, of the Virtual Multifrequency Spectrometer (VMS, http://dreamslab.sns.it/vms/) (Barone, 2016). VMS features two interconnected tools: VMS-Draw (Licari et al., 2015) and VMS-Comp (Barone et al., 2012), the former providing a user-friendly, graphical user interface to pilot the latter, which is in charge of computationally intensive tasks. VMS-Comp includes a wide set of algorithms and calculation options and it allows the user to predict with remarkable accuracy many types of spectroscopic data for a vast range of molecular systems and environments (Zerbetto et al., 2013; Licari et al., 2017; Presti et al., 2017). Further developments are, however, needed to deal with new and more sophisticated experimental techniques (Quack and Merkt, 2011; Lane, 2018). For small semi-rigid molecules, the accuracy of state-of-the-art quantum mechanical (QM) methods often rivals that of experimental techniques, but extension to large flexible systems (not to mention condensed phases) faces a number of difficulties (even within the Born-Oppenheimer approximation) ranging from the very unfavorable scaling of those methods with the number of active electrons to the proper description of flat potential energy surfaces (PESs) with stationary points that are ill defined (Puzzarini et al., 2019a).

A possible route to obtain accurate results, even for relatively large molecular systems (a few dozens of atoms), is provided by hybrid QM/QM′ models, which combine accurate quantum-mechanical (QM) calculations of “primary” properties (e.g., molecular structures or harmonic force fields) with cheaper yet reliable electronic structure approaches (QM′) for “secondary” properties (e.g., vibrational corrections or anharmonic effects). At the same time, computation of spectroscopic parameters often requires purposely tailored basis sets, whose selection must be based on extensive benchmarks. Finally, the customary rigid rotor (RR) / harmonic oscillator (HO) approximation is not sufficient for quantitative work and more advanced models must be employed. In the last few years, the second-order vibrational perturbation theory (VPT2) (Mills, 1972; Papoušek and Aliev, 1982; Aliev and Watson, 1985) has been exploited with considerable success for semi-rigid molecules of increasing dimensions. However, the identification of resonances in VPT2 treatments remains a daunting task due to the arbitrariness of their definition and their indirect influence on the energy and intensities. This leads to two distinct issues: finding the true resonances and then correcting them appropriately. The resonance conditions are strongly related to the quality of the electronic structure calculation, but they also depend on the coupling between the potentially resonant states. The most obvious solution is to combine perturbative and variational models, with the generalized (G)VPT2 model offering a very good accuracy-to-cost ratio not only for energies, but also for transition moments (Bloino et al., 2012). In recent years, we have implemented analytical second derivatives for double-hybrid functionals (Biczysko et al., 2010), general equations for Abelian and non-Abelian symmetry groups (Piccardo et al., 2015a), and full treatment of intensities for all conventional (IR, Raman) and chiral (VCD, ROA) vibrational spectroscopies up to three-quanta excitations (Bloino et al., 2015). These features compose, together with the particular care devoted to robustness, ease of use, and computational efficiency, the mandatory background for systematic evaluations of all the spectroscopic parameters beyond the rigid-rotor/harmonic-oscillator level for medium- to large-size semi-rigid molecules. Next, together with further developments for the treatment of flexible molecules (Puzzarini et al., 2019b), the selection of the most effective electronic structure models remains a central issue of computational spectroscopy. Here, the ongoing developments of methods rooted in the density functional theory (DFT) come into play. For microwave and vibrational spectroscopic applications, global hybrid density functionals (DFs), such as B3LYP (Lee et al., 1988; Becke, 1993) coupled to a polarized double-zeta basis set supplemented by a set of *sp* diffuse functions (hereafter B3) can deliver the accuracy required for the interpretation of the vibrational characteristics of medium and large molecules beyond the harmonic approximation for both transition frequencies and intensities. An increased accuracy, at the price of a more demanding computational loading, is brought by the double-hybrid B2PLYP functional (Grimme, 2006) in conjunction with partially augmented triple-ζ basis sets (hereafter B2). In fact, it has been shown that B2 calculations of vibrational frequencies and intensities can reach an average accuracy often within 10 cm^{−1} and a few km mol^{−1}, respectively (Biczysko et al., 2010; Puzzarini et al., 2019a; Boussessi et al., 2020a), thus performing equally to, or even better than, the CCSD(T)/cc-pVTZ approach. Concerning the prediction of rotational spectroscopic parameters, the same B2/B3 approach has been validated in several studies concerning rotational constants (Spada et al., 2017a; Li et al., 2018) together with quartic centrifugal distortion parameters (Tasinato, 2014; Boussessi et al., 2020a) and, very recently, also for sextic centrifugal distortion parameters (Pietropolli Charmet et al., 2017; Boussessi et al., 2020b). In the present work, hybrid and double-hybrid density functionals of the “last-generation” are analyzed to check if they provide improved performances with respect to the B3 and B2 paradigms. This is particularly significant since the improved performances for thermochemistry, kinetics and non-covalent interactions could have been obtained at the price of worsening other parameters of particular relevance for rotational and vibrational spectroscopy of closed- and open-shell species (Puzzarini et al., 2010). In parallel, the accuracy of different partially augmented correlation consistent basis sets is analyzed with the aim of defining the best choice in terms of the cost-to-accuracy ratio. We will consider, in particular, equilibrium geometries, ground state rotational constants and quartic centrifugal distortion constants, harmonic and anharmonic vibrational wavenumbers and IR intensities. The benchmark study is carried out on a set of ten molecules of atmospheric and astrochemical relevance, reported in Figure 1, which includes: difluoromethane (CH_{2}F_{2}) (Carlotti et al., 1988; Tasinato et al., 2012b; Piccardo et al., 2015b), chlorofluoromethane (CH_{2}ClF) (Blanco et al., 1995; Pietropolli Charmet et al., 2013), cis-1-chloro-2-fluoroethene (cis-ClHC=CHF) (Craig et al., 1970; Alonso et al., 1993; Gambi et al., 2002; Piccardo et al., 2015b), 1-chloro-1-fluoroethene (ClFC=CH_{2}) (Leung et al., 2009; Pietropolli Charmet et al., 2016; Gambi et al., 2019), chlorotrifluoroethene (F_{2}C=CFCl) (Hillig et al., 1988; Tasinato et al., 2012a), oxirane (cyc-C_{2}H_{4}O) (Russell and Wesendrup, 2003; Flaud et al., 2012; Medcraft et al., 2012; Lafferty et al., 2013; Puzzarini et al., 2014a; Piccardo et al., 2015b), glycolaldehyde (HOCH_{2}CHO) (Carroll et al., 2010; Johnson et al., 2013; Piccardo et al., 2015b; Boussessi et al., 2020a), E-ethanimine (CH_{3}CHNH) (Melli et al., 2018), sulfur dioxide (SO_{2}) (Flaud et al., 1993; Mller and Brnken, 2005; Tasinato et al., 2010; Boussessi et al., 2020a), and the gauche conformer of ethyl mercaptan (CH_{3}CH_{2}SH) (Smith et al., 1968; Wolff and Szydlowski, 1985; Miller et al., 2009; Kolesnikov et al., 2014; Puzzarini et al., 2014b; Hochlaf et al., 2015).

## 2. Computational Methodology

Quantum chemical calculations at the DFT level were carried out using some hybrid and double-hybrid density functional approximations of the last generation, which are considered the best performing and transferable according to a very recent benchmark (Peveratti, 2020). Among hybrid functionals, the PW6B95 meta exchange-correlation functional proposed by Zhao and Truhlar (2005) and the ωB97 family of long-range corrected functionals introduced by Chai and Head-Gordon (Chai and Head-Gordon, 2008a,b), namely ωB97, ωB97X and ωB97X-D, were considered in conjunction with aug-cc-pVDZ (Dunning, 1989; Kendall et al., 1992; Woon and Dunning, 1993) and jul-cc-pVDZ (Papajak et al., 2011) double-ζ basis sets. For the PW6B95 functional, calculations were also carried out by using the jun-cc-pVDZ basis set (Papajak et al., 2011). The rev-DSD-PBEP86 double-hybrid density functional, recently proposed by Martin and coworkers (Santra et al., 2019), was employed together with the aug-cc-pVTZ and jun-cc-pVTZ basis sets. Indeed, triple-ζ basis sets in conjunction with the B2PLYP double-hybrid functional (Grimme, 2006) have been demonstrated to provide accurate predictions of geometries, rotational spectroscopic parameters and vibrational properties (Biczysko et al., 2010; Penocchio et al., 2015; Spada et al., 2017b; Tasinato et al., 2017; Boussessi et al., 2020a,b). Both PW6B95 and rev-DSD-PBEP86 were augmented for dispersion correlation through the Grimme's DFT-D3 scheme (Grimme et al., 2010) with Becke-Johnson damping (Grimme et al., 2011), even if the bare PW6B95 functional can already provide a satisfactory description of dispersion forces (Tasinato and Grimme, 2015). Since tight *d* functions are important for a quantitative representation of the electronic structure of second-row elements, partially augmented basis sets, namely aug-/jul-/jun-cc-pV(*n*+*d*)Z, including an additional set of *d* functions, were employed for sulfur and chlorine atoms. These basis sets were downloaded from the Basis Set Exchange library (Pritchard et al., 2019). At each level of theory, geometries were first optimized and then harmonic vibrational frequencies were computed by means of analytical Hessian matrices. While details for the calculation of analytical second-order derivatives of double-hybrid density functionals can be found in Biczysko et al. (2010), here it is only mentioned that their evaluation requires the full derivatives of the correlation contribution to the one-particle density matrix, γ^{x}. Its occupied-occupied and virtual-virtual blocks depend on the products of second-order perturbation amplitudes and amplitude derivatives, whereas the occupied-virtual block can be found from the solution of the so-called derivative *Z*-vector equations, that involve the derivatives of the MP2 Lagrangian. The cubic and semi-diagonal quartic force constants, and the second and third derivatives of the dipole moment surface were calculated by numerical differentiation of analytic quadratic force constants and dipole moment first derivatives, respectively. These quantities were then employed for the computation of rotational and vibrational spectroscopic parameters beyond the RR/HO approximation. In particular, quartic centrifugal distortion constants, vibrational contributions to rotational constants and vibrational frequencies were derived in the framework of second-order vibrational perturbation theory (Mills, 1972; Papoušek and Aliev, 1982; Barone, 2005; Bloino et al., 2012). In order to tackle the problem of resonances plaguing the anharmonic vibrational frequencies, the generalized second-order vibrational perturbation theory (GVPT2) was adopted, in which the (near-) singular terms are removed from the perturbative summations of anharmonicity constants and transition dipole moments (leading to the so called deperturbed approach, DVPT2) and the energy levels coupled by the resonances are treated in a second step by a proper variational calculation of reduced dimensionality (Mills, 1972; Martin et al., 1995; Barone, 2005). All DFT calculations were carried out by using the Gaussian 16 suite of programs (Frisch et al., 2016), which was also employed for the perturbative treatment of the anharmonic force field in the framework of a general VPT2 engine. The latter, in addition to performing the calculation of anharmonic vibrational energies according to various flavors of VPT2 (namely, DVPT2, GVPT2 and the so-called hybrid-degeneracy corrected PT2, HDCPT2), allows the computation of transition integrals for a number of spectroscopic techniques (IR, Raman, VCD), from which the corresponding anharmonic transition intensities can be derived (Bloino et al., 2012, 2015). In addition, it should be noted that recently the VPT2 framework has been coupled to a 1-dimensional discrete-variable-representation (1D-DVR) approach for the treatment of molecular systems presenting one large-amplitude vibration (Baiardi et al., 2017), which has been applied to the simulation of the IR spectrum of the methyl-cyclopropenyl cation (Puzzarini et al., 2019b). The 1D-DVR method is currently implemented in the development version of the Gaussian code and it will be included in the next releases of the software, yet it is not required for simulating the spectroscopic properties of the molecules considered in this work, all being semi-rigid systems. Since rev-DSD-PBEP86 is not among the Gaussian built-in functionals, it has been defined by setting proper IOP flags on top of the DSD-PBEP86 functional.

## 3. Results and Discussion

The performance of last-generation density functionals (DFs) for applications in the field of rotational and vibrational spectroscopy is investigated in relation to (a) equilibrium geometry, (b) rotational spectroscopic parameters, i.e., ground state rotational constants and quartic centrifugal distortion parameters, (c) harmonic vibrational frequencies and IR intensities, and (d) fundamental anharmonic wavenumbers and IR integrated absorption cross sections. The computed data are benchmarked against values determined experimentally and, in addition, they are also compared to high-level CCSD(T)-based results taken from the literature. Statistical indicators, such as mean deviation (MD), mean absolute deviation (MAD) and mean absolute percentage deviation (MAD%) are used to assess the accuracy of the different model chemistries. For comparison purposes, B3LYP/SNSD results as well as those obtained by the B2PLYP functional in conjunction with cc-pVTZ, aug-cc-pVTZ and may′-cc-pVTZ [i.e., may-cc-pVTZ (Papajak et al., 2009) with *d* functions removed from hydrogen atoms] basis sets, obtained in a previous work (Boussessi et al., 2020a), are also reported.

### 3.1. Equilibrium Geometry

Theoretical equilibrium geometries have been benchmarked against semi-experimental equilibrium structures that, in view of their high accuracy, represent optimal reference values to test the predictive power of new computational approaches. The MDs and MADSs obtained for bond lengths and bond angles over the whole set of molecules are reported in Figures 2A,B, respectively, while the results obtained for each molecule can be found in Supplementary Tables 1–20. All the hybrid DFs well describe the equilibrium geometries, with MADs within 0.01 Å and 0.35° for bond lengths and bond angles, respectively. Interestingly, both the ωB97 family and the PW6B95 functional yield more accurate equilibrium geometries than the B3LYP/SNSD model, that during the last years, has been proposed as a good tradeoff between accuracy and computational cost to predict the structure and spectroscopic properties of medium- to large-size molecules. However, it should be noted that for the functionals belonging to the ωB97 family this improvement can be mainly attributed to the use of additional polarization functions on second-row atoms. In fact, for difluoromethane, oxirane, glycolaldehyde and ethanimine bond lengths obtained from ωB97, ωB97X and ωB97X-D functionals are in line with B3LYP/SNSD ones. Significantly more accurate bond lengths are obtained for molecules containing sulfur and chlorine atoms: in particular, it should be noted that the description of the C-S bond length of ethyl-mercaptan improves by one order of magnitude, and for SO_{2} the deviation from the semi-experimental equilibrium value for the S=O bond length decreases from 0.05 Å with the SNSD basis set to about 0.01 Å employing the ωB97(X-D) functional in conjunction with the aug-jul-cc-pV(D+d)Z basis sets. Similar conclusions can be drawn for the S-H and C-Cl bond lengths thus highlighting the importance of additional *d* functions for second-row elements. The improvement brought by PW6B95 over B3LYP is more systematic, indeed, in conjunction with aug- and jul-cc-pVDZ it attains lower deviations also for molecules containing only first-row elements. Moving to double-hybrid functionals, Figures 2A,B, demonstrate the good performance of B2PLYP and rev-DSD-PBEP86 approximations that, when coupled to triple-zeta basis sets, deliver an equivalent description of bond lengths (with MADs around 0.0035 Å), whereas for bond angles the rev-DSD-PBEP86 functional appears more accurate than B2PLYP, with excellent MDs and MADs of 0.03 and 0.15°, respectively.

**Figure 2**. Mean deviation (MD) and mean absolute deviation (MAD) from semi-experimental equilibrium structural parameters for **(A)** bond lengths and **(B)** bond angles.

### 3.2. Rotational Spectroscopic Parameters

Rotational constants are the leading terms for the prediction of rotational spectra as they rule the frequencies of the rotational transitions. While the rotational constants of the equilibrium configuration are straightforwardly derived from the equilibrium structure, for applications in the field of rotational spectroscopy, vibrational effects must be included in order to obtain the rotational constants of the molecule in a given (usually the ground) vibrational state. Even though vibrational corrections usually account only for ~1–3% of the total rotational constant value, their inclusion is mandatory for quantitative predictions of rotational spectra and, furthermore, it is necessary for the comparison with experimentally determined rotational constants. As well-known, ground state rotational constants, (*B*_{0}), are obtained by adding vibrational corrections, Δ*B*_{vib} to equilibrium rotational constants (*B*_{e}):

where the summation runs over all the normal modes of vibration and the vibration-rotation interaction constants, ${\alpha}_{r}^{i}$, are evaluated in the framework of VPT2 (Mills, 1972). Their computation requires to go beyond the RR/HO approximation as they depend on the semi-diagonal cubic force field, i.e., the third derivatives of the potential energy. Besides vibrational effects, also centrifugal distortions need to be accounted for, especially for predicting high-energy rotational transitions which are of interest in e.g., astrophysical and atmospheric applications of rotational and high-resolution IR spectroscopy. Given the different orders of magnitude that rotational constants and centrifugal distortion parameters can have, the performance of the different levels of theory considered in this work can be more conveniently evaluated by referring to percentage deviations.

#### 3.2.1. Ground State Rotational Constants

Mean percentage deviations and mean absolute percentage deviations from experimental ground state rotational constants are reported in Figure 3A and the full list of results can be found in Supplementary Tables 21–30. At first, it should be noted that CCSD(T)-based computations, B3LYP/SNSD and B2PLYP/triple-ζ levels of theory reproduce experimental values with MAD%s around 0.8, 2.3, and 1%, respectively. However, the average errors for CCSD(T) computations may be slightly overestimated due to the inaccuracy of the rotational constants of cis-ClHC=CHF computed by a scaled-CCSD(T) force field, in particular *A*_{0} which displays a deviation of −12% with respect to the experimental value. Indeed, by discarding this molecule, the global MD% and MAD% for CCSD(T)-based methods drop to −0.1 and 0.3%, respectively. The last-generation DFs examined in the present work provide excellent results in the prediction of ground state rotational constants: indeed, the ωB97 family shows MD%s and MAD%s around 1% and the PW6B95 functional performs even better, the MAD% being around 0.6% when coupled to the jun-cc-pVDZ and jul-cc-pVDZ basis sets and 0.8% in conjunction with the aug-cc-pVTZ basis set. A similar accuracy is delivered by the rev-DSD-PBEP86 double-hybrid functional, whose MAD% of 0.7% slightly improves the performance of the B2PLYP functional.

**Figure 3**. Mean percentage deviation (MD%) and mean absolute percentage deviation (MAD%) from experimental **(A)** ground state rotational constants and **(B)** quartic centrifugal distortion constants.

#### 3.2.2. Centrifugal Distortion Constants

As shown in Figure 3B, quartic centrifugal distortion constants are reproduced by hybrid DFs with MAD%s around 6–7% independently of the functional or basis set employed (the full list of results is reported in Supplementary Tables 31–40). Despite the similar performance, the last generation DFs considered in the present work perform slightly worse than the B3LYP/SNSD model, with the exception of the PW6B95/jul-cc-VDZ level of theory, which yields almost the same MAD% as B3LYP (5.5%) and a lower MD% (−2.0 vs. 3.2%), accompanied however by larger fluctuations (from −46 to 38% for PW6B95/jul-cc-pVDZ, in the −29 −19% range for B3LYP/SNSD). The rev-DSD-PBEP86 functional reproduces the quartic centrifugal distortion constants determined experimentally with a MD% and a MAD% around −2.2 and 3.9%, respectively, which are very similar to the scores of the B2PLYP functional in conjunction with augmented triple-ζ basis sets (MD% = −1.6% and MAD% = 3.6%). In passing, it should be noted that removal of diffuse functions worsens the accuracy of the results by about 1%. At this point a few remarks concerning the comparison between experimental and theoretical centrifugal distortion constants are deserved. First, it has to be noted that computed constants refer to the equilibrium configuration of the molecule whereas rotational spectroscopy measurements provide those of the ground vibrational state. Even if the vibrational dependence of centrifugal distortion parameters is expected to amount to a few percent, it is an experimentally measurable quantity. However, only a few studies have been devoted to the theoretical treatment of vibrational effects on the centrifugal distortion (Watson, 2005). Second, attention has to be paid in comparing theory and experiment because measured centrifugal distortion constants can be affected by both limited accuracy and shortcomings in the fitting procedure used for their determination, as already pointed out by Boussessi et al. (2020a,b) Indeed, during the fitting procedure, centrifugal distortion constants may absorb the effects of resonances not fully treated (in particular high-order Coriolis or anharmonic interactions) that can be difficult to describe properly within the ro-vibrational Hamiltonian employed for inverting the measured transitions.

### 3.3. Harmonic Vibrational Properties

Since experimental harmonic vibrational frequencies are not available except for a few very simple molecules, DFT predictions are here benchmarked against CCSD(T) computations performed either in conjunction with large basis sets or within composite schemes. This comparison can be justified *a posteriori* on the basis of the good agreement between high-level CCSD(T) computations and experimental anharmonic wavenumbers and integrated absorption cross sections, which imply an accurate underlying harmonic force field. In the following discussion, the gauche conformer of ethyl mercaptan has been excluded from the data set in view of the huge discrepancies, up to 113 cm^{−1}, between CCSD(T)-F12 and experimental fundamental wavenumbers as discussed in (Boussessi et al., 2020a). MDs and MADs for harmonic vibrational wavenumbers and IR intensities are shown in Figures 4A,B, respectively, with the results for individual molecules being listed in the Supplementary Tables 41–59. As it can be seen from Figure 4A, concerning harmonic wavenumbers, among the functionals belonging to the ωB97 family there is a steady improvement of the performance on moving from ωB97, to ωB97X up to ωB97X-D which, with a MAD of around 14 cm^{−1} and a MD very close to zero, is the only one reaching an accuracy better than the B3 model (MAD = 16 cm^{−1}). The performance of the PW6B95 functional with different basis sets is similar to that of B3LYP in conjunction with the SNSD basis set, with MADs between 15 and 17 cm^{−1} and MDs around 0.7 cm^{−1}. On the other hand, the rev-DSD-PBEP86 functional represents a slight improvement over the already notable predictive power of the B2 model, being able to reproduce CCSD(T)-based reference data with a MAD of about 5 cm^{−1} to be compared with about 8 cm^{−1} at the B2PYLP level. A different picture is obtained for harmonic IR intensities: in fact, as shown in Figure 4B, all the hybrid functionals show comparable accuracy with MADs around 5 km mol^{−1}, the only exception being PW6B95 in conjunction with the jun-cc-pVDZ basis set, which provides poorer results. This is probably related to the lack of diffuse *d*-functions in the basis set, whose role in the computation of IR intensities is well-known: the MAD of 9 km mol^{−1} is indeed very similar to that obtained in Boussessi et al. (2020a) for the B3LYP functional in conjunction with the pcs-1 basis set, also lacking diffuse functions. Interestingly, also for intensities the rev-DSD-PBEP86 functional, with a MAD of 2.2 km mol^{−1} performs better than B2PLYP, whose MAD in conjunction with the aug-cc-pVTZ basis set is around 3 km mol^{−1} from reference CCSD(T)-based results. From the above discussion, it can be concluded that, concerning the calculation of vibrational frequencies and IR intensities within the double-harmonic approximation, the ωB97X-D and PW6B95 hybrid functionals can represent good alternatives to the B3LYP/SNSD level of theory, and rev-DSD-PBEP86 in conjunction with triple-ζ basis sets including diffuse functions appears even more accurate than the already well-performing B2 model.

**Figure 4**. Mean deviation (MD) and mean absolute deviation (MAD) from reference CCSD(T)-based theoretical values for **(A)** harmonic vibrational frequencies, **(B)** harmonic infrared intensities.

### 3.4. Beyond the Double-Harmonic Approximation: Wave-Numbers and Absorption Cross Sections

The performance of the DFs of the last generation considered in this work for anharmonic fundamental frequencies are compared in Figure 5 (Supplementary Tables 60–69 report the results for each molecule). It is noted that, although CCSD(T)-based computations reach a MAD of 8 cm^{−1}, this result can be biased by the disagreement between experimental transition frequencies and CCSD(T)-F12/cc-pVTZ-F12 predictions previously reported for gauche-CH_{3}CH_{2}SH. As to the B3 and B2 models, they reproduce experimental observations with a MAD of 20 and 12 cm^{−1}, respectively. In comparison, all the model chemistries based on the ωB97 approximation, showing MADs between 26 and 36 cm^{−1} cannot compete with the B3LYP functional. Given the similar results obtained for harmonic properties, it is evident that the ωB97 DFs have problems with the computation of higher order derivatives of the potential energy (i.e., cubic and semi-diagonal quartic force constants). Some of the anharmonic contributions computed for glycolaldehyde at the ωB97X-D level are particularly disappointing: the harmonic frequencies of ν_{6} and ν_{8} normal modes are predicted at about 1,445 and 1,305 cm^{−1}, respectively, whereas at the ωB97X-D/jul-cc-pVDZ anharmonic level they are shifted at 904 and 719 cm^{−1} (1,012 and 832 cm^{−1} when using the aug-cc-pVTZ basis set). Even worse is the case of the ν_{12} vibration for which the anharmonic corrections evaluated by using jul-cc-pVDZ and aug-cc-pVDZ basis sets amount to −1,005 and −503 cm^{−1} thus resulting in a non-physical negative value of the ν_{12} fundamental transition frequency. Conversely, the PW6B95 DF turns out to be competitive with the B3 model, reaching a comparable MAD of 18 cm^{−1} when employed in conjunction with jul- and jun-cc-pVDZ basis sets and of 20 cm^{−1} together with the aug-cc-pVTZ basis set. Furthermore, recently Kreienborg and Merten (2019) pointed out that carbon-fluorine stretching vibrations are often strongly misplaced by common hybrid functionals with the possible exception of the M06-2X functional (restricting, however, the simulation to the double-harmonic approximation). Nevertheless, while this functional (and its predecessor M05-2X) predicts harmonic frequencies not far from experimental fundamentals, it becomes unreliable when anharmonic contributions are taken into the proper account (Puzzarini et al., 2010; Tasinato, 2014; Tasinato et al., 2018). Rather, it should be noted that the PW6B95 functional yields a consistent description of anharmonic C-F stretching frequencies, with deviations from experiment generally halved with respect to those of the B3 model. In fact, by focusing only on the C-F stretchings, the B3LYP/SNSD level of theory presents a MAD (computed over seven values) of 34 cm^{−1} and a maximum deviation of −50 cm^{−1}, while the PW6B95 DF (in conjunction with all the tested double-ζ basis sets) reproduces the experimental values with a MAD around 15 cm^{−1} and a maximum error of about −25 cm^{−1}. The rev-DSD-PBEP86 functional slightly improves over B2PLYP also for anharmonic fundamental wavenumbers: indeed it reproduces experimental data with a MD and a MAD around 2 and 8 cm^{−1}, respectively, in comparison to −1 and 12 cm^{−1} for the B2PLYP/aug-cc-pVTZ level of theory.

**Figure 5**. Mean deviation (MD) and mean absolute deviation (MAD) from experimental fundamental frequencies. MADs for the ωB97 family are out of range.

Some remarks are in order about the accuracy in the computation of IR integrated absorption cross sections (i.e., IR band intensities). It should be noted that their experimental determination is a daunting task prone to both systematic and random errors, hence a careful control of the experimental conditions and errors source needs to be performed during the measurements. Therefore, to assess the accuracy of DFT calculations, one must rely on precise experimental values, which are available only for a reduced number of molecules among those of the benchmark set, namely CH_{2}F_{2}, CH_{2}ClF, ClFC=CH_{2}, ClFC=CF_{2}, and SO_{2}. Moreover, in view of the poor results delivered by the ωB97-based functionals for anharmonic wave-numbers, the attention for IR anharmonic intensities is focused on PW6B95 and rev-DSD-PBEP86 DFs, whose performances are shown in Figure 6 together with those of B3 and B2 models (the full list of results is given in Supplementary Tables 70–74). In passing it should be stressed that the accuracy reached for integrated absorption cross sections depends on the reliability of both the anharmonic potential energy and dipole moment surfaces. In fact, when overlapping IR bands cannot be resolved at the experimental level, the integration required for determining band intensities is performed by considering all the absorptions within a given spectral interval. Here, the same approach has been mimicked at the theoretical level, i.e., the intensities of all the bands predicted in a given (experimental) integration range have been summed. As it can be seen, as for anharmonic wavenumbers, PW6B95 in conjunction with the aug-cc-pVDZ or jul-cc-pVDZ basis sets performs on par with B3LYP/SNSD, the MDs and MADs being around 9 and 10 km mol^{−1}, respectively. Furthermore, as already pointed out for harmonic IR intensities, the lack of diffuse *d* functions in the jun-cc-pVDZ basis set results in a worsening of the predictions, thus highlighting the necessity of diffuse polarization functions for a reliable description of the dipole moment surface. Moving to the double-hybrid DFs, Figure 6 shows that, also for integrated absorption cross sections, the rev-DSD-PBEP86 functional, with both MD and MAD around 5 km mol^{−1}, improves over B2PLYP that, in turn, reproduces experimental measurements with a mean absolute deviation ranging between 8 and 10 km mol^{−1} for the different triple-ζ basis sets.

**Figure 6**. Mean deviation (MD) and mean absolute deviation (MAD) from experimental integrated absorption cross sections for selected model chemistries.

Finally, anharmonic vibrational properties have been computed by a hybrid QM/QM′ approach featuring harmonic frequencies and intensities calculated at the rev-DSD-PBEP86/jul-cc-pVTZ level and anharmonic effects obtained from PW6B95/jun-cc-pVDZ computations. The hybrid force field obtained in this way simulates fundamental wavenumbers with a MAD of 11 cm^{−1}, thus showing a worsening of only 3 cm^{−1} in comparison to full rev-DSD-PBEP86 computations, but still performing similarly to full B2 anharmonic predictions. This result shows that the PW6B95/jun-cc-pVDZ level of theory represents a reliable and cost-effective approach to compute anharmonic corrections when employed in conjunction with harmonic force fields obtained from the rev-DSD-PBEP86 functional. Conversely, PW6B95 should be used together with the larger jul-cc-pVDZ basis set when performing full anharmonic computations in order to obtain reliable predictions of band intensities.

### 3.5. New Model Chemistries at Work

The previous sections have shown that, among the hybrid DFs of the last generation considered in this work, ωB97X-D and PW6B95 in conjunction with either aug- or jul-cc-pVDZ basis sets provide reliable predictions of equilibrium structures, rotational parameters and harmonic vibrational properties, sometimes even better than the well-tested B3 model. However, when anharmonic effects come into play, all the functionals of the ωB97 family yield unstable, sometimes disappointing, results. Conversely, the PW6B95 functional performs similarly to the B3 model for both fundamental transition frequencie and integrated absorption cross sections, provided that an additional set of *d* functions is employed on heavy elements, and in the case of intensities, at least a basis set of the jul-cc-pVDZ quality is used. Coming to double-hybrid functionals, the predictions of the rev-DSD-PBEP86 model are better than, or at least similar to, their B2 counterparts for both structural and rotational-vibrational spectroscopic properties.Given these premises, in this section the model chemistries able of rivaling with the B2 and B3 paradigms are applied to selected case studies in the field of structural determination and IR spectroscopy. A convenient case-study for the first subject is the complex formed by SO_{2} and dimethyl sulfide, (CH_{3})_{2}S (see Figure 7A), whose semi-experimental equilibrium structure has been recently determined (Obenchain et al., 2018). Then, the anharmonic IR spectrum of the cyclic-(CH)C_{3}${\text{H}}_{2}^{+}$ cation (Figure 7C) has been simulated because of the astrophysical interest of this molecule pointed out very recently by Westbook et al. (2020).

**Figure 7**. Structures of **(A)** SO_{2}⋯ (CH_{3})_{2}S complex, **(B)** dimethyl sulfide, and **(C)** cyc-(CH)C_{3}H_{2} cation.

The structural parameters of the equilibrium configuration of the SO_{2}⋯ (CH_{3})_{2}S 1:1 complex evaluated by different DFs and basis sets are compared in Table 1 with the semi-experimental structure obtained by Obenchain et al. (2018). Inspection of this table reveals that the B3 model reproduces bond lengths, valence and dihedral angles with average absolute errors of 0.03 Å, 0.6 and 4°, respectively. The PW6B95-D3/jul-cc-pV(D+*d*)Z model chemistry delivers significantly improved results, almost halving the deviations for bond lengths and valence angles and showing an absolute average error for dihedral angles around 0.6°. Concerning the B2PLYP-D3 (employed in conjunction with the may′-cc-pVTZ basis set) and rev-DSD-PBEP86-D3 (in conjunction with the jun-cc-pV(T+*d*) basis set) functionals, it can be seen that, on average, they perform equally well for both bond lengths (the mean absolute deviation over all the complex bond lengths are 0.014 and 0.013 Å, respectively) valence (0.2 vs. 0.1°) and dihedral (0.8 and 0.7°) angles. It is also noteworthy that the hybrid PW6B95-D3 functional reproduces the semi-experimental structure of the complex with an accuracy that, in spite of the considerably reduced computational cost, rivals that of the B2 and rev-DSD-PBEP86 double-hybrids. In order to understand whether the quite large deviation observed for the inter-molecular S-S distance is due to an intrinsic inaccuracy for S-S bonds or to an unbalanced treatment of inter-molecular interactions, the equilibrium geometry of hydrogen disulfide, HSSH (Figure 7B), has been computed at the PW6B95-D3/jul-cc-pV(D+*d*)Z level and compared with the semi-experimental structure recently reported by Ye et al. (2020). The obtained structural parameters, detailed in Table 2, show deviations of 8 and 9 mÅ for the S-S and S-H bond lengths, respectively, and of 0.3° for the HSS angle, while the HSSH dihedral is within the uncertainty of the semi-experimental value. Since the SS bond length of HSSH computed at PW6B95-D3 and rev-DSD-PBEP86-D3 levels is very close, the worse performance of the former functional for the SS distance in the complex is probably related to the description of non-covalent interactions.

Moving to the vibrational properties of the cyc-(CH)C_{3}${\text{H}}_{2}^{+}$ cation, harmonic frequencies are listed in Table 3 while anharmonic fundamental wavenumbers and IR intensities are reported in Table 4. In both Tables the results obtained by Westbrook et al. at the CCSD(T)-F12/aug-cc-pVTZ level (Westbook et al., 2020) are also reported for comparison purposes. In that work, the authors focused on the difficulties of post-Hartree-Fock methods in describing out-of-plane bending vibrations of molecules with C=C multiple bonds (especially, but not only, aromatic systems), an issue that should be of minor concern for DFT. For this reason, it is of some interest to compare their CCSD(T)-F12 predictions with the present DFT simulations. For hybrid DFs, the major differences with respect to CCSD(T)-F12/aug-cc-pVTZ harmonic vibrational frequencies, can be observed, for both B3LYP and PW6B95 functionals, for the ω_{7}, ω_{9}, ω_{14} and ω_{15} normal modes which correspond to the H-C=C-H out-of-plane bending, the ≡C-H out-of-plane vibration, the in-plane ring deformation and the in-plane bending vibration of the C-H groups, respectively. Concerning double-hybrids, the largest difference (+36 cm^{−1}) is observed for ω_{7} at B2PLYP/jun-cc-pVTZ level, while the worst rev-DSDPBEP86/jun-cc-pVTZ result concerns the ω_{3} C≡C stretching vibration (+27 cm^{−1}). Moving to the anharmonic vibrational wavenumbers, a huge anharmonic correction (−202 cm^{−1}) for the ν_{10} fundamental and an unusual positive contribution (20 cm^{−1}) for the ν_{15} vibration were obtained by Westbook et al. (2020) Conversely, according to the present calculations, the anharmonic correction for the ν_{15} vibration amounts to a few wave-numbers and that of ν_{10} ranges from about −45 cm^{−1} at the B3LYP level to −2 cm^{−1} at the B2PLYP/jun-cc-pVTZ level. While a number of Fermi resonances both of type 1 and 2 have been found, none of them strongly alters the spectral structure, just shifting the transitions by about 5 cm^{−1} from their unperturbed values. The only exception is the ν_{6} normal mode, whose excited *v*_{6} = 1 level is involved in a resonant triad together with *v*_{10} = 2 and *v*_{14} = *v*_{15} = 1.

**Table 3**. Harmonic frequencies (in cm^{−1}) and intensities (in km mol^{−1}) for *cyc*-CHC_{3}${\text{H}}_{2}^{+}$ computed at different levels of theory.

**Table 4**. Anharmonic frequencies (in cm^{−1}) and intensities (in km mol^{−1}) for *cyc*-CHC_{3}${\text{H}}_{2}^{+}$ obtained at different levels of theory.

Three different hybrid force fields have been also computed: in the B2/B3 hybrid force field, B2PLYP/jun-cc-pVTZ harmonic frequencies have been corrected by B3LYP/SNSD anharmonic contributions, while in the revDSD/B3 and revDSD/PW6 hybrid QM/QM′ approaches, the harmonic force field at the rev-DSD-PBEP86/jun-cc-pVTZ level has been coupled to cubic- and semi-diagonal quartic force constants evaluated at B3LYP/SNSD and PW6B95/jul-cc-pVDZ levels, respectively. The results collected in Table 4 show that the hybrid force fields deliver, with a few exceptions, the same results as the corresponding force field obtained by full anharmonic calculations at B2PLYP or rev-DSD-PBEP86 levels. The spectra of cyc-(CH)C_{3}${\text{H}}_{2}^{+}$ cation, simulated beyond the double-harmonic approximation at different levels of theory, which can be useful to guide experimental measurements on this species or even astronomical campaigns, are reported in Figure 8. It is quite apparent that the only major differences among the different theoretical models are found for the integrated absorption cross sections, whereas transition frequencies are only marginally affected.

**Figure 8**. Anharmonic infrared spectrum of cyc-CHC_{3}${\text{H}}_{2}^{+}$ simulated by using different levels of theory (computed spectral transitions have been convoluted with a Lorentzian function with an half-width at half-maximum of 2 cm^{−1}). B2: full B2PLYP/jun-cc-pVTZ anharmonic force field; B2/B3: B2PLYP/jun-cc-pVTZ harmonic force field and B3LYP/SNSD anharmonic effects; rev-DSD: full rev-DSDPBEP86/jun-cc-pVTZ anharmonic force field; rev-DSD/PW6: rev-DSDPBEP86/jun-cc-pVTZ harmonic force field and PW6B95/jul-cc-pVDZ anharmonic effects. Some traces have been displaced fpr clarity.

## 4. Conclusions

In the last decade we have witnessed the increasing accuracy of methods rooted in the density functional theory due to ongoing improvements in both methodological and numerical aspects. However, attention has mainly been focused on thermochemistry and kinetics, whereas theoretical support to rotational and vibrational spectroscopy requires accurate predictions of molecular geometries, harmonic force fields and leading anharmonic contributions, not to mention dipole moments and their derivatives. On these grounds, the present work has analyzed the performance of last-generation hybrid and double-hybrid functionals in conjunction with partially augmented correlation consistent basis sets for a benchmark set of 10 molecules of both atmospheric and astrochemical relevance. Equilibrium molecular geometries issued from DFT computations have been benchmarked against accurate semi-experimental equilibrium structures, while rotational constants, centrifugal distortion parameters, and vibrational frequencies have been compared to the experimental data available in the literature and with high level CCSD(T)-based *ab* *initio* calculations. The following conclusions can be drawn:(1) The jun-cc-pVDZ basis set performs remarkably well for hybrid functionals with the possible exception of IR intensities, which require diffuse *d*-functions, namely the jul-cc-pVDZ (or SNSD) basis set. In the case of double-hybrid functionals, the jun-cc-pVTZ basis set represents a nearly optimum cost/performance balance, but also the may−-cc-pVTZ basis set can be safely employed for larger systems. An additional set of *d*-functions is always mandatory for second-row atoms.(2) Among hybrid functionals, B3LYP-D3 (B3) is still very competitive, although the PW6B95 (PW6) model significantly improves equilibrium geometries.(3) Concerning double-hybrid functionals, the rev-DSD-PBEP86-D3 functional (rDSD) systematically improves the already reliable results delivered by the B2PLYP (B2) model, the enhancement being especially significant for non-covalent complexes. (4) Composite methods employing geometries and harmonic contributions evaluated by double-hybrid functionals coupled to anharmonic corrections, obtained with hybrid functionals, always lead to accurate results. In this connection the previously employed B2/B3 model remains very useful, but the new rDSD/PW6 variant seems capable of delivering even better results with the same cost.

In summary, even with further pending developments and validation, thanks to effective implementations in general electronic structure codes, last-generation hybrid and double-hybrid functionals provide unprecedented accuracy for all the parameters ruling rotational and vibrational spectroscopy with computer requirements well within current standards and, coupled to generalized second-order vibrational perturbation theory (GVPT2), can also be used by non-specialists to complement experimental studies of medium- and, even, large-size molecules of current fundamental and technological interest.

## Data Availability Statement

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

## Author Contributions

GC and MF performed the computations and analyzed the results. NT supervised the work, analyzed the results, and wrote part of the paper. VB defined the general strategy and wrote part of the paper. All authors contributed to the article and approved the submitted version.

## Funding

This work has been supported by MIUR PRIN 2015 (Grant Number 2015F59J3R), PRIN 2017 (Grant Number 2017A4XRCA), and by Scuola Nornale Superiore (Grant Number SNS18B-TASINATO).

## 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 technical staff of the SMART laboratory at Scuola Normale Superiore di Pisa is thanked for help with use of HPC resources.

## Supplementary Material

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

## References

Aliev, M. R., and Watson, J. K. G. (1985). “Higher-order effects in the vibration-rotation spectra of semirigid molecules,” in *Molecular Spectroscopy: Modern Research*, Vol. 3, ed K. N. Rao (New York, NY: Academic Press), 2–67. doi: 10.1016/B978-0-12-580643-5.50006-3

Alonso, J. L., Lesarri, A., Leal, L. A., and López, J. C. (1993). The millimeter-wave spectra of 1-chloro-1-fluoroethylene and cis-1-chloro-2-fluoroethylene. *J. Mol. Spectrosc.* 162, 4–9. doi: 10.1006/jmsp.1993.1265

Baiano, C., Lupi, J., Tasinto, N., Puzzarini, C., and Barone, V. (2020). The role of state-of-the-art quantum-chemical calculations in astrochemistry: formation route and spectroscopy of ethanimine as a paradigmatic case. *Molecules* 25:2873. doi: 10.3390/molecules25122873

Baiardi, A., Bloino, J., and Barone, V. (2017). Simulation of vibronic spectra of flexible systems: hybrid DVR-harmonic approaches. *J. Chem. Theory Comput.* 13, 2804–2822. doi: 10.1021/acs.jctc.7b00236

Barone, V. (2005). Anharmonic vibrational properties by a fully automated second-order perturbative approach. *J. Chem. Phys.* 122:014108. doi: 10.1063/1.1824881

Barone, V. (2012). *Computational Strategies for Spectroscopy: From Small Molecules to Nano Systems*. Hoboken, NJ: John Wiley & Sons.

Barone, V. (2016). The virtual multifrequency spectrometer: a new paradigm for spectroscopy. *WIREs Comput. Mol. Sci.* 6, 86–110. doi: 10.1002/wcms.1238

Barone, V., Baiardi, A., Biczysko, M., Bloino, J., Cappelli, C., and Lipparini, F. (2012). Implementation and validation of a multi-purpose virtual spectrometer for large systems in complex environments. *Phys. Chem. Chem. Phys.* 14, 12404–12422. doi: 10.1039/c2cp41006k

Becke, A. D. (1993). Density functional thermochemistry. III. The role of exact exchange. *J. Chem. Phys.* 98, 5648–5652. doi: 10.1063/1.464913

Berova, N., Polavarapu, P. L., Nakanishi, K., and Woody, R. W. (2012). *Comprehensive Chiroptical Spectroscopy: Instrumentation, Methodologies and Theoretical Simulations*. Hoboken, NJ: John Wiley & Sons.

Biczysko, M., 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. Theory Comput.* 6, 2115–2125. doi: 10.1021/ct100212p

Blanco, S., Lesarri, A., Alonso, J. L., and Guarnieri, A. (1995). The rotational spectrum of chlorofluoromethane. *J. Mol. Spectrosc.* 174, 397–416. doi: 10.1006/jmsp.1995.0011

Bloino, J., Biczysko, M., and Barone, V. (2012). General perturbative approach for spectroscopy, thermodynamics, and kinetics: methodological background and benchmark studies. *J. Chem. Theory Comput.* 8, 1015–1036. doi: 10.1021/ct200814m

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

Boussessi, R., Ceselin, G., Tasinato, N., and Barone, V. (2020a). DFT meets the segmented polarization consistent basis sets: Performances in the computation of molecular structures, rotational and vibrational spectroscopic properties. *J. Mol. Struct.* 1208:127886. doi: 10.1016/j.molstruc.2020.127886

Boussessi, R., Tasinato, N., Pietropolli Charmet, A., Stoppa, P., and Barone, V. (2020b). Sextic centrifugal distortion constants: interplay of density functional and basis set for accurate yet feasible computations. *Mol. Phys.* 118:e1734678. doi: 10.1080/00268976.2020.1734678

Carlotti, M., Nivellini, G., Tullini, F., and Carli, B. (1988). The far-infrared spectrum of methylene fluoride. *J. Mol. Spectrosc.* 132, 158–165. doi: 10.1016/0022-2852(88)90065-3

Carroll, P. B., Drouin, B. J., and Widicus-Weaver, S. L. (2010). The submillimeter spectrum of glycolaldehyde. *Astrophys. J.* 723, 845–849. doi: 10.1088/0004-637X/723/1/845

Chai, J.-D., and Head-Gordon, M. (2008a). Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. *Phys. Chem. Chem. Phys.* 10, 6615–6620. doi: 10.1039/b810189b

Chai, J.-D., and Head-Gordon, M. (2008b). Systematic optimization of long-range corrected hybrid density functionals. *J. Chem. Phys.* 128:084106. doi: 10.1063/1.2834918

Craig, N. C., Lo, Y. S., Piper, L. G., and Wheeler, J. C. (1970). Vibrational assignments and potential constants for cis- and trans-1-chloro-2-fluoroethylenes and their deuterated modifications. *J. Phys. Chem.* 74, 1712–1727. doi: 10.1021/j100720a010

Drummen, G. P. C. (2012). Fluorescent probes and fluorescence (microscopy) techniques-illuminating biological and biomedical research. *Molecules* 17, 14067–14090. doi: 10.3390/molecules171214067

Dunning, T. H. J. (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

Flaud, J.-M., Lafferty, W. J., Kwabia-Tchana, F., Perrin, A. M., and Landsheere, X. (2012). First high-resolution analysis of the ν_{1}5, ν_{1}2, ν_{5}, ν_{1}0 and ν_{2} bands of oxirane. *J. Mol. Spectrosc.* 271, 38–43. doi: 10.1016/j.jms.2011.11.005

Flaud, J.-M., Perrin, A. M., Salah, L. M., Lafferty, W. J., and Guelachvili, G. (1993). A reanalysis of the (010), (020), (100), and (001) rotational levels of ^{3}2S^{1}6O_{2}. *J. Mol. Spectrosc.* 160, 272–278. doi: 10.1006/jmsp.1993.1174

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

Gambi, A., Pietropolli Charmet, A., Stoppa, P., Tasinato, N., Ceselin, G., and Barone, V. (2019). Molecular synthons for accurate structural determinations: the equilibrium geometry of 1-chloro-1-fluoroethene. *Phys. Chem. Chem. Phys.* 21, 3615–3625. doi: 10.1039/C8CP04888F

Gambi, A., Puzzarini, C., Cazzoli, G., Dore, L., and Palmieri, P. (2002). The anharmonic force field of *cis*-1-chloro-2-fluoroethylene. *Mol. Phys.* 100, 3535–3543. doi: 10.1080/00268970210130155

Grimme, S. (2006). Semiempirical hybrid density functional with perturbative second-order correlation. *J. Chem. Phys.* 124:034108. doi: 10.1063/1.2148954

Grimme, S., Anthony, 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

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

He, Y., Tang, L., Wu, X., Hou, X., and Lee, Y. (2007). Spectroscopy: the best way toward green analytical chemistry? *Appl. Spectrosc. Rev.* 42, 119–138. doi: 10.1080/05704920601184259

Hillig, K. W., Bittner, E. R., Kuczkowski, R. L., Lewis-Bevan, W., and Gerry, M. C. L. (1988). The chlorine nuclear quadrupole coupling tensor in chlorotrifluoroethylene. *J. Mol. Spectrosc.* 132, 369–379. doi: 10.1016/0022-2852(88)90332-3

Hochlaf, M., Puzzarini, C., and Senent, M. L. (2015). Towards the computations of accurate spectroscopic parameters and vibrational spectra for organic compounds. *Mol. Phys.* 113, 1661–1673. doi: 10.1080/00268976.2014.1003986

Johnson, T. J., Sams, R. L., Profeta, L. T. M., Akagi, S. K., Burling, I. R., Yokelson, R. J., et al. (2013). Quantitative IR spectrum and vibrational assignments for glycolaldehyde vapor: glycolaldehyde measurements in biomass burning plumes. *J. Phys. Chem. A* 117, 4096–4107. doi: 10.1021/jp311945p

Kendall, R. A., Dunning, T. H. J., 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

Kolesnikov, L., Tercero, B., Cernicharo, J., Alonso, J. L., Daly, A. M., Gordon, B. P., et al. (2014). Spectroscopic characterization and detection of ethyl mercaptan in orion. *Astrophys. J. Lett.* 784:L7. doi: 10.1088/2041-8205/784/1/L7

Kreienborg, N. M., and Merten, C. (2019). How to treat C-F stretching vibrations? A vibrational CD study on chiral fluorinated molecules. *Phys. Chem. Chem. Phys.* 21, 3506–3511. doi: 10.1039/C8CP02395F

Lafferty, W. J., Flaud, J. M., Tchana, F. K., and Fernandez, J. M. (2013). Raman and infrared spectra of the ν_{1} band of oxirane. *Mol. Phys.* 111, 1983–1986. doi: 10.1080/00268976.2013.775516

Lee, C., Yang, W., and Parr, R. G. (1988). Development of the colle-salvetti correlation-energy formula into a functional of the electron density. *Phys. Rev. B* 37, 785–789. doi: 10.1103/PhysRevB.37.785

Leung, H. O., Marshall, M. D., Vasta, A. L., and Craig, N. C. (2009). Microwave spectra of eight isotopic modifications of 1-chloro-1-fluoroethylene. *J. Mol. Spectrosc.* 253, 116–121. doi: 10.1016/j.jms.2008.11.002

Li, W., Spada, L., Tasinato, N., Rampino, S., Evangelisti, L., Gualandi, A., et al. (2018). Theory meets experiment for noncovalent complexes: the puzzling case of pnicogen interactions. *Angew. Chem. Int. Ed.* 57, 13853–13857. doi: 10.1002/anie.201807751

Licari, D., Baiardi, A., Egidi, F., Latouche, C., and Barone, V. (2015). Implementation of a graphical user interface for the virtual multifrequency spectrometer: the VMS draw tool. *J. Comput. Chem.* 36, 321–334. doi: 10.1002/jcc.23785

Licari, D., Tasinato, N., Spada, L., Puzzarini, C., and Barone, V. (2017). VMS-ROT: a new module of the virtual multifrequency spectrometer for simulation, interpretation, and fitting of rotational spectra. *J. Chem. Theory Comput.* 13, 4382–4396. doi: 10.1021/acs.jctc.7b00533

Lindon, J., Tranter, G. E., and Koppenaal, D. (2017). *Encyclopedia of Spectroscopy and Spectrometry*. Oxford: Academic Press.

Martin, J. M. L., Lee, T. J., Taylor, P. R., and Françoise, J.-P. (1995). The anharmonic force field of ethylene, C_{2}H_{4}, by means of accurate *ab initio* calculations. *J. Chem. Phys.* 103:2589. doi: 10.1063/1.469681

Medcraft, C., Thompson, C. D., Robertson, E. G., Appadoo, D. R. T., and McNaughton, D. (2012). The far-infrared rotational spectrum of ethylene oxide. *Astrophys. J.* 753:18. doi: 10.1088/0004-637X/753/1/18

Melli, A., Melosso, M., Tasinato, N., Bosi, G., Spada, L., Bloino, J., et al. (2018). Rotational and infrared spectroscopy of ethanimine: a route toward its astrophysical and planetary detection. *Astrophys. J.* 855:123. doi: 10.3847/1538-4357/aaa899

Miller, B. J., Howard, D. L., Lane, J. R., Kjaergaard, H. G., Dunn, M. E., and Vaida, V. (2009). SH-stretching vibrational spectra of ethanethiol and *tert*-butylthiol. *J. Phys. Chem. A* 113, 7576–7583. doi: 10.1021/jp9017162

Mills, I. M. (1972). “Vibration-rotation structure in asymmetric- and symmetric-top molecules,” in *Molecular Spectroscopy: Modern Research*, eds K. N. Rao and C. W. Mathews (New York, NY: Academic Press), 115–140. doi: 10.1016/B978-0-12-580640-4.50013-3

Mller, H. S., and Brnken, S. (2005). Accurate rotational spectroscopy of sulfur dioxide, SO_{2}, in its ground vibrational and first excited bending states, *v*_{2} = 0, 1, up to 2 thz. *J. Mol. Spectrosc.* 232, 213–222. doi: 10.1016/j.jms.2005.04.010

Obenchain, D. A., Spada, L., Alessandrini, S., Rampino, S., Herbers, S., Tasinato, D. N., et al. (2018). Unveiling the sulfur-sulfur bridge: Accurate structural and energetic characterization of a homochalcogen intermolecular bond. *Angew. Chem. Int. Ed.* 57, 15822–15826. doi: 10.1002/anie.201810637

Papajak, E., Leverentz, H. R., Zheng, J., and Truhlar, D. G. (2009). Efficient diffuse basis sets: cc-pVxZ+ and maug-cc-pVxZ. *J. Chem. Theory Comput.* 5, 1197–1202. doi: 10.1021/ct800575z

Papajak, E., Zheng, J., Xu, X. R., Leverentz, H., and Truhlar, D. G. (2011). Perspectives on basis sets beautiful: seasonal plantings of diffuse basis functions. *J. Chem. Theory Comput.* 7, 3027–3034. doi: 10.1021/ct200106a

Papoušek, D., and Aliev, M. R. (1982). *Molecular Vibrational/Rotational Spectra*. Amsterdam: Elsevier.

Penocchio, E., Piccardo, M., and Barone, V. (2015). Semiexperimental equilibrium structures for building blocks of organic and biological molecules: the B2PLYP route. *J. Chem. Theory Comput.* 11, 4342–4363. doi: 10.1021/acs.jctc.5b00622

Peveratti, R. (2020). Fitting elephants in the density functional zoo: statistical criteria for the evaluation of density functional theory methods as a suitable replacement for counting parameters. *Int. J. Quant. Chem.* 120, 1–15. doi: 10.1002/qua.26379

Piccardo, M., Bloino, J., and Barone, V. (2015a). Generalized vibrational perturbation theory for rotovibrational energies of linear, symmetric, and asymmetric tops: theory, approximations, and authomated approaches to deal with medium-to-large molecular systems. *Int. J. Quant. Chem.* 115, 948–982. doi: 10.1002/qua.24931

Piccardo, M., Penocchio, E., Puzzarini, C., Biczysko, M., and Barone, V. (2015b). Semi-experimental equilibrium structure determinations by employing B3LYP/SNSD anharmonic force fields: validation and application to semirigid organic molecules. *J. Phys. Chem. A* 119, 2058–2082. doi: 10.1021/jp511432m

Pietropolli Charmet, A., Stoppa, P., Tasinato, N., and Giorgianni, S. (2017). Computing sextic centrifugal distortion constants by DFT: a benchmark analysis on halogenated compounds. *J. Mol. Spectrosc.* 335, 117–125. doi: 10.1016/j.jms.2017.02.006

Pietropolli Charmet, A., Stoppa, P., Tasinato, N., Giorgianni, S., Barone, V., Biczysko, M., et al. (2013). An integrated experimental and quantum-chemical investigation on the vibrational spectra of chlorofluoromethane. *J. Chem. Phys.* 139:164302. doi: 10.1063/1.4825380

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

Presti, D., Pedone, A., Licari, D., and Barone, V. (2017). A modular implementation for the simulation of 1D and 2D solid-state NMR spectra of quadrupolar nuclei in the virtual multifrequency spectrometer-draw graphical interface. *J. Chem. Theory Comput.* 13, 2215–2229. doi: 10.1021/acs.jctc.7b00154

Pritchard, B. P., Altarawy, D., Didier, B., Gibson, T. D., and Windus, T. L. (2019). New basis set exchange: an open, up-to-date resource for the molecular sciences community. *J. Chem. Inform. Model.* 59, 4814–4820. doi: 10.1021/acs.jcim.9b00725

Puzzarini, C., and Barone, V. (2018). Diving for accurate structures in the ocean of molecular systems with the help of spectroscopy and quantum chemistry. *Acc. Chem. Res.* 51, 548–556. doi: 10.1021/acs.accounts.7b00603

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. Theory Comput.* 6, 828–838. doi: 10.1021/ct900594h

Puzzarini, C., Biczysko, M., Bloino, J., and Barone, V. (2014a). Accurate spectroscopic characterization of oxirane: a valuable route to its identification in titan's atmosphere and the assignment of unidentified infrared bands. *Astrophys. J.* 785:107. doi: 10.1088/0004-637X/785/2/107

Puzzarini, C., Bloino, J., Tasinato, N., and Barone, V. (2019a). 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

Puzzarini, C., Senent, M. L., Domnguez-Gme, R., Carvajal, M., Hochlaf, M., and Al-Mogren, M. M. (2014b). Accurate spectroscopic characterization of ethyl mercaptan and dimethyl sulfide isotopologues: a route toward their astrophysical detection. *Astrophys. J.* 796:50. doi: 10.1088/0004-637X/796/1/50

Puzzarini, C., Tasinato, N., Bloino, J., Spada, L., and Barone, V. (2019b). State-of-the-art computation of the rotational and IR spectra of the methyl-cyclopropyl cation: hints on its detection in space. *Phys. Chem. Chem. Phys.* 21, 3615–3625. doi: 10.1039/C8CP04629H

Quack, M., and Merkt, F. (2011). *Handbook of High-Resolution Molecular Spectroscopy*. Chichester: John Wiley & Sons.

Russell, D. K., and Wesendrup, R. (2003). Tunable diode laser study of the ν_{3} band of oxirane. *J. Mol. Spectrosc.* 217, 59–71. doi: 10.1016/S0022-2852(02)00010-3

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

Smith, D., Devlin, J. P., and Scott, D. W. (1968). Conformational analysis of ethanethiol and 2-propanethiol. *J. Mol. Spectrosc.* 25, 174–184. doi: 10.1016/0022-2852(68)80004-9

Spada, L., Tasinato, N., Bosi, G., Vazart, F., Barone, V., and Puzzarini, C. (2017a). On the competition between weak OH⋯ F and CH⋯ F hydrogen bonds, in cooperation with CH⋯ O contacts, in the difluoromethane–tert-butyl alcohol cluster. *J. Mol. Spectrosc.* 337, 90–95. doi: 10.1016/j.jms.2017.04.001

Spada, L., Tasinato, N., Vazart, F., Barone, V., Caminati, W., and Puzzarini, C. (2017b). Noncovalent interactions and internal dynamics in pyridine-ammonia: a combined quantum-chemical and microwave spectroscopy study. *Chem. Eur. J.* 23, 4876–4883. doi: 10.1002/chem.201606014

Sugiki, T., Kobayashi, N., and Fujiwara, T. (2017). Modern technologies of solution nuclear magnetic resonance spectroscopy for three-dimensional structure determination of proteins open avenues for life scientists. *Comput. Struct. Biotechnol. J.* 15, 328–329. doi: 10.1016/j.csbj.2017.04.001

Tasinato, N. (2014). What are the spectroscopic properties of HFC 32? Answers from DFT. *Int. J. Quant. Chem.* 114, 1472–1485. doi: 10.1002/qua.24716

Tasinato, N., Ceselin, G., Stoppa, P., Pietropolli Charmet, A., and Giorgianni, S. (2018). A bit of sugar on TiO_{2}: quantum chemical insights on the interfacial interaction of glycolaldehyde over titanium dioxide. *J. Phys. Chem. C* 122, 6041–6051. doi: 10.1021/acs.jpcc.7b11911

Tasinato, N., and Grimme, S. (2015). Unveiling the non-covalent interactions of molecular homodimers by dispersion-corrected DFT calculations and collision-induced broadening of ro-vibrational transitions: application to (CH_{2}F_{2})_{2} and (SO_{2})_{2}. *Phys. Chem. Chem. Phys.* 17, 5659–5669. doi: 10.1039/C4CP05680A

Tasinato, N., Moro, D., Stoppa, P., Pietropolli Charmet, A., Toninello, P., and Giorgianni, S. (2015). Adsorption of F_{2}CCFCl on TiO_{2} nano-powder: structures, energetics and vibrational properties from DRIFT spectroscopy and periodic quantum chemical calculations. *Appl. Surf. Sci.* 353, 986–994. doi: 10.1016/j.apsusc.2015.07.006

Tasinato, N., Pietropolli Charmet, A., Stoppa, P., Giorgianni, S., and Buffa, G. (2010). Spectroscopic measurements of SO_{2} line parameters in the 9.2 μm atmospheric region and theoretical determination of self-broadening coefficients. *J. Chem. Phys.* 132:044315. doi: 10.1063/1.3299274

Tasinato, N., Pietropolli Charmet, A., Stoppa, P., Giorgianni, S., and Gambi, A. (2012a). Quantum-chemical ab initio investigation of the vibrational spectrum of halon 1113 and its anharmonic force field: a joint experimental and computational approach. *Chem. Phys.* 397, 55–64. doi: 10.1016/j.chemphys.2011.12.015

Tasinato, N., Puzzarini, C., and Barone, V. (2017). Correct modeling of cisplatin: a paradigmatic case. *Angew. Chem. Int. Ed.* 56, 13838–13841. doi: 10.1002/anie.201707683

Tasinato, N., Regini, G., Stoppa, P., Pietropolli Charmet, A., and Gambi, A. (2012b). Anharmonic force field and vibrational dynamics of CH_{2}F_{2} up to 5000 cm^{−}1 studied by Fourier transform infrared spectroscopy and state-of-the-art ab initio calculations. *J. Chem. Phys.* 136:214302. doi: 10.1063/1.4720502

Watson, J. K. (2005). The vibrational dependence of quartic centrifugal distortion. *J. Mol. Struct.* 742, 91–98. doi: 10.1016/j.molstruc.2005.01.009

Westbook, B. R., Rio, W. A. D., Lee, T. J., and Fortenberry, R. C. (2020). Overcoming the out-of-plane bending issue in an aromatic hydrocarbon: the anharmonic vibrational frequencies of c-(CH)C_{3}${\text{H}}_{2}^{+}$. *Phys. Chem. Chem. Phys.* 22, 12951–12958. doi: 10.1039/D0CP01889A

Wolff, H., and Szydlowski, J. (1985). Vibrational spectra and rotational isomerism of ethanethiol and ethanethiol-d_{1}. *Can. J. Chem.* 63, 1708–1712. doi: 10.1139/v85-287

Woon, D. E., and Dunning, T. H. J. (1993). Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. *J. Chem. Phys.* 98, 1358–1371. doi: 10.1063/1.464303

Ye, H., Mendolicchio, M., Kruse, H., Puzzarini, C., Biczysko, M., and Barone, V. (2020). The challenging equilibrium structure of HSSH: another success of the rotational spectroscopy/quantum chemistry synergism. *J. Mol. Struct.* 1211:127933. doi: 10.1016/j.molstruc.2020.127933

Zerbetto, M., Licari, D., Barone, V., and Polimeno, A. (2013). Computational tools for the interpretation of electron spin resonance spectra in solution. *Mol. Phys.* 111, 2746–2756. doi: 10.1080/00268976.2013.800602

Keywords: quantum chemistry, density functional theory, rotational spectroscopy, vibrational spectroscopy, benchmark, atmospheric molecules, astrochemical molecules

Citation: 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:584203. doi: 10.3389/fchem.2020.584203

Received: 16 July 2020; Accepted: 17 August 2020;

Published: 23 October 2020.

Edited by:

Wei Hu, Lawrence Berkeley National Laboratory, United StatesReviewed by:

Igor Ying Zhang, Fudan University, ChinaHonghui Shang, Chinese Academy of Sciences, China

Copyright © 2020 Barone, Ceselin, Fusè and Tasinato. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Vincenzo Barone, vincenzo.barone@sns.it