Empirical D3 Dispersion as a Replacement for ab Initio Dispersion

Subscriber access provided by University of Colorado Boulder

Article

Empirical D3 dispersion as a replacement for ab initio dispersion terms in the DFT-SAPT. Robert Sedlak, and Jan #ezá# J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01198 • Publication Date (Web): 07 Mar 2017 Downloaded from http://pubs.acs.org on March 8, 2017

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

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

Page 1 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Empirical D3 dispersion as a replacement for ab initio dispersion terms in the DFT-SAPT. Robert Sedlak1,2, Jan Řezáč1,* 1

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

Abstract In DFT-SAPT interaction energy calculations, the most demanding step is the calculation of the London dispersion term. To avoid this bottleneck and make DFT-SAPT applicable to larger systems, the ab initio dispersion terms can be replaced by one calculated empirically at an almost negligible cost (J. Phys. Chem. A 2011, 115, 11321–11330). We present an update of this approach that improves accuracy and makes the method applicable to a wider range of systems. It is based on Grimme’s D3 dispersion correction for DFT, where the damping function is changed to one suitable for the calculation of the complete dispersion energy. The best results have been achieved with the Tang-Toennies damping function. It has been parameterized on the S66x8 data set, for which we report DF-DFT-SAPT/aug-cc-pVTZ interaction energy decomposition. The method has been validated on a diverse set of noncovalent systems including difficult cases such as very compact noncovalent complexes of charge-transfer type. The root-mean-square error in the complete test set is 0.73 kcal.mol-1 and 0.42 kcal.mol-1 when charge-transfer complexes are excluded. The proposed empirical dispersion terms can also be used outside the DFT-SAPT framework, e.g. for the estimation of the amount of dispersion in a calculation where only the total interaction energy is known.

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Noncovalent interactions are responsible for very important phenomena in various areas of chemistry and physics and their understanding is necessary for the explanation of many processes at the molecular level. Both experiments and general computational methods can usually provide only the information on their energetics. Additional insight into the underlying mechanisms of the interaction is harder to obtain, but very valuable for their interpretation. On the side of theory, various energy decomposition analysis (EDA) methods have been developed to provide such information.1–5 Although the components obtained from EDA are not observable and their definition is more or less arbitrary, they are very useful tools that allow quantitative analysis and the categorization of noncovalent interactions, which are in line with general theoretical concepts. One of the most popular perturbative energy decomposition methods is the symmetry-adapted perturbation theory (SAPT).6 To obtain accurate results, the intermolecular perturbation has to be applied on top of a correlated calculation of the monomers and the resulting double perturbative variant of the SAPT is computationally extremely demanding.6 The local correlation can be, however, described well by the density functional theory (DFT). The DFT-based SAPT method (DFT-SAPT)7–9 can be further accelerated using the density fitting approximation (DF-DFT-SAPT), giving rise to an accurate, yet computationally efficient method.10 There exists an essentially identical approach to DFT-SAPT called SAPT(DFT)11–15 and its density-fitted version.16,17 The scaling of the DF-DFT-SAPT compared to the original SAPT method, with respect to the system size, is two orders of magnitude lower (i.e. N5 vs. N7, where N reflects the size of the studied system), which makes it routinely applicable to systems with dozens of atoms. In the most widely used form of DFT-SAPT, the total interaction energy is decomposed into the following terms: (1)

∆Eint = E(1)Pol + E(1)ex + E(2)ind + E(2)ex-ind + E(2)disp + E(2)ex-disp + δHF,

where E(1)Pol and E(1)ex are first-order electrostatic and exchange-repulsion terms; E(2)ind and E(2)ex-ind are second-order induction and exchange-induction energies, which can be grouped into the total induction term E(2)Ind; E(2)disp and E(2)ex-disp are second-order dispersion and exchange-dispersion terms, which can be grouped into the total dispersion E(2)Disp and δHF term, which represents higher-order corrections to the electrostatic and induction terms covered by the Hartree-Fock interaction energy but neglected in the second-order SAPT. For more details about this method, see Hesselmann et al.10 The DF-DFT-SAPT calculation performed in combination with a relatively large basis set, e.g. augcc-pVTZ,18,19 provides a reasonably accurate description for various types of noncovalent complexes.20 Even though the scaling of the DF-DFT-SAPT (N5) is not enormous, it is still relatively impractical to apply this method to larger complexes. The bottleneck of the method is the evaluation of the second-order exchange-dispersion term (E(2)ex-disp), which scales with the fifth power of the system size. The E(2)disp itself scales with the fourth power, but the evaluation of the overall dispersion requires the evaluation of both terms. Besides these, the evaluation of the


Page 2 of 23

Page 3 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


response function, which is a prerequisite for the calculation of the dispersion terms, scales also as (N5) in the DF-DFT-SAPT approach. This limitation can be avoided by calculating the whole dispersion component (E(2)disp + E(2)ex-disp) by more efficient means. The dispersion energy can be approximated rather accurately already by very simple pairwise formalism using fixed atomic dispersion coefficients. The cost of such a calculation is negligible compared to the calculation of the remaining DFT-SAPT terms. This approach has been widely used to calculate a dispersion correction for DFT in the DFT-D methodology.21–24 The same approach was adapted for DFT-SAPT by Hesselmann several years ago.25 Nevertheless, there is one important difference. In DFT, only the part of the dispersion energy not covered by the functional is calculated, while the calculation of E(2)Disp for DFT-SAPT requires the calculation of the complete dispersion energy. The formalism used by Hesselmann is, however, rather simple and the parameterization has only been limited to some elements. To exploit the full potential of this approach, improving both its applicability and accuracy, we have developed a new version of an empirical dispersion term for DFT-SAPT. It is based on the successful third-generation dispersion correction for DFT developed by Grimme (referred to as D3 correction).20 The D3 approach has several important advantages that result in higher accuracy. The details are discussed below; in brief, D3 can be viewed as a semiempirical methodology because it is based on precalculated dispersion coefficients and other atomic data and uses only a few adjustable parameters. To parameterize the empirical dispersion, a large set of high-quality DFT-SAPT data is needed. We have calculated the DFT-SAPT decomposition in the aug-cc-pVTZ basis set for the S66x8 data set, a set of dissociation curves of 66 non-covalent complexes of diverse nature.21,22 This is the largest basis set practically applicable to all the systems in the data set; the error caused by its use is discussed for smaller models where DFT-SAPT calculations in a series of basis sets up to aug-ccpV5Z were performed.

