Benchmarking the Effective Fragment Potential ... - ACS Publications

Mar 30, 2018 - Department of Chemistry, University of Colorado Denver, Denver, Colorado 80217, United States. ‡. Department of Chemistry and Ames La...
0 downloads 5 Views 600KB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

A: Molecular Structure, Quantum Chemistry, and General Theory

Benchmarking the Effective Fragment Potential Dispersion Correction on the S22 Test Set Shinae Kim, Chelsea M Kaliszewski, Emilie B Guidez, and Mark S. Gordon J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b11628 • Publication Date (Web): 30 Mar 2018 Downloaded from http://pubs.acs.org on March 30, 2018

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

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

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

The Journal of Physical Chemistry

Benchmarking the Effective Fragment Potential Dispersion Correction on the S22 Test Set Shinae Kim, Chelsea M. Kaliszewski, Emilie B. Guidez† and Mark S. Gordon* *Department of Chemistry and Ames Laboratory USDOE, Iowa State University, Ames, Iowa 50011, United States. Email: [email protected] †Chemistry Department, University of Colorado Denver, Denver, Colorado 80217, United States Email: [email protected] Abstract The usual modeling of dispersion interactions in density functional theory (DFT) is often limited by the use of empirically fitted parameters. In this study, the accuracies of the popular empirical dispersion corrections and the first principles derived effective fragment potential (EFP) dispersion correction are compared by computing the DFT-D and HF-D equilibrium interaction energies and intermolecular distances of the S22 test set dimers. Functionals based on the local density approximation (LDA) and generalized gradient approximation (GGA) as well as hybrid functionals are compared for the DFT-D calculations, using coupled cluster CCSD(T) at the complete basis set (CBS) limit as the reference method. In general, the HF-D(EFP) method provides accurate equilibrium dimerization energies and intermolecular distances of hydrogen bonded systems compared to the CCSD(T)/CBS reference data, without using any empirical parameters. For dispersion dominant and mixed systems, the structures and interaction energies obtained with the B3LYP-D(EFP) method are similar or better than those obtained with the other DFT-D and HF-D methods. Overall, the first-principles derived -D(EFP) correction presents a robust alternative to the empirical -D corrections when used with the B3LYP functional for dispersion dominant and mixed systems or with Hartree-Fock for hydrogen bonded systems.

1 ACS Paragon Plus Environment

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

Introduction Noncovalent interactions are ubiquitous in chemistry, biochemistry and biology. An obvious example is the mere existence of the liquid phase. In fact, noncovalent interactions are fundamental for many solvent properties and dynamics, and as a result, affect the rates of chemical reactions in the condensed phase. In addition, many vital biological features such as the stable folded structure of proteins,1-3 the double helix structure of DNA4 and selective enzymesubstrate reactions5,6 are the result of specific noncovalent interactions between molecular entities.3,7,8 Moreover, noncovalent interactions may be exploited for molecular assembly,9,10 catalyst design11-13 and drug delivery system design.14 Noncovalent interactions include electrostatic forces and dispersion forces. In DNA, for example, the hydrogen bonding between base pairs is highly directional, therefore contributing to the stability of the double helix 3D structure.4,15,16 While dispersion forces are weaker and less directional than the electrostatic interactions in hydrogen bonding, they facilitate the π-stacking of aromatic rings in base pairs and provide additional stability to the helix structure of DNA.15,17 Although noncovalent interactions are much weaker than covalent interactions (on the order of a few kcal/mol), a meticulous description is important to reach chemical accuracy. High level ab initio calculations such as coupled cluster and configuration interaction methods are able to accurately predict noncovalent interactions, but these methods are computationally expensive and can only be applied to small systems.18 Mean field methods such as density functional theory (DFT) are now commonly used to model large molecular systems (on the order of hundreds of atoms).18-20 However, local and semi-local functionals do not incorporate dispersion interactions (a long range electron correlation effect),21-23 leading to a poor estimation of the structures and

2 ACS Paragon Plus Environment

Page 2 of 27

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

The Journal of Physical Chemistry

energies of dispersion dominant complexes.24-26 Several approaches have been developed to overcome this important limitation. For example, the van der Waals-density functional (vdWDF)25,27-31 and the Andersson-Langreth-Lundqvist (ALL)32-34 theories have shown great success in modeling dispersion interactions within DFT. Truhlar and coworkers have developed parameterized density functionals such as MPWB1K35, M06-2X36, M08-HX37, M08-SO36, M1138 that account for dispersion in regions where intermolecular overlap is small but not negligible (near van der Waals minima). Semi-empirical dispersion corrections have also been applied to traditional DFT methods. For example, Grimme and coworkers developed empirical dispersion energy corrections to the total DFT energy (the DFT-Dn methods, n=1,2,3).39-41 Even though the DFT-Dn methods are more cost effective than coupled cluster theories, these methods have the disadvantage of potentially double counting some of the electron correlation.22,42 Unlike DFT, the Hartree-Fock (HF) method does not include electron correlation except for the Fermi hole and therefore, there is no double counting in HF-Dn calculations. Of course, HF-Dn does not include short-range correlation effects besides the Fermi hole. Nonetheless, the HF-Dn method, especially HF-D3, is reasonably accurate for dimerization energies of small, organic molecules compared to experimental data.42-45 Moreover, HF-D3 calculations performed on complexes with large π-π interactions predict interaction energies that are comparable to CCSD(T) benchmarks.43 However, the HF-D3 method contains empirical parameters. Other DFT-D methods compute dispersion energy corrections using quantities derived from first principles. For instance, the local response dispersion (LRD) correction by Sato46,47 and the Tkatchenko-Scheffler (TS) correction48,49 compute the C6 dispersion coefficients from the ground state electron density.

3 ACS Paragon Plus Environment

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

