Use of Low-Cost Quantum Chemistry Procedures for Geometry

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX-XXX

pubs.acs.org/JCTC

Use of Low-Cost Quantum Chemistry Procedures for Geometry Optimization and Vibrational Frequency Calculations: Determination of Frequency Scale Factors and Application to Reactions of Large Systems Bun Chan* Graduate School of Engineering, Nagasaki University, Bunkyo 1-14, Nagasaki 852-8521, Japan S Supporting Information *

ABSTRACT: We have assessed the performance of a variety of low-cost computational quantum chemistry procedures (semiempirical, pure-DFT, and screened-exchange DFT methods) for computing molecular geometries and thermochemical quantities associated with the vibrational frequencies. Frequency scale factors for zero-point vibrational energies and thermal corrections for 298 K enthalpies and 298 K entropies have been determined. In absolute terms, for small to mediumsized molecules, all procedures perform reasonably well. Semiempirical methods have mean absolute deviations (MADs) of ∼15 kJ mol−1 for total enthalpies and free energies. For DFT procedures, hybrid DFT generally performs better than pure DFT. Remarkably, the N12 pure functional shows very good performances (MADs ∼ 3 kJ mol−1) that are comparable to those for hybrid functionals. An examination of the basis set effect indicates N12/3-21G* and N12/6-31G(d) to be cost-effective for geometry optimization and vibrational frequency calculations, but the use of minimal basis sets leads to very large MADs for the calculated thermochemical quantities. Further testing with reaction energies of large systems shows that, by exploiting cancellation of systematic deviations, although the deviations can be very substantial in absolute terms (>100 kJ mol−1), those for relative energies are markedly reduced (∼10 kJ mol−1). This enables the use of even semiempirical procedures to obtain geometries and vibrational frequencies with reasonable accuracy in cases where the use of more expensive procedures is computationally prohibitive.

■

INTRODUCTION The calculation of thermochemical properties such as reaction energies and barriers is one of the most important utilities of computational quantum chemistry. To obtain such relative energies typically involves several steps. For each species involved, its geometry is optimized using a relatively economical method. This is followed by the calculation of vibrational frequencies at the same level to obtain quantities such as zero-point vibrational energies (ZPVEs) and thermal corrections to enthalpies (ΔHT). A refined single-point energy is usually computed with a higher-level procedure in order to improve the accuracy in the calculated thermochemical quantity. Among the various processes involved, obtaining accurate single-point energies is arguably the most difficult task. Much research effort has been dedicated to tackle this challenge with remarkable achievements.1,2 At one end of the scale, it has become routine nowadays to obtain thermochemical quantities with chemical accuracy for small to medium-sized molecules using high-level methodologies such as the Gn3 and Wn1,2 protocols. At the other end, the emergence of advanced density © XXXX American Chemical Society

functional theory (DFT) methods has enabled the calculations of sizable systems with adequate accuracy.4−6 Regardless of the method used for single-point-energy calculations, adequate geometries and quantities related to vibrational frequencies are also essential for overall reliability. For the treatment of small to medium-sized molecules, the use of hybrid-DFT or double-hybrid (DH) DFT geometries provides a good balance between accuracy and computational efficiency.7,8 Associated vibrational frequencies are typically computed within the harmonic oscillator approximation due to the high cost for explicit consideration of anharmonicity. Empirically determined scale factors are often applied to the harmonic vibrational frequencies in order to minimize systematic deviations related to anharmonicity and deficiencies in the method.8−14 The use of reliable hybrid-DFT or DH-DFT geometries and scaled vibrational frequencies is especially important for the use in combination with high-level singlepoint energies.15,16 This is because shortcomings associated Received: July 7, 2017

A

DOI: 10.1021/acs.jctc.7b00721 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 1. Scale Factors for Fundamentals, Low Frequencies (≤1000 cm−1), Zero-Point Vibrational Energies (ZPVE), and Thermal Corrections for 298 K Enthalpies (ΔH298) and Entropies (S298) for a Selection of Methods Examined AM1 PM3 PM6 B-LYP B-PW91 PBE B-LYP-D3BJ PBE-D3BJ N12 HSE HISS N12-SX B3-LYP B3-PW91

6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) VTZ+d VTZ+d

freq

low freq

ZPVE

ΔH298

S298

0.9542 0.9761 1.0611 0.9940 0.9865 0.9874 0.9937 0.9873 0.9625 0.9521 0.9262 0.9436 0.9681 0.9650

1.6921 2.4334 1.5976 1.0635 1.0483 1.0454 1.0584 1.0441 1.1654 0.9828 0.9484 0.9770 1.0028 0.9996

0.9772 0.9859 1.0851 1.0023 1.0083 1.0011 1.0024 1.0012 0.9849 0.9651 0.9391 0.9652 0.9886 0.9860

0.9851 1.0651 1.0654 1.0650 1.0459 1.0440 1.0632 1.0435 1.0205 0.9791 0.9390 0.9709 0.9926 0.9833

1.0615 1.1389 1.1719 1.0687 1.0506 1.0483 1.0662 1.0478 1.0351 0.9825 0.9438 0.9745 0.9970 0.9894

screened-exchange DFT procedures, namely HSE,34 N12-SX,5 and HISS.35 We include B3-LYP as it is fairly ubiquitous for molecular geometry optimization and vibrational frequency calculations. It is the prescribed method for these purposes within many widely used high-level composite protocols. 1,2 The B3-PW91 procedure is included as a closely related alternative to B3LYP due to some known shortcomings of the LYP correlation functional in some applications (see, for example ref 36). Screened-exchange methods can potentially offer the optimal compromise between accuracy and computational efficiency due to the exclusion of computationally intensive Hartree− Fock calculations in the long-range. For the basis sets employed in the present study, because we aim to search for computationally efficient methodologies, we focus mainly on minimal and double-ζ basis sets. These include the STO-3G,37 CEP-4G,38 LanL2MB,39 and MINI40 minimal basis sets and the 3-21G*, 3-21G(d), 6-31G*, and 6-31G(d) double-ζ basis sets.41,42 Note that our definitions for the “*” and “(d)” notations are somewhat different from those used in the Gaussian program. Thus, for both 3-21G and 6-31G basis sets, “*” indicates the inclusion of polarization functions only for second-row elements, whereas the “(d)” basis sets also include such functions for first-row elements (i.e., all nonhydrogen atoms). In some cases, we have used the cc-pVTZ+d basis set43 (abbreviated as VTZ+d) for benchmark purpose. In particular, B3-LYP/cc-pVTZ+d has been the recommended procedure for geometry optimization and vibrational frequency calculations within W144 and several WnX-type16,45−47 highlevel composite protocols. To obtain theoretical zero-point vibrational energies (ZPVEs), thermal corrections for enthalpies at 298 K (ΔH298), and 298 K entropies (S298), we used standard thermochemical formulas48 together with scaled vibrational Θ frequencies. Specifically, ZPVEs are obtained by ZPVE = ∑ 2

