Investigating the Calculation of Anharmonic Vibrational Frequencies

Apr 6, 2012 - *E-mail: [email protected]. .... James A. Calladine , Thomas S. Murphy , Michelle Hamilton , Ian P. Clark , Michael Towrie , ...
1 downloads 0 Views 926KB Size
Download PDF
Article pubs.acs.org/JPCA

Investigating the Calculation of Anharmonic Vibrational Frequencies Using Force Fields Derived from Density Functional Theory Magnus W. D. Hanson-Heine, Michael W. George, and Nicholas A. Besley* School of Chemistry, University of Nottingham, University Park, Nottingham NG7 2RD, U.K. S Supporting Information *

ABSTRACT: The calculation of anharmonic vibrational frequencies for a set of small molecules has been examined to explore the merit of applying such computationally expensive approaches for large molecules with density functional theory. The performance of different hybrid and gradient-corrected exchange-correlation functionals has been assessed for the calculation of anharmonic vibrational frequencies using second-order vibrational perturbation theory with two- and four-mode couplings and compared to the recently developed transition optimized shifted Hermite method. A range of exchange-correlation functionals (B3LYP, BLYP, EDF1, EDF2, B97-1, B97-2, HCTH-93, HCTH-120, HCTH-147, and HCTH-407) have been evaluated with reference to a large experimental data set comprising 88 species and 655 modes as well as a smaller set of shifts in frequency because of anharmonicity derived from experimental data. The anharmonic frequencies calculated using hybrid functionals provide the best agreement with experiment, and are not significantly improved by frequency scaling factors, indicating an absence of significant systematic error. For the molecules studied, the B97-1 and B97-2 functionals give the closest overall agreement with experiment, although the improvement over the best case for pure harmonic frequencies is modest. Predictions of the experimental anharmonic shifts are closest for the B3LYP and EDF2 functionals, with B97-1 performing well because of a good description of the harmonic force field. Investigations using modified hybrid functionals with increased fractions of Hartree−Fock exchange indicate that approximately 20% Hartree−Fock exchange is optimal.

■

INTRODUCTION With the continued development of both the modern computer and commercial quantum-chemical methods, the use of theoretical calculations has expanded for a wide range of applications, and the calculation of molecular vibrational properties has become an increasingly important spectroscopic tool.1 Such calculations are now performed routinely by both experimentalists and theoreticians, and are extremely useful in aiding the assignment and interpretation of experimental data. However, the vast majority of these calculations are based on the harmonic frequency approximation since to go beyond this approximation requires significantly greater sampling of the nuclear potential energy surface, and the calculations soon become prohibitive as the number of atoms increases. Of the electronic structure methods available, density functional theory (DFT) is potentially the most useful for the evaluation of nuclear force fields because of its relatively low computational cost making accurate calculations possible for larger systems, albeit with a wide choice of approximate exchange and correlation functionals.2 The use of DFT electronic structures in calculating anharmonic frequencies has therefore made a more rigorous treatment of molecular vibrations possible for molecules of continually increasing size,3 most commonly through second-order perturbation theory (VPT2).4 Infrared frequencies obtained from ab initio molecular dynamics © 2012 American Chemical Society

(AIMD) simulations using DFT have also been shown to include anharmonic effects by sampling regions of the potential energy surface accessed during classical trajectories. Recently AIMD frequencies have been reported for molecules as large as Ac-Ala15-LysH+ by Rossi et al.,5 Ala7H+ by Sediki et al.,6 and the sugars α-glucose, β-glucose, and sucrose by Gerber and coworkers,3f who have also recently employed vibrational selfconsistent field (VSCF) approaches to analyze the vibrational frequencies for a collection of large molecules with considerable success.7 In the harmonic approximation, the potential energy surface around an equilibrium structure is approximated by its secondorder Taylor series. Vibrational frequencies and normal modes can then be obtained simply from the diagonalization of the mass-weighted Hessian matrix, which requires only the second derivatives of the energy with respect to the nuclear coordinates to be evaluated. While this model provides an effective general treatment of vibrational structure, the specific anharmonicity in the fundamental frequencies is inherently ignored. Fortunately, the error in harmonic frequencies is quite systematic and results in the calculated frequencies being too high. This has led to the Received: February 20, 2012 Revised: April 2, 2012 Published: April 6, 2012 4417

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

Table 1. D3 and D4 Molecular Test Sets D3 Set: B2H6, C2Cl2, C2N2, cyclo-C2H4NH, cyclo-C2H4O, cyclo-C3H6, CH2CCH2, CH2CCHCl, CH2CCl2, CH2CH2, CH2CHCHO, CH2Cl2, 1CH2, 3CH2, CH3CH2Cl, CH3CH2CN, CH3CH2F, CH3CHO, CH3COF, CH3COOH, CH3NH2, CH, cis-CHClCHCl, cis-CHFCHF, Cl2O, ClCN, ClF3, ClNO2, ClNO, ClSN, CO2, COCl2, COClF, COF2, COS, CS2, CSF2, F2, F2NH, F2O, F2SO, FCN, FH, H2CO, H2, H2O2, H2O, H2S2, H2S, HCCCCH, HCCCH2Cl, HCCCH2F, HCCCl, HCCF, HCCH, HCN, HCO, HCOOCH3, HCOOH, HN3, HNCO, HNO3, HOCl, HOF, LiF, N2, N2F2, N2O, NCl2F, NClF2, NH, NO2, NSF, O2, O3, OCHCHO, OH, ONF, PH, S2F2, SCl2, SH, SiH2Cl2, SiH, SO2, SOCl2, trans-CHClCHCl, trans-CHFCHF. D4 Set: H2, CH, NH, OH, HF, N2, O2, F2, LiF, CO2, HCN, CH2Cl2, H2CO, HCCF, HCCH, H2O, H2S, CH2CH2.