The EFP method is an ab initio potential used to model intermolecular interactions.50-54 All of the EFP interaction energy terms are derived from first principles.50,55-57 The total EFP potential is composed of five terms: the electrostatic potential (Coulomb),50 polarization,50,58 dispersion,55,59,60 exchange-repulsion56,61 and charge transfer terms.57 The EFP method has been used to investigate the noncovalent interactions in dimers of benzene and benzene derivatives,24,62,63 styrene dimers,64 water-benzene,65 water-methanol,66 water-alanine,67,68 and DNA base pair complexes.16,24,52 In these studies, the EFP method demonstrated high accuracy compared to correlated ab initio methods. In order to avoid empirical parameters, the dispersion energy as implemented in the effective fragment potential (EFP) method can be used as a correction to the quantum mechanical energy (DFT or HF). This method also has the advantage that it is systematically improvable by, for instance, adding higher order terms and many-body terms. The EFP dispersion correction was previously added to the HF energy (HF-D(EFP) method) and the DFT energy using the B3LYP functional (DFT-D(EFP) method).42 These two methods were applied to a small subset of dimers from the S22 test set. The HF-D(EFP) and DFT-D(EFP) methods predicted interaction energies that are comparable to the CCSD(T) reference values. The most commonly used exchange-correlation functionals in DFT are modeled based on the local density approximation (LDA), generalized gradient approximation (GGA) or hybrids that contain HF exchange.19,20 However, the different density functionals exhibit contrasting performances in various chemical systems. For instance, GGAs offer systematic improvement over LDAs in hydrogen-bonded systems.69 On the other hand, both of these rungs on the DFT ladder fail to provide accurate descriptions for polymers, proteins and biomolecular surfaces.22,26 Hybrid functionals yield better accuracy than both LDA and GGA functionals for a wide range of chemical systems including hydrocarbons.22,70-72

4 ACS Paragon Plus Environment

Page 4 of 27

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

The Journal of Physical Chemistry

The objective of this study is to broadly assess the accuracy of the DFT–D(EFP) and HFD(EFP) methods by computing the intermolecular equilibrium distances and intermolecular interaction energies of the dimers in the S22 benchmark test set.73,74 LDA, GGA and hybrid functionals are used to compute the DFT-D(EFP) dimerization energies of the S22 test set dimers. The DFT-D(EFP) and HF-D(EFP) results are also compared with the DFT–Dn and HFDn methods of Grimme and co-workers.39-41

Methods The total energy E of the system may be written as,

E = Edisp + EQM .

(1)

Edisp is the dispersion energy calculated using the -D(EFP) method or the Grimme -Dn (n=1,2,3) corrections.39-41 Both methods are collectively referred to as –D methods. EQM is the quantum mechanical energy (Kohn-Sham or Hartree-Fock) of the system. In the -D(EFP) method, the dispersion energy is derived from second order perturbation theory. Within the isotropic and spherical approximations, the dispersion energy between two fragments (or monomers) A and B is written as:42,55,75 ∞

E

AB disp

ij In Eq.(2), C6 = −

=−

3h ∞

3h

π

∑∑

i

j

∫ α ( iω )α ( iω ) dω 0

Rij6

i∈ A j∈B

i

C6ij = ∑∑ 6 i∈ A j∈B R ij

(2)

j

α ( iω )α ( iω ) dω .The i and j labels represent the localized π ∫ 0

i

j

molecular orbitals (LMOs) of molecules A and B, respectively. α ( iω ) and α ( iω ) represent one third of the trace of the dynamic polarizability tensors computed at LMO centroids i and j, respectively. Rij is the distance between the two LMO centroids i and j of fragments A and B 5 ACS Paragon Plus Environment

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

Page 6 of 27

AB is multiplied by an overlap-based respectively. In order to avoid a singularity near Rij=0, Edisp



damping function f dmp ( S ij ) ,42,76 6

fdmp ( S ij ) = 1 − S ij2 ∑

( −2ln S ) ij

n!

n =0

n/2

.

(3)

AB Sij in Eq. (3) is the overlap between the LMOs i and j. In order to compute Edisp , the



LMOs of each fragment are generated from the HF wavefunction using the EdmistonRuedenberg scheme77 and the distributed dynamic dipole polarizabilities are obtained by solving the time-dependent Hartree-Fock equations.55,75 The -D(EFP) dispersion energy correction is therefore systematic without any empirically fitted parameters. Similar to the EFP method, monomers are currently kept internally rigid in all DFT-D(EFP) and HF-D(EFP) calculations. A more detailed derivation of the EFP dispersion energy was discussed by Guidez and Gordon42 and by Adamovic and coworkers.55 In the DFT-Dn and HF-Dn methods, dispersion interactions are computed between all atom pairs. In the -D1 and -D2 methods, the C6 coefficients and the damping functions contain empirically fitted parameters.39,40 In addition, the C6 coefficients are multiplied by an empirically derived, functional-dependent scaling factor.39,40 For the -D3 correction, the C6 coefficients are computed using time-dependent density functional theory (TDDFT). Recursion relations are used to compute the higher order dispersion energy terms C8/R8, which are scaled by an empirically derived factor. Two additional empirically derived parameters appear in the damping functions of the R-6 and R-8 dispersion terms.41 In order to be consistent with the DFTD(EFP) and HF-D(EFP) calculations, monomers are also kept rigid in all DFT-Dn and HF-Dn calculations. The DFT-D(EFP) dimerization energy is calculated using the following formula: 6 ACS Paragon Plus Environment

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

The Journal of Physical Chemistry

DFT − D( EFP ) DFT ∆E = Edim − ∑ Emonomer , er

(4)

DFT −D( EFP ) In Eq.(4), Edim is the energy of the dimer computed with the DFT-D(EFP) method and er

DFT Emonomer is the DFT energy of a monomer. Intramolecular dispersion interactions are not

computed in this method. An analogous equation can be written for the HF-D(EFP) method. The DFT-Dn dimerization energy is calculated using the following formula: DFT −Dn DFT − Dn ∆E = Edim − ∑ Emonomer , er

(5)

DFT − Dn DFT − Dn In Eq.(5), Edim er is the energy of the dimer computed with the DFT-Dn method. E monomer is the

DFT-Dn energy of a monomer. An analogous equation can be written for the HF-Dn method. The potential energy curves presented in this study show the dimerization energies as a function of the distance between the centers of mass of the monomers (also called the intermolecular distance). Using the reference geometries from the S22 database as a starting point,73 these potential energy curves were generated by displacing the monomers along the axis that connects their centers of mass. Only single-point energy calculations were performed for all DFT-D and HF-D methods investigated in this study since gradients are not yet implemented in the DFT-D(EFP) method. The LMOs and dynamic polarizability tensors were generated for each monomer at its geometry in the dimer. Monomers were kept internally rigid for all DFT-D and HF-D calculations. The equilibrium distances and energies are estimated to be at the minimum of the potential energy curves. The DFT-D(EFP) and HF-D(EFP) equilibrium distances and interaction energies are compared to the reference CCSD(T) and MP2 results.74,78

