Toward a Quantum-Chemical Benchmark Set for ... - ACS Publications

Subscriber access provided by UNIV OF SOUTHERN INDIANA

Article

Toward a Quantum-Chemical Benchmark Set for Enzymatically Catalyzed Reactions: Important Steps and Insights Dominique Ann Wappett, and Lars Goerigk J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b05088 • Publication Date (Web): 17 Jul 2019 Downloaded from pubs.acs.org on July 17, 2019

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 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 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.

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 62 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 Journal of Physical Chemistry

Toward a Quantum-Chemical Benchmark Set for Enzymatically Catalyzed Reactions: Important Steps and Insights Dominique A. Wappett and Lars Goerigk∗ School of Chemistry, The University of Melbourne, Victoria 3010, Australia E-mail: [email protected]

1 Environment ACS Paragon Plus

The Journal of Physical Chemistry 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

Dedicated to Prof. Leo Radom on the occasion of his upcoming 75th birthday. Abstract We present steps toward the construction of a quantum-chemical benchmark set for enzymatically-catalyzed reactions comprising different—quite sizeable—model systems for five reactions. Our work is inspired by a study of semi-empirical molecularorbital theory methods against B3LYP-based structures and energies [Kromann J. C.; Christensen, A. S.; Cui, Q.; Jensen, J. H. PeerJ 2016, 4, e1994.]. Based on the same model systems, we herein demonstrate first how adequate treatments of London dispersion and basis-set superposition error lead to structural differences compared to structures obtained at the popular B3LYP/6-31G(d,p) level. The changes observed by us have a significant influence on the final reaction barrier heights (BHs) and energies (REs), which are our main focus herein. We then proceed and carefully investigate different strategies to either obtain exact or estimated complete-basis-set (CBS) BHs and REs with the help of the DLPNO-CCSD(T) approach, which offers an exciting path into exploring larger biochemically relevant structures at the coupled cluster level. Based on our analysis, we present three strategies how to obtain reliable values. While these recommendations constitute our main findings and are meant to be used as useful guidelines for others in the field, we demonstrate their applicability in a preliminary benchmark set comprising 16 BHs and 12 REs on which we assess 35 density functional theory approximations. This preliminary benchmark study was used as an indicator for the validity of our approach toward generating such a benchmark set, based on the fact that we were able to reproduce major recommendations in the field of modern DFT, such as the reproduction of the famous Jacob’s Ladder scheme, while also pointing out subtle differences to previous benchmark studies. We end our study by critically assessing the still popular choice to gauge the accuracy of low-level methods against B3LYP data and we show the risks of such a strategy, adding to the ever-growing list of calls to move away from


Page 2 of 62

Page 3 of 62 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


B3LYP as a standard black-box approach.

Introduction Benchmarking is an essential tool in quantum-chemical research. Not only does it help method developers to gauge the accuracy of new methods, but more importantly it can provide fundamental insights that are of interest to the general chemist. Successful examples for the latter are carefully conducted, comprehensive benchmark studies in the field of Density Functional Theory (DFT) that have provided the general method user—experts and non-experts alike—with guidelines that allow them to navigate through the confusing ‘zoo’ of density functional approximations (DFAs), 1–6 address popular misconceptions in the field, 1,3,7–10 show that the popularity of a given method is not a scientifically reasonable justification for using it, 11,12 and attempt to close the communication gap between the developer and user communities. 2,3,12 Even more successful examples are benchmark studies in the field of London-dispersion interactions that showed how these influence chemical reactivity and structures, 1,3,13–20 thus inspiring experimentalists to synthesize compounds with unusual structural features 21–23 that are facilitated though ‘dispersion energy donors’. 16 Benchmark sets can in principle focus on any conceivable property—given that it is possible to establish an accurate reference against which lower-level methods can be compared—with single-property benchmark sets focussing on, for example, electron affinities and ionization potentials, 24 heats of formation and total atomization energies, 25–30 intermolecular and intramolecular noncovalent interactions, 13,15,31–42 reaction energies, 3,43–50 barrier heights, 3,51–55 vibrational frequencies, 56–61 magnetic properties, 62–64 and electronicexcitation energies, 65–74 where the list of examples given is by far not complete. Benchmark databases comprised of sets covering various properties have also been published, as



they are used to identify methods that perform well over a range of different properties, which can then subsequently be used for the application to previously unexplored problems. 1–3,5,6,75–80 On the occasion of Prof. Leo Radom’s 75th birthday we would also like to take the opportunity to point out his immense contributions to the field of benchmarking that covered areas such as vibrational frequencies and frequency scale factors, 56,57,61 bond dissociation energies and radical stabilization energies, 81–86 non-covalent interactions in peptide conformers, 40 thermochemistry of open shell systems, 87–91 and barrier heights (BHs) in proton-exchange, 55 proton-transfer 92 and hydrogen-abstraction 93 reactions, to name only a few of the many examples. In recent years, advances in hard- and software have made the study of large systems with one hundred or more atoms more achievable and some benchmark sets focussing particularly on biochemically relevant questions have gained interest in the community. 38,41,94 Our work is inspired by one of those benchmark studies, namely one that focussed on enzymatically catalyzed reactions entitled ‘Towards a barrier height benchmark set for biologically relevant systems’. 95 For that preliminary model study, Jensen and co-workers collected quantum-chemically optimized structures for five reactions separately published by Himo and co-workers (see figure 1), 96–100 and assessed various semi-empirical molecular orbital (MO) theories 101–103 with respect to their BHs and reaction energies (REs). As the treatment of the actual enzymes was not feasible, the authors presented for each reaction a set of different model systems ranging from 27 to 220 atoms, with the smallest models normally including the active sites and reactants, intermediates or products, while larger models also included more of the surrounding amino-acid residues. 95–100 The model structures from the literature were exclusively based on the popular B3LYP 104,105 /6-31G(d,p) 106 level of theory, while the reference BHs and REs were based on B3LYP/6-311+G(2d,2p) 107 single-point energies. Jensen et al. concluded their paper by encouraging others to use their computed structures and energies to use them toward building a new benchmark 4 Environment ACS Paragon Plus

Page 4 of 62

Page 5 of 62 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


set. 95

Figure 1: Reaction schemes for the five enzymes in the benchmark set. The acronyms are discussed in the following section. The scheme for the PTE reaction is simplified for clarity, with two histidine ligands on each zinc ion omitted.

In this paper, we accept Jensen et al.’s invitation and report the next significant steps toward such a set. However, besides the fact that enzymatically catalyzed reactions by themselves are important to study computationally, we also intend to raise awareness of a problem that to this date still exists in many DFT-based studies, namely the reliance on dispersion-uncorrected DFT approaches, such as B3LYP, and the usage of small atomicorbital (AO) basis sets; for a recent account addressing both issues for a non-expert user readership, see Ref. 12. Already in 2012 Kruse et al. pointed out why the B3LYP/6-31G* level of theory was inadequate for computational thermochemistry, 11 and in 2013 and 2014 Goerigk et al. followed up by demonstrating its inadequacy for gas-phase structures 5 Environment ACS Paragon Plus


of di- and tripeptides and for a protein-crystal fragment. 17,19 The main problem of just using B3LYP for a geometry optimization is the lack of an adequate treatment of London dispersion. That missing attractive force can be sometimes partially mimicked by basis-set superposition error (BSSE) stemming from small basis sets, however, this inherent error cancellation cannot be controlled. As a consequence, structures can be artificially distorted and chemically wrong. 17 Given the size of polypeptides and protein fragments, Grimme’s DFT-D3 108,109 dispersion correction and Kruse and Grimme’s geometrical counterpoise 110 (gCP) correction for BSSE are easy-to-use and time-efficient fixes to both problems that, combined with the right DFT method, provide the right structure for the right reasons, as pointed out by Goerigk and Reimers already in 2013. 17 When turning to the calculation of BHs and REs, 6-311+G(2d,2p) may be a larger AO basis set, however the lack of dispersion in the published B3LYP data must be pointed out. In addition, despite being one of the probably most popular DFT methods in the user community, our group’s recent large study of more than 325 dispersion-corrected and -uncorrected DFT methods for general main-group thermochemistry, kinetics and noncovalent interactions revealed that B3LYP itself ranked only in 197th position. 3,5,6,12 While the application of a dispersion correction improved its performance, even B3LYP-D3(BJ) still ranked in 72nd position and it was outperformed by a plethora of more accurate and robust dispersion-corrected DFAs. 3,5,6,12 This in itself does not only question common practice in the field, but more importantly how reliable B3LYP-based numbers are as a reference for the assessment of other quantumchemical methods; for other examples of how the quality of the reference benchmark affects the ranking of DFT methods, see Refs 3 and 53. Those familiar with the works by Radom and co-workers may also appreciate the thoroughness and detail of their studies. This work is inspired by that level of detail and we present an in-depth analysis of the importance of accurately handling dispersion and BSSE in geometry optimizations of the structures initially collated in Ref. 95. We then proceed 6 Environment ACS Paragon Plus

Page 6 of 62

Page 7 of 62 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


with a thorough analysis of different coupled-cluster based protocols and recommend a series of strategies to obtain reliable reference data for the same structures. We conclude our study by presenting a preliminary version of an updated version of the benchmark set from Ref. 95 with a brief study of DFT methods. As such, our study can be regarded as the next step toward an accurate benchmark set for enzymatically catalyzed reactions. Note that our main intention is not to present a final benchmark set and to focus on the DFT analysis, but on possible ways leading toward such a set. As such, our main focus will be on geometry optimizations and the generation of reference data, insights that can be used by others in a similar context in other studies.

The reactions and their various models The active-site models used in Ref. 95 and reused herein cover the five reactions catalyzed by the enzymes L-aspartate α-decarboxylase (AspDC), histone-lysine methyltransferase (HKMT), 4-oxalocrotonate tautomerase (4-OT), phosphotriesterase (PTE), and haloalcohol dehalogenase (HheC). The relevant reaction schemes are shown in figure 1. For all reactions, the models were created using the high-resolution crystal structures of the proteins. Details of the different model systems for each reaction are given in table S1 in the Supporting Information (SI). The first reaction, mediated by AspDC, involves the removal of carbon dioxide from L-aspartate, giving β-alanine. Work was initially done on this system by Liao, Yu and Himo. 96 Their work involved creating seven model systems of this reaction, ranging in size from 27 to 220 atoms. The smallest model system, called ‘model 0’, 95 consists of only the substrate with a pyruvoyl group covalently bonded to the nitrogen—no active site residues are included; this model is the only model of any of the five reactions that does not contain any parts of the enzyme. The substrate is, thus, protonated in model 0 to



reduce charge-delocalization issues, as the anionic group of the substrate forms salt bridges to a nearby arginine residue in the complete enzyme and also in the larger model systems, models 1-5, with the latter containing an increasing number of active site residues around the substrate. The second reaction is the methylation of the N -terminal histone tail of chromatin, which is catalyzed by HKMT. Initial studies were conducted by Georgieva and Himo, 97 and four model systems, models 0-3, of the reaction were created. In the smallest model, the active site is truncated leaving only two carbon atoms either side of the sulfur, as this is the smallest portion of the enzyme that must be included while still retaining the properties of the sulfur-carbon bonds. Again, larger structures include more residues around the active site. The third reaction is the isomerization of unconjugated alpha-keto acids to their conjugated isomers by 4-OT. This was taken from initial work by Sevastik and Himo, 98 in which two model systems for this reaction had been created. It has been shown 98 that three arginine residues are important for the reaction, as two bind to the carboxylate groups of the substrate, while a third stabilizes the negative charge formed in the first transition state and intermediate by hydrogen bonding to the oxygen in the carbonyl group. The smallest model of the system contains the substrate and these three arginine residues, as well as a catalytic proline residue and two water molecules. The fourth reaction is the hydrolysis of organophosphate triesters catalyzed by PTE. Initial computational analysis of this enzyme was conducted by Chen, Fang and Himo. 99 Only one model system, model 0, was created for this reaction, and it involves the substrate, the two ionic zinc centers and their first-shell ligands: the bridging hydroxide, four histidines, a carboxylated lysine and an aspartate. The ligands were truncated so that only the side chains were included. The final reaction is the dehalogenation of vicinal haloalcohols by HheC to form epox8 Environment ACS Paragon Plus

Page 8 of 62

Page 9 of 62 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


ides. Initial work on this enzyme, which catalyzes the reaction of (R)-2-chloro-1-phenylethanol to give (R)-styrene oxide, was done by Hopmann and Himo. 100 Three models of this system were constructed, with the smallest containing the side chains of arginine, tyrosine and serine residues, which act as catalysts for the reaction, as well as some neighboring residues which interact with the catalytic triad. The smallest model also uses a charge-separated structure, where arginine is protonated and tyrosine is deprotonated, as this is more stable than the neutral residues.