with geometry and vibrational frequencies are often major contributors to the overall uncertainty within a high-level thermochemical protocol.15,16 In comparison, for larger systems for which hybrid-DFT and DH-DFT are the only computationally viable options for the calculation of single-point energies, the use of more-economical methodologies for geometry optimization and vibrational frequency calculations is often unavoidable. In previous studies, we have used DH-DFT to investigate medium-sized fullerenes (up to C320).17,18 For the larger members within the series, the computation of vibrational frequencies was only viable with semiempirical procedures. In that investigation, we find that the use of custom scale factors for the semiempirical methods provides a reliable means for obtaining accurate scaled vibrational frequencies for the calculation of thermochemical quantities (ZPVEs and ΔH298) of fullerenes. In general, however, broadly applicable (i.e., nonspecific) frequency scale factors for thermochemical applications with lower-cost methods (compared to hybrid-DFT) are sparsely available. In the present study, we have determined frequency scale factors for a set of prototypical semiempirical and pureDFT procedures, in order to devise recommendations for their use for geometry optimization and vibrational frequency calculations for obtaining thermochemical quantities for larger systems.

■

COMPUTATIONAL DETAILS

Standard computational chemistry calculations were carried out with Gaussian 0919 and Spartan 16.20 In general, initial geometries were obtained from previous studies.8,17 For the tripeptide used in our assessment, we employed the hybrid systematic/Monte Carlo conformational search algorithm as implemented in Spartan, in conjunction with the MMFF force field,21 to generate an initial set of conformer structures. For geometry optimization and computation of harmonic vibrational frequencies, we have examined the use of a number of semiempirical and pure-DFT methods, as well as several hybrid-DFT procedures for comparison purposes. The semiempirical procedures are AM1,22 PM3,23 and PM6.24 We have investigated B-LYP,25,26 B-PW91,25,27 PBE,28 and N1229 as prototypical pure-DFT methods. For B-LYP and PBE, we have also examined their dispersion-corrected variants.30,31 The hybrid-DFT methods included in the present study are the global-hybrid B3-LYP32 and B3-PW9133 methods and several

where Θ is the vibrational temperature as defined by Θ =

hν kB

in

which h is the Planck constant, ν is the vibrational frequency, and kB is the Boltzmann constant. Vibrational enthalpy is

(

calculated with the equation ΔHT(vib) = R ∑ Θ

1 2

+

1 eΘ/T − 1

)

with R being the gas constant and T being the temperature. B


Article

Journal of Chemical Theory and Computation Table 2. MADs from Benchmark Values for the Various Properties for the LF10 Seta AM1 PM3 PM6 B-LYP B-PW91 PBE B-LYP-D3BJ PBE-D3BJ N12 HSE HISS N12-SX B3-LYP B3-PW91

6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) 6-31G(d) VTZ+d VTZ+d

ZPVE

ΔH298

S298(vib)

S(rot)

Eelec

H0

H298

G298

r

7.3 7.0 10.1 3.6 3.1 3.1 3.7 3.1 2.5 2.9 2.9 2.7 1.8 1.7

0.8 1.1 0.9 0.5 0.4 0.4 0.5 0.4 0.3 0.3 0.3 0.2 0.2 0.2

7.3 7.4 6.1 2.7 2.2 2.1 2.6 2.1 2.0 1.8 2.0 1.7 1.3 1.0

0.2 0.2 0.2 0.5 0.4 0.3 0.5 0.3 0.1 0.1 0.2 0.2 0.1 0.1

9.6 6.9 5.9 9.0 5.3 4.7 8.6 4.6 0.9 0.8 1.6 1.0 1.0 1.0

16.1 12.1 15.1 9.8 6.2 5.9 9.6 5.9 2.8 2.8 4.3 3.0 2.2 2.3

16.4 11.0 14.6 9.6 6.2 5.9 9.5 5.8 2.6 2.6 4.1 2.7 2.0 2.2

16.2 13.0 16.0 9.7 6.1 5.8 9.5 5.7 3.1 3.1 4.6 3.3 2.3 2.4

0.015 0.008 0.010 0.021 0.016 0.016 0.020 0.016 0.020 0.016 0.005 0.005 0.008 0.005

MADs for ZPVE and ΔH298 in kJ mol−1; S298(vib) and S(rot) = vibrational (298 K) and rotational entropy, respectively (J mol−1 K−1); Eelec, H0, H298, and G298 = DSD-PBE-P86/aug′-cc-pVTZ+d vibrationless energy, 0 K enthalpy, 298 K enthalpy, and 298 K free energy, respectively (kJ mol−1) and, r = bond length (Å).

a

previous study,50 it has been shown that the corrections for anharmonicity have comparable magnitudes for B-LYP (a pureDFT) and B3-LYP (a related hybrid-DFT). In that study, both the harmonic and anharmonic frequencies for B-LYP are found to be smaller than those for B3-LYP in a fairly consistent manner. This is in accord with the generally smaller scale factors for hybrid-DFT procedures determined in the present and other previous studies.9,10 General Performance for the LF10 Set of Larger Molecules. To assess the performance of the various theoretical procedures for obtaining geometries and (scaled) vibrational frequencies, we have in a previous study8 compiled the LF10 set of reference quantities taken from experimental data within CCCBDB.51 Because the LF10 set contains larger molecules than those in Z2 and F1, it “amplifies” the deviations and thus enables a more distinct comparison between the different methods. In the present study, we will use this set of systems as a major benchmark for our analysis. The mean absolute deviations (MADs) for the various thermochemical quantities computed with different methods are shown in Table 2. We have previously identified the use of low-level geometries for single-point-energy calculations and differences in ZPVEs to be major contributors to deviations from benchmark thermochemical quantities (e.g., H298 and G298).8 It can be seen that this continues to be the case, as shown by the variations among the more economical methods examined in the present study. Thus, the semiempirical methods have some of the largest MADs of 5.9−9.6 kJ mol−1 for single-point energies (Eelec). With geometries obtained using pure-DFT procedures, the MADs for Eelec also vary widely from 0.9 kJ mol−1 for N12/6-31G(d) to 9.0 kJ mol−1 for B-LYP/6-31G(d). The dispersion-corrected variants do not show significant differences from their DFT-only counterparts for this set of systems. The hybrid-DFT methods shown in Table 2 (both global-hybrid and screened-exchange) consistently yield reasonable MADs of ∼1−2 kJ mol−1. While the deviations in Eelec offer an indirect account for the deviations in geometries, it is of interest to interrogate the structures directly. To this end, we have examined the associated deviations in the bond lengths (Table 2, r). There is not a very strong correlation between the MADs for Eelec with