Computational details

7 ACS Paragon Plus Environment

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

The data reported in this paper were obtained using the GAMESS software package.79,80 The S22 test set molecules are divided into three categories based on the dominant type of interaction: hydrogen bonding dominant, dispersion dominant and mixed interactions.78 Popular DFT functionals with different approximations are tested in this study: B3LYP70,81 and PBE082 (hybrids), SVWN83 (LDA), and PBE84 (GGA). These functionals are used with the 6-311++G(d,p) basis set. The size of the basis set used may affect the calculated geometry and interaction energies. Therefore, for the HF-D(EFP) method, both the 6311++G(d,p) and 6-311G(d,p) basis sets are used to examine how the addition of diffuse functions85 affects the total energy. As explained in the Methods section, the reference geometries from the S22 database were used as a starting point to generate the potential energy curves. The reference structures of the water, ammonia, formic acid, ethene, ethene-ethyne and methane complexes were optimized with CCSD(T)/cc-pVQZ without counterpoise corrections, while the other complexes were optimized with MP2/cc-pVTZ CP.73,74 The equilibrium dimerization energy ∆Ee corresponds to the dimerization energy ∆E (eq.(4) and eq(5)) obtained at the minimum of the potential energy curve Re. The Re and ∆Ee values are reported for all of the methods and all of the dimers, and compared with high accuracy reference data.74,78 In order to examine the accuracy of the –D(EFP) and –Dn methods, the mean unsigned errors (MUE) and root mean square errors (RMS) were calculated for both the equilibrium dimerization energies and the equilibrium intermolecular distances. The MUE is calculated to measure how close DFT-D or HF-D values are to the reference values (CCSD(T) or MP2). For example, the MUE of the DFT-D equilibrium dimerization energy compared to the reference equilibrium dimerization energy is given in Eq. (6) below.

8 ACS Paragon Plus Environment

Page 8 of 27

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

The Journal of Physical Chemistry

MUE =

1 ∆Eereference − ∆EeDFT −D ∑ n

.

(6)

In Eq.(6), n is the total number of dimers of interest. ∆Eereference is the reference equilibrium dimerization energy obtained at the CCSD(T)/CBS CP level of theory for a given dimer. These values are reported on the S22 benchmark website.74,78 ∆EeDFT −D is the equilibrium dimerization energy obtained with the DFT-D method for the same dimer. The sum in Eq.(6) goes over all dimers of interest. An analogous expression can be written for the HF-D method. The root mean square (RMS) of the DFT-D equilibrium dimerization energy compared to the reference energy is calculated with the formula:

RMS =

∑ ( ∆E

reference e

− ∆E eDFT − D

)

2

n

.

(7)

Generally, the MUE offers a clear description of the average error in the test set. On the other hand, the RMS is a better statistical criterion for a data set that includes data points with significant deviations from reference values. Analogous expressions to Eq.(6) and Eq.(7) are used to compare the DFT-D and HF-D equilibrium intermolecular distances to the reference structures optimized at the MP2/cc-pVTZ CP or CCSD(T)/cc-pVQZ level of theory.73,78

Results and Discussion The equilibrium dimerization energies ∆Ee and equilibrium distances Re calculated using the DFT-Dn (n=1,2,3), DFT-D(EFP), HF-Dn (n=1,2,3), and HF-D(EFP) methods are compared with the reference values for all S22 complexes. The errors computed for the DFT-D and HF-D dimerization energies relative to the reference CCSD(T)/CBS energies are shown in Tables 1 and 9 ACS Paragon Plus Environment

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

2, respectively. The errors computed for the DFT-D and HF-D intermolecular distances relative to the reference intermolecular distances are shown in Tables 3 and 4, respectively.

Hydrogen bonded complexes The first type of complexes in the S22 test set is hydrogen bonded complexes. For hydrogen bonded complexes, the DFT-D methods overall tend to overestimate the interaction energies compared to CCSD(T) and MP2 calculations. The SVWN functional considerably overestimates interaction energies (MUE and RMS ≈10 kcal/mol) and underestimates equilibrium distances (MUE=0.16-0.3 Å and RMS=0.16-0.34 Å). This functional also displays inferior performance for dispersion dominant complexes and mixed complexes, and therefore will not be discussed further. With PBE, B3LYP and PBE0, the DFT-D(EFP) method overestimates the dimerization energy (MUE =3.51, 2.34 and 4.36 kcal/mol, respectively) and underestimates the intermolecular equilibrium distances (MUE=0.13, 0.10 and 0.15 Å, respectively) more than the DFT-Dn methods. Overall, the -D1 correction performs best for PBE and B3LYP, whereas all three -Dn corrections perform similarly for PBE0. The HF-D(EFP)/6-311++G(d,p) method yields more accurate equilibrium distances (MUE=0.04 Å and RMS=0.05 Å ) than the HF-Dn/6-311G++(d,p) methods (MUE and RMS >0.07 Å). Furthermore, HF-D(EFP)/6-311++G(d,p) yields interaction energies that are similar to those obtained with the HF-Dn methods and within chemical accuracy (within about 1 kcal/mol). It is also important to emphasize that the presence of diffuse functions in the basis set may be important for these hydrogen-bonded systems. The dimerization energies in Tables 2 and 4 show that the absence of diffuse functions in the basis set for the HF-D/6-311G(d,p) calculations leads to more negative dimerization energies than the HF-D/6-311++G(d,p) values. Potential energy curves for the water dimer (Figure 1A) and the ammonia complex (Figure 1B) illustrate this 10 ACS Paragon Plus Environment

Page 10 of 27

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

The Journal of Physical Chemistry

behavior. Potential energy curves for the other hydrogen bonded complexes represent similar behaviors as the water and the ammonia complex (Figure S1A and S2A of the SI). In general, D corrections are more appropriate when used with Hartree-Fock than with DFT, as there is no double counting of the electron correlation for HF. In fact, the ∆Ee and Re values computed with HF-D(EFP)/6-311++G(d,p) are closer to the CCSD(T) reference values (MUE=0.44 kcal/mol and 0.04 Å respectively) than all of the DFT-D(EFP) methods. In summary, for hydrogen bonded complexes, the HF-D(EFP)/6-311++G(d,p) results are in very good agreement with CCSD(T) and MP2 and within chemical accuracy.

