Assessment of the Performance of Density Functionals for Predicting Potential Energy Curves in Hydrogen Storage Applications
AAssessment of the Performance of DensityFunctionals for Predicting Potential EnergyCurves in Hydrogen Storage Applications
Srimukh Prasad Veccham and Martin Head-Gordon ∗ Department of Chemistry, University of California, Berkeley, California 94720, USAChemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California94720, USA
E-mail: [email protected]
Abstract
The availability of accurate computational tools for modeling and simulation is vitalto accelerate the discovery of materials capable of storing hydrogen (H ) under givenparameters of pressure swing and temperature. Previously, we compiled the H2Bind275dataset consisting of equilibrium geometries and assessed the performance of 55 densityfunctionals over this dataset (Veccham, S. P.; Head-Gordon, M. J. Chem. TheoryComput. , , , 4963–4982). As it is crucial for computational tools to accuratelymodel the entire potential energy curve (PEC), in addition to the equilibrium geometry,we have extended this dataset with 389 new data points to include two compressed andthree elongated geometries along 78 PECs for H binding, forming the H2Bind78 × × binding applications: PBE0-DH, ω B97X-V, ω B97M-V, andDSD-PBEPBE-D3(BJ). Addition of Hartree Fock exchange improves the performance a r X i v : . [ phy s i c s . c h e m - ph ] F e b f density functionals, albeit not uniformly throughout the PEC. We recommend theusage of ω B97X-V and ω B97M-V density functionals as they give good performancefor both geometries and energies In addition, we have also identified B97M-V andB97M-rV as the best semi-local density functionals for predicting H binding energyat its equilibrium geometry. Introduction
Hydrogen (H ) is a favorable substitute for fossil fuels as the only by-product of hydrogenfuel cell engines is water and the efficiency of a fuel cell is significantly higher than an internalcombustion engine. However, H is a light gas with low volumetric and gravimetric energydensities. This poses a significant hurdle to storage and transportation of H . Storing H reversibly in adsorbed form on porous materials is a promising solution to this problem. Ideally, such materials should adsorb H at high pressure and release it at low pressureso that the released H can be used for operating a fuel cell. Designing materials with thisproperty, while simultaneously not compromising on high volumetric and gravimetric storagecapacities, is an active area of research. While multiple porous materials like Metal-Organic Frameworks (MOFs), Covalent Or-ganic Frameworks (COFs), graphene, and other amorphous materials have been shown toadsorb H , none of these materials meet all the target criteria proposed by the U.S. Depart-ment of Energy for an ideal storage material. As experimental synthesis and characterizationof potential H storage materials is expensive and time-consuming, computational model-ing and screening of materials has emerged as a viable alternative to it. Computationaltechniques can be used in two different, potentially complementary ways. First, molecularmodeling can be used to understand the mechanism of H binding in different porous mate-rials and this understanding can be used to systematically tune materials to achieve targetproperties. Second, computational techniques can be used to screen materials in a high-throughput manner to select only a handful of potentially viable materials for synthesis and2haracterization.The ability of a material to store H is characterized by its usable capacity, which isdefined as the amount of H stored at the high operating pressure that is released when thepressure is reduced to the low operating pressure. Optimizing the usable capacity for typicalfixed operating pressures of 5 bar and 100 bar gives an optimal value for the Gibbs free energyof adsorption (∆ G ads ). Assuming a correlation between enthalpy and entropy of adsorptionin porous materials gives a range of −
15 to −
25 kJ/mol for the optimal value for enthalpyof adsorption (∆ H ads ). The internal energy of binding, which is the largest componentof ∆ H ads , can be computed using different quantum chemistry methods, including, butnot limited to, density functional theory (DFT), Møller-Plesset perturbation theory(MP2), and different variants of coupled-cluster theory.
Each of these methods havedifferent accuracies and computational costs associated with them.DFT, scaling as O ( N ) ( N is the number of basis functions in the system), can provide areasonable balance between cost and accuracy of computing H binding energy. However, asthe exact density functional remains unknown, different density functional approximations(DFAs), proposed in lieu of the exact density functional, provide varying accuracies fordifferent chemical systems and/or properties computed. In order to address this problem,we adopted a two-pronged approach. (1) We compiled the H2Bind275 dataset that consistsof H (s) interacting with binding motifs representative of different porous materials knownfor H adsorption. This dataset consists only of equilibrium geometries, that is, H (s) arelocated at the minimum of the potential energy curve (PEC) with respect to the binding site.We computed highly accurate reference interaction energies using coupled-cluster singles,doubles, and perturbative triples (CCSD(T)) extrapolated to the complete basis set limit forthis dataset. (2) We assessed the performance of 55 DFAs and identified the best performingdensity functionals for this dataset. In addition, we also identified inexpensive semi-localdensity functionals which give very good performance for low computational cost and aresuitable for in silico high-throughput screening purposes.3he H2Bind275 dataset, consisting of 275 data points, provides a balanced representa-tion of different H binding mechanisms like polarization, charge transfer interaction, anddispersion. It also captures the chemical diversity of binding motifs that H interactswith in porous frameworks. This dataset assesses the ability of density functionals to repro-duce H binding energies at the minima of the PEC. However, as DFAs are routinely usedfor geometry optimizations and molecular dynamics simulations either directly or indirectly(by generating reference data for training force fields), they should also be able to repro-duce the entire PEC which would ensure accurate nuclear gradients as required for geometryoptimization and molecular dynamics simulations. A strategy of assessing the performanceof DFAs for PECs has been previously employed for other non-covalent interaction energydatasets like S22, S66, and A24. The S22x5 dataset was created from the S22 dataset by including geometries that are shortened and elongated along a well-defined interactioncoordinate. Similarly, the S66x8 and A21x12 extended datasets were created from theS66 and A24 datasets. In order to address this issue for H storage, we have extended the H2Bind275 datasetto include geometries that are located at five different points on 78 separate PECs, not justthe minimum. This extended dataset, hereafter referred to as the H2Bind78 × and the binding motif. Thereference interaction energies were computed using CCSD(T) extrapolated to the completebasis set (CBS) limit using the same strategy outlined in Ref. 21. The performance of 55DFAs were assessed using regularized relative errors metrics by appropriately weighing theerror coming from different points on the PECs. We have analyzed the performance of theseDFAs for the extended dataset by comparing and contrasting it with the performance of theoriginal equilibrium H2Bind275 dataset.This paper is organized as follows. The H2Bind78 × interaction energies across the PEC is discussed andcontrasted with their performance for the previous H2Bind275 dataset. The performance ofDFAs for predicting equilibrium geometries and interaction energies at equilibrium geome-tries is explored. The best DFAs for predicting H binding energies are recommended whileconsidering their computational cost. Computational details
H2Bind78 × Table 1: Number of geometries and data points by chemical categories for the H2Bind78 × s-block ions salts organic ligands transition metals totalgeometries 19 13 5 41 78data points at PEC minimum 38 26 10 82 156data points not at PEC minimum 95 65 25 204 × The H2Bind275 dataset consists of 275 H interaction energies but only 86 unique ge-ometries as many of them have multiple H s. For example, the geometry of CaCl − (H ) has four hydrogen molecules bound to CaCl contributing four data points to the H2Bind275dataset. The H2Bind78 × , was chosen as itis closest to experimentally measurable values. For each minimum geometry, five additionalgeometries were generated by compressing and elongating the distance between the bindingmotif and the center of mass of H (denoted by r eq ). For geometries containing multiple One data point excluded due to convergence issues s bound to a single binding moiety, compressed and elongated geometries were generatedfor only one of the H s. This step was necessary in order to maintain low redundancy inthe dataset and make its size manageable. In addition to this, the interaction energy at theminimum of the PEC was also computed using vertical interaction energy method. In total,this dataset contains 78 adiabatic and 78 vertical interaction energies (a total of 156 datapoints) located at the PEC minimum.In this work, two compressed geometries (0 . r eq and 0 . r eq ) and three elongated ge-ometries (1 . r eq , 1 . r eq , and 1 . r eq ) were considered. These distances were chosen as theyare representative of the PEC in both the compressed and elongated regimes. In a porousmaterial, H interacts with not only its primary binding site but also has secondary interac-tions with other components of the framework. The binding distances of H to its secondaryinteraction sites of the porous material are often longer than their corresponding equilibriumdistances. As a consequence of this, when modeling H in a porous material, the elongatedportion of the PEC is sampled more than the compressed part. Additionally, the compressedportion of the PEC is usually significantly higher in energy (repulsive if compressed enough),and is sampled less often in a molecular dynamics or Monte Carlo simulation. Hence, DFAsshould be able to reproduce the elongated portion of the PEC more faithfully than thecompressed portion. We have included more data points in the elongated regime than thecompressed regime in order to underscore its relative importance. As shown in Table 1,the number of non-equilibrium data points is roughly 2 . × interaction energies with representative bindingmotifs.This dataset, like the H2Bind275 dataset, can also be divided into categories based on thechemical nature of the binding motif as shown in Table 2: (1) s-block ions: consisting of group1 and group 2 bare metal cations with unscreened charge binding one or multiple H s, (2)salts: consisting of small inorganic salts like AlF , CaCl , and MgF binding one or multiple6able 2: All 78 geometries in the H2Bind78 × s-block ions salts organic ligands transition metalsLi + − (H ) n , AlF − H benzene − H MX − H , X=H, F, Cl; M=Cu, Ag, Au n = 1 , , , , , − (H ) n , phenol − H CoF − H Na + − (H ) n , n = 1 , , , − H Cu(OMe) − H n = 1 , , − (H ) n , butene − H CuCN − H Mg − (H ) n , n = 1 , , , − H Sc + − (H ) n , V + − (H ) n , n = 3 , n = 1 , , , − (H ) n , Ti + − (H ) n , n = 2 , − (H ) n , n = 1 , , , + − (H ) n , Mn + − (H ) n , n = 1 , , , n = 1 , , , , , + − (H ) n , n = 1 , , , + − (H ) n , Ni + − (H ) n , n = 1 , + − (H ) n , n = 1 , , + − (H ) n , n = 1 , , , H s, (3) organic ligands: comprising of small aliphatic and aromatic molecules binding oneH , (4) transition metals: including small transition metal complexes and 3d transitionmetal cations binding one or multiple H s. Each of these categories is also representativeof various mechanisms of H binding found in porous materials. For example, H in theorganic ligands category is mostly dispersion-bound. The s-block metals category binds H using a combination of electrostatic and forward charge transfer (H → metal) interactions. This dataset captures both chemical and mechanistic diversity encountered in H bindingto porous materials. For a detailed discussion about different chemical categories in thisdataset, we refer readers to Ref. 21. Reference Binding Energies
Calculation of accurate reference interaction energies is an important task in compiling adataset. Reference interaction energies were computed using coupled-cluster theory withsingles, doubles, and perturbative triples (CCSD(T)) extrapolated to the complete basisset limit. Inspired by the success of composite extrapolation methods for computinghighly accurate reference values, we have developed our own composite extrapolation method7sing focal point analysis for computing accurate reference H binding energies: E ref = E HF/5Z + E MP2/QZ → + δE CCSD(T)/TZ + δE coreMP2/TZ (1) δE CCSD(T)/TZ = E CCSD(T)/TZ − E MP2/TZ (2) δE coreMP2/TZ = E core=0MP2/TZ − E core=nMP2/TZ (3)Here, E ref is the reference energy computed using the composite method, E HF/5Z is theHartree Fock energy computed using a basis set of quintuple-zeta (5Z) quality, E MP2/QZ → is the MP2 correlation energy extrapolated to the complete basis set limit with the 2-pointextrapolation formula using correlation energies computed with quadruple-zeta (QZ) andquintuple-zeta quality basis sets, and δE CCSD(T)/TZ is the difference between the CCSD(T)and MP2 correlation energies computed with a triple-zeta quality basis set. δE coreMP2/TZ is thecore-valence contribution to the correlation energy computed as the difference between MP2correlation energies with ( E core=nMP2/TZ ) and without ( E core=0MP2/TZ ) the frozen-core approximation.This composite method for computing reference H interaction energies ensures that the effectof higher-order excitations neglected in CCSD(T) are sufficiently small. It also ensures thatthe basis set incompleteness errors are small and that both HF and extrapolated correlationenergy components are of complete basis set limit quality. Further details of this scheme canbe found in Ref. 21.The cc-pVnZ ( n =T, Q, or 5) family of basis sets was used for all the HF andcorrelation energy calculations when core electrons were not included in the correlationcalculations. cc-pCVnZ ( n =T, Q, or 5) family of basis sets were employed when someor all of the core electrons were included in the correlation calculations. For transition metals,the cc-pwCVnZ ( n =T, Q, or 5) series of basis sets was used with a neon core excluded inall correlation energy computations. 8 ensity Functional Approximations
55 DFAs, including all the commonly used density functionals, were chosen to perform athorough assessment. We have also included DFAs that have previously shown very goodperformance for a range of non-covalent interaction energy prediction problems representedby multiple datasets. Based on the different quantities DFAs depend on, they are cate-gorized into rungs of the metaphorical Jacob’s ladder. From the first rung of the Jacob’sladder, in which DFAs depend only on electron density, SVWN5 and SPW92
DFAswere chosen. From the second rung called Generalized Gradient Approximation (GGA),12 different DFAs were chosen: the PBE family and its variants (PBE, PBE-D3(0), RPBE, revPBE, and revPBE-D3(op) ), BLYP and BLYP-D3(op) , dispersion-corrected variants of B97 (BLYP-D3(0) and BLYP-D3(BJ) ), PW91, and GAM. From the meta-GGA rung, the different variants of TPSS (TPSS, TPSS-D3(BJ), andrevTPSS ), SCAN and its dispersion-corrected version SCAN-D3(BJ) , MS2 and MS2-D3(op), the combinatorially-optimized B97M-V and B97M-rV were chosen. In addi-tion, mBEEF and the semi-local Minnesota functionals M06-L and MN15-L were alsoincluded in the assessment. Rung four DFAs, containing HF exchange, are generally more ac-curate than semi-local functionals as they partially alleviate the problem of self-interactionerror. In this work, global hybrid density functionals like B3LYP and B3LYP-D3(0), PBE0 and PBE0-D3(BJ), TPSSh and TPSSh-D3(BJ), the M06 family of densityfunctionals (M06, M06-2X, M06-2X-D3(0), and revM06 ), MVSh, and SCAN0 which is the hybrid variant of SCAN are included. Range-separated hybrids, which are hy-brid functionals containing DFT exchange and some HF exchange in the short-range andonly HF exchange in the long range, included in this study are ω B97X-D, ω B97X-D3, ω B97X-V, ω B97M-V, M11 and its revised version revM11. Two screened exchangedensity functionals (HSE-HJS and MN12-SX ), which contain DFT exchange in theshort range and attenuated HF exchange in the long range are also included. Double hy-brid density functionals, which are at the top the Jacob’s ladder classification, contain some9ercentage of correlation energy from wavefunction methods. These DFAs are characterizedby their superior accuracy and increased computational cost in comparison to semi-localand hybrid DFAs. We have included seven double hybrid density functionals in this study:B2PLYP-D3(BJ), XYG3, XYGJ-OS, PBE0-DH, , PTPSS-D3(0), DSD-PBEPBE-D3(BJ), and ω B97M(2). The def2-QZVPPD basis set was used for all DFA calculations with a quadraturegrid of 99 Euler-MacLaurin radial points and 590 Lebedev angular points for integrating theexchange-correlation contribution. SG-1 integration grid was used for integrating the VV10component. The choice of core for frozen core approximation and employment of densityfitting approximation for computing MP2 correlation energy in double hybrid density func-tionals is discussed in Table S1. All the PECs were interpolated using the one-dimensionalAkima interpolator. All computations were performed using Q-Chem 5. Results and Discussion
H2Bind78 × Typically, the coupled cluster reference H interaction energy with the binding motif isstrongest at equilibrium, that is at r eq . This implies that the geometries optimized using ω B97M-V/def2-TZVPD are also close to the CCSD(T)/CBS minima. Fig. 1 shows thedistribution of interaction energies for the entire H2Bind78 × . r eq , . r eq , 1 . r verteq , and1 . r eq ) have attractive interaction energies, with most of them smaller than 100 kJ/mol inmagnitude. Geometries that are stretched by 25% (1 . r eq ) are still attractive in nature, butmost interaction energies are smaller than 60 kJ/mol in magnitude. Geometries stretchedby 50% of their equilibrium distance are bound only weakly with a median binding energyof − . . r eq ) are mostly repulsive with a median interaction energy of +29 .
00 150 100 50 0 50 100Interaction energy (kJ/mol)020406080100120 F r equen cy r eq r eq r eq r eq r eq r eq r verteq Figure 1: Distribution of coupled-cluster reference interaction energies separated by loca-tion on the potential energy curve. The reference vertical interaction energy at equilibrium(1 . r verteq ) is also shown. 11eometry was also chosen in order to sample the repulsive part of the PEC and assess howaccurately different density functionals can reproduce it. The range of interaction energies covered by each PEC is also very large. The coinagemetal containing species are the strongest binders, as illustrated by the extreme example ofAuF which binds H with an interaction energy of − . with an energy of +87 . . r eq , thus spanning an interactionrange of 249 . interacting with a binding moiety has the shapeof a Morse potential. However, there is considerable variation in the well depth, well width,and decay in the long range for different chemical species. This variation can provide someclues into the dominant mechanism of interaction. For example, the AuCl binding motifinteracts with one H with an interaction energy of − . − . . r eq (82 .
7% decrease). This is a sharp decay in the interactionenergy in comparison to the Mg case. In the Mg interacting with one H case, theinteraction energy at equilibrium is − . − . . r eq (61 .
4% decrease). This suggested that the dominant mechanism of interaction in Mg caseis longer-ranged (like permanent electrostatics) than AuCl which is dominated by orbitalcontrolled short-ranged interactions like charge transfer. Performance of Density Functional Approximations on PECs
We will discuss the performance of DFAs using multiple error metrics. Each of these met-rics gives different weights to different aspects of the dataset. First, we will discuss theperformance of DFAs using the root mean square error (RMSE) metric which gives equalimportance to all data points in the H2Bind78 × × Rank DFA RMSE (kJ/mol) DFA RegMAPE (%)1 PBE0-DH 2.9 PBE0-DH 5.02 DSD-PBEPBE-D3(BJ) 3.7 ω B97X-V 5.43 ω B97X-V 4.0 ω B97M-V 6.34 ω B97X-D 4.1 DSD-PBEPBE-D3(BJ) 6.35 PBE0 4.2 XYGJ-OS 6.96 MVSh 4.3 ω B97M(2) 7.47 HSE-HJS 4.3 PBE0 7.68 ω B97M-V 4.5 HSE-HJS 7.69 XYGJ-OS 4.6 B2PLYP-D3(BJ) 8.210 ω B97X-D3 4.8 XYG3 8.511 XYG3 4.8 ω B97X-D 9.012 PTPSS-D3(0) 4.9 B97M-rV 9.013 ω B97M(2) 5.1 B97M-V 9.114 PBE0-D3(BJ) 5.2 SCAN0 9.115 MN15 5.7 PTPSS-D3(0) 9.316 B2PLYP-D3(BJ) 5.8 ω B97X-D3 9.817 SCAN0 5.8 MVSh 10.218 TPSSh 6.0 TPSSh 11.319 revM11 6.3 PBE0-D3(BJ) 11.520 mBEEF 6.7 M11 12.021 revM06 6.9 revTPSS 12.022 B3LYP 7.4 revM06 12.423 revTPSS 7.5 TPSS 13.624 B3LYP-D3(0) 7.5 TPSSh-D3(BJ) 13.825 B97M-V 7.6 B3LYP-D3(0) 14.026 B97M-rV 7.6 oTPSS-D3(BJ) 14.627 TPSS 7.7 TPSS-D3(BJ) 14.728 oTPSS-D3(BJ) 7.8 MN15 15.329 TPSSh-D3(BJ) 8.0 PBE 15.330 MN15-L 8.0 revM11 15.331 revPBE-D3(op) 8.5 BLYP-D3(op) 15.432 TPSS-D3(BJ) 8.5 B3LYP 15.533 revPBE 9.1 SCAN 15.634 MN12-SX 9.2 MN12-SX 16.235 RPBE 9.3 SCAN-D3(BJ) 16.336 M11 9.4 M06 16.937 B97-D3(BJ) 9.7 revPBE-D3(op) 17.038 BLYP-D3(op) 9.8 PW91 17.439 M06 9.9 MS2 17.440 BLYP 10.1 mBEEF 17.441 PBE 10.3 BP86-D3(BJ) 18.442 M06-L 10.5 M06-2X 19.343 BP86-D3(BJ) 10.9 M06-2X-D3(0) 19.844 B97-D3(0) 11.0 MN15-L 19.845 PBE-D3(0) 11.0 PBE-D3(0) 20.146 PW91 11.2 MS2-D3(op) 20.347 GAM 11.8 M06-L 20.748 MS2 12.1 RPBE 20.949 MS2-D3(op) 12.4 BLYP 22.150 SCAN 12.9 revPBE 23.451 SCAN-D3(BJ) 13.2 B97-D3(BJ) 24.252 M06-2X 13.4 GAM 24.753 M06-2X-D3(0) 13.4 B97-D3(0) 28.354 SPW92 32.6 SPW92 63.055 SVWN5 32.7 SVWN5 63.0 × This isclosely followed by ω B97X-V and ω B97X-D, both of which show a similar RMSEs of 4 . . × ω B97X-V in the hybrids rung, and PBE0-DH in the doublehybrids rung. The ranking of ω B97M-V deteriorates in the extended dataset in comparisonto the previous H2Bind275 dataset. Another interesting observation is the improvement inthe ranking of the MN15 density functional which is ranked 15 th in the H2Bind78 × . th with an RMSE of 6 . th and26 th ) in the H2Bind78 × × − . storage applications between 5 and 100 bar,interaction energies in the −
15 to −
25 kJ/mol range would be ideal.
A good error metricshould give more weight to data points in this range by considering the following aspects:1. The H2Bind78 × with different interaction energies at their corresponding equilibrium geometry.Binding motifs that bind H with an interaction energy in the range of −
15 to − BE B L YP - D ( op ) r e v PBE - D ( op ) P W BP - D ( B J ) PBE - D ( ) R PBE B L YP r e v PBE B - D ( B J ) G A M B - D ( ) B M -r V B M - V r e v T PSS T PSS o T PSS - D ( B J ) T PSS - D ( B J ) S C A N S C A N - D ( B J ) M S m BEE F M N - L M S - D ( op ) M - L B X - V PBE H SE - H J S B X - D B X - D PBE - D ( B J ) B YP - D ( ) B YP B M - V S C A N M VS h T PSS h M r e v M T PSS h - D ( B J ) M N r e v M M N - SX M M - X M - X - D ( ) PBE - DH D S D - PBEPBE - D ( B J ) XY G J O S B M ( ) B P L YP - D ( B J ) XY G P T PSS - D ( ) R eg M APE
GGA meta-GGA Hybrid GGA Hybrid meta-GGA Double hybrid 5%8%10%12%15%
Figure 2: Performance of density functional approximations for the H2Bind78 × −
15 to −
25 kJ/molrange, regularized percentage error (in order to avoid small denominators) for interactionenergies weaker than −
15 kJ/mol, and absolute error for interaction energies stronger than −
25 kJ/mol. The error metrics in neighboring ranges are also smoothly interpolated. For thesame amount of error, as percentage error is much larger in magnitude than absolute error,the RegMAPE error metric is able to satisfy criterion (1). For example, an error of 5 kJ/molfor a reference interaction energy of 100 kJ/mol will contribute 5 units to the total errorwhile the same error for a reference interaction energy of 20 kJ/mol will contribute 25 unitsto the total error. For a given binding motif, the equilibrium geometry should be given moreweight as it represents the H interaction with the primary binding site: the main lever totune while designing binding sites. As non-equilibrium geometries are higher in energy, theywould be encountered less frequently in a molecular dynamics or Monte Carlo simulation.Hence, lower weight for non-equilibrium geometries is achieved by using the equilibriumregularization value for non-equilibrium geometries as well. As the equilibrium geometryalways has a stronger interaction energy, the regularized value of its reference interactionenergy, ˜ E ( r eq ), will be larger in magnitude in comparison to the regularized values of non-equilibrium interaction energies ( ˜ E ( αr eq ) , α (cid:54) = 1 . E ( αr eq ) = E DFA ( αr eq ) − E ref ( αr eq )˜ E ( r eq ) , α ∈ { . , . , . , . , . , . } (4)where ∆ E ( αr eq ) is the RegMAPE, E DFA ( αr eq ) and E ref ( αr eq ) are the DFA and referenceinteraction energies at αr eq geometry, and ˜ E ( r eq ) is the regularized interaction energy for16he equilibrium geometry. As the vertical interaction energy is located at the minimum ofthe PEC, the error in this data point is regularized using the vertical reference interactionenergy.The performance of DFAs assessed by the RegMAPE error metric is shown in Fig. 2 andTable (3). While there are some similarities in the relative ordering of density functionals forRegMAPE and RMSE error metrics, there are also noteworthy differences. Again, PBE0-DHshows the best performance with the least RegMAPE of 5.0%. It is followed by the ω B97X-Vand ω B97M-V density functionals which have RegMAPEs of 5.4% and 6.3% respectively. TheDSD-PBEPBE-D3(BJ) density functional, which was the best performing density functionalin the H2Bind275 dataset with RegMAPE of 4.9%, is the fourth best performing DFA forthe H2Bind78 × th and 6 th ), and were recommended as inexpensive alternativesto the best performing and expensive density functionals. However, their performance inthe current H2Bind78 × th and 15 th ) respectively. While this reflects their lacklustreperformance for the non-equilibrium geometries, they still remain the best performing semi-local functionals. The next best performing semi-local density functional, revTPSS, is ranked21 st and shows a RegMAPE of 12.0%. The performance of PBE0 and B2PLYP-D3(BJ) DFAsshows significant improvement relative to their performance in the H2Bind275 dataset withboth density functionals entering the top 10 category for the H2Bind78 × .75 r eq r eq r eq r verteq r eq r eq r eq H distance from binding site2.55.07.510.012.515.0 R M SE ( kJ / m o l ) PBEPBE-D3(0)B3LYPB3LYP-D3(0)
Figure 3: Effect of addition of empirical dispersion corrections on the RMSE of overbinding(PBE) and underbinding (B3LYP) density functionals at different points on the potentialenergy curve.empirical dispersion correction can also be assessed using the mean signed error (MSE) andRegMAPE metrics. Addition of dispersion correction improves the performance only if theparent density functional has a systematic underbinding problem (characterized by a positivevalue of MSE). For example, B3LYP has an MSE of 2.7 kJ/mol and a RegMAPE of 15.5%and is systematically underbinding H (s). Addition of a dispersion correction to B3LYPleads to the B3LYP-D3(0) functional which overcomes this underbinding problem. B3LYP-D3(0) slightly overbinds with an MSE of − . with the bindingsite as dispersion corrections are usually damped in the short range. Dispersion correctionsimprove the performance of underbinding functionals in the short range. However, in thelong range, dispersion corrections overestimate its magnitude, causing BLYP, B3LYP, andrevPBE to overbind in the elongated regime. MSEs and RMSEs of all the DFAs containingdispersion corrections and their corresponding parent functionals is shown in Table S3. r eq r eq r eq r verteq r eq r eq r eq H distance from binding site5101520 R M SE ( kJ / m o l ) PBEPBE0SCANSCAN0
Figure 4: Performance of density functionals of the same family with and without HartreeFock exchange at different points on the potential energy curve.Addition of HF exact exchange is essential to ameliorate the effect of self interactionerror in density functionals. Comparing DFAs belonging to the same family, addition of HFexchange improves the performance of semi-local functionals for the H2Bind78 × th with a RegMAPE of 7.6%. Incontrast, the PBE functional is ranked 29 th with a RegMAPE for 15.3%, more than two timesthat of PBE0. HF exchange is a short-range effect and addition of HF exchange improvesthe performance of density functionals in the short range as shown in Fig. 4. While thehybrid functional performs better than its semi-local counterpart throughout the PEC, itseffect is more pronounced in the compressed region than the elongated region. The SCANand SCAN0 functionals show RMSEs of 21.6 and 9.5 kJ/mol (a difference of 12.2 kJ/mol) at19 . r eq of the PEC. Their RMSEs at 1 . r eq is 3.9 and 2.1 kJ/mol, with the hybrid functionalimproving on the semi-local one by only 1.8 kJ/mol.The RegMAPE error metric gives larger weights to data points whose reference inter-action energies are in the interesting range for H storage. Further, it gives more weightto the equilibrium data point than non-equilibrium data points. The relative weights ofdata points on the PES can be further tuned in order to assess the origin of errors of dif-ferent DFAs. Elongated geometries are encountered more often than compressed geometriesin porous material capable of storing H . Compressed geometries are also much higher inenergy (geometries compressed by 25% are almost always repulsive) and are encounteredless often in simulations. This would suggest retuning the weights of the regularized errorsby giving larger weights to equilibrium and elongated regions of the PEC. The weightedRegMAPE (denoted as wRegMAPE or ∆ E w ) is defined as:∆ E w = (cid:88) i w i ∆ E ( α i r eq )7 s.t. (cid:88) i w i = 7 (5)where ∆ E ( α i r eq ) is the RegMAPE at the point α i r eq defined in Eq. (4). Ensuring thatthe weights sum up to 7 would enable an apples-to-apples comparison of wRegMAPE andRegMAPE. In the case of RegMAPE, w i = 1, for all values of i . Reducing the weightsof the compressed geometries with the scheme shown in Table 4, the wRegMAPE can becomputed using Eq. (5). This wRegMAPE metric, shown in Table 5, gives more weight to theelongated geometries. As the vertical interaction is computed its respective PEC minimum,the 1.0 r verteq data point is assigned a weight equal to that of the adiabatic interaction energyat PEC minimum.Table 4: Weights for different points on the adiabatic PEC and vertical interaction energy forcalculating the weighted regularized mean absolute percentage error (wRegMAPE) metric.PEC location 0 . r eq r eq r eq r verteq r eq r eq r eq Weight ( w i ) 0.75 0.91 1.07 1.07 1.07 1.07 1.07Comparing the magnitude of the wRegMAPE (Table 5) of different DFAs to their cor-20able 5: Weighted regularized mean absolute percentage error (wRegMAPE) for selecteddensity functional approximations.Rank DFA wRegMAPE (%)1 PBE0-DH 4.82 ω B97X-V 5.23 DSD-PBEPBE-D3(BJ) 5.94 ω B97M-V 6.05 XYGJ-OS 6.56 ω B97M(2) 7.07 PBE0 7.48 HSE-HJS 7.49 XYG3 7.910 B2PLYP-D3(BJ) 8.011 B97M-rV 8.412 B97M-V 8.413 ω B97X-D 8.514 SCAN0 8.715 PTPSS-D3(0) 9.1responding RegMAPE (Table 3), we can notice that the wRegMAPEs are slightly smaller.Smaller wRegMAPEs suggest that density functionals perform better for the equilibrium andelongated geometries in comparison to the compressed ones. However, the relative orderingof density functionals remains more or less the same. The top five best performing DFAs(PBE0-DH, ω B97X-V, ω B97M-V, DSD-PBEPBE-D3(BJ), and XYGJ-OS) according to theRegMAPE metric are also the five best performing functionals according to the wRegMAPEmetric. We see that the recently parametrized ω B97M(2) double hybrid density functional,which uses ω B97M-V orbitals, is ranked sixth with wRegMAPE of 7.0%.
Performance of Density Functional Approximation for Geometries
All the equilibrium geometries for this dataset were obtained by geometry optimization usingthe ω B97M-V density functional in the def2-TZVPD basis set. With the exception of twochemical systems (AlF − H and Ti + − (H ) ), the CCSD(T)/CBS equilibrium geometry ofall other chemical systems coincides (up to sampling precision) with the ω B97M-V/def2-TZVPD equilibrium geometry. This further validates the choice of equilibrium geometries21
00 175 150 125 100 75 50 25 0Interaction energy at equilibrium (kJ/mol)12345 W e i gh t Figure 5: Weights of different chemical species as a function of their reference adiabaticinteraction energy at equilibrium.for the H2Bind78 × storage purposes, we have devised a weighting scheme that gives largerweights to more relevant data points. Data points with reference adiabatic interaction energyat equilibrium ( E ref (1 . r eq )) in the range of −
15 to −
25 kJ/mol are given a weight of 5 . −
25 kJ/mol are assigned a weight of 1 . storage (that is at −
20 kJ/mol) contributes 5 units to the totalerror as percentage error metric is used. The RegMAPE metric assigns a weight that is5 times larger to the species in the favorable regime in comparison to the strong binders,thus justifying the weights of 5 . . −
15 kJ/mol are mostly comprised of the organic ligands.These species are ubiquitously found in porous materials capable of adsorbing H (like MOFs)22nd form secondary binding sites for H . As this regime is not as important as the favorableone, it is assigned a weight of 4 .
0. This weighting scheme is used to form the weighted meansigned error (wMSE) and weighted mean unsigned error (wMAE) metrics in Fig. (5). Thereference adiabatic interaction energy at equilibrium decides the weight of the correspondingPEC. As the vertical interaction energy does not lie on the adiabatic PEC, those data pointswere not included in the analyses in this section. SV W N SP W M S - D ( op ) M S S C A N - D ( B J ) S C A N r e v M S C A N M PBE - D ( B J ) M - X M N M - X - D ( ) XY G J O S T PSS h - D ( B J ) B P L YP - D ( B J ) B M - V B M ( ) PBE - DH M - L B X - V T PSS - D ( B J ) P W BP - D ( B J ) PBE m BEE F H SE - H J S PBE M N - SX PBE - D ( ) r e v M G A M D S D - PBEPBE - D ( B J ) M VS h M XY G r e v PBE - D ( op ) B M -r V r e v T PSS B M - V P T PSS - D ( ) o T PSS - D ( B J ) T PSS h B X - D T PSS B L YP - D ( op ) B X - D M N - L B YP - D ( ) B - D ( B J ) B YP B - D ( ) B L YP R PBE r e v PBE E rr o r i n H d i s t an c e s ( i n A ng s t r o m s ) wMAE Negative wMSE Positive wMSE Figure 6: Weighted mean absolute error (wMAE) and weighted mean signed error (wMSE)of equilibrium H distances predicted by different DFAs.Most DFAs predict longer equilibrium binding motif – H distances which is shown as apositive value of wMSE in Fig. 6. With the exception of the LDA density functionals, we seethat all other density functionals which predict shorter equilibrium H distances have a verysmall negative wMSE ( > − . ω B97M(2) density functionalwith a wMAE of 0.09˚A. Both of these DFAs perform much better for equilibrium geometriesin comparison to their performance for PECs. It is also rather surprising to note the less23ood performance of PBE0-DH (wMAE of 0 . rd with wMAE of 0.10 ˚A) but its performance for PECs ismediocre (ranked 19 th with a RegMAPE of 11.5%). However, other top performing DFAsin the energetics category like XYGJ-OS, ω B97M-V, and ω B97X-V also perform well forgeometries giving low wMAEs of 0.10 ˚A, 0.11 ˚A, and 0.11˚A respectively. In particular, ω B97X-V shows no systematic error with virtually zero wMSE. It is also interesting to notethe good performance of some semi-local density functionals like mBEEF and B97M-rVwhich give very low errors despite having no HF exchange. In light of these observations,and given the enhanced computational cost of double hybrid DFA nuclear gradients, onecan use hybrid DFAs like PBE0-D3(BJ) or ω B97M-V to perform a geometry optimizationand then use the optimized geometry to perform a single point interaction energy calculationusing a hybrid or double hybrid functional.Another noticeable tread is the performance of DFAs upon the addition of some form ofempirical dispersion correction. Addition of empirical dispersion corrections to DFAs reducestheir errors for equilibrium H distance prediction when the parent functional overestimatesit. For example, the addition of the D3(op) correction to revPBE decreases its wMSE from0.68˚A to 0.09˚A (concurrently decreasing wMAE from 0.74˚A to 0.23˚A). The performance ofthe commonly used density functional B3LYP and M06-2X is quite poor with large wMSEand wMAEs.Typically, DFAs are used to optimize geometries of complexes containing an H boundto a binding motif. After the geometry optimization has converged to a minimum on thepotential energy surface, the binding energy of H is computed as the difference betweenthe energy of the complex at the minimum of the potential energy surface and energies ofthe binding motif and H in isolation. Alternatively, DFAs can also be used in moleculardynamics and Monte Carlo simulations either directly or indirectly (as reference energiesfor parametrizing force fields). In these typical use cases, error in equilibrium binding24able 6: Performance of density functional approximations (DFAs) for predicting H bindingenergy at equilibrium geometry. The adiabatic regularized mean absolute percentage error(RegMAPE ad ) for 15 best performing DFAs and some commonly used DFAs are shown.Rank DFA RegMAPE ad ω B97X-V 4.72 DSD-PBEPBE-D3(BJ) 4.73 PBE0-DH 5.34 ω B97M-V 6.05 XYGJ-OS 7.16 B97M-rV 7.27 B97M-V 7.28 ω B97M(2) 7.49 ω B97X-D 7.810 XYG3 7.911 B2PLYP-D3(BJ) 8.312 PBE0 8.413 HSE-HJS 8.414 ω B97X-D3 9.515 SCAN0 9.728 B3LYP 15.730 PBE 16.331 SCAN 16.433 revPBE-D3(op) 17.342 mBEEF 19.948 M06-2X 22.351 B97-D3(BJ) 24.452 GAM 26.553 B97-D3(0) 29.425nergy can be attributed to two sources: (1) Inaccurate prediction of equilibrium geometry(2) Incorrect prediction of binding energy for the equilibrium geometry. Using DFAs for mod-eling H binding materials, the DFA equilibrium geometry is typically used for computingthe equilibrium binding energy.We can assess the effect of relaxing the geometry along the PEC (defined by each DFA),the interaction coordinate of H with the binding site. Geometry optimization along thiscoordinate can either improve or deteriorate the performance of density functionals. Theperformance of selected density functionals for predicting the minimum energy geometry onits PEC and the H binding energy for this geometry assessed by the RegMAPE error metricis shown in Table 6 (the complete table showing the performance of all 55 DFAs is includedin Table S4).As both the reference and DFA geometries are relaxed along the potential energy curve,this error metric is adiabatic (ad) in nature as reflected in its subscript (RegMAPE ad ).The top five density functionals by the RegMAPE ad error metric are also the top five bestperformers according to their RegMAPE errors (Table 3), further emphasizing the superiorperformance of these DFAs for computing H interaction energies. While the top five densityfunctionals remain the same, it is interesting to note small changes in their order. ω B97X-V is the best performing DFA with a RegMAPE ad of 4.65% which is very closely followedby DSD-PBEPBE-D3(BJ) with a RegMAPE ad of 4.68%. Another noteworthy difference isthe performance of the B97M-rV and the B97M-V functionals which are ranked sixth andseventh with RegMAPE ad of 7.22% and 7.23% respectively. These DFAs were ranked 12 th and 13 th with RegMAPE of 9.0% and 9.1%. These functionals show a favorable cancellationof error in prediction of H equilibrium binding energies when the equilibrium geometry isalso optimized using the same functional. As these density functionals also do not haveany HF exchange, they are computationally less expensive making them well-suited forapplications in high-throughput material screening. The rVV10 non-local functional, whichis an approximation of the VV10 non-local functional, also allows for efficient evaluation26n a plane wave framework and can be useful for modeling periodic systems like MOFs. Athorough assessment of the DFA geometry relaxed on the entire potential energy surface, notjust along the one-dimensional PEC, is beyond the scope of this work and we refer interestedreaders to Ref. 28 for a detailed discussion of this topic. Conclusions
The H2Bind275 dataset published recently assesses the performance of density function-als for predicting the interaction energy of H with different model binding motifs at theequilibrium geometry. In this work, we have assessed the ability of DFAs to predict H interaction energies with binding motifs accurately throughout the PEC, not just at equilib-rium geometry. To that end, we have extended our previous H2Bind275 dataset by addingtwo compressed and three elongated geometries along the PEC to form the H2Bind78 × × . Reference interaction energies for alldata points were computed using CCSD(T) extrapolated to the complete basis set limit.The performance of 55 DFAs was assessed with the CCSD(T) reference interaction energiesusing multiple error metrics. The RMSE metric is democratic and gives equal important toall the 545 data points. The RegMAPE metric, on the other hand, gives more weight tobinding motifs with interaction energies in the range of −
15 to −
25 kJ/mol at equilibriumgeometry. For each binding motif, the RegMAPE metric is designed to give more weight tothe equilibrium than non-equilibrium data points as the latter are encountered less often inmodeling and simulation. DFAs are also assessed on the basis of their predicted equilibriumgeometry and binding energy at predicted equilibrium geometry.The CCSD(T) reference interaction energies for the H2Bind78 × × binding energy throughout the PEC. The ω B97X-V, ω B97M-V, and DSD-PBEPBE-D3(BJ) density functionals are also top performers. The semi-local density func-tionals, B97M-V and B97M-rV, show poorer performance for the H2Bind78 × × th and 12 th respectively withRegMAPEs of 9.1% and 9.0%. Previously in the H2Bind275 dataset, they were ranked5 th and 6 th with RegMAPE of 6.8%. In general, the good performance of the top den-sity functionals in the H2Bind275 dataset continues for the H2Bind78 × ω B97M-V and ω B97X-V give good performance for both ge-ometries and energies. Using the adiabatic RegMAPE metric (RegMAPE ad ) reveals that thesemi-local DFAs, B97M-V and B97M-rV, show very small errors. They benefit significantlyfrom cancellation between geometry-driven and energy-driven errors. ω B97M-V and ω B97X-V are the only density functionals that are not double hybrids which consistently show good28erformance for all the error metrics (energy and geometry-related) defined in this work. Asthese hybrid functionals have significantly lower computational cost in comparison to doublehybrids, we recommend their usage for H binding applications.The H2Bind78 × × ad ) can be usedfor assessment of other similar datasets with well-defined schemes for weighting differentdata points. As force field parametrization requires good reference energies throughout thepotential energy surface, the top performing density functionals in this work can be used forgenerating them. Acknowledgement
This work was supported by the Hydrogen Materials - Advanced Research Consortium (Hy-MARC), established as part of the Energy Materials Network under the U.S. Department ofEnergy, Office of Energy Efficiency and Renewable Energy, under Contract No. DE-AC02-05CH11231. The following author declares a competing financial interest. M. H. G. is a partowner of Q-Chem, Inc. 29 upporting Information Available
Additional information regarding the parameters used for double hybrid density function-als, performance of density functionals for the non-equilibrium subset of the H2Bind78 × ad error for all DFAs assessed in this work is included in the sup-plementary information. CCSD(T)/CBS and DFA interaction energies for all 55 functionalsare provided in the file supporting information.xlsx. The geometries of all complexes in theH2Bind78 × References (1) Takagi, H.; Hatori, H.; Soneda, Y.; Yoshizawa, N.; Yamada, Y. Adsorptive hydrogenstorage in carbon and porous materials.
Mater. Sci. Eng., B , , 143–147.(2) Thomas, K. M. Hydrogen adsorption and storage on porous materials. Catal. Today , , 389–398.(3) Murray, L. J.; Dinc˘a, M.; Long, J. R. Hydrogen storage in metal–organic frameworks. Chem. Soc. Rev. , , 1294–1314.(4) Park, N.; Choi, K.; Hwang, J.; Kim, D. W.; Kim, D. O.; Ihm, J. Progress on first-principles-based materials design for hydrogen storage. Proc. Natl. Acad. Sci. U. S.A. , , 19893–19899.(5) Allendorf, M. D.; Hulvey, Z.; Gennett, T.; Ahmed, A.; Autrey, T.; Camp, J.;Seon Cho, E.; Furukawa, H.; Haranczyk, M.; Head-Gordon, M.; Jeong, S.;Karkamkar, A.; Liu, D. J.; Long, J. R.; Meihaus, K. R.; Nayyar, I. H.; Nazarov, R.;Siegel, D. J.; Stavila, V.; Urban, J. J.; Veccham, S. P.; Wood, B. C. An assessment ofstrategies for the development of solid-state adsorbents for vehicular hydrogen storage. Energy Environ. Sci. , , 2784–2812.306) Colon, Y. J.; Fairen-Jimenez, D.; Wilmer, C. E.; Snurr, R. Q. High-throughput screen-ing of porous crystalline materials for hydrogen storage capacity near room tempera-ture. J. Phys. Chem. C , , 5383–5389.(7) Thornton, A. W.; Simon, C. M.; Kim, J.; Kwon, O.; Deeg, K. S.; Konstas, K.;Pas, S. J.; Hill, M. R.; Winkler, D. A.; Haranczyk, M.; Smit, B. Materials Genomein Action: Identifying the Performance Limits of Physical Hydrogen Storage. Chem.Mater. , , 2844–2854.(8) Ahmed, A.; Seth, S.; Purewal, J.; Wong-Foy, A. G.; Veenstra, M.; Matzger, A. J.;Siegel, D. J. Exceptional hydrogen storage achieved by screening nearly half a millionmetal-organic frameworks. Nat. Commun. , , 1–9.(9) Kapelewski, M. T.; Geier, S. J.; Hudson, M. R.; Stck, D.; Mason, J. A.; Nel-son, J. N.; Xiao, D. J.; Hulvey, Z.; Gilmour, E.; FitzGerald, S. A.; Head-Gordon, M.;Brown, C. M.; Long, J. R. M2(m-dobdc) (M = Mg, Mn, Fe, Co, Ni) MetalOrganicFrameworks Exhibiting Increased Charge Density and Enhanced H2 Binding at theOpen Metal Sites. J. Am. Chem. Soc. , , 12119–12129.(10) Tsivion, E.; Veccham, S. P.; Head-Gordon, M. High-Temperature Hydrogen Storage ofMultiple Molecules: Theoretical Insights from Metalated Catechols. ChemPhysChem , , 184–188.(11) Garrone, E.; Bonelli, B.; Otero Are´an, C. Enthalpy-entropy correlation for hydrogenadsorption on zeolites. Chem. Phys. Lett. , , 68–70.(12) Bhatia, S. K.; Myers, A. L. Optimum conditions for adsorptive storage. Langmuir , , 1688–1700.(13) Bae, Y. S.; Snurr, R. Q. Optimal isosteric heat of adsorption for hydrogen storage anddelivery using metal-organic frameworks. Microporous Mesoporous Mater. , ,300–303. 3114) Mueller, T.; Ceder, G. A density functional theory study of hydrogen adsorption inMOF-5. J. Phys. Chem. B , , 17974–17983.(15) Tsivion, E.; Long, J. R.; Head-Gordon, M. Hydrogen physisorption on metal-organicframework linkers and metalated linkers: A computational study of the factors thatcontrol binding strength. J. Am. Chem. Soc. , , 17827–17835.(16) Cabria, I.; L´opez, M.; Alonso, J. Hydrogen storage capacities of nanoporous carboncalculated by density functional and Møller-Plesset methods. Phys. Rev. B , ,075415.(17) Cabria, I.; L´opez, M.; Alonso, J. Simulation of the hydrogen storage in nanoporouscarbons with different pore shapes. Int. J. Hydrogen Energy , , 10748–10759.(18) Niaz, S.; Manzoor, T.; Islam, N.; Pandith, A. H. Theoretical investigations on C2H4Nbcomplex as a potential hydrogen storage system, using moller–plesset (MP2) and den-sity functional theory. Int. J. Quantum Chem. , , 449–457.(19) Kocman, M.; Jureˇcka, P.; Dubeck`y, M.; Otyepka, M.; Cho, Y.; Kim, K. S. Choosinga density functional for modeling adsorptive hydrogen storage: reference quantummechanical calculations and a comparison of dispersion-corrected density functionals. Phys. Chem. Chem. Phys. , , 6423–6432.(20) Ma, L.-J.; Jia, J.; Wu, H.-S. Computational investigation of hydrogen storage onscandium–acetylene system. Int. J. Hydrogen Energy , , 420–428.(21) Veccham, S. P.; Head-Gordon, M. Density Functionals for Hydrogen Storage: Definingthe H2Bind275 Test Set with Ab Initio Benchmarks and Assessment of 55 Functionals. J. Chem. Theory Comput. , , 4963–4982.(22) Mardirossian, N.; Head-Gordon, M. ω B97M-V: A combinatorially optimized, range-32eparated hybrid, meta-GGA density functional with VV10 nonlocal correlation.
J.Chem. Phys. , , 214110.(23) Sillar, K.; Hofmann, A.; Sauer, J. Ab Initio Study of Hydrogen Adsorption in MOF-5. J. Am. Chem. Soc. , , 4143–4150.(24) Koizumi, K.; Nobusada, K.; Boero, M. Hydrogen storage mechanism and diffusion inmetal–organic frameworks. Phys. Chem. Chem. Phys. , , 7756–7764.(25) Gr´afov´a, L.; Pitonak, M.; Rezac, J.; Hobza, P. Comparative study of selected wavefunction and density functional methods for noncovalent interaction energy calcula-tions using the extended S22 data set. J. Chem. Theory Comput. , , 2365–2376.(26) Jureˇcka, P.; ˇSponer, J.; ˇCern`y, J.; Hobza, P. Benchmark database of accurate (MP2and CCSD (T) complete basis set limit) interaction energies of small model complexes,DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. , , 1985–1993.(27) ˇRez´aˇc, J.; Riley, K. E.; Hobza, P. S66: A well-balanced database of benchmark inter-action energies relevant to biomolecular structures. J. Chem. Theory Comput. , , 2427–2438.(28) Witte, J.; Goldey, M.; Neaton, J. B.; Head-Gordon, M. Beyond energies: Geome-tries of nonbonded molecular complexes as metrics for assessing electronic structureapproaches. J. Chem. Theory Comput. , , 1481–1492.(29) Rezac, J.; Hobza, P. Describing noncovalent interactions beyond the common approx-imations: How accurate is the gold standard, CCSD (T) at the complete basis setlimit? J. Chem. Theory Comput. , , 2151–2155.(30) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth order com-parison of electron correlation theories. Chem. Phys. Lett. , , 479–483.3331) Tajti, A.; Szalay, P. G.; Cs´asz´ar, A. G.; K´allay, M.; Gauss, J.; Valeev, E. F.; Flow-ers, B. A.; V´azquez, J.; Stanton, J. F. HEAT: High accuracy extrapolated ab initiothermochemistry. J. Chem. Phys. , , 11599–11613.(32) DeYonker, N. J.; Cundari, T. R.; Wilson, A. K. The correlation consistent compositeapproach (ccCA): An alternative to the Gaussian-n methods. J. Chem. Phys. , , 84108.(33) Karton, A.; Rabinovich, E.; Martin, J. M.; Ruscic, B. W4 theory for computationalthermochemistry: In pursuit of confident sub-kJ/mol predictions. J. Chem. Phys. , , 11599.(34) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 theory. J. Chem. Phys. , , 084108.(35) East, A. L. L.; Allen, W. D. The heat of formation of NCO. J. Chem. Phys. , , 4638–4650.(36) Cs´asz´ar, A. G.; Allen, W. D.; Schaefer, H. F. In pursuit of the ab initio limit forconformational energy prototypes. J. Chem. Phys. , , 9751–9764.(37) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-set convergence of correlatedcalculations on water. J. Chem. Phys. , , 9639–9646.(38) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. Theatoms boron through neon and hydrogen. J. Chem. Phys. , , 1007–1023.(39) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecularcalculations. III. The atoms aluminum through argon. J. Chem. Phys. , , 1358–1371.(40) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular34alculations. V. Core-valence basis sets for boron through neon. J. Chem. Phys. , , 4572–4585.(41) Peterson, K. A.; Dunning, T. H. Accurate correlation consistent basis sets for molec-ular core-valence correlation effects: The second row atoms Al-Ar, and the first rowatoms B-Ne revisited. J. Chem. Phys. , , 10548–10560.(42) Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transitionmetals. I. All-electron correlation consistent basis sets for the 3d elements Sc-Zn. J.Chem. Phys. , , 64107.(43) Mardirossian, N.; Head-Gordon, M. Thirty years of density functional theory in com-putational chemistry: an overview and extensive assessment of 200 density functionals. Mol. Phys. , , 2315–2372.(44) Perdew, J. P.; Ruzsinszky, A.; Tao, J.; Staroverov, V. N.; Scuseria, G. E.; Csonka, G. I.Prescription for the design and selection of density functional approximations: Moreconstraint satisfaction with fewer fits. J. Chem. Phys. , , 62201.(45) Dirac, P. A. Note on exchange phenomena in the Thomas atom. Math. Proc. Cam-bridge Philos. Soc. 1930; pp 376–385.(46) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlationenergies for local spin density calculations: a critical analysis. Can. J. Phys. , , 1200–1211.(47) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B , , 13244.(48) Perdew, J.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Sim-ple. Phys. Rev. Lett. , , 3865–3868.3549) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initioparametrization of density functional dispersion correction (DFT-D) for the 94 ele-ments H-Pu. J. Chem. Phys. , , 154104.(50) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved adsorption energetics withindensity-functional theory using revised Perdew-Burke-Ernzerhof functionals. Phys.Rev. B , , 7413.(51) Zhang, Y.; Yang, W. Comment on Generalized gradient approximation made simple. Phys. Rev. Lett. , , 890.(52) Witte, J.; Mardirossian, N.; Neaton, J. B.; Head-Gordon, M. Assessing DFT-D3Damping Functions Across Widely Used Density Functionals: Can We Do Better? J. Chem. Theory Comput. , , 2043–2052.(53) Becke, A. D. Density-functional exchange-energy approximation with correct asymp-totic behavior. Phys. Rev. A , , 3098.(54) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energyformula into a functional of the electron density. Phys. Rev. B , , 785–789.(55) Becke, A. D. Density-functional thermochemistry. V. Systematic optimization ofexchange-correlation functionals. J. Chem. Phys. , , 8554–8560.(56) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersioncorrected density functional theory. J. Comput. Chem. , , 1456–1465.(57) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.;Singh, D. J.; Fiolhais, C. Atoms, molecules, solids, and surfaces: Applications ofthe generalized gradient approximation for exchange and correlation. Phys. Rev. B , , 6671–6687. 3658) Haoyu, S. Y.; Zhang, W.; Verma, P.; He, X.; Truhlar, D. G. Nonseparable exchange–correlation functional for molecules, including homogeneous catalysis involving tran-sition metals. Phys. Chem. Chem. Phys. , , 12146–12160.(59) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density func-tional ladder: Nonempirical meta–generalized gradient approximation designed formolecules and solids. Phys. Rev. Lett. , , 146401.(60) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Workhorsesemilocal density functional for condensed matter physics and quantum chemistry. Phys. Rev. Lett. , , 026403.(61) Sun, J.; Ruzsinszky, A.; Perdew, J. P. Strongly constrained and appropriately normedsemilocal density functional. Phys. Rev. Lett. , , 036402.(62) Brandenburg, J.; Bates, J.; Sun, J.; Perdew, J. Benchmark tests of a strongly con-strained semilocal functional with a long-range dispersion correction. Phys. Rev. B , , 115144.(63) Sun, J.; Haunschild, R.; Xiao, B.; Bulik, I. W.; Scuseria, G. E.; Perdew, J. P. Semilocaland hybrid meta-generalized gradient approximations based on the understanding ofthe kinetic-energy-density dependence. J. Chem. Phys. , , 044113.(64) Mardirossian, N.; Head-Gordon, M. Mapping the genome of meta-generalized gradientapproximation density functionals: The search for B97M-V. J. Chem. Phys. , , 074111.(65) Sabatini, R.; Gorni, T.; De Gironcoli, S. Nonlocal van der Waals density functionalmade simple and efficient. Phys. Rev. B , , 4–7.(66) Mardirossian, N.; Ruiz Pestana, L.; Womack, J. C.; Skylaris, C.-K.; Head-Gordon, T.;Head-Gordon, M. Use of the rVV10 nonlocal correlation functional in the B97M-V37ensity functional: defining B97M-rV and related functionals. J. Phys. Chem. Lett. , , 35–40.(67) Wellendorff, J.; Lundgaard, K. T.; Jacobsen, K. W.; Bligaard, T. mBEEF: An accuratesemi-local Bayesian error estimation density functional. J. Chem. Phys. , ,144107.(68) Zhao, Y.; Truhlar, D. G. A new local density functional for main-group thermochem-istry, transition metal bonding, thermochemical kinetics, and noncovalent interactions. J. Chem. Phys. , , 194101.(69) Yu, H. S.; He, X.; Truhlar, D. G. MN15-L: A new local exchange-correlation functionalfor Kohn–Sham density functional theory with broad accuracy for atoms, molecules,and solids. J. Chem. Theory Comput. , , 1280–1293.(70) Becke, A. D. Densityfunctional thermochemistry. III. The role of exact exchange. J.Chem. Phys. , , 5648–5652.(71) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustableparameters: The PBE0 model. J. Chem. Phys. , , 6158–6170.(72) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative assessment ofa new nonempirical density functional: Molecules and hydrogen-bonded complexes. J.Chem. Phys. , , 12129–12137.(73) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group ther-mochemistry, thermochemical kinetics, noncovalent interactions, excited states, andtransition elements: two new functionals and systematic testing of four M06-classfunctionals and 12 other functionals. Theor. Chem. Acc. , , 215–241.(74) Wang, Y.; Verma, P.; Jin, X.; Truhlar, D. G.; He, X. Revised M06 density functional38or main-group and transition-metal chemistry. Proc. Natl. Acad. Sci. U. S. A. , , 10257–10262.(75) Sun, J.; Perdew, J. P.; Ruzsinszky, A. Semilocal density functional obeying a stronglytightened bound for exchange. Proc. Natl. Acad. Sci. U. S. A. , , 685–689.(76) Hui, K.; Chai, J.-D. SCAN-based hybrid and double-hybrid density functionals frommodels without fitted parameters. J. Chem. Phys. , , 044114.(77) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals withdamped atom–atom dispersion corrections. Phys. Chem. Chem. Phys. , , 6615–6620.(78) Lin, Y.-S.; Li, G.-D.; Mao, S.-P.; Chai, J.-D. Long-range corrected hybrid densityfunctionals with improved dispersion corrections. J. Chem. Theory Comput. , ,263–272.(79) Mardirossian, N.; Head-Gordon, M. ω B97X-V: A 10-parameter, range-separated hy-brid, generalized gradient approximation density functional with nonlocal correlation,designed by a survival-of-the-fittest strategy.
Phys. Chem. Chem. Phys. , ,9904–9924.(80) Peverati, R.; Truhlar, D. G. Improving the accuracy of hybrid meta-GGA densityfunctionals by range separation. J. Phys. Chem. Lett. , , 2810–2817.(81) Verma, P.; Wang, Y.; Ghosh, S.; He, X.; Truhlar, D. G. Revised M11 Exchange-Correlation Functional for Electronic Excitation Energies and Ground-State Proper-ties. J. Phys. Chem. A , , 2966–2990.(82) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of theexchange screening parameter on the performance of screened hybrid functionals. J.Chem. Phys. , , 224106. 3983) Henderson, T. M.; Janesko, B. G.; Scuseria, G. E. Generalized gradient approximationmodel exchange holes for range-separated hybrids. J. Chem. Phys. , , 194105.(84) Peverati, R.; Truhlar, D. G. Screened-exchange density functionals with broad ac-curacy for chemistry and solid-state physics. Phys. Chem. Chem. Phys. , ,16187–16191.(85) Grimme, S. Semiempirical hybrid density functional with perturbative second-ordercorrelation. J. Chem. Phys. , , 6158.(86) Zhang, Y.; Xu, X.; Goddard, W. A. Doubly hybrid density functional for accuratedescriptions of nonbond interactions, thermochemistry, and thermochemical kinetics. Proc. Natl. Acad. Sci. U. S. A. , , 4963–4968.(87) Zhang, I. Y.; Xu, X.; Jung, Y.; Goddard, W. A. A fast doubly hybrid density functionalmethod close to chemical accuracy using a local opposite spin ansatz. Proc. Natl. Acad.Sci. U. S. A. , , 19896–19900.(88) Brmond, E.; Adamo, C. Seeking for parameter-free double-hybrid functionals: ThePBE0-DH model. J. Chem. Phys. , , 024106.(89) Goerigk, L.; Grimme, S. Efficient and accurate double-hybrid-meta-GGA densityfunctionals- evaluation with the extended GMTKN30 database for general main groupthermochemistry, kinetics, and noncovalent interactions. J. Chem. Theory Comput. , , 291–309.(90) Kozuch, S.; Martin, J. M. L. Spin-component-scaled double hybrids: An extensivesearch for the best fifth-rung functionals blending DFT and perturbation theory. J.Comput. Chem. , , 2327–2344.(91) Mardirossian, N.; Head-Gordon, M. Survival of the most transferable at the top of40acob’s ladder: Defining and testing the ω B97M(2) double hybrid density functional.