FULL PAPER

WWW.C-CHEM.ORG

Implementation of Nuclear Gradients of Range-Separated Hybrid Density Functionals and Benchmarking on Rotational Constants for Organic Molecules Tobias Risthaus,[a,b] Marc Steinmetz,[a] and Stefan Grimme*[a] We have implemented the nuclear gradient for several rangeseparated hybrid density functionals in the general quantum chemistry code ORCA. To benchmark the performance, we have used a recently published set of back-corrected gas phase rotational constants, which we extended by three molecules. In our evaluation, CAM-B3LYP-D3 and xB97X-D3 show great accuracy, and are surpassed by B2PLYP-D3 only. Lowercost alternatives to quadruple-f basis set-based calculations, among them a smaller basis set and the use of resolution-of-

the-identity approaches, are assessed and shown to yield acceptable deviations. In addition, the Hartree-Fock-based back-correction method is compared to a density functional theory alternative, which largely shows consistency between the two. A new, well-performing, spin-component scaled C 2014 Wiley MP2 variant is designed and discussed, as well. V Periodicals, Inc.

Introduction

optimization of the total energy of the system with respect to its nuclear coordinates. Very recently, some of us published a new test set called ROT25 aiming to provide a compilation of reliable gas phase rotational constants that are directly related to the entire geometry of the molecule, including conformational aspects, under consideration.[24] The main objective of this work is to extend that test set and apply RSH functionals to the resulting set. For context, some of the best performers of the previous investigation will be considered as well. That evaluation focused on the bare functional performance (which may include a dispersion correction) and contained efforts to exclude basis set effects and other sources of error. For that reason, no methods to reduce computational cost such as the “resolution-of-the-identity” (RI) approximation (see below) had been used. However, the use of these methods is nowadays fairly commonplace and as such, consideration in the current work seems warranted. Specifically, we investigate the effects of the RIJK[25] and RIJCOSX[26,27] variants of the RI approximation. In RIJK, appropriate auxiliary basis sets are used to substitute both Coulomb (J) and exchange integrals (K, X) ðlmjrkÞ as used in the KohnSham/Fock matrix by auxiliary three-center and two-center electron repulsion integrals ðlmjPÞ; ðPjQÞ:

Density functional theory (DFT) has emerged as a widely used tool in electronic structure theory.[1,2] It often outclasses more semiempirical approaches like NNDO (for instance PM6[3]) in terms of accuracy and while the advantage of semiempirical methods due to the lower cost is still appreciable in very large systems or molecular dynamics applications, advances in computational power make DFT the method of choice for the calculation of medium to large molecular equilibrium structures and their energetics. For this, hybrid density functionals (DF), which mix density-based and Fock exchange, have shown great success.[4] Range-separated hybrids (RSH) are DF that conceptually split the two-electron operator for the exchange to allow a different treatment at different interelectronic ranges. Several realizations and fields of applications have arisen, among them the use of the error function as the switching function[5,6] or a Gaussianbased two-electron operator[7] as well as the choice of having Fock exchange dominate either the short[8] or the long range[5] of electronic interaction with complementary (meta)-GGA exchange as needed. The extension of the concept of rangeseparation to DFT correlation has also been investigated.[9] In the following, we will focus on DFs that use the error function to yield up to 100% Fock exchange in the long-range limit. In previous years, these functionals, also dubbed longrange corrected (LC), have garnered increased attention, especially in the field of electronically excited states. For some recent applications, see Refs. 10–17 and for more indepth reviews, see Refs. 18,19. Truhlar and coworkers,[20,21] Chai and coworkers,[22] and Head-Gordon and coworkers[23] have each made efforts to design RSH functionals suitable for ground state properties, notably thermochemistry and weak intermolecular interactions. Properties often treated as an aside or completely neglected are structural equilibrium parameters (the “geometry” performance) as obtained by

DOI: 10.1002/jcc.23649

This article was published online on June 13, 2014. An error was subsequently identified. This notice is included in the online version to indicate that it has been corrected on June 25, 2014.

[a] T. Risthaus, M. Steinmetz, S. Grimme Mulliken Center for Theoretical Chemistry, Institut f€ ur Physikalische und Theoretische Chemie, Universit€ at Bonn, Beringstr. 4, D-53115 Bonn, Germany E-mail: [email protected] [b] T. Risthaus International NRW Graduate School of Chemistry, Wilhelm-Klemm-Str. 10, D-48149, M€ unster, Germany C 2014 Wiley Periodicals, Inc. V

Journal of Computational Chemistry 2014, 35, 1509–1516

1509

FULL PAPER

WWW.C-CHEM.ORG

ðlmjrkÞ 

X ðlmjPÞðPjQÞ21 ðQjkrÞ

(1)

P;Q

An additional error arises due to the finite size of the auxiliary basis set fPg, which is the trade-off for a significant drop in computational expense, in our experience a factor of 10 in small to medium-sized cases with standard auxiliary basis sets. In RIJCOSX (the RIJ chain-of-spheres exchange approach), the Coulomb integrals are approximated in the same way as in RIJK. The exchange integrals are treated differently in that only the integration over the spatial coordinates of one electron is performed analytically (up to the point where the Boys function is used), whereas the second one is done numerically on an appropriate grid indexed by g using grid point weights wg: ðlmjrkÞ 

X