Computational details Nearly all calculations were conducted in ORCA versions 4.0.0 and 4.1.2 111,112 with the default settings for self-consistent field (SCF) and geometry convergence. DFT calculations were carried out with ORCA’s multigrid option, using the “grid3 finalgrid5” settings. Some DFT calculations were carried out with a local version of ORCA 4 or with Turbomole 7.1.1, 113–115 as specified later. Turbomole calculations were carried out with the multigrid option “m4” and an SCF convergence criterion of 10−7 Eh , with the exception for the SCAN 116 DFA, for which the grid options “gridsize 4” and “radsize 40” were necessary, as documented in the literature. 3,116,117 Most DFT calculations were sped up with the resolution-of-the-identity (RI-J) approximation for Coulomb integrals. 118,119 An additional way to save computer time is to speed up the SCF step in an HF or hybrid-DFT calculation with the resolution of the identity for Coulomb integrals and the chain-of-spheres approximation for exchange integrals (RIJCOSX), 120 a standard setting in ORCA, which was sometimes applied (mostly to the ωB97X/M-type hybrids). All DFAs belonging to the generalized-gradient approximation (GGA) and meta-GGA classes, as well as most (double-)hybrids, were used together with the RI-J approximation. The frozen-core approximation was used for double-hybrid DFAs to prevent BSSE-effects in the treatment



of core-core electron correlation. 121 Specific details on tested DFAs, AO basis sets, and technical settings regarding wave-function calculations are provided in later sections.

Obtaining new geometries Technical details on geometry optimizations It has been established that any DFT or Hartree-Fock (HF) optimization of larger molecules requires both a dispersion and a BSSE correction if small AO basis sets are chosen. Works by Goerigk et al. recommend the use of the PW6B95 122 -D3(BJ) 108,109 /def2-TZVP 123 or HF-D3(BJ)-gCP/def2-SV(P) 123 levels of theory for calculating accurate polypeptide and protein structures. 17,19 The more cost-efficient HF-3c method, 124 which is HF-D3(BJ)gCP with a minimal basis set and a basis-set-incompleteness error (BSIE) correction, was recommended as a cost-effective alternative. Since the inception of HF-3c, two more methods have emerged, namely PBEh-3c 125 —a PBE-based 126 hybrid combined with a special double-ζ basis and enhanced with DFT-D3(BJ) and a modified gCP correction—as well as very recently the B97-3c 127 approach, which is a refit of the B97-D 128 GGA functional with a triple-ζ AO basis, DFT-D3(BJ) and a correction that addresses the known problem of systematically elongated bonds in GGA geometries. Based on the published literature at the time we started to undertake the present study, we initially chose PW6B95-D3(BJ)/def2TZVP as the level of theory for conducting the geometry optimizations of the smallest model, and then for efficiency reasons turned to PBEh-3c for all structures, including the larger models. In all models except the smallest model system of AspDC, the positions of some atoms had to be constrained. As the systems are active-site models designed to represent the mechanism of action of the enzyme as a whole, the model must remain as close to the Xray crystallographic structure used to create it. Most importantly, the outermost aminoACS Paragon 10 Plus Environment

Page 10 of 62

Page 11 of 62 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


acid residues and any parts of the peptide backbone that are included must not move and, thus, the original coordinates of certain atoms on these groups were fixed to retain their correct positions. Initially, we followed the guidelines from Ref. 95 regarding which atomic coordinates to constrain, but we found that the lists of coordinates published in the supporting material of Ref. 95 did not match the descriptions and diagrams in the original papers published by Himo and coworkers. 96–100 As such, we determined which atoms needed to be constrained in each model system from the details in those articles. Additionally, constraints were determined for the PTE model system for which no constraints had been published in either Ref. 95 or 99. The full list of atoms that were constrained in these calculations for each model system can be found in table S2 in the SI. The finally optimized structures can also be obtained as xyz files from the SI.

Validation of the new geometries We first compare the structures optimized at the PBEh-3c level of theory with the original B3LYP/6-31G(d,p) structures to see the effects of the dispersion and BSSE corrections. Representative comparisons can be seen in figure 2. In the smaller model systems and systems with a higher proportion of constrained atoms, the differences between the original structures are negligible. This is demonstrated in figure 2a, which compares the structures for HKMT model 1. This system is the third smallest model in the benchmark set, and contains only 46 atoms. As there are constraints applied to both molecules in the structure, there is little flexibility during the optimization, resulting in very little difference between the structures, with a root mean squared deviation (RMSD) of only 0.12 ˚ A. It is expected that with larger structures the effects of missing London dispersion and distortion due to BSSE become more apparent. Indeed, when looking at larger models— particularly those with proportionally fewer constraints—we observe more significant differences between the structures. In figure 2b, the effect of the dispersion correction can ACS Paragon 11 Plus Environment


Figure 2: Comparison between the original B3LYP structure (shown in black) and the reoptimized PBEh-3c structure (shown in yellow) of a) the HKMT model-1 reactant, b) the PTE model 0-reactant and c) the HheC model-0 reactant. One histidine ligand has been hidden in part (b) to make the reactant’s depiction clearer. be seen for PTE model 0. On the right hand side of the image there are two aromatic groups—a histidine ligand and a phenyl derivative in the phosphotriester substrate. In the PBEh-3c structure, these groups are closer to each other than in the B3LYP structure, which is due to the dispersion-induced stabilization. These differences in the structure result in an RMSD of 0.61 ˚ A. A similar effect of adding a dispersion term to the computational treatment can be seen in figure 2c, which shows the differences between the structures for HheC model 0, for which we observe a favorable interaction between a phenyl group from the substrate and a tyrosine side chain. In the PBEh-3c structure, the groups lie significantly closer together than in the original B3LYP structure without dispersion, with an arrangement that resembles off-center parallel π-stacking. This arrangement is stabilized by dispersion

ACS Paragon 12 Plus Environment

Page 12 of 62

Page 13 of 62 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


interactions, and thus is better represented by the dispersion-corrected method. The same model also contains a hydrogen bond between the part of the protein backbone connecting the Leu178 and Tyr177 residues and a water molecule, as seen in the top left of the structure. As approximately 25% of the strength in a hydrogen bond is due to dispersion interactions, 129 the hydrogen bond is represented better by the DFT-D3(BJ) model. This is demonstrated by the length of the hydrogen bond, which is 3.5 ˚ A in the original B3LYP structure, and 2.0 ˚ A in the PBEh-3c structure, with the latter being close to the the typical range for hydrogen bonds of this type. The RMSD between these PBEh-3c and B3LYP-based structures is 0.85 ˚ A. Visual structure comparisons are only included for smaller systems due to difficulties in obtaining clear pictures of the larger systems, so RMSDs between all obtained structures and the original B3LYP structures are given in table S4 in the SI to show that similar differences are observed for all model systems. In general, low RMSDs are seen for smaller systems and systems with many constraints where consequently the differences are only for a few small sections of the structure. Higher RMSDs are achieved in larger structures. For the model-0 systems for which PW6B95-D3(BJ)/def2-TZVP structures were also obtained, comparisons can be made with the PBEh-3c structures to determine whether they give similar results. A comparative analysis is shown in figure S1 (SI). While PBEh3c and PW6B95-D3(BJ) structures are not exactly identical, the differences between the two methods are small enough to justify exploiting the cost advantage PBEh-3c offers. Ultimately, describing the underlying physics better, the new structures are significantly more insightful than the original B3LYP ones. Moreover, studies on di- and tripeptide structures showed good performance for both PBEh-3c and PW6B95-D3. 17,125 As such, PBEh-3c is appropriate for the geometry optimizations conducted in this work.



Establishing reference values Coupled cluster singles doubles with quasi-perturbative triples 130 [CCSD(T)] is often labelled the “gold standard” of chemical accuracy (1 kcal/mol threshold for REs and BHs). When higher-order excitations are not computationally feasible, CCSD(T) complete-basisset (CBS) numbers are therefore ideal ab initio reference numbers for the assessment of lower-level methods. However, canonical CCSD(T) is often too expensive for chemically interesting systems, as the cost of the calculation scales as O(N 7 ), where N is the system size or number of AOs. Possible alternatives are explicitly correlated versions, 131–133 or methods that break down the entire space into separate local domains for which the electronic correction energy is calculated. For the present study, we considered Neese and co-workers’ Domain-based Local Pair Natural Orbital 134 (DLPNO) approximation in its latest linear-scaling implementation 135 in ORCA4 as ideally suited. Our aim is therefore to obtain DLPNO-CCSD(T)/CBS reference BHs and REs. However, despite having determined our reference method of choice, there are many factors that need to be considered when calculating DLPNO-CCSD(T) energies for the herein discussed systems. In particular, we need to address the effect of the choice of basis sets on the CBS extrapolations, the influence of three main parameters in the DLPNO approach, and the treatment of large systems, for which a conventional CBS extrapolation may go beyond computational resources. In case of the latter, we investigate various composite schemes that employ the DLPNO technique. All our tests were conducted using the HKMT model-0 structures optimized at the PBEh-3c level of theory, as this is one of the smallest model systems and the calculations were reasonably fast. This was followed up by assessing the AspDC model 0 to ensure that any trends observed for the first can also be transferred to the second system. Tests for both systems were also carried out with PW6B95-D3(BJ)/def2-TZVP and the original B3LYP/6-31G(d,p) structures


Page 14 of 62

Page 15 of 62 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


to determine the effect of the underlying geometry on the resulting BHs and REs. The guidelines we will end up determining in this part will not only be used in the later stage for our own preliminary benchmark study, but are important insights for others in the field.

DLPNO-CCSD(T)/CBS variants The first feature of the calculations tested was the choice of basis sets used in the CBS extrapolation. Herein, we followed the usual two-point schemes with individual extrapolations of HF energies 136 (ESCF ) and electron-correlation energies 137 (ECorr ), which are given as: (CBS)

ESCF

√ √ (X) (Y ) ESCF · exp(−α Y ) − ESCF · exp(−α X) √ √ = , exp(−α Y ) − exp(−α X)

(1)

and (X)

(CBS) ECorr

(Y )

X β · ECorr − Y β · ECorr , = Xβ − Y β

(2)

where X and Y are the cardinal numbers of the basis sets with X > Y , and α and β are optimized constants specific to the type of basis sets used (see table S5 for values). 138 To do the extrapolation, two successive basis sets of the same type are used. Ideally, this should be at least triple- and quadruple-ζ variants with diffuse functions, which are often Dunning’s augmented correlation-consistent basis sets (aug-cc-pVTZ/aug-ccpVQZ 139 ). The problem with those is that their diffuse functions go to relatively high angular momenta which may not in fact be relevant for the REs and BHs studied herein, as can also be seen in Refs. 1 and 3. Ahlrichs-type basis sets, even though initially developed to be used with DFT methods, have also been used in the context of obtaining ab initio reference values, and in fact they may provide a faster alternative given that they also rely on fewer primitive Gaussian-type orbitals. Herein, we investigate the minimally augmented variants ma-def2-TZVPP and ma-def2-QZVPP. 140 Finally, we also examine ACS Paragon 15 Plus Environment


Page 16 of 62

Table 1: DLPNO-CCSD(T)/CBS barrier heights (BHs) and reaction energies (REs) in kcal/mol for all three sets of structures for models 0 of HKMT and AspDC. The triple- and quadruple-ζ variants of each basis set family were used for the extrapolation to the complete basis set limit. RIJCOSX was used for all calculations except those using the aug-cc-pVnZ basis sets. B3LYP/6-31G(d,p) PW6B95-D3(BJ)/def2-TZVP BH RE BH RE

Model

PBEh-3c BH RE

HKMT

def2-nZVPP 26.2 ma-def2-nZVPP 26.2 cc-pVnZ 25.8 aug-cc-pVnZ 25.9

−5.7 −5.6 −5.9 −5.9

27.3 27.2 26.9 26.9

−5.0 −5.0 −5.3 −5.3

26.8 26.8 26.4 26.4

−5.2 −5.2 −5.6 −5.6

AspDC

def2-nZVPP ma-def2-nZVPP cc-pVnZ aug-cc-pVnZ

−6.4 −6.5 −6.6 −6.3

3.9 3.9 3.7 3.7

−5.7 −5.8 −5.9 −5.7

4.3 4.2 4.1 3.9

−5.6 −5.6 −5.8 −5.6

3.3 3.2 2.9 3.0

