Hard Numbers for Large Molecules: Toward Exact Energetics for

Letter pubs.acs.org/JPCL

Hard Numbers for Large Molecules: Toward Exact Energetics for Supramolecular Systems Alberto Ambrosetti,† Dario Alfè,‡ Robert A. DiStasio, Jr.,§ and Alexandre Tkatchenko† †

Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin, Germany Department of Earth Sciences and Department of Physics and Astronomy and London Centre for Nanotechnology and Thomas Young Centre@UCL, University College, London WC1E6BT, United Kingdom § Department of Chemistry, Princeton University, Princeton, New Jersey 08544, United States ‡

S Supporting Information *

ABSTRACT: Noncovalent interactions are ubiquitous in molecular and condensed-phase environments, and hence a reliable theoretical description of these fundamental interactions could pave the way toward a more complete understanding of the microscopic underpinnings for a diverse set of systems in chemistry and biology. In this work, we demonstrate that recent algorithmic advances coupled to the availability of large-scale computational resources make the stochastic quantum Monte Carlo approach to solving the Schrödinger equation an optimal contender for attaining “chemical accuracy” (1 kcal/mol) in the binding energies of supramolecular complexes of chemical relevance. To illustrate this point, we considered a select set of seven host−guest complexes, representing the spectrum of noncovalent interactions, including dispersion or van der Waals forces, π−π stacking, hydrogen bonding, hydrophobic interactions, and electrostatic (ion−dipole) attraction. A detailed analysis of the interaction energies reveals that a complete theoretical description necessitates treatment of terms well beyond the standard London and Axilrod−Teller contributions to the van der Waals dispersion energy. SECTION: Molecular Structure, Quantum Chemistry, and General Theory

S

A crucial step in the control and rational design of host− guest complexes is therefore an accurate theoretical description of this underlying set of noncovalent interactions in the absence of complex temperature and environment effects, that is, “cleanroom conditions”, which could provide direct access to the energetics of these supramolecular systems. However, the large size of most functional host−guest complexes of chemical relevance poses enormous challenges for current theoretical methodologies in terms of both accuracy and computational feasibility. In this regard, high-level quantum chemistry methods such as full configuration interaction (FCI) or coupled cluster theory with single, double, and perturbative triple excitations (CCSD(T)) could certainly provide highly accurate binding energies for such systems. In fact, a number of databases of binding energies computed at the CCSD(T) level of theory for small molecular dimers (containing up to a few dozen atoms in size) have recently become available.10−12 However, the steep associated computational cost (scaling as N7 for CCSD(T), where N is a measure of the system size), makes the application of such high-level quantum chemistry

upramolecular complexes have long been recognized for their remarkable versatility1−4 and have therefore become increasingly utilized in a vast array of practical applications, including molecular recognition, self-assembly, templatedirected synthesis, and biomimetics.4−8 Countless realizations of supramolecular complexes exist5,9 and typically consist of molecular assemblies stabilized by cooperative binding motifs, with energetic contributions arising from strong covalent and ionic bonds as well as weaker nonbonded intermolecular forces. Therefore, a central problem emergent in supramolecular chemistry is characterization and subsequent control of the delicate balance between the different underlying interactions that determine the relative stability of such systems. Of particular importance in supramolecular chemistry is the class of host−guest complexes, composed of a host molecule, such as the so-called molecular “tweezers” or “pincers,” and a guest molecule, typically a relatively smaller organic molecule, which are primarily stabilized by noncovalent interactions. Hence, host−guest complexes serve as prototypes for molecular recognition and transient binding eventsprocesses which are primarily dictated by this underlying set of noncovalent interactions. As such, noncovalent interactions play a central role in determining the functionality of host−guest complexes, with an influence that encompasses conformational energetics, entropic contributions, and solvation effects. © 2014 American Chemical Society

Received: December 9, 2013 Accepted: February 5, 2014 Published: February 5, 2014 849

dx.doi.org/10.1021/jz402663k | J. Phys. Chem. Lett. 2014, 5, 849−855

The Journal of Physical Chemistry Letters