Xlg Akr ðrg ÞXmg ;

(2)

g

where 1

Xlg 5ðwg Þ2  lðrg Þ; ð kðrÞrðrÞ dr: Akr ðrg Þ5 jr2rg j

(3) (4)

The finite grid size introduces an additional error but it also allows for spatial resolution and efficient screening of the exchange integrals leading to substantial savings in larger systems, especially when using basis sets with low maximum angular momentum. In the course of the currently presented work, we have implemented the nuclear gradient for several RSH functionals both within the framework of the RIJCOSX approximation and in the standard way (that is only subject to approximation in terms of prescreening by the Cauchy-Schwarz inequality, significance cut-offs etc.). This allows us to look into the numerical precision of the approximation. Since the usage together with valence quadruple-f basis sets is probably not the most common case, we also conduct investigations on the valence triple-f level. The remainder of the article is organized as follows. First, we give a detailed account of our used methodology. Second, we present the newly extended test set, which we dub ROT34, and the results for several DFs and computational settings. We then proceed to compare the back-corrections applied to rotational constants to equilibrium values via vibration-rotation interaction constants at LDA and HF level and recount the problems surrounding a molecule initially intended to be a member of the test set. We also investigate a new SCS-MP2 parametrization motivated by previous findings. Finally, we draw some conclusions.

Methodology Most optimizations were performed with a locally modified version of ORCA 2.9.[28] The nuclear gradients implemented (both standard and RIJCOSX) were tested against numerical derivatives as provided by the !NumGrad option. The function1510

Journal of Computational Chemistry 2014, 35, 1509–1516

als used are xB97, xB97X,[6] xB97X-D3,[22] CAM-B3LYP,[29] LCBLYP,[5] BHLYP (also known as Becke’s Half and Half functional or BHandHLYP),[30] B3LYP,[31] PW6B95,[32] and B2PLYP.[33,34] In some cases, the D3 dispersion correction[35–37] was applied using the Becke-Johnson damping function.[38,39] We used the def2-QZVP basis set,[40] and the def2-TZVP basis set,[41,42] which are abbreviated to QZ and TZ throughout. Unless noted otherwise, the very large DFT quadrature grids Grid7 in ORCA and 6 in TURBOMOLE were used. When the RIJCOSX approximation was used, the def2-QZVP/ J and def2-TZVP/J auxiliary basis sets were used as appropriate.[43,44] In addition, the def2-QZVP/C basis set was used for the RI-PT2 part of the B2PLYP calculations.[45] The moderately sized GridX5 was used for the semi-numerical integration of the Fock exchange, which also used the S-fitting option.[27] For these calculations, the DFT grid was lowered to Grid5. When the RIJK approximation was used, the def2-QZVP/JK and def2-TZVP/JK auxiliary basis sets were used as appropriate.[46] All RIJK-based calculations were performed using TURBOMOLE 6.5.[47] Following the same procedure as previously,[24] the effective rotational constant Bk0 along a principal rotation axis k 2 fA; B; Cg is related to the equilibrium rotational constant Bke and the vibration-rotation interaction constants akr for normal mode r as described by: Bk0 5Bke 2

1X k a : 2 r r

(5)

Each akr is obtained from perturbation theory,[48,49] using results obtained with CFOUR[50] at the HF level. The basis set used was Dunning’s DZ[51] as obtained from the EMSL Basis Set Exchange database.[52] For comparison, results were also obtained at the SVWN(V)/DZ level[53,54] using Gaussian.[55] The resulting correction will be denoted DBvib and called “anharmonic,” although it contains effects of the harmonic (zero-point) average over r12 instead of over r12 as well. Finite e 0 temperature effects in the experiment (supersonic jet expansion conditions corresponding to very low temperatures) do not affect the reference values due to the high resolution of the microwave experiments. For more discussion and a comparison to high-level computational results, see Ref 24. In general, we define the error as Bke 2Bkcalc , which means that a negative error corresponds to a rotational constant that is too large, i.e. the spatial extent of the structure is too small. The following error measurements are used: the mean deviation (MD) and the mean absolute deviation in MHz, the mean relative deviation (MRD), the relative error range (MINMAXR), and the estimated relative standard deviation SD in percent, defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 34  k u1 X Be 2Bkcalc SD 5t ; (6) 2MRD 33 k51 Bke where Bkcalc and Bke are the kth calculated and back-corrected rotational constant. Finally, we list the root mean square deviation (RMSD) between two structures as well as its mean WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Figure 1. Ball-and-stick representations of the molecules in the entire test set and the conformer under investigation. The color code is: carbon gray, hydrogen white, nitrogen blue, and oxygen red.

(RMSD /A˚) between the geometries obtained with certain ˚. methods for the whole test set in A

The rotational constants, the anharmonic corrections thereto, and the resulting reference values are listed in Table 1. Inspection reveals that the DBvib are between 20.8 and 21.3%

Results and Discussion Extension of the test set The test set is extended by the rotational constants of the molecules R-(1)-limonene 10,[56] (2)-lupinine 11,[57] and proline derivative Ac-Pro-NH2 12,[58] which are depicted along the members of the previous test set in Figure 1. From the restricted pool of available high-level experimental data, the molecules were chosen for size and the new structural motifs they provide. Limonene 10 is a hydrocarbon that features double bonds in a scaffold that is in between 1 and 4 in terms of flexibility. Lupinine 11 contains a very pronounced internal hydrogen bond with a nitrogen acceptor, a motif missing from the set thus far. The proline derivative 12 sports a fivemembered ring that is not more or less planar as in 6 and 9, but shows the puckering typical of saturated rings that size.