Empirical dispersion formalism Before moving to the empirical dispersion for DFT-SAPT, we will summarize the formalism as used in DFT and highlight the points where the calculation of the dispersion term to be used with DFT-SAPT differs. The dispersion is modeled via the atom–atom pairwise empirical potential (optionally, a three-body term can be added in a similar manner). The dispersion is a manifestation of a long-range electron correlation, which is not covered by the standard (LDA, GGA, metaGGA, hybrid) DFT functionals. As a consequence, the DFT provides a poor description for noncovalent complexes, especially those dominated by dispersion. Multiple other strategies to model dispersion in DFT have been developed, including, for example, truly nonlocal dispersion functionals,23–25 empirical re-parameterization of current functionals,31–33 and double–hybrid functionals.26,27 In the next paragraphs, attention will be exclusively paid to the atom–atom pairwise empirical dispersion correction. The first study to describe DFT-based interatomic potential augmented with atom pairwise dispersion correction was carried out by Cohen and Pack.36 Several papers about dispersioncorrected DFT were published by Gianturco et al. in the 1990s.37–40 The first attempts aimed at a more accurate description of noncovalent interactions in bio-organic complexes via DFT augmented ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

by dispersion correction go back to 2001–2002.41–43 A plenty of new, more advanced approaches have appeared since then20,28–30 and have recently been widely used in the field. The general expression for empirical atom–atom pairwise dispersion correction can be written as: Edisp [ρ ]= − ∑

∑

−n

ij

sn C n (ρ )Rij f damp (ρ )

i ≤ j n=6,8,10,. ..

(2)

,

where i and j are atom labels, Cijn(ρ) are density-dependent dispersion coefficients of particular atomic pairs ij, Rij is the inter-atomic distance and fdamp (ρ) is a density-dependent damping function. Optionally, each term in the expansion could be scaled by an empirical factor sn. The asymptotic expansion of Cijn(ρ) coefficients in Eq. 2 can be derived from the second-order perturbation theory expression for dispersion energy between isolated molecules.44 The specific form of Eq. 2 utilized by a vast majority of the DFT-D methods mostly differs in the following: (i) the length of the expansion in Eq. 2 (i.e. n=6,(8,10)), (ii) the type of the damping function utilized, and (iii) the degree of the approximation of the Cijn (ρ) coefficients. The consideration of the dispersion energy term only in the form of the C6/r6 + C8/r8 +... expansion, which is only valid for long intermolecular distances, would cause several problems. Firstly, the singularity problem of this expansion appears for RAB → 0. Secondly, as the DFT covers the local correlation, a part of the dispersion energy originating in the overlap region is included in the DFT energy, but it cannot be separated from it. To avoid the double-counting of the dispersion at these distances (including equilibrium intermolecular distances), the damping of the asymptotic expansion has to be introduced (see Eq. 2). This also solves the singularity problem for the RAB → 0. The use of an appropriate damping function is essential for obtaining accurate DFT-D energies. The most commonly used form of damping function is the Fermi function, whose modified version was used in Jurecka’s30 and Grimme’s D128 and D226 versions of empirical correction (cf. Eq. 3). This type of damping leads to zero contribution to dispersion energy as the intermolecular distance approaches zero: 1

f damp = −α

(3)

1+ exp

r ij −1 vdW vdW r +r i j

(

)

,

where α is an empirical parameter determining the slope of the damping function (e.g. a value of 23 was used in Grimme’s work21) and ri(j)vdW is the atomic van der Waals radius. The rational form of damping, proposed by Becke and Johnson, enters the dispersion energy in a different place (cf. Eq. 4). In the expansion up to n = 8, it reads ij

Edisp = −

(4)

∑ ∑ sn

n=6,8 i ≠ j

Cn

n

Rnij + (a1 Rij0 +a 2 )

,

where s6=1 and s8, a1 and a2 are function-dependent parameters. The R0ij is an empirical fixed cutoff radius for the atom pair ij. The Cnij are dispersion atomic pair parameters and the Rij are inter-atomic ACS Paragon Plus Environment

Page 4 of 23

Page 5 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distances. This type of damping, contrary to the Fermi (zero) type of damping (cf. Eq. 3), leads to a constant contribution to the dispersion energy for each spatially infinitely close atomic pair. In Grimme’s D3 method, it was possible to use both of these damping functions, the latter of which was found to be more accurate in most cases.31 In DFT-D3, the asymptotic expansion is truncated at n = 8, while the scaling of the C8/r8 term is freely optimized in order to account effectively also for higher-order terms. The D3 formalism optionally includes also the leading non-additive term,23 the three-body contribution modeled using the Axilrod-Teller-Muto approximation.46,47 While we do not include it in our method as it is not significant in the studied systems, it can be in principle added to improve the description of the dispersion in larger complexes. However, none of these approaches is suited well for use in DFT-SAPT. The dispersion correction for DFT should mostly account only for long-range correlation effects; in DFT-SAPT, on the other hand, the E(2)disp and E(2)ex-disp terms represent the complete dispersion energy. Hence, when designing the empirical dispersion term for DFT-SAPT, it is evident that the role of the damping is significantly different. In contrast to the DFT, there is no issue with double-counting of dispersion in the middle- and short-range region. On the other hand, the asymptotic expansion still has to be damped in order to account for the exchange and overlap effects at short distances. Such damping of the dispersion in this region does not have to be so extensive and different forms of the damping function are suited better for this purpose. The damping function used in Hesselmann’s original work is represented by an error function

(

f damp =erf a1

(5)

r ij vdW r i +r vdW j

)

,

where a1 is the optimized parameter (a1 = 1.087) and rivdW and rjvdW are the respective atomic van der Waals radii taken from Grimme’s study.34 To illustrate the difference between this damping function and those used in DFT-D, we have plotted the damped dispersion energies in He dimer along with DFT-SAPT dispersion in Figure 1. While the Fermi damping function damps the dispersion to zero for R→ 0, the Becke–Johnson damping provides a finite value. In the case of the damping via the error function, the dispersion energy is substantially larger than in the first two examples but still provides a finite value for small interatomic distances. However, it should be born in mind that at the shortest distances, the perturbation theory itself is no longer accurate.32–34 Another form of damping function suitable for the calculation of the overall dispersion energy was developed earlier by Tang and Toennies.51 This damping function, an incomplete gamma function of the order n+1, has been derived more rigorously as a smooth connection between long-range dispersion and short-range correlation in the H...H system. It takes the form n

