Construction of a Spin-Component Scaled Dual-Hybrid Random

Jan 4, 2017 - Among the local field corrections, exchange-only kernels have been proven very practical because this approach allows an analytic evalua...
2 downloads 6 Views 932KB Size
Subscriber access provided by Fudan University

Article

Construction of a spin-component scaled dual-hybrid random phase approximation Pál Dániel Mezei, Gabor I. Csonka, Adrienn Ruzsinszky, and Mihaly Kallay J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01140 • Publication Date (Web): 04 Jan 2017 Downloaded from http://pubs.acs.org on January 5, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Construction of a spin-component scaled dual-hybrid random phase approximation Pál D. Mezei,1* Gábor I. Csonka2, Adrienn Ruzsinszky,3 and Mihály Kállay1 MTA-BME Lendület Quantum Chemistry Research Group, Department of Physical Chemistry and Materials Science, Budapest University of Technology and Economics, H-1521 Budapest, Hungary 1

2

Department of Inorganic and Analytical Chemistry, Budapest University of Technology and Economics, H-1521 Budapest, Hungary 3

Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA

*

E-mail: [email protected]

Abstract Recently, we have constructed a dual-hybrid direct random phase approximation method, called dRPA75, and demonstrated its good performance on reaction energies, barrier heights, and non-covalent interactions of main-group elements. However, this method has also shown significant but quite systematic errors in the computed atomization energies. In this paper, we suggest a constrained spincomponent scaling formalism for the dRPA75 method (SCS-dRPA75) in order to overcome the large error in the computed atomization energies preserving the good performance of this method on spinunpolarized systems at the same time. The SCS-dRPA75 method with the aug-cc-pVTZ basis set results in an average error lower than 1.5 kcal mol-1 for the entire n-homodesmotic hierarchy of hydrocarbon reactions (RC0-5 test sets). The overall performance of this method is better than the related direct random phase approximation based double hybrid PWRB95 method on open-shell systems of maingroup elements (from the GMTKN30 database), and comparable to the best O(N4)-scaling opposite-spin second-order perturbation theory based double hybrid methods like PWPB95-D3, and to the O(N5)scaling RPAX2@PBEx method, which also includes exchange interactions. Furthermore, it gives a wellbalanced performance on many types of barrier heights similarly to the best O(N5)-scaling second-order perturbation theory based or spin-component scaled second-order perturbation theory based double hybrid methods such as XYG3 or DSD-PBEhB95. Finally, we show that the SCS-dRPA75 method has reduced self-interaction and delocalization errors compared to the parent dRPA75 method, and slightly smaller static correlation error than the related PWRB95 method. Introduction The direct random phase approximation (RPA or dRPA)1,2,3,4,5,6 is a promising non-local functional. The dRPA correlation seamlessly integrates dispersion interactions; therefore, it describes well (with a small underestimation of the long range correlation) the intra- and intermolecular dispersion interactions, the interaction energies of van der Waals bonded solids,7,8,9 as well as it might describe well the energetics of adsorption processes and interlayer interactions.10,11,12 However, it is well known that RPA overestimates the correlation energy at short range.13,14 This error leads to the difference between the RPA-parametrized, and the quantum Monte-Carlo (QMC) parametrized VWN or PW92 (local spin density approximation correlation functionals)15,16 correlation energies for the homogeneous electron gas. This overestimation is also responsible for the overestimation of the atomic energies, and thus the underestimation of the atomization energies of molecules and solids.17,18 Previously, we have demonstrated for hydrocarbons that this error is quite systematic and can be efficiently remedied by the corrected atomic energies technique.19 1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 2 of 17

Several approaches (such as density functional corrections, perturbative corrections, local field corrections, and optimized effective potential methods) were developed in order to diminish the atomization energy error and include beyond-RPA information.20 The semi-local RPA+ correction indeed improves the total atomic energies, but it works only for a well-localized exchange-correlation hole.21 Hence it worsens even the atomization energies, and as it has been recently shown, it also worsens the reaction energies for several classes of chemical reactions. The non-local RPA++ correction slightly improves the atomization energies, but it was concluded that a different kind of non-locality is required to overcome the problem.13 A more sophisticated version of RPA+ is the local-hybrid RPA+nlc correction, which performs somewhat better.14 The double hybrid (PBE+RPA+)/2 (mixture of the nonempirical PBE generalized gradient approximation functional and the RPA+ method) or the recent PWRB95 methods also provided possible solutions,13,22 although the double hybrid form needs some dispersion correction, which usually cannot capture the non-pairwise nature of the dispersion interaction between hollow molecules such as fullerenes.23 Recently, we have constructed the dRPA75 dual-hybrid direct random phase approximation method,19 and demonstrated its good performance on reaction energies, barrier heights, and non-covalent interactions of main-group elements. In this method, the orbitals are computed by the PBE75 global hybrid functional with 75% of exact exchange, and in the total energy, the PBE correlation energy is simply replaced by the non-self-consistent dRPA correlation energy. The high fraction of exact exchange both in the orbitals and in the total energy is beneficial for reducing the self-interaction error, which is important in the accurate description of reaction barrier heights and charge transfer interactions.24,25 The non-local correlation part captures the many-body nature of dispersion interactions at the level of double excitations and can be fine-tuned with an atom-pairwise correction to consider also some effect from the higher excitations as it was done by Brauer, Kesharwani, Kozuch, and Martin. 26 Despite the improved properties, the major problem with the dRPA75 method was its poor, RPA-like atomization energies. In this paper, we show that accurate atomization energies can be obtained by spin-component scaling in spite of the large static correlation error.22 We also compare our results to the dRPA-based double hybrid PWRB95 method developed by Grimme and Steinmetz,22 which uses pure semi-local reparametrized mPW91B95 orbitals,27,28,29,30 as well as 50% of exact exchange, 50% semi-local exchange, 35% of dRPA correlation, 65% of non-local van der Waals correction,31 and 71% of semi-local correlation in the total energy. Notice that the lower amount of exact exchange might lead to substantial self-interaction error, although this contribution was specifically chosen as a compromise between main-group and transition-metal thermochemistry similarly to the second-order perturbation theory-based (PT2-based) double hybrid PWPB95 method.32 Other fundamental problem with the PWRB95 method compared to the dRPA75 approach is that the semi-local correlation part also needs a non-local dispersion correction, which misses the many-body dispersion interaction terms. Furthermore, the spin-decomposition in the reparametrized B95 correlation functional (a large part of PWRB95) is based on the Stoll-Pavlidou-Preuss analysis,33 which is inaccurate for the spin-unpolarized homogeneous electron gas.34,35 Notice also that the perturbative second-order screened exchange (SOSEX) corrects only the same-spin component of the RPA correlation energy, thus improves the description of the short-range correlation and the ground-state energies upon RPA;36 however, this correction results in an underestimation of the long-range correlation energy, which is almost exact in RPA. Among the local field corrections, exchange-only kernels have been proven very practical since this approach allows an analytic evaluation of the coupling constant integral in the adiabatic-connection fluctuation-dissipation theorem. A specific group of exchange-only kernels is constructed to give the correct high-density limit of the spin-unpolarized uniform electron gas, although such corrections make somewhat larger discrepancy for lower densities.37,38,39 2 ACS Paragon Plus Environment