Letter

methods to large supramolecular systems very challenging, if not impossible, with the computational resources available today. Stochastic methods for solving the Schrödinger equation such as quantum Monte Carlo (QMC) could, in principle, be utilized to obtain the exact description of molecular groundstate wave functions and energies.13 In particular, diffusion quantum Monte Carlo (DQMC) represents an optimal approach for treating large systems14,15 because the DQMC method allows for a direct and highly accurate sampling of the ground-state electronic wave function, with a more favorable N3 computational scaling. While the DQMC approach has the ability to describe noncovalent interactions with benchmark (i.e., subchemical) accuracy in small molecular dimers,16 the applicability of DQMC to large supramolecular systems has not been systematically demonstrated to date. In this work, we considered a select set of six host−guest complexes (see Figure 1) from the recently proposed S12L database of Grimme,17,18 representing the spectrum of noncovalent interactions, including dispersion or van der Waals (vdW) forces, π−π stacking, hydrogen bonding, hydrophobic interactions, and electrostatic (ion−dipole)

attraction, and ascertained the quality of binding energies obtained via extrapolation of experimental association free energies17 with respect to benchmark binding energies computed at the DQMC level of theory. Although this comparison revealed a fair degree of overall fidelity between the extrapolated and DQMC binding energies, quantitative differences as large as 3.6 kcal/mol (e.g., for the case of the cucurbit[6]uril−butylammonium cationic complex, 6a in Figure 1), which are well above the chemical accuracy benchmark of 1 kcal/mol, persist and are indicative of the inherent limitations in the approximate corrections utilized to extract binding energies from experimentally determined association free energies. To further investigate the underlying noncovalent interactions determining the stability of host−guest complexes, we performed a many-body decomposition analysis of the longrange correlation energy in the aforementioned systems. Such an analysis is complementary to the DQMC methodology, which provides benchmark energetics for the systems considered herein, by yielding detailed physical insight into the fundamental role played by noncovalent interactions in governing supramolecular chemistry. As a result of this analysis, we found that the many-body expansion of the long-range correlation energy is slowly convergent and displays nontrivial behavior, depending on the symmetry and underlying topology of a given host−guest complex, strongly indicating that a chemically accurate theoretical description of supramolecular binding energies requires terms well beyond the standard London (two-) and Axilrod−Teller (three-) body contributions to the dispersion energy. We further investigated this point by extending our analysis to the long-range correlation energy of a double-walled carbon nanotube (DWCNT; see Figure 2). The

Figure 2. Graphical depiction of a periodic double-walled carbon nanotube (DWCNT), composed of coaxial (10,10) and (5,5) singlewalled carbon nanotubes. The supercell shown contains 900 carbon atoms.

