Accurate Intermolecular Potential for the C60 Dimer: The Performance

Article pubs.acs.org/JCTC

Accurate Intermolecular Potential for the C60 Dimer: The Performance of Different Levels of Quantum Theory Dmitry I. Sharapa,† Johannes T. Margraf,‡ Andreas Hesselmann,§ and Timothy Clark*,† †

Computer-Chemie-Centrum, Department Chemie und Pharmazie, Friedrich-Alexander-Universität Erlangen-Nürnberg, Nägelsbachstraße 25, 91052 Erlangen, Germany ‡ Quantum Theory Project, University of Florida, Gainesville, Florida 32611, United States § Lehrstuhl für Physikalische und Theoretische Chemie, Friedrich-Alexander-Universität Erlangen-Nürnberg, Egerlandstraße 3, 91058 Erlangen, Germany S Supporting Information *

ABSTRACT: The self-assembly of molecular building blocks is a promising route to low-cost nanoelectronic devices. It would be very appealing to use computer-aided design to identify suitable molecules. However, molecular self-assembly is guided by weak interactions, such as dispersion, which have long been notoriously difficult to describe with quantum chemical methods. In recent years, several viable techniques have emerged, ranging from empirical dispersion corrections for DFT to fast perturbation and coupled-cluster theories. In this work, we test these methods for the dimer of the prototypical building block for nanoelectronics, C60-fullerene. Benchmark quality data is obtained from DFT-based symmetryadapted perturbation theory (SAPT), the adiabatic-connection fluctuation dissipation (ACFD) theorem using an adiabatic LDA kernel, and domain-based local pair natural orbital (DLPNO) coupled-pair and coupled-cluster methods. These benchmarks are used to evaluate economical dispersion-corrected DFT methods, double-hybrid DFT functionals, and second-order Møller−Plesset theory. Furthermore, we provide analytical fits to the benchmark interaction curves, which can be used for a coarse-grain description of fullerene self-assembly. These analytical expressions differ significantly from those reported previously based on bulk data.

■

large for standard DFT calculations, so that classical force fields are still the most common methods used for computational life science. In this case, dispersion interactions are commonly described using Lennard-Jones24 or Buckingham25 potentials, which are perfectly adequate in the context of protein- and nucleic-acid force fields, although more specialized systems may require, for instance, anisotropic repulsion terms.26 The situation in the field of organic electronics is different; it is essential to determine the electronic structure of the system. For example, describing electron transport through self-assembled monolayers of functionalized fullerenes27−30 necessarily requires knowledge about the position and depth of the van der Waals minima between both neutral and charged fullerenes. Describing the intermolecular binding energy in a negatively charged C60-fullerene dimer represents a very challenging computational task because many techniques that are suitable for the problem cannot handle 120 carbon atoms.31 Ruzsinszky et al.32 pointed out that dispersion interactions between large N-atom fullerenes do not scale with N2, as would be the case if pairwise potentials were appropriate, but rather

INTRODUCTION In the past decade, quantum chemical methods for describing noncovalent (and in particular dispersion) interactions quantitatively have improved strongly.1−7 This is particularly true for density-functional theory (DFT)-based methods, since they often represent the best compromise between accuracy and cost but were originally unable to describe dispersion. Notable examples are symmetry-adapted perturbation theory based on DFT: DFT-SAPT,8−12 the exchange-hole dipole moment (XDM)approach,13,14 the dDsC density-dependent energy correction,15 the nonlocal van der Waals (NL-vdW) density functional in the most recent form derived by Vydrov and van Voorhis (VV10),16,17 and Grimme’s atom-pairwise D318−20 corrections. These methods allow fairly large systems of practical interest to be calculated.21 Their performance is typically evaluated using benchmark sets of molecular dimerization energies (e.g., S2222 and S12L23). DFT functionals with Grimme’s D3-correction have emerged as the most popular tools for describing noncovalent interactions because of the low cost of calculating the correction, reasonable agreement with experimental and ab initio results, and their wide availability in quantum-chemical software. Noncovalent interactions also play an important role in all biological processes, but biological systems are generally too © 2016 American Chemical Society

Received: September 1, 2016 Published: November 29, 2016 274

DOI: 10.1021/acs.jctc.6b00869 J. Chem. Theory Comput. 2017, 13, 274−285

Article

Journal of Chemical Theory and Computation with approximately N2.75, due to the increasing polarizability of the asymptotically (N → ∞) metallic fullerenes. Reilly and Tkatchenko33 have stressed that many-body (more precisely: many-atom) contributions to the correlation energy that are absent in standard pairwise dispersion corrections to DFT functionals are needed in many cases for molecular materials, including fullerenes, and have reviewed the techniques available. Such contributions were studied in detail for a set of 12 supramolecular complexes (the S12L set) containing systems up to about 800 electrons.21,34 Risthaus and Grimme21 observed that three-atom contributions to the interaction energy vary between +0.7 and +4.6 kcal mol−1 for the different complexes and therefore lead to a reduction of the interaction energy for all systems considered. Ambrosetti et al.34 have shown that accounting for many-atom contributions to all orders changes the second-order interaction energies (which only contain pairwise dispersion contributions) by up to 10% for the various complexes. These findings raise the question as to whether second-order approaches such as pairwise dispersion-corrected DFT methods or MP2 and SCS-MP2 are appropriate for describing the interaction energy of fullerene dimers. In this regard, note that the structure of the fullerene dimer is clearly different from the host−guest complexes of the S12L set because there is only one interaction site represented by two stacked C6 rings (see Figure 2), in contrast to the supramolecular complexes in which a (usually small) guest molecule is surrounded by a large host. The DFT-SAPT interaction energy decomposition presented below shows that the character of the interaction between the two C60 atoms is very similar to that of a stacked benzene dimer with the carbon atoms positioned at the same locations as the interaction site of the minimum structure of (C60)2. This indicates that many-atom contributions to the interaction energy do not play a dominant role in the stabilization of the dimer; otherwise a stronger variation of the percentage contribution of the dispersion energy to the total interaction energy would be expected in the two cases. Moreover, the three-atom Axilrod−Teller35 contribution to the dispersion energy as calculated using Grimme’s D3 program18 only amounts to +0.24 kcal mol−1 to the interaction energy of (C60)2 at equilibrium. This also indicates that manyatom terms only give small corrections to the interaction energy. Note, on the other hand, that early studies (from 1992) of the low-temperature (5K) crystal structure of C60-fullerene have shown that Lennard-Jones (two-center) potentials cannot reproduce the experimental structure,36 which indicates that many-body contributions may become important for reproducing the bulk properties of fullerene. The literature abounds with studies closely related to this one. We give a brief summary here: • π−π stacking has been discussed for coronene,37,38 corannulene, and sumanene.39−41 B3LYP+D3 was found to agree well with CCSD(T)38 and QCISD(T)37 benchmarks for the coronene dimer. B97-D was found to reproduce QCISD(T) and CCSD(T) results for corannulene well.42 A comparison between DFT functionals for planar and bowlshaped polyacenes is given in ref 43. The M06-2X binding energy for coronene obtained in this study with a cc-pVTZ basis set differs by a factor of 2.5 from results obtained in another study with the 6-31+G(d,p) basis set.44 • Large supramolecular test systems are collected in Grimme’s S12L21 and S30L45 data sets. The two complexes of “buckycatchers” with fullerenes are most relevant for the