Page 3 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

In this paper, we propose the spin-component scaling (SCS)40,41,42 of the dRPA correlation energy in our dRPA75 scheme to improve its performance for open-shell systems. We compare our new SCSdRPA75 method to O(N4)-scaling methods, such as PWPB95-D3, and PWRB95, as well as to O(N5)scaling methods such as XYG3, DSD-PBEhB95, and RPAX2. The PWPB95-D3 method (50% exact exchange; 26.9% opposite-spin second-order perturbative correlation, D3 dispersion correction with zero-damping),32,43 the XYG3 method (80.33% exact exchange; 32.11% second-order perturbative correlation, global hybrid B3LYP reference orbitals),44,45 and the DSD-PBEhB95 method (66% exact exchange; 47% opposite-spin, and 9% same-spin second-order perturbative correlation; D3 dispersion correction with Becke-Johnson damping)46,47 are the best performing PT2-based double hybrid methods for the GMTKN30 database.48,49 The RPAX2@PBEx method is a random phase approximation method which includes exchange interactions (RPAX2), and the RPAX2 energy is evaluated on PBE exchange reference orbitals.50 Theory Previously, we defined the dual-hybrid dRPA75 method,19 where the dRPA75 exchangecorrelation energy was calculated from the exact exchange(𝐸𝑋𝑒𝑥𝑎𝑐𝑡 ), semi-local PBE exchange (𝐸𝑋𝑃𝐵𝐸 ), and dRPA correlation energy(𝐸𝐶𝑑𝑅𝑃𝐴 ) evaluated on self-consistent PBE hybrid orbitals with 0.75 fraction of the exact exchange (PBE75) as 𝑑𝑅𝑃𝐴75 (1) 𝐸𝑋𝐶 = (0.75𝐸𝑋𝑒𝑥𝑎𝑐𝑡 + 0.25𝐸𝑋𝑃𝐵𝐸 + 𝐸𝐶𝑑𝑅𝑃𝐴 )@𝑃𝐵𝐸75. The dRPA correlation energy is obtained in the coupled cluster formalism from the B nonantisymmetrized two-electron integral matrix, and the T double excitation amplitude matrix20 as 1 (2) 𝐸𝐶𝑑𝑅𝑃𝐴 = 2 𝑡𝑟[𝐵𝑇], where the B matrix elements are defined by 𝐵𝑖𝑎,𝑗𝑏 = ⟨𝑖𝑗|𝑎𝑏⟩ for occupied orbitals i and j, and virtual orbitals a and b, and the matrix T with matrix elements 𝑇𝑖𝑎,𝑗𝑏 = 𝑇𝑖𝑗𝑎𝑏 is computed in an iterative update procedure initialized with T(0) = 0 as (3) 𝑇 (𝑛+1) = −∆ ∘ (𝐵 + 𝐵𝑇 (𝑛) + 𝑇 (𝑛) 𝐵 + 𝑇 (𝑛) 𝐵𝑇 (𝑛) ). The circle sign (∘ ) denotes the Hadamard product of matrices. The ∆orbital energy denominator matrix −1 elements are ∆𝑖𝑎,𝑗𝑏 = (𝜀𝑎 + 𝜀𝑏 − 𝜀𝑖 − 𝜀𝑗 ) with 𝜀𝑝 as the orbital energy of orbital p. The equations are solved in an O(N4)-scaling iterative procedure using density-fitting for the two-electron repulsion integrals, and Cholesky-decomposition for the orbital energy denominator matrix.51 We can rewrite the dRPA correlation energy as a sum of the opposite-spin 𝐸𝐶𝑑𝑅𝑃𝐴𝑜𝑠 , and samespin 𝐸𝐶𝑑𝑅𝑃𝐴𝑠𝑠 correlation energies as 𝑎𝛼 𝑏 (4) 𝐸 𝑑𝑅𝑃𝐴𝑜𝑠 = ∑ ∑ ⟨𝑖 𝑗 |𝑎 𝑏 ⟩ 𝑇 𝛽 𝛼 𝛽

𝐶

𝛼 𝛽

𝑖𝛼 𝑗𝛽

𝑖𝛼 𝑗𝛽 𝑎𝛼 𝑏𝛽

𝐸𝐶𝑑𝑅𝑃𝐴𝑠𝑠 =

1 1 𝑎 𝑏 𝑎 𝑏 ∑ ∑ ⟨𝑖𝛼 𝑗𝛼 |𝑎𝛼 𝑏𝛼 ⟩ 𝑇𝑖𝛼𝛼𝑗𝛼 𝛼 + ∑ ∑ ⟨𝑖𝛽 𝑗𝛽 |𝑎𝛽 𝑏𝛽 ⟩ 𝑇𝑖 𝛽𝑗 𝛽 𝛽 𝛽 2 2 𝑖𝛼 𝑗𝛼 𝑎𝛼 𝑏𝛼

(5) ,

𝑖𝛽 𝑗𝛽 𝑎𝛽 𝑏𝛽

where the subscript of the orbital index refers to its spin. At the derivation of eq (4), we utilized that matrix T is symmetric, which comes from eq (3) and the fact that matrix B is symmetric, and 𝑏 𝑎𝛼 𝑎𝛼 𝑏 consequently 𝑇𝑗 𝛽𝑖 = 𝑇𝑖 𝑗 𝛽 . For any restricted closed-shell calculation, the ⟨𝑖𝜎 𝑗𝜎′ |𝑎𝜎 𝑏𝜎′ ⟩ integrals are 𝛽 𝛼

𝛼 𝛽

3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

𝑎 𝑏

𝑎𝛽 𝑏𝛽 𝛽 𝑗𝛽

the same for each σσ' spin pair. The 𝑇𝑖𝛼𝛼𝑗𝛼 𝛼 and 𝑇𝑖

Page 4 of 17

amplitudes are trivially equal. Using the recursion 𝑎 𝑏

𝑎𝛼 𝑏𝛽

formula for the T matrix, it can be shown that also 𝑇𝑖𝛼𝛼𝑗𝛼 𝛼 = 𝑇𝑖𝛼 𝑗

𝛽

since the ∆𝑖𝜎 𝑎𝜎 ,𝑗𝜎′ 𝑏𝜎′ denominators

