Benchmark Calculations of Interaction Energies in Noncovalent

Review pubs.acs.org/CR

Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications Jan Ř ezác*̌ ,† and Pavel Hobza†,‡ †

Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, 166 10 Prague, Czech Republic Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Palacký University, 771 46 Olomouc, Czech Republic

‡

ABSTRACT: Data sets of benchmark interaction energies in noncovalent complexes are an important tool for quantifying the accuracy of computational methods used in this field, as well as for the development of new computational approaches. This review is intended as a guide to conscious use of these data sets. We discuss their construction and accuracy, list the data sets available in the literature, and demonstrate their application to validation and parametrization of quantum-mechanical computational methods. In practical model systems, the benchmark interaction energies are usually obtained using composite CCSD(T)/CBS schemes. To use these results as a benchmark, their accuracy should be estimated first. We analyze the errors of this methodology with respect to both the approximations involved and the basis set size. We list the most prominent data sets covering various aspects of the field, from general ones to sets focusing on specific types of interactions or systems. The benchmark data are then used to validate more efficient computational approaches, including those based on explicitly correlated methods. Special attention is paid to the transition to large systems, where accurate benchmarking is difficult or impossible, and to the importance of nonequilibrium geometries in parametrization of more approximate methods.

CONTENTS 1. Introduction 1.1. Noncovalent Interactions 1.2. The Scope of This Review 1.3. Benchmark Calculations of Noncovalent Interactions 1.4. Nature of Noncovalent Interactions 2. Methodology 2.1. Interaction Energy 2.2. Counterpoise Correction of the Basis Set Superposition Error 2.3. Basis Sets 2.4. Extrapolation to the Complete Basis Set Limit 2.5. Explicitly Correlated Methods 2.6. Composite CCSD(T)/CBS Schemes 2.7. Post-CCSD(T) Calculations 2.8. Approximate CCSD(T) Approaches 2.9. Quantum Monte Carlo 2.10. Wave Function Methods 2.11. DFT and Semiempirical Methods 2.12. Intermolecular Perturbation Theory and Interaction Energy Decomposition 3. Data Set Construction and Processing 3.1. Data Set Composition 3.2. Geometries 3.3. Nonequilibrium Geometries 3.4. Categorization of Interaction Types 3.5. Statistical Processing of the Results © 2016 American Chemical Society

3.6. Automated Data Set Processing 3.7. Versioning of Data Sets 4. Available Data Sets 4.1. Accurate Data Sets of the Smallest Model Systems 4.2. General Data Sets Covering Interactions of Organic Molecules 4.3. S66-Compatible Data Sets 4.4. Data Sets of Specific Interaction Types 4.5. Data Sets Targeting Specific Systems 4.6. Data Sets of Related Properties 4.7. Resources and Aggregated Data Sets 5. Important Results 5.1. The Accuracy of CCSD(T)/CBS 5.2. The Validation and Parameterization of Wave Function Methods 5.3. Transition to Large Systems 5.4. The Parametrization of Semiempirical QM Methods 6. A Summary 6.1. A State-of-the-Art Summary 6.2. Open Questions Author Information Corresponding Author Notes

5039 5039 5039 5039 5040 5040 5040 5041 5042 5042 5043 5044 5046 5047 5047 5047 5048 5049 5049 5049 5050 5050 5051 5051

5052 5052 5052 5052 5054 5057 5057 5057 5058 5059 5060 5060 5060 5062 5063 5064 5064 5064 5065 5065 5065

Special Issue: Noncovalent Interactions Received: September 8, 2015 Published: March 4, 2016 5038

DOI: 10.1021/acs.chemrev.5b00526 Chem. Rev. 2016, 116, 5038−5071

Chemical Reviews

Review

1. INTRODUCTION

by our group1−6 and by others.7,8 We will briefly mention related topics, such as more approximate calculations at semiempirical and density functional theory (DFT) levels, which are covered in depth by separate reviews in this issue. The methodology used in the paper is introduced starting from the basics to make it accessible to a broader audience.

1.1. Noncovalent Interactions

1.3. Benchmark Calculations of Noncovalent Interactions

While the structure and properties of small molecules are determined by covalent bonding, in any larger system, noncovalent interactions become equally important. They determine not only the physicochemical properties of noncovalent complexes, clusters, and condensed molecular matter but also the structure and properties of larger molecules, where noncovalent interactions occur between the building blocks constituting the molecule. Understanding noncovalent interactions is thus essential for many fields of chemistry. Noncovalent interactions are of great importance in biochemical processes, where they determine the structure, dynamics, and function of biomolecules. Here, life takes advantage of their reversibility and favorable energetics at ambient temperature; noncovalent interactions are strong enough to bind the molecules together for a sufficiently long time while being weak enough to be assembled and disassembled without too much energy expenditure. The same principles may be exploited in the design of new synthetic materials with properties and functions engineered for a specific purpose. In contrast to covalent bonding, noncovalent interactions are much subtler and affect the electronic structure of the interacting species only slightly. This must be reflected in the methods used in both experiment and theory. Isolating the effects of noncovalent interactions experimentally requires advanced, highly sensitive methods. In theoretical approaches, the principles of covalent bonding are already described in very approximate quantum mechanical (QM) methods, whereas noncovalent interactions can be qualitatively described only after long-range electron correlation is taken into account in post-Hartree−Fock methods. Quantitative accuracy can be reached by ab initio (nonempirical) means only when advanced correlated methods are used and the calculations are carried out in a large basis set. With the recent advances in both methodology and computer hardware, computational chemistry has become an important tool for rationalizing, modeling, and predicting chemical processes in systems of increasing size and complexity. In addition to reproducing the experiment, it provides insights at the atomic or electronic structure levels that would otherwise be difficult or impossible to obtain. Passing to larger systems, the ability of the computational methods to describe noncovalent interactions becomes more and more important. This is attested by the increasing attention to noncovalent interactions in the development of computational methods.

The ultimate test of a theoretical model is, of course, its comparison with experiment. Energetics of noncovalent interactions in smaller complexes can be measured directly in the gas phase using spectroscopic methods (there are multiple reviews covering the experiments in this issue) with accuracy that matches or exceeds the accuracy of the best computational methods. There are, however, two important limitations. First, the amount of available experimental data is rather small, and some important classes of noncovalent interactions have not been measured yet. Second, the measured quantity, interaction enthalpy at zero temperature ΔH00 (often called dissociation energy, D0), cannot be calculated directly, but it has to be composed from multiple components, some of which are difficult to calculate with the desired accuracy. The limiting step is usually the calculation of the change of zero-point vibrational energy, ΔZPVE, which requires the calculation of not only electronic but also vibrational properties of the system. To avoid these limitations, a common practice is to test the accuracy of approximate computational methods by comparing them with more accurate calculations. This can be done at the level of interaction energies, which are easy to calculate (see section 2.1). The conclusions drawn on the basis of such benchmarking studies are, however, not limited just to reproducing the dissociation energies. Theoretical prediction of more complex, experimentally observable quantities, such as rotational and vibrational spectra of the molecular complexes, relies on accurate description of the potential energy surface, which could be only as good as the interaction energy at a specific point. The recent trend of extending the benchmark interaction energy calculations also to nonequilibrium geometries provides further information in this direction. The question then is what computational methods to use as the benchmark and what their accuracy and computational demands are. Such a method should be accurate, robust, and reproducible. Among the methods that can be used, the coupled-cluster (CC) approach9 offers the possibility of improving the description of electron correlation systematically by adding excitation operators of increasing order (for more information, see, for example, the recent review by Bartlett and Musiał10). It has been found that for an accurate description of noncovalent interactions, triple excitations (triples) have to be used. These are fully covered by the coupled-cluster single, double, and triple excitation method (CCSDT),11 which is, however, very expensive, as it is iterative and scales with N8 (where N is an approximate measure of the size of the system combining the number of electrons and basis set size; for details, see section 2.7). This scaling is reduced to N7 when the triples are added to CCSD (iterative, scales with N6) by means of perturbation theory (up to fourth order, including also one term from the fifth order) in the CCSD(T) methodology.12 Here, the most demanding step, the calculation of the triples, is not iterative. At this level, the most important source of error is the incompleteness of the basis set. While the scaling limits the use of large basis sets, the energy at the complete basis set limit (CBS) can be estimated by extrapolation. The resulting

Biographies Acknowledgments References

5065 5065 5065

1.2. The Scope of This Review

The field of noncovalent interactions has a long history, which has been reviewed periodically, including thematic issues of Chemical Reviews in 1988, 1994, and 2000, edited by one of us (P.H.). The present issue is a continuation of this series, with the authors of this paper being also guest editors. In this review, we will summarize the recent advances in accurate calculations of noncovalent interactions by means of wave function theory (WFT) QM methods and their application in benchmarking more approximate QM approaches. It follows previous reviews 5039


Chemical Reviews

Review

elctrostatic and dispersion terms are of comparable size. Among the dispersion-dominated interactions, the π−π interactions (interactions involving contact between π orbitals, also called π−π stacking) are often set apart because of their specific properties. In these systems, the interaction is rather strong because the molecules are in very close contact via multiple heavier (non-hydrogen) atoms that have higher polarizability. Although π−π interactions are common in complexes of aromatic molecules, the aromaticity does not make the interaction stronger.20 Dispersion-dominated interactions represent an important test case for computational methods, as their accurate description requires either advanced correlated ab initio methods or empirical treatment. Among the interactions of mostly electrostatic nature, multiple classes had been defined on the basis of specific structural motifs. The most prominent category here is the hydrogen bond.21 The H-bond is a noncovalent interaction between an electron-deficient hydrogen covalently bound to an electronegative atom and a region of high electron density, most often a lone pair on another electronegative atom (for more precise definition, see ref 22). In addition to electrostatics, hydrogen bonding is associated with charge transfer accompanied by changes of electronic structure that can be described as partial covalent bonding.23 A special case of hydrogen bonding is the dihydogen bond,24 a H-bond where the acceptor atom is also hydrogen, in this case covalently bound to an electropositive atom such as alkali metal. Dihydrogen bonds are also obsreved in boron compounds.25 Another broader class of interactions is the σ-hole bonding,26 an interaction of a positively charged region (so-called σhole27,28) on an otherwise electronegative atom with another negatively charged site. The σ-hole lies in the extension of a covalent bond and it is caused by reduced occupation of the porbital involved in the bond compared to those perpendicular to it. This effect is most common in halogen compounds, where it is called the halogen bond.29−31 In most halogen bonds, the electrostatic component (which primarily defines the geometry of the complex) is complemented by equally important dispersion interaction.32 Due to the anisotropy of the charge density, also the repulsion and dispersion contributions to halogen bonding are orientationally dependent.33,34 The σ-hole is not limited only to halogens; the hole could also be found in other covalently bound atoms of groups IV−VI26 and, surprisingly, also VIII. In groups IV, V, VI, and VIII, this interaction is called tetrel, pnictogen,35 chalcogen,36 and aerogen bond.37 In this special issue, σ-hole bonding is addressed with a dedicated review by Kolár ̌ and Hobza.38

combination, CCSD(T)/CBS, yields an accurate and reliable description of noncovalent interactions, yet it is only applicable to molecular systems with no more than several tens of atoms. Because of this favorable accuracy-to-cost ratio, it is known as the “gold standard” of computational chemistry. It should be noted that there are also other approximate versions of CCSDT, both perturbative (CCSD[T]13) and iterative (CCSDT-1 to CCSDT-413−15). Despite being more demanding, the iterative approaches do not perform better than the perturbative ones.16 CCSD[T] has been shown to perform slightly better than CCSD(T) in small basis sets,16,17 but the difference does not justify its use instead of the widespread and reportedly more robust18 CCSD(T). In this review, we will focus closely on two aspects of the CCSD(T)/CBS calculations of noncovalent interactions. First, although CCSD(T)/CBS calculations are accurate, they are not exact. To handle the benchmark CCSD(T) calculations properly, it is necessary to know their accuracy first. Here, we will summarize the attempts in quantifying the accuracy of CCSD(T)/CBS interaction energies, mainly in comparison to higher-level CC calculations and with respect to more accurate estimates of the complete basis set limit. The effect of other approximations involved will be discussed as well. Second, CCSD(T)/CBS calculations are used not only to solve specific chemical problems but also for the benchmarking and parametrization of more approximate methods. Here, it is necessary to work with statistically relevant amounts of data, and it would be impractical to repeat the CCSD(T) calculations in each study. It is thus common to reuse previously published results; there exist benchmark data sets designed specifically for this purpose. The most important advantage of benchmarking computational methods against widely accepted data sets is that the results can be then compared across different studies. In this review, we will list and describe the commonly used data sets and some of their applications. 1.4. Nature of Noncovalent Interactions

Noncovalent interactions are often classified by their nature. Although all the interactions are governed by only few fundamental physical principles, various combinations of the individual contributions give rise to classes of noncovalent interactions of different character. This analysis can be done quantitatively by the means of computational methods that provide decomposition of the interaction energy to the individual terms. More detailed classification is then based on the presence of specific structural motifs or elements involved in the interactions. This classification is less rigorous, but there exist generally accepted definitions of the individual classes. While this classification itself is not the subject of this rewiew, it is important for the development of benchmark data sets because these sets are often designed to cover some specific classes of interactions. Here, we will briefly overview the usual categorization of noncovalent interactions; the practical use of such categorization in data set design is discussed later in section 3.4. To define the physical nature of a noncovalent interaction, London dispersion19 is usually put into contrast with the rest of electrostatic interactions (Coulombic interaction and induction). The interaction of nonpolar molecules is dominated by the dispersion; on the other end of the spectrum, the interaction of charged or highly polar molecules is mostly of electrostatic nature. Often, a third category of ”mixed character” interactions is inserted in between, covering systems where the

2. METHODOLOGY 2.1. Interaction Energy

The quantity central to this review is interaction energy, ΔEint, defined as the difference between the energy of the noncovalent complex and the sum of energies of the monomers constituting it. The energy of the monomers is evaluated at their geometry taken from the complex; i.e., deformation energy, ΔEdef, is not included in the ΔEint. In a dimer AB, consisting of molecules A and B, the interaction energy becomes ΔEint(AB) = E(AB) − E(A) − E(B)

(1)

In the so-called supramolecular approach, all the three energies on the right-hand side of eq 1 are calculated separately. This is the only possible option for the vast majority of the methods 5040


Chemical Reviews

Review

discussed in this work. Alternatively, in some methods, the interaction energy is obtained without the explicit evaluation of the energies of the individual subsystems, e.g., by means of perturbation theory. It must be noted that the ΔEint is not an observable but an artificial construct. In experiment, the only quantity accessible is the interaction enthalpy ΔH00, often referred to as dissociation energy, D0. It consists of three components: the interaction energy, ΔEint; the deformation energy, ΔEdef, describing the energy penalty associated with the change of the geometry of the monomers upon the formation of the complex; and the difference (calculated analogously to eq 1) of zero-point vibration energies, ΔZPVE. The reason why the interaction energy is the preferred quantity for benchmarking the computational method is the fact that it can be calculated from single-point energy calculations on a fixed geometry. Not only is it the most economical solution, but it is also a well-defined and easily reproducible procedure. For these reasons, it has become the approach used by default in the field. While the interaction energy is only one piece of the dissociation energy, understanding the accuracy of the calculations of ΔEint is the first step to accurate evaluation of the remaining components. If the deformation energy were added, costly geometry optimizations would be needed and some uncertainty would be introduced if a different setup of the optimization algorithm were used. On the other hand, there are no principal obstacles to evaluating the deformation energy with the same accuracy as the interaction energy. Finally, the calculation of ΔZPVE requires even more demanding calculations. To match the accuracy of the interaction energy, advanced anharmonic vibrational calculations are needed. Generally applicable, nonempirical approaches are based on complete characterization of the multidimensional potential energy surface (PES) that requires an enormous computational effort that grows exponentially with system size. In practice, it is not feasible to cover the whole PES at the same level as the calculation of ΔEint for all but simplest systems. For larger molecular complexes, such as those discussed in this review, further approximations and empirical assumptions have to be made. While it is possible to achieve accuracy approaching the experiment in an individual case, a general method applicable to routine processing of large sets of data is not available. The ΔZPVE, or the vibrational frequencies themselves, can be calulated at high level under the harmonic approximation and used as a benchmark, but the tight relation to experimentally observable D0 is lost. If ΔEdef and ΔZPVE are neglected in the benchmarking of computational methods, some information is evidently lost. This information characterizes the potential energy surface around the equilibrium geometry, and it can be recovered more economically by extending the ΔEint calculations to nonequilibrium geometries as described later. In our work, we use the sign convention implied by eq 1, where attractive interactions have negative values of ΔEint.

basis set functions are available. According to the variation principle, this leads to an artificial lowering of the energy of the complex (as compared to the dissociated state) and thus also of the interaction energy (in other words, the stability of the complex is overestimated). This issue can be mitigated by using the counterpoise (CP) correction introduced by Boys and Bernardi.39 In this approach, the monomers are calculated in the same basis as the complex, including the empty basis set functions located at the positions of the atoms of the other molecule(s). In a dimer, the interaction energy is calculated as follows

2.2. Counterpoise Correction of the Basis Set Superposition Error

The application of this scheme to the gradient in geometry optimizations is necessary for obtaining accurate geometries of noncovalent complexes when small- to medium-sized basis sets (up to ca. triple-ζ) are used. Despite making the calculations more complex, this approach is more economical than using a basis set large enough to make the effects of the BBSE negligible. In this review, all interaction energies listed are counterpoise-corrected unless explicitly noted otherwise.

int ΔECP (AB) = E(AB) − E(AAB) − E(BAB)

(2)

with the superscript AB denoting a calculation in a basis set of the whole dimer. This approach has been proven as very effective and it has become used widely. While some authors argue that it leads to overcorrecting the BSSE and may deteriorate the results, it happens only in specific cases, such as in atomic dimers where the basis set is optimized for a free atom calculation.40 The formal analysis of possible approaches to correcting the BSSE41 as well as numerical studies, including the recent ones based on very accurate reference data,42 prove that counterpoise correction is safe for general use. More detailed argumentation in favor of CP correction can be found in a review by van Duijneveldt.43 It must be admitted, however, that a CP-corrected result is not necessarily better (in the sense of being closer to the complete basis set limit) than the uncorrected one, because in any finite basis, the BSSE mixes with, and may accidentally be compensated by, other errors caused by the incompleteness of the basis set or the errors of the method itself.44 Some authors use a partial correction, most often taking the average of the corrected and uncorrected interaction energies (note that this requires two extra calculations). Although this empirical approach improves the results obtained in small basis sets, it is not effective in larger ones.42 Since the BSSE is not the only error caused by an insufficient basis set size, the combination of a sufficiently large basis set and counterpoise correction seems to be an effective and safe choice. As the CP correction allows for better isolation of the intrinsic error of the basis set, it makes the dependence of the results on basis set size smoother, which improves the accuracy of extrapolations to the complete basis set limit.45 The use of counterpoise correction can also be recommended for consistency with already published benchmark data, the majority of which are CP-corrected. The CP correction can also be applied to the calculation of the total energy and its derivatives.46,47 In contrast to the interaction energy calculation (eq 2), five calculations are needed to isolate and remove the BSSE: ECP(AB) = E(AB) − E(AA ) − E(BB) + E(AAB) + E(BAB) (3)

The use of finite basis sets in the supramolecular calculations of ΔEint produces the basis set superposition error (BSSE). It is caused by the fact that in the calculation of the complex, each molecule can use not only the basis set functions located on its own atoms but also the ones located on the other molecules, whereas in the calculation of the isolated monomers, no such 5041


Chemical Reviews

Review

2.3. Basis Sets

be mentioned specifically in the text. The details can be found in the original papers. To save space and increase the readability of the text, the following short notation will be used for the correlation-consistent basis sets: cc-pVXZ is abbreviated as XZ, aug-cc-pVXZ as aXZ, and heavy-aug-ccpVXZ as haXZ. The Dunning basis sets provide a reliable and systematic framework for the calculations of noncovalent interactions, but even the smallest basis sets from this series might be prohibitively large for application to very large systems or for methods scaling very steeply with the basis set size. The standard basis sets at this level (stepping down from double-ζ to split-valence basis sets) do not perform well for noncovalent interactions, but this can be improved by modifying them specifically for this purpose. The most successful approach is to use the 6-31G* or 6-31G** basis sets of Pople60 (the 6-31G basis with polarization functions on non-hydrogen and all atoms, respectively), where the polarization functions are made more diffuse.61 In particular, the exponents of the polarization functions of C, N, and O are reduced from 0.8 to 0.25; for hydrogen, the exponent is changed from 1.1 to 0.15. These modified basis sets are usually denoted as 6-31G*(0.25) and 631G**(0.25,0.15). These basis sets have been found to perform very well for noncovalent interactions of various nature.2,62 These basis sets are very useful for the calculation of the higherorder correction in composite schemes (see section 2.6).63,64

The importance of electron correlation in noncovalent interactions leads to specific requirements on the basis sets used in the calculations. Basis set functions with a high angular momentum are needed to model the effect of the correlation cusp on the wave function in correlated WFT methods. The most frequently used series of basis sets developed for this purpose are the correlation-consistent basis sets of Dunning,48 labeled cc-pVXZ, where X = D, T, Q, 5, and 6 for basis sets ranging from double- to sextuple-ζ. For an accurate description of the molecular properties that affect their noncovalent interactions, the basis sets should be augmented with diffuse functions, which leads to the aug-cc-pVXZ series of basis sets.49 However, the diffuse functions on hydrogen atoms can be neglected with only a small impact on the results. Therefore, the so-called heavy-augmented basis sets are often used. We denote these basis sets as heavy-aug-cc-pVXZ; recently, they have also been labeled as jul-cc-pVXZ.50 In these basis sets, the polarization functions are added only to the valence shell, as the majority of the correlated calculations are carried out within the frozen-core approximation (i.e., only the valence electrons are considered). The error introduced by this approximation can be neglected in most applications. If more accurate results are sought, calculations correlating all electrons can be carried out; in such a case, it is appropriate to use the polarized-core51 (ccpCVXZ) or weighted polarized-core (cc-pwCVXZ)52 versions of the basis sets. In the calculations of heavier elements, these basis sets can be combined with pseudopotentials modeling the relativistic effects. The appropriate combinations of correlation-consistent basis sets and pseudopotentials adjusted for use with these basis sets53 are labeled cc-pVXZ-PP. Many of the computational methods used in this work can be accelerated by evaluating the multicenter integrals using the resolution of identity (RI) or density fitting (DF) approximations. They can be utilized at both the self-consistent field (SCF) level in Hartree−Fock54 (HF) or DFT calculations and in the subsequent calculation of the correlation energy.55 Appropriate auxiliary basis sets are needed for the fitting; for the correlation-consistent basis sets, the corresponding auxiliary sets54,56 are used. In practical calculations, the error introduced by these approximations is negligible. The auxiliary basis sets are designed so that the error in the electronic energy due to the fitting is at least 1 order of magnitude smaller than the error caused by the incompleteness of the AO basis. In the interaction energies, part of the error cancels and the accuracy is even better. The default auxiliary basis sets were recommended even for very accurate calculations as long as diffuse functions are used.57 The error introduced by the RI approximation in calculations of interaction energies in nucleic acid base pairs in rather small basis sets was reported to be about 0.03 kcal/mol.58 In our calculations of interaction energies in the A24 data set,59 we found that the RI approximation in MP2/aug-cc-pVTZ calculations leads to an average error of 0.08% when the approximation is applied only in the calculation of the correlation energy and 0.10% when applied also at the HF level. The majority of the results presented in this review have been calculated using correlation-consistent basis sets within the frozen-core approximation. The exceptions will be listed individually in the discussion. A large part of the calculations was performed using RI or DF approximation, which will not