their counterparts without diffuse functions, which further reduces computational cost, i.e. the combinations cc-pVTZ/cc-pVQZ 141 and def2-TZVPP/def2-QZVPP. 123 Our results for the aug-cc-pVTZ/aug-cc-pVQZ pair will serve as our best-possible benchmark. However, when considering the other basis-set pairs, we additionally include the RIJCOSX application to further speed up the HF part. This approximation only introduces marginal differences, as discussed in the SI (table S6). Results with all four basis-set pairs for the HKMT and AspDC model 0 are presented in table 1. When considering the results for the PBEh-3c structures, which are the structures used for the benchmark study conducted herein, it can be seen that there is minimal difference between the def2-nZVPP and ma-def2-nZVPP basis-set pairs, with identical results for the HKMT model 0 and a difference in the BH of only 0.1 kcal/mol in the AspDC model 0. A similar result is seen with the cc-pVnZ and aug-cc-pVnZ basis-set pairs, with these also giving identical values for the HKMT model 0 BH and RE, and a difference of 0.2 kcal/mol for each of the BH and RE for AspDC model 0. While this suggests that the basis sets without diffuse functions are adequate for these test systems, in ACS Paragon 16 Plus Environment

Page 17 of 62 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


practice the largest basis sets that can be afforded should be used to ensure that the most accurate and reliable result is achieved. As such, the aug-cc-pVTZ/aug-cc-pVQZ basis-set pair is chosen as the first choice for the extrapolation to the CBS limit when calculating reference values for the remaining systems due to its size. The diffuse functions in the augmented basis sets may also become also important when dealing with anionic species. Our second, more cost-efficient choice is the ma-def2-TZVPP/ma-def2-QZVPP basis-set pair, which are also augmented with diffuse functions.

The effect of the geometries on the reference values The results presented in the previous section also allow us to examine the effect of the level of theory used for the geometry optimizations on the energies obtained (table 1). Overall, it can be seen that for each model system and basis set, the BHs and REs obtained using the PBEh-3c and PW6B95-D3(BJ) structures are closer to each other than to the values obtained using the original B3LYP structures. This corresponds to what was seen when visually comparing the different geometries obtained (figures 2 and S1), where there were smaller differences between the PBEh-3c and PW6B95-D3(BJ) structures than between PBEh-3c and the original B3LYP ones. While it may be tempting to just use the original structures obtained with dispersionuncorrected DFT and to then get an accurate reference value for them—something that is often still done in applications—our results clearly demonstrate that there is an effect of the underlying structure on the reference values. BHs for the B3LYP structures of HKMT model 0 are between 0.5 and 1 kcal/mol lower than for the other structures of the same reaction, and in the reaction energies we observe changes of up to 0.7 kcal/mol. For AspDC model 0, we observe a similar difference with the B3LYP structures giving BHs at least 0.6 kcal/mol lower than those of the other structures, and REs that are at least 0.6 kcal/mol lower. The trends are very similar for all four basis-set types used for the ACS Paragon 17 Plus Environment


Page 18 of 62

CBS extrapolation and they reinforce our decision to reoptimize the structures rather than using the original ones for this work, as the errors in the geometry optimizations would carry through to the rest of the results.

DLPNO thresholds The DLPNO scheme involves first selecting strong and weak pair correlations through a two-step approach. The first step is the use of an efficient estimate of the pair correlation energy through multipole analysis to prescreen for negligible contributions. The pairs that pass the prescreening step are then used in a semicanonical local MP2 (SC-LMP2) calculation, or full MP2 with the TightPNO setting, and the pairs with pair correlation energy lower than a threshold value dubbed TCutP airs are removed. Pairs above that threshold are carried through to the coupled-cluster calculation. Pair natural orbitals (PNOs) are created based on the results of the SC-LMP2 calculation, but only the PNOs that have an occupation number larger than a threshold called TCutP N O are retained. A third parameter, TCutDo , is used to truncate the size of a given domain in which the PNOs are expanded, while a fourth parameter TCutM KN is tied to the other three main cut-offs. 112,134 Three default sets of the truncation parameters have been created by Neese and coworkers. 142 These are labelled LoosePNO, NormalPNO, and TightPNO and listed in table 2. The lower the value for each truncation parameter, the lower the number of pairs cut from the calculation. This makes the results closer to canonical CCSD(T), and thus more accurate, but also increases the cost of the calculation, as more pairs must be treated in Table 2: Definitions of the default thresholds for DLPNO-CCSD(T). 112,142 TCutP airs LoosePNO 10−3 NormalPNO 10−4 TightPNO 10−5

TCutDO

TCutP N O

TCutM KN

2 × 10−2 1 × 10−2 5 × 10−3

1.00 × 10−6 3.33 × 10−7 1.00 × 10−7

10−3 10−3 10−3


Page 19 of 62 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 coupled-cluster calculation. NormalPNO is the default setting for DLPNO-CCSD(T) calculations in ORCA, and it was reported that the DLPNO error for REs falls below 1 kcal/mol. 142 TightPNO was recommended to get accurate noncovalent interaction energies, while LoosePNO could be used for potential-energy surface scans or other prescreening procedures. 142 In our work, we aim to test all three settings in order to ensure that our reference values are converged while trying to maintain computational efficiency. Therefore calculations and subsequent CBS extrapolations were performed on HKMT model 0 for each PNO setting using our two preferred basis-set types for the CBS extrapolation, namely the aug-cc-pVTZ/aug-cc-pVQZ and ma-def2-TZVPP/ma-def2-QZVPP basis-set pairs. All of these tests were conducted without the RIJCOSX approximation.

The results of our tests are shown in table 3. The energies calculated with the TightPNO thresholds were taken as the standard to emulate, as these are the most accurate energies due to the very low degree of truncation. As such, LoosePNO can be immediately discounted as an option, as the results are far too inaccurate. For instance the aug-ccpVTZ/QZ-based BH increased by 2.3 kcal/mol and the RE increased by 0.7 kcal/mol compared to TightPNO results. For the Dunning basis sets, the NormalPNO results were not significantly different from the TightPNO ones, with slight increases of 0.4 and 0.1 kcal/mol, respectively. Interestingly, we observe that the accuracy of the various settings seems to also depend on the type of basis set, something that to our knowledge has not been reported before. In particular, we observe a difference of 1.4 kcal/mol for the BH and 0.9 kcal/mol for the RE for the Ahlrichs-type basis sets. Based on our tests it was decided that the TightPNO settings had to be used for all DLPNO-CCSD(T) calculations to maintain the highest level of accuracy possible.



Page 20 of 62

Table 3: DLPNO-CCSD(T)/CBS barrier heights (BHs) and reaction energies (REs) in kcal/mol obtained for HKMT model 0 with three different PNO thresholds and two different basis-set types. aug-cc-pVTZ/QZ BH RE LoosePNO 28.4 NormalPNO 26.5 TightPNO 26.1 a

−5.0 −5.6 −5.7

ma-def2-TZVP/QZVP BH RE —a 27.9 26.5

—a −4.6 −5.5

left out because the calculation with ma-def2-QZVPP gave highly incorrect

results with CBS BH and RE showing large deviations.

Composite Schemes When a straightforward CBS extrapolation turns out to be too expensive, it has become common practice to generate benchmark reference values based on cost-efficient composite schemes that approximate the actual CBS limit for a given approach. Very accurate component-based strategies to get close to a CCSD(T)/CBS limit are the W1-(F12) and W2-(F12) approaches, 143,144 but if these are too expensive, a faster variant that is often used in the literature is to obtain a CBS value for a lower-level (LL) method and to correct for the insufficient description of electron correlation by adding the difference in correlation energy between a higher level (HL) and the low level of theory with a smaller basis set. In general, the energy expression for the estimated (est.) CBS value, thus becomes: 145,146

E(est. HL/CBS) = E(est. LL/CBS) + EC (HL/small basis) − EC (LL/small basis)

(3)

where the terms are self-explanatory. Given that our goal is to establish protocols to obtain reliable reference values for larger systems, the HL methodology in this work is DLPNO-CCSD(T), while the tested LL methods are MP2, 147 SCS-MP2, 148 and DLPNO-SCS-MP2, 148,149 with the latter again


Page 21 of 62 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


Table 4: Testing various composite schemes to reproduce the reference barrier height (BH) and reaction energy (RE) for HKMT model 0 (kcal/mol). The aug-cc-pVTZ/aug-cc-pVQZ basis-set pair was used to extrapolate the LL energies, and the cc-pVDZ, aug-cc-pVDZ, cc-pVTZ and aug-cc-pVTZ basis sets were used for the HL-correction term between higher and lower level. The TightPNO settings were used. BH

RE

26.1

−5.7

25.9 25.6 26.3 26.3

−5.0 −5.1 −5.6 −5.6

25.8 25.4 26.3 26.2

−5.1 −5.2 −5.7 −5.7

correction/cc-pVDZ 25.8 correction/aug-cc-pVDZ 25.4 correction/cc-pVTZ 26.4 correction/aug-cc-pVTZ 26.2

−5.0 −5.2 −5.6 −5.7

Reference MP2/CBS MP2/CBS MP2/CBS MP2/CBS

+ + + +

correction/cc-pVDZ correction/aug-cc-pVDZ correction/cc-pVTZ correction/aug-cc-pVTZ

SCS-MP2/CBS SCS-MP2/CBS SCS-MP2/CBS SCS-MP2/CBS

+ + + +

correction/cc-pVDZ correction/aug-cc-pVDZ correction/cc-pVTZ correction/aug-cc-pVTZ

DLPNO-SCS-MP2/CBS DLPNO-SCS-MP2/CBS DLPNO-SCS-MP2/CBS DLPNO-SCS-MP2/CBS

+ + + +

being particularly useful for larger systems. The focus of our initial tests is investigating the effect of the different second-order perturbative methods on the final result, followed by establishing the minimum requirements for the small basis set used to correct for missing electron correlation. We use the aug-cc-pVTZ/aug-cc-pVQZ basis-set pair for the CBSextrapolation of the perturbative methods, and cc-pVDZ, aug-cc-pVDZ, cc-pVTZ and aug-cc-pVTZ for the HL correction term. The values that are ideally to be replicated are the DLPNO-CCSD(T)/CBS energies extrapolated with aug-cc-pVTZ/aug-cc-pVQZ and TightPNO thresholds, which were the most accurate values established in the DLPNOCCSD(T) test section. Table 4 includes all relevant details of our tests. A comparison between the three main types of schemes shows little difference between the underlying low-level method—in contradiction to what has been hinted at for DLPNO-



MP2/CBS numbers in Ref. 150—while the HL correction clearly shows a basis-set dependence that is often not elaborated on in the literature (table 4). In fact, it is often believed that the inherent difference in how much electron correlation is recovered is basis-set independent. 145,146 Our results clearly indicate that this is not the case. The double-ζ based correction leads in all cases to an underestimation of the BH, with aug-cc-pVDZ numbers showing an even larger deviation from the reference than cc-pVDZ. The two triple-ζ sets slightly overestimate the BH between 0.1 [(DLPNO-)SCS-MP2/CBS + HL-correction/augcc-pVTZ] and 0.3 kcal/mol (DLPNO-SCS-MP2/CBS + HL-correction/cc-pVTZ). The absolute deviation from the reference RE for the double-ζ-based corrections is up to 0.7 kcal/mol, whereas the two triple-ζ basis sets show no deviations or maximum deviations of only 0.1 kcal/mol. Our results so far suggest that our preferred approach based on best possible accuracy and cost-efficiency is to use DLPNO-SCS-MP2/CBS numbers and correct with a term based on cc-pVTZ. Note that those numbers are based on DLPNO calculations with TightPNO thresholds in both DLPNO-SCS-MP2 and DLPNO-CCSD(T). Next, we explore if the assessed composite schemes can also be used with less tight PNO thresholds; note that the thresholds defined for DLPNO-(SCS-)MP2 (table 5) are not the same as those used for DLPNO-CCSD(T). In addition, it needs to be investigated if also the more efficient Ahlrichs-type basis sets can be useful in composite schemes. Both types of analyses are combined in table 6, where for all DLPNO-SCS-MP2/CBS extrapolations the ma-def2TZVPP/ma-def2-QZVPP basis-set pair was used. Note that double-ζ basis sets for the correction term again turned out to be insufficient, which is why those results are only reported in the SI (tables S7-S9). It can be seen that the different DLPNO thresholds have no significant effect on the DLPNO-SCS-MP2/CBS extrapolations and that even LoosePNO settings provide the same results as TightPNO ones, despite TightPNO calculations being significantly more expenACS Paragon 22 Plus Environment

Page 22 of 62

Page 23 of 62 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


Table 5: Definitions of the default DLPNO thresholds for DLPNO-MP2 for closed-shell systems. 112,149 TCutDO LoosePNO 2 × 10−2 NormalPNO 1 × 10−2 TightPNO 5 × 10−3

TCutP N O 10−7 10−8 10−9