are the same for each σσ' spin pair. Therefore, in the restricted closed-shell case, 𝐸𝐶𝑑𝑅𝑃𝐴𝑜𝑠 and 𝐸𝐶𝑑𝑅𝑃𝐴𝑠𝑠 are equal, and their sum is 2𝐸𝐶𝑑𝑅𝑃𝐴𝑜𝑠 . In the dRPA correlation energy as well as in the dRPA75 correlation energy, the exchange contribution to the same-spin correlation energy is missing, and consequently the same-spin correlation energy is too large in magnitude. To remedy this problem, it seems to be plausible to scale down this contribution. Thus, in this paper, we introduce a general spin-component scaled formula for the non-selfconsistent dRPA methods (SCS-dRPA@DFT), 𝑆𝐶𝑆−𝑑𝑅𝑃𝐴 (6) 𝐸𝑋𝐶 = (𝐸𝑋𝑒𝑥𝑎𝑐𝑡 + 𝑎𝑜𝑠 𝐸𝐶𝑑𝑅𝑃𝐴𝑜𝑠 + (2 − 𝑎𝑜𝑠 )𝐸𝐶𝑑𝑅𝑃𝐴𝑠𝑠 )@𝐷𝐹𝑇, and an analogous spin-component scaled formula for the dRPA75 approach (SCS-dRPA75), 𝑆𝐶𝑆−𝑑𝑅𝑃𝐴75 𝐸𝑋𝐶 = (0.75𝐸𝑋𝑒𝑥𝑎𝑐𝑡 + 0.25𝐸𝑋𝑃𝐵𝐸 + 𝑎𝑜𝑠 𝐸𝐶𝑑𝑅𝑃𝐴𝑜𝑠 + (2 − 𝑎𝑜𝑠 )𝐸𝐶𝑑𝑅𝑃𝐴𝑠𝑠 )@𝑃𝐵𝐸75,

(7)

which keep the original functionals intact for closed-shell systems with restricted orbitals. Although these energy components might be scaled independently, scaling-up the opposite-spin correlation energy leads to more accurate atomization energies than keeping it unchanged according to our preliminary results. Another purpose for the symmetrical scaling, besides it guarantees the unchanged results for closed-shell systems with restricted orbitals, that 𝑐𝑜𝑠 + 𝑐𝑠𝑠 = 2 should hold at large interelectronic distances for a general purpose SCS approach.41 Since the major problem that prevents the dRPA75 method from being a robust method for main-group thermochemistry is its poor atomization energies, we use atomization energy test sets (RC0 test set containing the atomization energies of 38 hydrocarbons with CCSD(T)/CBS reference energies; AE6 test set constituted of six difficult cases with CCSDT(Q)-full/CBS reference energies)52,53,54 for the optimization of the 𝑎𝑜𝑠 parameter, but we also give theoretical arguments for the determination of its value based on the fully spin-polarized homogeneous electron gas model. Similarly to our previous work, we combine our SCS-dRPA75 method with the aug-cc-pVTZ basis set for systematic error compensation. We define 𝑎𝑜𝑠 = 1.5 in eq (7) for the SCS-dRPA75/aug-cc-pVTZ method (vide infra the parameter determination). We note that the non-self-consistent treatment makes the formulation of the gradient and higher-order derivatives somewhat complicated, but the analytic gradients can be derived very similarly to the parent dRPA method.55 Computational details For the computations, the MRCC quantum chemistry software56,57 was applied (except for the PWPB95 computations, where we used the Orca quantum chemistry package),58 and the frozen core and density fitting approximations were invoked. For the non-self-consistent dRPA computations, we used different reference determinants obtained with the PBE, SCAN (strongly constrained and appropriately normed meta-GGA functional), and PBE0 (hybrid version of PBE with 25% of exact exchange) methods.59,60,61 For the dRPA computations, we applied the aug-cc-pVTZ basis set62,63,64,65 with the augcc-pVTZ-RI-JK,66 and aug-cc-pVTZ-RI67,68 density fitting basis sets for the SCF and dRPA correlation calculations, respectively if they are available for the element. Otherwise, we applied the def2-QZVPPRI-JK,69,70 and def2-QZVPP-RI68,70 density fitting basis sets. For the double hybrid XYG3 and DSD-PBEhB95 computations, we used the sufficiently large (aug)-def2-QZVPP basis set (def2-QZVP basis set augmented with aug-cc-pVQZ diffuse functions on the heavy atoms, and extra d and f polarization functions on the Na and Mg atoms)70,71,72,73,74 with the def2-QZVPP-RI-JK, and def2-QZVPP-RI auxiliary basis sets. We also performed double hybrid

4 ACS Paragon Plus Environment

Page 5 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

calculations with the smaller aug-cc-pVTZ basis set so that we could gain some experience what accuracy can be achieved with these methods at costs similar to our dRPA75/aug-cc-pVTZ approach. We have examined the basis set dependence of the RPAX2@PBEx energies with the aug-ccpVTZ, aug-cc-pVQZ, and aug-cc-pV5Z basis sets.62,63,64,65 For these basis sets, we applied the corresponding RI-JK and RI auxiliary basis sets if they are available for the element, otherwise the def2QZVPP-RI-JK and def2-QZVPP-RI density fitting basis sets. (For the Li, Be, Na, and Mg atoms, we used the cc-pV5Z basis set because the aug-cc-pV5Z basis set was unavailable.) In the following, we investigate the performance of the new SCS-dRPA75 method also on several open-shell test sets such as atomization energies (W4-08 containing 16 molecules with and 83 molecules without multireference character, with W4 reference energies), 75 self-interaction error related problems (SIE11 with CCSD(T)/CBS reference energies),76 ionization potentials (G21IP with back-corrected experimental reference values),77 electron affinities (G21EA with back-corrected experimental reference values),77 barrier heights (BH76 with W1 and theoretically estimated reference values; DBH24 with CCSDT(Q)-full/CBS reference energies for the BH6 subset, with W3.2 or W4 reference energies elsewhere),75,78,79,54 and radical stabilization energies (RSE43 with CCSD(T)/CBS reference energies).80 Furthermore, we provide the hydrogen molecular ion (H2+) and hydrogen (H2) dissociation curves to compare the delocalization and static correlation errors with those of other related methods.81,82 Parameter determination Homogeneous electron gas Similarly to the non-self-consistent dRPA@DFT methods, the reference energy and the exchange part of the dRPA75 total energy are exact for the homogeneous electron gas by construction. So all the errors of dRPA75 for the homogeneous electron gas come from the dRPA correlation energy. In this paper, we would like to correct the erroneous dRPA75/aug-cc-pVTZ atomization energies conserving the good performance of this method for reaction energies, barrier heights, and non-covalent interactions of main-group elements. Therefore, we require the total energy of spin-unpolarized systems not to change with the spin-component scaling (eq (7)). Here we investigate the effect of spin-scaling on the error of the RPA correlation energy on fully spin-polarized homogeneous electron gas with an 𝑟𝑆 = [3⁄(4𝜋𝑛)]1⁄3 Wigner-Seitz radius, where n is the electron density, and the spin-polarized correlation energy per particle 𝜀𝑐𝑝𝑜𝑙 can be well estimated within or beyond RPA using a simple expression by Perdew and Wang16 as ⁄ ⁄ (8) 𝜀𝑐𝑝𝑜𝑙 (𝑟𝑆 ) ≈ −2𝐴(1 + 𝛼1 𝑟𝑆 )𝑙𝑛(1 + 1⁄2𝐴(𝛽1 𝑟𝑆1 2 + 𝛽2 𝑟𝑆 + 𝛽3 𝑟𝑆3 2 + 𝛽4 𝑟𝑆𝑝+1 )). The parameter values for RPA (in a.u.) are: 𝐴 = 0.015545, 𝛼1 = 0.035374, 𝛽1 = 6.4869, 𝛽2 = 1.3083, 𝛽3 = 0.15180, 𝛽4 = 0.082349, and 𝑝 = 0.75. Whereas the QMC parameters are: 𝐴 = 0.015545, 𝛼1 = 0.20548, 𝛽1 = 14.1189, 𝛽2 = 6.1977, 𝛽3 = 3.3662, 𝛽4 = 0.62517, and 𝑝 = 1.00. The dRPA error for any fully spin-polarized density comes completely from the same-spin correlation energy. Figure 1 shows that this error can be minimized in a wide range of densities by scaling down the same-spin correlation energy approximately to 50%. Note that scaling down the samespin correlation energy to approximately 50% means 𝑎𝑜𝑠 ≈ 1.5 in the SCS-dRPA75 formula (eq (7)).