current study (as examples of systems with dominant nonpolar interactions) but also technically the most difficult. Thus, DLPNO-CCSD(T)/def2-TZVPP calculations successfully performed on other systems did not converge for these cases.4 Furthermore, for these two systems, a significant deviation of back-corrected empirical binding energies from values obtained by DFT-SAPT, MP2C, and NLDFT placed doubt on the accuracy of Grimme’s data.46 However, Grimme suggests PW6B95-D3/def2-QZVP as the benchmark combination of the functional and basis set.45 • The relevance of crystal structures is limited in our case; the effect of crystal structure should be taken into account when comparing data derived from fullerite parameters with gas-phase values. Tkatchenko showed that PBE performs well for the electronic properties of polyacenes, both in the gas and crystalline phases.47 An even more impressive demonstration of the reliability of dispersion-corrected PBE for modeling organic crystals was described for β-hematin.48 The interaction between neutral fullerenes has been studied many times, with differing results.49,50 Theoretical studies typically use the experimentally derived analytical potential of Girifalco51,52 as reference. However, experimental uncertainty and the neglect of effects such as the deformation and orientation of the fullerenes cast doubt on the accuracy of this potential, which uses spherical fullerenes as a model to derive a two-center interatomic potential. Two different studies53,54 have estimated that the difference between the maximum and rotationally averaged binding energies between two C60 molecules is more than 2 kcal mol−1. The Girifalco potential corresponds to a rotationally averaged interfullerene interaction and to a cubic lattice of fullerenes, rather than to a simple fullerene dimer. The aim of this work is to find appropriate methods to study interfullerene interactions and provide reliable benchmark numbers using the neutral C60-fullerene dimer as a test case. We have used the most accurate methods that we can apply to this problem, approximate coupled-cluster methods and intermolecular perturbation theory using the DFT-SAPT approach. The performance of more economical models is compared with these standards.

■

COMPUTATIONAL DETAILS All DFT, second-order Møller−Plesset (MP2), and DLPNO (domain based local pair-natural orbital) calculations were performed with ORCA.55 The LCB-OP-LRD (long-range exchange corrected BOP functional with local response dispersion) calculations were performed in Gamess.56,57 DFT-SAPT (symmetryadapted perturbation theory using a DFT monomer description)9,46 calculations were performed with MOLPRO.58 The DLPNO methods DLPNO-CCSD,59 DLPNO-CCSD(T)60 (referred to in the following as DLPNO-CC), and the significantly more efficient DLPNO-CEPA/1 methods were also used. The DLPNO-CC implementation in the currently publicly available version of ORCA (3.0.3) that was used in this work is not based on the concept of sparse maps.61 The TZV62,63 and cc-pVDZ64 basis sets were used to calculate intermolecular potential-energy curves with DLPNO-CC methods. The larger cc-pVTZ64 and def2-TZVP65 basis sets were used to calculate the binding energy at the equilibrium distance and at a distance of 14.17 Å (corresponding to the lattice cell parameter respectively the second-nearest neighbor in the bulk). The hierarchical Dunning basis sets64 were used to perform two-point extrapolations to the basis set limit (CBS) 275


Article

Journal of Chemical Theory and Computation

Figure 1. Interfullerene potential calculated with different approaches. Total energies are given in the Supporting Information.

corresponding distance (i.e., the lattice cell parameter divided by √2 for this crystal structure) is 10.0281−10.0836 Å in the orientationally disordered face-centered cubic phase at room temperature and 9.928 Å in the orientationally ordered lowtemperature simple cubic phase.36 The most stable orientation of the fullerene dimer as obtained from the BP86-D3 method is shown in Figure 2. The complex

with the DLPNO-CC values. The cc-pVDZ, cc-pVTZ, and def2-TZVP basis sets were used to calculate intermolecular potential-energy curves with DFT methods. We used auxiliary basis sets for correlation calculations (/C66−68) for the DLPNO part but without the resolution of the identity (RI)-approximation to accelerate the selfconsistent field (SCF) calculations. RI-MP2,69 RI-SCS-MP2,70 and RI-SOS-MP271 calculations used the/C-auxiliary basis set. DFT calculations were performed with the RIJCOSX (“chain-ofspheres exchange”) approximation72,73 and corresponding/J74 auxiliary basis set. Core electrons were kept frozen in all wave function correlation energy calculations. The Boys−Bernardi counterpoise correction75 was used to reduce the basis set superposition error in all ACFD(ALDA) calculations. All DFTSAPT calculations were performed in the full dimer basis set to reduce the basis set incompleteness error of the individual interaction energy terms.

■

RESULTS AND DISCUSSION Geometry. To the best of our knowledge, no experimental data is available for the geometry of the gas-phase C60-dimer. We therefore optimized initial geometries for the neutral C60fullerene monomer and dimer with the BP86 functional,76,77 TZV basis set,62,63 and D3BJ dispersion correction.18−20 Based on these geometries, rigid scans of the intermolecular distance were performed using different methods (in analogy to the methodology used in Hobza’s S66×8 data set).78−80 Energy profiles obtained from these scans are shown in Figure 1. DLPNO-CCSD/cc-pVDZ predicts an equilibrium distance for the dimer minimum of 9.748 Å, compared with 9.648 Å at the SOS-MP2/cc-pVDZ level. DFT-based methods give equilibrium distances in the range of 9.75−9.85 Å. Since the interfullerene interaction is mainly described by the dispersion correction, the influence of the functional and basis set is quite small for the DFT calculations. MP2/cc-pVDZ predicts a much shorter equilibrium distance that differs markedly from the results given by other methods. Given the good agreement between DFT and DLPNO-CC methods, we take 9.75 Å as the equilibrium distance of the gas-phase C60-dimer. Note that the

Figure 2. BP86-D3/TZV optimized geometry of the most stable orientation of the C60-dimer found.

exhibits C2‑symmetry and the interface can be described as “parallel, shifted 6-membered rings”. Several other possible orientations were optimized using DFT-D3 (fixed orientation, relaxed distances, and deformations). DLPNO-CEPA/1 calculations with the cc-pVDZ basis set showed that these alternative orientations (6-membered ring to 6-membered ring or bond-tobond) are significantly less favorable (up to 4 kcal mol−1). The calculated orientation differs strongly from the “pentagonal faces facing interpentagon bonds” orientation found in the lowest-temperature bulk crystalline phase,36 which we find to be approximately 0.7 kcal mol−1 less stable than our global minimum. Optimized geometries at the different interfullerene distances are given in the Supporting Information. Basis Set Influence. It is well-known that high angular momentum basis functions are necessary to describe electron correlation effects at short interelectronic distances in coupledcluster methods correctly.82 Indeed, DLPNO-CC calculations 276


Article


properties) allows the total intermolecular interaction energy to be decomposed into a sum of physically interpretable terms,8,9,11,12 as shown in eq 1