sive. Even when comparing the total energies in Hartree, the differing PNO thresholds only caused a change in the 5th or 6th decimal place, suggesting that the effect of the different thresholds on DLPNO-SCS-MP2 is not strong. As such, the choice of PNO setting is not important for this part of the procedure. On the other hand, for the calculation of the coupled-cluster term we see a significant dependence on the settings chosen for the DLPNO calculations. When using the LoosePNO settings for the correction term, we see large overestimation of the BH by about 2.5 kcal/mol, whereas the absolute RE is underestimated by 1 kcal/mol. While the RE seems to converge already for the NormalPNO setting (absolute deviation of 0.3 kcal/mol), the BH is lowered by another 0.5 kcal/mol compared to NormalPNO settings when the tightest thresholds are used. That being said, a comparison with our reference shows that the absolute differences are 0.3 kcal/mol for the NormalPNO option and 0.2 kcal/mol for the TightPNO. We can therefore conclude that the first is a viable solution if the latter setting poses computational-resource related problems. We also note that the results for the def2-TZVPP and ma-def2-TZVPP basis sets are practically identical and that therefore def2-TZVPP can be safely chosen to obtain the correction term. We conclude this section with an analysis of a second scheme that was suggested to approximate DLPNO-CCSD(T)/CBS results, namely the DLPNO-CCSD(T)/CBS* ap-



Page 24 of 62

Table 6: Composite scheme energies (kcal/mol) for HKMT model 0 using DLPNO-SCS-MP2/CBS with three different PNO thresholds, and two different AO basis sets for the HL correction term that was obtained with the three different PNO thresholds. The CBS-limit energies were obtained with the ma-def2-TZVPP/ma-def2-QZVPP basis-set pair. DLPNO-SCS-MP2 Correction Term Basis Correction Term PNO threshold def2-TZVPP LoosePNO ma-def2-TZVPP

def2-TZVPP NormalPNO ma-def2-TZVPP

def2-TZVPP TightPNO ma-def2-TZVPP

BH

RE

LoosePNO NormalPNO TightPNO LoosePNO NormalPNO TightPNO

28.6 26.4 25.9 28.5 26.4 25.9

−4.7 −5.4 −5.4 −4.7 −5.4 −5.4


28.6 26.4 25.9 28.5 26.4 25.9

−4.7 −5.4 −5.4 −4.7 −5.4 −5.4


28.5 26.4 25.9 28.5 26.4 25.9

−4.7 −5.4 −5.4 −4.7 −5.4 −5.4

26.1

−5.7

Reference proach, 38 which involves a multiplicative CBS extrapolation scheme of the form: DLPNO-CCSD(T)/def2-TZVP

E DLPNO-CCSD(T)/CBS* = E SCF/CBS(2,3) +

Ec

MP2/def2-TZVP Ec

! · fd · EcMP2/CBS(2,3)

(4) where fd = 1.08 is chosen to account for missing contributions from diffuse functions and CBS(2,3) denotes an extrapolation with the double and triple-ζ basis sets cc-pVDZ and cc-pVTZ. DLPNO-CCSD(T)/CBS* results have been used as reference values in recent bench-


Page 25 of 62 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


Table 7: DLPNO-CCSD(T)/CBS* results for HKMT model 0 (kcal/mol) obtained with the def2-SVP/def2-TZVP, ma-def2-SVP/ma-def2-TZVP, ccpVDZ/cc-pVTZ and aug-cc-pVDZ/aug-cc-pVTZ basis-set pairs. BH

RE

def2-SVP/def2-TZVP ma-def2-SVP/ma-def2-TZVP cc-pVDZ/cc-pVTZ aug-cc-pVDZ/aug-cc-pVTZ

22.5 23.6 24.8 25.4

−6.3 −5.8 −4.9 −5.4

Reference

26.1

−5.7

mark studies, including on other biochemical reactions, 38 and so it was of interest to test their accuracy on the enzymatically-catalyzed reactions of this benchmark set. Initially we tested this method as recommended, but after observing a severely underestimated BH and RE (table 7), further tests were conducted using the def2-SVP/def2-TZVP, madef2-SVP/ma-def2-TZVP and aug-cc-pVDZ/aug-cc-pVTZ basis-set pairs for the CBS(2,3) extrapolation terms in order to test if there was another viable alternative. The results with all four basis-set pairs are shown in table 7. When comparing the BHs in table 7 to the reference value of 26.1 kcal/mol, we see that the DLPNO-CCSD(T)/CBS* results significantly underestimate the barriers, with the recommended cc-pVDZ/cc-pVTZ basis-set pair giving a result that is 1.3 kcal/mol lower than the reference value. This basis-set pair also gives the worst result for the reaction energy, underestimating its absolute value by 0.8 kcal/mol. Of the additional basis-set pairs tested, the aug-cc-pVDZ/aug-cc-pVTZ basis-set pair gives the best result for the BH. However, even the best result is still 0.7 kcal/mol below the reference value, which is larger than for the other composite schemes. The REs calculated with the DLPNO-CCSD(T)/CBS* extrapolation scheme are closer to the −5.7 kcal/mol reference energy, particularly for the ma-def2-SVP/ma-def2-TZVP basis-set pair. However as all reference energies for a system must be calculated at the same level of theory, the consistent underestimation of the barrier



heights meant that the DLPNO-CCSD(T)/CBS* approach is not considered for reference energies for this work. We can thus conclude that for the present work the best strategy is a composite scheme of the form in Eq. 3. DLPNO-SCS-MP2 can be used as the lower level of theory and either extrapolated with the aug-cc-pVTZ/aug-cc-pVQZ basis set or the more efficient madef2-TZVPP/ma-def2-QZVPP pairs. A triple-ζ basis set without diffuse functions for the correction term is necessary (cc-pVTZ or def2-TZVPP). The PNO settings for the DLPNOSCS-MP2 part do not have to be as strict as for the DLPNO-CCSD(T) calculation, for which we recommend TightPNO when possible, and NormalPNO as an alternative. For these calculations, the RIJCOSX approximation can also be used to further decrease the cost of the calculation.

Preliminary benchmark set Reference values Based on the results of the DLPNO-CCSD(T) and composite scheme tests, we suggest a protocol for calculating reference values for related cases that can be adopted by others. In the present study, we apply our recommendations on a preliminary version of the updated benchmark set. All generated structures take into account London-dispersion effects and correct for potential problems arising from small AO basis sets; as outlined earlier, we made use of the PBEh-3c level of theory. For the calculation of benchmark REs as well as forward and reverse BHs, we suggest to carry out proper DLPNO-CCSD(T)/CBS extrapolations with TightPNO settings calculated with the aug-cc-pVTZ and aug-cc-PVQZ AO basis sets (strategy A). When this approach turns out to be too expensive, the second choice is DLPNO-CCSD(T)/CBS with TightPNO settings using the ma-def2-TZVPP and ma-def2-QZVPP basis sets for the extrapolation (strategy B). In this instance, we ACS Paragon 26 Plus Environment

Page 26 of 62

Page 27 of 62 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


also recommend the safe use of the RIJCOSX approximation to speed up the HF part. The third best option can be used when neither of the previous two strategies are practicable due to memory, disk space or other restrictions imposed by the high-performing computing system. In that instance, we recommend the composite scheme outlined in the previous section, namely DLPNO-SCS-MP2/CBS(ma-def2-TZVPP/ma-def2-QZVPP) using NormalPNO settings and an HL correction based on DLPNO-CCSD(T) and DLPNOSCS-MP2 with def2-TZVPP and TightPNO settings (strategy C). RIJCOSX can again be used for all HF calculations. The systems used in our small example benchmark study are listed in table 8 together with the strategy used to obtain reference values and the resulting REs, forward and reverse BHs. We intend to increase the number of data points continuously in the future.

Exemplified test of DFT methods The DFT methods used in this preliminary benchmark study were selected based upon popularity and recommendations based on various GMTKN55 studies. 3,5,6,12,174 The entire list of assessed DFAs, as well as the dispersion corrections and programs used for each DFA, is shown in table 9. The DFT-D3(BJ) correction was used for all DFAs except for ωB97X-D3(0), van der Waals (vdW) DFAs based on the VV10 176 kernel (suffix ‘-V’ or ‘NL’), and Minnesota DFAs (M06, M06L and M062X), for which the DFT-D3(0) correction is necessary. 1 Reasons behind using a dispersion correction for Minnesota functionals have been extensively provided in the literature. 1,3,7,9 Note that it has been established that vdW-DFT methods do not have to be treated in a fully-self-consistent way, and we follow the timely more efficient ‘post-SCF’ strategy, as explained in Ref. 6. The def2-QZVP basis set was used for all methods except for PBEh-3c, which uses its own special small basis set in conjunction with a BSSE correction, as stated earlier. The methods represent the four highest rungs in Perdew’s Jacob’s Ladder scheme 177 and we intend to determine ACS Paragon 27 Plus Environment


Page 28 of 62

Table 8: Summary of the obtained reference values (kcal/mol) and the strategies used to obtain them. Reaction Model Reference Strategy Reaction Stepa

BH

Reverse BH RE

AspDC

0 1 2

A B B

1 1 1

4.10 —b —

9.60 — —

HKMT

0 1 2

A A B

1 1 1

26.05 32.72 27.83 29.82 21.42 42.95

−5.67 1.99 −21.53

4-OT

0 0

B B

1 2

14.39 19.29 — —

−4.90 1.25

PTE

0 0

C C

1 2

6.48 9.44

4.51 10.90

1.97 −1.46

HheC

0 1

B C

1 1

— 17.64

— 3.58

16.05 14.06

−5.50 13.03 3.67

a

For AspDC, HKMT and HheC, which are one-step reactions, step 1 refers to the formation of the product. For 4-OT and PTE, which are two-step reactions, step 1 refers to the formation of the intermediate and step 2 refers to the formation of the product. b Left out of this preliminary study due to technical difficulties

whether our preliminary benchmark set already gives an indication if previously made recommendations can be applied to the systems included herein. 3,5,6,12 To analyze the results of the benchmark study, mean absolute deviations (MADs) are primarily used, as these indicate how the method performs against the reference values across all the models in the set. The MADs of each DFA with and without dispersion corrections calculated for all BHs and REs are shown in figure 3. Numerical values for all statistical quantities are also included in the SI (tables S13-S19). Our analysis herein combines BHs and REs to identify methods that are robust enough to follow an entire reaction pathway. Statistics for the individual BH and RE subsets are provided in the SI (tables S20 and S21).


Page 29 of 62 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


Table 9: List of DFAs used in the preliminary benchmark study. The dispersion corrections and programs used to run the calculations are also included alongside the references that determined dispersion-correction damping parameters for a given DFA. ‘ORCA*’ denotes that a local version of ORCA 4 was used. Name

Dispersion correction Program

Name

Dispersion correction Program

B97-D3(BJ) 109,128 BLYP 152–154 BP86 152,156,157 OLYP 153,154,158 PBE 126 revPBE 161 revPBE-NL 161,162

GGA D3(BJ) 109 D3(BJ) 109 D3(BJ) 109 D3(BJ) 1 D3(BJ) 109 D3(BJ) 109 VV10 162

Turbomole ORCA ORCA ORCA ORCA ORCA ORCA*

ωB97X-V 151 B3LYP 104,105 BHLYP 155 M06 75 M062X 75 PBE0 159,160 PBEh-3c 125 PW6B95 122

VV10 151 D3(BJ) 109 D3(BJ) 1 D3(0) 1 D3(0) 1 D3(BJ) 109 D3(BJ) 125 D3(BJ) 109

B97M-D3(BJ) 6,163 B97M-V 163 M06L 165 revTPSS 167 SCAN 116 TPSS 170

Meta-GGA D3(BJ) 6 VV10 163 D3(0) 1 D3(BJ) 3 D3(BJ) 117 D3(BJ) 109

ORCA* ORCA* ORCA ORCA Turbomole ORCA

ωB97M-D3(BJ) 6,173 ωB97M-V 173 ωB97X-D3(0) 175

Hybrid D3(BJ) 6 VV10 173 D3(0) 175

ORCA* ORCA* ORCA

Double Hybrid ωB97X-2 164 D3(BJ) 5 28 B2GP-PLYP D3(BJ) 1 166 B2NC-PLYP D3(BJ) 5 168 DSD-BLYP D3(BJ) 1 DSD-PBEP86 169 D3(BJ) 169 PWPB95 77 D3(BJ) 1 171 SOS0-PBE0-2 D3(BJ) 5 28,172 B2K-PLYP D3(BJ) 53 revDSD-PBEP86 174 D3(BJ) 174 revDOD-PBE 174 D3(BJ) 174 174 DOD-SCAN D3(BJ) 174