5 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

0.05

Error in the correlation energy per electron (a.u.)

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 6 of 17

0.04 0.03 0.02

0

0.01

0.25 0.5

0.00

0.75 -0.01

1

-0.02 -0.03 0

5

10

15

20

25

30

35

40

45

50

Wigner-Seitz radius (a.u.) Figure 1. Effect of scaling-down the same-spin correlation energy on the error of the RPA correlation energy (using the RPA parametrization of the PW92 correlation functional) for fully spin-polarized homogeneous electron gas. The reference method is the QMC parametrization of the PW92 correlation functional. The same-spin scaling factors are shown in the key. Atomization energies Previously, we have shown that the dRPA75/aug-cc-pVTZ method has 75.04 kcal mol-1 error on average for the hydrocarbon atomization energies of the RC0 test set.52 This serious RPA-like underestimation of the atomization energies is due to the overestimation of the short-range correlation in the atomic energies, which is hopefully corrected by the spin-scaling. In fact, the spin-component scaling has a quite dramatic effect on the computed dRPA75/aug-cc-pVTZ atomization energies for the RC0 test set (Figure 2). The endothermic atomization energy error linearly shifts towards the exothermic direction. According to the mean deviation (MD), mean absolute deviation (MAD), and standard deviation (SD), the optimal 𝑎𝑜𝑠 parameter in the SCS-dRPA75 formula (eq (7)) is 1.50 for this test set (resulting in MD = 0.13 kcal mol-1, MAD = 0.88 kcal mol-1, and SD = 1.09 kcal mol-1), which agrees well with our observation for the fully spin-polarized homogeneous electron gas. The more accurate atomization energies make the performance of the SCS-dRPA75/aug-cc-pVTZ method more balanced on the n-homodesmotic hierarchy of hydrocarbon reaction energies (MAD = 0.88, 1.45, 1.16, 0.42, 0.27, 0.19 kcal mol-1 for the RC0-RC5 test sets, in order) than the performance of the related dRPA-based double hybrid PWRB95/def2-QZVP method (MAD = 6.35, 1.98, 1.47, 0.33, 0.41, 0.33 kcal mol-1 for the RC0-RC5 test sets, in order).52,19,22

6 ACS Paragon Plus Environment

Page 7 of 17

80 60

Deviations (kcal mol-1)

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

40 20 MD

0 1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

-20

1.8

1.9

2

MAD SD

-40 -60 -80

Scaling factor

Figure 2. Mean (MD), mean absolute (MAD), and standard deviations (SD) from the benchmark atomization energies of the RC0 test set as a function of the scaling factor of the opposite-spin correlation energy in the SCS-dRPA75 method using the aug-cc-pVTZ basis set Figure 3 shows that 𝑎𝑜𝑠 = 1.51 is optimal (with MAD = 3.07 kcal mol-1) in eq (7) using the augcc-pVTZ basis set on the small but representative AE6 test set, which agrees well with the 𝑎𝑜𝑠 = 1.5 parameter value obtained for the RC0 test set (using 𝑎𝑜𝑠 = 1.50 the corresponding MAD is 3.18 kcal mol-1 for the AE6 test set). At this point, we have also tried to optimize the scaling factors of the oppositespin and same-spin correlation energy components separately (using 𝑎𝑜𝑠 and a linearly independent 𝑎𝑠𝑠 instead of 𝑎𝑠𝑠 = 2 − 𝑎𝑜𝑠 ) in eq (7) using the aug-cc-pVTZ basis set. Minimizing the root mean square deviation from the AE6 reference atomization energies, the obtained scaling factors are 𝑎𝑜𝑠 = 1.48 and 𝑎𝑠𝑠 = 0.58 (with MAD = 2.98 kcal mol-1). Therefore, 𝑎𝑜𝑠 + 𝑎𝑠𝑠 = 2.06 ≈ 2, which supports our previous assumption. For comparison, Grimme's SCS-MP2 scaling factors are somewhat smaller (𝑐𝑜𝑠 = 1.20 and 𝑐𝑠𝑠 = 0.33). (For other spin-component scaling approaches, the interested reader is referred to the works of Fink and Szabados.83,84) In contrast to the dRPA75 method, where the complete basis set extrapolated atomization energies are closer to the reference by 10.5 kcal mol -1 on average than those with the aug-cc-pVTZ basis set, the complete basis set extrapolation with the proposed parametrization worsens the calculated SCS-dRPA75 atomization energies (MAD = 14.9 kcal mol-1 using 𝑎𝑜𝑠 = 1.5 in eq (7)) since the optimal parametrization is different at the basis set limit (MAD = 2.4 kcal mol -1 using 𝑎𝑜𝑠 = 1.3 in eq (7)). (Here we follow our previous work:19 we extrapolate the hybrid PBE75 exchange and dRPA75 or SCS-dRPA75 correlation parts separately to the complete basis set limit with a two-point extrapolation scheme using aug-cc-pVTZ and aug-cc-pVQZ basis sets, and with the following extrapolation coefficients: 𝐶3𝑃𝐵𝐸75𝑥 = 0.312 and 𝐶3𝑑𝑅𝑃𝐴75𝑐 = 𝐶3𝑆𝐶𝑆−𝑑𝑅𝑃𝐴75𝑐 = 0.837.) ⁄4 ⁄4 ⁄4 The optimal 𝑎𝑜𝑠 parameter values are about 1.43, 1.45, and 1.44 for the SCS-dRPA@PBE/augcc-pVTZ, SCS-dRPA@PBE0/aug-cc-pVTZ, and SCS-dRPA@SCAN/aug-cc-pVTZ methods, respectively. However, the corresponding mean absolute deviations from the reference are considerably larger (~6 kcal mol-1) for the pure dRPA methods. It is also noticeable that the dependence of the dRPA 7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