f ndamp = 1− exp (− bR ij )∑

(6)

k=0

k

(bRij ) k!

,

where Rij is the interatomic distance and b is the distance scaling factor. Specifically, b is expressed as a linear function of the threshold radius b = a2.R0ij+a3, where a2 and a3 are optimized parameters.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

The empirical dispersion term for DFT-SAPT The empirical dispersion for DFT-SAPT by Hesselmann25 uses only the leading C6 term damped with the error damping function (Eq. 5). The complete formula is thus monomerA monomerB

Edisp = − s 6

∑

∑

i

j

(7)

(

erf a1

C ij6

r ij r vdW +r vdW i j

)( ) a2

r ij

.

The C6ij dispersion coefficients and ri(j)vdW values were adopted from Grimme’s D2 method.35 The parameters were optimized by fitting the total “SAPT+Dispersion” interaction energy towards the CCSD(T) interaction energies extrapolated to the complete basis set (CBS) limit.23 The fitting was performed on the S22+ data set,20 leading to the value of a1 = 1.087. The a2 exponent of the intermolecular distance (rij) was not fixed to 6 but optimized (a2 = 5.67). The optimization of this parameter was introduced as it is well known that a single-term expansion is not accurate enough in the van der Waals region, where the contribution to the dispersion originating from higher expansion terms is not negligible.36 Moreover, the scaling parameter s6, in analogy with Grimme’s D2 method, is introduced (s6 = 0.67) what leads to underestimated dispersion at larger distances. It is obvious that this formulation of the dispersion energy is rather empirical. In addition, the method is applicable only to a limited set of elements for which the C6 coefficients are available. Some of these issues can be addressed by starting from the D3 dispersion correction. Here, the C6 coefficients and threshold radii are available for a major part of the periodic table. Furthermore, to improve the accuracy, different valence states of each atom are considered as their C6 coefficients may differ. These data were precalculated for all pairs of atoms and not constructed from atomic coefficients. Finally, the inclusion of the higher-order C8 term should also improve the description of dispersion at short distances. The associated scaling factor s8 is made adjustable while s6 is set to 1 in order to conserve physical behavior at long distances. In the D3 correction, the C8 coefficients are derived from C6 using tabulated multipole expected-value ratios (see Eqs. 6–8 in Grimme’s D3 work23). This leads to the formula monomerA monomerB

Edisp = −

∑

∑

i

j

(8)

( f

6 damp

Cij6 r 6ij

+s8 f

8 damp

Cij8 r 8ij

)

.

In Eq. 8, two types of damping functions were tested. The first of them was, in analogy with Hesselmann’s original paper, the damping function based on an error function in the following form:

(

f ndamp =erf a1

(9)

r ij r vdW +r vdW i j

αn

)

.

Herein, an additional parameter α has been introduced. It affects the slope of the damping function. Furthermore, individual damping functions for the C6 and C8 terms have been utilized (cf. Eq. 9). The three options differing in the choice of the α6 and α8 are: a) the α8 parameter is not optimized but set to α8 = α6; b) the α8 parameter is not optimized but set to α8 = α6 + 2; c) the α8 parameter is ACS Paragon Plus Environment

Page 7 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


optimized independently. The sum of the van der Waals radii has been replaced by the Rij0 pairwise threshold radius from the D3 dispersion correction. Secondly, we used the Tang-Toennies (TT) damping function (cf. Eq. 6) as its theoretical background suggests that it is well suitable for this application.37 During the parameterization, we encountered an issue in the D3 dispersion that led to serious problems with the reproduction of the DFT-SAPT dispersion. In the original formulation, the C6 coefficients are geometry-dependent because they are interpolated from the tabulated data on the basis of a continuous valence number determined from the geometry (see Eqs. 15 and 16 in Grimme’s D3 work23). This valence number for a given atom thus differs slightly in the dimer and in the monomer when the intermolecular dispersion energy is calculated using the supramolecular approach (which we have adopted because of the simplicity of its implementation in an existing code which yields total dispersion energy). The resulting difference in the C6 coefficient is negligible, but it propagates into the C8 coefficient, which is derived from it (Eq. 6 in the Grimme’s D3 work23). Because of the steep slope of the C8/r8 term at short distances, even this small difference may lead to a large change in the interaction dispersion energy. This effect is hidden by the damping in DFT-D3 but becomes large when no or less aggressive damping is used, such as in this case. We have fixed this issue by using integer valence numbers which are the same in the dimer and in the monomers. This approximation is perfectly valid in the context of DFT-SAPT calculation where only the intermolecular dispersion between well-defined molecular species is calculated (i.e. no covalent bonds are being broken). In such cases, rounding the valence numbers to integers does not affect the accuracy of the C6 coefficients while it removes the unwanted artifacts. The calculation of the dispersion energy with both damping functions was implemented in the Cuby framework;38 practical information on how to perform the calculations is provided in the documentation.39