ORCA* ORCA ORCA ORCA ORCA ORCA ORCA ORCA

ORCA* ORCA ORCA ORCA ORCA ORCA ORCA ORCA ORCA ORCA ORCA

As many other studies have already shown the importance of the use of dispersion corrections, 1,3,13–20 we will only discuss the dispersion-uncorrected results briefly. When comparing the dispersion-corrected and uncorrected energies, it can be seen that the dispersion-corrected methods give lower MADs than the uncorrected ones. Only three DFAs show nearly negligible increase after adding a dispersion correction: revTPSS-D3(BJ) (0.06 kcal/mol), SCAN-D3(BJ) (0.19 kcal/mol), and B2NC-PLYP-D3(BJ) (0.01 kcal/mol). Note that when a dispersion-corrected method performs worse than its uncorrected counterpart, the dispersion correction is not necessarily to blame, but underlying problems in the DFA are revealed, as it may rely on error compensation when uncorrected which leads eventually to better-looking statistics. 3,54 This can be best seen for BHs and (meta-)GGAs and their mean deviations (MDs; see tables S20 and S21). revTPSS, for example, has an MD of −5.44 kcal/mol for BHs and −6.49 kcal/mol when corrected. This again demonACS Paragon 29 Plus Environment


Figure 3: Mean absolute deviations for all the functionals tested in the benchmark study with and without dispersion corrections. MADs are calculated for both barrier heights and reaction energies, and are given in kcal/mol. The vdW-DFT functionals and their D3(BJ) counterparts, as well as PBEh-3c, are presented only in their dispersion-corrected forms as their dispersion corrections are included in the functional itself. The results for the best functionals in each rung of Jacob’s Ladder are shown in red. strates the problem that a method that already underestimates a barrier will do so even more if the transition state is more stabilized by dispersion than the reactants. Coming back to the combined analysis of BHs and REs, we see there is a noticeable decrease in the MAD when using dispersion corrections (figure 3). The difference between the dispersion-corrected and uncorrected methods is very significant for OLYP and BHLYP, with reductions in the MADs of 3.45 and 1.92 kcal/mol, respectively—OLYP in particular gives the worst results of all the functionals studied when used without a dispersion correction, but OLYP-D3(BJ) is the best GGA functional (MAD=4.77 kcal/mol) and performs better than the two worst meta-GGA functionals, TPSS-D3(BJ) and revTPSSD3(BJ) (5.24 and 4.84 kcal/mol, respectively). Overall, our results reaffirm all the previous work into dispersion-corrected DFT, showing that DFAs should never be used without a dispersion correction.


Page 30 of 62

Page 31 of 62 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


Table 10: Best and worst dispersion-corrected DFAs in the benchmark study based on the MADs for the dispersion corrected methods. Results are given for each rung of Jacob’s Ladder. MAD values (kcal/mol) are shown in parentheses. Best Functional

Worst Functional

GGA OLYP-D3(BJ) (4.77) BLYP-D3(BJ) (6.50) Meta-GGA B97M-D3(BJ) (2.87) TPSS-D3(BJ) (5.24) Hybrid M062X-D3(0) (1.32) B3LYP-D3(BJ) (3.76) Double Hybrid SOS0-PBE0-2-D3(BJ) (1.07) DSD-BLYP-D3(BJ) (1.88) Now looking exclusively at the dispersion-corrected results, it is seen that GGA DFAs generally give the highest MADs, and the double hybrids give the lowest. This results in a downwards trend across the figure, which corresponds to the increasing accuracy that comes with each higher rung of Jacob’s Ladder. This is not a strict rule, however, as some methods outperform those in higher categories, such as B97M-D3(BJ) and B97M-V, two meta-GGA functionals that give results better than the two worst hybrids, B3LYP-D3(BJ) and M06-D3(0) (MADs of 2.87 and 2.97 kcal/mol vs 3.76 and 3.42 kcal/mol, respectively). M062X-D3(0) (MAD=1.32 kcal/mol), a hybrid, also outperforms many double hybrids except for SOS0-PBE0-2-D3(BJ) (MAD=1.07 kcal/mol), DOD-SCAN-D3(BJ) (MAD=1.22 kcal/mol), and revDOD-PBE-D3(BJ) (MAD=1.30 kcal/mol). MADs are, however, not the only measure of performance—the MDs and error spans can also be considered. When comparing M062X-D3(0) with the double hybrids, the first has an MD (−0.47 kcal/mol) which is not superior to any of the double hybrids, but comparable with those of B2NC-PLYPD3(BJ), and SOS0-PBE0-2-D3(BJ) (−0.46 and 0.51 kcal/mol, respectively). M062X-D3(0) is also outperformed by all the double hybrids except for ωB97X-2-D3(BJ) on error spans, with revDOD-PBE-D3(BJ) having the smallest error span (5.18 kcal/mol) compared to M062X-D3(0) (7.61 kcal/mol). Nevertheless, when cost is an important factor, M062XD3(0) is a reliable alternative that performs comparably well and we can comfortably recommend it.



Based on the MADs, the best and worst functionals for each rung of Jacob’s Ladder can be determined for this preliminary study. These results are given in table 10. Similar to what was seen in the original GMTKN55 study, OLYP-D3(BJ) performs very well. M062X-D3(0) also performs very well as mentioned previously; it was also the fourth best performing hybrid functional in the GMTKN55 study and therefore its performance here is not surprising. In Najibi and Goerigk’s work on van der Waals density functionals, 6 an updated list of the best functionals for the GMTKN55 benchmark set was developed—B97M-V was shown to be the best meta-GGA in the new list followed by its modification B97M-D3(BJ), with the latter being the best meta-GGA in our study (table 10). In fact, we observe that the two tested DFT-D3(BJ) variants of the B97M family of (hybrid-)meta-GGA DFAs perform as well as or slightly better than the original, more costly VV10-corrected versions. Our results for the double hybrids are close to previous recommendations, with ωB97X-2-D3(BJ) and B2NC-PLYP-D3(BJ) performing particularly well. 5 However the particularly good performance SOS0-PBE0-2-D3(BJ) was unexpected as it was only ninth-best double hybrid across the GMTKN55 database. 5,6 We can also recommend the fairly new and—with the correct code 174 —more time-efficient double-hybrids DOD-SCAN-D3(BJ) and revDOD-PBE-D3(BJ), which are the second and third best DFAs with MADs of 1.22 and 1.30 kcal/mol, respectively. So far, we have focussed mostly on DFT-D3 corrected results, as this correction is well established and accessible through a large number of quantum-chemical programs. However, for the interested reader we would also like to assess the impact of the new DFTD4 178,179 correction whose main difference to DFT-D3 is that its dispersion coefficients are atomic-charge dependent. Herein, we applied DFT-D4 with its ‘EEQ’ 180 parametrization and the three-body ‘ATM’-type 108 correction (see literature for more details). 178,179 We applied DFT-D4 to almost all tested DFAs with the exception of the B97 and B97M families and a handful of double hybrids for which the correction has not been parametrized ACS Paragon 32 Plus Environment

Page 32 of 62

Page 33 of 62 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


yet. We also applied the revDSDPBE-P86, DOD-SCAN and DOD-PBE versions that were specifically reparametrized for usage with DFT-D4. 174 The resulting statistical values for all methods are shown in tables S22-S24 in the SI. Overall, the new correction does not change any of our previous recommendations and most MADs are very close to the ones based on DFT-D3 with only minor absolute variations of 0.1 kcal/mol or less. The greatest observed difference is +0.20 kcal/mol for OLYP, followed by revDOD-PBE (+0.16 kcal/mol) and PW6B95 (−0.15 kcal/mol).

Effect of the reference values on the results The perceived performance of a DFA is entirely dependent on the reference values. Therefore, the level of theory used when calculating such reference data is of great importance, whether one is assessing the overall performance of a method or its performance relative to others. In order to show the actual accuracy and robustness of a DFA, the reference values must be calculated at the highest possible level of theory; for striking examples on how a ranking of DFAs changes with the quality of the reference, see Refs 3 and 53. In the original work on this benchmark set, 95 the reference values for the study had been based on the B3LYP/6-311+G(2d,2p) level of theory. As seen in our present study (figure 3), B3LYP performs particularly poorly for these systems, especially when used without a dispersion correction, and thus this level of theory should be completely inappropriate to be used as a reference value in a benchmark study, even when analyzing lower-level methods. To further emphasize the statement we made in the previous sentence, we present a brief analysis of GGA and meta-GGA DFAs benchmarked against B3LYP/6-311+G(2d,2p) reference values obtained for our new geometries (see table S25 for values). Because those alternative reference values were obtained without any dispersion corrections, the (meta)GGAs are also assessed without any such corrections, with the exception of vdW-DFAs ACS Paragon 33 Plus Environment


Page 34 of 62

Figure 4: Mean absolute deviations for the GGA and meta-GGA DFAs calculated for two sets of reference values: those detailed in table 8 and B3LYP/6-311+G(2d,2p) values. The MADs are calculated without dispersion corrections where possible for the B3LYP references and with dispersion corrections for the higher level references. The best DFA for each rung of Jacob’s Ladder is given in red (for the higher level references) or blue (for the B3LYP references). Table 11: Best and worst GGA and meta-GGA DFAs based on the MADs calculated with both the higher level reference values and B3LYP/6311+G(2d,2p) reference values. MAD values are given in kcal/mol. Higher level references Best Worst

B3LYP references Best Worst

GGA OLYP-D3(BJ) (4.77) BLYP-D3(BJ) (6.50) BP86 (2.49) revPBE-NL (5.09) Meta-GGA B97M-D3(BJ) (2.87) TPSS-D3(BJ) (5.24) TPSS (2.56) SCAN (3.65) or methods derived from them [rev-PBE and B97M-V/D3(BJ)]. Our results are presented in figure 4, along with the original MADs calculated with our higher level references. The best and worst DFAs for each set of reference values are also included in table 11. Note that the MD for B3LYP/6-311+G(2d,2p) compared to our high-level reference data is strongly negative (MD=−2.80 kcal/mol) and the MAD with 4.96 kcal/mol relatively high, which is an early indication that that level of theory should never be used as a reference to benchmark low-level methods. The most obvious and trivial result of this comparison is that the MADs for GGAs and meta-GGAs are (significantly) lower when based on the B3LYP references rather than ACS Paragon 34 Plus Environment

Page 35 of 62 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


on our high-level ones. This is likely due to the fact that any errors caused by the lack of dispersion correction in the DFA affect the reference method too, and thus there is a high potential for error cancellation that artificially improves the results. This also explains the poor performance of revPBE-NL with respect to the B3LYP references, as any error cancellation due to missing dispersion forces cannot occur with methods that are inherently dispersion-corrected. Another important result is that the best and worst DFAs are not consistent, with the greatest difference being observed for TPSS. When calculated without any dispersion correction, TPSS is the best performing meta-GGA compared with B3LYP references (MAD=2.56 kcal/mol), but when using the better reference values and dispersion corrections, TPSS-D3(BJ) becomes the worst meta-GGA, with an MAD of 5.24 kcal/mol. We, thus, demonstrated once more the importance of choosing the right level of theory as a benchmark for the assessment of lower-level methods.

Summary and Conclusions With ever-improving computational resources, benchmark sets need to be developed that allow the assessment of low-level methods—usually DFT approaches—for systems of chemically relevant size. Herein we continued work on a benchmark set on enzymatically catalyzed reactions proposed in Ref. 95. We presented an updated version of this set and assessed DFT methods for the first time. Whilst updating the test set, we first outlined the importance of using an adequate level of theory for geometry optimizations that properly addresses London dispersion and basis-set superposition error, as these two effects influence structural features and, as such, the resulting reaction barrier heights and energies in the final benchmark set. We could, thus, conclude that the still very popular B3LYP/6-31G(d,p) level of theory for geometry optimizations, on which the original set