errors on the reference determinant (GGA > meta-GGA > hyper-GGA) almost vanishes with the spincomponent scaling around the optimum similarly as it was observed earlier for PT2-based double hybrids.85 Note that the optimal 𝑎𝑜𝑠 parameter for the complete basis set extrapolated SCS-dRPA@PBE, SCS-dRPA@PBE0, and SCS-dRPA@SCAN methods would be about 1.26 with a mean absolute deviation of 3.4-3.2 kcal mol-1. (Here we follow our previous work:19 we extrapolate the exact exchange and dRPA correlation parts separately to the complete basis set limit with a two-point extrapolation scheme using aug-cc-pVTZ and aug-cc-pVQZ basis sets, and with the following extrapolation 𝑑𝑅𝑃𝐴𝑐@𝑃𝐵𝐸 coefficients: 𝐶3𝐸𝑋𝑋 = 𝐶3𝑑𝑅𝑃𝐴𝑐@𝑆𝐶𝐴𝑁 = 0.856, and 𝐶3𝑑𝑅𝑃𝐴𝑐@𝑃𝐵𝐸0 = 0.867.) ⁄4 = 0.274, 𝐶3⁄4 ⁄4 ⁄4

Mean absolute deviation (kcal mol-1)

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 8 of 17

SCS-dRPA@PBE/aug-cc-pVTZ

SCS-dRPA@SCAN/aug-cc-pVTZ

SCS-dRPA@PBE0/aug-cc-pVTZ

SCS-dRPA75/aug-cc-pVTZ

60 50 40 30 20 10 0 1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

Scaling factor Figure 3. Mean absolute deviations from the benchmark atomization energies of the AE6 test set as a function of the scaling factor of the opposite-spin correlation energy in the SCS-dRPA, and SCS-dRPA75 methods using the aug-cc-pVTZ basis set We observed for the AE6 test set that the MD = -2.85 kcal mol-1 and MAD = 3.14 kcal mol-1 values of the PWRB95/def2-QZVP method published in ref 22 disagree with the MD = -4.35 kcal mol-1 and MAD = 4.48 kcal mol-1 calculated from the results published in the supporting information of the same paper. Using this latter value, the MAD is 5.4 kcal mol-1 (instead of 4.7 kcal mol-1 as published in ref 22) for the AE6 and RC0 test sets. The SCS-dRPA75/aug-cc-pVTZ method yields a considerably smaller MAD for the same combined test set (2.0 kcal mol-1).

8 ACS Paragon Plus Environment

Page 9 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Results and discussion Open-shell systems For the larger W4-08 atomization energy test set, the SCS-dRPA75/aug-cc-pVTZ method (MAD = 4.8 kcal mol-1) shows moderately accurate results (Figure 4), similarly to the dRPA-based double hybrid PWRB95/(aug)-def2-QZVP method (MAD = 3.6 kcal mol-1). The slightly lower accuracy of the SCS-dRPA75/aug-cc-pVTZ results can be explained by the large basis set error since with the similarsized def2-TZVP(D) basis set the PWRB95 method gives larger errors (MAD = 5.0 kcal mol-1) as well.22 Note that the second-order perturbation theory-based double hybrid methods such as XYG3/(aug)-def2QZVP (MAD = 2.36 kcal mol-1), DSD-PBEhB95/(aug)-def2-QZVPP (MAD = 2.01 kcal mol-1), or PWPB95-D3/def2-QZVP (MAD = 1.89 kcal mol-1) perform relatively well, although large basis set dependence can be observed for the XYG3 and DSD-PBEhB95 atomization energies similarly to the XYG3 heats of formation in the literature.86 Notice that the PWPB95-D3 atomization energies have only a small basis set dependence since this is the only double hybrid method reported in the literature that can be used with a triple-zeta quality basis set without the need for reparametrization.49,87 The RPAX2@PBEx atomization energies also improve considerably with the increasing basis set size (MAD = 8.82, 4.02, 2.68 kcal mol-1 for the aug-cc-pVTZ, QZ, and 5Z basis sets, respectively). Figure 4 shows that the spin-component scaling somewhat also affects the self-interaction error related problems (SIE11), ionization potentials (G21IP), and electron affinities (G21EA); however, it has only a very small effect on forward/backward barrier heights (BH76), on the corresponding reaction energies (BH76RC: MAD = 1.30 kcal mol-1 for dRPA75/aug-cc-pVTZ, MAD = 1.42 kcal mol-1 for SCSdRPA75/aug-cc-pVTZ), and on the radical stabilization energies (RSE43). For these test sets, among the O(N4)-scaling methods, the SCS-dRPA75/aug-cc-pVTZ method performs better than the PWRB95/(aug)-def2-QZVP method, and it shows similar performance to the double hybrid PWPB95D3/def2-QZVP method. Note that the O(N5)-scaling methods such as DSD-PBEhB95, XYG3, and RPAX2 perform better for these test sets; however, they require larger basis sets. The DSDPBEhB95/(aug)-def2-QZVPP results are better for barrier heights and self-interaction related problems than the XYG3/(aug)-def2-QZVP results, but they are worse for ionization potentials, electron affinities, and radical stabilization energies.

9 ACS Paragon Plus Environment

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Mean absolute deviation (kcal mol-1)

Journal of Chemical Theory and Computation

9

Page 10 of 17

15.7

8 7 6 5 4 3 2 1

0 W4-08

SIE11

G21IP

G21EA

BH76

dRPA75/aug-cc-pVTZ

SCS-dRPA75/aug-cc-pVTZ

PWRB95/(aug)-def2-QZVP

PWPB95-D3/def2-QZVP

DSD-PBEhB95-D3(BJ)/(aug)-def2-QZVPP

XYG3/(aug)-def2-QZVPP

RSE43

RPAX2@PBEx/aug-cc-pVQZ