Vibrational entropy is obtained with the formula ST (vib) = R ∑

(

Θ/T eΘ/T − 1

)

− ln(1 − e−Θ / T ) .

Scale factors for ZPVE, ΔH298 and S298 were determined using literature methodologies.8,9 Thus, the root-mean-square deviation from benchmark values for contributions from each scaled vibrational mode is minimized. Unless otherwise noted, single-point energy calculations were carried out at the DSDPBE-P86 level49 with the aug′-cc-pVTZ+d basis set41 to determine the effect on energies for various geometries. This method is chosen as we consider DH-DFT to be the highestlevel viable for large systems for which the combination of hybrid-DFT and a triple-ζ basis set for optimization and frequency calculations is computationally prohibitive. Energies in the text are given in kJ mol−1, and entropies are given in J mol−1 K−1.

■

RESULTS AND DISCUSSION Optimized Scale Factors for the Various Low-Cost Procedures. The scale factors for a selection of methods examined in the present study are given in Table 1, while the full set of scale factors determined is provided in the Supporting Information. In accordance with previous practice,8,9 we have determined these scale factors using the Z2 and F1 sets of small molecules for which accurate ZPVEs and vibrational frequencies are available. They contain up to ten atoms, consisting of hydrogen and first- and second-row elements. The full set of molecules in these two sets is given in the Supporting Information. Separate scale factors are obtained for several distinct properties, namely fundamental frequencies, low frequencies (≤1000 cm−1), ZPVEs, ΔH298, and S298. In an overall sense, the scale factors for different classes of methods decrease in the order semiempirical > pure-DFT > hybrid-DFT. For example, the average scale factor for ZPVEs for semiempirical (AM1 to PM6) is 1.0161, whereas those for pure-DFT (B-LYP to N12) and hybrid-DFT methods (HSE to B3-PW91) are 0.9991 and 0.9688, respectively. We however note that there can be considerable variations within each class of methods. For instance, the scale factors for AM1 are smaller than those for many pure-DFT methods examined. The primary function of applying scale factors is to correct for anharmonicity and deficiencies in the methodology. In a C


Article

Journal of Chemical Theory and Computation Table 3. Deviations from Benchmark Values for N12/6-31G(d) for the LF10 Seta C6H6 PhF PhCl PhCN cubane CF3CN C2F6 N2O4 PF5 SF6 a

ZPVE

ΔH298

S298(vib)

S(rot)

Eelec

H0

H298

G298

−0.3 3.0 2.8 3.3 7.9 −0.6 −1.7 1.0 −1.6 −3.3

−0.2 −0.3 −0.2 −0.3 −0.7 0.0 0.2 −0.1 0.1 0.6

−1.7 −2.1 −1.8 −2.6 −3.9 −0.1 2.9 −1.7 0.2 2.8

−0.1 −0.1 0.0 −0.1 −0.1 −0.1 0.0 0.1 0.3 0.4

−0.1 0.1 −1.1 −1.1 −0.8 −0.9 −1.2 0.1 −1.9 −1.5

−0.4 3.1 1.7 2.1 7.1 −1.5 −2.9 1.1 −3.4 −4.7

−0.6 2.8 1.4 1.8 6.4 −1.4 −2.7 1.1 −3.3 −4.1

−0.1 3.5 2.0 2.6 7.6 −1.4 −3.6 1.5 −3.4 −5.0

ZPVE and ΔH298 in kJ mol−1, S298(vib) and S(rot) in J mol−1 K−1, Eelec, H0, H298, and G298 in kJ mol−1.

Table 4. MADs from Benchmark Values for the LF10 Set for N12 Used in Combination with a Variety of Basis Setsa,b STO-3G CEP-4G LanL2MB MINI 3-21G* 3-21G(d) 6-31G* 6-31G(d) VTZ+d

ZPVE

ΔH298

S298(vib)

S(rot)

Eelec

H0

H298

G298

7.8 11.4 9.9 9.5 4.9 3.4 5.0 2.5 2.6

1.7 1.9 2.1 1.4 0.4 0.2 0.6 0.3 0.2

12.1 11.9 12.7 8.6 2.4 1.8 3.2 2.0 1.2

1.2 3.1 1.5 2.0 0.2 0.2 0.3 0.1 0.1

76.6 403.8 123.5 217.5 4.1 2.4 5.7 0.9 1.2

73.0 398.0 119.7 220.5 6.8 5.7 6.9 2.8 3.6

73.5 397.5 119.7 220.6 6.6 5.7 6.7 2.6 3.4

71.7 397.6 119.3 219.9 7.0 5.7 7.0 3.1 3.7

a ZPVE and ΔH298 in kJ mol−1; S298(vib) and S(rot) in J mol−1 K−1; Eelec, H0, H298, and G298 in kJ mol−1. bNote that the definitions of the “*” and “(d)” notations are revised from those in the Gaussian program, see Computational Details for more information.

those for r, which indicates that other factors, such as bond angles and the type of bonds that are associated with significant deviations, might also be important contributors to deviations in energies. Nonetheless, we note that small MADs in Eelec tend to be associated with MADs in r that are ≤0.005 Å. The deviations for ZPVEs are most substantial for the semiempirical methods (MADs = 7.0−10.1 kJ mol−1), as one might intuitively anticipate. In contrast to the observations for Eelec, pure-DFT methods (2.5−3.6 kJ mol−1) perform only somewhat less well when compared with hybrid-DFT (1.7−2.9 kJ mol−1). When the deviations for both Eelec and ZPVEs are taken into account, the MADs for H0 show clearer distinctions between the three classes of methods, with ranges of 12.1−16.1, 2.8−9.8, and 2.2−4.3 kJ mol−1, respectively, for semiempirical, pure-DFT, and hybrid-DFT procedures. Going beyond H0, the MADs for H298 and G298 do not show large differences from the H0 values, indicating the relatively small impact of the level of theory on the additional quantities. Overall, among the three semiempirical methods, we do not see a “clear winner”. The use of this type of methods leads to MADs of ∼15 kJ mol−1 for the LF10 set of systems for H0, H298, and G298. At the other end of the scale, the corresponding MADs for the hybrid-DFT procedures in Table 2 are ∼2−5 kJ mol−1. In the middle of the pack, the variations between the pure-DFT methods are fairly substantial. Importantly, we find that the N12/6-31G(d) procedure performs exceptionally well for a pure-DFT method. Its MADs are generally quite comparable to more expensive hybrid-DFT procedures. Performance of N12/6-31G(d) for Individual Molecules in the LF10 Set. As we have pointed out in the previous section, the N12/6-31G(d) procedure performs particularly well among the pure-DFT procedures and is even competitive with those hybrid-DFT methods. We thus deem it to be a cost-