2.4. Extrapolation to the Complete Basis Set Limit

The steep scaling of the computational complexity of correlated quantum-mechanical methods imposes severe restrictions on the basis set size. In practical calculations with benchmarkquality methods such as CCSD(T), the incompleteness of the basis set becomes the main source of the error with respect not only to the complete basis set (CBS) limit but also to the exact solution.64 The same applies to higher-order post-CCSD(T) methods in the smallest model complexes. The convergence of the result with the basis set size can be divided into two contributions. Within the one-electron approximation at the SCF level, the basis set has to be large enough to provide a good description of the molecular orbitals and their changes caused by noncovalent interactions. This contribution, namely HF or DFT energy, converges relatively fast with the basis set size; a triple-ζ basis set is often sufficient and a quadruple-ζ basis set yields results close to the CBS limit. On the other hand, the convergence of correlation energy is much slower and basis sets larger than quadruple-ζ are needed. Such calculations can be extremely demanding or even intractable. There are two possible approaches to alleviating this problem. A series of basis sets of increasing size can be used to extrapolate the result to the CBS limit. Alternatively, methods covering the electron correlation explicitly can be used to accelerate the convergence, as discussed in the following section. The description of electron correlation depends on basis set functions with a higher angular momentum. In a series of otherwise saturated basis sets with the increasing highest angular momentum quantum number l, the correlation energy should improve systematically, making it possible to extrapolate the correlation energy. Basis sets specifically designed to meet these requirements should be used, such as the correlationconsistent basis sets described above. It can be derived that the correlation energy convergence follows approximately the inverse of the third power of the basis set size expressed by 5042


Chemical Reviews

Review

Density fitting is often employed in solving the additional integrals. The convergence of the correlation energy with the maximum angular number of the basis set l is improved from (l + 1)−3 to (l + 1)−7, and explicitly correlated calculations in a triple-ζ basis set yield results comparable to the ones in a quintuple-ζ basis without explicit correlation. For more technical details, we refer the reader to reviews on this topic.75,76 The explicit correlation is usually applied to second-order Møller−Plesset perturbation theory (MP2) and coupled-cluster calculations. In the coupled-cluster method, practical implementations up to CCSD-F12 are now available.77,78 These can serve as a basis of CCSD(T) calculations where the triples are calculated from the CCSD-F12 amplitudes but without the explicit correlation. These methods are generally referred to as CCSD(T)-F12, although the explicit correlation is not complete. The error due to the (T) term may limit the overall accuracy of the calculation.79 Full explicitly correlated CCSD(T) calculations have been also reported80 but are not generally used yet due to very high computational demands. Omitting the explicit correlation in the (T) term leads to analogous results as if it was calculated in smaller basis than the CCSD contribution. This imbalance introduces some error into the resulting interaction energies (this is demonstrated in section 2.6). An empirical correction has been devised, where the contribution of triples is scaled by a factor calculated as the ratio of MP2 to MP2-F12 correlation energies [denoted as CCSD(T*)-F12].79 Multiple versions of CCSD(T)-F12 involving different approximations are available,81−83 which, in combination with the optional scaling of the triples, leads to many variants of the method. These were systematically studied in the context of noncovalent interactions by Patkowski,57 who also reviewed previous applications of explicitly correlated methods to noncovalent interactions. He concluded that in large basis sets, the approximations involved in the CCSD-F12 and the inexact treatment of the triples limit the accuracy and that the conventional CCSD(T) yields similar or better results. In smaller basis sets, where the basis set size is by far the most limiting factor, the F12 approach is advantageous, but the same issue can be addressed better by basis set extrapolation (see section 2.4) and composite schemes (see section2.6). On the other hand, MP2-F12 calculations74 do not suffer from these problems, and the use of explicit correlation is a viable tool for approaching the MP2/CBS limit efficiently. While the use of explicit correlation makes the calculations more demanding, the scaling is not affected, and this overhead can be compensated for by using a smaller basis set. While the explicitly correlated calculations can theoretically be used with any basis set, the series of correlation-consistent basis sets was reoptimized specifically for this purpose, and the resulting cc-pVXZ-F12 basis sets84 are an obvious choice for these calculations. Although the name does not suggest that, these basis sets do contain diffuse functions and thus correspond to the aug-cc-pVXZ series. As all the practical implementations of explicitly correlated methods rely on density fitting, specialized auxiliary basis sets85 are also needed. It is further possible to extrapolate the explicitly correlated results to CBS. Since the convergence properties are different, a dedicated parametrization of the extrapolation formulas is necessary.86 However, the extrapolation does not have much practical value as in small (double-ζ) basis sets, the error is governed by other limitations of the basis rather than its ability

its cardinal number X (X = l + 1; X = 2 for double-ζ, 3 for triple-ζ, etc.). In this approach introduced by Helgaker et al.,65,66 the basis set dependence of the correlation energy on X can be expressed as corr ΔEXcorr = ECBS + aX −3

(4)

with a being a system-dependent parameter. To eliminate this parameter, the extrapolation has to be made from at least two points, usually from basis sets of sizes X and X + 1. The correlation energy at the CBS limit can then be expressed as corr ΔECBS =

3 corr (X + 1)3 EXcorr + 1 − X EX

(X + 1)3 − X3

(5)

This is the most widely used approach to the extrapolation to the CBS limit. It is used in a large part of the calculations discussed in this review. This extrapolation assumes that the limited description of electron correlation is the only error of the basis sets. In a small basis set, however, the other sources of error are equally important. This makes the extrapolation inaccurate when small basis sets are used. In general, extrapolation from double- to triple-ζ basis sets is not reliable and the result may be worse than the one from a triple-ζ basis set. There have been many attempts to modify the extrapolation procedure to account for the deficiency of double-ζ basis sets empirically,67−69 but for noncovalent interactions, none of them work better than the original Helgaker formula.64 It was shown earlier that in these small basis sets, the results of various possible extrapolation approaches depend strongly on the nature of the interaction.70 On the other end of the accuracy spectrum in smaller systems, the whole series of calculations in basis sets from triple-ζ to quintuple-ζ or even sextuple-ζ can be used in the extrapolation. One possibility is to fit the series with eq 4 to make the extrapolation more robust. However, we prefer the more flexible approach of making the exponent in the power law formula variable and then fitting it in a series of the three largest basis sets.59 This approach is better justified, as the theoretically derived scaling is close to eq 4 but does not follow it exactly. 2.5. Explicitly Correlated Methods

In standard correlated QM methods, the wave function is expressed as an expansion in terms of Slater determinants. The atomic orbital (AO) basis is not well-suited for the description of the effects of the correlation cusp (singularity of the Coulomb operator when the electron−electron distance r12 approaches zero) on the wave function. As a result, the correlation energy converges only slowly with the basis set size, and very large basis sets are needed to obtain accurate results. This can be avoided by taking the electron correlation into account explicitly. This is an approach that dates back to the origins of quantum theory.71 However, for n-electron systems, a complete formulation of the problem requires the evaluation of n-electron integrals, which makes the application of such methods to larger systems impossible. More recently, the explicit correlation effects have been added to standard methods as a correction added to a calculation in an AO basis. The R12 family of methods72 considers only terms linear in r12, covering well the most important short-range effects but becoming unphysical at larger r12. This is corrected in the F12 class of methods,73,74 where the correlation effects are modeled by exp(−γr12) functions (which are, in practice, represented by a Gaussian expansion), γ being a length-scale parameter. 5043


Chemical Reviews

Review

Figure 1. Error of the interaction correlation energy ECCSD(T) and its CCSD(T) correction component δECCSD(T) (a) in absolute value and (b) relative to the average magnitude of the term. The errors are evaluated against values extrapolated to the CBS and averaged over the A24 data set.

in all cases where the extrapolation cannot be performed directly at the CCSD(T) level. Traditionally, the EMP2 term is extrapolated to the CBS; extrapolation from aTZ and aQZ is accurate enough to prevent this step from being the one that limits the overall accuracy. Alternatively, explicitly correlated MP2 can be used to obtain a result close to the CBS limit without extrapolation. Despite the advantages of the explicitly correlated methods, we have found that as soon as the calculations use larger basis sets than double-ζ, both approaches yield similar results (i.e., MP2 extrapolation from aTZ and aQZ basis sets is comparable to MP2-F12 in the QZ-F12 basis) while the extrapolation may be computationally more efficient.64 In double-ζ calculations, however, MP2-F12 yields much more accurate results and MP2-F12/DZ-F12 can be recommended as the most accurate option applicable to large systems where the basis set size is severely limited. For a given size of the system, the step that limits the overall accuracy of the composite calculation is the δE CCSD(T) correction, where one cannot use as large basis sets as in EMP2. For systems with ca. 30 atoms, the calculations in a basis set larger than aTZ become practically impossible. Although it is not possible to eliminate this residual error in calculations of systems of this size, it is necessary to understand and quantify it in order to be able to handle the benchmark CCSD(T)/CBS calculations at this level properly. We have tested the accuracy of composite CCSD(T)/CBS schemes employing various combinations of basis sets on a set of smaller models, where a more accurate estimate of CCSD(T)/CBS interaction energies is possible.64 The main conclusions are presented in this review in section 5.1. In support of calculations of δECCSD(T) in a small basis set, it is often stated that the δECCSD(T) correction is less dependent on the basis set size than the total CCSD(T) correlation energy. The validity of this statement, however, depends on how this effect is evaluated. To clarify this issue, we provide here an analysis performed on interaction energies in the A24 data set59 (see section 4.1 for the description of the set). For a series of basis sets from aDZ to a5Z, we compare the errors of the contribution of the studied quantities, the correction δECCSD(T) and the total interaction correlation energy ECCSD(T), to the interaction energy, taking the value extrapolated from aQZ and a5Z as a reference. The absolute values of the errors are averaged over the data set and plotted in Figure 1a. In absolute terms, the errors of δECCSD(T) are smaller, which makes the curve look flatter, indicating a smaller dependence on the basis set. This is consistent with the first studies on this subject, where this behavior was demonstrated in hydrogen-bonded92

to describe the correlation, and in larger basis sets, the extrapolation does not bring any improvement that would justify the cost of additional calculations.64 2.6. Composite CCSD(T)/CBS Schemes

For benchmarking applications, CCSD(T) is the first method to offer the desired accuracy of the description of electron correlation. To match this accuracy with a comparably small error with respect to the complete basis set limit, the results should be extrapolated to the CBS. Because the first reliable extrapolation can be done from triple- and quadruple-ζ basis sets, the direct extrapolation of CCSD(T) correlation energies is limited to rather small systems. In larger systems, the best estimate of the CCSD(T)/CBS energy can be achieved using composite schemes where the final result is constructed in a stepwise manner from calculations in as large a basis set as possible. The most frequent approach used in the field of noncovalent interactions87−90 is to decompose the CCSD(T) energy ECCSD(T) into three terms: Hartree−Fock energy, EHF; MP2 correlation energy, EMP2; and a higher-order correction, δECCSD(T), often also labeled as ΔCCSD(T) E CCSD(T) = EHF + E MP2 + δE CCSD(T)

(6)

where each contribution can be calculated in a different basis set. The CCSD(T) correction is defined as δE CCSD(T) = E CCSD(T) − E MP2 CCSD(T)

(7)

MP2

where both E and E are calculated in the same basis. The Hartree−Fock energy converges faster with the basis set size, so this term is not a limiting factor. Since the convergence is not as systematic as in the case of correlation energy, extrapolation of this term provides little or no benefits; in practice, EHF is usually calculated in as large a basis set as possible and used without extrapolation. MP2 covers a large part of the correlation energy and the calculations can be performed in rather large basis sets, so that it is possible to estimate the MP2/CBS limit reliably. On the other hand, it is well-documented and theoretically explained that MP2 overestimates the dispersion energy (for a detailed analysis, see ref 91), although it neglects higher-order correlation contributions. This overestimation is the main source of error that has to be corrected by the δECCSD(T) term. Nevertheless, the main advantage of MP2, its computational efficiency, outweights these limitations. Composite schemes based on CBS extrapolation using more expensive methods such as MP3 or CCSD are not more accurate than the ones based on MP2 (see below), which makes MP2 the best choice 5044


Chemical Reviews

Review

and stacked90 complexes. On the other hand, when these errors are made relative to the average magnitude of the respective contribution in the data set (for quantity A calculated in the basis set X, with sharp brackets denoting averaging over the set) rel. error =

⟨|AX − A CBS|⟩ × 100% ⟨|A CBS|⟩

Table 1. Mean Unsigned Error (MUE, in kcal/mol) in the A24 Data Set Calculated for Multiple Composite Schemes Attempting To Estimate the CCSD(T)/CBS Result Extrapolated from aTZ and aQZ Basis Setsa

(8)

the trends, plotted in Figure 1b, look differently. The shape of the curve and the slopes of the segments are similar, indicating that both terms converge with the basis set size similarly (although some deviation might be observed in the aQZ basis). It is interesting to note that the relative errors are also very similar in size in each basis set. This analysis shows that the δECCSD(T) has no special convergence properties and the success of composite CCSD(T)/CBS calculations with δECCSD(T) in a small basis set should be attributed only to the fact that the correction term is about 1 order of magnitude smaller than the contribution of the correlation energy to the interaction energy. In larger molecular systems, the δECCSD(T) term can be calculated in only very small basis sets. In the series of correlation-consistent basis sets, these are DZ and aDZ. Another common choice is the 6-31G**(0.25,0.15) basis (defined in section 2.3), which is considered to be the smallest basis set to allow a reasonable description of noncovalent interactions.2,61 Marshall et al.93 has investigated the reliability of these calculations, showing that in some cases, especially in hydrogen bonds, the δECCSD(T) correction can have the opposite sign from that in larger basis sets and should not be used in these systems. Moreover, due to the errors in the double-ζ basis, extrapolating δECCSD(T) from the aDZ and aTZ basis sets in some systems introduces a larger error than the one obtained using δECCSD(T)/aTZ alone. It is possible to add one more step into eq 6, calculating separately ECCSD in a larger basis than the one used for EMP2 but smaller than the one in ECCSD(T) and possibly also extrapolating this term. It is closely related to a third option, the scheme built from EHF, ECCSD extrapolated to the CBS (or calculated with explicit correlation), and the contribution of triples, δE(T) = ECCSD(T) − ECCSD, calculated separately in a smaller basis (or without explicit correlation). Although these schemes are computationally more expensive, they are not necessarily better. The success of the scheme described above stems from the fact that the difference between EMP2 and ECCSD(T) is small, so any error in this term has only a small impact on the final result. Although CCSD is much more computationally demanding and recovers more correlation energy, it does not yield better interaction energies than MP2.94 The δE(T) correction is thus larger and the same relative inaccuracy in this term introduces a larger error into the composite result. This can be demonstrated by evaluating the errors of these schemes in the A24 data set. Taking CCSD(T)/CBS results extrapolated directly from CCSD(T) calculations in aTZ and aQZ basis sets as a reference, the errors of various less expensive composite schemes were evaluated. All the tested schemes are constructed from a baseline MP2 or CCSD calculation extrapolated to the CBS limit from aTZ and aQZ basis sets to which further corrections in smaller basis set (aDZ, aTZ) are added. The baseline level and basis sets used for the correction are listed in Table 1 along with the error with respect to the benchmark. Note that we intentionally omit calculations in larger basis sets, where the errors will be less

scheme

baseline

correction

MUE

1 2 3 4 5

MP2/CBS CCSD/CBS CCSD/CBS MP2/CBS CCSD/CBS

δECCSD(T)/aTZ δE(T)/aTZ δECCSD/aTZ + δE(T)/aDZ δECCSD(T)/aDZ δE(T)/aDZ

0.011 0.019 0.039 0.038 0.063

a

The same extrapolation is applied at the baseline level, which is subsequently corrected by adding further terms calculated in smaller basis sets.

significant. The schemes considered here are also the ones that would be applicable to larger systems. The data have been taken from the calculations presented in ref 59; the computational details are described therein. When a single correction is used, the schemes based on MP2/CBS with a δECCSD(T) correction (schemes 1 and 4) outperform their analogs based on CCSD/CBS with a δE(T) correction (schemes 2 and 5). The three-step scheme 3 has almost the same error as the cheaper, two-step one based on MP2/CBS and CCSD(T) calculations in the same basis. These observations can be further justified by comparing the average magnitude of the corrections, which are, in the aTZ basis set, 0.11 kcal/mol for δECCSD(T) and about twice as large, 0.23 kcal/mol, for δE(T). We also tested the related schemes based on explicitly correlated calculations. For this analysis, MP2-F12 and CCSD(T)-F12 calculations of the A24 data set were performed using the Molpro package,95 with the setup recommended for noncovalent interactions57 (using the CCSD-F12b approximation, 3CFIX ansatz, and a geminal exponent of 0.9 for DZ-F12 basis and 1.0 for larger basis sets). Since the -F12 basis sets had not been available for argon, the two complexes containing an argon atom were excluded from this analysis. All the calculations start with the MP2-F12 component calculated in the QZ-F12 basis set, which ensures results close to the CBS limit practically equivalent to extrapolation from aTZ and aQZ basis sets. The δECCSD(T) correction is then calculated either using the CCSD(T)-F12 method (optionally with the empirical scaling of the triples, denoted T*) or with conventional CCSD(T) in aDZ basis set. All calculations were counterpoisecorrected. CCSD(T)/CBS interaction energies extrapolated from aTZ, aQZ, and a5Z basis sets are used as a reference. The best results with an average relative unsigned error of 1.95% is achieved when the correction is calculated using the conventional CCSD(T)/aDZ method. The error grows considerably when we pass to the CCSD(T)-F12 method, where it amounts to 6.38%. This is again caused by breaking the balance between the CCSD contribution (that is calculated using explicit correlation, which is equivalent to the use of a larger basis set in a conventional calculation) and the triples (that are not explicitly correlated). Scaling the triples empirically in the CCSD(T*)-F12 approach improves the results only slightly to 3.66%. Despite their simplicity, the composite schemes based on a higher-order correction to MP2/CBS results represent the most accurate approach. This is caused by a favorable compensation of errors that makes the correction small and thus less sensitive to basis set size (in absolute terms). When this compensation is broken in approaches based on more 5045


Chemical Reviews

Review

basis in polar complexes such as hydrogen bonds. On the other hand, nonpolar, dispersion-bound complexes are described reasonably even in the 6-31G**(0.25,0.15) basis. This correction can also be calculated in two steps, first passing from CCSD(T) to CCSDT and then from CCSDT to CCSDT(Q), calculating CCSDT in an intermediate basis set. We have shown that the results improve slightly but the improvement is not significant enough to justify the additional computational expenses.105 (2) The frozen-core approximation is used in most correlated calculations of noncovalent interactions. In calculations of the interaction energy, a large part of the error in the correlation energy caused by this approximation cancels. The correlation of all electrons is straightforward but makes the calculation more expensive. To cover this effect fully, the basis sets should have a balanced description of the correlation for all the electrons, which is achieved by adding polarization functions to the core orbitals. The polarized-core correlation-consistent (ccpCVXZ) 51 and improved weighted polarized-core (ccpwCVXZ)52 basis sets are available for this purpose. (3) At this level of accuracy, relativistic effects cannot be neglected, even in light atoms. Among the computationally affordable relativistic methods, the Douglass−Kroll−Hess (DKH) Hamiltonian107−109 is used the most widely. It can be applied at the CCSD(T) level, and a specialized variant of the correlation-consistent basis sets is available for DKH calculations.110 This is the setup used for the A24 data set. The relativistic effects can be evaluated more rigorously using the four-component Dirac−Coulomb Hamiltonian, which can also be combined with CCSD(T), but these calculations are substantially more expensive. We have recently performed these calculations for the HF dimer.96 (4) Practically all electronic structure methods rely on the Born−Oppenheimer approximation. However, the errors associated with it are similar in size to the other corrections discussed here. The first correction covering the electron− nuclei coupling is the diagonal Born−Oppenheimer correction (DBOC).111 It is implemented for arbitrary level CC calculations.112 The calculation of the DBOC is based on the calculation of the second derivatives of the energy and is thus substantially more demanding than single-point calculations. The effect of DBOC on interaction energies has not been studied systematically in a larger number of complexes; in the HF dimer, it amounts to −0.012 kcal/mol when calculated at the CCSD level, compensating well for the relativistic correction, which has the value of 0.016 kcal/mol.96 This approach of correcting the baseline CCSD(T)/CBS result with all available corrections calculated at the best level and in the largest basis sets possible is closely related to the composite schemes used in computational thermochemistry.113 Among these, the focal point analysis,114 HEAT,115−117 W3,118 and W4119 protocols use multiple post-CCSD(T) terms. These schemes are, however, tuned for calculation of the total energies of isolated molecules where the relative importance of the individual contributions is different than in the calculations of interaction energies in noncovalent complexes. This is why the composite methods used for noncovalent interactions are usually constructed independently rather than by following the protocols used for thermochemical quantities. In the realm of noncovalent interactions, composite post-CCSD(T) calculations have been performed for small hydrogen bonds120 and the lithium−thiophene complex,121 as well as for small dimers, for which accurate experimental data are available.96,122

accurate CCSD calculations (either using larger basis set or explicitly correlated) complemented with triples calculated in smaller basis set, the error gets substantially larger. 2.7. Post-CCSD(T) Calculations

If an even higher accuracy is sought, it is possible to go beyond the CCSD(T)/CBS gold standard. There are two main reasons to perform such calculations. First, they become necessary for reproducing the most accurate experimental results (dissociation energies), which are available with errors as low as a few cm−1 (1 cm−1 = 2.859 cal/mol). Second, such calculations performed in model systems serve as a valuable benchmark that makes it possible to analyze and quantify the error in CCSD(T)/CBS calculations. In the standard CCSD(T)/CBS calculation, there are multiple approximations that can be lifted. Since all these improvements are associated with additional computational costs, which prevent the use of larger basis sets, these terms are most often evaluated individually and added as corrections to a well-converged CCSD(T)/CBS result. Here, we list these corrections sorted by decreasing importance (as evaluated in refs 59 and 96). (1) The most obvious approximation is the truncation of the coupled-cluster series. It is possible to add higher excitations, and modern computational techniques based on a computergenerated code now allow arbitrary-level CC calculations up to full configuration interaction97−99 (where the latest version developed by Kállay and Surján100 does not introduce any overhead in terms of scaling). In general, the computational complexity of CC calculations scales as onvn+2, with o being the number of occupied orbitals, v the number of virtual orbitals, and n being the order of the highest excitation operator. CCSDT thus scales with o3v5, CCSDTQ with o4v6, and so on. It is clear that these calculations can be performed only on small molecular systems and in small basis sets. In addition to the smallest noncovalent complexes, rare gas dimers,101−103 coupled cluster calculations up to CCSDTQP, and a full configuration interaction in a 6-31G**(0.25,0.15) basis were applied to small molecular complexes containing up to two first-row atoms saturated with several hydrogens.16 Both CCSDTQP and its approximate variant CCSDTQ(P) were found to be practically converged with respect to FCI. CCSDTQ(P) can be applied to larger systems (with up to about 20 electrons) and was used as a benchmark for CCSDT(Q)18 and CCSDTQ. Both of these methods exhibited only negligible errors below 1 cal/mol (0.1% of the interaction energy), while the errors of CCSDT and CCSD(T) were 0.5% and 1.6%, respectively. This is consistent with the previous findings that CCSDT(Q) recovers about 95% of correlation energy missing in CCSD(T).18 The triples and quadruples can both be treated perturbatively in the CCSD(TQ)104 method, but our results suggest that this approach does not provide sufficient improvement over CCSD(T) and yields results worse than CCSDT.16 This analysis shows that CCSDT(Q) is accurate enough for most applications, yet it can be applied to somewhat larger systems and combined with larger basis sets that are necessary to obtain quantitative results. We have applied CCSDT(Q)/ aDZ to the A24 data set,59,64 where the largest system calculated at this level is ethane dimer. The dependence of the results on the basis set size should be examined with special care. It has been shown105,106 that it is impossible to calculate reliably not only the values but also the sign of the correction calculated as ΔECCSDT(Q) − ΔECCSD(T) in the aDZ or smaller 5046