literature.12b,15 However, little exists by way of extensive comparison with experiment for a large set of molecules with a more varied electronic structure, and the newly developed TOSH method has not yet been thoroughly evaluated. In this paper, a range of widely used hybrid and gradient-corrected exchange-correlation functionals for the calculation of anharmonic vibrational frequencies are examined. Harmonic, 4MR VPT2, 2MR VPT2, and TOSH vibrational frequencies for the B3LYP,16 BLYP,17 EDF1,18 EDF2,19 B97-1,20 B97-2,21 HCTH93,20 HCTH-120,22 HCTH-147,22 and HCTH-40723 functionals with the 6-311+G(d,p) basis set, are compared to experiment for a large set of molecules and vibrational frequencies. The error in the anharmonic frequencies is also examined through comparison with a set of experimental anharmonic shifts, obtained by the difference between experimental frequencies and experimental projections toward infinitesimal amplitudes of vibration. The use of empirical linear frequency scaling as a further correction for these calculated anharmonic frequencies is examined, and the effect on the computed anharmonic frequencies resulting from varying fractions of Hartree−Fock exchange in the exchangecorrelation functional is investigated. The choice of an appropriate basis set that gives both sufficient accuracy and speed for application to larger systems is problematic since the calculations soon become intractable for very large basis sets. A study by Boese et al.15a indicates that polarized triple-ζ basis sets can be sufficiently complete for accurate DFT anharmonic vibrational frequency analysis, and separate works by Andersson and Uvdal,8a and Merrick et al.,8b suggest a general basis set convergence for scaled harmonic DFT frequencies by the 6-311+G(d,p) and 6-311+G(d) levels. For this study we have chosen the medium sized 6-311+G(d,p) basis set, and although basis set calibration studies by Jensen have shown an error of ∼6/7 cm−1 in harmonic frequencies calculated at this level,24 this error is expected to remain constant throughout.

incorporation of anharmonicity through empirical scaling factors, which has proved to be an effective correction, with a number of authors reporting scaling factors at various levels of theory.8 The main drawbacks of this approach are that the same scaling factor will not be appropriate for all vibrational modes. Although multiple scaling factors have been used to try and address this problem,8h this increased the level of empiricism and can also limit transferability of this approach. The mode specific nature of anharmonicity is important for band assignment where reordering of the vibrational modes occurs and where anharmonicity does not follow general molecular trends.3c,d Many comparisons between calculations and experimental data are based on a best fit as indicated by matching major bands and the overall spectral profile. The calculation of anharmonic frequencies may reveal subtle changes such as reordering of weaker vibrational modes. Insight into calculated anharmonic vibrational spectra also has the potential to provide a deeper understanding in a wide range of experimental techniques, including vibrational circular dichroism (VCD),9 ultrafast 2-dimensional infrared spectroscopy (2D-IR),10 and the assignment of infrared spectra for lowest triplet electronically excited states. For example, anharmonic effects are expected to be significantly more pronounced in the latter.11 Consequently, it remains desirable to have accurate calculations of vibrational frequencies incorporating anharmonicity from first principles to ensure their continued widespread applicability with increasing molecular size. Recently, hybrid schemes that combine accurate methods such as coupled-cluster theory for the harmonic force field with DFT for the anharmonic part have been demonstrated to give close agreement with experiment.12 Although the use of methods such as coupled cluster theory for ever larger systems in the future is limited by the rapid increase in computational cost compared to conventional DFT, the use of accurate harmonic force fields is clearly important.12 When calculating anharmonic frequencies with VPT2, the computational cost is reduced by using a quartic force field (QFF), and the coupling between modes reduced to an n-mode representation (nMR). However, VPT2 is also known to suffer from numerical failure in cases of normal mode degeneracy, which can be reduced by artificially altering the nuclear masses.13 An alternative approach to calculating anharmonic corrections called the transition optimized shifted Hermite (TOSH) method has been recently introduced,14 and is designed to correct for VPT2 failures in cases of spectral resonance and degeneracy by incorporating anharmonic information directly into the unperturbed wave function. This approach provides a slight reduction in anharmonic computational cost and subsequently has an expanded potential applicability within some truncations of the 2MR force field.14 The increased availability of anharmonic methods coupled with increased computational resources has made it possible to compare accurate computed vibrational frequencies from different electronic structure theories directly with experiment for a range of molecules, and a few studies have appeared in the

■

COMPUTATIONAL DETAILS