marked anisotropy of polarization interactions in this system leads to a reduction of the interwall dispersive binding, which amounts to ∼25% with respect to the (isotropic) pairwise vdW energy. To construct an accurate reference for the energetics in host−guest complexes, we performed DQMC calculations to determine the binding energies for a subset of six complexes from the S12L database17 (see Figure 1). This subset (namely,

Figure 1. Molecular geometries of the six host−guest complexes studied in this work following the original nomenclature of Grimme in ref 17. 2a: tetracyanoquinone−tweezer (TCNQ@tweezer), 2b: 1,4dicyanobenzene−tweezer (DCB@tweezer), 4a: buckyball−catcher (C60@catcher), 5a: glycine anhydride−macrocycle (GLH@mcyle), 6a: butylammonium−cucurbit[6]uril cation (BuNH4@CB6), and 7b: 1-hydroxyadamantane−cucurbit[7]uril (ADOH@CB7). 850



Letter

reliability of the fixed node approximation was also tested through the use of different trial wave functions. As discussed in greater detail later, both the time step and nodal errors fall within the statistical uncertainties reported in Table 1 (corresponding to ±σ). While DQMC provides direct and reliable access to benchmark energetics in supramolecular systems, an indirect empirical estimate of the binding energies for the S12L host− guest complexes can also be determined from experimental association free energies, as was recently done by Grimme.17 Apart from the binding energy, the measured free energies contain many other contributions, such as entropic and solvation effects, which are often of comparable magnitude and opposite sign to the binding energy. For these reasons, the experimental association free energies are roughly one order of magnitude smaller than the corresponding binding energies. Hence, the determination of reliable binding energies from experimentally determined association free energies is a delicate task, which is only further complicated by the need to introduce a number of approximations in the computation of both entropic and solvation contributions, the accuracy of which is often difficult to assess. For example, the solvation effects for the S12L database were computed utilizing a simplified continuum solvent model, while entropic contributions were treated in the rigid rotor-harmonic oscillator (RRHO) approximation, with further approximations required to avoid divergences in the low-frequency regime. Furthermore, the configurational entropy was neglected; that is, a nondynamical, single-structure approach was used during the computation of binding energies. To assess the reliability of this approach, we also present in Table 1 the comparison between binding energies computed at the DQMC level of theory and extrapolated utilizing the aforementioned prescription. From this data set, one can immediately observe absolute differences ranging from 1.4 to 3.6 kcal/mol, with the largest deviations for complexes 2a (2.4 kcal/mol), 2b (3.3 kcal/mol), and 6a (3.6 kcal/mol). In complexes 2a and 2b, the host−guest binding is predominantly due to dispersion or vdW forces. Therefore, the host−guest interaction is expected to be responsible for the appearance of soft vibrational modes, with further low-lying vibrational modes appearing that are related to the host system (i.e., the opening and closing of the “tweezer” moiety). In this regard, anharmonicity, which is neglected in the RRHO approximation, might play an important role in the qualitative and quantitative description of these modes. In addition, the empirical interpolation between the harmonic vibrational and rotational entropy contributions (i.e., to avoid divergences in the lowfrequency regime) might also contribute to such large deviations from the DQMC values. For the complex 6a, the electrostatic cation−dipolar binding is a likely source of inaccuracy for the continuum solvent model; the characterization of the solvent by a macroscopic dielectric function might not be well-suited here due to the polarization effects occurring within such a complex asymmetric system. We note in passing that the differences between the DQMC and extrapolated binding energies encountered for every host− guest complex considered herein (and, in particular, complexes 2a, 2b, and 6a) are well above the demanding “chemical accuracy” benchmark of 1 kcal/mol. Despite this fact, the qualitative agreement among these binding energies is certainly remarkable and constitutes an important “sanity check” between two radically different approaches for determining

2a 2b, 4a, 5a, 6a, and 7b, following the original nomenclature of Grimme17) was selected to represent the broad range of geometries and noncovalent interactions of primary relevance in supramolecular chemistry, thus preserving the general character of the full S12L database. The stochastic DQMC electronic structure method is a wellestablished ab initio, or first-principles, approach to solving the Schrödinger equation and can therefore be utilized in the computation of highly accurate ground-state energies and properties.14 Because the DQMC methodology intrinsically accounts for dynamical electron correlation effects at all interelectronic separations, DQMC can be considered as a natural benchmark reference for approximate density functionals and other perturbative approaches. Exhibiting a favorable computational scaling with system size (N3), the DQMC method can optimally utilize the computational resources afforded by high-performance massively parallel (super)computer architectures, thereby enabling the challenging large-scale applications carried out in this work. All DQMC calculations presented herein have been performed utilizing the CASINO suite of programs,19 employing Slater−Jastrow trial wave functions ΨT (R) = D↑D↓e J ↑

(1) ↓

in which D and D are Slater determinants assembled from single-particle spin orbitals representing the α (up) and β (down) electron-spin projections, respectively, and eJ is the socalled Jastrow factor, an exponential comprised of a sum over explicitly correlated one- (electron−nucleus), two- (electron− electron), and three-body (electron−electron−nucleus) terms. The computed DQMC binding energies for the six host− guest complexes considered in this work are provided in Table 1. The binding energy of the host (H) and guest (G) forming Table 1. Binding Energies for Each of the Considered S12L Host−Guest Complexes in kilocalories per mole binding energies complex 2a 2b 4a 5a 6a 7b