Table 1. Experimental rotational constants B0 as taken from the literature, calculated anharmonic corrections DBvib on HF/DZ level, and resulting equilibrium rotational constants Be, which serve as reference values for this work. Molecule

Principal axis

10

A B C A B C A B C

11

12

B0

DBvib

Be

3058.0 717.0 679.3 1414.1 811.7 671.5 1511.3 1060.2 714.5

228.2 26.7 25.7 218.0 28.8 27.9 211.9 210.3 25.4

3086.2 723.7 685.0 1432.1 820.5 679.4 1523.2 1070.5 719.9

All values in MHz.

Journal of Computational Chemistry 2014, 35, 1509–1516

1511

FULL PAPER

WWW.C-CHEM.ORG

Table 2. Statistical data for the results of the ROT34 test set for the methods tested without integral approximation using the def2-QZVP basis set. Method

MD/MHz MAD/MHz MRD/% MINMAXR/% SD/%

CAM-B3LYP-D3/QZ CAM-B3LYP/QZ LC-BLYP/QZ xB97/QZ xB97X/QZ xB97X-D3/QZ B3LYP-D3/QZ PW6B95-D3/QZ BHLYP-D3/QZ B2PLYP-D3/QZ B2PLYP/QZ HF-D3/QZ HF/QZ

0.1 2.6 27.8 6.0 3.4 0.9 10.9 28.6 210.4 2.6 5.7 230.2 5.1

3.5 5.4 8.0 6.1 4.8 3.5 10.9 8.7 10.4 3.5 6.2 30.2 11.3

0.06 0.41 20.54 0.55 0.39 0.06 0.83 20.61 20.69 0.21 0.58 22.19 0.73

1.25 2.96 1.92 1.34 1.75 1.40 1.41 1.89 1.57 1.09 2.11 1.83 5.48

0.29 0.59 0.33 0.36 0.41 0.31 0.36 0.37 0.32 0.23 0.46 0.46 1.15

See Methodology section for the definition of the error and the measures used.

of B0 with an average of 0.97%. Comparison with Table 2 of Ref. 24 shows that these corrections are of similar size compared to 1–9. For the entire ROT34 test set, the range of corrections is 20.2 to 21.6%. We reiterate that these corrections are comparable to those found on higher levels of wave function theory for smaller molecules.[59] Hence, we consider the backcorrected reference values to be trustworthy and the uncertainty at about 60.2% to be smaller than the deviations of most of the methods we test in this work.

Benchmark Results Comparison of RSHs to Standard Methods. In order to get a clear picture of the performance of RSH functionals, we have evaluated them first without integral approximations such as RIJCOSX or RIJK on the new test set. For context, we have also re-evaluated some of the standard and some of the best performing methods investigated in Ref. 24. The statistical data are listed in Table 2. Besides the emergence of the so-called Jacob’s Ladder,[60] one of the main findings of our previous work was that inclusion of a dispersion correction[36,37] significantly improved the results and that a deterioration of the error for a rotational constant was found only in very few cases. We now observe the same significant effect for the RSH functionals. In CAM-B3LYP-D3, the best RSH performer in the benchmark, the MD and MRD drops to almost zero when D3 is added. The error range and the SD are halved compared to CAM-B3LYP. Since the underlying functionals are similarly constructed, a comparable effect is found when matching up xB97X and xB97X-D3, in spite of their different parametrization. The latter and CAM-B3LYP-D3 show the lowest SD among the RSH functionals and are beaten only by B2PLYP-D3 in this regard. The difference in the other error measures is on the order of the estimated accuracy of the reference values. B3LYP-D3 and PW6B95-D3 are consistently beaten by both dispersion-corrected RSH functionals. The other four RSH functionals do not perform as well: while often showing similar precision measures (SD, MINMAXR) 1512

Journal of Computational Chemistry 2014, 35, 1509–1516

as the best performers, the increased MRD and MD values are disappointing. LC-BLYP shows a negative MD and MRD, meaning that the structures obtained are on average too small. In this regard, it is comparable to HF. However, the similarly constructed xB97 (which also varies the amount of Fock exchange from 0 to 100%) shows the opposite behavior. This is indicative of the fact that the form and parametrization of the complementary GGA functional matters in this application. This is different from applications to excited states, where the dominating influence is the amount of Fock exchange.[19,61] A reviewer suggested that the LYP functional may be responsible for the atypical behavior of LC-BLYP, which could be traced to its performance for weakly interacting systems.[62] This is substantiated by the similar performance of dispersionuncorrected BHLYP. However, plain CAM-B3LYP shows contrary behavior despite including 81% of LYP. The evidence is thus not conclusive on this matter. To sum up, no correlation between the RSH model and the performance could be discerned. Comparison of the statistical data for the methods tested both in this work and previously (B2PLYP(-D3), PW6B95-D3, B3LYP-D3, HF(-D3)) reveals little difference. This is perhaps not surprising, given that almost three quarters of the data is identical. The newly tested BHLYP-D3 functional shows a consistent performance with a SD of 0.32%, but a fairly large deviation from the reference with a MRD of 20.69%. Removal of the D3 correction (see Supporting Information) brings the MRD closer to zero, but at the cost of a doubled SD. Before moving on to choice of basis set and integral approximations, we shall briefly consider whether the inclusion of the dispersion correction D3 consistently improves the performance of a functional. Unfortunately for this investigation, the only new RSH case where inclusion of the D3 correction without changing the underlying functional parametrization is conceptually possible is CAM-B3LYP, since each functional from the xB97 family has a different parametrization, that is, xB97X-D3 is not simply xB97X plus a D3 correction. We thus plot the relative errors of CAM-B3LYP and CAM-B3LYP-D3 in Figure 2. While there are some cases of deterioration upon inclusion of the D3 correction, for instance for the problematic molecule 2 (which will be discussed in more detail in the second half of the following section), a general improvement is observed.