Geometry optimization and frequency calculations were performed using the Q-Chem quantum-chemical software suite25 using the B3LYP, BLYP, EDF1, EDF2, B97-1, B97-2, HCTH-93, HCTH-120, HCTH-147, and HCTH-407 functionals, a triple-ζ polarized 6-311+G(d,p) basis set, and the default SG-1 numerical integration grid,26 which was chosen for its applicability to larger systems. Prior to frequency evaluation, initial molecular structures were optimized to minimum energy geometries at each level of theory using tightened energy convergence criteria of 1 × 10−7 Eh, and a gradient convergence of 1 × 10−6 Eh a0−1 to ensure the efficacy of the harmonic and postharmonic approaches, and the resulting DFT force fields were used for respective harmonic, TOSH, and VPT2 2MR and 4MR vibrational calculations. Harmonic and 2MR TOSH calculations were also performed for hybrid modifications of the BLYP functional using the same procedure. 4418

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

we will briefly discuss the calculated harmonic frequencies since they provide a useful benchmark to assess the value of the anharmonic calculations. For gradient-corrected forms of the exchange-correlation functional the unscaled calculated frequencies are in surprisingly close agreement with the observed frequencies of the main test set, and show a level of accuracy similar to many of the anharmonic theoretical methods discussed later. Without the use of empirical scaling factors, this agreement is closest for the EDF1 functional, as has been noted previously for EDF1 in the case of amide dimers by Watson and Hirst.30 This agreement between harmonic frequencies and experiment, 36 cm−1 mean absolute deviation (MAD) for EDF1, is clearly derived from a fortuitous cancellation of errors arising from deficiencies in basis set, the form of the exchange-correlation functional, and the neglect of anharmonicity, which is particularly apparent for the gradientcorrected functionals. The EDF2 hybrid functional is an empirical reparameterization of the EDF1 and B3LYP functionals by Lin et al.,19 using a training set consisting purely of harmonic experimental frequencies mainly from diatomic species, with the intention of generating an accurate potential energy surface for second-derivatives of the energy with respect to nuclear position. It is therefore important to note that the harmonic frequencies obtained by this functional are not expected to agree well with the anharmonic experimental data of this test set, though they may provide a clearer indication of the harmonic approximations true performance with the second worst MAD from experiment of 49 cm−1. In line with previous studies, the harmonic frequencies of the D3 set are systematically overestimated using hybrid functionals, with between 25 and 38% of all frequencies showing a relative error inside the 98% confidence interval with respect to experimental values, as shown in Table 2. It is well-known that

The errors for common molecular benchmark sets have been evaluated, with errors determined using a set comprising 88 molecules across 655 frequencies, excluding atoms beyond the second periodic row, and having no more than 4 non-hydrogen atoms per molecule, with no more than 10 atoms in total per molecule. The set was derived from work by Shimanouchi27 and Johnson et al.28 and included only the most common isotopomer of each molecule in line with work by Lin et al.,14 ultimately incorporating the molecules shown in Table 1 which we have designated the D3 set. Errors in calculated anharmonic shifts were also evaluated using a subset of molecules, termed D4, for which experimental harmonic frequencies are available. Experimental harmonic frequencies were taken from the molecule specific references cited by Lin et al. in reference 19.28,29 Experimental anharmonic shifts were calculated from the difference between the harmonic and the anharmonic experimental frequencies of the D4 set, while calculated anharmonic shifts were taken to be the difference between harmonic and anharmonic DFT frequencies. Mean absolute experimental errors were also calculated for the full set of vibrational modes and are reported with the exclusion of all modes with absolute errors in excess of 200 cm−1, which were taken to be natural outliers arising from unphysical artifacts in the theoretical treatment, as well as removing spectral degeneracies from VPT2. Harmonic vibrational frequencies and normal modes were determined using analytic second derivatives of the energy with respect to nuclear displacement in all of the functionals studied with the exception of the HCTH family of functionals for which analytical second derivatives were not available. For the HCTH-93, HCTH-120, HCTH-147, and HCTH-407 functionals, the second derivatives were evaluated numerically from finite differences of the analytic gradients using the default step size of 0.189 a0. Third and fourth derivatives of the energy with respect to nuclear displacement used in the computation of anharmonic frequencies were also derived from finite differences of the analytic energy, gradient, and Hessian, where available, using the default step size of 0.1 a0,14 and using gradients only for the HCTH functionals. Linear frequency scaling factors for both harmonic and anharmonic frequencies were also considered, and were calculated following the methodology of least-squares fitting minimization, where scale factors at each level of theory are determined by minimization of the residual Δ=

∑ (λωi − νi)2

Table 2. Percentage of Theoretical Harmonic Frequencies That Lie within Specified Relative Error Ranges When Compared with Experimental Values

(1)

leading to the expression λ=

∑ ωiνi/∑ ωi 2

(2)

where the summations are over the N vibrational modes included in the optimization, ωi and νi are the wavenumber values of the ith theoretical and ith experimental fundamental frequencies, respectively, and λ is the scaling factor.8g Summary tables and figures giving overall set values are presented in this paper, and a complete listing of the molecule by molecule uncorrected mean unsigned frequency errors derived at each level of theory is available as Supporting Information, Tables S1−S25.