DFT-SAPT calculations The DF-DFT-SAPT reference calculations were performed in the Molpro2015 package.40 The monomer calculations were carried out in the dimer basis set. The asymptotically-corrected PBE0AC functional in combination with the aug-cc-pVTZ basis set was used.18,19,57,58 The asymptotic correction was calculated using PBE0 in the aug-cc-pVTZ basis set.41 All the individual DFT-SAPT terms as well as the overall DFT-SAPT and CCSD(T)/CBS interaction energies calculated for the S66x8 data set are provided in the Supplementary Information. We believe that these data are a valuable extension of the data set that can be used for other applications as well. To test the accuracy of the above mentioned reference data and the fitted dispersion with respect to basis set size, we use the A24 data set60 of small complexes for which we calculated the DFT-SAPT interaction energies in a series of basis sets up to aug-cc-pV5Z. The computational setup was the same as used for the S66x8 data set. The DFT-SAPT interaction energy at the complete basis set limit was then estimated as the sum of E(1)Pol, E(1)ex, E(2)ind, E(2)ex-ind and δHF terms calculated in aug-cc-pV5Z basis and second-order dispersion terms (E(2)disp and E(2)ex-disp) extrapolated using the Helgaker scheme61,62 from aug-cc-pVQZ and aug-cc-PV5Z.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Fitting procedure Unlike Hesselmann, we parameterize the empirical dispersion term solely to its DFT-SAPT counterpart, the sum of the DFT-SAPT second-order dispersion terms. Our aim is to provide a dispersion term fully compatible with the DFT-SAPT calculation rather than reproducing the most accurate interaction energy. The S66x8 data set26,27 was used as the training set for the parameterization. This data set consists of potential energy curves for 66 noncovalent complexes, which can be divided into three categories: hydrogen-bonded, dispersion-dominated complexes and mixed ones. Each potential curve is represented by 8 structures with the following multiplication of the equilibrium intermolecular distance: 0.9, 0.95, 1.0, 1.05, 1.10, 1.25, 1.50 and 2.0. For more details about this dataset, see Řezáč et al.22 The coverage of the non-equilibrium distance is of key importance as it is well known that the reliability of the perturbative DFT-SAPT method for really short intermolecular distances is relatively poor. This is a problem of the perturbative expansion, which becomes illdefined when the perturbation becomes large. Hence, we repeated the parameterization also with the data points (i.e. structures) with the intermonomer distance smaller than the equilibrium removed. Henceforth, the optimization towards the whole data set (i.e. S66x8) and its non-repulsive (truncated) part will be noted as “S66x8” and “S66x6”, respectively. The following combinations of parameters were optimized. In the α8=α6 variant (cf. Eqs. 8 and 9) of the error function, these were: a1; a1 and α6; a1 and s8; a1, α6 and s8; in the case of the α8=α6 + 2 variant (cf. Eqs. 8 and 9): a1; a1 and α6; a1 and s8; a1, α6 and s8; a1, α6 and α8; a1, α6, α8 and s8. The Tang–Toennies damping function was tested for the following combinations of parameters: a2 and a3; a2, a3 and s8. The starting values of all the considered parameters before the fitting procedure were the following: a1=1.0, α6=6.0, α8=8.0 and s8=1.0, and in the case of the Tang–Toennies damping function, they were: a2=-0.2, a3=4.0, s8=1.0. Gradient optimization based on the numerically calculated gradient of the parameters was used to minimize the root-mean-square-error (RMSE) in the training set. The convergence threshold for the optimization procedure was set to 0.0001 kcal.mol-1 for all the fitting calculations.

Validation data sets We utilized several data sets for the testing of the developed methods. In all of these validation data sets, the sum of the second-order dispersion terms, calculated at the (DF-)DFT-SAPT (PBE0AC/aug-cc-pVTZ) level, was used as a reference. The charge-transfer (CT) data set. A set of 11 complexes with weak, medium and strong chargetransfer (CT) character.63 These systems make tests difficult because of their very short equilibrium intermolecular distances. The structures of the complexes as well as the DFT-SAPT results have been taken from the original reference without any modifications. The halogen-bonded (X40-XB) data set. To test the transferability of the parameters to different elements, a subset of 18 complexes from the original X40 halogen bond data set was chosen.42 These complexes contain a direct Cl, Br, I….N, O contact which is involved in the interaction.


Page 8 of 23

Page 9 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The mixed (11x6-Mix) data set. A set of potential curves for 11 noncovalent complexes was constructed, including the C2H2 and C2H4 dimers in the stacking conformation as representatives of the π...π interaction; the C2H2 and C2H4 dimers in the T-shape conformation as models of the CH...π interaction; CH3Br...OCH2 and CHF3...CH4 complexes as halogen-bond representatives; aliphatic interaction was captured by the CH4 dimer; the CH2O dimer, H2O...CH3OH, H2O dimer and H2O...NH2CH3 were all representatives of hydrogen-bonded complexes. The xyz-coordinates for all of these complexes can be found in the Supplementary Information. The potential curves are represented by 6 structures, with the following multiplication of the main intermolecular coordinate: 0.7, 0.8, 0.9, 1.0, 1.25 and 1.50. The definition of the main intermolecular coordinate depends on the type of noncovalent interaction, e.g. for hydrogen bonds it is defined by the donor and acceptor atoms, in most other cases it is defined by the centers of mass of the monomers. (For details see Supporting Information of ref. 26). The empirical dispersion term constructed in the above-mentioned manner can be used not only within the DFT-SAPT framework, but also for the quantification of the dispersion part of the total interaction energy. The simplified version of such an empirical dispersion model has already been used for the investigation of bonding in Lappert’s stannylene complexes and their analogs.65