Figure 2. Plot of the relative errors of CAM-B3LYP and CAM-B3LYP-D3 for each rotational constant considered in this work. Within each group of columns the largest rotational constant is on the left, followed by the medium one. See Methodology section for the definition of the error. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

Table 3. Selected statistical data for the results of the ROT34 test set comparing the def2-QZVP and def2-TZVP basis set.

Method

QZ MRD/%

SD/%

TZ MRD/%

SD/%

CAM-B3LYP-D3 CAM-B3LYP LC-BLYP xB97 xB97X xB97X-D3 B3LYP-D3 PW6B95-D3 BHLYP-D3 HF-D3 HF B2PLYP-D3

0.06 0.41 20.54 0.55 0.39 0.06 0.83 20.61 20.69 22.19 0.73 0.21

0.29 0.59 0.33 0.36 0.41 0.31 0.36 0.37 0.32 0.46 1.15 0.23

0.07 0.42 20.55 0.52 0.39 0.06 0.82 20.60 20.66 22.15 0.75 0.41

0.26 0.51 0.30 0.31 0.35 0.34 0.36 0.33 0.30 0.46 1.12 0.23

See Methodology section for the definition of the error and the measures used.

Results for Computational Levels with Diminished Cost. One of the easiest and most commonly used ways to lower the computational cost associated with either a WFT or DFT method is to use a smaller basis set. Also, a commonly held wisdom is that geometries are not very sensitive to the basis set. To test this statement more broadly than was done in our preceding work, the methods evaluated were also used together with the def2-TZVP basis set. The MRD and SD data for this evaluation are presented in Table 3. The differences observed for all methods evaluated are negligibly small with the exception of B2PLYP-D3, for which a sharp rise in the MRD is observed. This corroborates the earlier finding that larger, that is, quadruplef, basis sets must be recommended for standard use with perturbative methods such as this double-hybrid DF.[4] Surprisingly, the SD stays low. Across all methods, the error range (see Supporting Information) does not change significantly nor increase consistently. We thus conclude that for the SCF methods considered herein, the error introduced by using the TZ instead of the QZ basis set is very small. An alternative and nowadays common way to reduce the computational cost is the use of integral approximations to

the two-electron integrals used within the SCF and gradient calculation. We have successfully implemented the RSH functionals for use with the RIJCOSX approximation and list the statistical data obtained from the geometry optimizations performed with them in Table 4. For context, the same data on the other methods is listed as well as data pertaining to the RIJK approximation, where available. In general, when moving from NoRI to RIJCOSX, the MRD changes by about 0.1% in any direction. However, the actual changes are larger, as indicated by the mostly sizable increases in the SD. The strongest increases are found for xB97X-D3 (from 0.31% to 0.73%) and HF-D3 (from 0.46% to 1.06%). The error ranges MINMAXR, showing the effect on individual rotational constants, are typically increased to twice the NoRI size (see Supporting Information). An exception to this behavior is B2PLYP(-D3), wherein the RI-PT2 gradient may be dominating. This is also reflected by its low RMSD . For the global hybrids and HF, one finds a positive correlation between the amount of Fock exchange (the quantity approximated by COSX) and the RMSD . On that scale, the RSH functionals are ranked a bit below BHLYP-D3. The results for the triple-f basis set are very similar; we thus forego discussion here and list the results in the Supporting Information. In order to investigate the influence of the COSX grid, we have looked to the most affected molecule 2 to provide insight (Table 5, Fig. 3). Inspection revealed that the driving force behind the RMSD is a certain dihedral angle. It thus becomes clear that one needs to use rather large COSX grid setups to obtain a good agreement with the unapproximated calculation if the potential energy surface is rather flat, as in this case. We have investigated xB97X-D3 in the same manner and found virtually identical results (not shown here). The errors introduced by the RIJK approximation are significantly lower. The RMSD is negligibly small compared to both typical fluctuations caused by different starting structures and the less restrictive convergence criteria commonly used. Accordingly, the changes in MRD and SD are very small. By cross-comparing the performance of RIJK and RIJCOSX, we can conclude that the RIJ approximation shares the characteristics

˚ for the RI-approximated calculation with respect to the Table 4. Selected statistical data for the results of the ROT34 test set, including the RMSD in A unapproximated calculation (NoRI).

Method

NoRI MRD/%

SD/%

RIJCOSX MRD/%

SD/%

RMSD /A˚

RIJK MRD/%