TCNQ@tweezer DCB@tweezer C60@catcher GLH@mcyle BuNH4@CB6 ADOH@CB7

DQMCa −27.2 −17.2 −25.8 −33.4 −81.0 −24.1

(0.3)d (1.0) (1.5) (1.0) (1.6) (1.8)

extrap. expt.b

PBE +MBD*c

−29.9 −20.5 −27.5 −34.8 −77.4 −22.6

−29.0 −18.8 −28.3 −33.8 −82.1 −27.4

a

DQMC computed binding energies with associated statistical sampling uncertainties given in parentheses. bBinding energies extrapolated from experimentally measured association free energies via approximate solvation and entropic corrections.18 cBinding energies for the PBE functional including long-range correlation at the MBD* level (see the text for details). dAs visible from the lower statistical error, a longer DQMC sampling was performed for complex 2a in order to test the actual importance of statistical noise (see Supporting Information for details).

the host−guest complex (H−G) was defined as ΔE = E(H−G) − E(H) − E(G). All DQMC calculations were performed utilizing molecular geometries optimized with dispersioncorrected density functional theory (DFT) in ref 18. Convergence tests were performed to verify the dependence of the computed binding energies on the imaginary timepropagation step (see Supporting Information). In addition, the 851



Letter

considered in this work are also provided in Table 1. (All DFT calculations were performed with the FHI-aims code.25) In general, we observed very good performance across the entire S12L database, with a mean absolute relative error (MARE) of 5.5% computed with respect to the extrapolated experimental values, which is similar to the MARE obtained for smaller gasphase molecular dimers.21 In this regard, it should be emphasized that the PBE+MBD* method, with a corresponding mean absolute error (MAE) of 1.6 kcal/mol over the entire S12L database, provides an accuracy comparable to that of the reference data. In fact, the binding energies computed at the PBE+MBD* level of theory are consistently in closer agreement with the benchmark DQMC results than the extrapolated values: by performing a statistical analysis restricted to the subset of six host−guest complexes considered in this work, we found a MAE of 1.7 kcal/mol with respect to the DQMC results, which should be compared to the MAE of 2.3 kcal/mol obtained when comparing the extrapolated binding energies to the same reference DQMC values. Through this analysis, we further confirm the importance of the longrange correlation energy, the contribution of which can amount to more than 90% of the total binding,26 which is clearly an integral component of the binding energy that is not captured at the underlying DFT level of theory. The agreement between the PBE+MBD* and DQMC binding energies is essentially comparable to the statistical error of the stochastic DQMC method, except for the complexes 4a (C60@catcher) and 7b (ADOH@CB7), which can be explained by the underlying approximations employed in the PBE+MBD* method. For instance, the approximation of localized QHOs might not provide the flexibility to adequately describe the response of the delocalized electrons present in the C60 guest of complex 4a, leading to a moderate deviation (i.e., 1.0 kcal/mol outside the error bar). Regarding complex 7b, the elevated number of hydrogens in the guest pointing toward the host causes the combined appearance of weak hydrogen bonds and Pauli exchange−repulsion effects. Here the approximate treatment of exchange at the semilocal DFT level of theory combined with the inexact range separation of the Coulomb interaction certainly represent a limitation in this context. A systematic analysis of the influence of these effects on the binding energies of the entire S12L database is currently under investigation. To better analyze the role of the many-body (many-atom) interactions in the binding energies of the S12L complexes, we have carried out a many-body decomposition of the infiniteorder MBD* energy into pairwise (two-body), three-body, and higher-order terms. This many-body decomposition is facili*) tated by the fact that the MBD* correlation energy (EMBD c can be expanded in powers of the product of the bare response function χ0 with the interaction v as:20

the gas-phase binding energetics of large supramolecular systems. The binding energies computed at the DQMC level of theory and presented herein provide the most reliable benchmarks available to date for the energetics of large supramolecular systems of chemical relevance and can therefore serve as reliable references for the development (and subsequent validation) of computationally efficient approximate electronic structure methods. In fact, a main advantage of utilizing the DQMC method lies in its ability to accurately treat long-range correlation effectseffects that are inherently quantum mechanical, many-body, and nonlocal in character which pose quite a challenge for many approximate electronic structure methods. However, the DQMC method alone does not directly allow for a detailed analysis of the various electron correlation contributions to the binding energies in question. Because of the stochastic nature of the sampling of the manybody ground-state wave function, the DQMC method computes ground-state energies in a nonperturbative fashion, making a distinction among the different energy components impractical. Therefore, to gain direct physical insight into the role played by the long-range correlation energy in the stabilization of host−guest complexes, we will make combined use of an alternative and complementary approach, which allows for a detailed analysis of the pairwise and many-body contributions to the long-range correlation energy within the framework of DFT. Semilocal DFT is a self-consistent quantum-mechanical electronic structure method that accurately describes electrostatics, induction, and hybridization effects but does not include long-range electron correlation and therefore fails to account for dispersion or vdW interactions. As discussed in greater detail later, we explicitly treat the long-range correlation energy within DFT by utilizing the random-phase approximation (RPA) in the dipole limit20 (obtained through the MBD* methodsee Supporting Information) based on a rangeseparation of the interelectronic Coulomb potential. The RPA approach seamlessly includes many-body effects in the correlation energy to all orders and is a very accurate theory for the long-range correlation energy, provided that correct polarizabilities are utilized as an input. To efficiently compute the long-range RPA correlation energy, we map the molecular system onto a set of atom-centered quantum harmonic oscillators (QHOs) and utilize an effective oscillator Hamiltonian.21 With respect to the recently published MBD method,21 MBD* (called MBD@rsSCS in ref 22) offers an improved description of highly anisotropic systems, primarily due to range-separation of the Coulomb operator and a correct resultant treatment of the long-range electrodynamic response. Although the MBD* method is not expected to reach the same degree of accuracy as DQMC, performance beyond chemical accuracy has been demonstrated utilizing this methodology in both small molecules21 and extended systems.23 In this regard, the MBD* method can be regarded as complementary to DQMC, which not only provides a detailed many-body decomposition analysis of the long-range correlation energy but also allows for even more challenging large-scale applications due to its high computational efficiency. Before proceeding any further, a comparison against reliable benchmark data remains essential to assess the accuracy of the MBD* method and its predictivity for supramolecular systems. The computed PBE+MBD* (MBD* coupled to the PBE functional24) binding energies for the six host−guest complexes

EcM BD * = −

1 2π

∫0

∞

∞

dω ∑ n=2

1 Tr[(χ0 v)n ] n

(2)

The present procedure for the perturbative expansion of the MBD* long-range correlation energy differs from that followed in ref 27, which was based on an averaging of the free QHO characteristic frequencies, all of which were set to a single value, chosen to preserve the total correlation energy. In contrast, the many-body expansion utilized herein is based on a straightforward Taylor series decomposition of the logarithm term that naturally arises in the RPA correlation energy; as such, this 852



Letter

expansion allows for a clear diagrammatic interpretation of each perturbative term, now expressed in powers of χ0v. In addition, we also note that an alternative range-separation of the Coulomb interaction, with smoother and isotropic shortrange behavior, has been employed in this work.22 The second-order truncation of the series contained in eq 2 leads to the well-known C6/R6 London pairwise summation, the third-order term contains the so-called Axilrod−Teller− Muto28 (ATM) energy contribution, while the summation to infinite order corresponds to the full MBD* energy. As an expected general trend, the terms corresponding to the even powers in eq 2 are typically negative (attractive), while odd powers provide contributions of positive sign (repulsive). The alternating sign behavior was observed in all complexes, the only exception being the 2a host−guest complex (which loses this alternating trend after fifth order). An analysis of the data contained in Table 2 confirms that the many-body terms provide substantial decreases in the binding Table 2. Truncation Errors in the PBE+MBD* Long-Range Correlation Energies for the Host−Guest Complexes Considered in This Work correlation energy analysis 2a 2b 4a 5a 6a 7b

complex

2nd−∞a

2nd−3rdb

TCNQ@tweezer DCB@tweezer C60@catcher GLH@mcyle BuNH4@CB6 ADOH@CB7

12.2% 11.9% 12.9% 7.2% 12.5% 12.3%

20.7% 19.6% 22.6% 13.2% 21.0% 21.3%

a

Differences between the second- and infinite-order MBD* correlation energies (percentage-wise, with respect to the infinite-order value). b Differences between the second- and third-order MBD* correlation energies (percentage-wise, with respect to the infinite-order value).

energies. The magnitude of these decreases depends on the particular system and exceeds 10% for all systems considered, apart from 5a (7.2%). The reason for this smaller effect (encountered also in the 5b S12L complex, with a 7.5% decrease) can be attributed to a combination of effects, such as the rather sparse and symmetric conformation and the relatively low impact of dispersion forces on the overall binding in a complex that is primarily bound by hydrogen bonding. Interestingly, with the mere inclusion of the ATM (three-body) term, the binding decrease is exaggerated, as the higher-order terms provide a further reduction by about a factor of two. Hence, the neglect of higher-order (n > 3) terms often amounts to a few kilocalories per mole and will lead to errors much larger than the highly desired chemical accuracy threshold. As visible in Figure 3, the progressive decrease in the absolute values of the perturbative contributions with n only converges at relatively higher orders. The rate of convergence again depends on the system, and deviations from MBD* below 1% are usually achieved only after n ≈ 6. We stress that a deviation of 1% from MBD* can correspond to ∼1 kcal/mol for such supramolecular systems, that is, only slightly below the accuracy of the reference data. To further illustrate this point, we analyzed the MBD* longrange correlation energy of an infinite, periodic, double-walled carbon nanotube (DWCNT). This system is composed of two coaxial single-walled nanotubes, namely, of the (10,10) and (5,5) type, as seen in Figure 2. Despite the metallicity inherent

Figure 3. Many-body decomposition of the long-range MBD* correlation energy, as defined in eq 2. Results are reported for the complexes showing the slowest and fastest convergence (4a and 5a, respectively) and for an intermediate case (6a). The additional case of a periodic double-walled carbon nanotube (DWCNT) is also considered for comparison. (Upper Panel) Ratio of the single nthorder energy contribution with respect to the full infinite-order MBD* binding energy. (Lower Panel) Deviation of the cumulative summation of the many-body contributions up to nth order with respect to the full infinite-order MBD* binding energy.

to the DWCNT, this system possesses only two graphene-like crossing linear bands29 at the Fermi surface. Hence, the contribution of these delocalized states to the overall correlation energy will be rather limited. To ensure convergence with respect to finite-size effects, we adopted a supercell of 35.7 Å along the DWCNT longitudinal axis, corresponding to 900 carbon atoms. The MBD* contribution to the binding energy is defined as previously described, in which the inner and outer nanotubes represent the two fragments of the host−guest complex. The apparent dimension of the system does not represent a limitation for the MBD* method, which can easily be applied to systems containing thousands of atoms. By virtue of the spatial extension and highly anisotropic character of the DWCNT, long-range manybody effects appear to be strongly enhanced in this complex: the long-range MBD* contribution to the binding energy is 853



Letter

semilocal density functional approximation and many-body correlation effects beyond the dipole approximation.

reduced by 24.7% from second to infinite order, as illustrated in Figure 3. In addition, a strikingly slow percentage-wise convergence of the perturbative MBD* series is also observed, leading to a deviation of 5.7% from full infinite-order MBD* at sixth-order, corresponding to a remarkable ∼32 kcal/mol per supercell. A comparable behavior is also expected in finiteextension nanotubes with length scales comparable to the present supercell, although the convergence of the many-body expansion was shown to be very slow in anisotropic lowdimensional systems.30 The sensitive dependence of the many-body effects on the structure of the host−guest complex and the nature of the binding strongly indicate that a simple renormalization of the standard C6/R6 pairwise summation is unable to provide the same degree of accuracy as a full many-body treatment of the long-range correlation energy in all types of supramolecular systems. In the same breath, a perturbative approach limited to a fixed finite-order also appears to be inaccurate due to the relatively slow convergence of the perturbative series. In particular, second-order Moeller−Plesset perturbation theory (MP2) completely neglects many-body dispersion effects, as they arise from higher-order correlation contributions. The ad hoc addition of the ATM term to empirical dispersioncorrected DFT approaches31,32 is questionable in this regard, further complicated by the fact that unphysical damping functions have to be utilized. In fact, employing different empirical damping functions for two- and three-atom interactions is inconsistent, as the attenuation of this fundamental interaction at any perturbative order stems from the same origin. In the MBD* approach, instead, the shortrange attenuation of the Coulomb interaction is seamlessly achieved, requiring only the attenuation of the interaction between two atoms at a time, which is consistent with the presence of terms up to two-particle interactions (at most) in the full nucleo-electronic Hamiltonian. By making use of state-of-the-art DQMC algorithms and large-scale computational resources, we provided benchmark binding energies for a set of six host−guest complexes from the S12L database. The DQMC data represent the first accurate benchmark for large supramolecular systems, with estimated errors not far from chemical accuracy. Very close agreement is found between DQMC and the PBE+MBD* method, and a perturbative many-body decomposition analysis of the longrange correlation energy in these host−guest complexes clearly demonstrated the need for an accurate description of manybody correlation effects. The influence of these many-body interactions was found to have a significant dependence on the symmetry and underlying topology of the host−guest complexes. As a consequence, these effects cannot be recovered by an effective pairwise approach as high perturbative orders are required to converge to the full infinite-order MBD* limit for the long-range correlation energy. Moreover, the mere inclusion of the three-body Axilrod−Teller−Muto energy contribution was shown to provide an overestimated reduction of the binding with respect to the full infinite-order energy and cannot be used to reproduce the infinite-order long-range correlation energy with high fidelity. The relatively successful application of the PBE+MBD* method to host−guest complexes of chemical relevance demonstrates that this is a promising approach for the challenging investigation of largescale supramolecular systems. The remaining issues to be addressed include a systematic analysis of the underlying

■

ASSOCIATED CONTENT

* Supporting Information S

Description of the MBD* computational method along with technical details regarding the DQMC calculations. Further analysis concerning the role of fractional exact exchange and a comparison among pairwise dispersion methods is also given for completeness. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.T. and A.A. received support from the European Research Council (ERC Starting Grant VDW-CMAT). This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy (DOE) under Contract No. DEAC05-00OR22725. R.A.D. received funding from the Department of Energy under Grant No. DE-SC0005180.

■

REFERENCES

(1) Pedersen, C. J. Cyclic Polyethers and Their Complexes with Metal Salts. J. Am. Chem. Soc. 1967, 89, 7017−7036. (2) Pedersen, C. J. The Discovery of Crown Ethers (Nobel Lecture). Angew. Chem., Int. Ed. Engl. 1988, 27, 1021−1027. (3) Cram, D. J. The Design of Molecular Hosts, Guests, and Their Complexes (Nobel Lecture). Angew. Chem., Int. Ed. Engl. 1988, 27, 1009−1020. (4) Lehn, J.-M. Supramolecular Chemistry−Scope and Perspectives Molecules, Supermolecules, and Molecular Devices (Nobel Lecture). Angew. Chem., Int. Ed. Engl. 1988, 27, 89−112. (5) Atwood, J. L.; Steed, J. W. Supramolecular Chemistry; Wiley: Weinheim, Germany, 2000. (6) Lehn, J.-M. Supramolecular Chemistry. Science 1993, 260, 1762− 1763. (7) Zhang, S. G. Fabrication of Novel Biomaterials through Molecular Self-Assembly. Nat. Biotechnol. 2003, 21, 1171−1178. (8) Balzani, V.; Gomez-Lopez, M.; Stoddart, J. F. Molecular Machines. Acc. Chem. Res. 1998, 31, 405−414. (9) Waller, M. P.; Kruse, H.; Mueck-Lichtenfeld, C.; Grimme, S. Investigating Inclusion Complexes Using Quantum Chemical Methods. Chem. Soc. Rev. 2012, 41, 3119−3128. (10) Rezac, 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. (11) Jurecka, P.; Sponer, J.; Cerny, 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. (12) Rezac, J.; Riley, K. E.; Hobza, P. Benchmark Calculations of Noncovalent Interactions of Halogenated Molecules. J. Chem. Theory Comput. 2012, 8, 4285−4292. (13) Booth, G. H.; Grüneis, A.; Kresse, G.; Alavi, A. Towards an Exact Description of Electronic Wavefunctions in Real Solids. Nature 2012, 493, 365−370. (14) Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, J. Quantum Monte Carlo Simulations of Solids. Rev. Mod. Phys. 2001, 73, 33−83. 854



Letter

(15) Diedrich, C.; Lüchow, A.; Grimme, S. Weak Intermolecular Interactions Calculated with Diffusion Monte Carlo. J. Chem. Phys. 2005, 123, 184106. (16) Dubecký, M.; Jurečka, P.; Derian, R.; Hobza, P.; Otyepka, M.; Mitas, L. Quantum Monte Carlo Methods Describe Noncovalent Interactions with Subchemical Accuracy. J. Chem. Theory Comput. 2013, 9, 4287−4292. (17) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem.Eur. J. 2012, 18, 9955−9964. (18) Risthaus, T.; Grimme, S. Benchmarking of London DispersionAccounting Density Functional Theory Methods on Very Large Molecular Complexes. J. Chem. Theory Comput. 2013, 9, 1580−1591. (19) Needs, R. J.; Towler, M. D.; Drummond, N. D.; López Ríos, P. Continuum Variational and Diffusion Quantum Monte Carlo Calculations. J. Phys.: Condens. Matter 2010, 22, 023201−023215. (20) Tkatchenko, A.; Ambrosetti, A.; DiStasio, R. A., Jr. Interatomic Methods for the Dispersion Energy Derived from the Adiabatic Connection Fluctuation-Dissipation Theorem. J. Chem. Phys. 2013, 138, 074106. (21) Tkatchenko, A.; DiStasio, R. A., Jr.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body van der Waals Interactions. Phys. Rev. Lett. 2012, 108, 236402. (22) Ambrosetti, A.; Reilly, A. M.; DiStasio, R. A., Jr.; Tkatchenko, A. Long-range correlation energy calculated from coupled atomic response functions. J. Chem. Phys. 2014, accepted for publication. (23) Reilly, A. M.; Tkatchenko, A. Seamless and Accurate Modeling of Organic Molecular Materials. J. Phys. Chem. Lett. 2013, 4, 1028− 1033. (24) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (25) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. Comput. Phys. Commun. 2009, 180, 2175−2196. (26) Tkatchenko, A.; Alfè, D.; Kim, K. S. First-Principles Modeling of Non-Covalent Interactions in Supramolecular Systems: The Role of Many-Body Effects. J. Chem. Theory Comput. 2012, 8, 4317−4322. (27) DiStasio, R. A., Jr.; von Lilienfeld, O. A.; Tkatchenko, A. Collective Many-Body van der Waals Interactions in Molecular Systems. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 14791−14795. (28) Axilrod, B. M.; Teller, E. Interaction of The van der Waals Type Between Three Atoms. J. Chem. Phys. 1943, 11, 299−300. (29) Kwon, Y.-K.; Tománek, D. Electronic and Structural Properties of Multiwall Carbon Nanotubes. Phys. Rev. B 1998, 58, 16001(R). (30) Gobre, V. V.; Tkatchenko, A. Scaling Laws for van der Waals Interactions in Nanostructured Materials. Nat. Commun. 2013, 4, 2341. (31) 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. (32) Otero-de-la-Roza, A.; Johnson, E. R. Many-Body Dispersion Interactions from the Exchange-Hole Dipole Moment Model. J. Chem. Phys. 2013, 138, 054103.

855


Hard Numbers for Large Molecules: Toward Exact Energetics for

Recommend Documents