Results and discussion The RMSE and relative RMSE (rRMSE, calculated as RMSE divided by the average magnitude of the reference value, in percent) of the three best- and worst-performing methods for the validation sets are listed in Table 1. The results for all the twelve variants of the damping function are provided in the Supplementary Information. The optimized values of the fitted parameters for all the 12 tested combinations of parameters together with a brief discussion can be found in the Supplementary Information. Here, we will discuss only the three variants best performing in the validation data sets. We list the parameters optimized for the complete S66x8 as well as the truncated S66x6 set (in parentheses). The former are more general but the latter make it possible to achieve better accuracy in systems where close contacts are not present. It is clear that the Tang–Toennies damping function works better than the error-function damping. The two best variants use the TT damping with two and three optimized parameters: a2, a3 and a2, a3, s8, respectively. The corresponding optimal values of fitted parameters are -0.752 (-0.825), 5.357 (5.576) and -0.436 (-0.365), 4.757 (4.656), 0.869 (0.882). The third best performing variant of the empirical term utilizes the error function as damping with only one fitted parameter, a1, α6 and α8 have been set to 6. The optimized values of a1 parameter are 1.234 (1.241). The RMSE and rRMSE values for model as well as validation data sets are noticeably lower for the S66x6 optimization (cf. Table 1, part B) than for the S66x8 optimization (cf. Table 1, part A). In the case of the model data set, the average RMSE and rRMSE (listed in parentheses) for all the 12 tested variants of empirical dispersion are 0.185 (7.2) as compared to 0.287 (7.3) kcal.mol-1 (%) for the S66x8 parameterization (cf. Table S2). The same trend holds for each validation data set, but not so significantly. The respective average RMSE and rRMSE (listed in parentheses) values for the CT, X40-XB and 11x6-Mix data sets are 2.956 (29.1), 0.465 (12.8), 0.996 (13.5) and 3.090 (32.1),



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.477 (13.1), 1.044 (14.3) kcal.mol-1 (%) for the S66x6 and S66x8 optimization, respectively (cf. Table S2). The slightly larger errors in the case of S66x8 parameterization have two sources that cannot be separated. Firstly, the coverage of the larger interval of distances places higher demands on the damping function and the empirical dispersion formalism in general (e.g. the use of higher-order terms and their accuracy). Secondly, it is well known that the accuracy of the perturbative DFTSAPT method itself deteriorates in the region of small intermonomer distances (both due to the limitation of the truncated expansion in the perturbation theory as well as due to an approximation used in the evaluation of the exchange-dispersion term, which can be avoided).66 Thus, fitting a model with limited flexibility (with fixed form of the damping function and only few parameters) to these points may compromise the accuracy even at larger distances. First, we discuss the performance of empirical dispersion optimized towards the whole S66x8 set. Among all the variants tested, the best performance is exhibited by the one that utilizes TT function for damping with three optimized parameters (a2, a3 and s8), with an average RMSE (rRMSE) of 0.727 kcal.mol-1 (15.2 %). The second best variant also uses the TT damping function with two optimized parameters and it yields a RMSE (rRMSE) of 0.971 kcal.mol-1 (15.8 %). The third one utilizes error function for damping with fixed exponents α6 = α8 = 6.000, with the only optimized parameter being a1. It yields a substantially larger RMSE of 1.450 kcal.mol-1 (19.3 %). Surprisingly, the same model where an additional parameter (s8) is optimized performs the worst, with the RMSE (rRMSE) of 3.142 kcal.mol-1 (35.8 %). The situation is analogical for models optimized on the S66x6 data set. Here again, the performance of two forms utilizing the TT damping function is superior. The respective RMSE (rRMSE) values are 0.675 kcal.mol-1 (11.4 %) (optimized parameters: a2, a3 and s8) and 0.956 kcal.mol-1 (11.0 %) (optimized parameters: a2 and a3), followed by a model employing an error function for damping with the absolute (relative) error of 1.415 kcal.mol-1 (19.0 %) (optimized parameter: a1, α6 = α8 = 6.000). The improvement over the S66x8 parameterization is small (only 0.05 kcal.mol-1); for practical use, we recommend the variant optimized on the S66x8 set because it is more general. Closer inspection of Table 1 and Tables S2 and S3 in the Supplementary Information reveals that the largest errors are observed for CT validation data sets. This is understandable as the CT complexes are very often characterized by their close inter-monomer distance (for more details, see the discussion in previous sections). The exclusion of the CT data set from statistics brings a significant reduction of errors. Specifically, the RMSE (rRMSE) of the originally best threeparameter TT variant is decreased from 0.727 to 0.458 kcal.mol-1 (15.2 → 10.3 %) and from 0.675 to 0.567 kcal.mol-1 (11.4 → 11.1 %) for S66x8 and S66x6 optimization. A similar reduction of errors is also observed in the other variants, most notably in the case of error-function damping. This only confirms that the error function for damping is not as good as the TT function, especially in the case of short intermolecular distances. To test the quality of the fit across various distances, we calculated the errors of the most accurate variant of the empirical dispersion (which employs TT damping function with 3 optimized parameters) separately for each of the relative distances used in the S66x8 curves (defined by factors scaling the closest intramolecular distance in equilibroium). The Table 2 list the absolute


Page 10 of 23

Page 11 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


RMSE (2nd column) and relative RMSE (rRMSE, 3rd column) values. From these data, it is evident that absolute error is systematically decreased with increasing intermolecular distance from 0.41 kcal.mol-1 to 0.086 kcal.mol-1 which is expected as the absolute value of the dispersion energy also decreases. The relative errors are a more important measure as they are independent of the magnitude of the dispersion energy. They vary only negligibly between 6.5 and 8.6 % which indicates that the form of the damping function is appropriate and the fit is good. Finally, we compare our method to Hesselmann’s original empirical dispersion25 as implemented in MOLPRO.56 It should be, however, kept in mind that such a comparison is not fully appropriate, as Hesselmann’s empirical term was fitted to reproduce overall interaction energies, whereas we fit our dispersion only to its DFT-SAPT counterpart. Nevertheless, this difference in the reference data is smaller than the errors observed in the results. We compare these methods for the X40-XB data set of halogen bonds, for which the results of Hesselmann’s method have been taken from Kolář et al.67 This is a good test of the transferability of the method as both approaches have been fitted on a set of organic molecules not containing halogens. The RMSE (rRMSE) of our best-performing method is 0.401 kcal.mol-1, while the original approach yields an error of 0.605 kcal.mol-1. To illustrate the computational efficiency of the empirical dispersion, we look at the largest system in the S66x8 set, the stacked uracil dimer in equilibrium geometry (run at 8-core Xeon server with 48 GB RAM). Here, about half of the CPU time is spend in DFT-SAPT (the rest are the calculations of the monomers and the evaluation of the δHF term). Out of that, 82% of the time is spent in the evaluation of the dispersion and exchange-dispersion terms which can be replaced by the empirical dispersion calculated at practically no cost. This ratio will be even more favorable in larger systems because the complexity of the calculation of the dispersion terms has a steeper scaling than the rest of the calculation.

Accuracy with respect to basis set size The largest basis set applicable consistently to all the systems in the S66x8 data set, aug-cc-pVTZ, was used as a reference for fitting the dispersion term. It is, however, well known that dispersion in this basis set is somewhat underestimated. Hence, the extrapolation of the dispersion energy to the complete basis set limit is commonly used.68–70 Here, this error is discussed in comparison with both DFT-SAPT and CCSD(T) interaction energies extrapolated to the CBS limit. For this purpose, we use the A24 data set of small complexes60 for which accurate CCSD(T) results (extrapolated from a series of basis sets up to aug-cc-pV5Z) are available. DFT-SAPT/CBS interaction energies in this set were calculated for this work. The results are summarized in Table 3. To put the accuracy of DFT-SAPT with D3 dispersion into proper context, we list also the errors of complete DFT-SAPT interaction energies calculated in different basis sets as well as estimates of DFT-SAPT/CBS limit obtained by both extrapolation61,62 or scaling of the dispersion terms.20 The error with respect to DFT-SAPT/CBS (2nd column in Table 3) represents the error due to the incompleteness of the basis set, whereas the error calculated against the CCSD(T)/CBS reference (3rd column) includes also the inherent error of the DFT-SAPT method itself. The error of DFT-SAPT in finite basis sets with respect to DFT-SAPT/CBS reference illustrates the smooth but rather slow convergence towards the CBS limit. Here, the calculations in the aTZ basis ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(which was used for the S66x8 data set) yield a RMSE of 0.12 kcal/mol. When these results are compared to the CCSD(T)/CBS reference, the RMSE of DFT-SAPT/aug-cc-pVTZ is twice as large, 0.24 kcal/mol. It is an improvement compared to aDZ calculations but the error stays about the same (and increases slightly) in larger basis sets. It is thus evident that, except for the aug-cc-pVDZ basis set, the error originating in the method is larger than the error due to the basis set incompleteness. Using DFT-SAPT/aug-cc-pVTZ for fitting the empirical dispersion is thus a safe choice, using larger basis set would not improve the absolute accuracy of the method. The second part of Table 3 lists the RMSE for estimated DFT-SAPT/CBS interaction energies, which were obtained by extrapolation or scaling of the dispersion energy terms. The scaling factors were taken from ref. 20; our present calculations in the A24 set confirm their validity. The extrapolation is more accurate than the simple scaling. However, the scaled aug-cc-pVTZ calculations still yield an error smaller than unscaled calculations in aug-cc-pVQZ what is a remarkable result given the simplicity of the method. Again, when these results are compared to CCSD(T)/CBS, the limitation is the accuracy of DFT-SAPT itself. Finally, DFT-SAPT (in aug-cc-pVTZ) with the D3 dispersion yields errors of about 0.3 kcal/mol with respect to both DFT/SAPT and CCSD(T) reference (see the bottom of Table 3). The error with respect to DFT-SAPT/CBS is larger than the error of the aug-cc-pVTZ reference which is caused by the limited accuracy of the empirical dispersion formalism and parameters. On the other hand, the error with respect to CCSD(T)/CBS is not much larger than the error of complete DFT-SAPT calculations, including those in very large basis set. We have tried scaling the empirical dispersion using the scaling factor derived for the aug-cc-pVTZ it was fitted to but the improvement is not significant. Analogously, fitting the D3 dispersion to scaled DFT-SAPT dispersion did not bring any improvement.