% error

B3LYP

EDF2

B97-1

B97-2

HCTH-407

0−2 2−4 4−6 6−8 8−10 10−12 12−14 14−16 16−18 18−20 20−22 22−24 24−26 26−28 28−30 >30

34.7 34.0 16.3 7.2 2.7 1.7 1.1 0.5 0.6 0.2 0.2 0.3 0.5 0.0 0.0 0.2

35.7 29.9 19.7 6.7 3.8 0.5 0.8 1.1 0.6 0.2 0.2 0.3 0.3 0.0 0.0 0.3

37.7 34.2 15.7 6.1 2.4 0.8 0.8 0.9 0.2 0.3 0.2 0.2 0.2 0.0 0.0 0.5

25.0 27.0 30.7 7.6 4.3 2.6 0.5 0.5 0.3 0.5 0.3 0.2 0.3 0.0 0.0 0.3

41.8 28.7 13.1 5.5 3.1 2.4 1.4 1.4 0.3 0.2 1.1 0.2 0.0 0.0 0.0 0.9

these frequencies are significantly improved by uniform scaling, and here the MAD is reduced to 26 cm−1 for B3LYP and B971, surpassing the accuracy of the EDF1 functional. In this way, scaling factors can be used to effectively reduce both the extent and the spread of errors, and although empirical scaling fails to account for the mode specific anharmonicity that can prove crucial for accurate vibrational assignment, matching its global

■

RESULTS AND DISCUSSION The main focus of this paper is the calculation of anharmonic IR frequencies for the D3 and D4 test sets listed above. Initially 4419

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

accuracy for a wide range of molecules represents a continuing challenge for first principles anharmonic methods. The incorporation of anharmonicity directly through the 4MR VPT2 method leads to a significant improvement in accuracy for the hybrid functionals, with MADs decreasing by 15−27 cm−1. A total of 56−57% of frequencies are within the 98% confidence interval for relative errors, and a more symmetrical distribution of signed errors is observed, as shown in Table 3 and more generally for all the methods in Table 3. Percentage of Theoretical 4MR VPT2 Frequencies That Lie within Specified Relative Error Ranges When Compared with Experimental Values % error

B3LYP

EDF2

B97-1

B97-2

HCTH-407

0−2 2−4 4−6 6−8 8−10 10−12 12−14 14−16 16−18 18−20 20−22 22−24 24−26 26−28 28−30 >30

55.7 20.5 6.9 5.5 4.6 2.6 1.4 0.2 0.6 0.2 0.3 0.2 0.0 0.0 0.2 1.4

57.4 19.4 7.2 5.8 3.7 2.1 1.5 0.5 0.2 0.2 0.3 0.3 0.0 0.0 0.0 1.5

56.5 19.4 7.6 5.3 3.7 1.5 2.1 0.5 0.8 0.2 0.2 0.0 0.0 0.2 0.0 2.1

57.1 20.0 8.2 4.7 4.1 1.4 0.8 0.5 0.3 0.8 0.2 0.0 0.0 0.2 0.0 1.8

40.3 24.7 11.5 7.5 3.1 3.8 2.0 1.8 1.2 0.3 0.8 0.2 0.5 0.2 0.2 2.1

Figures 1−4. For the gradient-corrected functionals, the previously accurate BLYP and EDF1 functionals now underestimate experimental frequencies, with MAD values increasing by 33 cm−1 and 12 cm−1 respectively, suggesting that the previous cancellation of errors is less significant. There is a smaller yet still significant net deterioration for HCTH-93, HCTH-120, and HCTH-147 and a slight improvement is observed for the HCTH-407 functional. Overall, B97-2 performs best with the 4MR VPT2 method showing a MAD of 28 cm−1, closely followed by B97-1 and EDF2, then B3LYP. The significant improvement in performance seen for all four hybrids in comparison to the gradientcorrected methods, together with the relatively small differences in performance between hybrid functionals further illustrates the importance of exact exchange for the accurate description of vibrational properties. The hybrid functionals B3LYP, EDF2, B97-1, and B97-2 contain 20%, 17%, 21%, and 21% exact exchange, respectively, and it is interesting to note that the spectroscopically optimized EDF2 functional yields a lower fraction than the B3LYP functional from which it originated. The previously reported superiority of the B97-1 functional over other hybrids15a also diminishes significantly when applied to a wide range of functional groups, even though it still performs relatively well. Incorporating anharmonicity through the 2MR VPT2 method shows a similar ordering of functional performance to 4MR, with B97-1 marginally better than B97-2, and a similar distribution of signed errors across the functionals. The MAD for the 2MR force field is between 5 and 4 cm−1 greater than 4MR for the hybrid functionals and between 0 and 5 cm−1 lower for the gradient-corrected case. Relative errors in Table 4

Figure 1. Signed error distributions for calculated harmonic frequencies.