CAM-B3LYP-D3/QZ CAM-B3LYP/QZ LC-BLYP/QZ xB97/QZ xB97X/QZ xB97X-D3/QZ B3LYP-D3/QZ PW6B95-D3/QZ BHLYP-D3/QZ B2PLYP-D3/QZ B2PLYP/QZ HF-D3/QZ HF/QZ

0.06 0.41 20.54 0.55 0.39 0.06 0.83 20.61 20.69 0.21 0.58 22.19 0.73

0.29 0.59 0.33 0.36 0.41 0.31 0.36 0.37 0.32 0.23 0.46 0.46 1.15

20.03 0.31 20.65 0.52 0.33 0.01 0.77 20.77 20.79 0.21 0.59 22.30 0.50

0.49 0.44 0.47 0.41 0.45 0.73 0.51 0.58 0.67 0.30 0.42 1.06 1.11

0.030 0.030 0.035 0.033 0.034 0.032 0.021 0.025 0.037 0.010 0.008 0.067 0.066

– – – – – – 0.82 20.63 20.69 – – 22.19 0.72

SD/%

RMSD /A˚

– – – – – – 0.36 0.39 0.33

– – – – – – 0.002 0.002 0.001 – – 0.001 0.003

– 0.46 1.16

See Methodology section for the definition of the error and the measures used.

Journal of Computational Chemistry 2014, 35, 1509–1516

1513

FULL PAPER

WWW.C-CHEM.ORG

Table 5. Results of the COSX grid size investigation for 2, performed on the B3LYP-D3/QZ level. Gridsize

RMSD/A˚

3602U=

NoRI GridX7 GridX6 GridX5 GridX4 GridXS1 GridXS2

– 0.03 0.09 0.08 0.09 0.37 0.37

57.9 58.4 60.1 59.9 59.9 69.7 69.8

Table 6. Experimental rotational constants B0 in MHz, taken from the literature,[56–58,63–72] calculated DBvib corrections in MHz on the HF/DZ and SVWN/DZ level. 

The NoRI result serves as reference. The investigated dihedral angle U is illustrated in Figure 3.

Molecule

Principal axis

B0/MHz

1

A B C A B C A B C A B C A A B C A B C A B C A B C A B C A B C A B C A B C

4246.2 1386.5 1121.7 3280.9 713.0 690.1 3043.1 1274.7 1235.7 2730.1 2650.5 2631.7 2314.9 1455.7 760.5 575.4 1163.2 650.6 450.1 1156.1 762.6 509.0 855.7 745.5 508.5 3058.0 717.0 679.3 1414.1 811.7 671.5 1511.3 1060.2 714.5 4437.2 612.2 595.8

2

3

of the RIK approximation. In summary, an additional and significant error introduced by COSX due to its grid is clearly apparent, whereas RIJ and RIK are not a source of significant error.

4

Comparison of Corrections at the HF and SVWN Levels

7

In order to corroborate the anharmonic corrections beyond the favorable comparison with high-level coupled-cluster results of other groups referenced in our previous work,[24] we have calculated them also at the SVWN/DZ DFT level. This level was chosen due to its relatively low cost and the difference afforded compared to HF, since SVWN contains some LDA correlation through its functional form. The resulting corrections are listed along the used HF/DZ corrections in Table 6. It is heartening to see that two methods, that are theoretically rather different, yield similar results. Nevertheless, it is also apparent that if differences in the structures are obtained (as indicated by the RMSD), different corrections DBvib will arise. In the case of 5, a slightly different conformer is found with SVWN/DZ. Still, we are confident that the correct conformer was chosen for the test set. Otherwise, no qualitative deviations are found for the molecules included in the test set, and we thus are confident that the DBvib values are correct within the margins mentioned before (60.2% of Be). For the reason of qualitative disagreement, we have excluded the ethyl valerate molecule in the conformation depicted in Figure 4 (EtVl). The absolute difference between the HF and SVWN

Figure 3. Ball-and-stick representations of 2 and the dihedral angle U in blue listed in Table 5. The color code is: carbon gray, hydrogen white, and oxygen red.

1514

HF/DZ

Journal of Computational Chemistry 2014, 35, 1509–1516

5 6

8

9

10

11

12

EtVl

SVWN/DZ

DBvib/MHz 247.7 29.4 28.5 241.6 26.8 27.9 228.0 210.3 213.0 225.8 225.1 221.6 222.0 28.5 27.7 25.2 22.5 210.6 23.9 210.2 25.0 24.0 26.8 28.7 25.2 228.2 26.7 25.7 218.0 28.8 27.9 211.9 210.3 25.4 10.5 26.6 26.4

246.7 29.4 28.4 233.8 29.1 26.6 229.5 214.3 211.0 226.8 225.8 222.5 245.1 216.4 24.9 25.2 29.6 25.8 23.9 28.8 25.2 23.7 26.3 28.6 24.9 242.8 24.2 27.3 217.8 27.5 26.1 212.6 210.0 25.3 214.1 27.6 26.7

RMSD/A˚ 0.032

0.359

0.086

0.023

0.529 0.123

0.126

0.152

0.051

0.048

0.116

0.130

0.100

We also list the individual RMSD in A˚ for each system between the HF and the SVWN geometry.