Conclusions We have developed an enhanced version of the empirical dispersion term for the DFT-SAPT framework based on Grimme’s D3 dispersion correction for DFT. Two different damping functions: A/ the Tang–Toennies (TT) function and B/ the error function, with different numbers and combinations of freely optimized parameters were tested. The fitting was done towards the robust S66x8 data set in two variants: A/ the whole S66x8 dataset and B/ S66x6 (a subset of S66x8 where points with inter-monomer distances smaller than equilibrium were excluded). Overall, the errors of empirical dispersion terms that employ the TT function for damping are substantially smaller than those with the error function. Moreover, the stronger physical background of the TT provides another argument supporting its use. Among the multiple parameterizations we have tested, the three-parameter TT version parameterized to the whole S66x8 data set is the most general, yet very accurate with a RMSE (rRMSE) of 0.727 kcal.mol-1 (15.2 %), and it can be recommended for general use. The same approach parameterized to the S66x6 set performs marginally better, but the improvement is too small to justify. The exclusion of problematic, very compact CT complexes from validation data sets generally brings a significant reduction of errors. The two-parameter TT parameterization towards the S66x8 data set performs overall best with a RMSE (rRMSE) of only 0.419 kcal.mol-1 (10.0 %).


Page 12 of 23

Page 13 of 23

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We believe that the proposed empirical dispersion terms can be used in the DFT-SAPT framework as a reasonably accurate substitution for the original dispersion terms in the cases where the full DFT-SAPT calculation would be too demanding. It is also applicable outside the DFT-SAPT framework, e.g. for the estimation of the amount of dispersion in a calculation where only the total interaction energy is known. The implementation of the method is readily available in the Cuby framework.38 The DF-DFT-SAPT decomposition for the whole S66x8 data set reported in this work is a valuable extension of this data set, which may be useful in other applications.

Acknowledgements This work has been supported by the Czech Science Foundation by grant No. P208/16-11321Y and is a part of the Research Project RVO 61388963 of the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic.

Supporting Information The optimized values of the fitted parameters for all the 12 tested combinations of parameters withing the empirical dispersion models are presented. The RMSE and relative RMSE (rRMSE, calculated as RMSE divided by the average magnitude of the reference value, in percent) results for all the twelve variants of the damping function are provided for all considered validation data sets. All the individual DFT-SAPT terms as well as the overall DFT-SAPT and CCSD(T)/CBS interaction energies calculated for the S66x8 data set are provided.



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

Page 14 of 23

Table 1. The RMSE (in kcal.mol-1) and relative RMSE (rRMSE in %, listed in parentheses) values of the three best- and worst-performing variants of the empirical dispersion term with respect to the reference (DF-)DFT-SAPT dispersion data. The RMSE and rRMSE values are listed for model as well as all three validation data sets. Part A lists data corresponding to the fitting that was performed for the whole S66x8 data set, whereas Part B only for its non-repulsive part. A/ Type of damping function TTa TTa erfb

Optimized parameters a2, a3, s8 a2, a3 a1

RMSE (rRMSE) model CT 0.259 (7.3) 1.266 (25.0) 0.265 (7.3) 2.076 (27.3) 0.302 (7.3) 3.065 (32.3)

X40-XB 0.401 (10.8) 0.385 (10.6) 0.505 (14.1)

11x6-Mix 0.515 (9.7) 0.453 (9.5) 0.781 (11.4)

aver 0.727 (15.2) 0.971 (15.8) 1.450 (19.3)

avernoCTc 0.458 (10.3) 0.419 (10.0) 0.643 (12.8)

erfb

a1, s8

0.291 (7.5)

4.167 (39.1)

0.664 (18.4)

4.595 (50.0)

3.142 (35.8)

2.630 (34.2)

Type of damping function TTa TTa erfb

Optimized parameters a2, a3, s8 a2, a3 a1

RMSE (rRMSE) model CT 0.173 (7.3) 0.893 (12.1) 0.177 (7.3) 1.986 (12.3) 0.189 (7.3) 2.967 (31.7)

X40-XB 0.415 (10.6) 0.378 (10.4) 0.489 (13.7)

11x6-Mix 0.718 (11.5) 0.503 (10.3) 0.790 (11.6)

aver 0.675 (11.4) 0.956 (11.0) 1.415 (19.0)

avernoCTc 0.567 (11.1) 0.441 (10.3) 0.640 (12.7)

erfb

a1, s8

0.187 (7.4)

0.472 (13.2)

4.003 (40.8)

2.516 (28.7)

2.238 (27.0)

B/

a b c

3.074 (32.0)

stands for the Tang–Toennies damping function (cf. Eq. 6) stands for the error function as a damping function with the α8 parameter set up as: α8=α6 (cf. Eq. 9) This average value of RMSE (rRMSE) is calculated as follows: (RMSEX40-XB + RMSE11x6-Mix)/2, i.e. the CT data set has been excluded.