Figure 4. Performance of the SCS-dRPA75/aug-cc-pVTZ method compared to other approaches on several test sets containing open-shell systems (W4-08: atomization energies; SIE11: self-interaction error related problems; G21IP: ionization potentials; G21EA: electron affinities; BH76: barrier heights; RSE43: radical stabilization energies). The error bars denote the PWRB95/def2-TZVP(D), PWPB95D3/aug-cc-pVTZ, DSD-PBEhB95/aug-cc-pVTZ, XYG3/aug-cc-pVTZ, and RPAX2@PBEx/aug-cc-pVTZ results. (The PWRB95/def2-TZVP(D) and PWRB95/def2-QZVP results are taken from ref 22. The PWPB95-D3/def2-QZVP results are taken from ref 32.) In the BH76 test set, the hydrogen atom transfer barrier heights are overrepresented, and the reference energies are not sufficiently accurate to distinguish the performance of the dRPA75 and SCSdRPA75 methods. For a better comparison, we have selected the more balanced DBH24 barrier height test set, which is composed of four subsets (HAT6: six heavy atom transfer barriers; NS6: six nucleophilic substitution reaction barriers; BH6: six hydrogen atom transfer barriers; UR6: six unimolecular and recombination reaction barriers), and the reference energies are also more accurate than those of the BH76 test set. Figure 5 shows that the largest effect of the spin-component scaling is on the heavy atom transfer barrier heights, and there is no effect on the nucleophilic substitution barrier heights since each component has closed shell electronic structure in the reactions. The SCSdRPA75/aug-cc-pVTZ approach performs significantly better than the PWRB95/def2-QZVP method on hydrogen transfer, heavy atom transfer, and nucleophilic substitution barrier heights, and only slightly worse on the barrier heights of unimolecular and recombination reactions, where the spin-component scaling worsens the dRPA75/aug-cc-pVTZ results. (Note that the lack of diffuse functions in the def2QZVP basis set partially explains the larger errors of the PWRB95/def2-QZVP method, especially for the NS6 test set, which contains negatively charged species. Nevertheless the PWRB95/(aug)-def2QZVP results for the BH76 test set suggest that the PWRB95 method is generally worse for barrier heights than the SCS-dRPA75/aug-cc-pVTZ method.) Except for the unimolecular and recombination reaction barrier heights, the SCS-dRPA75/aug-cc-pVTZ results are more accurate than the PWPB95D3/def2-QZVP results. The SCS-dRPA75/aug-cc-pVTZ barrier heights are generally only a bit less 10 ACS Paragon Plus Environment

Page 11 of 17

accurate than the XYG3/(aug)-def2-QZVP barrier heights. Furthermore, the RPAX2@PBEx/aug-ccpVQZ method shows a less balanced performance for barrier heights than the SCS-dRPA75/aug-ccpVTZ scheme. Finally, the DSD-PBEhB95/(aug)-def2-QZVPP barrier heights are generally more accurate than the SCS-dRPA75/aug-cc-pVTZ ones, although the lack of diffuse functions in the def2QZVPP basis set can deteriorate this good performance. (The basis set dependence of the dRPA75 and SCS-dRPA75 barrier heights is illustrated on the BH6 test set, and the results are presented in the supporting information.)

Mean absolute deviation (kcal mol-1)

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 HAT6

NS6

BH6

dRPA75/aug-cc-pVTZ

SCS-dRPA75/aug-cc-pVTZ

PWRB95/def2-QZVP

PWPB95-D3/def2-QZVP

DSD-PBEhB95-D3(BJ)/(aug)-def2-QZVPP

XYG3/(aug)-def2-QZVPP

UR6

RPAX2@PBEx/aug-cc-pVQZ

Figure 5. Performance of the SCS-dRPA75/aug-cc-pVTZ method compared to the other approaches on the different barrier heights of the DBH24 test set (HAT6: heavy atom transfer; NS6: nucleophilic substitution; BH6: hydrogen atom transfer; UR6: unimolecular and recombination reactions). The error bars denote the DSD-PBEhB95/def2-QZVP, and XYG3/def2-QZVP results. (The PWRB95/def2-QZVP results are taken from ref 22, the PWPB95-D3/def2-QZVP, and XYG3/def2-QZVP results are taken from ref 32.)

11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 12 of 17

Self-interaction, delocalization, and static correlation errors In the following, we investigate the self-interation, delocalization, and static correlation errors of the SCS-dRPA75 method. The error in the energy of the hydrogen atom gives a very simple indication of the self-interaction error. Furthermore, we can compare the magnitude of the delocalization error of the density functional methods on the energy of the stretched H2+ molecule, and the magnitude of the static correlation error on the energy of the stretched H2 molecule. (Note that the basis set error is expected to be small for hydrogen-atom-containing systems, so we can use the aug-cc-pVTZ and def2QZVP basis sets henceforward.) The dRPA@HF method (with the aug-cc-pVTZ basis set) shows 8.32 kcal mol-1 self-interaction error on the hydrogen atom, which comes completely from the dRPA correlation. The dRPA75 method (with the aug-cc-pVTZ basis set) suffers almost from the same amount of self-interaction error (-7.96 kcal mol-1). Scaling down the same-spin correlation energy, the SCS-dRPA75 method halves this error (-3.41 kcal mol-1). For comparison, the PWRB95 method (with the def2-QZVP basis set) shows larger (-4.89 kcal mol-1) self-interaction error on the hydrogen atom. In the asymptotic region of the H2+ binding energy curve (Figure 6), the spin-component scaling has larger effect than near the equilibrium bond length. The SCS-dRPA75 binding energy curve shows a delocalization error similar to the dRPA@HF binding energy curve. The PWRB95 binding energy curve runs below the dRPA75 curve, which supports our previous observation for the SIE11 test set that the PWRB95 method is loaded by larger delocalization error. The spin-component scaling has no effect on the static correlation error since the hydrogen molecule has a closed shell electronic structure. Although it affects the energy of the hydrogen atom in the H2 binding energy, which results in a constant exothermic shift in the binding energy curve (Figure 7). In conclusion, the self-interaction error was responsible for the deviation of the dRPA75 binding energy curve from the reference near the equilibrium distance. (The other binding energy curves run close to the reference near the equilibrium distance.) Furthermore, the dRPA75 and SCS-dRPA75 methods are loaded by slightly smaller static correlation error than the PWRB95 method, although they show still larger static correlation error than the dRPA@PBE method.

12 ACS Paragon Plus Environment

Page 13 of 17

Binding energy (kcal mol-1)

0 0

1

2

3

4

5

6

-10 -20 HF/CBS

-30

SCS-dRPA75/aug-cc-pVTZ

-40

dRPA@HF/aug-cc-pVTZ

-50

dRPA75/aug-cc-pVTZ PWRB95/def2-QZVP

-60 -70 -80

Bond length (Å)