effective means for geometry optimization and vibrational frequency calculations for medium-sized to large molecules. In this regard, we will further inspect its performance for each species within the LF10 set, with the aim to identify its strengths and shortcomings. The deviations from benchmark values are shown in Table 3. We can see that, for the ten molecules in the LF10 set, N12/ 6-31G(d) performs reasonably well for most of them. One notable exception is for the ZPVE of cubane, where the calculated value deviates from benchmark by +7.9 kJ mol−1. For other species, the deviations for ZPVE are ∼±3 kJ mol−1. For other quantities, the deviations are also quite modest, with values of ∼±1 kJ mol−1 for ΔH298, ±4 J mol−1 K−1 for S (corresponds to ∼±1 kJ mol−1 at 298 K), and ∼±2 kJ mol−1 for Eelec. The combinations of these thermochemical quantities lead to deviations of ∼±3 kJ mol−1 for H0, H298, and G298 in general. Larger deviations occur in two instances, specifically for cubane (∼+7 kJ mol−1) and SF6 (∼−5 kJ mol−1). Overall, N12/6-31G(d) appears to be fairly adequate for the purpose of obtaining geometries and vibrational frequencies for “typical” molecules, but one may encounter larger deviations for strained species such as cubane. Effect of Basis Set on the Performance of N12. While N12/6-31G(d) appears to offer a cost-effective means for optimizing geometries and calculating vibrational frequencies for applications to medium-sized systems, for larger systems, even more economical methods are desirable. To this end, it is of interest to explore the possibility of further reducing the size of the basis set in order to extract untapped computational efficiency. We have thus examined the combination of N12 with several smaller double-ζ basis sets, namely 6-31G*, 321G(d), and 3-21G* (see Computational Details for their definitions in the present study), and the minimal basis sets D


Article


expect the deficiencies in the computed geometries and vibrational frequencies would also cancel to a considerable extent. As we have noted in the Introduction, we have previously employed a series of isodesmic-type reactions to accurately compute the heats of formation of medium-sized fullerenes.17 We used semiempirical methods to obtain the vibrational frequencies for some of the larger species, together with validated custom scale factors for fullerenes. How adequate are the generic scale factors determined in the present study when applied to this kind of reactions? To seek the answer to this question, we have examined several isodesmic-type reactions of fullerenes with a number of methodologies for geometry optimization and vibrational frequency calculations, with the highest-level employed being B3-LYP/6-31G(d). We obtained single-point energies at the DSD-PBE-P86/cc-pVTZ (for C60 and C90), B3-PW91/cc-pVTZ (for C60, C90, and C180), and B3PW91/cc-pVDZ (for C60, C90, C180, and C320) levels in order to assess the energetic effect of variations in geometries. We will first look at the absolute deviations in Eelec and ZPVE for the various fullerenes for the few selected methods (Table 5), in order to establish a baseline for comparison when we next discuss the results for the isodesmic-type reactions. Before we compare the different methods for obtaining these two

STO-3G, CEP-4G, LanL2MB, and MINI. For comparison, we have also obtained values for the larger cc-pVTZ+d basis set. The results for the LF10 set are shown in Table 4. Before we discuss the performance of the smaller basis sets, let us first compare the results obtained with 6-31G(d) to those with the larger cc-pVTZ+d basis set. We can see that the MADs for 6-31G(d) are comparable to or are sometimes (perhaps fortuitously) smaller than those for cc-pVTZ+d. These results further support the use of N12/6-31G(d) as a cost-effective means for geometry optimization and vibrational frequency calculations. Among the smaller basis sets, 6-31G* and 3-21G(d) can be compared directly to 6-31G(d). The 6-31G* basis set contains polarization functions only for second-row and later nonhydrogen elements, whereas these functions are included in 631G(d) for all supported elements starting from Li. It is apparent from the notable differences between the two basis sets that polarization functions play an important role for geometry optimizations as well as vibrational frequency calculations. Thus, the MAD for Eelec is 0.9 kJ mol−1 for 631G(d), but it is significantly larger for 6-31G* (5.7 kJ mol−1). Likewise, the exclusion of polarization functions for first-row elements in 6-31G* also leads to a considerable deterioration in the performance for computation of ZPVEs [MAD = 5.0 kJ mol−1, versus 2.5 kJ mol−1 for 6-31G(d)]. The 3-21G(d) basis set differs from 6-31G(d) in that the 631G component is replaced by 3-21G, but the same set of polarization functions is employed. We can see that the use of a more modest split-valence basis set leads to a smaller deterioration in the accuracy (e.g., from 0.9 to 2.4 kJ mol−1 for Eelec). In a similar manner, the performance of N12/3-21G* is quite comparable to that for N12/6-31G*. Overall, even though the MADs for the smaller double-ζ basis sets are larger than those for 6-31G(d), the worsening in performance is byno-means excessive. Thus, we deem the use of, e.g., 3-21G* to be appropriate when the calculation of a large system necessitates the use of a smaller basis set than 6-31G(d). In contrast, the use of minimal basis sets leads to very dramatic changes in the accuracy. Even for the best performing one, namely STO-3G, the MAD of 76.6 kJ mol−1 for Eelec is massive. The largest MAD for Eelec is associated with the CEP4G basis set (403.8 kJ mol−1). Interestingly, despite very large deviations in the optimized geometries (as indicated by the MADs for Eelec), the deviations for other properties are not overly excessive (MADs ∼ 10 kJ mol−1 for ZPVE). At any rate, the MADs for N12 with minimal basis sets are comparable to or larger than those for semiempirical methods (Table 2). They thus offer no apparent benefit for the purpose of geometry optimization and vibrational frequency calculations when compared with computationally less-demanding semiempirical procedures. Evaluation with Fullerenes and the Utility of Isodesmic-Type Reactions. So far we have examined the deviations of the various quantities in absolute terms, i.e., deviations associated with each molecule on its own. In some cases, such deviations would be carried over to calculated thermochemical quantities unchanged. These include, for example, calculated atomization energies and associated heats of formation, where the atomic species do not contribute to, e.g., the Eelec and the ZPVE components. In other cases, one may expect cancellation of deviations to different degrees. In particular, for isodesmic-type reactions52,53 that are formulated with the aim of maximizing such benefit, one may reasonably