corrections is not larger than in the case of the first principal axis of 10 and, accordingly, is not the reason for exclusion. Instead, EtVl was dropped because a positive DBvib value on the HF level has not been observed by us before. In addition, it was determined that the test set would be sufficiently large to yield statistically meaningful data without inclusion of this strange case.

Figure 4. Ball-and-stick representation of ethyl valerate (EtVl). The color code is: carbon gray, hydrogen white, and oxygen red. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

WWW.CHEMISTRYVIEWS.COM

FULL PAPER

WWW.C-CHEM.ORG

A SCS-MP2 Parametrization for Gas-Phase Structures When examining the results of the ROT25 benchmark, we observed that the SD for both MP2 and SCS-MP2 was low and of similar size (0.36% and 0.25%, respectively) and that the MD were very close, but of opposite sign (20.44% and 0.40%, respectively). For this reason, we decided to investigate the performance of a SCS-MP2 variant with a spin-component parametrization halfway between MP2 and SCS-MP2. The resulting parameters are cOS 5 1.1 and cSS 5 23. For a review of other SCS-MP2 variants and an in-depth discussion of the theoretical foundations of spin-component scaling, see Ref. 73. The deviations found for each rotational constant using MP2, SCS-MP2 and the new variant are plotted in Figure 5. It is perhaps surprising that the errors of the new variant lie very consistently halfway between MP2 and SCS-MP2. Accordingly, its statistical data, as listed in Table 7, is very good. The MD and MRD are close to zero and the SD is extremely small. Since the conformational energetics of both MP2 and SCSMP2 are decent, as evidenced by the results for the conformational subsets in the GMTKN30 database,[4] we believe the new variant to be a useful tool for the field of assigning conformations to experimental results for rotational constants. However, it must be stressed that due to the high basis set dependence of all MP2-based methods, the use of this variant with a smaller basis set is not recommended since significantly worse results are to be expected. In this regard, we refer the reader to the large shift in MD when going from the def2QZVP to the def2-TZVPP basis set listed in Table 4 of Ref. 24.

Conclusions In this work, we have extended our previous test set for gas phase rotational constants useful to test the performance of methods for isolated molecular structures to 34 data points, yielding the ROT34 set. To that end, we have calculated corrections covering anharmonic and zero-point effects for three molecules. We consider this undertaking to be successful because the corrections fall within the limits of our expectations and compare well with corrections obtained for the other cases. Additionally, we have investigated corrections obtained on two different computational levels which largely showed agreement and thus raised our confidence in the ref-

Table 7. Statistical data for the results of the ROT34 test set for the three MP2 variants using the def2-QZVP basis set. Method

MD/MHz MAD/MHz MRD/% MINMAXR/% SD/%

MP2 SCS-MP2 SCS(1:1; 0: 6)-MP2

25.0 5.2 20.3

5.5 5.4 1.8

20.45 0.40 20.05

1.45 1.38 1.12

0.32 0.24 0.22

See Methodology section for a definition of the error and the measures used.

erence values. This comparison also led to the discarding of one molecule we considered for inclusion in the test set because a qualitative disagreement was found. We then used the resulting test set to extensively evaluate the performance of several RSH DFs. Similarly to our previous work, we found that a dispersion correction is vital for a good agreement with the reference values. xB97X-D3 and CAMB3LYP-D3 are shown to be the best performers on the hybrid DF level and are only surpassed by B2PLYP-D3 in terms of consistency. We have also looked at common methodologies to lower computational cost. Among them, the smaller def2-TZVP basis set instead of the near-complete def2-QZVP basis set introduced negligible uncertainty for the SCF methods. The RIJCOSX method, which combines a standard RI approach for the Coulomb terms with a seminumerical integration for the Fock exchange, yields errors that are tolerable, but show a clear dependence on the amount of Fock exchange included. RIJK on the other hand introduces only minimal errors at somewhat larger computational demands. Finally, we have tested a spin-component scaled variant of MP2, that shows a minimal mean deviation and a tiny standard deviation comparable to MP2 and SCS-MP2, and we thus recommend it for conformer assignments if use of WFT methods over DFT is desired.

Acknowledgments T.R. would like to thank the International NRW Graduate School of Chemistry for support. S.G. and M.S. express their gratitude to the Deutsche Forschungsgemeinschaft for continued support. Keywords: density functional theory  benchmark  molecular structure  range-separated hybrids  spin-component scaled MP2

How to cite this article: T. Risthaus, M. Steinmetz, S. Grimme. J. Comput. Chem. 2014, 35, 1509–1516. DOI: 10.1002/jcc.23649

]

 Figure 5. Plot of the relative errors of MP2, SCS-MP2 and SCS(1:1; 0:6)-MP2 for each rotational constant considered in this work. Within each group of columns the largest rotational constant is on the left, followed by the medium one. See Methodology section for the definition of the error. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

Additional Supporting Information may be found in the online version of this article.

[1] F. Jensen, Introduction to Computational Chemistry; Wiley: Chicester, UK, 2006. [2] W. Koch, M. C. Holthausen, A Chemist’s Guide to Density Functional Theory; Wiley-VCH: Weinheim, Germany, 2001. [3] J. J. P. Stewart, J. Mol. Mod. 2007, 13, 1173. [4] L. Goerigk, S. Grimme, Phys. Chem. Chem. Phys. 2011, 13, 6670. [5] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 2001, 115, 3540.