Figure 6. Binding energy curves with different methods for the stretched positively charged hydrogen molecular ion (The reference method is the HF/CBS.)

60

Binding energy (kcal mol-1)

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

40 20 0 0

1

2

3

4

5

6

dRPA75/aug-cc-pVTZ

-20

PWRB95/def2-QZVP

-40

SCS-dRPA75/aug-cc-pVTZ dRPA@PBE/aug-cc-pVTZ

-60

CISD/CBS -80

-100 -120

Bond length (Å)

Figure 7. Binding energy curves with different methods for the stretched hydrogen molecule (The reference method is the CISD/CBS.)

13 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 14 of 17

Conclusions For the fully spin-polarized homogeneous electron gas, the (same-spin) random phase approximation correlation energy is roughly twice as large as the exact one. Due to the lack of exchange terms in the correlation energy, similar trend can be observed for molecular systems. To overcome this problem, we have proposed a constrained spin-component scaling scheme for our recent dRPA-based dual-hybrid dRPA75 method, which preserves its good performance for spin-unpolarized systems. Similarly to our previous work, we take advantage of the systematic error of the aug-cc-pVTZ basis set in contrast to the XYG3 and DSD-PBEhB95 double hybrid density functionals, which require larger basis sets. We have demonstrated that the RPA-like atomization energy error in the dRPA75/aug-ccpVTZ method can be significantly reduced by scaling-down the same-spin correlation energy to 50% and simultaneously scaling-up the opposite-spin correlation energy to 150%. Using this new method, called SCS-dRPA75/aug-cc-pVTZ, the entire n-homodesmotic hierarchy of hydrocarbon reactions (RC0-5 test sets) can be solved with an average error lower than 1.5 kcal mol-1. The spin-component scaling has similar effect also on the atomization energies of other molecules constituted of main-group elements (AE6 and W4-08 test sets). We have tested the SCS-dRPA75/aug-cc-pVTZ method also on other open-shell systems from the GMTKN30 database (SIE11, G21IP, G21EA, BH76, and RSE43 test sets) excluding atomization. These benchmarks test self-interaction error, ionization potentials, electron affinities, barrier heights, and radical stabilization energies. For these test sets, the spin-component scaling has a moderate effect. Among the O(N4)-scaling dRPA-based methods, the SCS-dRPA75/aug-cc-pVTZ method shows an overall better performance than the related double hybrid PWRB95 method. Its overall accuracy is comparable to the most successful O(N4)-scaling second-order perturbation theory-based double hybrid PWPB95-D3 method. Furthermore, similarly to the best O(N5)-scaling double hybrid density functionals, such as XYG3 and DSD-PBEhB95, it has a more balanced performance for many kinds of barrier heights (i.e., heavy atom transfer, nucleophilic substitution, hydrogen atom transfer, unimolecular, and recombination reactions) than the other O(N4)-scaling methods, such as PWRB95 and PWPB95-D3. Finally, we have estimated the magnitude of the self-interaction, delocalization, and static correlation errors of the SCS-dRPA75 method using the one-electron hydrogen atom, the potential energy curve of the positively charged hydrogen molecular ion, and the potential energy curve of the neutral hydrogen molecule. The spin-component scaling has a large effect on the energy of the openshell hydrogen atom, and the hydrogen molecular cation; however, no effect on the energy of the closedshell hydrogen molecule. The SCS-dRPA75 method shows considerably smaller self-interaction and delocalization errors than the parent dRPA75 method, and the related dRPA-based PWRB95 method. Furthermore, the new approach also exhibits a somewhat smaller static correlation error than the PWRB95 method, although it is still larger than the static correlation error in the dRPA@PBE method. Acknowledgments This paper is dedicated to Professor József Nagy on the occasion of his 90th birthday. The research work has been accomplished in the framework of the “BME R+D+I project,” supported by Grant No. TÁMOP 4.2.1/B-09/1/KMR-2010-0002. A.R. acknowledges the support of the National Science Foundation under Grant No. DMR-1553022. The computing time granted on the Hungarian HPC Infrastructure at NIIF Institute, Hungary, is gratefully acknowledged. Supporting information Further information about the parameter optimization, and detailed results with statistics for the discussed test sets. This information is available free of charge via the Internet at http://pubs.acs.org/. 14 ACS Paragon Plus Environment