Chemical Reviews

Review

2.8. Approximate CCSD(T) Approaches

The practical value of these approaches in calculations of noncovalent interactions is difficult to establish. The local methods target at larger systems where no reliable benchmark is available. To bridge the gap between small model systems and the large ones, we have analyzed the trends in accuracy in a series of systems of systematically increasing size. The results are discussed in section 5.3.

To increase the efficiency of CCSD(T) calculations, various approximations can be employed at this level. They range from the mild ones, which provide some acceleration without compromising the results, to very approximate calculations, which are less accurate but applicable to much larger systems than the canonical CCSD(T). In the latter case, however, it has to be carefully considered whether some other methodology with more favorable scaling may yield results superior to the approximate CCSD(T). There are two main approaches to approximating full CCSD(T) calculations. First, it is possible to truncate the virtual space (the set of unoccupied MOs), lowering the number of excitations used in the CC calculation. Although this does not change the unfavorable scaling, the computational cost can be substantially reduced. As CCSD(T) scales with a number of virtual orbitals v as v4, the computational effort can be reduced by (vtruncated/ vfull)4. To make this procedure efficient while conserving as much accuracy as possible, the active space has to be defined carefully. It has been shown that passing to MP2 natural orbitals allows better separation of the virtual orbitals that are important from the ones that can be discarded. This is the basis of the frozen natural orbital (FNO) approximation.123 Here, the orbitals are selected on the basis of their occupancy, with a threshold of 10−5 being suitable for noncovalent interactions.124 A more advanced approach is the optimized virtual orbital space (OVOS) method,125 where the truncated wave function is optimized iteratively. The threshold for discarding VOs can be varied automatically to achieve a specified target accuracy.126 Given the nature of these approximations, they are most efficient in large basis sets, while in small (double-ζ) basis sets, practically no virtual orbitals fall below the threshold. On the other hand, the error associated with the FNO approximation slightly grows with the basis set size, and we have shown that in composite CCSD(T)/CBS calculations of counterpoisecorrected interaction energies, it becomes equally important as the error caused by the basis set incompleteness when the CCSD(T) correction is evaluated in the aQZ basis.64 These approximations do not extend the applicability of CCSD(T) to significantly larger systems; nevertheless, since it is possible to achieve very small errors, they can be used to accelerate the calculations in most benchmarking applications. Additionally, in the calculations of noncovalent interactions, these methods provide one very important advantage. When the counterpoise correction is used, the virtual space originating from the ghost basis functions in the calculations of the monomers can be substantially reduced without compromising the accuracy.127 Second, it is possible to introduce the local approximation.128,129 Many variants of this approximations have been developed for CC methods,130−135 but they share common principles: the orbitals are localized and the correlation is evaluated only between those that are spatially close. While this approach can be successful in recovering a major part of the total correlation energy, its application to noncovalent interactions is more difficult. By definition, noncovalent interactions are not local, and they represent only a small fraction of the total correlation energy. Therefore, very tight thresholds have to be used, and even then, errors of about 0.1 kcal/mol cannot be avoided.134,136 On the other hand, these methods do actually reduce the scaling, in some cases up to almost linear in very large systems (where each local domain has only a fixed number of neighbors).

2.9. Quantum Monte Carlo

Recent advances in the Quantum Monte Carlo (QMC) methodology based on the fixed-node diffusion Monte Carlo (FN-DMC) have made this method useful also for benchmark calculations. Details can be found in a separate review in this issue by Dubecký et al.;137 here we will only summarize the most important applications of the method related to this review. When compared to post-HF correlated methods, the QMC approach has both advantages and limitations. While the statistical error inherent to Monte Carlo can be converged to a sufficiently small value, there are systematic errors originating from the approximate form of the trial wave function, most importantly from the use of fixed-node approximation. On the other hand, the QMC describes the electron correlation at the full CI limit. The calculations are demanding, but their scaling is about O(N)3 to O(N)4, which makes it feasible (but timeconsuming) to apply them to large systems. Another advantage is that the QMC algorithms can be parallelized efficiently even in massively parallel environments. Two interesting applications of QMC methodology to noncovalent interactions have recently emerged. First, it is the development of the methodology toward increased accuracy. Smaller systems can now be calculated with a degree of accuracy comparable to benchmark WFT calculations.138 The observed agreement between FN-DMC and higher-order coupled cluster calculations provides valuable verification of the accuracy of both approaches.64 Second, the QMC methodology has been applied to some very large noncovalent complexes with more than 100 atoms.139 In systems of this size, the QMC may be the least approximate method applicable. 2.10. Wave Function Methods

We will demonstrate the application of some of the benchmark data sets to the evaluation of the accuracy of more approximate correlated wave function methods. We list them here, putting them into the context of related methods not covered in the Important Results section. As the simplest correlated method, MP2 is widely used in the field. Although MP2 describes polar complexes and hydrogen bonds very well, it is notorious for overestimating London dispersion, yielding rather large errors. This is caused by the treatment of polarization and dispersion at the uncoupled Hartree−Fock level.91 This can be partially compensated for by using smaller basis sets, where the overestimated dispersion is compensated for by the error caused by the basis set incompleteness. Various possible setups have been investigated; for general use, the MP2/TZ combination can be recommended.62 In addition to canonical MP2, we test its variants aimed at improving the description of noncovalent interactions. Since the main issue in MP2 calculations is the overestimation of dispersion, the attempts to improve MP2 are based on scaling these effects down more or less empirically. Splitting the MP2 correlation energy into its spin components provides more room for empirical scaling in the 5047


Chemical Reviews

Review

MP2.5 is very robust and performs equally well in a wide range of systems,151 including three-body interactions.152 It can be recommended for all applications where high accuracy is needed but more expensive CCSD(T)- or CCSD-based calculations are not affordable. The computational time required for the rate-limiting step, the MP3 calculation, corresponds roughly to one CCSD iteration (about 15 iterations are usually needed to achieve the convergence in CCSD calculations). Moreover, the MP3 calculation can be substantially accelerated by using the RI approximation. More empiricism can be introduced by optimizing the ratio of MP2 and MP3 contributions, resulting in the MP2.X method.63 In larger basis sets, however, the optimal ratio is very close to 0.5 and the negligible gain in accuracy does not justify abandoning the elegance of MP2.5. In composite calculations where the correction to MP2/CBS is calculated in a small basis set (roughly smaller than aDZ), this optimization can substantially improve the accuracy and provides a more balanced description of different classes of interactions.63,151 The last step in the series before reaching CCSD(T) are methods based on CCSD. The scaling of CCSD, N6, is the same as of MP3, but CCSD is iterative and thus substantially more expensive. CCSD itself performs rather poorly, with errors often larger than MP2. This can be improved by the scaling of the spin components. The general parametrization, SCS-CCSD153 with cos = 1.27 and css = 1.13, improves the results significantly, but it does not outperform the cheaper MP2.5. On the other hand, after the scaling coefficients were optimized for molecular interactions in the SCS-MI-CCSD method,154 excellent results were achieved. The optimization of the scaling coefficients was performed on the S22 data set, leading to the values cos = 1.11 and css = 1.28. SCS-MI-CCSD/ CBS results closely follow CCSD(T)/CBS; e.g., in the S66 data set the RMSE is below 0.1 kcal/mol. Again, the high accuracy of SCS-MI-CCSD was also confirmed in calculations of threebody interactions.152

so-called spin-component-scaled (SCS) approaches. The original SCS-MP2 method140 has been developed to reproduce reaction energies, but it improves also the description of noncovalent interactions. The scaling coefficients applied to the other-spin (singlet) and same-spin (triplet) terms are cos = 6/5 and css = 1/3. This approach has been reoptimized specifically for interaction energies (in the S22 data set), forming the SCSMI-MP2 method141 (where MI stands for molecular interactions). The coefficients have been optimized for particular basis set combinations; in this work, we use the one developed for results extrapolated from TZ and QZ basis sets (cos = 0.4 and css = 1.29) while we extrapolate from aTZ and aQZ (after validating that the difference is negligible94). Furthermore, some computational time can be saved when only one spin component is evaluated, as in the SOS-MP2 method (where cos = 1.3).142 All of these approaches have recently been reviewed by Grimme et al.143 The spin-component scaling has been further developed into the dispersion-weighted MP2 (DWMP2),144 which interpolates between MP2 and SCS-MP2 using a system-dependent factor calculated from the ratio between the MP2 and HF contributions. The empirical nature of these approaches, combined with the inherent limitations of MP2 theory, lead to rather unbalanced results. None of the SCSbased methods yields a satisfactory description of all kinds of noncovalent interactions, and even in the case of SCS-MI-MP2, the error remains rather high, even in the case of simple organic systems it was fitted on. Alternatively, the MP2 energy can be divided into short- and long-range contributions by means of introducing a rangeseparated Coulomb operator. By keeping only the short-range component and optimizing the range-separation function, the attenuated MP2 method has been obtained.145,146 It has been developed specifically for calculation in small basis sets (aDZ and aTZ), where this approximation compensates for the error of the basis set surprisingly well. A more rigorous approach is to remove the problematic component, uncoupled HF dispersion, from MP2, replacing it with the dispersion calculated using time-dependent DFT (TDDFT). This is used in the MP2C (C stands for “coupled”) method,147,148 which performs the best among the MP2-based approaches. Another deficiency of MP2 is observed in larger π−π complexes, where MP2 diverges because of the narrowing HOMO−LUMO gap. This has been addressed by the restrained-denominator MP2 (RD-MP2) approach,149 which also halves the error of MP2 in the S66 data set. The next step in perturbation theory, MP3, does not bring much improvement over MP2 when interaction energies are considered. The small improvement in accuracy cannot justify the increase of scaling from N5 to N6. In contrast to MP2, MP3 underestimates the strength of the interaction. Interestingly, the magnitude of these errors correlates well across systems of different nature. This observation has led to the construction of the MP2.5 method,150 where the correlation energy is taken as the average of MP2 and MP3 energies. To achieve the CBS limit, the MP3 contribution can be added to MP2/CBS analogously to other composite schemes: E MP2.5 = EHF + E MP2/CBS + 0.5(E MP3 − E MP2)

2.11. DFT and Semiempirical Methods

While both of these topics are covered by separate reviews in this issue, we are including several widely used methods to put the WFT methodology into a broader perspective. It is not intended as a complete survey of the methods available. For more details on the DFT and semiempirical methods, we refer to the respective reviews in this issue. The main deficiency of the classical DFT methodology in the field of noncovalent interactions is its inability to describe London dispersion. While DFT functionals cover some part of correlation energy, they do not cover the long-range effects responsible for dispersion. Here, we will discuss only one of the possible solutions, adding the dispersion by means of an empirical correction independent of the DFT calculation in the so-called DFT-D methods.155−157 This approach is very successful for two reasons: First, the dispersion is rather well separated from the DFT energy, and, second, even a simple pairwise empirical model can describe the missing dispersion with a high degree of accuracy. Among the DFT-D approaches, the DFT-D3 method of Grimme et al.158,159 is accurate, robust, and applicable to a wide range of systems. When combined with well-chosen DFT functionals and calculated in a large def2-QZVP basis set, this method yields very good agreement with the CCSD(T)/CBS benchmark data.160 It has recently become available in many DFT codes, which makes it the first choice for general-purpose applications. The DFT-D method-

(9)

Even though this is only an empirical approach (since no theory guarantees that the errors of MP2 and MP3 are of similar size), it works very well and it has no competition at this complexity level. In subsequent studies, it has been found that 5048


Chemical Reviews

Review

accurate description of the exchange repulsion. For a recent review of this topic, see ref 181. The main advantage of the SAPT methodology, compared to the variational methods, is its ability to yield not only the total interaction energy but also its decomposition into physically well-defined terms. DFT-SAPT is usually applied up to the (1) second order, providing separate electrostatic (Eelst ) and (1) exchange (Eexch) terms in the first order and induction (E(2) ind ) and dispersion (E(2) disp) in the second. The second-order terms are both associated with complementary exchange contribu(2) tions (E(2) exch−ind) and (Eexch−disp). To reach quantitative accuracy, the higher-order induction terms are included in the form of the δHF correction calculated from the difference between Hartree−Fock interaction energy and the corresponding SAPT terms. DFT-SAPT calculations in large basis set can be very accurate, approaching the accuracy of supramolecular CCSD(T) calculations.182,183 In the S22 data set, DFT-SAPT calculations extrapolated to CBS from aDZ and aTZ basis sets have a RMSE of 0.49 kcal/mol with respect to the CSCD(T)/CBS reference;6 part of this error could be attributed to the relatively small basis sets used. From a computational point of view, the DFT-based SAPT calculations are applicable to rather large systems; SAPT(DFT) calculations in aDZ basis set augmented with midbond functions were recently reported for coronene dimer.184

ology and related methods are reviewed separately by Grimme et al.161 There exist further approaches to correcting DFT methods, including nonempirical corrections and explicit treatment of electron correlation within the DFT functional. Among the latter, the random phase approximation (RPA) improves the description of noncovalent interactions substantially while very favorable scaling is retained.162 The situation is similar in semiempirical quantum-mechanical (SQM) methods (see the review by Elstner and co-workers163 in this issue), where the dispersion is missing but can be rather easily added empirically. It was done first in the self-consistent density-functional tight-binding (SCC-DFTB) method164 and later applied also to the modified neglect of diatomic overlap (MNDO) type methods.165−170 The latest development was porting the successful D3 dispersion correction to SQM methods.170 Additionally, the combination of the approximations involved in SQM methods leads to a severe underestimation of hydrogen bonds. We have pioneered the use of an additional empirical correction addressing this problem,168 and the last version, combined with the dispersion correction in the D3H4 methods,170 yields excellent results with errors below 1 kcal/mol for all types of interactions in organic molecules. The parametrization of these corrections is an important example of the use of large databases of benchmark data and is described in detail in section 5.4. In some applications, these advanced SQM methods are becoming a viable replacement for empirical force fields. Recent developments in linear scaling algorithms make it possible to calculate thousands of atoms at the time scale of seconds or minutes.171 Because these methods are based on the proper QM description of the system, they can naturally describe effects that are difficult to model in force fields. This applies not only to noncovalent interactions (e.g., charge ansiotropy, polarization and charge transfer) but also to chemical reactions. Unlike force fields, SQM methods do not require any systemspecific parametrization what makes the preparation of the calculations straightforward.

3. DATA SET CONSTRUCTION AND PROCESSING 3.1. Data Set Composition

In this part of the review, we would like to summarize our experience with building data sets of benchmark interaction energies. Constructing a useful data set is not just a matter of running expensive calculations; special care must be paid to the composition of the set. We will discuss the construction of a data set with a broader scope, targeted not only the validation but also the parametrization of approximate QM methods. The prime example of this approach is our S66 data set,94 aimed at covering the most common interaction motifs in organic molecules and biomolecular building blocks. This information is important not only for the construction of new data sets but also for a conscious interpretation of the results obtained with existing data sets. The first important parameter is the size of the data set. Here, two main requirements oppose each other, and a suitable compromise must be found. A general-purpose data set must be large enough to enable reliable statistical processing of the results. This means that for each type of interaction, several model systems should be included. These model systems should not be too similar so as to be redundant. On the other hand, the data set should not be excessively large. The first reason for that is to make the use of the data set practical and inexpensive. This applies not only to the computational side but also to working with the results. Using a very large, broad set of data in a parametrization has some benefits, but isolated, yet important, errors might be lost in the amounts of data that make it difficult to manually inspect and understand the processes involved. The second reason is that at some point, the coverage of the desired chemical space becomes sufficient and the addition of more similar model systems does not bring any further information. We have found that for most applications, the lower limit is about 20 or 30 systems and the upper limit is 100. Specific conditions and the focus on one particular issue may warrant the development of smaller data

2.12. Intermolecular Perturbation Theory and Interaction Energy Decomposition

Another approach to calculation of intermolecular noncovalent interactions is the application of perturbation theory where the interaction is represented as a perturbation to a Hamiltonian constructed from the Hamiltonians of the noninteracting monomers. This is a valid approach at large intermolecular distances but it cannot be applied directly in the overlap region where the wave function constructed from the fragments violates the Pauli exclusion principle. To address this, projections that enforce the antisymmetry of the wave function had been introduced in the development of symmetry-adapted perturbation theory (SAPT) (for review of this topic, see refs 172 and 173). To describe noncovalent interactions accurately, this approach has to be based on a Hamiltonian covering correlation energy, which was achieved by using the coupledcluster formalism.174 This leads to very complicated tripleperturbation formalism where the dispersion energy term scales with N7, making this method applicable only to small systems. It is also possible to start from the DFT Hamiltonian, which covers the correlation effects in the subsystems but is much simpler to handle.175 This method was developed independently as DFT-SAPT176−178 and SAPT(DFT).179,180 These approaches are practically identical, and in both cases, the use of an asymptotically corrected DFT functional is necessary for 5049


Chemical Reviews

Review

sets. If multiple data are needed, it is better to split them into more separate, well-defined data sets that can be also used independently. Each data set should have a well-defined scope, to provide a reasonably complete coverage of a specific chemical subspace. This has to be realized at three levels: First, the monomers are chosen that comprise the complexes. Second, the complexes should be built from these monomers so that all combinations resulting in meaningful interactions are sampled. Third, and most importantly, the interactions should be classified by their nature (see section 3.4) and the data set should cover all classes of interactions in a balanced way. Attention should be paid not only to the counts of systems representing some class of interactions but also to their energetic contributions. If possible, the interaction energies in the set should lie in a not too wide range so that measuring the error across the whole data set provides reliable and balanced information on all the systems. Systems interacting too weakly do not contribute to the overall statistics. This can be partially remedied by using relative errors, but when the interaction energy is too close to zero, the relative error becomes unstable and might be too high, even though the absolute error is insignificant (see also section 3.5). On the other hand, systems with too strong interaction may override the remainder of the set. This balance is, however, difficult to achieve in practice. For example, in the S66 data set, about a third of the data set is devoted to hydrogen bonds. These interactions are substantially stronger than the dispersion-dominated and nonspecific interactions in the remainder of the set. It is something to be kept in mind, but in this particular case, it is not a problem, because in most applications, the errors are distributed more uniformly than the interaction energies themselves. This is because most methods describe hydrogen bonding better than dispersion. Additionally, when we calculate the energy balance not as sums of interaction energies in the groups but more physically as a sum of the interaction energy components in the whole set, the results are more balanced. Using the DFT-based symmetry-adapted perturbation theory (DFT-SAPT) to decompose and sum electrostatic and dispersion contributions across the set, their ratio is 1:0.86.94

out in a smaller (cc-pVTZ) basis set, using counterpoise correction. Recent advances in DFT-D methodology make it a viable tool for this purpose as well; geometries obtained with DFT-D3 in a large basis set are on par with or even better than the ones from MP2.160 In a smaller system, it is possible to apply post-MP2 methods including CCSD(T)/CBS also to geometry optimizations. This has already been applied to some data sets of small systems.59,185 A data set of geometries complemented by harmonic vibrational frequencies calculated at the CCSD(T) level is also available.186 The harmonic vibrational frequencies may serve for benchmarking other methods; for accurate reproduction of experimental results, anharmonic calculations are needed. Noncovalent interactions, and dispersion in particular, are, among the intermolecular degrees of freedom, the most sensitive to the distance. This can be exploited for refining the geometry at a higher level, performing a scan in the distance coordinate, and interpolating a minimum on that curve. This is applicable even in larger systems, where direct optimization at the post-MP2 level would not be possible. This approach was used in the construction of the S66 and X40 data sets,94,187 where the geometries were constructed as a minimum of the distance scan at the CCSD(T)/CBS level. The same approach was also used for analyzing the optimal intermolecular distances obtained by interpolation in various post-HF methods and comparing them to CCSD(T)/CBS.188 This work suggests that MP2.5/CBS is a promising method for obtaining accurate geometries at a not too high cost, with the error being reduced to one-third as compared to the best MP2 results. Another possible source of geometries is experiment. In the gas phase, the experimental data can be used only indirectly in modeling the geometries. A more useful source of geometries of noncovalent complexes are crystallographic data. Multiple dimer or cluster geometries with different properties can be extracted from the crystal. Although these do not correspond to the gas-phase equilibrium, they provide a realistic sampling of the interactions in the condensed phase. For the recent 3B-69 data set of three-body interactions,152 we have used geometries extracted from DFT-D-optimized crystal geometries, extracting three different conformations from each crystal.

3.2. Geometries

3.3. Nonequilibrium Geometries

When only interaction energies are concerned, the choice of the geometries of the complexes is not very important. The only requirement is that the covalent structure of the monomers should not be too distorted for their electronic structure to represent well their ground state. On the other hand, using nonarbitrary geometries is often preferred, because the choice of specific geometries can add further value. The first option is to use optimized geometries, either global minima or the local ones featuring the desired interaction motif. This has several advantages. It makes it possible to generate consistent and reproducible geometries, and if the method used for the geometry optimization is accurate enough, the geometries themselves can serve as a benchmark. Additionally, it maximizes the strength of the interaction. This is the approach used in the majority of the data sets available. For noncovalently interacting systems, it is important to cover London dispersion properly even at the level of the geometry optimization. Among correlated WFN methods, the cheapest option, MP2, is used most often and it is the only option for larger systems. To remediate the problems of MP2 in strong dispersion-bound complexes, the optimization is often carried

While single-point interaction energy calculations in equilibrium geometries are sufficient for benchmarking more accurate methods that can be expected to behave similarly over the entire potential energy surface, working with more approximate methods requires extending the benchmarking also to nonequilibrium geometries. This is especially important in the parametrization of highly empirical methods, such as SQM methods, and the corresponding corrections for noncovalent interactions. These corrections are often completely empirical and contain a rather large number of parameters. To ensure that the correction is valid for any geometry, extended sampling of the potential energy surface is required. The use of a benchmark data set in the parametrization of corrections for SQM methods is discussed in detail in section 5.4. A complete characterization of the intermolecular PES at the benchmark level is very expensive and thus applicable only to very small systems. Moreover, the majority of the generated points would be redundant or have a low relevance for the effects studied. The PES can be sampled randomly to reduce the overall number of points, but the other disadvantages remain the same. Adding an energetic criterion to the random 5050