Dispersion dominant complexes The second type of complexes in the S22 test set is dispersion dominant complexes (see Tables 1-4). Comparing the errors reported for ∆Ee and Re with DFT-D and HF-D, the B3LYP-D methods have overall the lowest errors among all calculations. B3LYP-D3 gives a MUE for the dimerization energy of 0.29 kcal/mol and a MUE for the equilibrium distance of 0.07 Å. The B3LYP-D(EFP) method demonstrates a similar performance, with a MUE of 0.49 kcal/mol for ∆Ee and a MUE of 0.04 Å for Re. The B3LYP-D(EFP) method shows smaller ∆Ee and Re errors than all of the other DFTD(EFP) methods and the HF-D(EFP) method. The -D(EFP) correction overall underestimates equilibrium intermolecular distances and overestimates interaction energies compared to the CCSD(T)/CBS CP values when used with the PBE and PBE0 functionals. With these two functionals, the -Dn corrections are more accurate. Moreover, although both PBE0 and B3LYP are hybrid functionals, PBE0-D overall tends to exhibit larger errors than B3LYP-D for both ∆Ee and Re.

11 ACS Paragon Plus Environment

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

Previous studies have shown that MP2 exhibits a tendency to overbind dispersion dominant complexes such as benzene dimers.39,86-89 For instance, Hobza, Selzle and Schlag showed that MP2 greatly overbinds polarized or mixed interaction systems such as the T-shaped benzene dimer, even with a significant increase in basis set size.89 Reference [89] also suggests that “floppy systems” like the benzene dimers (T-shaped, as well as sandwich structures and parallel-displaced benzene dimers) cannot rely on MP2 geometries. In addition, Conrad and Gordon also discussed in reference [43] that the HF-D3 interaction energies are better than MP2 interaction energies for large dispersion dominant complexes, in which MP2 is known to overbind. The results presented in this study tend to support these observations. There are eight dispersion dominant complexes in the S22 test set: the adenine-thymine stack, the benzenemethane complex, the benzene stack, the ethene dimer, the indole benzene complex stack, the methane dimer, and the pyrazine and uracil dimer stacks. All were optimized at the MP2/ccpVTZ CP level of theory except the ethene and methane dimers, which were optimized with CCSD(T)/cc-pVQZ and CCSD(T)/cc-pVTZ, respectively. MP2 intermolecular distances are consistently smaller than those obtained with HF-D. In addition, MP2 dimerization energies are systematically larger than those obtained with HF-D for all systems except the methane and ethyne dimers (cf. Tables S1 and S2 of the SI). It is noteworthy that the methane and ethyne dimers are the only two dispersion-dominant systems optimized with CCSD(T), not MP2. The HF-D method in principle accounts for the correlation in dispersion dominant systems, without double counting, and so one could assert that HF-D might be more reliable than MP2 for these systems. The HF-D(EFP) method predicts equilibrium distances and interaction energies that are slightly better than HF-D1 but less accurate than HF-D2 and HF-D3, with both the 6-

12 ACS Paragon Plus Environment

Page 12 of 27

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

The Journal of Physical Chemistry

311++G(d,p) basis set and the 6-311G(d,p) basis set. The absence of diffuse functions in HFD/6-311G(d,p) yields slightly larger intermolecular interaction energies and smaller equilibrium distances than the HF-D/6-311++G(d,p) method for the dispersion dominant complexes (Figure S1B and Figure S2B of the SI). The HF-D(EFP)/6-311G(d,p) intermolecular interaction energies and equilibrium distances are then closer to the reference values (MUE=0.81 kcal/mol and 0.16 Å respectively) than those obtained with HF-D(EFP)/6-311++G(d,p) (MUE=1.15 kcal/mol and 0.23 Å respectively). The equilibrium dimerization energies and intermolecular distances obtained with the larger 6-311++G(d,p) basis set should be more accurate than those obtained with the smaller 6-311G(d,p) basis set. For these dispersion dominant systems, the overbinding tendency of MP2 is cancelled out by the inadequacy of the smaller basis set. In general, the calculated errors for all of the –D methods for the dispersion dominant complexes embed the potential errors resulting from the overbinding tendency of MP2 for the dispersion dominant systems.

Mixed complexes: For mixed systems, the PBE-D1 method shows the lowest error for all of the dimerization energies calculated for mixed complexes and a comparably low Re error. (MUE=0.24 kcal/mol and 0.06 Å respectively). However, except for PBE-D1, the errors obtained with the B3LYP-D methods tend to be in general smaller than the errors obtained with the other DFT-D methods. Among all –D(EFP) results, B3LYP-D(EFP) has the lowest error for ∆Ee and Re (MUE=0.37 kcal/mol and 0.05 Å, respectively). The errors obtained with B3LYP-D(EFP) for the equilibrium intermolecular distances and dimerization energies are similar to or lower than those obtained with the B3LYP-Dn methods. When used with the PBE and PBE0 functionals, the DFT-D(EFP)

13 ACS Paragon Plus Environment

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

method largely overestimates ∆Ee and underestimates Re compared to CCSD(T) and MP2 (Figure S1C and Figure S2C of the SI). With these functionals, the -Dn corrections are more accurate. The equilibrium intermolecular distances obtained with the HF-D methods tend to be much larger than the reference values and the MUEs tend to be large (between 0.07 and 0.21 Å). The MP2 intermolecular distance is used as a reference for all the mixed systems except the ethene-ethyne dimer, which was optimized with CCSD(T). The known tendency of the MP2 method to overbind some of these systems (such as the T-shaped benzene dimer) and some missing correlation in the HF-D method might contribute to the large MUE values. As observed for the dispersion dominant complexes, the equilibrium distances obtained with HF-D tend to be smaller and the interaction energies tend to be larger when diffuse functions are not included. The HF-D(EFP)/6-311G(d,p) level of theory thus yields equilibrium geometries and interaction energies that are closer to the reference values. This behavior might again result from a cancellation of error between the overbinding tendency of MP2 and the inadequacy of the smaller basis set. Missing correlation in HF might also contribute to this phenomenon. Overall, the results demonstrate that the ∆Ee and Re values obtained with PBE-D1 are closest to the reference energies and geometries and within chemical accuracy. B3LYP-D(EFP) performs the best among –D(EFP) methods. However, it should again be emphasized that just like dispersion dominant systems, the reference MP2 geometries for mixed systems might not be reliable since studies have shown that for benzene dimers, including the T-shaped dimer, dimerization energies calculated with MP2 are overestimated and equilibrium distances are underestimated compared to CCSD(T).90 In addition, DNA base pair studies by Hobza and