Journal of Computational Chemistry 2014, 35, 1509–1516

1515

FULL PAPER [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47]

[48] [49] [50]

1516

WWW.C-CHEM.ORG

J.-D. Chai, M. Head-Gordon, J. Chem. Phys. 2008, 128, 084106. J.-W. Song, K. Yamashita, K. Hirao, J. Chem. Phys. 2011, 135, 071103. J. Heyd, G. E. Scuseria, M. Ernzerhof, J. Chem. Phys. 2003, 118, 8207. A. Stoyanova, A. M. Teale, J. Toulouse, T. Helgaker, E. Fromager, J. Chem. Phys. 2013, 139, 134113. P. Dev, S. Agrawal, N. J. English, J. Chem. Phys. 2012, 136, 224301. K. A. Nguyen, P. N. Day, R. Pachter, J. Chem. Phys. 2011, 135, 074109. N. Holmgaard List, J. M. Olsen, T. Rocha-Rinza, O. Christiansen, J. Kongsted, Int. J. Quantum Chem. 2012, 112, 789. D. W. Silverstein, N. Govind, H. J. J. van Dam, L. Jensen, J. Chem. Theory Comput. 2013, 9, 5490. R. D. Eithiraj, K. R. Geethalakshmi, Chem. Phys. Lett. 2013, 585, 138. B. Moore II, J. Autschbach, J. Chem. Theory Comput. 2013, 9, 4991. N. De Mitri, S. Monti, G. Prampolini, V. Barone, J. Chem. Theory Comput. 2013, 9, 4507. N. Shao, Q. Wu, Phys. Chem. Chem. Phys. 2014, 16, 6700. D. Jacquemin, B. Mennucci, C. Adamo, Phys. Chem. Chem. Phys. 2011, 13, 16987. R. Baer, E. Livshits, U. Salzner, Annu. Rev. Phys. Chem. 2010, 61, 85. R. Peverati, D. G. Truhlar, J. Phys. Chem. Lett. 2011, 2, 2810. R. Peverati, D. G. Truhlar, Phys. Chem. Chem. Phys. 2012, 14, 16187. Y.-S. Lin, G.-D. Li, S.-P. Mao, J.-D. Chai, J. Chem. Theory Comput. 2013, 9, 263. N. Mardirossian, M. Head-Gordon, Phys. Chem. Chem. Phys. 2014, 16, 9904. S. Grimme, M. Steinmetz, Phys. Chem. Chem. Phys. 2013, 15, 16031. F. Weigend, Phys. Chem. Chem. Phys. 2002, 4, 4285. F. Neese, F. Wennmohs, A. Hansen, U. Becker, Chem. Phys. 2009, 356, 98. R. Izsak, F. Neese, J. Chem. Phys. 2011, 135, 144105. F. Neese, WIREs Comput. Mol. Sci. 2012, 2, 73. T. Yanai, D. P. Tew, N. C. Handy, Chem. Phys. Lett 2004, 393, 51. A. D. Becke, J. Chem. Phys. 1993, 98, 1372. A. D. Becke, J. Chem. Phys. 1993, 98, 5648. Y. Zhao, D. G. Truhlar, J. Phys. Chem. A 2005, 109, 5656. S. Grimme, J. Chem. Phys. 2006, 124, 034108. F. Neese, T. Schwabe, S. Grimme, J. Chem. Phys. 2007, 126, 124115. S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132, 154104. J. Klimes, A. Michaelides, J. Chem. Phys. 2012, 137, 120901. S. Grimme, WIREs Comput. Mol. Sci. 2011, 1, 211. S. Grimme, S. Ehrlich, L. Goerigk, J. Comput. Chem. 2011, 32, 1456. E. R. Johnson, A. D. Becke, J. Chem. Phys. 2006, 124, 174104. F. Weigend, F. Furche, R. Ahlrichs, J. Chem. Phys. 2003, 119, 12753. F. Weigend, R. Ahlrichs, Phys. Chem. Chem. Phys. 2005, 7, 3297. F. Weigend, M. H€aser, H. Patzelt, R. Ahlrichs, Chem. Phys. Lett. 1997, 294, 143. F. Weigend, Phys. Chem. Chem. Phys. 2006, 8, 1057. K. Eichkorn, F. Weigend, O. Treutler, R. Ahlrichs, Theor. Chem. Acc. 1997, 97, 119. C. H€attig, Phys. Chem. Chem. Phys. 2005, 7, 59. F. Weigend, J. Comput. Chem. 2008, 29, 167. R. Ahlrichs, M. K. Armbruster, R. A. Bachorz, M. B€ar, H.-P. Baron, R. Bauernschmitt, F. A. Bischoff, S. B€ ocker, N. Crawford, P. Deglmann, F. D. Sala, M. Diedenhofen, M. Ehrig, K. Eichkorn, S. Elliott, F. Furche, A. Gl€ oß, F. Haase, M. H€aser, C. H€attig, A. Hellweg, S. H€ ofener, H. Horn, C. Huber, U. Huniar, M. Kattannek, W. Klopper, A. K€ ohn, C. K€ olmel, M. € Kollwitz, K. May, P. Nava, C. Ochsenfeld, H. Ohm, M. Pabst, H. Patzelt, D. Rappoport, O. Rubner, A. Sch€afer, U. Schneider, M. Sierka, D. P. Tew, O. Treutler, B. Unterreiner, M. von Arnim, F. Weigend, P. Weis, H. Weiss, N. Winter, TURBOMOLE 6.5; Universit€at Karlsruhe, 2009. Available at: http://www.turbomole.com. D. Papousek, M. R. Aliev, Molecular Vibrational-Rotational Spectra; Elsevier: Amsterdam, 1982. I. M. Mills, in Molecular Spectroscopy: Modern Research, K. N. Rao, C. W. Mathews, Eds.; Academic Press: New York, 1972; Vol. 1, pp. 115–140. J. F. Stanton, J. Gauss, M. E. Harding, P.G. Szalay with contributions from A. A. Auer, R. J. Bartlett, U. Benedikt, C. Berger, D. E. Bernholdt, Y.