Page 15 of 23

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


Table 2. The absolute RMSE (in kcal.mol-1) and relative RMSE (rRMSE, in %) values of empirical dispersion correction for different intermolecular distance. Values are listed for subsets of structures from the S66x8 data set, belonging to specific multiplication of main intermolecular coordinate. Presented numbers correspond to the most accurate empirical dispersion employing TT damping function (three parameter version) with respect to the reference DFT-SAPT values. Distancea

RMSE

rRMSE

0.9 0.95 1 1.05 1.1 1.25 1.5 2

0.409 0.349 0.257 0.232 0.200 0.141 0.086 0.122

6.6 6.5 6.6 6.7 6.9 7.6 8.3 8.6

a

values represent multiplication of the equilibrium value of the main intermolecular coordinate



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

Page 16 of 23

Table 3. The RMSE (in kcal.mol-1) values of DFT-SAPT in different basis sets with respect to the DFT-SAPT/CBS and CCSD(T)/CBS reference in the A24 data set. Method

RMSEDFT-SAPT/CBS

RMSECCSD(T)/CBS

Finite basis aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVQZ aug-cc-pV5Z

0.389 0.118 0.045 0.023

0.363 0.239 0.250 0.259

Extrapolation of dispersion aug-cc-pV(D→T)Z 0.022 aug-cc-pV(T→Q)Z 0.007 aug-cc-pV(Q→5)Z 0.000

0.264 0.271 0.269

Scaled dispersion aug-cc-pVDZ*1.19 aug-cc-pVTZ*1.05

0.117 0.038

0.305 0.272

Empirical dispersion aug-cc-pVTZ + D3 aug-cc-pVTZ + D3*1.05

0.339 0.314

0.313 0.329


Page 17 of 23

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


Figure 1. The dispersion energy evaluated at five different levels of theory for the He...He system. The EE2dispSAPT energy (purple) corresponds to the SAPT second-order dispersion term;71 the -EE2ex-dispSAPT (green) corresponds to the the SAPT second-order exchange dispersion term;71 the EempirDFTSAPT (light blue) corresponds to the Hesselmann-like SAPT empirical dispersion with the error function as a damping function; the EzeroG3/B3-LYP energies (yellow) to Grimme’s D3 dispersion with zero damping, and the EzeroG3/B3-LYP energies (dark blue) to the D3 dispersion with the Becke– Johnson damping. The vertical dotted black line represents the equilibrium distance for the He dimer (req=2.96 Å).



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

Page 18 of 23

References (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16)

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. te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931–967. Su, P.; Li, H. Energy decomposition analysis of covalent bonds and intermolecular interactions. J. Chem. Phys. 2009, 131, 014102. Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment molecular orbital method: application to polypeptides. Chem. Phys. Lett. 2000, 318, 614–618. Horn, P. R.; Mao, Y.; Head-Gordon, M. Probing non-covalent interactions with a second generation energy decomposition analysis using absolutely localized molecular orbitals. Phys Chem Chem Phys 2016, 18, 23067–23079. 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. Heßelmann, A.; Jansen, G. First-order intermolecular interaction energies from Kohn–Sham orbitals. Chem. Phys. Lett. 2002, 357, 464–470. 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. Heßelmann, A.; Jansen, G. Intermolecular dispersion energies from time-dependent density functional theory. Chem. Phys. Lett. 2003, 367, 778– 784. Heßelmann, A.; Jansen, G.; Schütz, M. Density-functional theory-symmetry-adapted intermolecular perturbation theory with density fitting: A new efficient method to study intermolecular interaction energies. J. Chem. Phys. 2005, 122, 014103. 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. Misquitta, A. J.; Szalewicz, K. Intermolecular forces from asymptotically corrected density functional description of monomers. Chem. Phys. Lett. 2002, 357, 301–306. Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Dispersion Energy from Density-Functional Theory Description of Monomers. Phys. Rev. Lett. 2003, 91, 033201. Misquitta, A. J.; Szalewicz, K. Symmetry-adapted perturbation-theory calculations of intermolecular forces employing density-functional description of monomers. J. Chem. Phys. 2005, 122, 214109. Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. Intermolecular potentials based on symmetry-adapted perturbation theory with dispersion energies from time-dependent density-functional calculations. J. Chem. Phys. 2005, 123, 214103. Bukowski, R.; Podeszwa, R.; Szalewicz, K. Efficient calculation of coupled Kohn–Sham dynamic susceptibility functions and dispersion energies with density fitting. Chem. Phys. Lett. 2005, 414, 111–116.


Page 19 of 23

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


(17) 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. (18) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions aug-cc-pVnZ_1992.pdf http://www.chem.ccu.edu.tw/~hu/Web_Lib/literature/quantum/aug-cc-pVnZ_1992.pdf (accessed Dec 6, 2016). (19) Woon, D. E.; Jr, T. H. D. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358–1371. (20) Řezáč, J.; Hobza, P. Extrapolation and Scaling of the DFT-SAPT Interaction Energies toward the Basis Set Limit. J. Chem. Theory Comput. 2011, 7, 685–689. (21) Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463–1473. (22) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction - Grimme - 2006 - Journal of Computational Chemistry - Wiley Online Library http://onlinelibrary.wiley.com/doi/10.1002/jcc.20495/abstract (accessed Dec 6, 2016). (23) 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. (24) Jurečka, P.; Černý, J.; Hobza, P.; Salahub, D. R. 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. (25) Heßelmann, A. Comparison of Intermolecular Interaction Energies from SAPT and DFT Including Empirical Dispersion Contributions. J. Phys. Chem. A 2011, 115, 11321–11330. (26) Řezáč, 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. (27) Řezáč, 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. (28) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92. (29) Sato, T.; Tsuneda, T.; Hirao *, K. Van der Waals interactions studied by density functional theory. Mol. Phys. 2005, 103, 1151–1164. (30) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103. (31) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157–167. (32) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Exchange-correlation functional with broad accuracy for metallic and nonmetallic compounds, kinetics, and noncovalent interactions. J. Chem. Phys. 2005, 123, 161103. (33) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364–382. (34) Grimme, S. Semiempirical hybrid density functional with perturbative second-order correlation. J. Chem. Phys. 2006, 124, 034108.



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

Page 20 of 23