Table 5. Deviations (kJ mol−1) from Eelec [DSD-PBE-P86/ cc-pVTZ, (B3-PW91/cc-pVTZ) and B3-PW91/cc-pVDZ] and ZPVE Values Obtained with Benchmark [B3-LYP/631G(d)] Geometries for Selected Medium-Sized Fullerenes and Computer Time Consumed (s) for the Vibrational Frequency Calculations C60 AM1

PM6

N12

3-21G*

N12

6-31G(d)

PBE

6-31G(d)

AM1 PM6 N12 N12 PBE AM1 PM6 N12 N12 PBE B3-LYP E

3-21G* 6-31G(d) 6-31G(d)

3-21G* 6-31G(d) 6-31G(d) 6-31G(d)

C90

C180

deviations for Eelec (kJ mol−1) 24.7 35.9 (−21.7) (−27.3) (−66.6) 20.7 31.9 75.5 41.7 58.1 (−44.0) (−60.8) (−137.6) 34.2 45.7 108.1 1.5 2.9 (−56.1) (−79.4) (−177.7) 6.6 12.8 31.3 2.4 3.4 (−19.8) (−23.2) (−64.5) 4.3 9.1 19.8 14.1 24.4 (−57.0) (−77.5) (−172.6) 9.6 14.8 28.6 deviations for ZPVE (kJ mol−1) 120.1 172.5 124.8 177.3 −16.8 −28.4 −7.1 −12.4 −0.2 −0.3 time for frequency calculation (s) 190 615 1850 138 287 2025 16490 26518 125477 71858 119748 106662 168556 181265 322766

C320

129.5

179.6

60.9

36.2

49.5

67142 77851


Article

Journal of Chemical Theory and Computation quantities, we note that, in cases for which DSD-PBE-P86 Eelec values are available, among the two lower-level methods, B3PW91/cc-pVDZ in fact yields deviations that are in better agreement with the DSD-PBE-P86 values. Thus, for this part of our discussion, we will focus on the deviations for B3-PW91/ cc-pVDZ Eelec for which values for all fullerenes are accessible. For Eelec, there are notable deviations from the values obtained with benchmark geometries. The deviations with various DFT geometries are reasonably close to one another for C60, but the gap widens with increasing size of the system. Thus, for C320, the deviations for N12/6-31G(d), PBE/631G(d), and N12/3-21G* are 36.2, 49.5, and 60.9 kJ mol−1, respectively. In comparison, the deviations for Eelec for values obtained with semiempirical geometries are considerably larger. For PM6, they vary from 34.2 kJ mol−1 for C60 to 179.6 kJ mol−1 for C320. The computations of vibrational frequencies with the full set of methodologies were only viable for C60 and C90 with our current computational resources. Nonetheless, similar qualitative trends can also be observed for deviations for ZPVE. In the face of the challenge for carrying out vibrational frequency calculations for the larger fullerenes, it is instructive to briefly digress and examine the resources required for these computations. The total computer time taken is also shown in Table 5. For one of the lowest-cost methods, namely AM1, the computation for C60 was completed in 190 s. This figure grew to just over 67,000 s for C320. With the least-costly DFT method (N12/3-21G*), the vibrational frequency calculation for C60 already took over 16,000 s, and it extended to ∼125,000 s for C180. Further improvements in the basis set and the methodology from N12/3-21G* to N12/6-31G(d) [and PBE/ 6-31G(d)] to B3-LYP/6-31G(d) led to notable (but perhaps somewhat less dramatic) increases in the computer time consumed, with the most expensive calculation taking ∼323,000 s. Let us now return to the subject of accuracy. The large deviations shown in Table 5 speak for the importance of cancellation of errors if one were to employ the computationally more efficient procedures to obtain geometries and frequencies. To this end, we have examined four isodesmictype reactions of the fullerenes (Table 6). Remarkably, none of

on the two sides of each equation because of larger differences in the sizes of the fullerenes involved. These results provide a glimpse of the utility for the use of low-cost procedures to obtain geometries and thermochemical quantities for chemical reactions of medium-sized to large systems, yielding reasonable accuracies. Relative Energies for Peptide Isomers. Fullerenes and the small systems in the Z2, F1, and LF10 sets are structurally fairly rigid. For geometry optimizations (and associated vibrational frequency calculations), molecules with substantial conformational freedoms may represent a different challenge for relatively low-level theoretical procedures. To this end, we have assessed the performance of AM1, PM6 N12/3-21G*, N12/6-31G(d), and PBE/6-31G(d) using the tripeptide Nformyl-L-Met-L-Leu-L-Phe-OH (PDB code 1Q7O). In this case, we also employ B3-LYP/6-31G(d) geometries and thermochemical quantities as benchmarks, as is the case for the above assessment using fullerenes. DSD-PBE-P86/cc-pVTZ singlepoint energies are employed to evaluate the energetic effects of using different underlying geometries. Our preliminary conformational search yields over 300 conformers, among which over 90% of the population is represented by the 16 lowest-energy structures at the MMFF level (i.e., the method used in the conformational search). For a more thorough assessment, we employ the 50 lowest-energy conformers in subsequent benchmark. The lowest-energy structure obtained with DSD-PBE-P86/cc-pVTZ//B3-LYP/631G(d) is shown in Figure 1, which involves several hydrogen bonds. We note that, importantly, all other procedures also give the same conformer as the lowest-energy one.

Table 6. Deviations (kJ mol−1) from ΔEelec and ΔZPVE (in Parentheses) Values Obtained with Benchmark Geometries for Isodesmic-Type Reactions of Medium-Sized Fullerenes 1.5 C60 → C90 AM1 PM6 N12 N12 PBE

3-21G* 6-31G(d) 6-31G(d)

0.8 −5.7 2.9 2.6 0.4

(−7.7) (−9.9) (−3.2) (−1.8) (−0.1)

3 C60 → C180

2 C90 → C180

2 C180 → C320 + 2/3 C60

13.3 5.5 11.4 6.9 −0.2

11.6 16.8 5.6 1.7 −1.1

−7.6 −13.8 2.7 −0.6 −1.3

Figure 1. B3-LYP/6-31G(d)-optimized structure for the tripeptide 1Q7O. Non-hydrogen-bonded hydrogen atoms are omitted for clarity.

Table 7 shows the deviations for DSD-PBE-P86 Eelec and B3LYP ZPVE for this conformer. The qualitative observations are generally in accord with those in Table 5. Thus, the deviations between the N12/6-31G(d) values and those for B3-LYP/631G(d) are reasonably small (−5.1 and −1.3 kJ mol−1, respectively, for Eelec and ZPVE). The deviations are somewhat larger when the smaller 3-21G* basis set is used or when N12 is replaced with PBE. For instance, there is a 16.5 kJ mol−1 deviation for ZPVE for N12/3-21G* and a 11.6 kJ mol−1 deviation for Eelec for PBE/6-31G(d). For AM1 and PM6, the deviations in Eelec are quite large (85.7 and 64.7 kJ mol−1, respectively), which reflects considerable structural differences from the B3-LYP/6-31G(d) geometry. Aligning the semi-