14 ACS Paragon Plus Environment

Page 14 of 27

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

The Journal of Physical Chemistry

coworkers showed that there is a substantial discrepancy in dimerization energy between MP2/CBS CP and CCSD(T)/CBS CP, which demonstrates that the MP2 level of theory is not sufficient when a significant dispersion contribution is expected.74 It is therefore important to correct for higher order correlation effects that are not provided by the MP2 level of theory, when non-negligible dispersion is expected, as is the case for dispersion dominant and mixed systems.74,90 The reference dimerization energies and equilibrium distances obtained for the MP2 optimized dispersion dominant and mixed complexes are expected to be less accurate than those for the CCSD(T) optimized complexes.

Conclusions The effective fragment potential (EFP) dispersion energy has been applied as a correction to the DFT (DFT-D(EFP)) and Hartree-Fock (HF-D(EFP)) energies. The DFT-D(EFP) and HFD(EFP) methods are compared to the semi-empirical DFT-Dn (n=1,2,3) and HF-Dn (n=1,2,3) methods and to the MP2 and high accuracy CCSD(T) method for the dimers of the S22 test set. For hydrogen bonded systems, the dimerization energies calculated with HF-D(EFP)/6311++G(d,p) and HF-Dn (n=2,3) have good accuracy (even higher than MP2) compared to CCSD(T)/CBS CP results. Importantly, contrary to the HF-Dn methods, there are no fitted parameters in HF-D(EFP). For dispersion dominant and mixed complexes, the B3LYP-D(EFP) dimerization energies and intermolecular distances are comparable to or better than those obtained with the DFT-Dn methods, without having to resort to fitted parameters. However, MP2 tends to overbind these systems, making the MP2 reference geometries unreliable. The absence of diffuse functions in the basis set for the HF-D(EFP) calculations results in larger equilibrium interaction energies and smaller intermolecular distances. The equilibrium distance calculations

15 ACS Paragon Plus Environment

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

performed in this study demonstrate the good accuracy of the -D(EFP) correction compared to the popular semi-empirical –Dn (n = 1,2,3) dispersion correction approaches. Finally, while it is true that the -D(EFP) correction is more expensive than empirical corrections, the computational cost is trivial compared to the ab initio part of the calculation.

Supporting Information: The supporting information is available free of charge on the ACS Publications website. Figure S1. DFT-D3 and HF-D(EFP) potential energy curves for the S22 test set. Potential curves generated by DFT-D3 (DFT: SVWN, B3LYP, PBE, PBE0), HF-D(EFP) /6-311G(d,p), HF-D(EFP)/6-311++G(d,p), CCSD(T)/6-311++G(d,p) and MP2/6-311++G(d,p).

Figure S2. DFT-D(EFP) and HF-D(EFP)potential energy curves for the S22 test set. Potential curves generated by DFT-D(EFP) (DFT: SVWN, B3LYP, PBE, PBE0), HFD(EFP)/6-311G(d,p), HF-D(EFP)/6-311++G(d,p), CCSD(T)/6-311++G(d,p) and MP2/6311++G(d,p).

Table S1 Intermolecular dimerization energies (∆Ee) in kcal/mol Intermolecular dimerization energies (∆Ee) in kcal/mol are calculated with the DFT-D and HF-D methods for all 22 complexes.

Table S2. Equilibrium distances (Re) in Å for the S22 test set. Equilibrium distances (Re) in Å are calculated with the DFT-D and HF-D methods for all 22 complexes.

16 ACS Paragon Plus Environment

Page 16 of 27

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

The Journal of Physical Chemistry

Acknowledgements. This work was supported in part by a National Science Foundation Software Infrastructure (SI2) grant (ACI-1047772) and in part by a grant from the Air Force Office of Scientific Research (FA9550-14-1-0306). C.K. participated in this research as part of the Iowa State university freshman honors program.

References (1) Kumar, S.; Nussinov, R. Close-Range Electrostatic Interactions in Proteins. ChemBioChem 2002, 3, 604-617. (2) Koumanov, A.; Ladenstein, R.; Karshikoff, A. Electrostatic Interactions in Proteins: Contribution to Structure-Function Relationships and Stability. Recent Res. Dev. Protein Eng. 2001, 1, 123-148. (3) Vondrasek, J.; Kubar, T.; Jenney, F. E.; Adams, M. W. W.; Kozisek, M.; Cerny, J.; Sklenar, V.; Hobza, P. Dispersion Interactions Govern the Strong Thermal Stability of a Protein. Chem. Eur. J. 2007, 13, 9022-9027. (4) Černý, J.; Hobza, P. Non-Covalent Interactions in Biomacromolecules. Phys. Chem. Chem. Phys. 2007, 9, 5291-5303. (5) Hargis, J. C.; Vankayala, S. L.; White, J. K.; Woodcock, H. L. Identification and Characterization of Noncovalent Interactions That Drive Binding and Specificity in DDPeptidases and β-Lactamases. J. Chem. Theory Comput. 2014, 10, 855-864. (6) Kollman, P. A.; Kuhn, B.; Peräkylä, M. Computational Studies of EnzymeCatalyzed Reactions:  Where Are We in Predicting Mechanisms and in Understanding the Nature of Enzyme Catalysis? J. Phys. Chem. B 2002, 106, 1537-1542. (7) Hunter, C. A.; Singh, J.; Thornton, J. M. Π-Π Interactions: The Geometry and Energetics of Phenylalanine-Phenylalanine Interactions in Proteins. J. Mol. Biol. 1991, 218, 837846. (8) Wheeler, S. E.; Bloom, J. W. G. Toward a More Complete Understanding of Noncovalent Interactions Involving Aromatic Rings. J. Phys. Chem. A 2014, 118, 6133-6147. (9) Aakeröy, C. B.; Schultheiss, N. In Making Crystals by Design: Methods, Techniques and Applications; Wiley-VCH Verlag GmbH & Co. KGaA: 2006.. (10) Gazit, E. Self-Assembled Peptide Nanostructures: The Design of Molecular Building Blocks and Their Technological Utilization. Chem. Soc. Rev. 2007, 36, 1263-1269. (11) Wheeler, S. E.; Seguin, T. J.; Guan, Y.; Doney, A. C. Noncovalent Interactions in Organocatalysis and the Prospect of Computational Catalyst Design. Acc. Chem. Res. 2016, 49, 1061-1069. (12) Neel, A. J.; Hilton, M. J.; Sigman, M. S.; Toste, F. D. Exploiting Non-Covalent Π Interactions for Catalyst Design. Nature (London, U. K.) 2017, 543, 637-646.