Journal of Computational Chemistry 2014, 35, 1509–1516

[51] [52] [53] [54] [55]

[56] [57] [58] [59] [60]

[61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73]

J. Bomble, L. Cheng, O. Christiansen, M. Heckert, O. Heun, C. Huber, T.C. Jagau, D. Jonsson, J. Jus elius, K. Klein, W. J. Lauderdale, D. A. Matthews, T. Metzroth, L. A. M€ uck, D. P. O’Neill, D. R. Price, E. Prochnow, C. Puzzarini, K. Ruud, F. Schiffmann, W. Schwalbach, S. Stopkowicz, A. Tajti, J. Vazquez, F. Wang, J. D. Watts, the integral packages MOLECULE (J. Alml€ of and P. R. Taylor), PROPS (P. R. Taylor), ABACUS (T. Helgaker, H. J. Aa. Jensen, P. Jïrgensen, J. Olsen), and ECP routines by A. V. Mitin, C. van W€ ullen. CFOUR, Coupled-Cluster techniques for Computational Chemistry, a quantum-chemical program package. For the current version, see http://www.cfour.de. T. H. Dunning, Jr., J. Chem. Phys. 1970, 53, 2823. K. L. Schuchardt, B. T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, T. L. Windus, J. Chem. Inf. Model. 2007, 47, 1045. J. C. Slater, Phys. Rev. 1951, 81, 385. S. H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 1980, 58, 1200. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, J. V. Dannenberg, S. Dapprich, A. D. Daniels, O. Ortiz, J. Cioslowski, D. J. Fox, Gaussian 09, Rev. A.02; Gaussian Inc.: Pittsburgh, PA. J. R. A. Moreno, T. R. Huet, J. J. L. Gonzalez, Struct. Chem. 2013, 24, 1163. M. K. Jahn, D. Dewald, M. Vallejo-L opez, E. J. Cocinero, A. Lesarri, J.-U. Grabow, J. Phys. Chem. A 2013, 117, 13673. C. Cabezas, M. Varela, J. L. Alonso, ChemPhysChem 2013, 14, 2539. C. Puzzarini, M. Heckert, . Gauss, J. Chem. Phys. 2008, 128, 194108. J. P. Perdew, K. Schmidt, in Density Functional Theory and Its Applications to Materials, AIP Conference Proceedings, Vol. 577; V. E. Van Doren, C. Van Alsenoy, P. Geerlings, Eds.; American Institute of Physics, Melville, NY, USA, 2001. T. Risthaus, A. Hansen, S. Grimme, Phys. Chem. Chem. Phys. (in press), DOI: 10.1039/C3CP54517B. M. Kamiya, T. Tsuneda, K. Hirao, J. Chem. Phys. 2002, 117, 6010. F. J. Lovas, R. D. Suenram, J. Phys. Chem. Ref. Data 1989, 18, 1245. L. W. Sutikdja, D. Jelisavac, W. Stahl, I. Kleiner, Mol. Phys. 2012, 110, 2883. Y. Zhao, H. Mouhib, W. Stahl, J. Phys. Chem. A 2013, 117, 311. R. N. Nandi, M. D. Harmony, A. E. Howard, S. W. Staley, J. Chem. Phys. 1983, 78, 3560. H. V. L. Nguyen, R. Kannengießer, W. Stahl, Phys. Chem. Chem. Phys. 2012, 14, 11753. I. Pe~ na, A. M. Daly, C. Cabezas, S. Mata, C. Berm udez, A. Ni~ no, J. C. L opez, J.-U. Grabow, J. L. Alonso, J. Phys. Chem. Lett. 2013, 4, 65. C. Cabezas, M. Varela, I. Pe~ na, J. C. L opez, J. L. Alonso, Phys. Chem. Chem. Phys. 2012, 14, 13618. C. Cabezas, J. L. Alonso, J. C. L opez, S. Mata, Angew. Chem. Int. Ed. 2012, 51, 1375. H. Mouhib, W. Stahl, M. L€ uthy, M. B€ uchel, P. Kraft, Angew. Chem. Int. Ed. 2011, 50, 5576. H. Mouhib, W. Stahl, ChemPhysChem 2012, 13, 1297. S. Grimme, L. Goerigk, R. F. Fink, WIREs Comput. Mol. Sci. 2012, 2, 886.

Received: 2 April 2014 Revised: 16 May 2014 Accepted: 21 May 2014 Published online on 13 June 2014

WWW.CHEMISTRYVIEWS.COM