the deviations for the contributions to reaction energies is very large, especially when one considers the excessive deviations in absolute terms (Table 5). For the 1.5 C60 → C90 reaction, the deviations in ΔEelec are just a few kJ mol−1, while the corresponding deviations for ΔZPVE are larger. The deviations for the other three reactions are larger. This may be in part due to a stronger reliance on cancellations of deviations associated with larger deviations in absolute terms and in part due to poorer cancellations resulting from less-well balanced chemistry F


Article

Journal of Chemical Theory and Computation Table 7. Deviations (kJ mol−1) from Benchmark [DSD-PBEP86/cc-pVTZ//B3-LYP/6-31G(d)] Eelec and ZPVE Values for the Lowest-Energy Conformer of Tripeptide 1Q7O, Mean Absolute Deviations (kJ mol−1) for Relative Energies (ΔE) of Its 50 Lowest-Energy Conformers, and Deviations (kJ mol−1) for the Correction for Boltzmann Distributions AM1 PM6 N12 N12 PBE

3-21G* 6-31G(d) 6-31G(d)

Eelec

ZPVE

MAD(ΔE)

EBoltzmann

85.7 64.7 16.0 −5.1 11.6

2.6 3.7 16.5 −1.3 −3.1

5.9 6.5 3.3 1.0 1.6

0.4 1.1 0.5 0.0 −0.2

DSD-PBE-P86/cc-pVTZ//B3-LYP/6-31G(d) energies and the associated B3-LYP ZPVEs as our benchmark. The results are shown in Table 8. Table 8. Mean Absolute Deviations (kJ mol−1) from Benchmark [DSD-PBE-P86/cc-pVTZ//B3-LYP/6-31G(d)] Eelec and ZPVE Values for the W4-11 Set of Molecules (mol) and for the Associated W4-11-RE Set of Relative Energies (RE) AM1 PM6 N12 N12 PBE

empirical structures with the B3-LYP one shows that there are relatively small variations in strongly hydrogen-bonded regions, but larger structural differences occur in other locations, in particular in the side chain of Phe. Curiously, despite the large deviations in Eelec, for AM1 and PM6 the corresponding deviations in ZPVE are just 2.6 and 3.7 kJ mol−1, respectively. As we have seen earlier in our assessment with fullerenes, the calculation of relative energies can benefit substantially from cancellation of deviations in the absolute energies. This is also the case for the peptide conformational energies in Table 7. For example, while the deviation in Eelec for AM1 is 85.7 kJ mol−1 for the lowest-energy conformer (the corresponding MAD for all 50 conformers is 82.8 kJ mol−1), the MAD for relative energies (ΔE) between the conformers is just 5.9 kJ mol−1. In a similar manner, the MADs for ΔE for N12/3-21G*, N12/631G(d), and PBE/6-31G(d) are smaller (3.3, 1.0, and 1.6 kJ mol−1, respectively) in comparison with the larger deviations in Eelec. From a thermochemical point of view, the significance of having a large number of conformers is the associated Boltzmann contribution (EBoltzmann) to the energy. At the B3LYP/6-31G(d) level, this contribution amounts to 2.3 kJ mol−1 with the 50 low-energy conformers of the tripeptide. When this quantity is obtained with N12/6-31G(d), there is virtually no deviation from the B3-LYP value. The deviations for other methods are somewhat larger. Overall, the results of the assessments with fullerenes and the small peptide further substantiate the use of N12/6-31G(d) as a robust method for geometry optimization and the calculation of vibrational frequencies for thermochemical applications. In addition, they show that even-lower-cost methods may also be useful for these purposes in a somewhat limited fashion. Performance for an Extensive Set of Reactions. The two sets of systems considered above, namely isodesmic reactions of fullerenes and relative energies of peptide isomers, do not involve substantial changes in the bonding structures between reactants and products. These examples fittingly illustrate the utility of low-cost methods for geometry optimizations and vibrational frequency calculations when used in combination with reaction schemes that maximize cancellation of systematic deviations. Nonetheless, it is also of importance to assess more general situations in which such cancellations are less extensive. Thus, we have examined the performance of the various semiempirical and DFT methods for obtaining geometries and ZPVEs for the W4-11-RE set54 of reaction energies. This compilation comprises over 10,000 reactions, which are generated systematically using the diverse set of 140 molecules in the W4-11 set.55 As with our assessments using fullerenes and the peptide, we employ

3-21G* 6-31G(d) 6-31G(d)

Eelec(mol)

ZPVE(mol)

Eelec(RE)

ZPVE(RE)

24.0 12.1 3.8 0.6 2.4

2.0 2.7 1.2 0.3 0.3

21.6 21.2 6.2 1.6 4.6

3.0 4.4 2.0 0.5 0.7

Let us first examine the MADs for Eelec and ZPVE for the molecules themselves, i.e., Eelec(mol) and ZPVE(mol). For the semiempirical methods, the MADs in Eelec(mol) are quite substantial, with values of 24.0 and 12.1 kJ mol−1, respectively, for AM1 and PM6. The corresponding MADs for the DFT methods are much smaller, ranging from 0.6 kJ mol−1 for N12/ 6-31G(d) to 3.8 kJ mol−1 for N12/3-21G*. In a similar manner, we find that the MADs for ZPVEs are smaller for the DFT methods (∼1 kJ mol−1) than those for the semiempirical methods (∼2−3 kJ mol−1), but the differences are much less pronounced. These observations are consistent with other results in the present study. How do the deviations in individual molecules combine in the wide range of reactions in the W4-11-RE set? For Eelec(RE), the MADs for the semiempirical methods are ∼21 kJ mol−1. Thus, while the Eelec(RE) MAD with AM1 geometries is somewhat smaller than the corresponding Eelec(mol) value, for PM6, the MAD for Eelec(RE) is considerably larger than that for Eelec(mol). We can also see that, for the three DFT procedures examined, the Eelec(RE) MADs are all larger than the corresponding Eelec(mol) values. The larger MADs for Eelec(RE) than those for Eelec(mol) can be readily rationalized, in qualitative terms, by considering standard error-propagation principles and the general form of A + B → C + D reactions in the W4-11-RE set. In this regard, the comparable Eelec(mol) and Eelec(RE) MADs for AM1 represent an especially fortuitous outcome. For ZPVEs, the MADs for the reaction energies are also larger than those for the molecules themselves. Overall, N12/6-31G(d) appears to remain as a fairly reasonable method for geometry optimization and vibrational frequency calculations in a more general sense. Our results show that the deviations in individual molecular energies resulting from the use of approximate geometries and frequencies would accumulate in chemical reactions, which is more-or-less in accordance with general error-propagation principles. Nonetheless, if one were able to exploit systematic cancellation of deviations, it would enable the use of very lowcost methods without excessive deterioration in the accuracy.

■

CONCLUDING REMARKS In the present study, we have determined frequency scale factors for obtaining thermochemical quantities [zero-point vibrational energies (ZPVE), thermal corrections for 298 K enthalpies (ΔH298), and 298 K entropies (S298)] using a variety of low-cost quantum chemistry procedures. These include semiempirical, pure-DFT, and screened-exchange DFT methG


Article

Journal of Chemical Theory and Computation Funding

ods. We have assessed their performance for computing molecular geometries and thermochemical quantities associated with vibrational frequencies. In absolute terms, for the LF10 set of small to medium-sized molecules (e.g., benzene and cubane), all procedures perform reasonably well. Major deviations are found for calculated ZPVEs and single-point energies (Eelec) associated with the use of approximate geometries. These then become main contributors to the deviations for total enthalpies (H0 and H298) and free energies (G298). The semiempirical methods show similar accuracy to one another, with mean absolute deviations (MADs) of ∼15 kJ mol−1 for H0, H298, and G298. For the DFT procedures, hybrid DFT generally performs better than pure DFT. In spite of this, the N12 pure functional shows exceptionally good performance (MADs for H0, H298, and G298 ∼ 3 kJ mol−1) that is comparable to hybrid functionals. An examination of the effect of the basis set on the performance of N12 shows that the moderate 6-31G(d) basis set provides geometries and ZPVEs that are comparable to those obtained with the larger cc-pVTZ+d basis set. While the use of smaller ones leads to worsening of the accuracy, the deteriorations are not excessive for typical double-ζ basis sets. On the other hand, using minimal basis sets leads to very large MADs for the calculated thermochemical quantities. Overall, we deem the use of N12/3-21G* and N12/6-31G(d) to be a cost-effective means for geometry optimization and vibrational frequency calculations. Further testing with isodesmic-type reactions of mediumsized to large fullerenes and conformational energies of a small peptide shows that, although the deviations for Eelec and ZPVEs can be large in absolute terms (up to >100 kJ mol−1), those for contributions to reaction energies are markedly reduced (∼10 kJ mol−1). This enables the use of low-cost methodologies such as the semiempirical AM1 procedure to obtain geometries and vibrational frequencies (and hence ΔZPVE, ΔH298, etc.) with reasonable accuracy for large systems for which the use of more expensive procedures is computationally prohibitive. For reactions that do not involve substantial cancellation of deviations, the deviations for individual species accumulate in accordance with error-propagation principle. In these situations, N12/6-31G(d) would be preferred over lower-cost methods.

■

We gratefully acknowledge research funding from the Japan Society for the Promotion of Science (JSPS) (Grant number 16H07074001) and generous grants of computer time from the RIKEN Advanced Center for Computing and Communication (ACCC), Japan, Institute for Molecular Science (IMS), Japan, and National Computational Infrastructure (NCI), Australia. Notes

The author declares no competing financial interest.

■

(1) Karton, A. A Computational Chemist’s Guide to Accurate Thermochemistry for Organic Molecules. WIREs Comput. Mol. Sci. 2016, 6, 292−310. (2) Chan, B. How to Computationally Calculate Thermochemical Properties Objectively, Accurately, and as Economically as Possible. Pure Appl. Chem. 2017, 89, 699−713. (3) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gn Theory. WIREs Comput. Mol. Sci. 2011, 1, 810−825. (4) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (5) Peverati, R.; Truhlar, D. G. Screened-Exchange Density Functionals with Broad Accuracy for Chemistry and Solid-State Physics. Phys. Chem. Chem. Phys. 2012, 14, 16187−16191. (6) Mardirossian, N.; Head-Gordon, M. Thirty Years of Density Functional Theory in Computational Chemistry: An Overview and Extensive Assessment of 200 Density Functionals. Mol. Phys. 2017, 115, 2315−2372. (7) Brémond, É.; Savarese, M.; Su, N. Q.; Pérez-Jiménez, Á . J.; Xu, X.; Sancho-García, J. C.; Adamo, C. Benchmarking Density Functionals on Structural Parameters of Small/Medium-Sized Organic Molecules. J. Chem. Theory Comput. 2016, 12, 459−465. (8) Chan, B.; Radom, L. Frequency Scale Factors for Some DoubleHybrid Density Functional Theory Procedures: Accurate Thermochemical Components for High-Level Composite Protocols. J. Chem. Theory Comput. 2016, 12, 3774−3780. (9) Scott, A. P.; Radom, L. Harmonic Vibrational Frequencies: An Evaluation Of Hartree−Fock, Møller−Plesset, Quadratic Configuration Interaction, Density Functional Theory, and Semiempirical Scale Factors. J. Phys. Chem. 1996, 100, 16502−16513. (10) Merrick, J. P.; Moran, D.; Radom, L. An Evaluation of Harmonic Vibrational Frequency Scale Factors. J. Phys. Chem. A 2007, 111, 11683−11700. (11) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. Computational Thermochemistry: Scale Factor Databases and Scale Factors for Vibrational Frequencies Obtained from Electronic Model Chemistries. J. Chem. Theory Comput. 2010, 6, 2872−2887. (12) Laury, M. L.; Boesch, S. E.; Haken, I.; Sinha, P.; Wheeler, R. A.; Wilson, A. K. Harmonic Vibrational Frequencies: Scale Factors for Pure, Hybrid, Hybrid Meta, and Double-Hybrid Functionals in Conjunction with Correlation Consistent Basis Sets. J. Comput. Chem. 2011, 32, 2339−2347. (13) Laury, M. L.; Carlson, M. J.; Wilson, A. K. Vibrational Frequency Scale Factors for Density Functional Theory and the Polarization Consistent Basis Sets. J. Comput. Chem. 2012, 33, 2380− 2387. (14) Kesharwani, M. K.; Brauer, B.; Martin, J. M. L. Frequency and Zero-Point Vibrational Energy Scale Factors for Double-Hybrid Density Functionals (and Other Selected Methods): Can Anharmonic Force Fields be Avoided? J. Phys. Chem. A 2015, 119, 1701−1714. (15) Chan, B.; Deng, J.; Radom, L. G4(MP2)-6X: A Cost-Effective Improvement to G4(MP2). J. Chem. Theory Comput. 2011, 7, 112− 120.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00721. Complete sets of frequency scale factors, deviations from benchmark values for ZPVEs and DSD-PBE-P86 singlepoint energies, and calculated energies and thermochemical quantities for fullerenes, the tripeptide and W4-11RE set of molecules (Tables S1−S6) (XLSX) B3-LYP/6-31G(d)-optimized structures for the tripeptide N-formyl-L-Met-L-Leu-L-Phe-OH in XYZ format, incorporated into a single compressed archive (ZIP)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Bun Chan: 0000-0002-0082-5497 H


Article


(35) Henderson, T. M.; Izmaylov, A. F.; Scuseria, G. E.; Savin, A. J. Chem. Theory Comput. 2008, 4, 1254−1262. (36) Chan, B.; Gill, P. M. W.; Radom, L. Performance of GradientCorrected and Hybrid Density Functional Theory: Role of the Underlying Local Density Approximation and the Gradient Correction. J. Chem. Theory Comput. 2012, 8, 4899−4906. (37) Hehre, W. J.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of SlaterType Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657−2664. (38) Stevens, W. J.; Basch, H.; Krauss, M. Compact Effective Potentials and Efficient Shared-Exponent Basis Sets for the First- and Second-Row Atoms. J. Chem. Phys. 1984, 81, 6026−6033. (39) Hay, P. J.; Wadt, W. R. Ab Initio Effective Core Potentials for Molecular Calculations. Potentials for the Transition Metal Atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270−283. (40) Huzinaga, S.; Andzelm, J.; Klobukowski, M.; Radzio-Andzelm, E.; Sakai, Y.; Tatewaki, H. Gaussian Basis Sets for Molecular Calculations; Elsevier: Amsterdam, 1984. (41) Binkley, J. S.; Pople, J. A.; Hehre, W. J. Self-Consistent Molecular Orbital Methods. 21. Small Split-Valence Basis Sets for First-Row Elements. J. Am. Chem. Soc. 1980, 102, 939−947. (42) Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules. J. Chem. Phys. 1971, 54, 724−728. (43) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (44) Martin, J. M. L.; De Oliveira, G. Towards Standard Methods for Benchmark Quality Ab Initio Thermochemistry − W1 and W2 Theory. J. Chem. Phys. 1999, 111, 1843−1856. (45) Chan, B.; Radom, L. W1X-1 and W1X-2: W1-Quality Accuracy with an Order of Magnitude Reduction in Computational Cost. J. Chem. Theory Comput. 2012, 8, 4259−4269. (46) Chan, B.; Radom, L. W3X: A Cost-Effective Post-CCSD(T) Composite Procedure. J. Chem. Theory Comput. 2013, 9, 4769−4778. (47) Chan, B. Unification of the W1X and G4 (MP2)-6X Composite Protocols. J. Chem. Theory Comput. 2017, 13, 2642−2649. (48) McQuarrie, D. A.; Simon, J. D. Molecular Thermodynamics; University Science Books: Sausalito, CA, 1999. (49) Kozuch, S.; Martin, J. M. L. DSD-PBEP86: In Search of the Best Double-Hybrid DFT with Spin-Component Scaled MP2 and Dispersion Corrections. Phys. Chem. Chem. Phys. 2011, 13, 20104− 20107. (50) Ringholm, M.; Jonsson, D.; Bast, R.; Gao, B.; Thorvaldsen, A. J.; Ekström, U.; Helgaker, T.; Ruud, K. Analytic Cubic and Quartic Force Fields Using Density-Functional Theory. J. Chem. Phys. 2014, 140, 034103. (51) NIST Computational Chemistry Comparison and Benchmark Database; Johnson, R. D., III, Ed.; NIST Standard Reference Database Number 101 Release 18; National Institute of Standards and Technology: Gaithersburg, MD, 2016. (52) Wheeler, S. E.; Houk, K. N.; Schleyer, P. v. R.; Allen, W. D. A Hierarchy of Homodesmotic Reactions for Thermochemistry. J. Am. Chem. Soc. 2009, 131, 2547−2560. (53) Ramabhadran, R. O.; Raghavachari, K. Theoretical Thermochemistry for Organic Molecules: Development of the Generalized Connectivity-Based Hierarchy. J. Chem. Theory Comput. 2011, 7, 2094−2103. (54) Margraf, J. T.; Ranasinghe, D.; Bartlett, R. J. Automatic Generation of Reaction Energy Databases from Highly Accurate Atomization Energy Benchmark Sets. Phys. Chem. Chem. Phys. 2017, 19, 9798−9805. (55) Karton, A.; Daon, S.; Martin, J. M. L. W4-11: A HighConfidence Benchmark Dataset for Computational Thermochemistry Derived from First-Principles W4 Data. Chem. Phys. Lett. 2011, 510, 165−178.

(16) Chan, B.; Radom, L. W2X and W3X-L: Cost-Effective Approximations to W2 and W4 with kJ mol−1 Accuracy. J. Chem. Theory Comput. 2015, 11, 2109−2119. (17) Karton, A.; Chan, B.; Raghavachari, K.; Radom, L. Evaluation of the Heats of Formation of Corannulene and C60 by Means of HighLevel Theoretical Procedures. J. Phys. Chem. A 2013, 117, 1834−1842. (18) Chan, B.; Kawashima, Y.; Katouda, M.; Nakajima, T.; Hirao, K. From C60 to Infinity: Large-Scale Quantum Chemistry Calculations of the Heats of Formation of Higher Fullerenes. J. Am. Chem. Soc. 2016, 138, 1420−1429. (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision E.01; Gaussian, Inc.: Wallingford, CT, 2009. (20) Hehre, W. J. Spartan 16; Wavefunction, Inc.: Irvine, CA, 2016. (21) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and Performance of MMFF94. J. Comput. Chem. 1996, 17, 490−519. (22) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (23) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods I. Method. J. Comput. Chem. 1989, 10, 209−220. (24) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods V: Modification of NDDO Approximations and Application to 70 Elements. J. Mol. Model. 2007, 13, 1173−1213. (25) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behaviour. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (26) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (27) Perdew, J. P. In Electronic Structure of Solids ’91; Ziesche, P., Eschrig, H., Eds.; Akademie Verlag: Berlin, 1991; p 11. (28) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (29) Peverati, R.; Truhlar, D. G. Exchange−Correlation Functional with Good Accuracy for Both Structural And Energetic Properties while Depending Only on the Density and its Gradient. J. Chem. Theory Comput. 2012, 8, 2310−2319. (30) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104−1−19. (31) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465. (32) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (33) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (34) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid Functionals Based on a Screened Coulomb Potential. J. Chem. Phys. 2003, 118, 8207−8215. I


Use of Low-Cost Quantum Chemistry Procedures for Geometry

Recommend Documents