with the TZV basis set, 11s6p [62111/411], strongly underestimate the binding energy (by a factor of more than three in comparison to experimental estimates (−6.34 ± 1.84 kcal mol−1)83). Also, removing the f-functions from the def2-TZVP basis set (def2-TZVP(-f), 11s6p2d [62111/411/11]), reduces the binding energy by approximately 6.5 kcal mol−1 (in comparison to employing the unmodified def2-TZVP basis set). A similar strong underestimation of the binding energy with the TZV basis set (in comparison to Dunning basis sets) was also found at the MP2 and SOS- and SCS-MP2 levels. In an independent test, we also found that post-HF methods with the TZV basis set give an anthracene dimerization energy half as large as with more flexible basis sets. However, the TZV basis set is acceptable for standard DFT calculations based on the generalized gradient approximation (GGA) because these methods do not model the interelectronic cusp explicitly. Furthermore, note that standard pairwise dispersion corrections are independent of the basis set. DLPNO-CC. A series of test calculations was performed in order to determine the effects of the basis set, the DLPNO method, and cutoffs on the calculated dimerization energy. Details of the individual results are given in the Supporting Information. We outline the major findings only here. Generally, DLPNO-CEPA/1 gives results within 0.5 kcal mol−1 of DLPNO-CCSD with any given basis set. The accuracy of DLPNO-CCSD(T) for noncovalent interactions has, for example, been demonstrated in ref 84. Because the three DLPNO-CC methods agreed closely, the least computationally expensive DLPNO-CEPA/1 method was used for the largest calculations. All discussions in this section refer to DLPNO-CEPA/1 unless otherwise stated. The cc-pVDZ basis set is not sufficient for describing the interfullerene interaction; Liakos has shown for the S66 data set78,79 that the CCSD(T) and LPNO-CEPA/1 methods only recover 70−75% of the CBS correlation energy with the cc-pVDZ basis set.85 To improve the theoretical description, we tried enlarging the atom-centered basis set, adding bondcenter basis functions at the geometrical center of the dimer and tightening the PNO cutoffs. The first two approaches result in a significantly higher (more negative) binding energy, while tightening the PNO thresholds has more complex effects due to significant changes in the counterpoise correction. Adding bond-center functions is the most economical way to improve the basis set, whereas tightening the PNO cutoffs makes calculations significantly more expensive. We have proposed an approximate scheme to reduce the error caused by the low cutoffs in the PNO calculations. These results can be found in the Supporting Information. To obtain benchmark data, DLPNO-CEPA/1 calculations with the cc-pVTZ basis set and TightPNO cutoffs were performed. The counterpoise-corrected binding energy extrapolated to a complete basis set, free of any scaling coefficients, is −8.1 kcal mol−1 at a center-to-center distance of 9.75 Å. We found that using TightPNO cutoffs with the cc-pVDZ basis set gives the position of minimum at a slightly larger center-to-center distance. The difference between binding energies (cc-pVDZ basis set) at 9.75 and 9.85 Å is smaller than 0.2 kcal mol−1. Thus, we can estimate the minimum binding energy of (C60)2 to be −8.3 kcal mol−1 with an uncertainty of ±0.1 kcal mol−1. DFT-SAPT. The DFT-SAPT technique (symmetry-adapted perturbation theory employing intramolecular correlation corrections from static and time-dependent DFT monomer

(1) (1) (2) (2) (2) (2) Eint = Eelst + Eexch + Eind + Eexch − ind + Edisp + Eexch − disp + δ(HF )

(1) (2) where E(1) elst denotes the electrostatic interaction energy, Eind denotes the induction energy, E(2) denotes the dispersion disp (2) (2) energy, and the terms E(1) exch and Eexch−ind and Eexch−disp are the corresponding short-range exchange-repulsion contributions, which stem from a tunneling of the electrons between the two monomers. As indicated in eq 1, electron correlation terms higher than second-order are not accounted for in the standard DFT-SAPT method. The question arises as to what extent the second-order dispersion energy is able to capture many-atom correlation interactions that are absent in pairwise DFT+D dispersion corrections or in second-order correlation methods such as MP2. If the dispersion energy is calculated using infinite-order response functions from TDDFT calculations, as in DFTSAPT, it can be shown that the dispersion energy will be complete through all three-atom interaction energy terms (i.e, captures the Axilrod−Teller interaction energy contribution),35 while it misses those fourth- and higher-order terms in which the intermolecular interaction is of third and higher orders. Note, of course, that this also means that E(2) disp cannot describe simultaneous electron correlation interactions between three or more molecules. The accuracy of the individual interaction-energy terms can depend critically on the shape of the underlying exchangecorrelation (xc) potentials used in the DFT monomer calculations, in particular their asymptotic behavior.8,9 In this work, the PBE0AC xc potential, which uses the Becke-RousselJohnson exchange potential to approximate the exact exchange potential,86,87 was used to calculate the Kohn−Sham orbitals and orbital energies. All other technical details of the calculations are described in ref 46. The first-order terms and the second-order (exchange-) induction energies converge quickly toward the basis set limit upon increasing the basis set size. They can generally be estimated well with nonaugmented triple-ζ basis sets, see ref 46. However, the intermolecular correlation contributions (i.e., the (2) E(2) disp and Eexch−dispenergies) depend strongly on the quality of the basis set, especially on the diffuseness of the basis functions. Since, however, on the other hand, diffuse functions can lead to poor convergence of the SCF calculations for large systems due to almost linearly dependent basis functions, a modified basis-set extrapolation technique for DFT-SAPT calculations was introduced,46 shown in eq 2

∞ VDZ VTZ E int (corr ) ≈ α[E int (corr ) → E int (corr )]

(2)

where corr here denotes the correlation contribution to the VTZ interaction energy and [EVDZ int (corr) → Eint (corr)] is the extrapolated correlation energy using the vdz → vtz two-point extrapolation scheme. VDZ and VTZ denote the (nonaugmented) double and triple-ξ Dunning basis sets64 and α is a constant factor (>1) fitted46 to basis-set converged dispersion and exchange-dispersion interaction-energy terms. It has been found46 that an optimized value of α (1.08) reduces the deviation between the modified extrapolation scheme of the above equation and aug-cc-pVTZ → aug-cc-pVQZ extrapolated results to approximately ±5% for the S22 dimers. The results in 277


Article

Journal of Chemical Theory and Computation ref 46 showed, however, that this value of α can lead to an overestimation of the interaction energies for dispersion-bound complexes. Because of this, the above recipe of eq 2 was tested for a substructure of the full (C60)2 dimer for which large aug-cc-pVTZ basis set calculations can still be performed. For this, only the interaction between the two π-stacked benzene rings (see Figure 2) was considered. Figure 3 shows the interaction potential obtained from the cc-pVDZ → cc-pVTZ and the aug-cc-pVTZ → aug-cc-pVQZ

interaction energy terms to the total interaction energy of the benzene and C60 dimers at their equilibrium distances (3.87 and 9.84 Å, respectively). As can be seen in the figure, the interaction energy decomposition is very similar for the two dimers, and, clearly, in both cases the dispersion interaction energy yields the strongest contribution to the total interaction energy. A notable difference between the two systems revealed by Figure 4 is that, in the case of (C60)2, both the contribution of the first-order exchange interaction energy and that of the dispersion energy are reduced by about 50% as compared to the benzene dimer. This can clearly be attributed to the larger distance (of the centers of mass) of the C60 dimer and shows that the minimum-energy structure of the fullerene dimer is determined by a delicate balance between the strongly repulsive exchange interaction energy and the dispersion interaction energy, see also below. Apart from this distance-dependent feature, however, the characters of the interactions in the full C60 dimer and the model benzene-dimer substructure are very similar. We can therefore expect that the basis-set dependences of electron correlation contributions to the interaction energy also behave similarly for the two dimers. Figure 5 shows the individual interaction-energy contributions obtained from the DFT-SAPT calculations along with the

Figure 3. Benzene dimer: basis set extrapolations of the DFT-SAPT interaction energy.

extrapolations. As can be seen, at the equilibrium distance of the two stacked benzene rings the cc-pVDZ → cc-pVTZ extrapolated interaction energy underestimates the interaction energy from the two-point extrapolation using the augmented double and triple-ζ basis sets by almost 0.4 kcal mol−1 (12% of the total interaction energy). A scaling of the cc-pVDZ → cc-pVTZ extrapolated results by a factor of α = 1.06 (red curve in Figure 3) corrects this error, while the value of α = 1.08 proposed in ref 46 would lead to an overestimation of the interaction energy of about 0.13 kcal mol−1 at the geometry of the minimum. Because of this, in this work the scaling factor of α = 1.06 was used in eq 2 to extrapolate the DFT-SAPT correlation contributions. In order to justify the extrapolation formula described above, Figure 4 shows the percentage contributions of the individual

Figure 5. Interaction energy contributions and total intermolecular interaction energy (Eint) calculated using the DFT-SAPT method.

total interaction potential (the dark green line). It shows that the minimum of the interaction energy potential is located at about 9.8 Å and has a depth of −8.0 kcal mol−1. The attractive first-order electrostatic interaction energy is completely quenched by the first-order exchange contribution, which is by far the most repulsive interaction-energy term at short distances between the monomers. At second order, a fairly strong induction energy is also observed, almost twice as large as the electrostatic interaction at a distance of 9 Å, see Figure 5. Again, however, this contribution almost completely cancels with its exchange-induction counterpart, so that in total one can argue that polarization effects do not play a role in the stabilization of the dimer. This is also supported by the very weak contribution of the δ(HF) energy, which sums higherorder polarization contributions on the supermolecular Hartree−Fock level. Because of this, as expected, the main attractive contribution to the total interaction energy of the C60 dimer originates from the dispersion interaction, which at the

Figure 4. Percentage contributions of individual interaction energy terms yielded by the DFT-SAPT method to the total interaction energy for the benzene and C60 dimer. 278


Article

Journal of Chemical Theory and Computation equilibrium distance is more than twice as large as the total interaction energy (Figure 5). The exchange-dispersion interaction energy is far smaller at short intermolecular distances and quickly decays to zero at larger distances. It should be noted, however, that the E(2) exch−disp term cannot be neglected at intermediate distances because it amounts to almost +4 kcal mol−1 at the equilibrium distance. ACFD(ALDA). The adiabatic-connection fluctuation dissipation theorem (ACFD)88,89 provides a route to obtain electron correlation energies from the Kohn−Sham response functions and the generally nonadiabatic and nonlocal exchangecorrelation (xc) kernel. Provided the latter were known exactly, the ACFD method would yield the exact correlation energy of a many-electron system. In practice, however, the xc kernel must be approximated. The most crucial approximation is to neglect any xc contributions to the kernel completely, which results in the random-phase approximation (RPA).90 It is known, however, that the RPA method underestimates the interaction energies of dispersion-dominated complexes,91 and thus this method will not yield accurate results for the interaction energy of (C60)2. In this work the xc kernel used within the ACFD approach has been approximated on the adiabatic local density approximation (ALDA) level using the approach described in ref 92 in order to resolve the singularities of the pair correlation function at short interelectronic distances, see also refs 93 and 94 for two related ACFD approaches. In ref 92 this method has also been tested for its performance for intermolecular interactions using the S22 benchmark set by Hobza et al., and it was found that it reproduces the coupledcluster reference values of this data set with a fairly high accuracy yielding an absolute error of only 0.12 kcal mol−1 on average for the 22 systems. The ACFD(ALDA) method scales formally with N4 with respect to the molecular size only and therefore can even be applied to extended systems using large basis sets. However, as noted above in the DFT-SAPT discussion, the use of extended basis sets containing diffuse functions becomes technically problematic for the C60 molecule due to the almost linearly dependent basis functions. Because of this, as for DFT-SAPT, the complete basis set limit for the ACFD(ALDA) interaction energy was estimated using eq 2. It should be taken into account, however, that now the scaling applies to the complete correlation contribution to the interaction energy (containing also intramolecular correlation contributions), and it can thus be expected that the optimal scaling factor of α will be different for the DFT-SAPT method and generally supermolecular approaches. Figure 6 again shows the two-point extrapolated interaction energy curves for the benzene dimer substructure of (C60)2 (see DFT-SAPT section above) using the cc-pVDZ → ccpVTZ and the aug-cc-pVDZ → aug-cc-pVTZ extrapolation schemes. Compared to the DFT-SAPT method, see Figure 3, now the cc-pVDZ → cc-pVTZ extrapolated interaction energy at the minimum only underestimates the extrapolated result that takes into account diffuse basis functions by 0.16 kcal mol−1. The correlation contribution obtained from the cc-pVDZ → cc-pVTZ here therefore needs to be scaled by a factor of α that is smaller than the optimized value of α = 1.06 found for the DFT-SAPT method. A minimization of the average percentage deviation between the aug-cc-pVDZ → aug-cc-pVTZ extrapolated results and eq 2 for the calculated points of the benzene dimer potential yield α = 1.03. With this, the scaled extrapolated potential reproduces the potential from the

Figure 6. Benzene dimer: basis set extrapolations of the ACFD(ALDA) interaction energy.

aug-cc-pVDZ → aug-cc-pVTZ extrapolation with a reasonable accuracy around the equilibrium; see the red curve in Figure 6. The counterpoise-corrected interaction energy potential for (C60)2 of the ACFD(ALDA) method using the corrected basis set extrapolation technique of eq 2 is presented in Figure 7

Figure 7. C60 dimer interaction energy potential for the DFT-SAPT and ACFD(ALDA) methods obtained from the scaled complete basis set extrapolation technique. Total energies are given in the Supporting Information.

along with the extrapolated DFT-SAPT potential. As can be seen in the figure, the two curves are very close to each other for all distances plotted in the diagram. Both methods yield a minimum at a distance of 9.85 Å with depths of −8.1 (ACFD) and −8.0 kcal mol−1 (DFT-SAPT). Since the interaction energy is determined quite differently in the two methods, this result indicates strongly that the predicted minimum obtained by the two methods is close to the true one. Moreover, the small differences between the extrapolated DFT-SAPT and ACFD(ALDA) interaction energies also indicate that third- and higher-order intermolecular electron correlation contributions, which are not accounted for in the DFT-SAPT method (see above), turn out to be negligible. These estimates also agree with the result obtained using DLPNO-CEPA/1 with TightPNO cut-offs and CBS extrapolation (−8.1 kcal mol−1 at 9.75 Å, see above). Finally, these values lie at the upper bound of experimental interval suggested by Branz for the gas-phase dimer (4.5−8.2 kcal mol−1).83 279


Article


Three different types of dispersion correction were used: atompairwise dispersion corrections with Becke-Johnson damping (D3BJ),18−20 local response dispersion (LRD),114,115 and the nonlocal van der Waals density functional (NL-vdW DF) in the form by Vydrov and Van Voorhis (VV10),16,17 with and without self-consistent (SC) treatment. The VV10 and LRD methods116 allow the dispersion energy to be calculated from the ground-state electron density only, providing an alternative to Grimme’s D3-correction. A direct comparison (performed with the B3LYP functional) shows that the difference between results calculated with the D3BJ and SCNL corrections is less than 1 kcal mol−1 and thus smaller than the effect of the basis set (up to 2.5 kcal mol−1). Small differences between the results obtained with the D3BJ and NL-vdW corrections combined with different DFT functionals have been pointed out by Grimme.99 Overall, the binding energies calculated with the dispersioncorrected DFT-functionals range from −6 to −12 kcal mol−1. The effect of three-body corrections was investigated for the BP86, M06-2X, and B3LYP functionals (all with D3BJ correction) with the cc-pVDZ, cc-pVTZ, and def2-TZVP basis sets. In all cases binding energy increases by only 0.24− 0.25 kcal mol−1 Of the double-hybrid functionals, canonical B2PLYP and mPW2LYP with the D3BJ correction give results between −11 and −12.5 kcal mol−1 and PWPB95 between −8.8 and −10.8 kcal mol−1. The DSD- and DOD-PBEP86 functionals with the SCNL correction116 give values between 0 and 5 kcal mol−1. A similar strong lowering of the binding energy was observed for the spin-scaled MP2 methods (see below). These results deviate strongly from those found with other functionals and cast doubt on the suitability of DSD- and DOD-PBEP86-SCNL for interfullerene and similar interactions containing π-stacking motifs. The canonical, spin-component-scaled and spin-oppositescaled MP2 techniques give differing results. Canonical MP2 consistently gives the lowest binding energies (−22.5 to −23.5 kcal mol−1), while the SCS-MP2 results lie between −16.5 and −17.5 kcal mol−1 and the SOS-MP2 results lie between −13.5 and −14.5. The latter two methods give an interfullerene distance close to the reference one, whereas canonical MP2 gives a far shorter contact (less than 9.5 Å) which can be ascribed to the strong overestimation of dispersion interaction by MP2. In contrast to DLPNO-CC and DFT-SAPT, the difference between cc-pVDZ and cc-pVTZ basis sets is less than 0.5 kcal mol−1 for all three MP2 methods. Interaction energies for dispersion corrected DFT scatter less strongly than for MP2-derived methods but still depend significantly on the basis set and functional. Recommended Analytical Potential. Figure 7 shows that DLPNO-CC, DFT-SAPT, and ACFD methods agree well on the position of the minimum and the binding energy. Dispersion-corrected DFT methods also reproduce this interfullerene distance, and the empirical dispersion corrections can reproduce the binding energy reasonably well. However, these curves differ significantly from the published analytical potentials because these neglected deformational and orientational effects and were fitted to bulk data obtained from crystals. We fitted the Girifalco51 (eq 3), Guerin117 (eq 4), Smith-Thakkar118 (eq 5), and Lim119 (eq 6) potentials to the DFT-SAPT/CBS data points (Figure 9). The fitting coefficients together with the squared correlation coefficient are shown in Table 2. The form of the potential used is given by

Table 1. Density Functionals Used in This Work functional

ref

Functional and Dispersion Correction LDA 96 Generalized Gradient Approximation Functionals (GGAs) BP86-D3BJ 76, 77 rPW86PBE-SCNL 97−99 Meta-GGAs TPSS-D3BJ 100−102 Hybrid and Long-Range Corrected Functionals B3LYP-D3BJ or −NL 103−106 LC-BOP-LRD 107 CAM-B3LYP-D3BJ 108 wB97X-D3 109, 110 M06-2X-D3BJ 111

DFT and MP2. Nine different DFT functionals from different rungs of Jacob’s ladder95 were tested with and without dispersion corrections (Table 1). Double hybrid functionals and MP2 techniques will be discussed separately below. Functionals without dispersion corrections gave a wide range of the interaction energies from clearly repulsive to slightly attractive, as shown in Figure 8. As might have been expected from its performance for graphite,112,113 LDA without dispersion-correction gives the correct form of the interaction potential (and position of the minimum), with a binding energy of around −6 kcal mol−1.

Figure 8. Schematic representation of the C60-dimerization energies obtained at different levels of theory. The benchmark value (−8.2 kcal mol−1) shown as a black line for comparison. Total and relative energies are given in the Supporting Information. 280


Article

Journal of Chemical Theory and Computation m ⎧ ⎡m⎛ ⎪ n ⎜⎛ L ⎟⎞ 2 r ⎞⎤ U (r ) = D ⎨ exp⎢ ⎜1 − ⎟⎥ ⎪ ⎣2⎝ L ⎠⎦ ⎩m − n⎝ r ⎠

−

n ⎫ ⎡n⎛ m ⎜⎛ L ⎟⎞ 2 r ⎞⎤⎪ ⎬ exp⎢ ⎜1 − ⎟⎥⎪ ⎣2⎝ m − n⎝ r ⎠ L ⎠⎦⎭

(6)

We emphasize that the calculated potential corresponds to the electronic energy of two optimally oriented molecules in the gas phase. These results differ from crystal calculations or empirical potentials built based on data for bulk phases.

■

CONCLUSIONS In this work, we have studied the noncovalent binding of the C60 dimer using various quantum chemical methods. It has been found that high-level coupled-cluster and coupled electron pair wave function correlation methods, which are known to provide accurate interaction energies for small dimer systems, become impractical for studying systems of this size (containing about 700 electrons) due to the extensive computer resources required. This holds even true for the recently developed DLPNO variants of these methods, which can drastically reduce the computational cost of the corresponding canonical approaches by using a transformation of the wave function amplitudes to a local orbital representation. In the case of the latter, it has been found that the interaction energy of the C60 dimer depends strongly on the cutoff value used to discard pair natural orbitals with small eigenvalues. Using tight cutoff thresholds and performing a CBS extrapolation, we found the counterpoise-corrected binding energy of (C60)2 to be −8.1 kcal mol−1 at a center-to-center distance of 9.75 Å and estimate minimum binding energy at 9.85 Å to be −8.3 kcal mol−1 with an uncertainty of ±0.1 kcal mol−1. We have also performed DFT-SAPT and ACFD(ALDA) calculations for computing the interaction energy profile of (C60)2. Both DFT-SAPT and ACFD(ALDA) possess a lower scaling behavior with respect to the molecular size than the CC and CEPA/1 methods and thus can also be used for the C60

Figure 9. Fitted analytical potentials.

N 2A ⎛ 1 1 2⎞ + − ⎟ ⎜ s4 ⎠ s(s + 1)3 12(2R )6 ⎝ s(s − 1)3 N 2B ⎛ 1 1 2 ⎞ + + − ⎜ ⎟ s(s + 1)9 90(2R )12 ⎝ s(s − 1)9 s10 ⎠

U (r ) = −

(3) r − 2R ,

where s = r is the distance between the centers of the fullerenes, R is the radius of the fullerenes, taken by Girifalco to be 3.55 Å, and A, B, d, L, m, and n are fitting parameters. N 2A ⎛ 1 1 2⎞ + − 4⎟ ⎜ 3 3 6 s ⎠ s(s + 1) 12(2R ) ⎝ s(s − 1) 2 ⎛ CR ⎞⎤ N B −2CRs⎡ 2 e sinh(2CR )⎟⎥ + ⎢⎣(CRs + 1)sinh (CR ) − ⎜⎝ ⎠⎦ 2 CR3s

U (r ) = −

(4)

U (r ) =

A B − 2 (r − d 2)3 (r 2 − d 2)6

(5)

Table 2. Coefficients of Fitted Expressions Girifalco potential (eq 3) A (kcal Å mol−1)

B (kcal mol−1)

6

51

original this work

459.3636 442.373

800632.4 443061.4 Guerin potential (eq 4)

A (kcal Å6 mol−1) 117,120

original this work

421.13 477.66 A (kcal Å12 mol−1)

118

original this work

119

original this work this work this work a

8.47 × 10 13.50 × 1010 10

B (kcal mol−1)

R̅ 2 0.92232 C (Å−1)

R̅ 2

71625.72 13892.91 Smith-Thakkar potential (eq 5)

3.68 3.18

0.9998

B (kcal Å6 mol−1)

d (Å)

R̅ 2

7.258 6.755

0.99731

1.50 × 10 2.07 × 106 Lim potential (Eq 6) 6

m/n

m*n = mn

4 4 (constant) 2.88 (fitted parameter) 3 (constant)

17.743 *4.436 = 314.804 32*8b = 256 28.078*9.745 = 273.62 30*10b = 300 a

a

D 6.93 7.92 8.02 8.06

R̅ 2

L a

10.074 9.832 9.835 9.802

a

0.99833 0.99988 0.99755

Not listed in the original paper. Fitted to the original Smith-Thakkar potential. bConstrained to integers. 281


Article


GAMESS calculations. A.H. gratefully acknowledges financial support of this work through the DFG priority Program No. SPP1807 (“Control of London dispersion interactions in molecular chemistry”). J.T.M. was supported by the Alexander von Humboldt Foundation via a Feodor Lynen Research Fellowship.

dimer without local approximations and with fairly large basis sets that also enable proper extrapolation techniques to reduce the basis set incompleteness error. The DFT-SAPT and ACFD(ALDA) energy profiles shown in Figure 7, which are rather close to each other, are therefore recommended as a reference to derive analytical forms of the potential. We note, however, that the dimerization energy depends quite strongly on the relative orientation of the monomers, so that, in contrast to the original Girifalco potential, Figure 7 takes the specific relative orientation of the monomers into account. The recommended potential has a Born−Oppenheimer dimer dimerization energy of −8.1 ± 0.1 kcal mol−1 at a center-to-center distance of 9.85 Å. The best-fitting analytical equation is the Guerin potential shown in Table 2. Of the DFT-techniques tested, B3LYP-NL, TPSS-D3BJ, B3LYP-D3BJ, or VV10-SCNL with the cc-pVTZ or def2-TZVP basis sets give reasonable agreement with the reference curve. The M06-2X-D3BJ, BP86-D3BJ, and the double-hybrid techniques are less suitable. In the case of the double hybrid density functionals, the reason for the strong overestimation of the reference interaction energy of the C60 dimer is likely the general failure of second-order correlation methods to describe the interaction energy of extended systems.121 This failure is also observable in the large overestimation of the interaction energies yielded by the MP2, SCS-MP2, and SOS-MP2 methods, as found in this work. This shows that, for the description of the interaction energy of large systems, in particular the dispersion energy contribution, one must take intermolecular correlation contributions higher than second order into account. The results of this work indicate that these contributions can be captured well by the recently developed ACFD method, which extends the random-phase approximation both by high order xc-kernel contributions and by symmetry-adapted perturbation theory using a static and timedependent description of the monomer properties (DFTSAPT). The results of this work may serve as a general guide for the investigation of the interaction energy of large complexes using quantum chemistry methods.

■

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00869. Total and relative energies that correspond to Figures 1, 7, 8 and Cartesian coordinates for each geometry (PDF)

■

REFERENCES

(1) Hobza, P. Calculations on Noncovalent Interactions and Databases of Benchmark Interaction Energies. Acc. Chem. Res. 2012, 45, 663−672. (2) Wagner, J. P.; Schreiner, P. R. London Dispersion in Molecular Chemistry-Reconsidering Steric Effects. Angew. Chem., Int. Ed. 2015, 54, 12274−12296. (3) Mahadevi, A. S.; Sastry, G. N. Cooperativity in Noncovalent Interactions. Chem. Rev. (Washington, DC, U. S.) 2016, 116, 2775− 2825. (4) Calbo, J.; Orti, E.; Sancho-Garcia, J. C.; Arago, J. Accurate Treatment of Large Supramolecular Complexes by Double-Hybrid Density Functionals Coupled with Nonlocal Van Der Waals Corrections. J. Chem. Theory Comput. 2015, 11, 932−939. (5) Swart, M.; Sola, M.; Bickelhaupt, F. M. Inter- and Intramolecular Dispersion Interactions. J. Comput. Chem. 2011, 32, 1117−1127. (6) Rezac, J.; Riley, K. E.; Hobza, P. Evaluation of the Performance of Post-Hartree-Fock Methods in Terms of Intermolecular Distance in Noncovalent Complexes. J. Comput. Chem. 2012, 33, 691−694. (7) Grimme, S.; Diedrich, C.; Korth, M. The Importance of Interand Intramolecular Van Der Waals Interactions in Organic Reactions: The Dimerization of Anthracene Revisited. Angew. Chem., Int. Ed. 2006, 45, 625−629. (8) Hesselmann, A.; Jansen, G. The Helium Dimer Potential from a Combined Density Functional Theory and Symmetry-Adapted Perturbation Theory Approach Using an Exact Exchange-Correlation Potential. Phys. Chem. Chem. Phys. 2003, 5, 5010−5014. (9) Hesselmann, A.; Jansen, G.; Schutz, M. Density-Functional Theory-Symmetry-Adapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. J. Chem. Phys. 2005, 122, 014103. (10) Misquitta, A. J.; Szalewicz, K. Intermolecular Forces from Asymptotically Corrected Density Functional Description of Monomers. Chem. Phys. Lett. 2002, 357, 301−306. (11) Podeszwa, R.; Bukowski, R.; Szalewicz, K. Density-Fitting Method in Symmetry-Adapted Perturbation Theory Based on KohnSham Description of Monomers. J. Chem. Theory Comput. 2006, 2, 400−412. (12) Jansen, G. Symmetry-Adapted Perturbation Theory Based on Density Functional Theory for Noncovalent Interactions. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 127−144. (13) Becke, A. D.; Johnson, E. R. Exchange-Hole Dipole Moment and the Dispersion Interaction Revisited. J. Chem. Phys. 2007, 127, 154108. (14) Becke, A. D.; Johnson, E. R. A Unified Density-Functional Treatment of Dynamical, Nondynamical, and Dispersion Correlations. J. Chem. Phys. 2007, 127, 124108. (15) Bremond, E.; Golubev, N.; Steinmann, S. N.; Corminboeuf, C. How Important Is Self-Consistency for the Ddsc Density Dependent Dispersion Correction? J. Chem. Phys. 2014, 140, 18A516. (16) Vydrov, O. A.; Van Voorhis, T. Nonlocal Van Der Waals Density Functional: The Simpler the Better. J. Chem. Phys. 2010, 133, 244103. (17) Hujo, W.; Grimme, S. Comparison of the Performance of Dispersion-Corrected Density Functional Theory for Weak Hydrogen Bonds. Phys. Chem. Chem. Phys. 2011, 13, 13942−13950. (18) 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.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Timothy Clark: 0000-0001-7931-4659 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank the Deutsche Forschungsgemeinschaft (DFG) for financial support within CRC953 “Synthetic Carbon Allotropes” and the Regionales Rechenzentrum Erlangen for computer time. We thank Tobias Klöffel and Bernd Meyer for DFT-MD calculations and Denis Tikhonov for help with 282


Article


(40) Chen, X.; Bai, F.-Q.; Tang, Y.; Zhang, H.-X. How the Substituents in Corannulene and Sumanene Derivatives Alter Their Molecular Assemblings and Charge Transport Properties? - a Theoretical Study with a Dimer Model. J. Comput. Chem. 2016, 37, 813−824. (41) Josa, D.; Rodriguez-Otero, J.; Cabaleiro-Lago, E. M.; dos Santos, L. A.; Ramalho, T. d. C. MDPI AG 2013, e014. (42) Janowski, T.; Pulay, P.; Sasith Karunarathna, A. A.; Sygula, A.; Saebo, S. Convex-Concave Stacking of Curved Conjugated Networks: Benchmark Calculations on the Corannulene Dimer. Chem. Phys. Lett. 2011, 512, 155−160. (43) Josa, D.; Rodriguez-Otero, J.; Cabaleiro-Lago, E. M.; RellanPineiro, M. Analysis of the Performance of Dft-D, M05−2x and M06− 2x Functionals for Studying Π···Π Interactions. Chem. Phys. Lett. 2013, 557, 170−175. (44) Zhao, Y.; Truhlar, D. G. A Prototype for Graphene Material Simulation: Structures and Interaction Potentials of Coronene Dimers. J. Phys. Chem. C 2008, 112, 4061−4067. (45) Sure, R.; Grimme, S. Comprehensive Benchmark of Association (Free) Energies of Realistic Host-Guest Complexes. J. Chem. Theory Comput. 2015, 11, 3785−3801. (46) Hesselmann, A.; Korona, T. Intermolecular Symmetry-Adapted Perturbation Theory Study of Large Organic Complexes. J. Chem. Phys. 2014, 141, 094107. (47) Schatschneider, B.; Monaco, S.; Liang, J.-J.; Tkatchenko, A. High-Throughput Investigation of the Geometry and Electronic Structures of Gas-Phase and Crystalline Polycyclic Aromatic Hydrocarbons. J. Phys. Chem. C 2014, 118, 19964−19974. (48) Marom, N.; Tkatchenko, A.; Kapishnikov, S.; Kronik, L.; Leiserowitz, L. Structure and Formation of Synthetic Hemozoin: Insights from First-Principles Calculations. Cryst. Growth Des. 2011, 11, 3332−3341. (49) Nakamura, M. In Handbook of Nanophysics: Clusters and Fullerenes; Sattler, K. D., Ed.; CRC Press: 2011; Vol. 2, pp 37/1−37/ 15. (50) Anatole von Lilienfeld, O.; Tkatchenko, A. Two- and ThreeBody Interatomic Dispersion Energy Contributions to Binding in Molecules and Solids. J. Chem. Phys. 2010, 132, 234109. (51) Girifalco, L. A. Molecular Properties of Fullerene in the Gas and Solid Phases. J. Phys. Chem. 1992, 96, 858−61. (52) Girifalco, L. A. Interaction Potential for Carbon (C60) Molecules. J. Phys. Chem. 1991, 95, 5370−1. (53) Sweetman, A.; Rashid, M. A.; Jarvis, S. P.; Dunn, J. L.; Rahe, P.; Moriarty, P. Visualizing the Orientational Dependence of an Intermolecular Potential. Nat. Commun. 2016, 7, 10621. (54) Jakubov, T.; Mainwaring, D. Molecular Interactions in Fullerite Nanostructures. In Physics, Chemistry and Applications of Nanostructures - Proceedings of International Conference Nanomeeting - 2011, Borisenko, V. E., Gaponenko, S. V., Gurin, V. S., Kam, C. H., Eds.; pp 102−105. (55) Neese, F. The Orca Program System. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73−78. (56) 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.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−63. (57) Gordon, M. S.; Schmidt, M. W. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier B.V.: Amsterdam, 2005; pp 1167−1189. (58) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: A General-Purpose Quantum Chemistry Program Package. Wiley Interdisciplinary Reviews: Computational Molecular Science 2012, 2, 242−253. (59) Riplinger, C.; Neese, F. An Efficient and near Linear Scaling Pair Natural Orbital Based Local Coupled Cluster Method. J. Chem. Phys. 2013, 138, 034106. (60) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural Triple Excitations in Local Coupled Cluster Calculations with Pair Natural Orbitals. J. Chem. Phys. 2013, 139, 134101.

(19) Johnson, E. R.; Becke, A. D. A Post-Hartree-Fock Model of Intermolecular Interactions: Inclusion of Higher-Order Corrections. J. Chem. Phys. 2006, 124, 174104. (20) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465. (21) 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. (22) 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. (23) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - Eur. J. 2012, 18, 9955−9964 S9955/1-S9955/53. (24) Jones, J. R. The Determination of Molecular Fields (Ii) from the Equation of State of a Gas. Proc. R. Soc. London, Ser. A 1924, 106, 463−77. (25) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proc. R. Soc. London, Ser. A 1938, 168, 264− 83. (26) El Kerdawy, A.; Murray, J. S.; Politzer, P.; Bleiziffer, P.; Heßelmann, A.; Görling, A.; Clark, T. Directional Noncovalent Interactions: Repulsion and Dispersion. J. Chem. Theory Comput. 2013, 9, 2264−2275. (27) Leitherer, S.; Jaeger, C. M.; Halik, M.; Clark, T.; Thoss, M. Modeling Charge Transport in C60-Based Self-Assembled Monolayers for Applications in Field-Effect Transistors. J. Chem. Phys. 2014, 140, 204702. (28) Bauer, T.; Jaeger, C. M.; Jordan, M. J. T.; Clark, T. A MultiAgent Quantum Monte Carlo Model for Charge Transport: Application to Organic Field-Effect Transistors. J. Chem. Phys. 2015, 143, 044114. (29) Jaeger, C. M.; Schmaltz, T.; Novak, M.; Khassanov, A.; Vorobiev, A.; Hennemann, M.; Krause, A.; Dietrich, H.; Zahn, D.; Hirsch, A.; Halik, M.; Clark, T. Improving the Charge Transport in Self-Assembled Monolayer Field-Effect Transistors: From Theory to Devices. J. Am. Chem. Soc. 2013, 135, 4893−4900. (30) Clark, T. Simulating Charge Transport in Flexible Systems. Perspectives in Science 2015, 6, 58−65. (31) Shubina, T. E.; Sharapa, D. I.; Schubert, C.; Zahn, D.; Halik, M.; Keller, P. A.; Pyne, S. G.; Jennepalli, S.; Guldi, D. M.; Clark, T. Fullerene Van Der Waals Oligomers as Electron Traps. J. Am. Chem. Soc. 2014, 136, 10890−10893. (32) Ruzsinszky, A.; Perdew, J. P.; Tao, J.; Csonka, G. I.; Pitarke, J. M. Van Der Waals Coefficients for Nanostructures: Fullerenes Defy Conventional Wisdom. Phys. Rev. Lett. 2012, 109, 233203. (33) Reilly, A. M.; Tkatchenko, A. Van Der Waals Dispersion Interactions in Molecular Materials: Beyond Pairwise Additivity. Chem. Sci. 2015, 6, 3289−3301. (34) Ambrosetti, A.; Alfe, D.; DiStasio, R. A.; Tkatchenko, A. Hard Numbers for Large Molecules: Toward Exact Energetics for Supramolecular Systems. J. Phys. Chem. Lett. 2014, 5, 849−855. (35) Axilrod, B. M.; Teller, E. Interaction of the Van Der Waals Type among Three Atoms. J. Chem. Phys. 1943, 11, 299−300. (36) David, W. I. F.; Ibberson, R. M.; Dennis, T. J. S.; Hare, J. P.; Prassides, K. Structural Phase Transitions in a Solid of Carbon SixtyAtom Fullerene Molecules. Europhys. Lett. 1992, 18, 219−25. (37) Janowski, T.; Ford, A. R.; Pulay, P. Accurate Correlated Calculation of the Intermolecular Potential Surface in the Coronene Dimer. Mol. Phys. 2010, 108, 249−257. (38) Janowski, T.; Pulay, P. A Benchmark Comparison of Σ/Σ and Π/Π Dispersion: The Dimers of Naphthalene and Decalin, and Coronene and Perhydrocoronene. J. Am. Chem. Soc. 2012, 134, 17520−17525. (39) Karunarathna, A. A. S.; Saebo, S. Computational Studies of the Intermolecular Interactions in Dimers of the Bowl-Shaped Sumanene Molecule. Struct. Chem. 2014, 25, 1831−1836. 283


Article

Journal of Chemical Theory and Computation (61) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse Maps - a Systematic Infrastructure for Reduced-Scaling Electronic Structure Methods. Ii. Linear Scaling Domain Based Pair Natural Orbital Coupled Cluster Theory. J. Chem. Phys. 2016, 144, 024109. (62) Schaefer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Lithium to Krypton. J. Chem. Phys. 1992, 97, 2571−7. (63) Schaefer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−35. (64) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−23. (65) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (66) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Ri-Mp2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143−152. (67) Weigend, F.; Kohn, A.; Hattig, C. Efficient Use of the Correlation Consistent Basis Sets in Resolution of the Identity Mp2 Calculations. J. Chem. Phys. 2002, 116, 3175−3183. (68) Hättig, C. Optimization of Auxiliary Basis Sets for Ri-Mp2 and Ri-Cc2 Calculations: Core-Valence and Quintuple-Z Basis Sets for H to Ar and Qzvpp Basis Sets for Li to Kr. Phys. Chem. Chem. Phys. 2005, 7, 59−66. (69) Weigend, F.; Häser, M. Ri-Mp2. First Derivatives and Global Consistency. Theor. Chem. Acc. 1997, 97, 331−340. (70) Grimme, S. Improved Second-Order Moller-Plesset Perturbation Theory by Separate Scaling of Parallel- and Antiparallel-Spin Pair Correlation Energies. J. Chem. Phys. 2003, 118, 9095−9102. (71) Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled Opposite-Spin Second Order M?ller-Plesset Correlation Energy: An Economical Electronic Structure Method. J. Chem. Phys. 2004, 121, 9793−9802. (72) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, Approximate and Parallel Hartree-Fock and Hybrid Dft Calculations. A ’Chain-of-Spheres’ Algorithm for the Hartree-Fock Exchange. Chem. Phys. 2009, 356, 98−109. (73) Izsak, R.; Neese, F. An Overlap Fitted Chain of Spheres Exchange Method. J. Chem. Phys. 2011, 135, 144105. (74) Weigend, F. Accurate Coulomb-Fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (75) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (76) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−100. (77) Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (78) 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. (79) Rezac, J.; Riley, K. E.; Hobza, P. Extensions of the S66 Data Set: More Accurate Interaction Energies and Angular-Displaced Nonequilibrium Geometries. J. Chem. Theory Comput. 2011, 7, 3466−3470. (80) Brauer, B.; Kesharwani, M. K.; Kozuch, S.; Martin, J. M. L. The S66×8 Benchmark for Noncovalent Interactions Revisited: Explicitly Correlated Ab Initio Methods and Density Functional Theory. Phys. Chem. Chem. Phys. 2016, 18, 20905. (81) Heiney, P. A.; Fischer, J. E.; McGhie, A. R.; Romanow, W. J.; Denenstein, A. M.; McCauley, J. P., Jr.; Smith, A. B., III; Cox, D. E. Orientational Ordering Transition in a Solid of Carbon Sixty-Atom Molecules. Phys. Rev. Lett. 1991, 66, 2911−14.

(82) Halkier, A.; Helgaker, T.; Jorgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Basis-Set Convergence in Correlated Calculations on Ne, N2, and H2o. Chem. Phys. Lett. 1998, 286, 243−252. (83) Branz, W.; Malinowski, N.; Enders, A.; Martin, T. P. Structural Transition in (C60)N Clusters. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 094107. (84) Minenkov, Y.; Chermak, E.; Cavallo, L. Accuracy of DlpnoCcsd(T) Method for Noncovalent Bond Dissociation Enthalpies from Coinage Metal Cation Complexes. J. Chem. Theory Comput. 2015, 11, 4664−4676. (85) Liakos, D. G.; Neese, F. Improved Correlation Energy Extrapolation Schemes Based on Local Pair Natural Orbital Methods. J. Phys. Chem. A 2012, 116, 4801−4816. (86) Becke, A. D.; Roussel, M. R. Exchange Holes in Inhomogeneous Systems: A Coordinate-Space Model. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 39, 3761−7. (87) Becke, A. D.; Johnson, E. R. A Simple Effective Potential for Exchange. J. Chem. Phys. 2006, 124, 221101. (88) Langreth, D. C.; Perdew, J. P. The Exchange-Correlation Energy of a Metallic Surface. Solid State Commun. 1975, 17, 1425−1429. (89) Langreth, D. C.; Perdew, J. P. Exchange-Correlation Energy of a Metallic Surface: Wave-Vector Analysis. Phys. Rev. B 1977, 15, 2884− 2901. (90) Pines, D.; Nozieres, P. In Phys. Today; Benjamin: 1966; Vol. 9, pp 99−100. (91) Bleiziffer, P.; Hesselmann, A.; Goerling, A. Efficient SelfConsistent Treatment of Electron Correlation within the Random Phase Approximation. J. Chem. Phys. 2013, 139, 084113. (92) Hesselmann, A.; Goerling, A. On the Short-Range Behavior of Correlated Pair Functions from the Adiabatic-Connection FluctuationDissipation Theorem of Density-Functional Theory. J. Chem. Theory Comput. 2013, 9, 4382−4395. (93) Furche, F.; Van Voorhis, T. Fluctuation-Dissipation Theorem Density-Functional Theory. J. Chem. Phys. 2005, 122, 164106. (94) Olsen, T.; Thygesen, K. S. Extending the Random-Phase Approximation for Electronic Correlation Energies: The Renormalized Adiabatic Local Density Approximation. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 081103. (95) Perdew, J. P.; Schmidt, K. Jacob’s Ladder of Density Functional Approximations for the Exchange-Correlation Energy. AIP Conf. Proc. 2000, 577, 1−20. (96) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133− A1138. (97) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (98) Murray, E. D.; Lee, K.; Langreth, D. C. Investigation of Exchange Energy Density Functional Accuracy for Interacting Molecules. J. Chem. Theory Comput. 2009, 5, 2754−2762. (99) Hujo, W.; Grimme, S. Performance of the Van Der Waals Density Functional Vv10 and (Hybrid)Gga Variants for Thermochemistry and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 3866−3871. (100) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical MetaGeneralized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. (101) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative Assessment of a New Nonempirical Density Functional: Molecules and Hydrogen-Bonded Complexes. J. Chem. Phys. 2003, 119, 12129−12137. (102) Perdew, J. P.; Tao, J.; Staroverov, V. N.; Scuseria, G. E. MetaGeneralized Gradient Approximation: Explanation of a Realistic Nonempirical Density Functional. J. Chem. Phys. 2004, 120, 6898− 6911. (103) Becke, A. D. Density-Functional Thermochemistry. Iii. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−52. 284


Article

Journal of Chemical Theory and Computation (104) Lee, C.; Yang, W.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−9. (105) 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−11. (106) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−7. (107) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-Range-Corrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425−8433. (108) Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid ExchangeCorrelation Functional Using the Coulomb-Attenuating Method (Cam-B3lyp). Chem. Phys. Lett. 2004, 393, 51−57. (109) Chai, J.-D.; Head-Gordon, M. Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. J. Chem. Phys. 2008, 128, 084106. (110) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (111) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (112) Kolmogorov, A. N.; Crespi, V. H. Registry-Dependent Interlayer Potential for Graphitic Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 235415. (113) 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: Condens. Matter Mater. Phys. 2006, 73, 205101. (114) Sato, T.; Nakai, H. Density Functional Method Including Weak Interactions: Dispersion Coefficients Based on the Local Response Approximation. J. Chem. Phys. 2009, 131, 224104. (115) Sato, T.; Nakai, H. Local Response Dispersion Method. Ii. Generalized Multicenter Interactions. J. Chem. Phys. 2010, 133, 194101. (116) Klimes, J.; Michaelides, A. Perspective: Advances and Challenges in Treating Van Der Waals Dispersion Forces in Density Functional Theory. J. Chem. Phys. 2012, 137, 120901. (117) Guerin, H. Analytical Expressions for Atom-Fullerene, Fullerene-Fullerene and Fullerene-Graphite-Surface Interaction Energies Using the Surface Continuum Approximation with an AtomAtom Van Der Waals Buckingham Potential. J. Chim. Phys. Phys.-Chim. Biol. 1998, 95, 561−573. (118) Liu, F.-L.; Wang, J. The Intermolecular Potential Function of Smith-Thakkar Type for C60. J. Mol. Struct.: THEOCHEM 2006, 778, 105−109. (119) Lim, T.-C. Preliminary Assessment of a Multifunctional Potential Energy Function. Mol. Phys. 2010, 108, 1589−1597. (120) Gavezzotti, A. Calculations on Packing Energies, Packing Efficiencies and Rotational Freedom for Molecular Crystals. Nouv. J. Chim. 1982, 6, 443−50. (121) Cybulski, S. M.; Lytle, M. L. The Origin of Deficiency of the Supermolecule Second-Order Møller-Plesset Approach for Evaluating Interaction Energies. J. Chem. Phys. 2007, 127, 141102.

285


Accurate Intermolecular Potential for the C60 Dimer: The Performance

Recommend Documents