was based, is inadequate, adding to an ever growing list of evidence and calls to not rely on it as a black-box method for geometry optimizations. 11,17,19 Next, we proceeded with a detailed analysis on how to obtain DLPNO-CCSD(T)/CBS reference values with the best-possible accuracy and investigated different basis sets for the extrapolation to the complete-basis-set (CBS) limit, different technical thresholds of the DLPNO approach and various composite approaches that can be used if a conventional extrapolation is technically too demanding. From these tests we developed a protocol comprising three strategies that can be applied to systems of similar or larger size. DLPNOCCSD(T) should ideally be chosen with the “TightPNO” setting and if possible we advocated to carry out two-point extrapolations to the CBS limit with the aug-cc-pVTZ/QZ basis-set pair. Timely more efficient, with almost similar accuracy, is an alternative extrapolation with the ma-def2-TZ/QZVPP basis-set pair and the RIJCOSX approximation to speed up HF calculations. The best and most-efficient composite approach that can be used when neither of the first two strategies are feasible is to obtain DLPNO-SCS-MP2/CBS numbers with the NormalPNO setting and the ma-def2-TZ/QZVPP basis-set pair followed by a higher-level correction that represents the difference between DLPNO-CCSD(T) and DLPNO-SCS-MP2 with TightPNO settings and the def2-TZVPP basis set. We note that using a double-ζ basis set for the higher-level correction turned out to be insufficient, indicating that there is indeed a basis-set dependence for such corrections contrary to popular belief. While our main focus has been to establish strategies for the generation of accurate reference values that others can use in their future works, we also apply our recommendations and presented data for a preliminary benchmark set that we used in an assessment of 35 selected DFT methods. Our preliminary test set contains 10 model systems for enzymatically catalyzed reactions, which range in size from 27 to 112 atoms, in which we analyzed 16 barrier heights and 12 reaction energies. The validity of our preliminary benchmark set and our ACS Paragon 36 Plus Environment

Page 36 of 62

Page 37 of 62 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


protocol to generate reference values can be seen from the fact that our findings generally agree with more detailed comprehensive benchmark studies, 3,5,6,12 namely that Londondispersion corrections always have to be taken into consideration in DFT treatments of such systems and that the Jacob’s Ladder scheme is generally reproduced with double hybrids being generally more accurate than hybrids, followed by meta-GGA and GGA methods. However, taken into consideration the well-known self-interaction-error problem for barrier heights, we can particularly recommend the SOS0-PBE0-2-D3(BJ), DOD-SCAN-D3(BJ), revDOD-PBE-D3(BJ), and M062X-D3(0) approaches, as well as other standard doublehybrid DFAs. We also took the opportunity to shed light on the importance of using the best-possible reference values to assess lower-level methods with dramatic changes in rankings and recommendations when an inadequate level of theory is chosen as the benchmark. In particular, such an inadequate level would be the popular B3LYP DFA in conjunction with a triple-ζ basis set, rendering that approach obsolete in both the generation of benchmark data and applications of a similar kind. We are confident that with our findings we have done the next important step toward the establishment of an accurate and representative set of benchmark data for biologically relevant reactions. In the future, our recommended guidelines could not only be used to further enhance our preliminary benchmark set, but we would like to encourage the reader to apply our guidelines in their own work.

Acknowledgement LG is deeply grateful to Prof. Leo Radom for his continuing support and mentorship over the past 12 years. We thank Marina Jansen for performing some preliminary calculations on the AspDC system. We also thank for the allocation of computing resources by Research Platform Services (ResPlat) at The University of Melbourne (Project ID:



punim0094), Melbourne Bioinformatics (Project ID: RA0005), and the National Computational Infrastructure (NCI) National Facility within the National Computational Merit Allocation Scheme (Project ID: fk5). DAW acknowledges an Australian Government Research Training Program Scholarship.

Notes The authors declare no competing financial interest.

Supporting Information Available Information on the geometry optimizations and xyz files of all PBEh-3c optimized structures; further details on the determination of reference values; detailed information on the results for every tested DFT method.

References (1) Goerigk, L.; Grimme, S. A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions. Phys. Chem. Chem. Phys. 2011, 13, 6670–6688. (2) Mardirossian, N.; Head-Gordon, M. Thirty years of density functional theory in computational chemistry: An overview and extensive assessment of 200 density functionals. Mol. Phys. 2017, 115, 2315–2372. (3) Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Grimme, S. A look at the density functional theory zoo with the advanced GMTKN55 database for


Page 38 of 62

Page 39 of 62 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


general main group thermochemistry, kinetics and noncovalent interactions. Phys. Chem. Chem. Phys. 2017, 19, 32184–32215. (4) Mardirossian, N.; Head-Gordon, M. Survival of the most transferable at the top of Jacob’s ladder: Defining and testing the ωB97M(2) double hybrid density functional. J. Chem. Phys. 2018, 148, 241736. (5) Mehta, N.; Casanova-Páez, M.; Goerigk, L. Semi-empirical or non-empirical doublehybrid density functionals: Which are more robust? Phys. Chem. Chem. Phys. 2018, 20, 23175–23194. (6) Najibi, A.; Goerigk, L. The non-local kernel in van-der-Waals density functionals as an additive correction — an extensive analysis with special emphasis on the B97M-V and ωB97M-V approaches. J. Chem. Theory Comput. 2018, 14, 5725–5738. (7) Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking density functional methods against the S66 and S66x8 datasets for non-covalent interactions. Chem. Phys. Chem. 2011, 12, 3421–3433. (8) Goerigk, L. How do DFT-DCP, DFT-NL, and DFT-D3 compare for the description of London-dispersion effects in conformers and general thermochemistry? J. Chem. Theory Comput 2014, 10, 968–980. (9) Goerigk, L. Treating London-dispersion effects with the latest Minnesota density functionals: Problems and possible solutions. J. Phys. Chem. Lett. 2015, 6, 3891– 3896. (10) Mardirossian, N.; Head-Gordon, M. How accurate are the Minnesota density functionals for noncovalent interactions, isomerization energies, thermochemistry, and barrier heights involving molecules composed of main-group elements? J. Chem. Theory Comput. 2016, 12, 4303–4325. ACS Paragon 39 Plus Environment