Chemical Reviews

Review

groups: (1) hydrogen bonds; (2) dispersion-dominated complexes; and (3) “mixed” or “other” systems, where the dispersion is complemented by an electrostatic contribution that is not as localized as in H-bonds; it is usually a dipole− dipole interaction. Among the dispersion-dominated complexes, it is useful to separate the π−π (stacking) interactions because of their specific properties. Similarly, more specific classes of interactions can be identified in a broader chemical space. In data sets mixing neutral and charged systems, it is useful to separate them into groups by charge, as their description by approximate methods may also differ. Moreover, the systems with significant intermolecular charge transfer are often set apart. This categorization of the interactions can also be done more rigorously by means of a decomposition of the interaction energy. For example, the categories described above in the S66 data set were compared with the sorting of the systems by the ratio of dispersion and electrostatic (the sum of Coulombic and induction) contributions from DFT-SAPT calculations. Using a single empirically adjusted threshold, it was possible to reproduce these three categories rather well.94 The interaction energy decomposition and various relative measures based on it are also useful in the analysis and interpretation of the results without sorting them into the groups. The errors can be correlated with and plotted against these measures, revealing important trends. For obtaining these measures, various interaction energy decomposition methods can be used. Among them, the SAPT and especially its version based on Kohn−Sham orbitals, DFT-SAPT or SAPT(DFT) (discussed in section 2.12), provides a well-defined decomposition and accurate results. However, even simpler approaches might be useful: in the study of three-body interactions where rigorous decomposition is not tractable, we have used polarization (Epol) and dispersion (Edisp) contributions from a QM-derived force field.152 However, any relative measures used for this quantification have to be designed carefully. If the values of some of the components are close to zero, a simple ratio of the values is not reliable. Instead, for robust quantification of the dispersionelectrostatic scale, we have defined152 the fraction of the dispersion as

sampling, e.g., using snapshots from molecular dynamics, improves the coverage of the important regions of the PES. Such a random sampling is useful for some specific purposes, including very thorough testing of the robustness of a method. For most uses, however, it is more useful to sample the PES systematically in a more controlled way, isolating the main features. The most important coordinate is usually the intermolecular distance. For our S66 and X40 data sets, we have developed the following protocol to ensure consistent results: We scale the closest intermolecular distance by a series of fixed factors, moving the molecules along the axis of the interaction. This axis is chosen manually in order to conserve the character of the interaction; e.g., for hydrogen-bonded systems, it is the vector corresponding to the hydrogen bond, and for nonspecific interactions, the vector connecting molecular centers of mass can be used. For reliable parametrization of approximate methods, we found it necessary to start the sampling already at repulsive distances. In the newer X40x10 data set,187 we start with a scaling factor of 0.8; in S66x8,94 the dissociation curves start at 0.9. The region around the expected minimum is sampled more densely and the curves end at 2 times the closest distance in equilibrium. Using 8 or 10 points along the curve provides more than enough data for accurate interpolation. Some interactions, such as hydrogen and halogen bonds, have specific angular dependence that should be examined in method validation and reproduced in parametrization. Covering the angular dependence in detail would produce a very large data set with a high level of redundancy. Since these effects are usually well transferable between similar systems, in most cases we opt for working with a few specific model systems. For less detailed verification, we have extended the S66 data set with angular-displaced geometries, building the S66a8 data set.189 Here, each of the molecules in the complex is rotated around two axes perpendicular to the intermolecular axis by 30° in each direction, which gives eight points for each system. In general, the nonequilibrium geometries sampling the surroundings of the energy minima as described above (a combination of radial and angular displacements) provide all the information necessary to locate the minimum. If the tested method reproduces all these points well, it should also be able to predict the benchmark geometry with good accuracy. Although this approach cannot fully replace the benchmarking carried out on geometries optimized by the method tested, it can provide a quick overview at a substantially lower cost. Additionally, parts of the PES not visited by the optimizations but important for modeling realistic systems (such as the sampling of the short distances that are crucial in compressed condensed systems) are covered as well.

%dispersion =

|Edisp| |Edisp| + |Epol|

× 100 (10)

3.5. Statistical Processing of the Results

Multiple statistical measures can be used for the evaluation of an error in a data set or its subsets. Here, we will overview the widely used options. Each of them carries different information. For a thorough comparison of the benchmarked methods, more of them should be listed. If a single error measure has to be chosen, the root-meansquare error (RMSE) is the best candidate (it is also known as the root-mean-square deviation, RMSD). In the case of energies E compared to a reference Eref in a data set of size N, it is defined as

3.4. Categorization of Interaction Types

For interpreting the results in data sets featuring different types of interactions, it is useful to split the data set into groups of systems of similar character. A comparison of the errors between the groups provides valuable information on the strengths and weaknesses of the method tested. Additionally, quantifying the nature of the interactions in the set makes it possible to correlate the results with the character of the interaction. The simplest approach is grouping the systems by their character assigned on the basis of “chemical intuition”. In neutral, organic molecules, such as in the S22 and S66 data sets, the complexes have been divided into three easily recognizable

RMSE =

1 N

N

∑ (Ei − Eiref )2 i=1

(11)

Since the RMSE is based on the square of the differences, it combines two important aspects of the error: its average size and its extreme values. It is thus a good measure of both the 5051


Chemical Reviews

Review

is the Psi4 package,190 which contains multiple data sets in a form that makes it possible to run a calculation over the data set simply by selecting the method and the data set to be used. For our work, one of the authors (J.Ř .) has developed a more universal piece of software, which is now available to the public. The Cuby framework191,192 provides an interface to various computational chemistry software packages, giving the user access to a wide range of methods from molecular mechanics to advanced quantum chemistry. Although it is not intended solely for this purpose, automation of complex workflows on data sets is one of the important goals of the software. It contains a library of data sets (and the option to use user-defined ones); the calculations on a predefined data set can be called very easily just by specifying the method and the name of the set. The list of the data sets available can be found in the online documentation.191 The individual calculations can run in parallel either on a single machine or over the network in a cluster. After the calculation, the basic statistical analysis of the results is performed automatically, including individual statistics for each predefined group of systems. Moreover, the framework also automates basis set extrapolations and the composite calculations based on them.

accuracy and robustness of a method. The RMSE is also the preferred error function to be minimized in the parametrization of empirical methods. Another widespread error measure is the mean unsigned error (MUE, also called the mean absolute error, MAE). It is the simple average of the absolute values of the residuals: MUE =

1 N

N

∑ |Ei − Eiref | i=1

(12)

It is a very natural way of expressing the error, but it is less sensitive to problematic points. Although many authors report it as the only error, this approach cannot be recommended, because the RMSE is a better overall error measure. The systematic part of the error can be easily isolated by means of the mean signed error (MSE), the simple average of the errors: MSE =

1 N

N

∑ (Ei − Eiref ) i=1

(13)

It is not a measure of the overall error but a tool for determining whether the results are systematically overestimated or underestimated. In addition to the average measures, it is also useful to list the largest (unsigned) error. A comparison of this value to the average error provides a good clue to the performance of the tested method in the worst-case scenario. All the error measures are absolute; they have the unit of energy and are tied to the magnitudes of the studied variable in the set. They can be used for the discussion of a particular data set or applied to very similar systems. For more universal generalization of the results, the relative error should be used. Most importantly, it allows the comparison of the results between different sets of data with different ranges of energies. It is usually expressed in percent. However, the use of relative errors has an important pitfall of which one should be aware. If the data set contains energies close to zero, the relative error in such systems can be very high, although the absolute one is negligible. As a result, calculating a relative error in the data set as an average of the relative errors in the individual systems makes sense only if the energies in the set are sufficiently large. Despite this problem, this definition of relative error is often used because it is a variable easy to describe and understand. To avoid this problem, it is possible to express the relative error as the ratio of the average error and the average magnitude (i.e., the average of the absolute values) of the variable itself. When applied to the RMSE, a robust relative error that combines the advantages of the RMSE and the transferability of relative errors can be obtained. In statistics, it is called the coefficient of the variation of the RMSE. It is the error measure that we have used for the comparison of errors in different types of interactions in the S66 and X40 data sets.94,187

3.7. Versioning of Data Sets

Some of the data sets evolve even after they have been published, most often by updating the benchmark results with more accurate calculations. In addition to citing the source of the data, it is useful to introduce some simplified labeling to distinguish between the individual versions of the benchmark. Sherrill and co-workers, who recalculated several existing data sets in larger basis sets, label the updated reference data by appending a letter (A, B, ...) to the name of the set.93,193 We propose the use of the standard decimal version numbers in the format X.Y. The major version number X denotes the update of the geometries while the minor version number Y indicates the method used for the calculation of the benchmark energies. This versioning scheme is used in the data set library of our Cuby software.191 In this work, it is applied to the S66 data set, for which three versions exist: version 1.0 corresponds to the original paper,94 version 1.1 was published soon afterward,189 and the most recent version discussed here, 1.2, was later constructed from the available data, improving the accuracy further (see section 4.2 below). It must be noted, however, that updating an existing data set has to be considered carefully and a new version should be introduced only if substantial improvement is achieved.

4. AVAILABLE DATA SETS 4.1. Accurate Data Sets of the Smallest Model Systems

Small model systems, where highly accurate calculations are possible, can serve as tools for estimating the accuracy of the methods used as a benchmark in larger systems. On the other hand, only the molecules above a certain size can be used if we seek realistic modeling of the interactions that occur in larger systems. The development of a data set for this purpose is thus a compromise between these two requirements on the size of the systems. Noncovalent complexes built from molecules with one or two atoms from the second period, saturated with hydrogens, strike this balance rather well. Smaller neutral complexes exhibit too weak interactions. The nature of the interaction of individual atoms, such as in rare-gas dimers, which are often used as a model, does not represent well the complexity of the interactions of molecules. In larger systems, it

3.6. Automated Data Set Processing

Working with data sets requires the management of large amounts of input and output data. For this work to be efficient, it should be automated. This is usually done by means of scripts automating the preparation of the input for the calculations, written individually by each author. However, due to the growing importance of the benchmark data sets, automated calculations of them are becoming supported by some quantum chemistry software packages. The most advanced in this respect 5052


Chemical Reviews

Review

Figure 2. A24 data set of small noncovalent complexes. The color of each box denotes one of the following groups of interactions: blue, hydrogen bonds; green, mixed electrostatic/dispersion; red, dispersion; and purple, π−π dispersion.

suggested that this can lead to a rather large error in this term, especially in polar systems. The values of the CCSDT(Q) corrections in the A24 set, calculated in a series of basis sets, are plotted in Figure 3. It is clear that for all the nonpolar systems,

would be difficult to estimate the CCSD(T)/CBS limit with the required accuracy, and higher-order calculations may not be possible at all. Although there are very accurate calculations in small noncovalent complexes reported in the literature, a systematically built data set of such systems had been lacking. This has led us to the creation of the A24 data set,59 a set of 24 complexes that includes 5 hydrogen bonds, 8 dispersiondominated systems (out of which 3 are models of π−π interactions, covering all combinations of C−C double and triple bonds) and the remaining 11 systems representing various interactions of mixed electrostatic/dispersion nature (see Figure 2). For these, accurate CCSD(T)/CBS energies have been obtained using a three-point extrapolation from aTZ, aQZ, and a5Z data sets on geometries optimized by composite CCSD(T)/CBS gradients. Later, these benchmark interaction energies were improved further by Burns et al., who extended the series of basis sets to a6Z in some systems and extrapolated the CBS limit as a weighted average of counterpoise-corrected and uncorrected results.42 Further corrections to the CCSD(T)/CBS energies were evaluated in order to obtain the best estimate of the interaction energy: (1) correction to CCSDT(Q), recently updated with calculations in larger (mostly aDZ) basis sets;64 (2) the core-correlation contribution extrapolated from aug-cc-pCVTZ and aug-cc-pCVQZ basis sets at the CCSD(T) level; and (3) the contribution of relativistic effects calculated using the Douglas−Kroll−Hess Hamiltonian at the CCSD(T) level. Special attention should be paid to the CCSDT(Q) correction, which is calculated in a rather small basis set. Recent studies105,106 (discussed above in section 2.7) have

Figure 3. CCSDT(Q) correction (calculated as ΔECCSDT(Q) − ΔECCSD(T)) in the A24 data set, calculated in a series of basis sets and extrapolated to the CBS limit from aDZ and aTZ basis sets where applicable.

the accuracy is sufficient and this correction should be applied. Moreover, in the polar systems, the values of the correction are reasonably converged in the aTZ basis. In the polar complexes where aTZ calculations are not available, this correction may even have an incorrect sign. This is why we do not include it in the best estimate of interaction energies in HCN and formaldehyde dimers, where the correction is positive in the 5053


Chemical Reviews

Review

Figure 4. S66 data set. The complexes are divided into groups according to the nature of the interaction. The groups are represented here by colored boxes: blue, hydrogen bonds; orange, red, and purple, dispersion (π−π, aliphatic hydrocarbons, and hydrocarbons−π, respectively); and green, mixed electrostatic/dispersion interactions.

calculations, apart from the error caused by the basis set incompleteness. The A24 data set was used for the extensive benchmarking of composite CCSD(T)/CBS schemes.64 Because of the high importance of these results, they are discussed in detail below in section 5.1.

largest basis set applicable, aDZ. Among the post-CCSD(T) corrections, the higher-order correlation covered by the CCSDT(Q) calculations was found to be of the highest importance [1.14% of the CCSD(T)/CBS interaction energy], followed by the core correlation (0.57%) and relativistic effects (0.14%). These terms partially compensate for each other, and the overall difference between the best estimate and CCSD(T)/CBS is 1.72% [this analysis is based on the updated CCSDT(Q) contribution from ref 59]. This can be considered as an estimate of the intrinsic error of CCSD(T)/CBS

4.2. General Data Sets Covering Interactions of Organic Molecules

The first step in the testing or parametrization of any generally applicable computational method is usually done on 5054


Chemical Reviews

Review

Table 2. Benchmark CCSD(T)/CBS Interaction Energies in the S66 Data Set (in kcal/mol) and Updated Values with the CCSD(T) Correction Calculated in the haTZ Basis Set Are Presented (S66 version 1.2)a interaction 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

Hydrogen Bonds water···water water···methanol water···methylamine water···peptide methanol···methanol methanol···methylamine methanol···peptide methanol···water methylamine···methanol methylamine···methylamine methylamine···peptide methylamine···water peptide···methanol peptide···methylamine peptide···peptide peptide···water uracil···uracil (BP) water···pyridine methanol···pyridine acetic acid···acetic acid acetamide···acetamide acetic acid···uracil acetamide···uracil Dispersion benzene···benzene (π−π) pyridine···pyridine (π−π) uracil···uracil (π−π) benzene···pyridine (π−π) benzene···uracil (π−π) pyridine···uracil (π−π) benzene···ethene uracil···ethene uracil···ethyne pyridine···ethene pentane···pentane neopentane···pentane

ΔE

ΔE

interaction

−4.966 −5.653 −6.984 −8.169 −5.811 −7.619 −8.292 −5.052 −3.085 −4.193 −5.451 −7.347 −6.242 −7.517 −8.681 −5.167 −17.356 −6.931 −7.475 −19.301 −16.434 −19.681 −19.370

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 61 62 63 64 65 66

−2.758 −3.831 −9.771 −3.375 −5.628 −6.734 −1.388 −3.341 −3.707 −1.827 −3.749 −2.591

Dispersion neopentane···neopentane cyclopentane···neopentane cyclopentane···cyclopentane benzene···cyclopentane benzene···neopentane uracil···pentane uracil···cyclopentane uracil···neopentane ethene···pentane ethyne···pentane peptide···pentane Other benzene···benzene (T-shape) pyridine···pyridine (T-shape) benzene···pyridine (T-shape) benzene···ethyne (CH···π) ethyne···ethyne (T-shape) benzene···acetic acid (OH···π) benzene···acetamide (NH···π) benzene···water (OH···π) benzene···methanol (OH···π) benzene···methylamine (NH···π) benzene···peptide (NH···π) pyridine···pyridine (CH···N) ethyne···water (CH···O) ethyne···acetic acid (OH···π) pentane···acetic acid pentane···acetamide benzene···acetic acid peptide···ethene pyridine···ethyne methylamine···pyridine

−1.756 −2.386 −2.977 −3.533 −2.863 −4.811 −4.098 −3.689 −1.988 −1.727 −4.243 −2.843 −3.513 −3.304 −2.860 −1.536 −4.715 −4.387 −3.277 −4.168 −3.204 −5.258 −4.205 −2.899 −4.933 −2.894 −3.516 −3.760 −2.995 −4.068 −3.963

The molecule labeled “peptide” is N-methylacetamide, a peptidebond model.

a

energies in nucleic acid base pairs and amino acid complexes. The JSCH2005 set is rather large but with too narrow a focus, which makes it impractical for general use. Furthermore, due to the size of the systems, which required the use of smaller basis sets, the accuracy was limited. On the other hand, the S22 set was more universal, small enough for routine use, and featured the most accurate calculations available at that time. Because of these advantages, it became a de facto standard for the testing and parametrization of various computational methods. The data set consists of 22 complexes of a size ranging from 6 (water dimer) to 30 atoms (the adenine−thymine base pair), divided into three groups of roughly the same size, each covering one type of interaction: hydrogen bonds, dispersion, and mixed interactions. The S22 set was intended mainly as a model for interactions of nucleic acid bases, which affected the choice of the model systems in these groups. Among the dispersion-dominated interactions, π−π stacking is prevalent and the only representative of another type of dispersion-bound complex is the methane dimer, which is very weak. In the group of hydrogen bonds, the majority of the systems feature multiple very strong H-bonds modeling the Watson−Crick base pairing.

interactions of simple organic molecules, namely, for two reasons: First, these interactions are well-understood, which makes it possible to analyze the results thoroughly. Second, they are a sufficient model for many applications, including the study of biomolecules. Nevertheless, it should be kept in mind that this class of model systems represents only the minimum requirements, and for broader applicability of the method, it has to be validated in other systems as well. The first comprehensive data set of interaction energies was published by Zhao and Truhlar in the testing and development of DFT methods.194,195 It was named NCCE31/05 and consists of several subsets covering dispersion, hydrogen bonds, and charge-transfer complexes. It features the smallest possible models for these interactions, calculated using the W1 theory,196 a composite method including correlation up to CCSD(T), core correlation, and relativistic correction. The history of the wide acceptance of the benchmark data sets of noncovalent interactions started with another set developed at about the same time in our laboratory, the S22 set by Jurečka et al.197 Originally, the S22 data set was designed as a supplement to the JSCH2005 database of the interaction 5055


Chemical Reviews

Review

Figure 5. X40 data set of noncovalent complexes of halogenated molecules. The color of each box corresponds to one of the following groups of interactions: yellow, dispersion; green, induction and dipole−dipole interactions; cyan, π−π stacking; purple, halogen bonds; violet, halogen−π interaction; and orange, hydrogen bonds. The halogens F, Cl, Br, and I are visualized in cyan, green, brown, and purple, respectively.

This makes the set rather unbalanced. To alleviate the issue of hydrogen bonds, four complexes with single H-bonds were added later, forming the S26 set.62 The benchmark interaction energies were calculated by means of the CCSD(T)/CBS scheme using the largest basis set applicable to each complex at that time. The component limiting the accuracy is the CCSD(T) correction, which was calculated only in the cc-pVDZ basis in the case of the largest complexes (a recent estimate of the accuracy of this setup is about 5%). Later, the data set was repeatedly recalculated with higher accuracy. In 2010, Takatani et al. recalculated S22 with the CCSD(T) correction in the aTZ basis,193 bringing the estimated accuracy close to 1%. The updated interaction energies, named S22A, use the correction extrapolated from aDZ to aTZ basis sets. It was later shown that this extrapolation could lead to the deterioration of the results; the use of the aTZ results without extrapolation is now preferred.93 In the same year, S22 was also recalculated by Podeszwa et al. in basis sets augmented with midbond functions.198 Last, Marshall et al. have compiled the best results from the previous works listed here and recalculated some complexes, building the S22B version of the data set.93 The realization of the deficiencies of the S22 set, most importantly its small size and the unbalanced coverage of some classes of the interactions, led us to the development of the S66 data set.94,189 It is intended as a more general replacement of S22, targeting a broader spectrum of interactions. Biomolecular interactions are represented by not only nucleic acid bases (uracil) but also by various interactions of a simple peptidebond model (N-methylacetamide). Special attention was paid to cover diverse interactions in a balanced way. The hydrogenbond group covers various combinations of C and N donors

and acceptors, featuring both single and cyclic double bonds. The dispersion group covers both π−π interactions, interactions in various isomers of aliphatic hydrocarbons, and their mutual combinations. The 66 complexes forming the set are shown in Figure 4. Unlike in S22, all the systems are calculated using the same setup. Although it limits the basis set size to the one applicable to the largest systems in the set, this consistency makes it easier to isolate the errors of the studied methods from basis set size effects. The benchmark interaction energies in the S66 data set were obtained using the composite CCSD(T)/CBS scheme based on MP2/CBS calculation built from HF/aQZ and MP2 correlation extrapolated from aTZ and aQZ. In the original work,94 the CCSD(T) correction was calculated in aDZ (S66 version 1.0). In the A24 data set, this composite scheme has an average error of 2.14%. Later it was recalculated in haTZ, and extrapolation from haDZ and haTZ was used189 (version 1.1; the average error of the scheme in A24 set is 1.45%). Since this scheme was found to be inappropriate for hydrogen bonds, where the aDZ results are not reliable, it is better to drop the extrapolation and use the CCSD(T) correction from the haTZ basis only (the error in the A24 set is 1.06%). This latest version (version 1.2) of the set has not been published per se yet, so we take the opportunity to list the benchmark interaction energies here in Table 2. Especially in the use of S22 for the parametrization of empirical methods, it was soon realized that a data set of equilibrium geometries is not sufficient to obtain a robust fit. This led to the extension of the S22 set with dissociation curves simultaneously by Fusti-Molnar et al.199 and our group,200 where the S22x5 set was developed. The importance of nonequilibrium geometries was taken into consideration 5056


Chemical Reviews

Review

already in the design of the S66 data set, which was published along with a set of dissociation curves for all 66 complexes. This data set of eight points per each system, obtained by scaling the closest intermolecular distance by 0.9, 0.95, 1.0, 1.05, 1.1, 1.25, 1.5 and 2.0, was named S66x8.94 Soon afterward, the coverage of nonequilibrium geometries of the S66 was extended to angular-displaced ones, forming the S66a8 set.189 When combined, the S66 set can be complemented with 1056 nonequilibrium geometries. The interaction energies in the nonequilibrium geometries are calculated at a slightly lower level, using the same MP2/CBS setup but with the CCSD(T) correction calculated in the aDZ basis set.