17 ACS Paragon Plus Environment

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

(13) Knowles, R. R.; Jacobsen, E. N. Attractive Noncovalent Interactions in Asymmetric Catalysis: Links between Enzymes and Small Molecule Catalysts. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 20678-20685. (14) Zhou, P.; Huang, J.; Tian, F. Specific Noncovalent Interactions at Protein-Ligand Interface: Implications for Rational Drug Design. Curr. Med. Chem. 2012, 19, 226-238. (15) Černý, J.; Kabeláč, M.; Hobza, P. Double-Helical → Ladder Structural Transition in the B-DNA Is Induced by a Loss of Dispersion Energy. J. Am. Chem. Soc. 2008, 130, 16055-16059. (16) Smith, Q. A.; Gordon, M. S.; Slipchenko, L. V. Effective Fragment Potential Study of the Interaction of DNA Bases. J. Phys. Chem. A 2011, 115, 11269-11276. (17) Riley, K. E.; Hobza, P. Noncovalent Interactions in Biochemistry. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 3-17. (18) Jensen, F. Introduction to Computational Chemistry; Wiley, 2007. (19) Levine, I. N. Quantum Chemistry; Pearson Prentice Hall Upper Saddle River, NJ, 2009; Vol. 6. (20) Koch, W.; Holthausen, M. C. Front Matter and Index in A Chemist's Guide to Density Functional Theory; Wiley-VCH Verlag GmbH: 2001. (21) Kristyán, S.; Pulay, P. Can (Semi)Local Density Functional Theory Account for the London Dispersion Forces? Chem. Phys. Lett. 1994, 229, 175-180. (22) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289-320. (23) Allen, M. J.; Tozer, D. J. Helium Dimer Dispersion Forces and Correlation Potentials in Density Functional Theory. J. Chem. Phys. 2002, 117, 11113-11120. (24) Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. Accurate Prediction of Noncovalent Interaction Energies with the Effective Fragment Potential Method: Comparison of Energy Components to Symmetry-Adapted Perturbation Theory for the S22 Test Set. J. Chem. Theory Comput. 2012, 8, 2835-2843. (25) 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, 246401. (26) Ortmann, F.; Bechstedt, F.; Schmidt, W. G. Semiempirical Van Der Waals Correction to the Density Functional Description of Solids and Molecular Structures. Phys. Rev. B 2006, 73, 205101. (27) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy Van Der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (28) Jiří, K.; David, R. B.; Angelos, M. Chemical Accuracy for the Van Der Waals Density Functional. J. Phys.: Condens. Matter 2010, 22, 022201. (29) Kristian, B.; Valentino, R. C.; Kyuho, L.; Elsebeth, S.; Thonhauser, T.; Per, H.; Bengt, I. L. Van Der Waals Forces in Density Functional Theory: A Review of the VDW-DF Method. Rep. Prog. Phys. 2015, 78, 066501. (30) Vydrov, O. A.; Van Voorhis, T. Nonlocal Van Der Waals Density Functional Made Simple. Phys. Rev. Lett. 2009, 103, 063004. (31) Vydrov, O. A.; Van Voorhis, T. Improving the Accuracy of the Nonlocal Van Der Waals Density Functional with Minimal Empiricism. J. Chem. Phys. 2009, 130, 104105. (32) Andersson, Y.; Langreth, D. C.; Lundqvist, B. I. Van Der Waals Interactions in Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 102-105. 18 ACS Paragon Plus Environment

Page 18 of 27

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

The Journal of Physical Chemistry

(33) Sato, T.; Tsuneda, T.; Hirao, K. A Density-Functional Study on Π-Aromatic Interaction: Benzene Dimer and Naphthalene Dimer. J. Chem. Phys. 2005, 123, 104307. (34) Sato, T.; Tsuneda, T.; Hirao, K. Long-Range Corrected Density Functional Study on Weakly Bound Systems: Balanced Descriptions of Various Types of Molecular Interactions. J. Chem. Phys. 2007, 126, 234114. (35) Zhao, Y.; Truhlar, D. G. Exploring the Limit of Accuracy of the Global Hybrid Meta Density Functional for Main-Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2008, 4, 1849-1868. (36) Zhao, Y.; Truhlar, D. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215-241. (37) 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, 364382. (38) Peverati, R.; Truhlar, D. G. Improving the Accuracy of Hybrid Meta-GGGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810-2817. (39) Grimme, S. Accurate Description of Van Der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463-1473. (40) Grimme, S. Semiempirical Gga-Type Density Functional Constructed with a LongRange Dispersion Correction. J. Comput. Chem. 2006, 27, 1787-1799. (41) 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. (42) Guidez, E. B.; Gordon, M. S. Dispersion Correction Derived from First Principles for Density Functional Theory and Hartree–Fock Theory. J. Phys. Chem. A 2015, 119, 2161-2168. (43) Conrad, J. A.; Gordon, M. S. Modeling Systems with Π–Π Interactions Using the Hartree–Fock Method with an Empirical Dispersion Correction. J. Phys. Chem. A 2015, 119, 5377-5385. (44) Goerigk, L.; Collyer, C. A.; Reimers, J. R. Recommending Hartree–Fock Theory with London-Dispersion and Basis-Set-Superposition Corrections for the Optimization or Quantum Refinement of Protein Structures. J. Phys. Chem. B 2014, 118, 14612-14626. (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) Sato, T.; Nakai, H. Density Functional Method Including Weak Interactions: Dispersion Coefficients Based on the Local Response Approximation. J. Chem. Phys. 2009, 131, 224104. (47) Sato, T.; Nakai, H. Local Response Dispersion Method. II. Generalized Multicenter Interactions. J. Chem. Phys. 2010, 133, 194101. (48) Tkatchenko, A.; Scheffler, M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. 19 ACS Paragon Plus Environment

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