are also found to be similarly worse. While there is a small decrease in accuracy for the 2MR method, it does represent a significant decrease in computational time, and for the largest molecule studied in this work, CH3CH2CN, the 2MR calculation is approximately a factor of 7 times less expensive than the 4MR calculation. While low frequency vibrations are known to experience a higher degree of mode coupling than high frequency modes, separating the vibrational modes for the B97-1 functional into frequency ranges of low frequency (0− 1000 cm−1), mid range (1000−2000 cm−1), and high frequency (2000−4000 cm−1) showed an improved prediction with increased mode coupling of 5 cm−1 MAD for the low frequency modes, a deterioration of 1 cm−1 in the mid range, and a 9 cm−1 4420

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

Figure 2. Signed error distributions for calculated TOSH frequencies.

Figure 3. Signed error distributions for calculated 2MR VPT2 frequencies.

improvement for high frequency MAD where the overall anharmonicity is much larger. Incorporation of anharmonicity using the TOSH method with its current 2MR implementation14 gives a similar ordering of the functionals in terms of performance with similar error distributions shown in Table 5 as well as MAD values, however with the advantage of a natural treatment of the vibrational degeneracies in VPT2. Furthermore, despite its approximation of VPT2 theory, the TOSH method still performs better than the analogous 2MR VPT2 by ∼1 cm−1, even with degeneracies accounted for through error cancellation. VPT2 was estimated by Lin et al. to be at least 2-fold slower with the implementation of method specific force field reductions.14

The specific total MAD values of the different functionals for each vibrational method is illustrated by Figure 5, and the relative accuracy is summarized below. Harmonic Functionals EDF1 = BLYP > HCTH-93 = HCTH-120 = HCTH-147 = HCTH-407 > B97-1 > B3LYP > EDF2 > B97-2 Best MAD: 36 cm−1 VPT2 2MR Functionals B97-1 = B97-2 > EDF2 > B3LYP > HCTH-407 > HCTH-93 > HCTH-147 > HCTH-120 > EDF1 > BLYP Best MAD: 33 cm−1 TOSH Functionals 4421

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

Table 4. Percentage of Theoretical 2MR VPT2 Frequencies That Lie within Specified Relative Error Ranges When Compared with Experimental Values % error

B3LYP

EDF2

B97-1

B97-2

HCTH-407

0−2 2−4 4−6 6−8 8−10 10−12 12−14 14−16 16−18 18−20 20−22 22−24 24−26 26−28 28−30 >30

48.4 24.0 9.6 5.5 3.1 3.5 2.3 0.5 0.5 0.5 0.5 0.0 0.2 0.0 0.3 1.4

50.1 22.8 8.6 6.4 2.9 3.4 2.0 0.8 0.3 0.0 0.8 0.3 0.0 0.2 0.2 1.5

52.1 20.9 8.1 6.6 2.4 2.3 1.2 2.0 0.8 0.6 0.5 0.2 0.2 0.0 0.2 2.1

48.2 25.3 9.6 6.7 3.5 1.5 1.2 0.5 0.2 0.5 0.5 0.2 0.0 0.3 0.2 1.7

41.4 26.4 11.2 6.7 3.1 2.8 1.7 1.7 0.9 0.6 0.5 0.0 0.2 0.2 0.6 2.3

Table 5. Percentage of Theoretical TOSH Frequencies That Lie within Specified Relative Error Ranges When Compared with Experimental Values % error

B3LYP

EDF2

B97-1

B97-2

HCTH-407

0−2 2−4 4−6 6−8 8−10 10−12 12−14 14−16 16−18 18−20 20−22 22−24 24−26 26−28 28−30 >30

48.7 24.3 9.3 5.7 2.9 3.4 2.3 0.5 0.3 0.5 0.5 0.0 0.2 0.0 0.3 1.4

50.1 22.8 8.9 6.7 2.4 3.4 2.0 0.8 0.3 0.0 0.8 0.3 0.0 0.2 0.2 1.4

53.7 20.0 7.9 5.8 3.1 2.3 1.1 2.0 0.6 0.6 0.5 0.2 0.2 0.0 0.2 2.0

47.6 26.4 9.8 7.3 3.1 1.4 1.2 0.3 0.3 0.5 0.3 0.2 0.0 0.2 0.2 1.4

40.6 27.9 10.1 6.4 3.5 2.6 1.8 2.1 0.9 0.6 0.6 0.0 0.2 0.2 0.5 2.0

error, and suggesting that more accurate electronic structure theories are required for significant improvements in accuracy. The theoretical anharmonic corrections were further tested against a smaller set of experimental anharmonic corrections, and the calculated errors are reported in Table 7. Although the B97-1 functional performed better in reproducing observable frequencies, anharmonic shifts were best represented by B3LYP, giving a frequency shift MAD of 26, 24, and 19 cm−1 for the TOSH, 2MR VPT2, and 4MR VPT2 methods, respectively. Calculated anharmonic frequency shifts are in part dependent on the computation of accurate harmonic frequencies and hence second derivatives of the potential energy surface, which are worse for the B3LYP functional than B97-1. This importance of the harmonic part of the force field provides the basis for hybrid methods that combine high level second derivatives from coupled cluster theory with third and fourth derivatives from DFT, which achieve good agreement with experiment. Further uncertainty in the experimental values of these shifts owing to their being gathered from a variety of

Figure 4. Signed error distributions for calculated 4MR VPT2 frequencies.