charge density on the halogen atom), but empirical corrections can be added to fix this issue. We have developed such corrections for PM6202,203 and SCC-DFTB.204 There, the availability of the extended dissociation curves in the X40x10 set is indispensable. For the parametrization of the hydrogen-bonding corrections of SQM methods, we have complemented the S66x8 data with a set of 15 dissociation curves of ionic hydrogen bonds.170 The curves are constructed and calculated using exactly the same setup as used in the S66x8 set. The systems in the set have been built specifically as models of the side chains of charged amino acids.

4.3. S66-Compatible Data Sets

4.4. Data Sets of Specific Interaction Types

To extend the coverage to other elements while keeping the compatibility with the S66 and S66x8 data sets, the following sets were developed along the same guidelines and using the same methods of geometry optimization and benchmark calculations. A small data set dedicated to the compounds of sulfur was developed by Mintz and Parks.201 The data set is biased toward hydrogen bonds, containing 10 complexes with H-bonds of sulfane, methanethiol, and thioacetone with each other or with water. This is complemented by only three complexes of these molecules with dispersion-dominated interaction. Eight-point dissociation curves for all of these complexes were generated analogously to S66x8. To cover the interactions of halogen-containing molecules, we have developed the X40 data set.187 The 40 noncovalent complexes it consists of are shown in Figure 5. The data set focuses mainly on halogen bonding, featuring 14 complexes covering the interactions of Cl, Br, and I with N, O, and S. This is an important feature, because the halogen bond is a prototypical example of the so-called σ-hole bonding, a class of interactions not covered by the S66 data set. The rest of the set covers dispersion, dipole−dipole, and induction interactions; the effects of fluorination on benzene stacking; halogen−π interactions; and hydrogen bonds. The methodology used is exactly the same as in the S66 data set (using the CCSD(T) correction in the haTZ basis, which is equivalent to S66 version 1.2). The core set of equilibrium geometries is complemented by dissociation curves forming the X40x10 set. Because of the importance of the close contacts in halogen bonding, these curves contain two extra points at 0.8 and 0.85 of the equilibrium distance in addition to the protocol used for S66x8. This data set has two important uses. The first is in method testing, where it can be used as a cross-check of transferability from organic molecules to other elements. Especially for methods including a small number of empirical parameters that are usually parametrized on organic systems (e.g., S66x8), comparing the relative error between those and the halogenated systems is a good measure of the robustness of the method. We have found that correlated methods, including the parametrized ones such as spin-component scaled approaches and MP2.5, conserve their performance rather well (with MP2.5 and SCSMI-MP2 being the most accurate at the given complexity level). The error of DFT-D grows slightly, and SQM methods have serious difficulties describing all the interactions well.187 The second important use is the parametrization of more empirical methods toward the description of halogen bonding. Semiempirical QM methods cannot describe this phenomenon well (due to their inability to reproduce the anisotropy of the

There exist multiple smaller data sets targeting interactions of specific nature. Many of these interactions are covered by universal data sets; in this case, these extra data sets might be used as additional validation data (after checking for the overlap of the sets). Some classes of interactions, however, are not included in the universal sets. From the sets listed here, this concerns pnictogen and chalcogen bonds (σ-hole bonding of group 15 and 16 elements) and charge-transfer complexes. All of these data sets are briefly introduced here and listed in Table 3. Table 3. Data Sets of Benchmark Interaction Energies Dedicated to Specific Interaction Types featured interactions

set name

set size

refsa

dispersion dispersion in hydrocarbons dispersion in heavy-atom hydrides hydrogen bonds halogen, chalcogen and pnictogen bonds charge transfer

NBC10  HEAVY28 HBC6 

10 curves 12 28 6 curves 30

93, 205 206 158 93, 207 208



11

209, 210

a

The additional references point to updated versions of the data sets.

Dispersion interactions are covered by the following data sets: The NBC10 data set of Sherrill et al.205 contains 10 dissociation curves, focusing mainly on complexes of benzene. The coverage of dispersion in both saturated and unsaturated hydrocarbons is more balanced in the set of Granatier et al.206 For the validation of the methods applicable to complexes of heavier atoms, Grimme et al. have built the HEAVY28 data set, featuring hydrides of Pb, N, Sb, Bi, O, S, Te, Cl, Br, and I.158 Hydrogen bonds are covered by the HBC6 set of Thanthiriwatte et al., featuring dissociation curves of six complexes.207 A data set of σ-hole interactions, namely halogen, chalcogen, and pnictogen bonds in very small model systems (both neutral and charged), has recently been built by Bauzá et al.208 Another class of interactions often studied separately are complexes with strong charge-transfer character. Small model systems are covered by the set of Karthikeyan et al.,209 recently recalculated with higher accuracy by Ř ezáč and de la Lande.210 4.5. Data Sets Targeting Specific Systems

In this section, we list data sets designed as models of particular systems (for a summary, see Table 4). Among them, special attention is paid to data sets relevant to biochemistry. Large biomolecules are often studied by approximate computational methods, including MM force fields. The benchmark data covering interactions of the building blocks comprising these 5057


Chemical Reviews

Review

Another option to obtain benchmark data for large systems is to derive them from experiment, starting from free energies measured in solution. From these, computed solvation and thermodynamic contributions are removed to acquire the gasphase interaction energies. These calculations have to be carried out using approximate and empirical methods, which limits the accuracy of the final interaction energies. Grimme applied this approach to a set of 12 supramolecular complexes, building the S12L data set.218,219 Some of these systems were later recalculated using quantum Monte Carlo.139 This set has recently been extended to 30 systems and named S30L.220

Table 4. Data Sets of Benchmark Interaction Energies Focusing on Specific Systems featured interactions

set name

set size

DNA and RNA base pairs, amino acid dimers amino acid side chains fragments of the protein−ligand complex nucleic acid−protein contacts water clusters water clusters large noncovalent complexes large supramolecular complexes

JSCH2005

143

197

SCAI HSG

24 11

211 93, 213

 WATER27  L7 S12L

272 27 61 7 12

S30L

30

212 214 216 217 139, 218, 219 220

large supramolecular complexes

refs

4.6. Data Sets of Related Properties

Here we will briefly review the available data sets of other properties relevant to noncovalent interactions other than interaction energies. The geometry of smaller noncovalent complexes can now be optimized at a very high level and used as a benchmark. The geometric properties are, however, less sensitive to the quality of the methods tested than interaction energies, and the interpretation of the results is less straightforward. Efficient implementations of analytical gradients of CCSD(T)221 are now available, which makes it possible to obtain benchmark geometries at this level. Apart from applications to individual systems, several data sets of such geometries are already available. The A24 data set59 features geometries optimized at the CCSD(T)/CBS level, employing a composite scheme based on MP2 extrapolated from aTZ and aQZ basis sets and CCSD(T) correction in the aDZ basis. These systems, along with larger complexes for which CCSD(T)-quality geometries were interpolated from single-point calculations, have recently been used for the benchmarking of geometries optimized with multiple cheaper methods.185 Using geometries as a benchmark is more difficult because it is necessary to run the geometry optimization employing the method tested. While a reasonable estimate of errors in the geometry can also be obtained by analyzing single-point energies in nonequilibrium geometries, more direct information is provided by a detailed comparison of optimized geometries with the benchmark ones. It is also possible to go beyond two-body interactions and calculate many-body interaction energies. Among them, threebody interaction energies, which are energetically the most important, can be calculated with benchmark accuracy in medium-sized molecules. The importance of benchmarking computational methods on the many-body energies is relevant not only to this property itself but also, more generally, to how the methods work in larger systems. The accuracy of the description of three-body energies is a good indicator of how the accuracy measured on pairwise interaction energies in smaller model systems translates to large systems composed of a large number of such units, either covalently or noncovalently bound. We have published the first comprehensive benchmark data set of three-body energies.152 It contains 69 trimers, obtained by taking three different conformers from 23 molecular crystals. It covers both stabilizing and destabilizing three-body interactions of both polarization and dispersion nature. The three-body energies are calculated at the CCSD(T)/CBS level, using a composite scheme with CCSD(T) correction in the aDZ basis set. Another quantity useful for benchmarking is conformational energies of flexible molecules, the structure of which is determined by intramolecular noncovalent interactions. Here,

molecules are a valuable tool for the testing and parametrization of these methods, and multiple such data sets are available. The JSCH2005 data set of Jurečka et al. was the first data set of this type.197 It covers nucleic acid base pairs and amino acid dimers. Interactions of amino acid side chains in geometries representing the most common motifs in proteins form the SCAI data set of Berka et al.211 Recently, a large data set of nucleobase−amino acid complexes mapping representative nucleic acid−protein contacts has been developed by Hostaš et al.212 As a benchmark model for protein−ligand interactions, Faver et al. have built a data set of complexes obtained by the fragmentation of both the ligand and the surrounding protein chain in an experimental geometry (PDB ID 1HSG).213 For the modeling of the interactions in water, the WATER27 data set of water clusters of Bryantsev et al.214 can be used. The larger clusters from the set have recently been recalculated with higher accuracy.215 Another data set of water clusters has been published by Temelso et al.,216 who used it for the study of the basis set dependence of the CCSD(T) correction specifically in water. Another criterion for sorting noncovalent complexes is their size. Small- and medium-sized systems (up to several tens of atoms) are covered well by all the data sets discussed so far, but accurate interaction energies in larger systems are scarce. These large complexes (with a size of ca. 50 atoms and larger) are highly important for the validation of approximate methods developed using the benchmark data sets of smaller systems. Nevertheless, it is very difficult to obtain accurate and reliable interaction energies in systems of this size. It is possible to extend the standard CCSD(T)/CBS methodology to systems of this size, but it is necessary to both reduce the basis set size and use very large computational resources. The limitations of the basis set size can be partially compensated for by fine-tuning the basis set and extrapolation to the CBS limit. Although it is not possible to achieve accuracy as high as in smaller systems, this approach has the advantage of making it possible to compare various computational methods calculated in the same basis sets, eliminating a large part of the basis set effects. We used this approach in the development of the L7 data set,217 a set of 7 mostly dispersion-bound complexes with 48−112 atoms. The benchmark energies were obtained with either the CCSD(T) or very similar QCISD(T) (quadratic configuration interaction up to perturbative triples) method and extrapolated to CBS using a composite scheme. 5058


Chemical Reviews

Review

Table 5. Errors (the relative mean unsigned error, MUE%, and the root mean square error, RMSE, in kcal/mol) of Composite CCSD(T) Calculations in the A24 Data Seta basis sets HF QZ-F12 QZ-F12 aQZ aQZ QZ-F12 QZ-F12 QZ-F12 QZ-F12 QZ-F12 aQZ aQZ QZ-F12 QZ-F12 aQZ QZ-F12 aQZ QZ-F12 aQZ aQZ aQZ QZ-F12 aQZ aQZ aQZ aQZ QZ-F12 QZ-F12 QZ-F12 aQZ QZ-F12 QZ-F12 aQZ aQZ QZ-F12 aQZ aQZ aQZ QZ-F12 QZ-F12 aQZ TZ-F12 TZ-F12

MP2

basis sets

errors ΔCCSD(T)

Quadruple-ζ QZ-F12 haQZ QZ-F12 aQZ aTZ→aQZ QZ aTZ→aQZ haQZ TZ-F12→QZ-F12 QZ QZ-F12 haTZ QZ-F12 QZ QZ-F12 aTZ TZ-F12→QZ-F12 haQZ aTZ→aQZ haTZ aTZ→aQZ aQZ TZ-F12→QZ-F12 aQZ TZ-F12→QZ-F12 haTZ aTZ→aQZ aTZ TZ-F12→QZ-F12 aTZ aQZ aTZ QZ-F12 aDZ aTZ→aQZ aDZ aQZ aQZ aQZ 6-31G**′b TZ-F12→QZ-F12 aDZ aQZ haQZ aQZ haTZ aQZ aDZ aTZ→aQZ TZ TZ-F12→QZ-F12 TZ QZ-F12 TZ QZ-F12 6-31G**b aTZ→aQZ haDZ QZ-F12 haDZ TZ-F12→QZ-F12 haDZ aQZ QZ aTZ→aQZ 6-31G**b TZ-F12→QZ-F12 6-31G**b aQZ haDZ aQZ TZ aTZ→aQZ DZ TZ-F12→QZ-F12 DZ QZ-F12 DZ aQZ DZ Triple-ζ DZ-F12→TZ-F12 aTZ TZ-F12 aTZ

MUE%

RMSE

HF

0.53 0.57 0.84 0.91 0.94 0.94 0.95 1.00 1.02 1.06 1.07 1.18 1.21 1.53 1.62 1.72 1.95 2.14 2.16 2.20 2.26 2.31 2.47 2.63 2.71 2.80 2.92 3.13 3.20 3.29 3.31 3.53 3.54 3.60 4.45 4.89 5.88 6.01 6.02 7.08

0.012 0.015 0.018 0.016 0.016 0.014 0.017 0.018 0.020 0.018 0.018 0.022 0.021 0.023 0.026 0.048 0.041 0.044 0.047 0.074 0.044 0.050 0.059 0.085 0.054 0.052 0.055 0.085 0.061 0.061 0.060 0.072 0.087 0.091 0.110 0.108 0.125 0.125 0.125 0.168

TZ-F12 aTZ TZ-F12 aTZ TZ-F12 TZ-F12 aTZ TZ-F12 TZ-F12 aTZ TZ-F12 TZ-F12 aTZ aTZ TZ-F12 TZ-F12 aTZ TZ-F12 aTZ TZ-F12 aTZ aTZ aTZ aTZ aTZ aTZ

0.62 0.70

0.014 0.012

DZ-F12 DZ-F12 DZ-F12 DZ-F12 DZ-F12 aDZ aDZ aDZ aDZ aDZ

MP2

errors ΔCCSD(T)

Triple-ζ DZ-F12→TZ-F12 haTZ aDZ→aTZ haTZ TZ-F12 haTZ aDZ→aTZ aTZ DZ-F12→TZ-F12 aDZ TZ-F12 aDZ aDZ→aTZ aDZ TZ-F12 6-31G**b DZ-F12→TZ-F12 6-31G**b aDZ→aTZ 6-31G**b DZ-F12→TZ-F12 TZ DZ-F12→TZ-F12 haDZ aDZ→aTZ haDZ aDZ→aTZ TZ TZ-F12 haDZ TZ-F12 TZ aTZ 6-31G**b DZ-F12→TZ-F12 DZ aDZ→aTZ DZ TZ-F12 DZ aTZ aTZ aTZ aDZ aTZ haTZ aTZ haDZ aTZ TZ aTZ DZ Double-ζ DZ-F12 6-31G**b DZ-F12 aDZ DZ-F12 haDZ DZ-F12 TZ DZ-F12 DZ aDZ 6-31G**b aDZ aDZ aDZ haDZ aDZ TZ aDZ DZ

MUE%

RMSE

0.77 1.24 1.32 1.34 1.58 1.85 1.91 2.22 2.54 2.86 2.89 3.05 3.16 3.37 3.60 3.74 4.85 5.73 5.93 6.24 6.39 6.86 7.14 8.77 9.56 9.72

0.012 0.060 0.020 0.052 0.037 0.045 0.085 0.074 0.082 0.089 0.056 0.060 0.105 0.106 0.070 0.068 0.139 0.123 0.166 0.131 0.145 0.178 0.156 0.203 0.205 0.254

2.35 3.29 5.14 5.83 7.54 20.20 22.21 24.12 24.91 25.07

0.067 0.074 0.100 0.101 0.155 0.399 0.449 0.475 0.476 0.517

a

The calculations in the -F12 basis sets use explicit correlation. The calculations are grouped by the basis set used in the MP2 term, ordered by increasing MUE%. The arrows denote extrapolation to the complete basis set limit. CCSD(T) calculations in a single basis are highlighted in bold. b6-31G**′ stands for the 6-31G**(0.25,0.15) basis set.

4.7. Resources and Aggregated Data Sets

the resulting energies are defined by the interplay of these interactions and the covalent structure, mainly torsional potentials of rotatable bonds. These conformational energies are very sensitive to the computational method used. There are several sets of benchmark-quality CCSD(T)/CBS conformational energies for biomolecules from amino acids222−225 to DNA backbone fragments226 and small peptides.227 The latter systems contain a very large number of the local minima, and many conformers have energy very close to the global minimum. These data are directly relevant to the gas-phase spectroscopy of these molecules, where only limited information on their structure is available, but they can also serve as a universal benchmark.

The benchmark data and the geometries of the systems are usually published as Supporting Information accompanying the papers introducing them. This is, however, not a very convenient way of working with the data. To simplify the use of the benchmark data sets, several groups have created repositories of these data, which allow easier access to them and provide additional functions. Some of these collections are also intended as more general benchmark suites for method testing. To make ours as well as contributed benchmark data more accessible, we have built the Benchmark Energy and Geometry Database (BEGDB), a Web-based repository of data sets.228,229 5059


Chemical Reviews

Review

extrapolation is possible. The composite schemes based on MP2-F12/DZ-F12 calculations can reach surprisingly low errors below 3%, whereas those without explicit correlation have errors an order of magnitude larger. Second, the accuracy at a given level of MP2 calculations is determined by the basis set used for the CCSD(T) correction. Here, the importance of diffuse functions is clearly visible. With the exception of the QZ basis set, which is large enough to compete with the augmented basis sets, the basis sets without those functions perform substantially worse. Even the TZ basis set is outperformed by a basis as small as 6-31G**(0.25,0.15). The heavy-augmented basis sets, however, yield results about the same as (and sometimes even better than) their fully augmented counterparts. The use of these basis sets thus makes it possible to save some computational resources without compromising the accuracy. Third, it is interesting to compare the CCSD(T) calculations in a single basis set (highlighted in Table 5) with the composite schemes. In the aQZ basis, the error is rather small (2.2%), which makes this method a viable alternative in applications where composite calculations are not practical. Passing to aTZ, the error increases to 6.4%, which is substantially more than in many cheaper composite calculations. In the aDZ basis, the error is as large as 22.2%. So far, we have discussed the errors against the CBS limit obtained from frozen-core CCSD(T) calculations. It is also possible to evaluate the error against the best estimate that contains the most important effects neglected at this level. This analysis is discussed in detail in ref 64, where updated values of the best estimate are available. The only difference is observed in the most accurate schemes, where this small change of reference plays some role. In some cases, a favorable error compensation between the error of the CCSD(T)/CBS estimate and the missing corrections occurs and the overall results are thus very similar. Especially the scheme based on MP2 extrapolated from aTZ and aQZ basis sets and CCSD(T) correction in the haTZ basis (the one used as a benchmark in the S66 and X40 data sets) has similar errors of slightly over 1% against both references. As a summary, we would like to highlight two setups that provide an exceptional accuracy-to-cost ratio. For accurate calculations of medium-sized systems, we can recommend using either MP2-F12/QZ-F12 or extrapolation from MP2/aTZ and MP2/aQZ with CCSD(T) correction in the haTZ basis. With this scheme, accuracy close to 1% can be reached. If necessary, the aDZ basis can be used for calculating the correction, and the errors grow to about 2%. For calculations of larger systems, MP2-F12 calculation in the DZ-F12 basis with CCSD(T) correction in 6-31G**(0.25,0.15) uses basis sets as small as possible while it retains an excellent accuracy of 2.35% due to (presumably rather large) error compensation. The accuracy of this scheme was further validated in the S66 data set, and the results are equally good.64

It makes it possible to browse the data sets with details on each calculation, visualize molecular geometries, plot the errors of the methods in the data sets, and download tables of the results and geometries. The database can also be searched across the data sets. Other authors have published collections of previously published data sets. Goerigk and Grimme have published the GMTKN24 database,230 a collection of 24 data sets of various properties including also noncovalent interactions. Later, it was extended by six more data sets to GMTKN30.231 The geometries and reference data are available at the author’s Web site.232 The data sets of Zhao and Truhlar194,195 are also available online.233

5. IMPORTANT RESULTS 5.1. The Accuracy of CCSD(T)/CBS

When using CCSD(T)/CBS calculations as a benchmark, one should be aware of their accuracy. The errors in interaction energies caused by the neglected contributions (described above in section 2.7) are about 2% in organic molecules (see section 4.1). The second source of error, more important in practice, is the incompleteness of the basis set. It can be partially corrected by using extrapolation to the CBS limit and composite CCSD(T) schemes. This provides many possible combinations of techniques and basis sets. We have recently studied the error of a large number of composite schemes applicable to larger molecules in the A24 data set, which offers an accurate estimate of the CBS limit as a reference.64 To select meaningful composite schemes (based on eqs 5 and 6), we explore the combinations of basis sets up to aQZ satisfying the following rules: (1) HF energy is calculated in the same basis set as MP2; if the MP2 term is extrapolated, the HF energy in the larger basis set is used without extrapolation. (2) MP2 correlation is either calculated in a single basis set or extrapolated to the CBS limit and calculated with or without explicit correlation. (3) The ΔCCSD(T) term is calculated in a single basis set not larger than the one used for MP2. The resulting errors are summarized in Table 5. Calculations based on MP2 calculations in quadruple-, triple-, and double-ζ basis sets are listed separately. The errors are measured against the accurate CCSD(T)/CBS from the original paper,59 which were extrapolated from CCSD(T) calculations in aTZ, aQZ, and a5Z basis sets using an adaptive three-point fit. Besides the obvious trend of improving the results with basis set size, there are several interesting features to be discussed. First, we will focus on the comparison of schemes based on canonical and explicitly correlated MP2 calculations. Of course, the use of canonical MP2 without extrapolation leads to results substantially worse than the other options. On the other hand, the use of extrapolation in the explicitly correlated calculations (parametrized specifically for this method86) does not bring any improvement. Only the two remaining options, the extrapolation from canonical MP2 calculations and explicitly correlated MP2-F12 calculations without extrapolation, are useful in practice. The results show that extrapolation from a(X-1)Z and aXZ basis sets leads to errors comparable to explicitly correlated calculation in the XZ-F12 basis. Since the explicitly correlated calculations are more demanding, these two approaches are comparably expensive.234 The most important advantage of the explicitly correlated calculations thus lies in the calculations in double-ζ basis sets, where no

5.2. The Validation and Parameterization of Wave Function Methods

The two most important applications of the data sets of benchmark data are the evaluation of the accuracy of existing computational methods and the parametrization of newly developed ones. Here, we will demonstrate the use of our best data sets in the testing of mostly high level correlated wave function methods and the use of benchmark data sets in their parametrization. 5060


Chemical Reviews

Review

Table 6. Errors (RMSE in kcal/mol and the relative RMSEa in percent) in S66 and X40 Data Sets S66

X40

method

RMSE

RMSE, %

SCS-MI-CCSD/CBS SCS-CCSD/CBS CCSD/CBS MP2.5/CBS MP3/CBS MP2C SCS-MI-MP2/CBS SCS-MP2/CBS DW-MP2/CBS MP2/CBS MP2/TZ MP2/aDZ BLYP-D3(BJ)/def2-QZVP PM6-D3H4X PM6

0.12 0.31 0.70 0.18 0.60 0.14 0.38 0.90 0.40 0.72 0.76 0.85 0.24 0.67 3.05

2.2 5.7 12.8 3.3 10.9 2.6 6.9 16.4 7.3 13.1 13.9 15.5 4.4 12.2 55.6

RMSE

RMSE, %

0.06 0.19 0.48 0.18 0.45

1.7 4.9 12.7 4.7 12.0

0.39 0.40 0.32 0.70 0.58 0.55 0.39 2.31 (0.62)b 2.80

10.4 10.7 8.5 18.6 15.4 14.6 10.4 61.4 (16.5)b 74.3

a

Calculated as the RMSE divided by the average magnitude of the interaction energy in the set. bThe PM6-D3H4X results in parentheses were evaluated in the X40 set without the hydrogen bonds of halogens, for which no correction is available.

case, the (T) term scaling with N7 is the most time-consuming step of the calculation. At the MP3 level, MP3 itself considerably underestimates the strength of the interactions. When the MP3 term is scaled down in the MP2.5 method, the results improve significantly (the relative RMSE is 3.3% in S66 and 4.7% in X40). Although there is a slight increase of the relative error in the X40 set, MP2.5 is still the most accurate method at this or lower levels of computational complexity. Here, the limiting step is the calculation of the MP3 correlation, which scales with N6. The scaling is the same as in CCSD, but because MP3 is not iterative, it is about 15−20 times faster. At the MP2 level, the only method to have accuracy comparable to higher-level approaches is MP2C, where the deficiencies of MP2 are corrected by physically sound means of replacing the ill-behaved component of the dispersion with a more rigorous calculation. It yields a very low error of 3.3% in the S66 set. Unfortunately, we have not been able to test the transferability of this method to the X40 set because the current implementation of the method is not compatible with the use of the pseudopotentials, which are required for the description of the heavier halogens. Spin-component scaling, even in SCSMI-MP2 or DW-MP2 versions, developed specifically for noncovalent interactions, does not improve the results of MP2 much with an error of 7−10%. A comparison of SCS-MI-MP2 in S66 and X40 sets shows that the parametrization carried out on organic molecules is not well-transferable to other elements and types of interaction. The worst results are obtained with plain MP2 extrapolated to the CBS, where the dispersion interaction is the most exaggerated. The error amounts to 13% in S66 and 19% in X40. This overestimation can be partially compensated for by using a smaller basis set; according to our survey of the options, the use of the TZ basis can be recommended.62 Even though the overall error does not change much, it is better balanced; the description of dispersion and especially π−π interactions is improved while the error in hydrogen bonds (described well at the CBS limit) increases. For more details on the decomposition of the results in the groups of complexes with different natures of interactions, see the original papers on S6694 and X40187 data sets.

Using both the S66 set, which covers organic molecules, and X40, which adds halogenated compounds, makes it possible to test the transferability of the methods to different systems. For the first time, the benchmark data in such a comparison were obtained using exactly the same setup, the haTZ basis set for the calculation of CCSD(T) correction (i.e., using the updated S66 interaction energies available in this paper in Table 2). Furthermore, some computational methods were added to the tests published previously in the paper on the X40 seta set.187 To allow a fair comparison, the relative RMSE is used (see section 3.5 for definition). The tested methods are described above in section 2.10. The results are summarized in Table 6, sorted by decreasing computational requirements. In addition to correlated wave function methods, one representative of the DFT-D methodology and one corrected semiempirical method are added for comparison. Both were selected as the best performers in their class. Unless otherwise indicated, the results of all the methods tested were obtained using a composite scheme based on the same HF/aQZ and extrapolated MP2/(aTZ,aQZ) terms as utilized in the benchmark. The same extrapolation was used in MP2/CBS calculations and all their variants. There is, however, a difference in the higher-order correction in CCSD- and MP3based methods, which is calculated in the aDZ basis (only these data were available in the literature for all the methods). When these results are compared to the benchmark using CCSD(T) correction calculated in the haTZ basis set, some error caused by this difference is added to the intrinsic error of the method. This error is, however, very small (about 1%), so that its effect on the overall results can be neglected. Starting from the most expensive methods, plain CCSD substantially underestimates the interaction energies, with the error in both sets being 13%. This can be improved by spin-component scaling, and its variant optimized for noncovalent interactions, SCS-MI-CCSD, achieves excellent results with errors around 2%. It should be highlighted that this method has been parametrized on the S22 data set154 and the improvement observed upon transition to the halogenated complexes in the X40 set is an important proof of its reliability and transferability. It should be noted that while CCSD calculations are still very demanding, they are substantially more economic than CCSD(T). In the latter 5061


Chemical Reviews

Review

many-body interactions. These results should be applicable also to many-body interactions of individual blocks of large molecules. Using CCSD(T)/CBS three-body energies as a reference, we studied the accuracy of a wide range of methods from the accurate correlated ones to DFT.152 In general, the high-level methods describe both two- and three-body interactions similarly. This is a valuable proof of the robustness of the parametrized QM methods, such as SCS-MI-CCSD or MP2.5, which both worked very well here. On the other hand, the difficulties in describing the three-body energies with DFTD (including an empirical three-body dispersion term in the correction) indicate that the error of DFT may grow when passing to larger systems, where the three-body interactions play an increasingly important role. More direct conclusions can be drawn by studying interactions in the dimers of molecules of a systematically increasing size. We are currently working on such a project focused on the performance of local and approximate CCSD(T) methods. Since we do not have other data allowing a comparison of these methods, we will include some of the results here. Namely, we will discuss two of the studied systems, a dimer of polyacetylene chains in CH−π and π−π arrangements. Chains of two, four, six, and eight carbon atoms have been used. Including hydrogens, the largest dimer has 36 atoms, which makes it possible to use full CCSD(T)/aDZ calculations as a benchmark. In addition to the methods tested in smaller systems, we have also investigated the trends in the errors caused by approximations at the CCSD(T) level, namely, the use of frozen natural orbitals124 in FNO−CCSD(T) and the local methods utilizing domain-localized pair natural orbitals, DLPNO−CCSD(T),135 and local natural orbitals, LNO−CCSD(T).134 The CCSD(T) benchmark was calculated in Turbomole.235 In all the approximate methods, the most accurate setup recommended by the authors for noncovalent interactions was used: FNO−CCSD(T) in Psi4190 (occupancy cutoff 10−6), DLPNO−CCSD(T) in Orca236 (using the TightPNO set of thresholds136), and LNO−CCSD(T) in the MRCC package237 (using the parameters recommended in ref 134 and the MP2 correction described therein). The errors discussed here are unsigned relative errors of the contribution of the correlation to the interaction energy; this is done for consistency between the systems, where the HF interaction energy is negative in the CH−π geometry and positive in the π−π one. The results are summarized in Figure 6.

As a representative DFT-based method, we have selected the BLYP/def2-QZVP setup with the D3 correction using the Becke−Johnson damping.159 The performance of DFT-D depends quite strongly on the functional used, and this one was selected as the most accurate among those that can be calculated efficiently.160 The use of a rather large basis set is needed to obtain this level of accuracy as well. In organic systems, the results are excellent and DFT-D3 clearly outperforms all MP2-based approaches, which are much more expensive. In the X40 set, the results are slightly worse. One hypothesis was that this can be caused by charge-transfer effects incorrectly described by the GGA functional, but no improvement has been observed by passing to an analogous hybrid one, B3LYP, which yields the same error. It may be caused by other errors of DFT or by the inaccuracy of the dispersion correction, which has been parametrized on the S66x8 data set and thus performs extremely well in S66. The PM6-D3H4X method serves as an example of an advanced SQM method with corrections for dispersion, hydrogen bonds, and halogen bonds. However, hydrogen bonds involving halogens are not corrected. This leads to a large error in the X40 set, which reduces to a value comparable to S66 only when these systems are removed from the set (the values in parentheses in Table 6). The errors of uncorrected PM6 are much higher, above 50% in both sets. These results, in combination with studies confirming the quality of the CCSD(T)/CBS benchmark used, allow us to quantify confidently the accuracy of various approximate methods and select those with the best properties. Among the most complex methods based on CCSD, SCS-MI-CCSD yields the best results, which closely reproduce the CCSD(T) benchmark. It works well in different types of systems and describes even three-body interactions with equally high accuracy.152 One step down in the computational complexity are methods based on MP3, among which MP2.5 works best. Again, it performs well in a wide range of systems,151 including three-body interactions.152 With the exception of MP2C, methods based on MP2 cannot be recommended for most applications, as they are less accurate than the substantially less expensive dispersion-corrected DFT. 5.3. Transition to Large Systems

The performance of approximate computational methods in small- and medium-sized model systems is now welldocumented, as demonstrated in the previous two sections. The situation is less clear in larger molecular complexes, where accurate benchmark calculations are difficult or impossible. It is possible to extend the same benchmarking methodology to larger systems, such as was done in the L7 data set,217 although the accuracy of the benchmark is somewhat lower because of the use of smaller basis sets. Another option is to resort to experimental data, such as in the S12L and S30L data sets.218,220 These data sets are described above in section 4.5. Here, it is possible to study larger systems, but the accuracy of the interaction energies is even lower. In both cases, the results confirm the observations made in smaller systems, but the limited accuracy of the benchmark does not allow for such accurate quantification of the errors. An alternative approach to this problem is to study the transition from small to larger systems in the range where accurate benchmarking is possible and to extrapolate this information to larger systems. For condensed systems consisting of smaller molecules, this means investigating

Figure 6. Relative error of the studied methods as a function of the size of the system. The points represent results for the two geometric arrangements, and the line represents the average error. 5062


Chemical Reviews

Review

5.4. The Parametrization of Semiempirical QM Methods

The trends in the errors of MP2 and CCSD seem to converge to a finite value in systems larger than those studied. The rate convergence of the MP2 error is especially slow. MP2.5 behaves much better and the errors in all but the smallest systems are very similar, with values of about 4.4%. This indicates that the error cancellation between MP2 and MP3, on which the MP2.5 method is built, gives the method very favorable properties and makes it a good choice for studying large systems. Among the approximate CCSD(T) methods, FNO−CCSD(T) introduces only a negligible error (as the approximation is very gentle and the scaling remains the same). Importantly, this error of about 0.1% is independent of the system size. The local DLPNO−CCSD(T) method has a larger error, as expected. In the three larger systems, it is practically constant, which implies that this method will behave well in even larger systems. This is an important finding, as this method, with its near-linear scaling, is applicable to much larger systems, where no benchmarking is possible. The observed error of about 5% is slightly higher than the error evaluated earlier in the S66 data set, where it amounts to 1.8% (0.1 kcal/mol),136 possibly because of the difficulties specific to conjugated π-electron systems. Only the errors of LNO−CCSD(T) seem to grow with no signs of convergence, which is an unfavorable behavior, although they are comparable to the errors of DLPNO− CCSD(T) in the studied range of the system sizes. This magnitude of errors is consistent with the tests performed by the authors of the methods on a subset of the S22 data set, where an average error of 0.12 kcal/mol was reported.134 It should be noted, however, that these results are only preliminary and the study will be extended to more systems in order to provide more exhaustive coverage of different classes of interactions. The accuracy of these methods should be also discussed in the context of the time needed to perform the computation. Here, we provide the timings for the calculation of the interaction energy in the largest CH−π dimer running in parallel on an eight-core Intel Xeon server (this comparison is, however, only approximate, as the calculations were run on similar but not identical nodes of a heterogeneous cluster and are performed with different software). The benchmark CCSD(T) calculation took 187 h. The FNO approximation decreases this time to 142 h, 75% of the full CCSD(T) calculation. Here, the savings come from the calculations of the monomers with counterpoise correction, where the virtual space can be substantially reduced; the calculation of the dimer is almost identical to a full calculation. The local methods are much faster with DLPNO−CCSD(T) at 9 h (5%) and LNO− CCSD(T) at 5 h (3%). These times are comparable to the MP2.5 calculation, which took 6 h. The reduced scaling of the local method must be also considered; the methods should be applicable to even larger systems where the full or FNO− CCSD(T) calculations would not be possible at all. The results presented here illustrate the difficulty of obtaining reliable interaction energies in large noncovalent systems. On the other hand, substantial savings in the computational time can be expected. Some of the recent developments are very promising, although there is a great deal of work to be done to characterize the behavior and quantify the accuracy of these methods in large systems.

The approximations involved in semiempirical QM methods severely limit their ability to describe noncovalent interactions. Additionally, the performance of these methods depends not only on their formulation but also on their parametrization, in which noncovalent interactions were, until very recently, often neglected or underrepresented. There are four main issues that stem from both these limitations: (1) As in all one-electron methods, London dispersion is missing, although some shortrange part of it may be covered by the empirical treatment of the correlation. (2) The use of the subminimal basis set, the limited polarizability of atoms, and a simplified electrostatic model lead to the underestimation of H-bonds. (3) Intermolecular Pauli repulsion at short-range is often incorrect due to the insufficient parametrization of the core−core potentials. (4) Less common interactions, such as σ-hole bonds, are described poorly both because of the limitations of the theory and their neglect in parametrization. This topic is covered by a separate review in this issue; here we will discuss only the role of benchmark data sets in the development of corrections to these issues. The missing dispersion is a well-defined issue separated from the rest of the method and it is easy to correct by stand-alone empirical correction analogous to DFT-D. In fact, such a correction was first applied to the semiempirical self-consistent density-functional tight-binding (SCC-DFTB) method by Elstner et al.164 in 2001 and only later applied to the DFT. Subsequently, it was applied also to a wide range of MNDOtype methods.165−169,238 We have recently adapted the advanced D3 correction of Grimme to many SQM methods.170 The underestimation of hydrogen bonds is more difficult to solve, as it has origins in multiple approximations that make the SQM methods so efficient. In the DFTB method, the description of H-bonds has been improved, although not completely fixed, by modifying the short-range electrostatic interaction of all hydrogen atoms.239 In MNDO-type methods, we have introduced and developed the idea of empirical corrections for hydrogen bonding.168,169 These corrections involve the application of a fully empirical, local force field to all possible H-bonds in the system, calculated independently and added to unmodified SQM calculation. The last version of the corrections, named H4 (to be used with the D3 dispersion as the D3H4 method), makes it possible to achieve an accuracy below 1 kcal/mol in a wide range of organic systems and is simple and robust.170 It has been parametrized for a wide range of SQM methods; see Table 7. Another empirical correction was necessary for halogen bonds. Here, the SQM methods we used (SCC-DFTB and PM6) underestimate the repulsion, which leads to substantial overbinding at the short distance found in halogen bonds. The correction is thus repulsive.202 To develop and parametrize it, it was necessary to use benchmark data at very short distances. Therefore, we included geometries up 80% of the equilibrium distances in the X40x10 data set,187 which were used for the reparameterization of the halogen-bonding correction for PM6.203 The same data have recently been used for the correction of the description of halogen bonding in the SCCDFTB method.204 When all these corrections are used together, the resulting method can describe a wide range of noncovalent interactions, and because of the efficiency of the SQM methods, it can be routinely applied to very large systems (up to tens of thousands of atoms171). 5063


Chemical Reviews

Review

6. A SUMMARY

Table 7. Errors of Corrected SQM Methods in S66 and S22 Data Sets (RMSE in kcal/mol)a method

S66

S22

PM6 PM6-DH2 PM6-DH+ PM6-D3H4 DFTB DFTB-D DFTB-D,γ DFTB-DH2 DFTB-D3H4 RM1 RM1-D3H4 AM1 AM1-DH2 AM1-D3H4 PM3 PM3-D3H4 B3LYP/6-31G* TPSS/TZVP-D BLYP/def2QZVP-D3 MP2/cc-pVTZ

3.02 0.91 0.82 0.65 2.88 1.50 1.17 1.44 0.67 5.39 0.92 6.24 1.93 1.35 5.08 1.40 2.68 0.69 0.25 0.70

4.16 0.54 0.80 0.78 3.45 1.63 1.21 1.86 0.97 7.15 1.03 8.66 0.85 1.76 7.64 2.51 3.63 0.58 0.33 1.85

6.1. A State-of-the-Art Summary

After a period of fast development, the field of benchmark calculations of noncovalent interactions has reached a rather stable state and is close to saturation. The accuracy of the method that most often serves as the benchmark, CCSD(T)/ CBS, has been well-characterized and it is clear that interaction energies in systems of a few tens of atoms can be routinely evaluated with an accuracy of a few percent. Even more accurate calculations are now possible for smaller systems. Although there is still some room for improving the benchmark data by using more demanding calculations, the accuracy of the existing data sets is sufficient for the vast majority of applications. A large body of data, including data sets designed specifically for benchmarking, is now available. Moreover, it has been made freely accessible in online repositories, which make the use of benchmark data sets very easy. At present, these data sets are also being integrated directly into software packages, which makes their use even more straightforward. A very positive trend is the widespread use of a small number of established data sets by many authors. Most notably, the S66 data set has recently become a de facto standard in the field. This allows for an easy comparison of methods across various studies. Another important development is the appreciation of the importance of nonequilibrium geometries. Nowadays, many of the data sets contain them or are accompanied by them, most often in the form of dissociation curves. These nonequilibrium geometries are a crucial tool for the parametrization of empirical methods; recent progress in methods such as DFTD and empirically corrected semiempirical methods would not be possible without these data. Large-scale benchmarking also enabled very thorough quantification of the accuracy of more approximate methods. A hierarchy of methods with an optimal accuracy-to-cost ratio can be established and it is now easy to select the most suitable method for a particular application. Among the methods available, we would like to highlight the parametrized, correlated methods, such as SCS-MI-CCSD and MP2.5, which yield results close to the benchmark at a substantially reduced cost. Among more efficient methods, the modern versions of DFT-D can be recommended. At the end of the spectrum, empirically corrected semiempirical QM methods can now describe noncovalent interactions with reasonable accuracy at a very low cost.

a

Several DFT methods and MP2 are listed for comparison. Note that the D3H4 corrections were parameterized on the S66x8 data set and the DH2 and DH+ corrections on data including the S22 set.

The development and parametrization of these corrections would not be possible without a large body of benchmark data covering also nonequilibrium geometries. Since the SQM methods can have large and systematic errors, the corrections are very empirical and it is necessary to parametrize them on the whole relevant PES. The D3 correction has been parametrized on the dispersion-bound systems from the S66x8 set. The use of the dissociation curves in the parametrization of a pairwise potential ensures consistent accuracy in any geometry. The use of geometries closer than in equilibrium revealed insufficient hydrogen−hydrogen repulsion in the vast majority of the SQM methods. This effect was too strong to be compensated for by adjusting the dispersion correction. An extra repulsive term had to be added to the standard form of the D3 dispersion correction. It was later revised by improving the potential at even shorter distances;240 while this does not affect the description of realistic geometries, it is crucial for work with geometries obtained by computer modeling, where such close contacts are common artifacts (e.g., in docking or when hydrogens are added to X-ray structures). In the case of H-bond correction, there are further requirements on the benchmark data. Here, the correction is a function of both H-bond distance and angle, and its form is fully empirical. Therefore, not only the dissociation curves from S66x8 but also angular-displaced geometries from S66a8 have to be used. It is clear that the parametrization of these corrections puts very high demands on the benchmark data sets used and the success of the corrections, especially the latest D3H4 version, indicates that the data sets used (S66x8, S66a8, X40x10) provide sufficient coverage for the noncovalent interactions considered.

6.2. Open Questions

There are, however, still many open questions. In our opinion, the most important issues are now related to the coverage of the chemical space by the benchmark data. Noncovalent interactions in organic molecules are mapped exceptionally well and there is not much to add. On the other hand, interactions involving other elements and especially all their combinations are under-represented or missing. This certainly slows down the development of universally applicable empirical methods, where these data are necessary for parametrization. Addressing this issue is within current technical means, but it will involve an enormous amount of work. The second issue is related to the size of the model systems used in the common benchmark data sets. The coverage of larger systems is sparse and the accuracy of the benchmark data available for them is questionable. This issue is of high 5064


Chemical Reviews

Review

REFERENCES

importance in many applications, where it is necessary to work with large molecular systems accessible only with approximate methods. If the performance of these methods can be reliably quantified beforehand using accurate benchmark data relevant to the problem, the interpretation of the results will be far easier. Solving this issue is not easy, as the accurate calculations in large systems are technically demanding and may require additional approximations. We are, however, optimistic and we believe that the ongoing progress in the methodology will gradually make larger and larger systems accessible even to benchmark-quality calculations.

(1) Hobza, P.; Zahradník, R. Intermolecular interactions between medium-sized systems. Nonempirical and empirical calculations of interaction energies. Successes and failures. Chem. Rev. 1988, 88, 871− 897. (2) Hobza, P.; Šponer, J. Structure, Energetics, and Dynamics of the Nucleic Acid Base Pairs: Nonempirical Ab Initio Calculations. Chem. Rev. 1999, 99, 3247−3276. (3) Müller-Dethlefs, K.; Hobza, P. Noncovalent Interactions: A Challenge for Experiment and Theory. Chem. Rev. 2000, 100, 143− 168. (4) Hobza, P.; Zahradník, R.; Müller-Dethlefs, K. The world of noncovalent interactions: 2006. Collect. Czech. Chem. Commun. 2006, 71, 443−531. (5) Č erný, J.; Hobza, P. Non-covalent interactions in biomacromolecules. Phys. Chem. Chem. Phys. 2007, 9, 5291−5303. (6) Riley, K. E.; Pitoňaḱ , M.; Jurečka, P.; Hobza, P. Stabilization and Structure Calculations for Noncovalent Interactions in Extended Molecular Systems Based on Wave Function and Density Functional Theories. Chem. Rev. 2010, 110, 5023−5063. (7) Sherrill, C. D. In Reviews in Computational Chemistry; Lipkowitz, K. B., Cundari, T. R., Eds.; John Wiley & Sons, Inc., 2008; pp 1−38. (8) Hohenstein, E. G.; Sherrill, C. D. Wavefunction methods for noncovalent interactions. WIREs Comput. Mol. Sci. 2012, 2, 304−326. (9) Č ížek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256−4266. (10) Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. Rev. Mod. Phys. 2007, 79, 291−352. (11) Noga, J.; Bartlett, R. J. The full CCSDT model for molecular electronic structure. J. Chem. Phys. 1987, 86, 7041−7050. (12) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A. A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479−483. (13) Urban, M.; Noga, J.; Cole, S. J.; Bartlett, R. J. Towards a full CCSDT model for electron correlation. J. Chem. Phys. 1985, 83, 4041−4046. (14) Lee, Y. S.; Kucharski, S. A.; Bartlett, R. J. A coupled cluster approach with triple excitations. J. Chem. Phys. 1984, 81, 5906−5912. (15) Noga, J.; Bartlett, R. J.; Urban, M. Towards a full CCSDT model for electron correlation. CCSDT-n models. Chem. Phys. Lett. 1987, 134, 126−132. (16) Šimová, L.; Ř ezác,̌ J.; Hobza, P. Convergence of the Interaction Energies in Noncovalent Complexes in the Coupled-Cluster Methods Up to Full Configuration Interaction. J. Chem. Theory Comput. 2013, 9, 3420−3428. (17) Ř ezáč, J.; Šimová, L.; Hobza, P. CCSD[T] Describes Noncovalent Interactions Better than the CCSD(T), CCSD(TQ), and CCSDT Methods. J. Chem. Theory Comput. 2013, 9, 364−369. (18) Bomble, Y. J.; Stanton, J. F.; Kállay, M.; Gauss, J. Coupledcluster methods including noniterative corrections for quadruple excitations. J. Chem. Phys. 2005, 123, 054101. (19) London, F. Some characteristics and uses of molecular force. Z. Phys. Chem. B-Chem. Elem. Aufbau. Mater. 1930, 11, 222−251. (20) Bloom, J. W. G.; Wheeler, S. E. Taking the Aromaticity out of Aromatic Interactions. Angew. Chem. 2011, 123, 7993−7995. (21) Pimentel, G.; McClellan, A. The Hydrogen Bond; W. H. Freeman & Co., 1960. (22) Arunan, E.; Desiraju, G. R.; Klein, R. A.; Sadlej, J.; Scheiner, S.; Alkorta, I.; Clary, D. C.; Crabtree, R. H.; Dannenberg, J. J.; Hobza, P.; et al. Definition of the hydrogen bond (IUPAC Recommendations 2011). Pure Appl. Chem. 2011, 83, 1637−1641. (23) Isaacs, E. D.; Shukla, A.; Platzman, P. M.; Hamann, D. R.; Barbiellini, B.; Tulk, C. A. Covalency of the Hydrogen Bond in Ice: A Direct X-Ray Measurement. Phys. Rev. Lett. 1999, 82, 600−603. (24) Liu, Q.; Hoffmann, R. Theoretical Aspects of a Novel Mode of Hydrogen-Hydrogen Bonding. J. Am. Chem. Soc. 1995, 117, 10108− 10112.

AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. Biographies Jan Ř ezáč is senior researcher at the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Prague, Czech Republic, and associate professor at Palacky University in Olomouc, Czech Republic. He studied physical chemistry at Charles University in Prague, where he obtained his Ph.D. with Pavel Hobza in 2008. Then, he spent his postdoctoral study with D. R. Salahub at the University of Calgary, Calgary, Canada, and two shorter stays as a visiting professor in Puerto Rico and Paris. His main research interests are benchmark calculations of noncovalent interactions, improvement of their description in semiempirical QM methods, and software development. He has published more than 65 papers on these and related topics. Pavel Hobza is holder of the Distinguished Chair at Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Prague, Czech Republic, and is professor of chemistry at Charles University in Prague and Palacky University in Olomouc, Czech Republic. He obtained his Ph.D. degree with Rudolf Zahradnı ́k in 1974 at Academy of Sciences in Prague. After postdoctorial studies in Montreal, Canada (with Camille Sandorfy); Erlangen, Germany (with Paul von Rague Schleyer); and Munich, Germany (with Edward W. Schlag), he spent several periods as a visiting professor in Montreal, Munich, and Pohang, Korea. He has authored or coauthored more than 500 papers and 3 books. These studies focus mainly on noncovalent interactions and their role in chemistry, biodisciplines, and nanosciences. According Thomson Reuters, P.H. is ranked among the top 1% of researchers for most cited documents, in the field of chemistry in 2014 and 2015.

ACKNOWLEDGMENTS The work is part of research project RVO:61388963 of the IOCB AS CR. Funding for this work from the Czech Science Foundation (P208/12/G016) and from the Ministry of Education, Youth and Sports of the Czech Republic (project LO1305) is gratefully acknowledged. Some of the calculations presented in this work were carried out in the T4Innovations National Supercomputing Center, the Ministry of Education, Youth and Sports Large Infrastructures for Research, Experimental Development and Innovations project LM2015070. 5065


Chemical Reviews

Review

(25) Klooster, W. T.; Koetzle, T. F.; Siegbahn, P. E. M.; Richardson, T. B.; Crabtree, R. H. Study of the N-H ··· H-B dihydrogen bond including the crystal structure of BH3NH3 by neutron diffraction. J. Am. Chem. Soc. 1999, 121, 6337−6343. (26) Politzer, P.; Murray, J. S.; Clark, T. Halogen bonding and other σ-hole interactions: a perspective. Phys. Chem. Chem. Phys. 2013, 15, 11178−11189. (27) Brinck, T.; Murray, J.; Politzer, P. Surface Electrostatic Potentials of Halogenated Methanes as Indicators. Int. J. Quantum Chem. 1992, 44, 57−64. (28) Clark, T.; Hennemann, M.; Murray, J. S.; Politzer, P. Halogen bonding: the σ-hole. J. Mol. Model. 2007, 13, 291−296. (29) Metrangolo, P.; Resnati, G. Halogen Bonding: A Paradigm in Supramolecular Chemistry. Chem. - Eur. J. 2001, 7, 2511−2519. (30) Auffinger, P.; Hays, F. A.; Westhof, E.; Ho, P. S. Halogen bonds in biological molecules. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 16789−16794. (31) Politzer, P.; Lane, P.; Concha, M. C.; Ma, Y.; Murray, J. S. An overview of halogen bonding. J. Mol. Model. 2007, 13, 305−311. (32) Hobza, P.; Riley, K. The Relative Roles of Electrostatics and Dispersion in the Stabilization of Halogen Bonds. Phys. Chem. Chem. Phys. 2013, 15, 17742. (33) El Kerdawy, A.; Murray, J. S.; Politzer, P.; Bleiziffer, P.; Heßelmann, A.; Görling, A.; Clark, T. Directional Noncovalent Interactions: Repulsion and Dispersion. J. Chem. Theory Comput. 2013, 9, 2264. (34) Sedlák, R.; Kolár,̌ M. H.; Hobza, P. Polar Flattening and the Strength of Halogen Bonding. J. Chem. Theory Comput. 2015, 11, 4727−4732. (35) Scheiner, S. The Pnicogen Bond: Its Relation to Hydrogen, Halogen, and Other Noncovalent Bonds. Acc. Chem. Res. 2013, 46, 280−288. (36) Wang, W.; Ji, B.; Zhang, Y. Chalcogen Bond: A Sister Noncovalent Bond to Halogen Bond. J. Phys. Chem. A 2009, 113, 8132−8135. (37) Bauzá, A.; Frontera, A. Aerogen Bonding Interaction: A New Supramolecular Force? Angew. Chem., Int. Ed. 2015, 54, 7340−7343. (38) Kolár,̌ M. H.; Hobza, P. Computer Modeling of Halogen Bonds and Other σ-Hole Interactions. Chem. Rev. 2016, DOI: 10.1021/ acs.chemrev.5b00560. (39) Boys, S.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (40) Alvarez-Idaboy, J. R.; Galano, A. Counterpoise corrected interaction energies are not systematically better than uncorrected ones: comparison with CCSD(T) CBS extrapolated values. Theor. Chem. Acc. 2010, 126, 75−85. (41) Gutowski, M.; Van Duijneveldt, F. B.; Chałasiński, G.; Piela, L. Proper correction for the basis set superposition error in SCF calculations of intermolecular interactions. Mol. Phys. 1987, 61, 233− 247. (42) Burns, L. A.; Marshall, M. S.; Sherrill, C. D. Comparing Counterpoise-Corrected, Uncorrected, and Averaged Binding Energies for Benchmarking Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10, 49−57. (43) van Duijneveldt, F. B.; van Duijneveldt - van de Rijdt, J. G. C. M.; van Lenthe, J. H. State of the Art in Counterpoise Theory. Chem. Rev. 1994, 94, 1873−1885. (44) Mentel, Ł. M.; Baerends, E. J. Can the Counterpoise Correction for Basis Set Superposition Effect Be Justified? J. Chem. Theory Comput. 2014, 10, 252−267. (45) Halkier, A.; Klopper, W.; Helgaker, T.; Jørgensen, P.; Taylor, P. R. Basis set convergence of the interaction energy of hydrogen-bonded complexes. J. Chem. Phys. 1999, 111, 9157−9167. (46) Simon, S.; Duran, M.; Dannenberg, J. J. How does basis set superposition error change the potential surfaces for hydrogen-bonded dimers? J. Chem. Phys. 1996, 105, 11024−11031.

(47) Hobza, P.; Havlas, Z. Counterpoise-corrected potential energy surfaces of simple H-bonded systems. Theor. Chem. Acc. 1998, 99, 372−377. (48) 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. (49) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties. J. Chem. Phys. 1994, 100, 2975. (50) Papajak, E.; Zheng, J.; Xu, X.; Leverentz, H. R.; Truhlar, D. G. Perspectives on Basis Sets Beautiful: Seasonal Plantings of Diffuse Basis Functions. J. Chem. Theory Comput. 2011, 7, 3027−3034. (51) Woon, D. E.; Dunning, T. H., Jr. Gaussian basis sets for use in correlated molecular calculations. V. Core-valence basis sets for boron through neon. J. Chem. Phys. 1995, 103, 4572−4585. (52) Peterson, K. A.; Dunning, T. H. Accurate correlation consistent basis sets for molecular core-valence correlation effects: The second row atoms Al-Ar, and the first row atoms B-Ne revisited. J. Chem. Phys. 2002, 117, 10548−10560. (53) Peterson, K. A.; Figgen, D.; Goll, E.; Stoll, H.; Dolg, M. Systematically convergent basis sets with relativistic pseudopotentials. II. Small-core pseudopotentials and correlation consistent basis sets for the post-d group 16−18 elements. J. Chem. Phys. 2003, 119, 11113− 11123. (54) Weigend, F. A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285−4291. (55) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (56) Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116, 3175−3183. (57) Patkowski, K. On the accuracy of explicitly correlated coupledcluster interaction energies  have orbital results been beaten yet? J. Chem. Phys. 2012, 137, 034103. (58) Jurečka, P.; Nachtigall, P.; Hobza, P. RI-MP2 calculations with extended basis setsa promising tool for study of H-bonded and stacked DNA base pairs. Phys. Chem. Chem. Phys. 2001, 3, 4578−4582. (59) Ř ezác,̌ J.; Hobza, P. Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the “Gold Standard,” CCSD(T) at the Complete Basis Set Limit? J. Chem. Theory Comput. 2013, 9, 2151−2155. (60) Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theoret. Chim. Acta 1973, 28, 213−222. (61) Kroon-Batenburg, L.; van Duijneveldt, F. The use of a momentoptimized DZP basis set for describing the interaction in the water dimer. J. Mol. Struct.: THEOCHEM 1985, 121, 185−199. (62) Riley, K. E.; Hobza, P. Assessment of the MP2 Method, along with Several Basis Sets, for the Computation of Interaction Energies of Biologically Relevant Hydrogen Bonded and Dispersion Bound Complexes. J. Phys. Chem. A 2007, 111, 8257−8263. (63) Riley, K. E.; Ř ezác,̌ J.; Hobza, P. MP2.X: a generalized MP2.5 method that produces improved binding energies with smaller basis sets. Phys. Chem. Chem. Phys. 2011, 13, 21121−21125. (64) Ř ezác,̌ J.; Dubecký, M.; Jurečka, P.; Hobza, P. Extensions and applications of the A24 data set of accurate interaction energies. Phys. Chem. Chem. Phys. 2015, 17, 19268−19277. (65) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-set convergence of correlated calculations on water. J. Chem. Phys. 1997, 106, 9639−9646. (66) 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. (67) Truhlar, D. G. Basis-set extrapolation. Chem. Phys. Lett. 1998, 294, 45−48. (68) Bakowies, D. Accurate extrapolation of electron correlation energies from small basis sets. J. Chem. Phys. 2007, 127, 164109. 5066


Chemical Reviews

Review

(69) Okoshi, M.; Atsumi, T.; Nakai, H. Revisiting the extrapolation of correlation energies to complete basis set limit. J. Comput. Chem. 2015, 36, 1075−1082. (70) Min, S. K.; Lee, E. C.; Lee, H. M.; Kim, D. Y.; Kim, D.; Kim, K. S. Complete basis set limit of Ab initio binding energies and geometrical parameters for various typical types of complexes. J. Comput. Chem. 2008, 29, 1208−1221. (71) Hylleraas, E. A. Neue Berechnung der Energie des Heliums im Grundzustande, sowie des tiefsten Terms von Ortho-Helium. Eur. Phys. J. A 1929, 54, 347−366. (72) Kutzelnigg, W.; Klopper, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General theory. J. Chem. Phys. 1991, 94, 1985−2001. (73) Ten-no, S. Initiation of explicitly correlated Slater-type geminal theory. Chem. Phys. Lett. 2004, 398, 56−61. (74) May, A. J.; Manby, F. R. An explicitly correlated second order Møller-Plesset theory using a frozen Gaussian geminal. J. Chem. Phys. 2004, 121, 4479−4485. (75) 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. (76) Ten-no, S.; Noga, J. Explicitly correlated electronic structure theory from R12/F12 ansätze. WIREs Comput. Mol. Sci. 2012, 2, 114− 125. (77) Köhn, A.; Richings, G. W.; Tew, D. P. Implementation of the full explicitly correlated coupled-cluster singles and doubles model CCSD-F12 with optimally reduced auxiliary basis dependence. J. Chem. Phys. 2008, 129, 201103. (78) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Explicitly correlated coupled-cluster singles and doubles method based on complete diagrammatic equations. J. Chem. Phys. 2008, 129, 071101. (79) Knizia, G.; Adler, T. B.; Werner, H.-J. Simplified CCSD(T)-F12 methods: Theory and benchmarks. J. Chem. Phys. 2009, 130, 054104. (80) Köhn, A. Explicitly correlated connected triple excitations in coupled-cluster theory. J. Chem. Phys. 2009, 130, 131101. (81) Tew, D. P.; Klopper, W.; Neiss, C.; Hättig, C. Quintuple-ζ quality coupled-cluster correlation energies with triple-ζ basis sets. Phys. Chem. Chem. Phys. 2007, 9, 1921−1930. (82) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (83) Hättig, C.; Tew, D. P.; Köhn, A. Communications: Accurate and efficient approximations to explicitly correlated coupled-cluster singles and doubles, CCSD-F12. J. Chem. Phys. 2010, 132, 231102. (84) Peterson, K. A.; Adler, T. B.; Werner, H.-J. Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, B−Ne, and Al−Ar. J. Chem. Phys. 2008, 128, 084102. (85) Yousaf, K. E.; Peterson, K. A. Optimized auxiliary basis sets for explicitly correlated methods. J. Chem. Phys. 2008, 129, 184108. (86) Hill, J. G.; Peterson, K. A.; Knizia, G.; Werner, H.-J. Extrapolating MP2 and CCSD explicitly correlated correlation energies to the complete basis set limit with first and second row correlation consistent basis sets. J. Chem. Phys. 2009, 131, 194105. (87) Koch, H.; Fernández, B.; Christiansen, O. The benzene−argon complex: A ground and excited state ab initio study. J. Chem. Phys. 1998, 108, 2784. (88) Sinnokrot, M. O.; Valeev, E. F.; Sherrill, C. D. Estimates of the Ab Initio Limit for π-π Interactions: The Benzene Dimer. J. Am. Chem. Soc. 2002, 124, 10887−10893. (89) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. Origin of Attraction and Directionality of the π/π Interaction: Model Chemistry Calculations of Benzene Dimer Interaction. J. Am. Chem. Soc. 2002, 124, 104−112. (90) Hobza, P.; Šponer, J. Toward True DNA Base-Stacking Energies: MP2, CCSD(T), and Complete Basis Set Calculations. J. Am. Chem. Soc. 2002, 124, 11802−11808. (91) Cybulski, S. M.; Lytle, M. L. The origin of deficiency of the supermolecule second-order Møller-Plesset approach for evaluating interaction energies. J. Chem. Phys. 2007, 127, 141102.

(92) Jurečka, P.; Hobza, P. On the convergence of the (ΔECCSD(T)-ΔEMP2) term for complexes with multiple H-bonds. Chem. Phys. Lett. 2002, 365, 89−94. (93) Marshall, M. S.; Burns, L. A.; Sherrill, C. D. Basis set convergence of the coupled-cluster correction, δMP2CCSD(T): Best practices for benchmarking non-covalent interactions and the attendant revision of the S22, NBC10, HBC6, and HSG databases. J. Chem. Phys. 2011, 135, 194102. (94) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. S66: A Well-balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures. J. Chem. Theory Comput. 2011, 7, 2427−2438. (95) Werner, H.-J.; Knowles, P. J.; Manby, F. R.; Schütz, M.; et al. MOLPRO, a package of ab initio programs, version 2010.1; 2010; http://www.molpro.net (accessed Feb 10, 2016). (96) Ř ezác,̌ J.; Hobza, P. Ab Initio Quantum Mechanical Description of Noncovalent Interactions at Its Limits: Approaching the Experimental Dissociation Energy of the HF Dimer. J. Chem. Theory Comput. 2014, 10, 3066−3073. (97) Kállay, M.; Surján, P. R. Computing coupled-cluster wave functions with arbitrary excitations. J. Chem. Phys. 2000, 113, 1359− 1365. (98) Hirata, S.; Bartlett, R. J. High-order coupled-cluster calculations through connected octuple excitations. Chem. Phys. Lett. 2000, 321, 216−224. (99) Olsen, J. The initial implementation and applications of a general active space coupled cluster method. J. Chem. Phys. 2000, 113, 7140−7148. (100) Kállay, M.; Surján, P. R. Higher excitations in coupled-cluster theory. J. Chem. Phys. 2001, 115, 2945−2954. (101) Jäger, B.; Hellmann, R.; Bich, E.; Vogel, E. Ab initio pair potential energy curve for the argon atom pair and thermophysical properties of the dilute argon gas. I. Argon−argon interatomic potential and rovibrational spectra. Mol. Phys. 2009, 107, 2181−2188. (102) Hellmann, R.; Bich, E.; Vogel, E. Ab initio potential energy curve for the neon atom pair and thermophysical properties of the dilute neon gas. I. Neon−neon interatomic potential and rovibrational spectra. Mol. Phys. 2008, 106, 133−140. (103) Patkowski, K.; Szalewicz, K. Argon pair potential at basis set and excitation limits. J. Chem. Phys. 2010, 133, 094304. (104) Raghavachari, K.; Pople, J. A.; Replogle, E. S.; Head-Gordon, M. Fifth order Moeller-Plesset perturbation theory: comparison of existing correlation methods and implementation of new methods correct to fifth order. J. Phys. Chem. 1990, 94, 5579−5586. (105) Demovičová, L.; Hobza, P.; Ř ezác,̌ J. Evaluation of composite schemes for CCSDT(Q) calculations of interaction energies of noncovalent complexes. Phys. Chem. Chem. Phys. 2014, 16, 19115− 19121. (106) Smith, D. G. A.; Jankowski, P.; Slawik, M.; Witek, H. A.; Patkowski, K. Basis Set Convergence of the Post-CCSD(T) Contribution to Noncovalent Interaction Energies. J. Chem. Theory Comput. 2014, 10, 3140−3150. (107) Douglas, M.; Kroll, N. M. Quantum electrodynamical corrections to the fine structure of helium. Ann. Phys. 1974, 82, 89− 155. (108) Hess, B. A. Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33, 3742−3748. (109) Reiher, M. Relativistic Douglas−Kroll−Hess theory. WIREs Comput. Mol. Sci. 2012, 2, 139−149. (110) de Jong, W. A.; Harrison, R. J.; Dixon, D. A. Parallel Douglas− Kroll energy and gradients in NWChem: Estimating scalar relativistic effects using Douglas−Kroll contracted basis sets. J. Chem. Phys. 2001, 114, 48−53. (111) Kutzelnigg, W. The adiabatic approximation I. The physical background of the Born-Handy ansatz. Mol. Phys. 1997, 90, 909−916. (112) Gauss, J.; Tajti, A.; Kállay, M.; Stanton, J. F.; Szalay, P. G. Analytic calculation of the diagonal Born-Oppenheimer correction 5067


Chemical Reviews

Review

within configuration-interaction and coupled-cluster theory. J. Chem. Phys. 2006, 125, 144111. (113) Feller, D.; Peterson, K. A.; Dixon, D. A. A survey of factors contributing to accurate theoretical predictions of atomization energies and molecular structures. J. Chem. Phys. 2008, 129, 204105. (114) Császár, A. G.; Allen, W. D.; Schaefer, H. F. In pursuit of the ab initio limit for conformational energy prototypes. J. Chem. Phys. 1998, 108, 9751−9764. (115) Tajti, A.; Szalay, P. G.; Császár, A. G.; Kállay, M.; Gauss, J.; Valeev, E. F.; Flowers, B. A.; Vázquez, J.; Stanton, J. F. HEAT: High accuracy extrapolated ab initio thermochemistry. J. Chem. Phys. 2004, 121, 11599−11613. (116) Bomble, Y. J.; Vázquez, J.; Kállay, M.; Michauk, C.; Szalay, P. G.; Császár, A. G.; Gauss, J.; Stanton, J. F. High-accuracy extrapolated ab initio thermochemistry. II. Minor improvements to the protocol and a vital simplification. J. Chem. Phys. 2006, 125, 064108. (117) Harding, M. E.; Vázquez, J.; Ruscic, B.; Wilson, A. K.; Gauss, J.; Stanton, J. F. High-accuracy extrapolated ab initio thermochemistry. III. Additional improvements and overview. J. Chem. Phys. 2008, 128, 114111. (118) Boese, A. D.; Oren, M.; Atasoylu, O.; Martin, J. M. L.; Kállay, M.; Gauss, J. W3 theory: Robust computational thermochemistry in the kJ/mol accuracy range. J. Chem. Phys. 2004, 120, 4129−4141. (119) Karton, A.; Rabinovich, E.; Martin, J. M. L.; Ruscic, B. W4 theory for computational thermochemistry: In pursuit of confident sub-kJ/mol predictions. J. Chem. Phys. 2006, 125, 144108. (120) Boese, A. D. Assessment of Coupled Cluster Theory and more Approximate Methods for Hydrogen Bonded Systems. J. Chem. Theory Comput. 2013, 9, 4403−4413. (121) Harding, M. E.; Klopper, W. Benchmarking the Lithium− Thiophene Complex. ChemPhysChem 2013, 14, 708−715. (122) Jankowski, P.; McKellar, A. R. W.; Szalewicz, K. Theory Untangles the High-Resolution Infrared Spectrum of the ortho-H2CO van der Waals Complex. Science 2012, 336, 1147−1150. (123) Taube, A. G.; Bartlett, R. J. Frozen natural orbitals: Systematic basis set truncation for coupled-cluster theory. Collect. Czech. Chem. Commun. 2005, 70, 837−850. (124) DePrince, A. E.; Sherrill, C. D. Accurate Noncovalent Interaction Energies Using Truncated Basis Sets Based on Frozen Natural Orbitals. J. Chem. Theory Comput. 2013, 9, 293−299. (125) Adamowicz, L.; Bartlett, R. J. Optimized virtual orbital space for high-level correlated calculations. J. Chem. Phys. 1987, 86, 6314− 6324. (126) Pitoňaḱ , M.; Neogrády, P.; Kellö, V.; Urban, M. Optimized virtual orbitals for relativistic calculations: an alternative approach to the basis set construction for correlation calculations. Mol. Phys. 2006, 104, 2277−2292. (127) Kraus, M.; Pitoňaḱ , M.; Hobza, P.; Urban, M.; Neogrády, P. Highly correlated calculations using optimized virtual orbital space with controlled accuracy. Application to counterpoise corrected interaction energy calculations. Int. J. Quantum Chem. 2012, 112, 948−959. (128) Sæbø, S.; Pulay, P. Local configuration interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13−18. (129) Sæbø, S.; Pulay, P. Local Treatment of Electron Correlation. Annu. Rev. Phys. Chem. 1993, 44, 213−236. (130) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286−6297. (131) Scuseria, G. E.; Ayala, P. Y. Linear scaling coupled cluster and perturbation theories in the atomic orbital basis. J. Chem. Phys. 1999, 111, 8330−8343. (132) Schütz, M.; Werner, H.-J. Local perturbative triples correction (T) with linear cost scaling. Chem. Phys. Lett. 2000, 318, 370−378. (133) Rolik, Z.; Kállay, M. A general-order local coupled-cluster method based on the cluster-in-molecule approach. J. Chem. Phys. 2011, 135, 104111.

(134) Rolik, Z.; Szegedy, L.; Ladjánszki, I.; Ladóczki, B.; Kállay, M. An efficient linear-scaling CCSD(T) method based on local natural orbitals. J. Chem. Phys. 2013, 139, 094105. (135) Sparta, M.; Neese, F. Chemical applications carried out by local pair natural orbital based coupled-cluster methods. Chem. Soc. Rev. 2014, 43, 5032. (136) 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. (137) Dubecký, S.; Mitas, L.; Jurecka, P. Noncovalent Interactions by Quantum Monte Carlo. Chem. Rev. 2016, in this issue. (138) Dubecký, M.; Jurečka, P.; Derian, R.; Hobza, P.; Otyepka, M.; Mitas, L. Quantum Monte Carlo Methods Describe Noncovalent Interactions with Subchemical Accuracy. J. Chem. Theory Comput. 2013, 9, 4287−4292. (139) Ambrosetti, A.; Alfè, D.; DiStasio, R. A.; Tkatchenko, A. Hard Numbers for Large Molecules: Toward Exact Energetics for Supramolecular Systems. J. Phys. Chem. Lett. 2014, 5, 849−855. (140) 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. (141) DiStasio, R. A., Jr.; Head-Gordon, M. Optimized spincomponent scaled second-order M?ller-Plesset perturbation theory for intermolecular interaction energies. Mol. Phys. 2007, 105, 1073− 1083. (142) Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled opposite-spin second order Møller−Plesset correlation energy: An economical electronic structure method. J. Chem. Phys. 2004, 121, 9793−9802. (143) Grimme, S.; Goerigk, L.; Fink, R. F. Spin-component-scaled electron correlation methods. WIREs Comput. Mol. Sci. 2012, 2, 886− 906. (144) Marchetti, O.; Werner, H.-J. Accurate Calculations of Intermolecular Interaction Energies Using Explicitly Correlated Coupled Cluster Wave Functions and a Dispersion-Weighted MP2 Method†. J. Phys. Chem. A 2009, 113, 11580−11585. (145) Goldey, M.; Head-Gordon, M. Attenuating Away the Errors in Inter- and Intramolecular Interactions from Second-Order Møller− Plesset Calculations in the Small Aug-cc-pVDZ Basis Set. J. Phys. Chem. Lett. 2012, 3, 3592−3598. (146) Goldey, M.; Dutoi, A.; Head-Gordon, M. Attenuated secondorder Møller−Plesset perturbation theory: performance within the aug-cc-pVTZ basis. Phys. Chem. Chem. Phys. 2013, 15, 15869−15875. (147) Heßelmann, A. Improved supermolecular second order Møller−Plesset intermolecular interaction energies using time-dependent density functional response theory. J. Chem. Phys. 2008, 128, 144112. (148) Pitoň aḱ , M.; Heßelmann, A. Accurate Intermolecular Interaction Energies from a Combination of MP2 and TDDFT Response Theory. J. Chem. Theory Comput. 2010, 6, 168−178. (149) Ohnishi, Y.-y.; Ishimura, K.; Ten-no, S. Interaction Energy of Large Molecules from Restrained Denominator MP2-F12. J. Chem. Theory Comput. 2014, 10, 4857−4861. (150) Pitoňaḱ , M.; Neogrády, P.; Č erný, J.; Grimme, S.; Hobza, P. Scaled MP3 Non-Covalent Interaction Energies Agree Closely with Accurate CCSD(T) Benchmark Data. ChemPhysChem 2009, 10, 282− 289. (151) Sedlák, R.; Riley, K. E.; Ř ezác,̌ J.; Pitoňaḱ , M.; Hobza, P. MP2.5 and MP2.X: Approaching CCSD(T) Quality Description of Noncovalent Interaction at the Cost of a Single CCSD Iteration. ChemPhysChem 2013, 14, 698−707. (152) Ř ezác,̌ J.; Huang, Y.; Hobza, P.; Beran, G. J. O. Benchmark Calculations of Three-Body Intermolecular Interactions and the Performance of Low-Cost Electronic Structure Methods. J. Chem. Theory Comput. 2015, 11, 3065−3079. (153) Takatani, T.; Hohenstein, E. G.; Sherrill, C. D. Improvement of the coupled-cluster singles and doubles method via scaling same5068


Chemical Reviews

Review

and opposite-spin components of the double excitation correlation energy. J. Chem. Phys. 2008, 128, 124111. (154) Pitoňaḱ , M.; Ř ezác,̌ J.; Hobza, P. Spin-component scaled coupled-clusters singles and doubles optimized towards calculation of noncovalent interactions. Phys. Chem. Chem. Phys. 2010, 12, 9611. (155) Wu, Q.; Yang, W. Empirical correction to density functional theory for van der Waals interactions. J. Chem. Phys. 2002, 116, 515− 524. (156) Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463−1473. (157) Jurečka, P.; Č erný, J.; Hobza, P.; Salahub, D. Density functional theory augmented with an empirical dispersion term. Interaction energies and geometries of 80 noncovalent complexes compared with ab initio quantum mechanics calculations. J. Comput. Chem. 2007, 28, 555−569. (158) 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. (159) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456−1465. (160) Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking Density Functional Methods against the S66 and S66x8 Datasets for NonCovalent Interactions. ChemPhysChem 2011, 12, 3421−3433. (161) Grimme, S.; Hansen, A.; Brandenburg, J.; Bannwarth, C. Dispersion-Corrected Mean-Field Electronic Structure Methods. Chem. Rev. 2016, in this issue. (162) Eshuis, H.; Bates, J. E.; Furche, F. Electron correlation methods based on the random phase approximation. Theor. Chem. Acc. 2012, 131, 1084. (163) Christensen, A.; Kubar, T.; Cui, Q.; Elstner, M. Semi-Empirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications. Chem. Rev. 2016, in this issue. (164) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen bonding and stacking interactions of nucleic acid base pairs: A density-functional-theory based treatment. J. Chem. Phys. 2001, 114, 5149. (165) Martin, B.; Clark, T. Dispersion treatment for NDDO-based semiempirical MO techniques. Int. J. Quantum Chem. 2006, 106, 1208−1216. (166) McNamara, J. P.; Hillier, I. H. Semi-empirical molecular orbital methods including dispersion corrections for the accurate prediction of the full range of intermolecular interactions in biomolecules. Phys. Chem. Chem. Phys. 2007, 9, 2362. (167) Morgado, C. A.; McNamara, J. P.; Hillier, I. H.; Burton, N. A.; Vincent, M. A. Density Functional and Semiempirical Molecular Orbital Methods Including Dispersion Corrections for the Accurate Description of Noncovalent Interactions Involving Sulfur-Containing Molecules. J. Chem. Theory Comput. 2007, 3, 1656−1664. (168) Ř ezác,̌ J.; Fanfrlík, J.; Salahub, D.; Hobza, P. Semiempirical Quantum Chemical PM6 Method Augmented by Dispersion and HBonding Correction Terms Reliably Describes Various Types of Noncovalent Complexes. J. Chem. Theory Comput. 2009, 5, 1749− 1760. (169) Korth, M.; Pitoňaḱ , M.; Ř ezác,̌ J.; Hobza, P. A Transferable HBonding Correction for Semiempirical Quantum-Chemical Methods. J. Chem. Theory Comput. 2010, 6, 344−352. (170) Ř ezác,̌ J.; Hobza, P. Advanced Corrections of Hydrogen Bonding and Dispersion for Semiempirical Quantum Mechanical Methods. J. Chem. Theory Comput. 2012, 8, 141−151. (171) Stewart, J. J. P. Application of the PM6 method to modeling proteins. J. Mol. Model. 2009, 15, 765−805. (172) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (173) Stone, A. The Theory of Intermolecular Forces, 2nd ed.; Oxford University Press: Oxford, 2013.

(174) Monkhorst, H. J.; Jeziorski, B.; Harris, F. E. Recursive scheme for order-by-order many-body perturbation theory. Phys. Rev. A: At., Mol., Opt. Phys. 1981, 23, 1639−1644. (175) Williams, H. L.; Chabalowski, C. F. Using Kohn-Sham Orbitals in Symmetry-Adapted Perturbation Theory to Investigate Intermolecular Interactions. J. Phys. Chem. A 2001, 105, 646−659. (176) Heßelmann, A.; Jansen, G. First-order intermolecular interaction energies from Kohn-Sham orbitals. Chem. Phys. Lett. 2002, 357, 464−470. (177) Heßelmann, A.; Jansen, G. Intermolecular induction and exchange-induction energies from coupled-perturbed Kohn-Sham density functional theory. Chem. Phys. Lett. 2002, 362, 319−325. (178) Heßelmann, A.; Jansen, G. Intermolecular dispersion energies from time-dependent density functional theory. Chem. Phys. Lett. 2003, 367, 778−784. (179) Misquitta, A. J.; Szalewicz, K. Intermolecular forces from asymptotically corrected density functional description of monomers. Chem. Phys. Lett. 2002, 357, 301−306. (180) Misquitta, A. J.; Szalewicz, K. Symmetry-adapted perturbationtheory calculations of intermolecular forces employing densityfunctional description of monomers. J. Chem. Phys. 2005, 122, 214109. (181) Szalewicz, K. Symmetry-adapted perturbation theory of intermolecular forces. WIREs Comput. Mol. Sci. 2012, 2, 254−272. (182) Podeszwa, R.; Bukowski, R.; Szalewicz, K. Potential Energy Surface for the Benzene Dimer and Perturbational Analysis of π-π Interactions. J. Phys. Chem. A 2006, 110, 10345−10354. (183) Tekin, A.; Jansen, G. How accurate is the density functional theory combined with symmetry-adapted perturbation theory approach for CH−π and π−π interactions? A comparison to supermolecular calculations for the acetylene−benzene dimer. Phys. Chem. Chem. Phys. 2007, 9, 1680−1687. (184) Podeszwa, R. Interactions of graphene sheets deduced from properties of polycyclic aromatic hydrocarbons. J. Chem. Phys. 2010, 132, 044704. (185) Witte, J.; Goldey, M.; Neaton, J. B.; Head-Gordon, M. Beyond Energies: Geometries of Nonbonded Molecular Complexes as Metrics for Assessing Electronic Structure Approaches. J. Chem. Theory Comput. 2015, 11, 1481−1492. (186) Howard, J. C.; Tschumper, G. S. Benchmark Structures and Harmonic Vibrational Frequencies Near the CCSD(T) Complete Basis Set Limit for Small Water Clusters: (H2O)n = 2, 3, 4, 5, 6. J. Chem. Theory Comput. 2015, 11, 2126−2136. (187) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. Benchmark calculations of noncovalent interactions of halogenated molecules. J. Chem. Theory Comput. 2012, 8, 4285−4292. (188) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. Evaluation of the performance of post-Hartree-Fock methods in terms of intermolecular distance in noncovalent complexes. J. Comput. Chem. 2012, 33, 691−694. (189) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. Extensions of the S66 Data Set: More Accurate Interaction Energies and Angular-Displaced Nonequilibrium Geometries. J. Chem. Theory Comput. 2011, 7, 3466−3470. (190) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; et al. Psi4: an open-source ab initio electronic structure program. WIREs Comput. Mol. Sci. 2012, 2, 556−565. (191) Ř ezác,̌ J. CubyRuby Framework for Computational Chemistry, version 4; http://cuby4.molecular.cz/ (accessed Feb 10, 2016). (192) Ř ezác,̌ J. Cuby: An integrative framework for computational chemistry. J. Comput. Chemm 2016, DOI: 10.1002/jcc.24312. (193) Takatani, T.; Hohenstein, E. G.; Malagoli, M.; Marshall, M. S.; Sherrill, C. D. Basis set consistent revision of the S22 test set of noncovalent interaction energies. J. Chem. Phys. 2010, 132, 144104. (194) 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. (195) Zhao, Y.; Truhlar, D. G. Benchmark Databases for Nonbonded Interactions and Their Use To Test Density Functional Theory. J. Chem. Theory Comput. 2005, 1, 415−432. 5069


Chemical Reviews

Review

(196) Martin, J. M. L.; de Oliveira, G. Towards standard methods for benchmark quality ab initio thermochemistryW1 and W2 theory. J. Chem. Phys. 1999, 111, 1843−1856. (197) Jurečka, P.; Šponer, J.; Č erný, J.; Hobza, P. Benchmark database of accurate (MP2 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. (198) Podeszwa, R.; Patkowski, K.; Szalewicz, K. Improved interaction energy benchmarks for dimers of biological relevance. Phys. Chem. Chem. Phys. 2010, 12, 5974. (199) Molnar, L. F.; He, X.; Wang, B.; Merz, K. M., Jr Further analysis and comparative study of intermolecular interactions using dimers from the S22 database. J. Chem. Phys. 2009, 131, 065102. (200) Gráfová, L.; Pitoňaḱ , M.; Ř ezác,̌ J.; Hobza, P. Comparative Study of Selected Wave Function and Density Functional Methods for Noncovalent Interaction Energy Calculations Using the Extended S22 Data Set. J. Chem. Theory Comput. 2010, 6, 2365−2376. (201) Mintz, B. J.; Parks, J. M. Benchmark Interaction Energies for Biologically Relevant Noncovalent Complexes Containing Divalent Sulfur. J. Phys. Chem. A 2012, 116, 1086−1092. (202) Ř ezác,̌ J.; Hobza, P. A halogen-bonding correction for the semiempirical PM6 method. Chem. Phys. Lett. 2011, 506, 286−289. (203) Brahmkshatriya, P. S.; Dobeš, P.; Fanfrlík, J.; Ř ezác,̌ J.; Paruch, K.; Bronowska, A.; Lepšík, M.; Hobza, P. Quantum Mechanical Scoring: Structural and Energetic Insights into Cyclin-Dependent Kinase 2 Inhibition by Pyrazolo [1,5-a]pyrimidines. Curr. Comput.Aided Drug Des. 2013, 9, 118−129. (204) Kubillus, M.; Kubař, T.; Gaus, M.; Ř ezác,̌ J.; Elstner, M. Parameterization of the DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332−342. (205) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. An Assessment of Theoretical Methods for Nonbonded Interactions: Comparison to Complete Basis Set Limit Coupled-Cluster Potential Energy Curves for the Benzene Dimer, the Methane Dimer, Benzene-Methane, and Benzene-H2S†. J. Phys. Chem. A 2009, 113, 10146−10159. (206) Granatier, J.; Pitoňaḱ , M.; Hobza, P. Accuracy of Several Wave Function and Density Functional Theory Methods for Description of Noncovalent Interaction of Saturated and Unsaturated Hydrocarbon Dimers. J. Chem. Theory Comput. 2012, 8, 2282−2292. (207) Thanthiriwatte, K. S.; Hohenstein, E. G.; Burns, L. A.; Sherrill, C. D. Assessment of the Performance of DFT and DFT-D Methods for Describing Distance Dependence of Hydrogen-Bonded Interactions. J. Chem. Theory Comput. 2011, 7, 88−96. (208) Bauzá, A.; Alkorta, I.; Frontera, A.; Elguero, J. On the Reliability of Pure and Hybrid DFT Methods for the Evaluation of Halogen, Chalcogen, and Pnicogen Bonds Involving Anionic and Neutral Electron Donors. J. Chem. Theory Comput. 2013, 9, 5201− 5210. (209) Karthikeyan, S.; Sedlák, R.; Hobza, P. On the Nature of Stabilization in Weak, Medium, and Strong Charge-Transfer Complexes: CCSD(T)/CBS and SAPT Calculations. J. Phys. Chem. A 2011, 115, 9422−9428. (210) Ř ezác,̌ J.; de la Lande, A. Robust, Basis-Set Independent Method for the Evaluation of Charge-Transfer Energy in Noncovalent Complexes. J. Chem. Theory Comput. 2015, 11, 528−537. (211) Berka, K.; Laskowski, R.; Riley, K. E.; Hobza, P.; Vondrásě k, J. Representative Amino Acid Side Chain Interactions in Proteins. A Comparison of Highly Accurate Correlated ab Initio Quantum Chemical and Empirical Potential Procedures. J. Chem. Theory Comput. 2009, 5, 982−992. (212) Hostaš, J.; Jakubec, D.; Laskowski, R. A.; Gnanasekaran, R.; Ř ezác,̌ J.; Vondrásě k, J.; Hobza, P. Representative Amino Acid SideChain Interactions in Protein−DNA Complexes: A Comparison of Highly Accurate Correlated Ab Initio Quantum Mechanical Calculations and Efficient Approaches for Applications to Large Systems. J. Chem. Theory Comput. 2015, 11, 4086. (213) Faver, J. C.; Benson, M. L.; He, X.; Roberts, B. P.; Wang, B.; Marshall, M. S.; Kennedy, M. R.; Sherrill, C. D.; Merz, K. M. Formal

Estimation of Errors in Computed Absolute Interaction Energies of Protein-Ligand Complexes. J. Chem. Theory Comput. 2011, 7, 790− 797. (214) Bryantsev, V. S.; Diallo, M. S.; van Duin, A. C. T.; Goddard, W. A. Evaluation of B3LYP, X3LYP, and M06-Class Density Functionals for Predicting the Binding Energies of Neutral, Protonated, and Deprotonated Water Clusters. J. Chem. Theory Comput. 2009, 5, 1016−1026. (215) Anacker, T.; Friedrich, J. New accurate benchmark energies for large water clusters: DFT is better than expected. J. Comput. Chem. 2014, 35, 634−643. (216) Temelso, B.; Renner, C. R.; Shields, G. C. Importance and Reliability of Small Basis Set CCSD(T) Corrections to MP2 Binding and Relative Energies of Water Clusters. J. Chem. Theory Comput. 2015, 11, 1439−1448. (217) Sedlák, R.; Janowski, T.; Pitoňaḱ , M.; Ř ezác,̌ J.; Pulay, P.; Hobza, P. Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes. J. Chem. Theory Comput. 2013, 9, 3364− 3374. (218) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - Eur. J. 2012, 18, 9955−9964. (219) Risthaus, T.; Grimme, S. Benchmarking of London DispersionAccounting Density Functional Theory Methods on Very Large Molecular Complexes. J. Chem. Theory Comput. 2013, 9, 1580−1591. (220) Sure, R.; Grimme, S. Comprehensive Benchmark of Association (Free) Energies of Realistic Host−Guest Complexes. J. Chem. Theory Comput. 2015, 11, 3785−3801. (221) Scuseria, G. E. Analytic evaluation of energy gradients for the singles and doubles coupled cluster method including perturbative triple excitations: Theory and applications to FOOF and Cr2. J. Chem. Phys. 1991, 94, 442−447. (222) Barone, V.; Biczysko, M.; Bloino, J.; Puzzarini, C. Characterization of the Elusive Conformers of Glycine from State-of-the-Art Structural, Thermodynamic, and Spectroscopic Computations: Theory Complements Experiment. J. Chem. Theory Comput. 2013, 9, 1533− 1547. (223) Barone, V.; Biczysko, M.; Bloino, J.; Puzzarini, C. Accurate structure, thermodynamic and spectroscopic parameters from CC and CC/DFT schemes: the challenge of the conformational equilibrium in glycine. Phys. Chem. Chem. Phys. 2013, 15, 10094−10111. (224) Riffet, V.; Bouchoux, G. Gas-phase structures and thermochemistry of neutral histidine and its conjugated acid and base. Phys. Chem. Chem. Phys. 2013, 15, 6097−6106. (225) Fogueri, U. R.; Kozuch, S.; Karton, A.; Martin, J. M. 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. (226) Mládek, A.; Krepl, M.; Svozil, D.; Č ech, P.; Otyepka, M.; Banás,̌ P.; Zgarbová, M.; Jurečka, P.; Šponer, J. Benchmark quantumchemical calculations on a complete set of rotameric families of the DNA sugar−phosphate backbone and their comparison with modern density functional theory. Phys. Chem. Chem. Phys. 2013, 15, 7295− 7310. (227) Valdes, H.; Pluhácǩ ová, K.; Pitonák, M.; Ř ezác,̌ J.; Hobza, P. Benchmark database on isolated small peptides containing an aromatic side chain: comparison between wave function and density functional theory methods and empirical force field. Phys. Chem. Chem. Phys. 2008, 10, 2747−2757. (228) BEGDB: Benchmark Energy and Geometry Database. http:// www.begdb.com/ (accessed Nov 20, 2013). (229) Ř ezác,̌ J.; Jurečka, P.; Riley, K. E.; Č erný, J.; Valdes, H.; Pluhácǩ ová, K.; Berka, K.; Ř ezác,̌ T.; Pitoňaḱ , M.; Vondrásě k, J. 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, www.begdb.com,. (230) Goerigk, L.; Grimme, S. A General Database for Main Group Thermochemistry, Kinetics, and Noncovalent Interactions - Assess5070


Chemical Reviews

Review

ment of Common and Reparameterized (meta-)GGA Density Functionals. J. Chem. Theory Comput. 2010, 6, 107−126. (231) Goerigk, L.; Grimme, S. Efficient and Accurate Double-HybridMeta-GGA Density FunctionalsEvaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 291−309. (232) Grimme, S. GMTKN30 database. http://www.thch.uni-bonn. de/tc/index.php?section=downloads&subsection=GMTKN30&lang= english (accessed Feb 10, 2016). (233) Minnesota Databases 2.0. http://comp.chem.umn.edu/db/ (accessed Feb 10, 2016). (234) Liakos, D. G.; Izsák, R.; Valeev, E. F.; Neese, F. What is the most efficient way to reach the canonical MP2 basis set limit? Mol. Phys. 2013, 111, 2653−2662. (235) TURBOMOLE, v6.5.; 2011; http://www.turbomole.com (accessed Feb 10, 2016). (236) Neese, F. The ORCA program system. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (237) Kállay, M.; Rolik, Z.; Csontos, J.; Ladjánszki, I.; Szegedy, L.; Ladóczki, B.; Samu, G. MRCC, a quantum chemical program suite; http://www.mrcc.hu (accessed Feb 10, 2016). (238) Korth, M. Third-Generation Hydrogen-Bonding Corrections for Semiempirical QM Methods and Force Fields. J. Chem. Theory Comput. 2010, 6, 3808−3816. (239) Yang; York, D.; Cui, Q.; Elstner, M.; Yu, H. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (240) Vorlová, B.; Nachtigallová, D.; Jirásková-Vaníčková, J.; Ajani, H.; Jansa, P.; Ř ezác,̌ J.; Fanfrlík, J.; Otyepka, M.; Hobza, P.; Konvalinka, J.; et al. Malonate-based inhibitors of mammalian serine racemase: Kinetic characterization and structure-based computational study. Eur. J. Med. Chem. 2015, 89, 189−197.

5071


Benchmark Calculations of Interaction Energies in Noncovalent

Recommend Documents