(49) Kronik, L.; Tkatchenko, A. Understanding Molecular Crystals with DispersionInclusive Density Functional Theory: Pairwise Corrections and Beyond. Acc. Chem. Res. 2014, 47, 3208-3216. (50) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An Effective Fragment Method for Modeling Solvent Effects in Quantum Mechanical Calculations. J. Chem. Phys. 1996, 105, 1968-1986. (51) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646-9663. (52) Ghosh, D.; Kosenkov, D.; Vanovschi, V.; Williams, C. F.; Herbert, J. M.; Gordon, M. S.; Schmidt, M. W.; Slipchenko, L. V.; Krylov, A. I. Noncovalent Interactions in Extended Systems Described by the Effective Fragment Potential Method: Theory and Application to Nucleobase Oligomers. J. Phys. Chem. A 2010, 114, 12739-12754. (53) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632-672. (54) Pruitt, S. R.; Bertoni, C.; Brorsen, K. R.; Gordon, M. S. Efficient and Accurate Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2786-2794. (55) Adamovic, I.; Gordon, M. S. Dynamic Polarizability, Dispersion Coefficient C6 and Dispersion Energy in the Effective Fragment Potential Method. Mol. Phys. 2005, 103, 379-387. (56) Jensen, J. H.; Gordon, M. S. An Approximate Formula for the Intermolecular Pauli Repulsion between Closed Shell Molecules. Mol. Phys. 1996, 89, 1313-1325. (57) Li, H.; Gordon, M. S.; Jensen, J. H. Charge Transfer Interaction in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 124, 214108. (58) Li, H.; Netzloff, H. M.; Gordon, M. S. Gradients of the Polarization Energy in the Effective Fragment Potential Method. J. Chem. Phys. 2006, 125, 194103. (59) Xu, P.; Zahariev, F.; Gordon, M. S. The R–7 Dispersion Interaction in the General Effective Fragment Potential Method. J. Chem. Theory Comput. 2014, 10, 1576-1587. (60) Guidez, E. B.; Xu, P.; Gordon, M. S. Derivation and Implementation of the Gradient of the R–7 Dispersion Interaction in the Effective Fragment Potential Method. J. Phys. Chem. A 2016, 120, 639-647. (61) Li, H.; Gordon, M. Gradients of the Exchange-Repulsion Energy in the General Effective Fragment Potential Method. Theor. Chem. Acc. 2006, 115, 385-390. (62) Smith, Q. A.; Gordon, M. S.; Slipchenko, L. V. Benzene−Pyridine InteracVons Predicted by the Effective Fragment Potential Method. J. Phys. Chem. A 2011, 115, 4598-4609. (63) Slipchenko, L. V.; Gordon, M. S. Electrostatic Energy in the Effective Fragment Potential Method: Theory and Application to Benzene Dimer. J. Comput. Chem. 2007, 28, 276291. (64) Adamovic, I.; Li, H.; Lamm, M. H.; Gordon, M. S. Modeling Styrene−Styrene Interactions. J. Phys. Chem. A 2006, 110, 519-525. (65) Slipchenko, L. V.; Gordon, M. S. Water−Benzene InteracVons: An EffecVve Fragment Potential and Correlated Quantum Chemistry Study. J. Phys. Chem. A 2009, 113, 2092-2102. (66) Adamovic, I.; Gordon, M. S. Methanol−Water Mixtures:  A MicrosolvaVon Study Using the Effective Fragment Potential Method. J. Phys. Chem. A 2006, 110, 10267-10273. 20 ACS Paragon Plus Environment

Page 20 of 27

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

The Journal of Physical Chemistry

(67) Mullin, J. M.; Gordon, M. S. Water and Alanine: From Puddles(32) to Ponds(49). J. Phys. Chem. B 2009, 113, 14413-14420. (68) Mullin, J. M.; Gordon, M. S. Alanine: Then There Was Water. J. Phys. Chem. B 2009, 113, 8657-8669. (69) Hamann, D. R. H2o Hydrogen Bonding in Density-Functional Theory. Phys. Rev. B 1997, 55, R10157-R10160. (70) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098-3100. (71) Becke, A. D. Density‐Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648-5652. (72) Zhou, Y.; Wu, J.; Xu, X. How Well Can B3LYP Heats of Formation Be Improved by Dispersion Correction Models? Theor. Chem. Acc. 2016, 135, 44. (73) 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. (74) Jurečka, P.; Šponer, J.; Č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-1993. (75) Stone, A. The Theory of Intermolecular Forces; 2nd ed.; Oxford University Press: Oxford, 2013. (76) Slipchenko, L. V.; Gordon, M. S. Damping Functions in the Effective Fragment Potential Method. Mol. Phys. 2009, 107, 999-1016. (77) Edmiston, C.; Ruedenberg, K. Localized Atomic and Molecular Orbitals. Rev. Mod. Phys. 1963, 35, 457-464. (78) Řezáč, J. J., P.; Riley, K. E.; Cerny, J.; Valdes, H.; Pluhackova, K.; Berka, K.; Rezac, T.; Pitonak, M.; Vondrasek, J.; Hobza, P. 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. (Accessed Sep. 26th, 2017) (79) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347-1363. (80) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: Gamess a Decade Later in Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005. (81) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti CorrelationEnergy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785-789. (82) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for Mixing Exact Exchange with Density Functional Approximations. J. Chem. Phys. 1996, 105, 9982-9985. (83) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J. Phys. 1980, 58, 1200-1211. (84) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865-3868. 21 ACS Paragon Plus Environment

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

(85) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent Molecular Orbital Methods 25. Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80, 3265-3269. (86) Burns, L. A.; Marshall, M. S.; Sherrill, C. D. Appointing Silver and Bronze Standards for Noncovalent Interactions: A Comparison of Spin-Component-Scaled (SCS), Explicitly Correlated (F12), and Specialized Wavefunction Approaches. J. Chem. Phys. 2014, 141, 234111. (87) Takatani, T.; Hohenstein, E. G.; Sherrill, C. D. Improvement of the CoupledCluster Singles and Doubles Method Via Scaling Same- and Opposite-Spin Components of the Double Excitation Correlation Energy. J. Chem. Phys. 2008, 128, 124111. (88) Hopkins, B. W.; Tschumper, G. S. Ab Initio Studies of Π···Π Interactions:  The Effects of Quadruple Excitations. J. Phys. Chem. A 2004, 108, 2941-2948. (89) Hobza, P.; Selzle, H. L.; Schlag, E. W. Potential Energy Surface for the Benzene Dimer. Results of Ab Initio CCSD(T) Calculations Show Two Nearly Isoenergetic Structures:  TShaped and Parallel-Displaced. J. Phys. Chem. 1996, 100, 18790-18794. (90) 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.