B97-1 > B97-2 > EDF2 = B3LYP > HCTH-407 > HCTH-93 > HCTH-147 > HCTH-120 > EDF1 > BLYP Best MAD: 32 cm−1 VPT2 4MR Functionals B97-2 > B97-1 = EDF2 > B3LYP > HCTH-407 > HCTH-93 > HCTH-147 > HCTH-120 > EDF1 > BLYP Best MAD: 28 cm−1 Table 6 shows empirical frequency scaling factors derived for the 2MR TOSH and 4MR VPT2 anharmonic methods. For the gradient-corrected functionals BLYP and EDF1 this provides an empirical correction for the frequency underestimation. For both VPT2 and TOSH, the hybrid functionals have scaling factors close to unity, indicating that there is no systematic 4422

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

Figure 5. Mean absolute errors (cm−1) between experimental and calculated frequencies for the different DFT functionals, using the harmonic, 2MR, and 4MR VPT2 methods (columns from left to right) shown for each functional excluding modes with error in excess of 200 cm−1.

Table 6. Optimal Frequency Scaling Factorsa Harm. B3LYP EDF2 B97-1 B97-2 BLYP EDF1 HCTH-93 HCTH-120 HCTH-147 HCTH-407

0.9665 0.9642 0.9669 0.9555 0.9993 0.9844 0.9788 0.9801 0.9791 0.9747

(−22) (−23) (−20) (−34) (1) (−3) (−6) (−5) (−5) (−8)

TOSH 1.0066 1.0042 1.0000 0.9947 1.0436 1.0261 1.0195 1.0210 1.0201 0.9992

B88 exchange functional with HF exchange in increments of 10% between BLYP (0%) and HFLYP (100%). The results of this modification are shown in Table 8 and graphically in Figure

VPT2

(0) (0) (0) (−2) (−28) (−12) (−7) (−9) (−8) (0)

1.0039 1.0032 0.9954 0.9939 1.0418 1.0268 1.0119 1.0193 1.0180 0.9925

(−1) (−1) (1) (0) (−35) (−18) (−7) (−11) (−10) (6)

Table 8. Observed MAD Values for Hybrid Modifications of the BLYP Functional in cm−1

Derived from a least-squares fitting procedure using modes calculated for the specific molecular set and vibrational approximations. The changes in mean absolute error in cm−1 arising from scaling are quoted in parentheses. a

Table 7. Experimental Anharmonic Frequency Shift Errors in cm−1 B3LYP EDF2 B97-1 B97-2 BLYP EDF1 HCTH-93 HCTH-120 HCTH-147 HCTH-407

Harm.a

TOSHb

2MRb

4MRb

26 27 23 37 60 36 30 31 30 29

26 26 31 27 31 31 29 30 31 34

24 24 31 26 31 30 28 29 30 32

19 20 27 29 26 46 59 36 38 35

HF (%)

Harm.

TOSH

scaled

scale factor

0 10 20 30 40 50 60 70 80 90 100

37 36 47 63 83 103 123 142 161 179 198

72 48 36 38 48 67 86 107 127 145 164

37 30 26 25 25 27 29 32 36 39 45

0.9992 0.9833 0.9685 0.9546 0.9417 0.9296 0.9183 0.9076 0.8975 0.8879 0.8788

6. Compared with observed frequencies, the accuracy of the harmonic calculations deteriorates in a linear fashion as the faction of exact exchange increases beyond ∼10%. TOSH errors are minimized with the inclusion of ∼20% exact exchange, increasing linearly above 30−40%. Errors in the scaled harmonic frequencies were also minimized by the inclusion of 30−40% HF exchange. This result for scaled frequencies is in line with the scaling factor work of Merrick et al.8b using a modification of the B3LYP functional. A key consideration in any calculation is the computational cost. As stated above, the cost in performing anharmonic calculations is much greater than for the corresponding harmonic predictions. To provide indicative timings for the different methods, we have performed specific calculations for CH3CH2CN on a standalone Intel Xeon E5506 2.13 GHz CPU. As expected, the order of computational cost is harmonic (≪ 1 h) < 2MR VPT2 (∼3 h) < 4MR VPT2 (∼21 h). This illustrates the dramatic increase in cost of 4MR compared to 2MR, and also the increased cost of anharmonic calculations compared to harmonic. In the current Q-Chem implementation, the 2MR VPT2 and TOSH methods have equal costs. Comparison of functionals for CH3CH2CN show that hybrid functionals are about 30% more expensive than the pure gradient-corrected functionals, taking ∼2.8 h for B3LYP and ∼2.2 h for BLYP in the 2MR implementation.

a

The harmonic MAD values are the differences between the experimental and theoretical harmonic frequencies in cm−1. bThe anharmonic MAD values are the differences between the experimental and theoretical anharmonic shifts in cm−1.