Page 15 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1 Langreth, D.; Perdew, J. P. Phys. Rev. B 1977, 15, 2884–2901. 2 Langreth, D.; Perdew, J. P. Phys. Rev. B 1980, 21, 5469–5493. 3 Yan, Z.; Perdew, J. P.; Kurth, S. Phys. Rev. B 2000, 61, 16430-16439. 4 Furche, F. J. Chem. Phys. 2008, 129, 114105. 5 Kubo, R. Rep. Prog. Phys. 2002, 29, 255–284. 6 Scuseria, G. E.; Henderson, T. M.; Sorensen, D. C. J. Chem. Phys. 2008, 129, 231101. 7 Eshuis, H.; Furche, F. J. Phys. Chem. Lett. 2011, 2, 983–989. 8 Zhu, W.; Toulouse, J.; Savin, A.; Angyán, J. G. J. Chem. Phys. 2010, 132, 244108. 9 Harl, J.; Kresse, G. Phys. Rev. B 2008, 77, 1–8. 10 Ren, X.; Rinke, P.; Scheffler, M. Phys. Rev. B 2009, 80, 1–8. 11 Schimka, L.; Harl, J.; Stroppa, A.; Grüneis, A.; Marsman, M.; Mittendorfer, F.; Kresse, G. Nat. Mater. 2010, 9, 741– 744. 12 Lebègue, S.; Harl, J.; Gould, T.; Ángyán, J. G.; Kresse, G.; Dobson, J. F. Phys. Rev. Lett. 2010, 105, 196401. 13 Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I. J. Chem. Theory Comput. 2010, 6, 127–134. 14 Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I. J. Chem. Phys. 2011, 134, 114110. 15 Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200–1211. 16 Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45, 13244–13249. 17 Furche, F. Phys. Rev. B 2001, 64, 195120. 18 Harl, J.; Schimka, L.; Kresse, G. Phys. Rev. B 2010, 81, 115126. 19 Mezei, P. D.; Csonka, G. I.; Ruzsinszky, A.; Kállay, M. J. Chem. Theory Comput. 2015, 11, 4615-4626. 20 Eshuis, H.; Bates, J. E.; Furche, F. Theor. Chem. Acc. 2012, 131, 1084. 21 Yan, Z.; Perdew, J. P.; Kurth, S. Phys. Rev. B 2000, 61, 16430-16439. 22 Grimme, S.; Steinmetz, M. Phys. Chem. Chem. Phys. 2016, 18, 20926-20937. 23 Ruzsinszky, A.; Perdew, J. P.; Tao, J.; Csonka, G. I. Pitarke, J. M. Phys. Rev. Lett. 2012, 109, 233203. 24 Sai, N.; Barbara, P. F.; Leung, K. Phys. Rev. Lett. 2011, 106, 226403. 25 Atalla, V.; Yoon, M.; Caruso, F.; Rinke, P.; Scheffler, M. Phys. Rev. B 2013, 88, 165122. 26 Brauer, B.; Kesharwani, M. K.; Kozuch, S.; Martin, J. M. L. Phys. Chem. Chem. Phys. 2016, 18, 20905-20925. 27 Perdew, J. P. In Proceedings of the 21st Annual International Symposium on the Electronic Structure of Solids; Ziesche, P., Eschrig, H., Akademie Verlag: Berlin, Germany, 1991, p 11. 28 Adamo, C.; Barone, V. J. Chem. Phys. 1998, 108, 664-675. 29 Becke, A. D. J Chem. Phys. 1996, 104, 1040-1046. 30 Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 5656-5667. 31 Vydrov, O. A.; Van Voorhis, T. J. Chem. Phys. 2010, 133, 244103. 32 Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2011, 7, 291–309. 33 Stoll, H.; Pavlidou, C.M.E.; Preuss, H. Theor. Chim. Acta 1978, 49, 143-149. 34 Gori-Giorgi, P.; Sacchetti, F.; Bachelet, G. B. Phys. Rev. B 2000, 61, 7353-7363. 35 Gori-Giorgi, P.; Perdew, J. P. Phys. Rev. B 2004, 69, 041103. 36 Grüneis, A.; Marsman, M.; Harl, J.; Schimka, L.; Kresse, G. E. J. Chem. Phys. 2009, 131, 154115. 37 Dobson, J. F.; Wang, J. Phys. Rev. B 2000, 62, 10038-10045. 38 Jung, J.; García-González, P.; Dobson, J. F.; Godby, R.W. Phys. Rev. B 2004, 70, 205107. 39 Bates, J. E.; Laricchia, S.; Ruzsinszky, A. Phys. Rev. B 2016, 93, 045119. 40 Grimme, S. J. Chem. Phys. 2003, 118, 9095-9102. 41 Grimme, S.; Goerigk, L.; Fink, R. F. WIREs Comput. Mol. Sci. 2012, 2, 886. 42 Chai, J. D.; Head-Gordon, M. J. Chem. Phys. 2008, 128, 084106. 43 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. 44 Zhang, Y.; Xu, X.; Goddard III, W. A. Proc. Natl. Acad. Sci. USA 2009, 106, 4963-4968. 45 Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652. 46 Kozuch, S.; Martin, J. M. J. Comput. Chem. 2013, 34, 2327-2344. 47 Grimme, S.; Ehrlich, S.; Goerigk, L. J. Comput. Chem. 2011, 32, 1456-1465. 48 Goerigk, L.; Grimme, S. Phys. Chem. Chem. Phys. 2011, 13, 6670-6688. 49 Goerigk, L.; Grimme, S. WIREs Comput. Mol. Sci. 2014, 4, 576-600. 50 Heßelmann, A. Phys. Rev. A 2012, 85, 012517. 51 Kállay, M. J. Chem. Phys. 2014, 141, 244113. 52 Wheeler, S. E.; Houk, K. N.; Schleyer, P. V. R.; Allen, W. D. J. Am. Chem. Soc. 2009, 131, 2547-2560. 53 Lynch, B. J.; Truhlar, D. G. J. Phys. Chem. A 2003, 107, 8996-8999.

15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Page 16 of 17

54 Haunschild, R.; Klopper, W. Theor. Chem. Acc. 2012, 131, 1112. 55 Rekkedal, J.; Coriani, S.; Iozzi, M. F.; Teale, A. M.; Helgaker, T.; Pedersen, T. B. J. Chem. Phys. 2013, 139, 081101. 56 Mrcc, a quantum chemical program suite written by M. Kállay, Z. Rolik, J. Csontos, I. Ladjánszki, L. Szegedy, B. Ladóczki, G. Samu, K. Petrov, M. Farkas, P. Nagy, D. Mester, and B. Hégely. See also: www.mrcc.hu. 57 Rolik, Z.; Szegedy, L.; Ladjánszki, I.; Ladóczki, B.; Kállay, M. J. Chem. Phys. 2013, 139, 094105. 58 Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2010, 2, 73-78. 59 Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 18, 3865-3868. 60 Sun, J.; Ruzsinszky, A.; Perdew, J. P. Phys. Rev. Lett. 2015, 115, 036402. 61 Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158-6170. 62 Dunning Jr., T. H.; J. Chem. Phys. 1989, 90, 1007-1023. 63 Woon, D. E.; Dunning Jr., T. H. J. Chem. Phys. 1994, 100, 2975-2988. 64 Kendall, R. A.; Dunning Jr., T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796-6806. 65 Woon, D. E.; Dunning Jr., T. H. J. Chem. Phys. 1993, 98, 1358-1371. 66 Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285-4291. 67 Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175-3183. 68 Hättig, C. Phys. Chem. Chem. Phys. 2005, 7, 59-66. 69 Weigend, F. J. Comput. Chem. 2008, 29, 167-175. 70 Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297-3305. 71 Schäfer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571-2577. 72 Schäfer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829-5835. 73 Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119-124. 74 Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003, 119, 12753-12762. 75 Karton, A.; Tarnopolsky, A.; Lamère, J.-F.; Schatz, G. C.; Martin, J. M. L. J. Phys. Chem. A 2008, 112, 12868–12886. 76 Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2010, 6, 107. 77 Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J. Chem. Phys. 1991, 94, 7221–7230. 78 Zhao, Y.; Lynch, B. J.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 2715–2719. 79 Zhao, Y.; González-García, N.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 2012–2018. 80 Neese, F.; Schwabe, T.; Kossmann, S.; Schirmer, B.; Grimme, S. J. Chem. Theor. Comput. 2009, 5, 3060–3073. 81 Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E. J. Chem. Phys. 2006, 125, 194112. 82 Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 792–794. 83 Fink, R. F. J. Chem. Phys. 2010, 133, 174113. 84 Szabados, Á.; Nagy, P. J. Phys. Chem. A 2011, 115, 523-534. 85 Kesharwani, M. K.; Kozuch, S.; Martin, J. M. L. J. Chem. Phys. 2015, 143, 187101. 86 Zhang, I. Y.; Luo, Y.; Xu, X. J. Chem. Phys. 2010, 133, 104105. 87 Chan, B.; Radom, L. J. Chem. Theory Comput. 2011, 7, 2852-2863.

16 ACS Paragon Plus Environment

Page 17 of 17

1 2 3 4 5 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 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

TOC graphics 85x47mm (300 x 300 DPI)

ACS Paragon Plus Environment