(11) Kruse, H.; Goerigk, L.; Grimme, S. Why the standard B3LYP/6-31G* model chemistry should not be used in DFT calculations of molecular thermochemistry: Understanding and correcting the problem. J. Org. Chem. 2012, 77, 10824–10834. (12) Goerigk, L.; Mehta, N. A trip to the Density Functional Theory zoo: Warnings and recommendations for the user. Aust. J. Chem. 2019, published online, DOI: 10.1071/CH19023. ˇ ˇ y, J.; Hobza, P. Benchmark database of accurate (MP2 (13) Jureˇcka, P.; Sponer, J.; Cern´ and CCSD(T) complete basis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985–1993. (14) Antony, J.; Grimme, S. Density functional theory including dispersion corrections for intermolecular interactions in a large benchmark set of biologically relevant molecules. Phys. Chem. Chem. Phys. 2006, 8, 5287–5293. ˇ aˇc, J.; Riley, K. E.; Hobza, P. S66: A well-balanced database of benchmark (15) Rez´ interaction energies relevant to biomolecular structures. J. Chem. Theory Comput 2011, 7, 2427–2438. (16) Grimme, S.; Huenerbein, R.; Ehrlich, S. On the importance of the dispersion energy for the thermodynamic stability of molecules. ChemPhysChem 2011, 12, 1258–1261. (17) Goerigk, L.; Reimers, J. R. Efficient methods for the quantum chemical treatment of protein structures: The effects of London-dispersion and basis-set incompleteness on peptide and water-cluster geometries. J. Chem. Theory Comput. 2013, 9, 3240– 3251. (18) Grimme, S.; Steinmetz, M. Effects of London dispersion correction in density func-


Page 40 of 62

Page 41 of 62 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


tional theory on the structures of organic molecules in the gas phase. Phys. Chem. Chem. Phys. 2013, 15, 16031–16042. (19) Goerigk, L.; Collyer, C. A.; Reimers, J. R. Recommending Hartree-Fock theory with London-dispersion and basis-set-superposition corrections for the optimization or quantum refinement of protein structures. J. Phys. Chem. B 2014, 118, 14612– 14626. (20) Reimers, J. R.; Panduwinata, D.; Visser, J.; Chin, Y.; Tang, C.; Goerigk, L.; Ford, M. J.; Sintic, M.; Sum, T.-J.; Coenen, M. J. J. et al. A priori calculations of the free energy of formation from solution of polymorphic self-assembled monolayers. Proc. Natl. Acad. Sci. USA 2015, 112, E6101–E6110. (21) Grimme, S.; Schreiner, P. R. Steric crowding can stabilize a labile molecule: Solving the hexaphenylethane riddle. Angew. Chem. Int. Ed. 2011, 50, 12639–12642. (22) Wagner, J. P.; Schreiner, P. R. London dispersion in molecular chemistry— reconsidering steric effects. Angew. Chem. Int. Ed. 2015, 54, 12274–12296. (23) Rösel, S.; Quanz, H.; Logemann, C.; Becker, J.; Mossou, E.; Ca˜ nadillas-Delgado, L.; Caldeweyher, E.; Grimme, S.; Schreiner, P. R. London dispersion enables the shortest intermolecular hydrocarbon H···H contact. J. Amer. Chem. Soc. 2017, 139, 7428– 7431. (24) Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. Gaussian-2 theory for molecular energies of first- and second-row compounds. J. Chem. Phys. 1991, 94, 7221–7230. (25) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. J. Chem. Phys. 1997, 106, 1063–1079. ACS Paragon 41 Plus Environment


(26) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. Gaussian3 (G3) theory for molecules containing first and second-row atoms. J. Chem. Phys. 1998, 109, 7764–7776. (27) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Assessment of Gaussian-3 and density-functional theories on the G3/05 test set of experimental energies. J. Chem. Phys. 2005, 123, 124107. (28) Karton, A.; Tarnopolsky, A.; Lame`re, J.-F.; Schatz, G. C.; Martin, J. M. L. Highly accurate first-principles benchmark data sets for the parametrization and validation of density functional and other approximate methods. Derivation of a robust, generally applicable, double-hybrid functional for thermochemistry and thermochemical kinetics. J. Phys. Chem. A 2008, 112, 12868–12886. (29) Karton, A.; Daon, S.; Martin, J. M. L.; Ruscic, B. W4-11: A high-confidence benchmark dataset for computational thermochemistry derived from first-principles W4 data. Chem. Phys. Lett. 2011, 510, 165–178. (30) Karton, A.; Sylvetsky, N.; Martin, J. M. L. W4-17: A diverse and high-confidence dataset of atomization energies for benchmarking high-level electronic structure methods. J. Comp. Chem. 2017, 38, 2063–2075. (31) Rezac, J.; Jurecka, P.; Riley, K. E.; Cerny, J.; Valdes, H.; Pluhackova, K.; Berka, K.; Rezac, T.; Pitonak, M.; Vondrasek, J. et al. Quantum chemical benchmark energy and geometry database for molecular clusters and complex molecular systems (www.begdb.com): A users manual and examples. Collect. Czech. Chem. Commun. 2008, 73, 1261–1270. ˇ (32) Reha, D.; Valdes, H.; Vondrasek, J.; Hobza, P.; Abu-Riziq, A.; Crews, B.; de Vries, M. S. Structure and IR spectrum of phenylalanyl-glycyl-glycine tripetide ACS Paragon 42 Plus Environment

Page 42 of 62

Page 43 of 62 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


in the gas-phase: IR/UV experiments, ab initio quantum chemical calculations, and molecular dynamic simulations. Chem. Eur. J. 2005, 11, 6803–6817. (33) Gruzman, D.; Karton, A.; Martin, J. M. L. Performance of ab initio and density functional methods for conformational equilibria of Cn H2n+2 Alkane Isomers (n= 4-8). J. Phys. Chem. A 2009, 113, 11974–11983. (34) Kesharwani, M. K.; Karton, A.; Martin, J. M. L. Benchmark ab initio conformational energies for the proteinogenic amino acids through explicitly correlated methods. assessment of density functional methods. J. Chem. Theory Comput. 2016, 12, 444– 454. (35) Fogueri, U. R.; Kozuch, S.; Karton, A.; Martin, J. M. L. The melatonin conformer space: Benchmark and assessment of wave function and DFT methods for a paradigmatic biological and pharmacological molecule. J. Phys. Chem. A 2013, 117, 2269– 2277. (36) Lao, K. U.; Schäffer, R.; Jansen, G.; Herbert, J. M. Accurate description of intermolecular interactions involving ions using symmetry-adapted perturbation theory. J. Chem. Theory Comput. 2015, 11, 2473–2486. (37) Kozuch, S.; Bachrach, S. M.; Martin, J. M. Conformational equilibria in butane-1,4diol: A benchmark of a prototypical system with strong intramolecular H-bonds. J. Phys. Chem. A 2014, 118, 293–303. ˇ (38) Kruse, H.; Mladek, A.; Gkionis, K.; Hansen, A.; Grimme, S.; Sponer, J. Quantum chemical benchmark study on 46 RNA backbone families using a dinucleotide unit. J. Chem. Theory Comput. 2015, 11, 4972–4991. (39) Kozuch, S.; Martin, J. M. L. Halogen bonds: Benchmarks and theoretical analysis. J. Chem. Theory Comput. 2013, 9, 1918–1931. ACS Paragon 43 Plus Environment


(40) Goerigk, L.; Karton, A.; Martin, J. M. L.; Radom, L. Accurate quantum chemical energies for tetrapeptide conformations: Why MP2 data with an insufficient basis set should be handled with caution. Phys. Chem. Chem. Phys. 2013, 15, 7028–7031. (41) Beckett, D.; El-Baba, T. J.; Clemmer, D. E.; Raghavachari, K. Electronic energies are not enough: An ion mobility-aided, quantum chemical benchmark analysis of H+GPGG conformers. J. Chem. Theory Comput. 2018, 14, 5406–5418. (42) Prasad, V. K.; de-la Roza, A. O.; DiLabio, G. A. PEPCONF, a diverse data set of peptide conformational energies. Sci. Data 2019, 6, 180310. (43) Grimme, S.; Steinmetz, M.; Korth, M. How to compute isomerization energies of organic molecules with quantum chemical methods. J. Org. Chem. 2007, 72, 2118– 2126. (44) Margraf, J. T.; Ranasinghe, D. S.; Bartlett, R. J. Automatic generation of reaction energy databases from highly accurate atomization energy benchmark sets. Phys. Chem. Chem. Phys. 2017, 19, 9798–9805. (45) Yu, H.; Truhlar, D. G. Components of the bond energy in polar diatomic molecules, radicals, and ions formed by group-1 and group-2 metal atoms. J. Chem. Theory Comput. 2015, 11, 2968–2983. (46) Zhao, Y.; Ng, H. T.; Peverati, R.; Truhlar, D. G. Benchmark database for ylidic bond dissociation energies and its use for assessments of electronic structure methods. J. Chem. Theory Comput. 2012, 8, 2824–2834. (47) Yu, L.-J.; Karton, A. Assessment of theoretical procedures for a diverse set of isomerization reactions involving double-bond migration in conjugated dienes. Chem. Phys. 2014, 441, 166 – 177.


Page 44 of 62

Page 45 of 62 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


(48) Karton, A.; Martin, J. M. Explicitly correlated benchmark calculations on C8H8 isomer energy separations: How accurate are DFT, double-hybrid, and composite ab initio procedures? Mol. Phys. 2012, 110, 2477–2491. (49) Friedrich, J.; Hänchen, J. Incremental CCSD(T)(F12*)—MP2: A black box method to obtain highly accurate reaction energies. J. Chem. Theory Comput. 2013, 9, 5381–5394. (50) Huenerbein, R.; Schirmer, B.; Moellmann, J.; Grimme, S. Effects of London dispersion on the isomerization reactions of large organic molecules: A density functional benchmark study. Phys. Chem. Chem. Phys. 2010, 12, 6940–6948. (51) Zhao, Y.; González-Garc´ıa, N.; Truhlar, D. G. Benchmark database of barrier heights for heavy atom transfer, nucleophilic substitution, association, and unimolecular reactions and its use to test theoretical methods. J. Phys. Chem. A 2005, 109, 2012–2018. (52) Zhao, Y.; Lynch, B. J.; Truhlar, D. G. Multi-coefficient extrapolated density functional theory for thermochemistry and thermochemical kinetics. Phys. Chem. Chem. Phys. 2005, 17, 43–52. (53) Karton, A.; Goerigk, L. Accurate reaction barrier heights of pericyclic reactions: Surprisingly large deviations for the CBS-QB3 composite method and their consequences in DFT benchmark studies. J. Comput. Chem. 2015, 36, 622–632. (54) Goerigk, L.; Sharma, R. The INV24 test set: How well do quantum-chemical methods describe inversion and racemization barriers? Can. J. Chem. 2016, 94, 1133– 1143. (55) Karton, A.; O’Reilly, R. J.; Chan, B.; Radom, L. Determination of barrier heights for proton exchange in small water, ammonia, and hydrogen fluoride clusters with ACS Paragon 45 Plus Environment


G4(MP2)-Type, MPn, and SCS-MPn procedures: A caveat. J. Chem. Theory Comput. 2012, 8, 3128–3136. (56) Scott, A. P.; Radom, L. Harmonic vibrational frequencies: An evaluation of HartreeFock, Møller-Plesset, quadratic configuration interaction, density functional theory, and semiempirical scale factors. J. Phys. Chem. 1996, 100, 16502–16513. (57) Merrick, J. P.; Moran, D.; Radom, L. An evaluation of harmonic vibrational frequency scale factors. J. Phys. Chem. A 2007, 111, 11683–11700. (58) Biczysko, M.; Panek, P.; Scalmani, G.; Bloino, J.; Barone, V. Harmonic and anharmonic vibrational frequency calculations with the double-hybrid B2PLYP method: Analytic second derivatives and benchmark studies. J. Chem. Theory Comput. 2010, 6, 2115–2125. (59) Puzzarini, C.; Biczysko, M.; Barone, V. Accurate harmonic/anharmonic vibrational frequencies for open-shell systems: Performances of the B3LYP/N07D model for semirigid free radicals benchmarked by CCSD(T) computations. J. Chem. Theory Comput. 2010, 6, 828–838. (60) Bloino, J.; Biczysko, M.; Barone, V. General perturbative approach for spectroscopy, thermodynamics, and kinetics: Methodological background and benchmark studies. J. Chem. Theory Comput. 2012, 8, 1015–1036. (61) Chan, B.; Radom, L. Frequency scale factors for some double-hybrid density functional theory procedures: Accurate thermochemical components for high-level composite protocols. J. Chem. Theory Comput. 2016, 12, 3774–3780. (62) Wilson, D. J. D.; Mohn, C. E.; Helgaker, T. The rotational g tensor as a benchmark for density-functional theory calculations of molecular magnetic properties. J. Chem. Theory Comput. 2005, 1, 877–888. ACS Paragon 46 Plus Environment

Page 46 of 62

Page 47 of 62 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


(63) Lutnæs, O. B.; Teale, A. M.; Helgaker, T.; Tozer, D. J.; Ruud, K.; Gauss, J. Benchmarking density-functional-theory calculations of rotational g tensors and magnetizabilities using accurate coupled-cluster calculations. J. Chem. Phys. 2009, 131, 144104. (64) Teale, A. M.; Lutnæs, O. B.; Helgaker, T.; Tozer, D. J.; Gauss, J. Benchmarking density-functional theory calculations of NMR shielding constants and spin–rotation constants using accurate coupled-cluster calculations. J. Chem. Phys. 2013, 138, 024111. (65) Diedrich, C.; Grimme, S. Systematic investigation of modern quantum chemical methods to predict electronic circular dichroism spectra. J. Phys. Chem. A 2003, 107, 2524–2539. (66) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, 134110. (67) Goerigk, L.; Grimme, S. Calculation of electronic circular dichroism spectra with time-dependent double-hybrid density functional theory. J. Phys. Chem. A 2009, 113, 767–776. (68) Goerigk, L.; Moellmann, J.; Grimme, S. Computation of accurate excitation energies for large organic molecules with double-hybrid density functionals. Phys. Chem. Chem. Phys. 2009, 11, 4611–4620. (69) Goerigk, L.; Grimme, S. Assessment of TD-DFT methods and of various spin scaled CIS(D) and CC2 versions for the treatment of low-lying valence excitations of large organic dyes. J. Chem. Phys. 2010, 132, 184103.



(70) Leang, S. S.; Zahariev, F.; Gordon, M. S. Benchmarking the performance of timedependent density functional methods. J. Chem. Phys. 2012, 136, 104101. (71) Jacquemin, D.; Duchemin, I.; Blase, X. Is the Bethe-Salpeter formalism accurate for excitation energies? Comparisons with TD-DFT, CASPT2, and EOM-CCSD. J. Phys. Chem. Lett. 2017, 8, 1524–1529. (72) Schwabe, T.; Goerigk, L. Time-dependent double-hybrid density functionals with spin-component and spin-opposite scaling. J. Chem. Theory Comput. 2017, 13, 4307–4323. (73) Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. A mountaineering strategy to excited states: Highly accurate reference energies and benchmarks. J. Chem. Theory Comput. 2018, 14, 4360–4379. (74) Casanova-Páez, M.; Dardis, M. B.; Goerigk, L. ωB2PLYP & ωB2GPPLYP: The first two double-hybrid density functionals with long-range correction optimized for excitation energies. J. Chem. Theory Comput. 2019, published online, DOI: 10.1021/acs.jctc.9b00013. (75) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215–241. (76) Goerigk, L.; Grimme, S. A general database for main group thermochemistry, kinetics, and noncovalent interactions — Assessment of common and reparameterized (meta-)GGA density functionals. J. Chem. Theory Comput. 2010, 6, 107–126. (77) Goerigk, L.; Grimme, S. Efficient and accurate double-hybrid-meta-GGA density functionals — Evaluation with the extended GMTKN30 database for general main ACS Paragon 48 Plus Environment

Page 48 of 62

Page 49 of 62 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


group thermochemistry, kinetics and noncovalent interactions. J. Chem. Theory Comput. 2011, 7, 291–309. (78) Peverati, R.; Truhlar, D. G. Quest for a universal density functional: The accuracy of density functionals across a broad spectrum of databases in chemistry and physics. Philos. Trans. R. Soc., A 2014, 372, 20120476. (79) Haoyu, S. Y.; He, X.; Li, S. L.; Truhlar, D. G. MN15: A Kohn–Sham global-hybrid exchange–correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions. Chem. Sci. 2016, 7, 5032– 5051. (80) Morgante, P.; Peverati, R. ACCDB: A collection of chemistry databases for broad computational purposes. J. Comput. Chem. 2019, 40, 839–848. (81) O’Reilly, R. J.; Karton, A.; Radom, L. N–H and N–Cl homolytic bond dissociation energies and radical stabilization energies: An assessment of theoretical procedures through comparison with benchmark-quality W2w data. Int. J. Quantum Chem. 2011, 112, 1862–1878. (82) Parkinson, C. J.; Mayer, P. M.; Radom, L. Cyanovinyl radical: An illustration of the poor performance of unrestricted perturbation theory and density functional theory procedures in calculating radical stabilization energies. Theor. Chem. Acc. 1999, 102, 92–96. (83) Henry, D. J.; Parkinson, C. J.; Mayer, P. M.; Radom, L. Bond dissociation energies and radical stabilization energies associated with substituted methyl radicals. J. Phys. Chem. A 2001, 105, 6750–6756. (84) Wood, G. P. F.; Henry, D. J.; Radom, L. Performance of the RB3-LYP, RMP2,



and UCCSD(T) procedures in calculating radical stabilization energies for •NHX radicals. J. Phys. Chem. A 2003, 107, 7985–7990. (85) Wood, G. P. F.; Moran, D.; Jacob, R.; Radom, L. Bond dissociation energies and radical stabilization energies associated with model peptide-backbone radicals. J. Phys. Chem. A 2005, 109, 6318–6325. (86) Chan, B.; Radom, L. BDE261: A comprehensive set of high-level theoretical bond dissociation enthalpies. J. Phys. Chem. A 2012, 116, 4975–4986. (87) Mayer, P. M.; Parkinson, C. J.; Smith, D. M.; Radom, L. An assessment of theoretical procedures for the calculation of reliable free radical thermochemistry: A recommended new procedure. J. Chem. Phys. 1998, 108, 604–615. (88) Graham, D. C.; Menon, A. S.; Goerigk, L.; Grimme, S.; Radom, L. Optimization and basis-set dependence of a restricted-open-shell form of B2-PLYP double-hybrid density functional theory. J. Phys. Chem. A 2009, 113, 9861–9873. (89) Chan, B.; Deng, J.; Radom, L. G4(MP2)-6X: A cost-effective improvement to G4(MP2). J. Chem. Theory Comput. 2011, 7, 112–120. (90) Henry, D. J.; Sullivan, M. B.; Radom, L. G3-RAD and G3X-RAD: Modified Gaussian-3 (G3) and Gaussian-3X (G3X) procedures for radical thermochemistry. J. Chem. Phys. 2003, 118, 4849–4860. (91) Chan, B.; Coote, M. L.; Radom, L. G4-SP, G4(MP2)-SP, G4-sc, and G4(MP2)-sc: Modifications to G4 and G4(MP2) for the treatment of medium-sized radicals. J. Chem. Theory Comput. 2010, 6, 2647–2653. (92) Karton, A.; O’Reilly, R. J.; Radom, L. Assessment of theoretical procedures for cal-


Page 50 of 62

Page 51 of 62 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


culating barrier heights for a diverse set of water-catalyzed proton-transfer reactions. J. Phys. Chem. A 2012, 116, 4211–4221. (93) Chan, B.; Radom, L. Assessment of theoretical procedures for hydrogen-atom abstraction by chlorine, and related reactions. Theor. Chem. Acc. 2011, 130, 251–260. ˇ (94) Szabla, R.; Havrila, M.; Kruse, H.; Sponer, J. Comparative assessment of different RNA tetranucleotides from the DFT-D3 and force field perspective. J. Phys. Chem. B 2016, 120, 10635–10648. (95) Kromann, J. C.; Christensen, A. S.; Cui, Q.; Jensen, J. H. Towards a barrier height benchmark set for biologically relevant systems. PeerJ 2016, 4, e1994. (96) Liao, R.-Z.; Yu, J.-G.; Himo, F. Quantum chemical modeling of enzymatic reactions: The case of decarboxylation. J. Chem. Theory Comput. 2011, 7, 1494–1501. (97) Georgieva, P.; Himo, F. Quantum chemical modeling of enzymatic reactions: The case of histone lysine methyltransferase. J. Comput. Chem. 2010, 31, 1707–1714. (98) Sevastik, R.; Himo, F. Quantum chemical modeling of enzymatic reactions: The case of 4-oxalocrotonate tautomerase. Bioorg. Chem. 2007, 35, 444–457. (99) Chen, S.-L.; Fang, W.-H.; Himo, F. Theoretical study of the phosphotriesterase reaction mechanism. J. Phys. Chem. B 2007, 111, 1253–1255. (100) Hopmann, K. H.; Himo, F. Quantum chemical modeling of the dehalogenation reaction of haloalcohol dehalogenase. J. Chem. Theory Comput. 2008, 4, 1129–1137. (101) Stewart, J. J. P. Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173–1213.



(102) Stewart, J. J. P. Optimization of parameters for semiempirical methods VI: More modifications to the NDDO approximations and re-optimization of parameters. J. Mol. Model. 2012, 19, 1–32. (103) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the self-consistent-charge density-functional tight-binding method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931–948. (104) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (105) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 1994, 98, 11623–11627. (106) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-consistent molecular orbital methods. XII. Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules. J. Chem. Phys. 1972, 56, 2257–2261. (107) Hehre, W. J. Ab initio molecular orbital theory. Acc. Chem. Res. 1976, 9, 399–406. (108) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (109) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (110) Kruse, H.; Grimme, S. A geometrical correction for the inter- and intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems. J. Chem. Phys. 2012, 136, 04B613. ACS Paragon 52 Plus Environment

Page 52 of 62

Page 53 of 62 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


(111) Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 2, 73–78. (112) Neese, F. ORCA, versions 4.0.0 and 4.1.2; , MPI f¨ ur Kohleforschung: M¨ ulheim a. d. Ruhr, Germany, 2017. (113) TURBOMOLE, version 7.1.1; COSMOlogic: Leverkusen, Germany, 2016. (114) Häser, M.; Ahlrichs, R. Improvements on the direct SCF method. J. Comput. Chem. 1989, 10, 104–111. (115) Bauernschmitt, R.; Ahlrichs, R. Stability analysis for solutions of the closed shell Kohn–Sham equation. J. Chem. Phys. 1996, 104, 9047–9052. (116) Sun, J.; Ruzsinszky, A.; Perdew, J. Strongly constrained and appropriately normed semilocal density functional. Phys. Rev. Letters 2015, 115, 036402. (117) Brandenburg, J. G.; Bates, J. E.; Sun, J.; Perdew, J. P. Benchmark tests of a strongly constrained semilocal functional with a long-range dispersion correction. Phys. Rev. B 2016, 94, 115144. ¨ (118) Eichkorn, K.; Treutler, O.; Ohm, H.; Häser, M.; Ahlrichs, R. Auxiliary basis sets to approximate Coulomb potentials. Chem. Phys. Lett. 1995, 242, 652–660. (119) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials. Theor. Chem. Acc. 1997, 97, 119–124. (120) Izsák, R.; Neese, F. An overlap fitted chain of spheres exchange method. J. Chem. Phys. 2011, 135, 144105.



(121) Schwabe, T. Systematic study of the basis set superposition error in core-electron correlation effects. J. Phys. Chem. A 2013, 117, 2879–2883. (122) Zhao, Y.; Truhlar, D. G. Design of density functionals that are broadly accurate for thermochemistry, thermochemical kinetics, and nonbonded interactions. J. Phys. Chem. A 2005, 109, 5656–5667. (123) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (124) Sure, R.; Grimme, S. Corrected small basis set Hartree-Fock method for large systems. J. Comput. Chem. 2013, 34, 1672–1685. (125) Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A. Consistent structures and interactions by density functional theory with small atomic orbital basis sets. J. Chem. Phys. 2015, 143, 054107. (126) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (127) Brandenburg, J. G.; Bannwarth, C.; Hansen, A.; Grimme, S. B97-3c: A revised lowcost variant of the B97-D density functional method. J. Chem. Phys. 2018, 148, 064104. (128) Grimme, S. Semiempirical GGA-type density functional constructed with a longrange dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. (129) Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. Accurate prediction of noncovalent interaction energies with the effective fragment


Page 54 of 62

Page 55 of 62 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


potential method: Comparison of energy components to symmetry-adapted perturbation theory for the S22 test set. J. Chem. Theory Comput. 2012, 8, 2835–2843. (130) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479–483. (131) Klopper, W.; Manby, F. R.; Ten-No, S.; Valeev, E. F. R12 methods in explicitly correlated molecular electronic structure theory. Int. Rev. Phys. Chem. 2006, 25, 427–468. (132) Ma, Q.; Werner, H.-J. Explicitly correlated local coupled-cluster methods using pair natural orbitals. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2018, 8, e1371. (133) Karton, A. A computational chemists guide to accurate thermochemistry for organic molecules. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2016, 6, 292–310. (134) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. J. Chem. Phys. 2013, 139, 134101. (135) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. (136) Karton, A.; Martin, J. M. L. Comment on: “Estimating the Hartree–Fock limit from finite basis set calculations” [Jensen F (2005) Theor Chem Acc 113:267]. Theor. Chem. Acc. 2005, 115, 330–333.



(137) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Basis-set convergence in correlated calculations on Ne, N2, and H2O. Chem. Phys. Lett. 1998, 286, 243–252. (138) Neese, F.; Valeev, E. F. Revisiting the atomic natural orbital approach for basis sets: Robust systematic basis sets for explicitly correlated and conventional correlated ab initio methods? J. Chem. Theory Comput. 2011, 7, 33–43. (139) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (140) Zheng, J.; Xu, X.; Truhlar, D. G. Minimally augmented Karlsruhe basis sets. Theor. Chem. Acc. 2010, 128, 295–305. (141) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (142) Liakos, D. G.; Sparta, M.; Kesharwani, M. K.; Martin, J. M. L.; Neese, F. Exploring the accuracy limits of local pair natural orbital coupled-cluster theory. J. Chem. Theory Comput. 2015, 11, 1525–1539. (143) Martin, J. M. L.; Oliveira, G. D. Towards standard methods for benchmark quality ab initio thermochemistry — W1 and W2 theory. J. Chem. Phys. 1999, 111, 1843– 1856. (144) Karton, A.; Martin, J. M. L. Explicitly correlated Wn theory: W1-F12 and W2-F12. J. Chem. Phys. 2012, 136, 124114. (145) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. The magnitude


Page 56 of 62

Page 57 of 62 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


of the CH/π interaction between benzene and some model hydrocarbons. J. Am. Chem. Soc. 2000, 122, 3746–3753. (146) Jurecka, P.; Hobza, P. On the convergence of the (∆E CCSD(T ) -∆E M P 2 ) term for complexes with multiple H-bonds. Chem. Phys. Lett. 2002, 365, 89–94. (147) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 618. (148) Grimme, S. Improved second-order Møller–Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies. J. Chem. Phys. 2003, 118, 9095–9102. (149) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. J. Chem. Phys. 2015, 143, 034108. (150) Rezac, J.; Bim, D.; Gutten, O.; Rulisek, L. Toward accurate conformational energies of smaller peptides and medium-sized macrocycles: MPCONF196 benchmark energy data set. J. Chem. Theory Comput. 2018, 14, 1254–1266. (151) Mardirossian, N.; Head-Gordon, M. ωB97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-the-fittest strategy. Phys. Chem. Chem. Phys. 2014, 16, 9904–9924. (152) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098.



(153) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785. (154) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Results obtained with the correlation energy density functionals of Becke and Lee, Yang and Parr. Chem. Phys. Lett. 1989, 157, 200–206. (155) Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372–1377. (156) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 33, 8822. (157) Perdew, J. P. Erratum: Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 34, 7406. (158) Handy, N. C.; Cohen, A. J. Left-right correlation energy. Mol. Phys. 2001, 99, 403– 412. (159) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. (160) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew–Burke–Ernzerhof exchangecorrelation functional. J. Chem. Phys. 1999, 110, 5029–5036. (161) Zhang, Y.; Yang, W. Comment on “Generalized gradient approximation made simple”. Phys. Rev. Letters 1998, 80, 890. (162) Hujo, W.; Grimme, S. Performance of the van der Waals Density Functional VV10 and (hybrid)GGA Variants for Thermochemistry and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 3866–3871.


Page 58 of 62

Page 59 of 62 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


(163) Mardirossian, N.; Head-Gordon, M. Mapping the genome of meta-generalized gradient approximation density functionals: The search for B97M-V. J. Chem. Phys. 2015, 142, 074111. (164) Chai, J.-D.; Head-Gordon, M. Long-range corrected double-hybrid density functionals. J. Chem. Phys. 2009, 131, 174105. (165) Zhao, Y.; Truhlar, D. G. A new local density functional for main-group thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions. J. Chem. Phys. 2006, 125, 194101. (166) Yu, F. Double-hybrid density functionals free of dispersion and counterpoise corrections for non-covalent interactions. J. Phys. Chem. A 2014, 118, 3175–3182. (167) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Workhorse semilocal density functional for condensed matter physics and quantum chemistry. Phys. Rev. Lett. 2009, 103, 026403. (168) Kozuch, S.; Gruzman, D.; Martin, J. M. L. DSD-BLYP: A general purpose double hybrid density functional including spin component scaling and dispersion correction. J. Phys. Chem. C 2010, 114, 20801–20808. (169) Kozuch, S.; Martin, J. M. L. DSD-PBEP86: In search of the best double-hybrid DFT with spin-component scaled MP2 and dispersion corrections. Phys. Chem. Chem. Phys. 2011, 13, 20104–20107. (170) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta–generalized gradient approximation designed for molecules and solids. Phys. Rev. Letters 2003, 91, 146401.



(171) Alipour, M. Seeking for spin-opposite-scaled double-hybrid models free of fitted parameters. J. Phys. Chem. A 2016, 120, 3726–3730. (172) Tarnopolsky, A.; Karton, A.; Sertchook, R.; Vuzman, D.; Martin, J. M. L. Doublehybrid functionals for thermochemical kinetics. J. Phys. Chem. A 2008, 112, 3–8. (173) Mardirossian, N.; Head-Gordon, M. ωB97M-V: A combinatorially optimized, rangeseparated hybrid, meta-GGA density functional with VV10 nonlocal correlation. J. Chem. Phys. 2016, 144, 214110. (174) Santra, G.; Sylvetsky, N.; Martin, J. M. L. Minimally empirical double-hybrid functionals trained against the GMTKN55 database: revDSD-PBEP86-D4, revDODPBE-D4, and DOD-SCAN-D4. J. Phys. Chem. A 2019, 123, 5129–5143. (175) Lin, Y.-S.; Li, G.-D.; Mao, S.-P.; Chai, J.-D. Long-range corrected hybrid density functionals with improved dispersion corrections. J. Chem. Theory Comput. 2012, 9, 263–272. (176) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals density functional: The simpler the better. J. Chem. Phys. 2010, 133, 244103. (177) Perdew, J. P.; Schmidt, K. Jacob’s Ladder of Density Functional Approximations for the Exchange-Correlation Energy. AIP Conf. Proc. 2001, 577, 1–20. (178) Caldeweyher, E.; Bannwarth, C.; Grimme, S. Extension of the D3 dispersion coefficient model. J. Chem. Phys. 2017, 147, 034112. (179) Caldeweyher, E.; Ehlert, S.; Hansen, A.; Neugebauer, H.; Spicher, S.; Bannwarth, C.; Grimme, S. A generally applicable atomic-charge dependent London dispersion correction. J. Chem. Phys. 2019, 150, 154122.


Page 60 of 62

Page 61 of 62 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


(180) Ghasemi, S. A.; Hofstetter, A.; Saha, S.; Goedecker, S. Interatomic potentials for ionic systems with density functional accuracy based on charge densities obtained by a neural network. Phys. Rev. B 2015, 92, 045131.




Page 62 of 62

Toward a Quantum-Chemical Benchmark Set for ... - ACS Publications

Recommend Documents