sources and extrapolated from experimental spectra, as well as the reduced size and variation of the D4 subset also limits wider conclusion. However, the results indicate that the superior performance of the B97-1 functional in reproducing observable vibrations is largely due to accuracy in the harmonic component, which showed an improvement over both B3LYP and EDF2, in spite of the EDF2 parametrization. The anharmonic shifts of the gradient-corrected functionals were of lower quality, and all, with the exception of BLYP, become worse moving from 2MR to 4MR mode-coupling. The effects of increasing the percentage of exact Hartree− Fock exchange were also examined for the D3 set. A modified BLYP functional was constructed by replacing a portion of the

■

CONCLUSIONS We have reported an extensive assessment of hybrid and gradient-corrected exchange-correlation functionals for a large 4423

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

Figure 6. Mean absolute errors between experimental and calculated frequencies as a function of the percentage of exact exchange mixed with the BLYP functional, calculated using (a) the harmonic approximation and (b) the 2MR TOSH approximation.

■

ACKNOWLEDGMENTS We thank The University of Nottingham for funding. M.W.G. gratefully acknowledges receipt of a Royal Society Wolfson Merit Award.

set of molecules comprising many different functional groups. This work both supplements and consolidates previous studies into the performance of DFT force fields for anharmonic frequency calculations that have employed much smaller test sets of vibrational modes. Furthermore, this work benchmarks the newer TOSH method, which shows a marginal improvement in both accuracy and speed for this test set. Our results emphasize the importance of ∼20% exact exchange for the calculation of anharmonic frequencies, and show that the B97-1 and B97-2 exchange-correlation functionals perform best out of those tested within a standard triple-ζ basis set. However, all the hybrid functionals tested perform to a similar level, and while each functional's anharmonic correction is significant, their improvements over the best harmonic values are less pronounced. Accurate description of the harmonic potential energy surface also still seems to be a dominant factor in correctly representing observable vibrations even within the hybrid schemes tested. For the gradient-corrected functionals, the incorporation of anharmonicity also does not result in an improved agreement with experiment, although they outperform hybrid functionals when their harmonic frequencies are compared directly with observables. The results for the 2MR VPT2 and TOSH methods show only a slightly larger deviation from experiment than 4MR VPT2, and the lower computational cost and implicit inclusion of vibrational degeneracies make the 2MR TOSH method particularly suitable for the study of larger systems where time is often a limiting factor and degeneracies are more likely to occur. Despite these findings, the failure of VPT2 coupled with DFT to provide a significant improvement over the best case of harmonic frequencies with a wider test set of vibrational modes indicates that there is still an unmet demand for an accurate and theoretically rigorous treatment of anharmonicity that can be applied to larger scale chemistry.

■

■

(1) (a) Koch, W.; Holthausen, C. A Chemist’s Guide to Density Functional Theory; Wiley-VCH: Weinheim, Germany, 2000; (b) Meier, R. J. Vib. Spectrosc. 2007, 43, 26−37. (c) Jonas, V.; Thiel, W. J. Chem. Phys. 1996, 105, 3636−3648. (d) Besley, N. A.; Bryan, J. A. J. Phys. Chem. C 2008, 112, 4308−4314. (e) Gaigeot, M.-P. Phys. Chem. Chem. Phys. 2010, 12, 3336−3359. (f) Langhoff, S. R.; Bauschlicher, C. W.; Hudgins, D. M.; Sandford, S. A.; Allamandola, L. J. J. Phys. Chem. A 1998, 102, 1632−1646. (2) Kohn, W.; Becke, A. D.; Parr, R. G. J. Phys. Chem. 1996, 100, 12974−12980. (3) (a) Barone, V. J. Phys. Chem. A 2004, 108, 4146−4150. (b) Boese, A. D.; Martin, J. M. L. J. Phys. Chem. A 2004, 108, 3085− 3096. (c) Cane, E.; Trombetti, A. Phys. Chem. Chem. Phys. 2009, 11, 2428−2432. (d) Cane, E.; Miani, A.; Trombetti, A. J. Phys. Chem. A 2007, 111, 8218−8222. (e) Miani, A.; Cane, E.; Palmieri, P.; Trombetti, A.; Handy, N. C. J. Chem. Phys. 2000, 112, 248−259. (f) Brauer, B.; Pincu, M.; Buch, V.; Bar, I.; Simons, J. P.; Gerber, R. B. J. Phys. Chem. A 2011, 115, 5859−5872. (g) Dreyer, J. J. Chem. Phys. 2007, 127, 054309. (h) Gaigeot, M.-P.; Besley, N. A.; Hirst, J. D. J. Phys. Chem. B 2011, 115, 5526−5535. (4) Lounila, J.; Wasser, R.; Diehl, P. Mol. Phys. 1987, 62, 19−31. (5) Rossi, M.; Blum, V.; Kupser, P.; von Helden, G.; Bierau, F.; Pagel, K.; Meijer, G.; Scheffler, M. J. Phys. Chem. Lett. 2010, 1, 3465−3470. (6) Sediki, A.; Snoek, L. C.; Gaigeot, M. P. Int. J. Mass Spectrom. 2011, 308, 281−288. (7) (a) Schweke, D.; Brauer, B.; Gerber, R. B.; Haas, Y. Chem. Phys. 2007, 333, 168−178. (b) Pele, L.; Gerber, R. B. J. Phys. Chem. C 2010, 114, 20603−20608. (c) Xie, H.-b.; Pincu, M.; Brauer, B.; Gerber, R. B.; Bar, I. Chem. Phys. Lett. 2011, 514, 284−290. (8) (a) Andersson, M. P.; Uvdal, P. J. Phys. Chem. A 2005, 109, 2937−2941. (b) Merrick, J. P.; Moran, D.; Radom, L. J. Phys. Chem. A 2007, 111, 11683−11700. (c) Tantirungrotechai, Y.; Phanasant, K.; Roddecha, S.; Surawatanawong, P.; Sutthikhum, V.; Limtrakul, J. J. Mol. Struct. (Theochem) 2006, 760, 189−192. (d) Wong, M. W. Chem. Phys. Lett. 1996, 256, 391−399. (e) Defrees, D. J.; Mclean, A. D. J. Chem. Phys. 1985, 82, 333−341. (f) Pople, J. A.; Scott, A. P.; Wong, M. W.; Radom, L. Isr. J. Chem. 1993, 33, 345−350. (g) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502−16513. (h) Rauhut, G.; Pulay, P. J. Phys. Chem. 1995, 99, 3093−3100. (9) Nafie, L. A. Theory of Vibrational Circular Dichroism. In Vibrational Optical Activity: Principles and Applications; John Wiley & Sons, Ltd: Chichester, U.K., 2011. (10) Hamm, P.; Zanni, M. T. Concepts and Methods of 2D Infrared Spectroscopy; Cambridge Press: Cambridge, U.K., 2011. (11) Toscano, J. P. Structure and Reactivity of Organic Intermediates as Revealed by Time-Resolved Infrared Spectroscopy. In Advances in Photochemistry; Neckers, D. C., Bünau, G. V., Jenks, W. S., Eds.; John Wiley & Sons, Inc.: Hoboken, NJ, 2007; Vol. 26.