22 ACS Paragon Plus Environment

Page 22 of 27

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

The Journal of Physical Chemistry

Table 1. Errors for intermolecular dimerization energy (∆Ee) in kcal/mol for S22 test set using DFT-D methods relative to CCSD(T)/CBS CP reference energies SVWN PBE B3LYP PBE0 D2 D3 D(EFP) D1 D2 D3 D(EFP) D1 D2 D3 D1 D2 D3 D(EFP) D1

D(EFP)

Hydrogen bonding MUE RMS

8.21 8.59

9.67 10.1

10.02 10.55

11.94 12.6

0.73 0.92

1.59 1.68

1.2 1.33

3.51 3.76

0.37 0.54

2.31 3.34

1.16 1.22

2.34 2.51

1.67 1.74

1.61 1.64

1.64 1.71

4.36 4.64

Dispersion Dominant 5.71 6.38

10.13 11.58

10.23 11.7

9.15 11.48

0.69 0.91

0.5 0.55

0.3 0.33

1.54 1.81

0.97 1.11

0.55 0.87

0.29 0.44

0.49 0.51

1.07 1.2

0.3 0.41

0.54 0.66

2.05 2.37

MUE

4.5

6.2

6.23

5.93

0.24

0.95

0.64

0.99

0.35

0.85

0.62

0.37

0.89

0.7

0.88

1.46

RMS

4.75

6.50

6.55

6.26

0.30

0.99

0.68

1.04

0.47

0.93

0.71

0.43

0.94

0.76

0.92

1.54

MUE RMS

Mixed

Table 2. Errors for intermolecular dimerization energy (∆Ee) in kcal/mol for S22 test set using HF-D methods relative to CCSD(T)/CBS CP reference energies

D1

Hydrogen bonding 1.21 MUE 1.38 RMS Dispersion Dominant 1.3 MUE 1.54 RMS Mixed 0.51 MUE 0.62 RMS

HF/6-311++G(d,p) D2 D3

D(EFP)

D1

HF/6-311G(d,p) D2 D3

D(EFP)

0.5 0.57

0.33 0.44

0.44 0.6

0.82 0.95

1.12 1.21

1.29 1.41

1.27 1.47

0.63 0.88

0.71 1.03

1.15 1.28

0.9 1.17

0.73 0.9

1.46 2.34

0.81 0.83

0.27 0.39

0.85 0.92

0.5 0.6

0.31 0.37

0.78 0.88

1.17 1.25

0.39 0.38

23 ACS Paragon Plus Environment

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

Page 24 of 27

Table 3. Errors for equilibrium distances (Re) in Å for S22 test set using the DFT-D methods. SVWN PBE B3LYP D1 D2 D3 D(EFP) D1 D2 D3 D(EFP) D1 D2 D3 D(EFP)

D1

PBE0 D2 D3

D(EFP)

Hydrogen bonding MUE RMS

0.16 0.16

0.16 0.17

0.17 0.18

0.3 0.34

0.02 0.03

0.03 0.04

0.02 0.03

0.13 0.14

0.02 0.03

0.02 0.02

0.15 0.37

0.1 0.1

0.03 0.03

0.04 0.04

0.03 0.04

0.15 0.16

Dispersion Dominant 0.28 0.3

0.42 0.43

0.34 0.35

0.52 0.52

0.09 0.1

0.08 0.1

0.11 0.12

0.16 0.16

0.11 0.12

0.09 0.13

0.07 0.08

0.04 0.06

0.27 0.53

0.25 0.5

0.23 0.53

0.34 0.49

MUE

0.27

0.32

0.3

0.46

0.06

0.08

0.05

0.15

0.08

0.08

0.05

0.05

0.05

0.07

0.05

0.19

RMS

0.28

0.33

0.31

0.46

0.08

0.1

0.07

0.16

0.1

0.09

0.06

0.06

0.06

0.08

0.05

0.19

MUE RMS

Mixed

Table 4. Errors for equilibrium distances (Re) in Å for S22 test set using the HF-D methods. D1

Hydrogen bonding 0.11 MUE 0.12 RMS Dispersion Dominant 0.26 MUE 0.3 RMS Mixed 0.21 MUE 0.22 RMS

HF/6-311++G(d,p) D2 D3

D(EFP)

D1

HF/6-311G(d,p) D2 D3

D(EFP)

0.07 0.07

0.1 0.1

0.04 0.05

0.09 0.09

0.05 0.05

0.06 0.07

0.03 0.05

0.16 0.2

0.17 0.19

0.23 0.25

0.22 0.25

0.13 0.16

0.13 0.15

0.16 0.18

0.07 0.08

0.14 0.15

0.18 0.19

0.16 0.17

0.2 0.41

0.09 0.11

0.13 0.15

24 ACS Paragon Plus Environment

Page 25 of 27

A

Water 0 2.5

3

3.5

4

4.5

5

5.5

-2 Energy (kcal/mol)

SVWN-D3 B3LYP-D3

-4

B3LYP-D(EFP) HF-D(EFP)/6-311++G(d,p)

-6

HF-D(EFP)/6-311G(d,p) -8

PBE-D3 PBE0-D3

-10

CCSD(T)/6-311++G(d,p)

-12

Intermoledular Distance (Å)

Ammonia

B 10 8

Energy (kcal/mol)

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

The Journal of Physical Chemistry

6

SVWN-D3

4

B3LYP-D3

2

B3LYP-D(EFP)

0

HF-D(EFP)/6-311++G(d,p)

-2

2.5

3

3.5

4

4.5

5

5.5

HF-D(EFP)/6-311G(d,p) PBE-D3

-4

PBE0-D3

-6

CCSD(T)/6-311++G(d,p) -8 -10 Intermoledular Distance (Å)

Figure 1. Potential energy curve of hydrogen-bonded complex for water (Figure 1A) and ammonia (Figure 1B). SVWN, B3LYP, PBE, and PBE0 are used for DFT calculations and HF calculations are compared. 6-311++G(d,p) basis sets are used for all DFTs and HF, and 6311G(d,p) basis set is used additionally to HF.

25 ACS Paragon Plus Environment

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

TOC GRAPHIC

26 ACS Paragon Plus Environment

Page 26 of 27

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

The Journal of Physical Chemistry

27 ACS Paragon Plus Environment