(35) Schwabe, T.; Grimme, S. Double-hybrid density functionals with long-range dispersion corrections: higher accuracy and extended applicability. Phys. Chem. Chem. Phys. 2007, 9, 3397. (36) Cohen, J. S.; Pack, R. T. Modified statistical method for intermolecular potentials. Combining rules for higher van der Waals coefficients. J. Chem. Phys. 1974, 61, 2372–2382. (37) Gianturco, F. A.; Paesani, F. The rovibrational structure of the Ar–CO complex from a model interaction potential. J. Chem. Phys. 2001, 115, 249–256. (38) GIANTURCO, F. A.; PAESANI, F.; LARANJEIRA, M. F.; VASSILENKO, V.; CUNHA, M. A.; SHASHKOV, A. G.; ZOLOTOUKHINA, A. F. Computed and measured thermal diffusion factor for CO-He mixtures: a test of recent interaction potentials. Mol. Phys. 1997, 92, 957–972. (39) GIANTURCO, F. A.; PAESANI, F.; LARANJEIRA, M. F.; VASSILENKO, V.; CUNHA, M. A.; SHASHKOV, A. G.; ZOLOTOUKHINA, A. F. Computed and measured transport coefficients for CO-He mixtures: testing a density functional approach. Mol. Phys. 1998, 94, 605–622. (40) Gianturco, F. A.; Paesani, F.; Laranjeira, M. F.; Vassilenko, V.; Cunha, M. A. Intermolecular forces from density functional theory. III. A multiproperty analysis for the Ar(1S)-CO(1Σ) interaction. J. Chem. Phys. 1999, 110, 7832–7845. (41) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen bonding and stacking interactions of nucleic acid base pairs: A densityfunctional-theory based treatment. J. Chem. Phys. 2001, 114, 5149–5155. (42) Wu, Q.; Yang, W. Empirical correction to density functional theory for van der Waals interactions. J. Chem. Phys. 2002, 116, 515–524. (43) Wu, X.; Vargas, M. C.; Nayak, S.; Lotrich, V.; Scoles, G. Towards extending the applicability of density functional theory to weakly bound systems. J. Chem. Phys. 2001, 115, 8748–8757. (44) Kaplan, I. G. Intermolecular Interactions: Physical Picture, Computational Methods and Model Potentials; John Wiley & Sons, Ltd: Chichester, UK, 2006. (45) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (46) Axilrod, M.; Teller, E. Interaction of the van der Waals Type Between Three Atoms. J. Chem. Phys. 1943, 11, 299–300. (47) Muto, Y. Force between nonpolar molecules. Proc. Phys.-Math. Soc. Jpn. 1943, 17, 629. (48) Jeziorski, B.; Szalewicz, K.; Chałasiński, G. Symmetry forcing and convergence properties of perturbation expansions for molecular interaction energies. Int. J. Quantum Chem. 1978, 14, 271–287. (49) Korona, T.; Moszynski, R.; Jeziorski, B. Convergence of Symmetry-Adapted Perturbation Theory for the Interaction between Helium Atoms and between a Hydrogen Molecule and a Helium Atom. In Advances in Quantum Chemistry; Per-Olov Löwdin, J. R. S., Michael C.Zerner, Jacek Karwowski and Mati Karelson, Ed.; Academic Press, 1997; Vol. 28, pp 171–188. (50) Patkowski, K.; Jeziorski, B.; Szalewicz, K. Symmetry-adapted perturbation theory with regularized Coulomb potential. J. Mol. Struct. THEOCHEM 2001, 547, 293–307. (51) Tang, K. T.; Toennies, J. P. An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients. J. Chem. Phys. 1984, 80, 3726.


Page 21 of 23

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


(52) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787– 1799. (53) Adamovic, I.; Gordon, M. S. Dynamic polarizability, dispersion coefficient C 6 and dispersion energy in the effective fragment potential method. Mol. Phys. 2005, 103, 379–387. (54) Řezáč, J. Cuby: An integrative framework for computational chemistry. J. Comput. Chem. 2016, 37, 1230–1237. (55) Cuby 4 http://cuby4.molecular.cz/ (accessed Dec 6, 2016). (56) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Györffy, W.; Kats, D.; Korona, T.; Lindh, R.; et al. MOLPRO, version 2015.1, a package of ab initio programs; molpro, 2015. (57) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. (58) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys. 1996, 105, 9982–9985. (59) Grüning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. Shape corrections to exchange-correlation potentials by gradientregulated seamless connection of model potentials for inner and outer region. J. Chem. Phys. 2001, 114, 652. (60) Řezáč, 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. (61) 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. (62) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Olsen, J. Basis-set convergence of the energy in molecular Hartree–Fock calculations. Chem. Phys. Lett. 1999, 302, 437–446. (63) Karthikeyan, S.; Sedlá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. (64) Řezáč, J.; Riley, K. E.; Hobza, P. Benchmark Calculations of Noncovalent Interactions of Halogenated Molecules. J. Chem. Theory Comput. 2012, 8, 4285–4292. (65) Sedlak, R.; Stasyuk, O. A.; Fonseca Guerra, C.; Řezáč, J.; Růžička, A.; Hobza, P. New Insight into the Nature of Bonding in the Dimers of Lappert’s Stannylene and Its Ge Analogs: A Quantum Mechanical Study. J. Chem. Theory Comput. 2016, 12, 1696–1704. (66) Schäffer, R.; Jansen, G. Single-determinant-based symmetry-adapted perturbation theory without single-exchange approximation. Mol. Phys. 2013, 111, 2570–2584. (67) Kolář, M. H.; Deepa, P.; Ajani, H.; Pecina, A.; Hobza, P. Characteristics of a σ-Hole and the Nature of a Halogen Bond. In Halogen Bonding II; Metrangolo, P., Resnati, G., Eds.; Springer International Publishing: Cham, 2014; Vol. 359, pp 1–25. (68) Leforestier, C.; Tekin, A.; Jansen, G.; Herman, M. First principles potential for the acetylene dimer and refinement by fitting to experiments. J. Chem. Phys. 2011, 135, 234306.



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

Page 22 of 23

(69) 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. (70) Jansen, G. Symmetry-adapted perturbation theory based on density functional theory for noncovalent interactions. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 127–144. (71) Korona, T.; Williams, H. L.; Bukowski, R.; Jeziorski, B.; Szalewicz, K. Helium dimer potential from symmetry-adapted perturbation theory calculations using large Gaussian geminal and orbital basis sets. J. Chem. Phys. 1997, 106, 5109–5122.


Page 23 of 23

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


For Table of Contents Only


Empirical D3 Dispersion as a Replacement for ab Initio Dispersion

Recommend Documents