ASSOCIATED CONTENT

S Supporting Information *

Tables listing molecule specific MAD values associated with the prediction of vibrations at each level of electronic structure theory (Tables S1−S25). This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 4424

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425

The Journal of Physical Chemistry A

Article

(12) (a) Puzzarini, C.; Biczysko, M.; Barone, V. J. Chem. Theory Comput. 2011, 7, 3702−3710. (b) Biczysko, M.; Panek, P.; Scalmani, G.; Bloino, J.; Barone, V. J. Chem. Theory Comput. 2010, 6, 2115− 2125. (13) Barone, V. J. Chem. Phys. 2005, 122, 10. (14) Lin, C. Y.; Gilbert, A. T. B.; Gill, P. M. W. Theor. Chem. Acc. 2008, 120, 23−35. (15) (a) Boese, A. D.; Klopper, W.; Martin, J. M. L. Int. J. Quantum Chem. 2005, 104, 830−845. (b) Carbonniere, P.; Barone, V. Chem. Phys. Lett. 2004, 399, 226−229. (c) Carbonniere, P.; Lucca, T.; Pouchan, C.; Rega, N.; Barone, V. J. Comput. Chem. 2005, 26, 384− 388. (d) Rauhut, G.; Knizia, G.; Werner, H. J. J. Chem. Phys. 2009, 130. (16) (a) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (b) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (17) (a) Becke, A. D. Phys. Rev. A 1988, 38, 3098−3100. (b) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37, 785−789. (18) Adamson, R. D.; Gill, P. M. W.; Pople, J. A. Chem. Phys. Lett. 1998, 284, 6−11. (19) Lin, C. Y.; George, M. W.; Gill, P. M. W. Aust. J. Chem. 2004, 57, 365−370. (20) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 1998, 109, 6264−6271. (21) Wilson, P. J.; Bradley, T. J.; Tozer, D. J. J. Chem. Phys. 2001, 115, 9233−9242. (22) Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M. J. Chem. Phys. 2000, 112, 1670−1678. (23) Boese, A. D.; Handy, N. C. J. Chem. Phys. 2001, 114, 5497− 5503. (24) Jensen, F. J. Chem. Phys. 2003, 118, 2459−2463. (25) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (26) Gill, P. M. W.; Johnson, B. G.; Pople, J. A. Chem. Phys. Lett. 1993, 209, 506−512. (27) (a) Shimanouchi, T. Tables of Molecular Vibrational Frequencies Consolidated Vol. I. National Bureau of Standards: Washington, DC, 1972; p 164; (b) Shimanouchi, T. J. Phys. Chem. Ref. Data 1977, 6, 993−1102. (28) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1993, 98, 5612−5626. (29) (a) Harris, N. J. J. Phys. Chem. 1995, 99, 14689−14699. (b) Thomas, J. R.; Deleeuw, B. J.; Vacek, G.; Crawford, T. D.; Yamaguchi, Y.; Schaefer, H. F. J. Chem. Phys. 1993, 99, 403−416. (c) Miller, K. J.; Gandakesuma, F. S. J. Mol. Spectrosc. 1991, 145, 429− 447. (30) Watson, T. M.; Hirst, J. D. J. Phys. Chem. A 2002, 106, 7858− 7867.

4425

dx.doi.org/10.1021/jp301670f | J. Phys. Chem. A 2012, 116, 4417−4425