Dispersion-Corrected Spin-Component-Scaled ... - ACS Publications

May 24, 2017 - scaled double-hybrid (DSD) density functional theory (DFT) ... of the perturbative treatment of correlation effects is shown to signifi...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/JCTC

Dispersion-Corrected Spin-Component-Scaled Double-Hybrid Density Functional Theory: Implementation and Performance for Non-covalent Interactions Loïc M. Roch† and Kim K. Baldridge*,†,‡,§ †

Department of Chemistry, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland School of Pharmaceutical Science and Technology, University of Tianjin, 92 Weijin Road, Nankai District, Tianjin 3000072, P. R. China



S Supporting Information *

ABSTRACT: The implementation of 300 combinations of generalized gradient approximation/local density approximation exchange−correlation dispersion-corrected spin-componentscaled double-hybrid (DSD) density functional theory (DFT) methods has been carried out and the performance assessed against several DFT and post-Hartree−Fock methods, enabling further advancements toward the long-standing challenge of accurate prediction of interaction energies and associated properties. The resulting framework is flexible and has been further extended to include the resolution of identity (RI) approximation for solving the critical four-center two-electron repulsion integrals in the basis of the Kohn−Sham orbitals for cost effectiveness. To evaluate the performance of this set of new cost-effective methods, denoted as RI-DSD-DFTs, seven validation data sets were designed to cover a broad range of non-covalent interactions with characteristic stabilizing contributions. Inclusion of the perturbative treatment of correlation effects is shown to significantly improve the description of weak interactions. The set of DSD-DFTs provide interaction energies with root-mean-square deviations and mean absolute errors within 0.5 kcal/mol. The cost-effective RI-DSD-DFT counterparts deviate by less than 0.18 kcal/mol on average with only 2% of the computational cost.



INTRODUCTION Considerable efforts in the past decade have been extended toward achieving chemical accuracy for predictions of interaction energies and associated phenomena in systems where weak interactions dominate.1−5 Despite this, there are still important gaps in our present array of first-principles methodologies that enable predictions within 0.5 kcal/mol of experimental values. The issue can be further complicated by the type as well as size of system being investigated. For example, investigations into largescale polynuclear aromatic hydrocarbons for materials design pose challenges for effective treatment of electron correlation, dispersion, and polarization in a manner appropriate to the size of the system and materials property being investigated.6,7 Delocalization lies at the heart of the aromatic character in these systems, where large surface areas can amplify the importance of van der Waals (vdW) interactions and polarizable electron distributions. While dispersion interactions involve pure electron correlation effects,8,9 delocalization comprises both short and long-range effects, all of which require appropriate methodology to capture. Among the variety of quantum-chemical methodologies, challenges of system size are easily reached, making Kohn−Sham density functional theory (DFT) as the theoretical method of choice. However, the appropriate choice of functional becomes © XXXX American Chemical Society

another challenge, complicated by the effectiveness of the treatment of dispersion, which parallels the sophistication of the exchange−correlation potential. In the particular class of polynuclear aromatic hydrocarbon systems, adequate treatment of delocalization requires a proper balance of exchange (accounting for antisymmetry due to Pauli exclusion) and correlation (accounting for many-body effects). Despite the success of standard functionals over the last 20 years, it is widely known that most do not provide an adequate description of nonlocal dispersion forces.1,10 Recent attempts to overcome insufficiencies in the long-range electron correlation in generalized gradient approximation (GGA) functionals have led to several important strategies, including empirical dispersioncorrected functionals,11−13,12−15 specialized functionals,16−20 and the dispersion-corrected spin-component-scaled doublehybrid (DSD) method.21,22 Within that last category of methods is a promising double-hybrid scheme introduced by Martin with demonstrated high accuracy for correlated chemical systems of small to medium size.21−23 In that scheme, the incorrect longrange asymptotic behavior of standard functionals is treated by mixing the approximate DFT exchange−correlation functional Received: March 2, 2017

A

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

reference values is reported in Tables S1−S7 in the Supporting Information. Single-point energetic analysis for structures in these sets was carried out using the Def2-QZVPD47 basis set with representative DFT48,49 types, including 10 GGA DFs (BLYP,50,51 BPBE,50,52,53 BP86,50,54 BPW91,50,55 PBELYP,51−53 PBE,52,53 PBEP86,52−54 PBEPW91,52,53,55 OLYP,51,56 and B97-D11), one meta-GGA DF (tHCTH57), two hybrid DFs (B3LYP58 and PBE059), one meta-hybrid-GGA DF (tHCTHhyb57), two range-separated DFs (CAMB3LYP60 and ωB9761), and 10 double-hybrid DFs (B2PLYP,17 DSD-BLYP,21 DSD-BPBE,22 DSD-BP86,22 DSD-BPW91,22 DSD-PBELYP,22 DSD-PBEPBE,22 DSD-PBEP86,22 DSD-PBEPW91,22 and DSD-OLYP22), MP2,34 SCS-MP2,35 and CCSD(T).36 All of the DSD-DFT and RI-DSD-DFT implementations were carried out with a local version of GAMESS, which will be made available in the community version. The new logical keyword XCDH ensures the branching between standard GGA and DSD-DFTs. The cost-effective RI-DSD-DFTs are selected by setting the new keyword RIDH to TRUE in the $DFT group. This enables the RI approximation for solving the critical four-center two-electron repulsion integrals in the basis of the Kohn−Sham orbitals. The choice of basis set was assessed from an exhaustive study comparing the interaction energies obtained using the HF, MP2, and DSD-PBEPBE wave function types with their corresponding HF/complete basis set (CBS), MP2/CBS, and DSD-PBEPBE/ CBS limits. Among the series of double-, triple-, quadruple-, and quintuple-ζ basis sets investigated, Def2-QZVPD provides the best accuracy, on par with quintuple-ζ basis sets, and additionally outperforms the corresponding Dunning-style quadruple-ζ set at a lower computational cost (further details on the accuracy can be found in the Supporting Information). At the DFT level, a Lebedev62 grid with 155 radial points in the Euler−MacLaurin quadrature and with 1202 angular points (NRAD = 155, NLEB = 1202) was used to solve the integrals. Both the D211 and D313 empirical corrections of Grimme were compared as the standard DFT energy to correct for the missing long-range attraction. It is worth mentioning that the damping function from Grimme was used when the DFT-D2 correction scheme was applied. At the MP2 level, second-order perturbation energies were obtained on the unoccupied orbital space only (i.e., omitting the chemical core orbitals) and were shown to differ by 0.02 to 0.3 kcal/mol from the calculations on the full orbital space (i.e., including the chemical core orbitals). For the CCSD(T) calculations, the direct SCF density convergence criterion was lowered to 2.5 × 10−7 (CONV = 2.5D-07), and only integrals with small exponents were dropped out for the SCF (ICUT = 11). The default values for the twoelectron integral transformation were lowered to CUTTRF = 1D-11. In all cases, because the systems contain a large set of diffuse functions, all of the two-electron contributions to the Fock matrices were computed at each SCF step (FDIFF = FALSE). The interaction energy Eint for a complex was computed according to eq 1:

with a pertubative correction term within the basis of the Kohn−Sham orbitals. One drawback to this strategy is extension to large-scale molecular systems because of the computationally intensive cost of the perturbative correction component. Consideration of the performance of any methodology requires adequate benchmarking for validation and assessment. This creates the need for appropriate validation data sets consisting of highly reliable data appropriate for the phenomenon being investigated. While standardized validation sets are valuable for uniform treatment and assessment of methodology, it is also important that they enable extensibility and transferability of the resulting model beyond the validation set. To address each of the aforementioned issues, efforts in the present work focused on (i) performance assessment of single-reference methods for accurate description of electron correlation effects across a broad class of system types, (ii) assembly of reliable validation sets for comparison of new and existing methods, and (iii) development and implementation of a strategy enabling cost-effective treatment of large structures without compromise in accuracy. Toward this end, a generalized implementation of the double-hybrid scheme was used to enable the availability of ca. 300 combinations of GGA/LDA exchange− correlation functionals, which can be run either as stand-alone GGAs or as DSDs. For the sake of cost-effectiveness, the implementation was extended to enable the resolution of identity (RI) approximation24−29 to facilitate solving the critical fourcenter two-electron repulsion integrals in the basis of the Kohn−Sham orbitals, generating the corresponding set of RI-DSD-DFTs, which are particularly well suited to address large-sized systems. On the basis of the previous work initiated by Zhao and Truhlar30,31 and Mikami and co-workers32 in the search for small benchmark data sets, six new validation data sets were introduced for benchmarking and comparison across the methods as well as comparison with conventional wave function theory (WFT), including Hartree−Fock,33 Møller−Plesset second-order perturbation theory (MP2),34 spin-componentscaled MP2 (SCS-MP2),35 and coupled-cluster CCSD(T).36 Application of the methods to several key studies analyzing molecular recognition properties for curved polynuclear aromatic hydrocarbon materials demonstrates the capabilities with a high degree of accuracy compared to experiment. Systems in these investigations typically reach sizes that make them prohibitively expensive using DSD-DFT or conventional wave function strategies (>4000 Cartesian Gaussian basis functions) and therefore provide a good “test bed” for the RI-DSD-DFT methods. The results demonstrate predictive strengths of the methods as well as suggested avenues of further development. All of the implementations have been carried out in the computational chemistry suite GAMESS,37 an open-source quantum chemistry software.



COMPUTATIONAL SETUP All of the calculations reported here were carried out using a locally modified version of the GAMESS37 (2014R1) electronic structure program. Validation studies were carried out using six extended existing and newly designed data sets covering hydrogen-bonding interactions [HB9 (extended HB6/0424)], internal dispersion (IDISP4), alkane dimerization interactions (ADIM4), rare-gas dimer/trimer interactions (RG21), dipolar interactions [DI9 (extended DI6/0424)], and π−π interactions (PPS11) together with the charge-transfer interactions data set of Zhao and Truhlar (CT7/04).30,31 Geometries in these sets were obtained from existing sets;30−32,38−46 a complete tabulation of

E int = E AB − E A − EB

(1)

where EAB is the total energy of the complex AB and EA and EB are the total energies of the isolated monomers A and B, respectively. B

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. Seven validation sets used to assess the performance of the computational methods. The data sets cover (A) hydrogen bonds (HB9), (B) charge-transfer interactions (CT7/04, (C) dipole interactions (DI9), (D) alkane dimerization processes (ADIM5), (E) internal dispersion forces (IDISP4), (F) π−π interactions (PPS11), and (G) rare-gas dimers and trimers (RG21).



VALIDATION SETS Despite the large number of existing validation sets, for this study it was necessary to explore more specific sets related to applications in the area of polynuclear aromatic hydrocarbon nanostructures. In particular, interest lies in the area of curved aromatic building blocks, where the significant curvature confines the electronic states anisotropically in one or more of the

nanoscale directions. Derivatives based on the smallest curved aromatic bowl structure, C20H10 (1), have large surface areas and significant dipole and polarizability properties, manifesting significant variations in the associated material properties. Accurate prediction of properties in structural arrays resulting from supramolecular ordering of such building blocks requires consideration of dispersive effects, electronic correlation, and C

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation quantum confinement and thereby naturally falls to firstprinciples approaches. With this in mind, validation sets were designed with a philosophy not unlike those developed by Zhao and Truhlar,30 which explore the performance of a key methodology as a function of the degree and type of non-covalent interactions, in this case across a broad array of curved aromatic nanostructures. The data sets were categorized into classes in accord with key types of non-covalent contributions, with associated interaction energies that range over ∼20 kcal/mol. The charge-transfer data set of Zhao and Truhlar, CT7/04,30,31 was additionally considered an important data set to include with the six designed/extended data sets. The resulting seven validation sets (Figure 1) include hydrogen-bonding interactions, chargetransfer interactions, dipole−dipole interactions, vdW interactions via π effects, interaction between rare-gas atoms, and interaction between alkane dimers. Overall, 66 chemical systems with interaction energies ranging from 18.6 to 0.02 kcal/mol were involved in the validation studies. The benchmark data collected from the literature are of MP2/CBS (counterpoise (CP) where available) and CCSD(T)/CBS (CP where available) quality. It has to be mentioned that both the IDISP4 and RG21 sets contain experimental interaction energies. The complete tabulation of reference values is reported in Tables S1−S7, and the XYZ coordinates of the 66 complexes are listed in section VI in the Supporting Information.

convergence rates. The reported ab initio energies have uncertainties of up to 0.1 kcal/mol, DSD-DFs have uncertainties of up to 0.04 kcal/mol, and DFTs have uncertainties of up to 0.01 kcal/mol (see the Supporting Information). Figure 2 provides a perspective of the overall performance across all of the interaction types represented in the validation

Figure 2. RMSD results for the 38 methodologies compared for the seven validation sets depicted in Figure 1. Light blue: standard DFTs (i.e., no dispersion correction). Red: DFT+D2. Yellow: DFT+D3. Green: DSD-DFTs. Solid black line: CCSD(T). Dashed black line: SCS-MP2. Dashed/dotted black line: HF. For clarity, the MP2 result (RMSD = 0.673 kcal/mol) is not displayed, as it is essentially the same as the SCS-MP2 result.



RESULTS AND DISCUSSION From the plethora of dispersion-corrected density functionals,12−15 this work focuses on dispersion-corrected spincomponent-scaled double-hybrid (DSD) methods.21−23 In the original formulation of Martin,21 the exchange−correlation term (EXC) for the DSD-DFT type is represented as

sets (66 complexes) in terms of RMSD (in kcal/mol) across all wave function types with the Def2-QZVPD basis set. Of the conventional wave function techniques, CCSD(T) shows the overall best performance, as expected, with an RMSD of 0.246 kcal/mol, significantly exceeding our targeted benchmark value of 0.5 kcal/mol. On the low-performance end, as also might be anticipated, is HF theory, with an RMSD of 2.598 kcal/mol. However, although HF must set the upper RMSD boundary, it appears to outperform three cases of standard GGA DFs without dispersion corrections, in particular the DFs with the PW91 correlation functionalBPW91 and PBEPW91as well as OLYP. PBE without empirical dispersion corrections in several cases outperforms its respective D3-dispersion-enabled counterparts. Among the standard GGAs, PBELYP provides the best results with an RMSD of 1.494 kcal/mol. On average, non-GGA functionals show an across-the-board improvement over standard GGAs. The ωB97 functional with no empirical dispersion shows the best performance with an RMSD of 0.4557 kcal/mol. Also, the popular B3LYP hybrid functional with the inclusion of either the D2 or D3 dispersion correction shows good performance, competing with SCS-MP2. With the exception of the few cases already mentioned, inclusion of an empirical correction generally improves the performance of the different DFs (D2 correction in red and D3 correction in yellow in Figure 2). Of the cases where the empirical correction fails to provide improved results, the D3 correction tends to overestimate the nonbonding interactions, although with decreased variance. Figure 3 shows the trends in the computed interaction energies against their reference values (either MP2/CBS or CCSD(T)/CBS as detailed in Tables S8−S12). Panel A displays standard GGAs, panel B the GGAs with the

E XC = C XE XGGA + (1 − C X)E XHF + CCECGGA + COEOMP2 + CSESMP2 + S6E D

(2)

where CX (CC) is the amount of DFT exchange (correlation) at the GGA level, CO (CS) the opposite (same)-spin MP2 contribution, and S6 is the dispersion energy (ED) contribution. The latter term is obtained via the D2 scheme, where both the exponential parameter to scale the van der Waals radii and the linear parameter for scaling the C6 term are taken from the corresponding GGA functional and then scaled down according to S6 in eq 2. The computationally intensive component of this strategy is attributed to the perturbative MP2 component. The goal in this work is to enable extension of such methods efficiently and economically to large-scale aromatic systems with a chemical accuracy within our target of 0.5 kcal/mol. Initial performance analysis was carried out on 38 wave function types spanning a variety of standard density functional (DF) types, double-hybrid DFs, and conventional wave function theory (WFT), as detailed in Computational Setup. The computations were carried out with the Def2-QZVPD basis set.47 Whenever possible, the DF types were analyzed without and with dispersion corrections, the latter comparing the D211 and D313 empirical dispersion corrections. The detailed data, including root-mean-square deviations (RMSDs) and mean absolute errors (MAEs) for the 38 methodologies, are summarized in Tables S8−S12. I. Overall Performance of Standard DFTs, DSD-DFs, and WFTs. Prior to the discussion of the performance of the methodologies, it is prudent to discuss the difference in D

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 3. Climbing the stairway to heaven of DFs. Shown are correlation graphs of reference values vs computed values for (A) GGAs, (B) GGA+D2, (C) GGA+D3, and (D) DSD-DFs. Stars refer to the B88 exchange functional and squares to PBE; blue refers to the LYP correlation functional, red to PBE, cyan to P86, and green to PW91.

D2 correction, panel C the GGAs with the D3 correction, and panel D the double-hybrid DSD-DFTs. Climbing the stairway to heaven1 in accord with DF type, in going from panel A to panel B one finds a general improvement in the description of the weak component of the interaction energies. Panels B and C appear to be rather similar, but a closer look shows the trend for D3 correction is toward overbinding, as evidenced by the slight left shift of the data points along the horizontal axis. Inclusion of a perturbative treatment of the correlation energy with the DSD-DFs results in a considerable improvement in the performance of the GGAs (green bars in Figure 2 and panel D in Figure 3). In terms of RMSD, all of the double hybrids (DHs) are bracketed between the predicted SCS-MP2 value (dashed line in Figure 2) and the state-of-the-art CCSD(T) value. MP2 and SCSMP2 have very similar RMSD values (0.673 and 0.665 kcal/mol, respectively). Of the DHs, DSD-OLYP shows the highest RMSD and is on par with the MP2 and SCS-MP2 results (RMSD = 0.686, 0.673, and 0.665 kcal/mol, respectively). On the other hand, the lowest RMSD is found for DSD-PBEPBE, which even challenges the performance of coupled cluster (RMSD = 0.299 kcal/mol vs 0.246 kcal/mol). The overall substantial improvement in performance with inclusion of the perturbative treatment is highlighted in panel D of Figure 3, where an essentially perfect alignment of the 594 interaction energies with that of the correlation line (depicted in black) is observed. Out of all the methods tested, the

set of methods that hit our target benchmark of 0.5 kcal/mol includes all of the DSD DFs, B3LYP with dispersion correction, and ωB97 without any correction. II. Performance of DSD-DFs for Specific Types of Stabilizing Interactions. To better understand the more global behavior profile discussed in the last section, it is of interest to consider the performance with respect to each of the seven individual validation data sets presented above. Figure 4 shows a summary of results for each of the seven sets, and the individual sets of raw results are given in the Supporting Information. HB9. Hydrogen bonds represent an important type of interaction in the construction of functional architectures and supramolecular materials. For example, bioconjugated corannulene cores can exploit hydrogen-bonding networks to form functional supramolecular structures in solution, as in dimeric or polymeric columnar assemblies of the corannulene core stabilized by hydrogen bonding between oligopeptide chains.63 The hydrogen-bonding structure is typically complex because of the low mass and therefore quantum nature of the proton. Moreover, with increasing hydrogen-bond strength the proton typically becomes more delocalized, commensurate with a decrease in the frequency of the bond with the atom to which it is attached. The performance results for this set are shown in Figure 4A. The poorest performance is obtained with BPW91, PBEPW91, and OLYP. Without any dispersion correction, PBELYP, PBEPBE, E

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 4. RMSDs of the 38 methodologies for validation sets (A) HB9, (B) CT7/04, (C) DI9, (D) ADIM5, (E) IDISP4, (F) PPS11, and (G) RG21 and (H) for all sets combined. Values in white are for the standard DFTs (i.e., no dispersion correction), those in red for DFT+D2, those in yellow for DFT +D3, and those in green for the DSD-DFTs. The solid black line is for CCSD(T), the dashed black line for SCS-MP2, and the dashed/dotted black line for HF. For clarity, MP2 is not displayed. Values from Tables S8−S12 represent data displayed in panels A−H.

provide some of the best performance, with the interesting exception of BLYP+D2(D3). CCSD(T) provides the overall best performance with an RMSD of 0.091 kcal/mol, whereas SCS-MP2 exceeds our target threshold with an RMSD of 0.802 kcal/mol.

PBE0, CAMB3LYP, and tHCTH-hyb all perform within our target range. Interestingly, the empirical dispersion correction leads to more accurate prediction for only about half of the methods, with the other half being overestimated significantly compared with the uncorrected DF. The DSD-DFTs again F

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation CT7/04. The recent strong interest in curved aromatic materials in general and organic semiconductors in particular is accentuated by the rich and diverse array of new concepts being highlighted relating their curvature and π character to properties especially significant for use in device technology.6,7 In this category of functional components are structures that show charge-transfer behavior, particularly in the crystal environment. Thus, the ability of methods to capture charge-transfer character effectively is quite important. The validation set CT7/04 represents a set of charge-transfer complexes assembled by Zhao and Truhlar.30,31 The performance results for this validation set across the series of methods are depicted in Figure 4B. It is perhaps surprising that the majority of the DFs have RMSD values considerably above the target accuracy of 0.5 kcal/mol. Even SCS-MP2 has an RMSD of 0.571 kcal/mol, which just exceeds the target threshold. Even inclusion of the empirical dispersion correction does not significantly help the situation. In fact, inclusion of empirical dispersion in many cases results in an increase in RMSD. However, addition of the perturbative correction for correlation effects in the DSD-DFTs has a large impact. The exceptions in this category are again DSD-PBELYP, DSD-PBEPBE, and DSD-OLYP, for which the RMSD values exceed the target value of 0.5 kcal/mol. In addition, CAMB3LYP and ωB97 without any dispersion correction are also within the target threshold. The much more expensive CCSD(T) treatment also enables accuracy within the target threshold, with an RMSD of 0.349 kcal/mol. DI9. Another strong feature of curved aromatic bowl structures based on corannulene is the ability to fine-tune the molecular structure and therefore the electronic and dynamic properties through functionalization. A consequence of the local curvature of such structures is their rather substantial intrinsic dipole moment. In the series of structures with increasing bowl depth (C20H10, C30H10, C40H10, and C50H10), the intrinsic dipole moment ranges from 2.1 to 6.0 D.7 Both intrinsic and local dipole features of these types of curved aromatic structures manifest important static and field-oriented properties that are important in materials design, such as crystal packing and transport properties.6,64,65 The performance results for the set of dipole−dipole interaction systems are depicted in Figure 4C. In this case, of the standard DFs, only PBELYP and ωB97 are within the target threshold of 0.5 kcal/mol. Inclusion of an empirical dispersion correction brings B3LYP+D2 and B3LYP+D3 to within the target threshold. All of the remaining DFs are well over the target. One sees a general improvement in the non-GGA DFs over the GGA DFs, but the results of all but the above-mentioned ones

exceed the 0.5 kcal/mol target. With the exception (again) of OLYP and PBELYP, the DSD-DFTs provide good accuracy within the target threshold, with DSD-PBEPW91 displaying the lowest RMSD (0.320 kcal/mol). SCS-MP2 also provides results within the threshold at RMSD = 0.375 kcal/mol, and CCSD(T) shows the overall best performance with an RMSD of 0.192 kcal/mol. ADIM5. An important property of the curved aromatic class of structures is their bowl-to-bowl interconversion ability, as controlled by the variations in structure through functionalization of the core. Corannulene itself has a rather shallow bowl structure (see Figure 9), resulting in a rapid bowl-to-bowl interconversion rate. Predictable trends in the bowl-inversion barrier with bowl depth and curvature have been determined among corannulene and associated derivatives, providing a way to gauge steric or electronic aspects of molecular recognition in supramolecular processes.66 It has been shown that functionalized derivatives with transannular interactions of substituents play a role in the dynamic process.67 The ability to understand the structure−dynamic relationships in these derivatives requires computational methods that can pick up these effects, in this case the alkane interactions. Figure 4D depicts the RMSD results for this alkane interaction validation set. In this case, only one of the DFs, the non-GGA DF ωB97, provides results within the target accuracy of 0.5 kcal/mol. Across the set, the standard GGAs fail to describe the attractive interactions between alkane monomers. This observation is additionally illustrated in Figure 5, which compares correlation graphs for the ADIM5 test set at the (A) GGA and (B) non-GGA levels (Figure 5A is an enlargement of Figure 3A but restricted to ADIM5). Here one can also see that the non-GGA functionals other than ωB97 also provide very poor results (e.g., Figure 5B). The ωB97 DF actually outperforms MP2 and SCS-MP2 in this case, with an RMSD of 0.208 kcal/mol. For the other DFs, inclusion of empirical dispersion with either D2 or D3 recovers the missing dispersion and lowers the RMSD; however, all but BP86+D2, B3LYP+D3, and PBE0+D2(D3) are still considerably over our target threshold. Inclusion of the perturbative term via the DSD-DFT ansatz results in RMSDs well within the target of 0.5 kcal/mol and very close to the accuracy of CCSD(T). The clear exceptions are OLYP and PBELYP, as observed for the previous validation sets discussed. It should be mentioned that the MP2 results are within 0.5 kcal/mol (RMSD = 0.242 kcal/mol) and that the SCS-MP2 results are outside the target threshold with an RMSD of 0.814 kcal/mol.

Figure 5. Correlation graphs for the ADIM5 test set at the (A) GGA and (B) non-GGA levels. (A) is an enlargement of Figure 3A restricted to ADIM5. G

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation IDISP4. The facile functionalization of corannulene enables rapid access to a myriad of derivatives, providing a path toward creative tuning of the electronic properties for specific materials, for example, ones suitable for organic light-emitting diodes or organic field-effect transistors.68−75 Many derivatives have internal dispersive interactions among the substituents that control the dynamic behavior,66 in particular when the magnitude and importance of such interactions are modulated from one bowl conformation through a flat transition state structure to another bowl structure.66,76 High-symmetry functionalized derivatives can have more than one isomeric form, the variation of which can come down to differences in internal dispersive interactions. Therefore, suitable computational methodology must be able to pick up such interactions in a predictable way to distinguish isomeric forms. The benchmarking results for the IDISP4 set are shown in Figure 4E. In this case, unlike the previous results involving the other validation sets, the RMSD values across all methods are under 1.6 kcal/mol. However, also unlike the previous results, none of the DF methods provide RMSDs within the target threshold of 0.5 kcal/mol, including the DSD-DFs. Unfortunately, inclusion of either the empirical dispersion correction or correlation via MP2 perturbation (via D2, D3, or DH) does not lead to noticeable improvements. In this case, HF, BPW91, and OLYP have the largest RMSD values at 1.462, 1.373, and 1.576 kcal/mol, respectively. The best agreement is obtained with ωB97, with an RMSD of 0.703 kcal/mol. Unfortunately, even the post-HF methods, i.e., MP2, SCS-MP2, and CCSD(T), exceed the target threshold, with RMSDs of 0.858, 0.829, and 0.744 kcal/mol, respectively. Notably, MP2 and SCS-MP2 provide very similar results, suggesting that this set of systems is governed by non-covalent interactions that take place in the strong regime. These rather confusing results could be related to the fact that the reference values are all from experimental results and therefore may have additional factors that have not been included in these gas-phase benchmark studies. PPS11. Curved aromatic structures can serve as mimics of fullerene structures or open cavities for host/guest chemistry. In the former case, corannulene has been used as a base component toward the creation of several types of corannulene cyclophane motifs.77 In the designed “basketball” cyclophane motif, for example, depending on the orientation of the opposing corannulene fragment, one finds close contact (∼3.3 Å) between the aromatic isomeric structures, mimicking the arrangement of corannulene fragments in C60. Unlike C60, however, these cyclophane structures have an open cavity, the size of which can be tailored for particular functions in host−guest chemistry. Several other examples in the literature show the effects of curvature on π−π stacking interactions for prototype dimeric systems.78 Benchmarking results for the set of 11 π−π systems, PPS11, are shown in Figure 4F. Of the group of DFs, only ωB97 reaches the target threshold with an RMSD of 0.311 kcal/mol. However, it is important to note that the interaction energy is of the wrong sign for all of the GGA DFs. This is illustrated in Figure 6 (an enlargement of Figure 3A restricted to PPS11), which shows that all 88 points predict repulsion between the monomers. The non-GGA functionals recover incorrect long-range asymptotic behavior and predict a favorable interaction between the monomers in each case (e.g., Figure 3B). Inclusion of an empirical dispersion correction improves the RMSD in all cases. However, only in the cases B3LYP+D2(D3) and PBE0+D2(D3)

Figure 6. Correlation graph for the PPS11 test set. Stars refer to the B88 exchange functional and squares to PBE; blue refers to the LYP correlation functional, red to PBE, cyan to P86, and green to PW91. This figure is an enlargement of Figure 3A restricted to PPS11.

is the correction enough to bring the RMSD within the target range of 1 kcal/mol for all of the DSD-DF methods. These data points correspond to the small cluster of data points in the 0−4 kcal/mol range that stray off the correlation line in Figure 8 (see the Supporting Information for more a detailed comment). Application of the RI approximation shows a specific limitation for these two structures in particular, resulting in an error of ∼50% in the prediction of the interaction energy. Table 2 compares the computational times for the DSD-DF and RI-DSD-DF sets of methods for the seven validation sets. The CPU time in each category represents the sum of CPU time for each member of the validation set for each of the nine DSD-DFT types. The ratio is then determined as the total CPU time for all of the DSD-DF methods divided by the total CPU time for all of the RI-DSD methods for that particular validation set. The overall CPU time is taken as the sum across all validation sets of both the DSD-DFs and RI-DSD-DFs. The total CPU time required to compute all 594 interaction energies was 155 days across all of the DSD-DFT methods. Implementation of the RI approximation brought this down to only 3 days, giving an average speedup factor of 51.5. Given the considerable savings in computational time together with the overall excellent agreement to the standard DSD-DFT method (with the exception of the two outlier complexes), one can conclude that the RI-DSD-DFT method is particularly well suited as an accurate cost-effective substitute for the treatment of systems governed by correlation

fail to describe non-covalent interactions. Non-GGA functionals, however, show an across-the-board improvement, with ωB97 giving an overall accuracy competing with that of the DSD-DFTs. The ωB97 hybrid functional includes 100% long-range exact exchange partitioned into short- and longrange interactions through the error function. Additionally, both the exchange and correlation contributions to the total energy are decomposed into same-spin and opposite-spin components. The high accuracy of the range separation is further observed by comparing B3LYP to its Coulomb-attenuated version, CAMB3LYP. Indeed, the RMSD of CAMB3LYP (1.150 kcal/mol) is on average 0.542 kcal/mol lower than the RMSD of B3LYP (1.692 kcal/mol). Addition of a dispersion correction scheme on top of the DFT energy has generally been seen to improve the description of weak interactions, as depicted in Figure 3 and already reported in the literature.11,13,84−88 The main difference between DFT+D2 and DFT+D3 lies in the dispersion scheme yielding a reasonable −1/r6 asymptotic behavior for the interactions of particles in the gas phase.1 In the DFT+D2 approach, the dispersion coefficients are calculated through the combination of ionization potentials and static polarizabilities of isolated atoms. The DFT+D3 approach captures the chemical environmental dependence of the dispersion coefficients by considering the number of neighbors each atom has. In the current analysis, addition of D2 and D3 corrections showed mixed performance with regard to the validation studies. To capture this, values of ΔD3 and ΔD2 have been calculated as the deviations of DFT+D3 and DFT+D2 values from the reference values for the respective complexes in the seven validation sets, as determined in the previous section. Figure 7 displays ΔD3 as a function of ΔD2. For the 66

Figure 7. DFT+D3 deviation (ΔD3) as a function of the DFT+D2 deviation (ΔD2) for the 66 complexes of the seven validation sets. Green: improvement of DFT+D3 over DFT+D2. Red: overcorrection of DFT+D3 over DFT+D2.

complexes making up the seven validation sets, this graph shows improvement of DFT+D3 over DFT+D2 (green) and overcorrection of DFT+D3 over DFT+D2 (red). Notably, among the interactions overcorrected by DFT+D2, ca. 50% are further overcorrected by the use of DFT+D3 as highlighted with a red background in Figure 7. For the set of DSD-DFTs, the dependence of the accuracy on the functional is less pronounced than in the respective GGAs. From the worst-case performance of GGA, BPW91 (RMSD = 4.488 kcal/mol), to the best performance, PBELYP (RMSD = 1.494 kcal/mol), I

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 1. Root-Mean-Square Deviations and Mean Absolute Errors (in Italics, in Parentheses) of RI-DSD-DFTs method

HB9

CT7/04

DI9

ADIM5

IDISP4

PPS11

RG21

overall

RI-DSD-BLYP RI-DSD-BPBE RI-DSD-BP86 RI-DSD-BPW91 RI-DSD-PBELYP RI-DSD-PBEPBE RI-DSD-PBEP86 RI-DSD-PBEPW91 RI-DSD-OLYP

0.103 (0.057) 0.076 (0.038) 0.084 (0.043) 0.077 (0.039) 0.106 (0.060) 0.078 (0.041) 0.086 (0.046) 0.078 (0.041) 0.099 (0.056)

0.079 (0.044) 0.064 (0.035) 0.069 (0.038) 0.064 (0.035) 0.083 (0.046) 0.067 (0.036) 0.072 (0.040) 0.066 (0.036) 0.078 (0.043)

0.532 (0.203) 0.468 (0.173) 0.491 (0.185) 0.466 (0.173) 0.557 (0.213) 0.484 (0.180) 0.507 (0.191) 0.474 (0.177) 0.520 (0.198)

0.498 (0.377) 0.325 (0.232) 0.376 (0.278) 0.332 (0.238) 0.488 (0.378) 0.319 (0.235) 0.370 (0.279) 0.320 (0.237) 0.405 (0.307)

0.217 (0.181) 0.167 (0.131) 0.166 (0.141) 0.165 (0.127) 0.205 (0.176) 0.171 (0.143) 0.165 (0.137) 0.164 (0.136) 0.197 (0.172)

1.120 (0.516) 1.012 (0.471) 1.036 (0.471) 1.007 (0.467) 1.178 (0.537) 1.048 (0.481) 1.073 (0.485) 1.025 (0.470) 1.119 (0.521)

0.001 (0.001) 0.001 (0.001) 0.002 (0.001) 0.203 (0.145) 0.280 (0.205) 0.001 (0.001) 0.001 (0.001) 0.208 (0.152) 0.264 (0.192)

0.533 (0.174) 0.471 (0.143) 0.486 (0.150) 0.481 (0.184) 0.577 (0.238) 0.486 (0.148) 0.502 (0.154) 0.489 (0.188) 0.543 (0.222)

reference is made to the highest available theoretical result in accord with the benchmarking results in this paper, and in several cases serves as prediction for future studies. Application I: Dynamics and Complexations Involving Corannulene. Although corannulene (1) (Figure 9) is a structural

Figure 8. Correlation of interaction energies between the set of DSD-DFs and the respective RI-DSD-DFs for the seven validation data sets (HB9, CT7/04, DI9, ADIM5, IDISP4, PPS11, and RG21).

Figure 9. Structures and associated bowl depths of corannulene (1), pentaindenocorannulene (2), and C60 (3).

effects. In particular, RI-DSD-DFT methods hold particular promise for large-scale molecular complexes. V. Illustrative Applications. The extended family of curved polynuclear aromatic hydrocarbon fragments based on the parent corannulene91 provides a particularly suitable set of systems to illustrate performance capabilities across methods discussed in this article. The large structural diversification of this motif and the availability of comparable experimental data offer a means of testing a variety of structural, electronic, and dynamic properties of single molecules as well as molecular recognition complexes. Elucidation of accurate methodology for applications involving curved aromatics is fundamental to the design of new classes of functional building blocks as well as better understanding of the mechanistic aspects of device technologies.6,7 Each of the examples below enables analysis of the predictability of the various classes of functionals with regard to dynamics and molecular recognition properties of curved aromatics, which involve various types of weak interactions. Results are compared across the functional classes evaluated in the previous sections and in comparison to experiment where known. In cases where no experimental data are yet available,

subunit of C60 (3) and the smallest bowl-shaped polynuclear aromatic hydrocarbon, the surface curvature of the bowl and the bowl depth are significantly less than those of C60 (8.9° and 0.92 Å for 1 vs 12.4° and 1.43 Å for 3), enabling a dynamic bowl-to-bowl interconversion process. The C5v bowl structure converts to another C5v bowl through a D5h flat transition state structure, with an experimental estimate of between 11.1 and 11.5 kcal/mol for the interconversion barrier (EBtB).66 Previous investigations across a large variety of corannulene derivatives enabled the demonstration of a strong structure−energy correlation between the bowl depth and EBtB, which follows a quartic function. The observed differences in activation energies of corannulene derivatives are attributed to ground-state perturbations by destabilization/stabilization of the bowl-shaped structure through repulsive/attractive interactions. As a result, this property is highly method-dependent and represents a good test of the methodology discussed in this paper. Table 3 summarizes the values of EBtB for corannulene as a function of wave function type in conjunction with the Def2-TZVP basis set. The results can be compared to the experimental estimate as well as the theoretical reference value determined

Table 2. Comparison of CPU Times for the Nine DSD-DFs and Corresponding RI-DSD-DFs (Values Are Expressed as DaysHours:Minutes) Methods

HB9

CT7/04

DI9

ADIM5

IDISP4

PPS11

RG21

overall

DSD RI-DSD ratio

6-22:39 0-2:51 58.3

1-13:18 0-0:47 47.1

5-3:27 0-2:13 55.5

32-9:9 1-2:51 28.9

10-6:48 0-5:49 42.4

98-1:37 1-9:21 70.6

1-1:31 0-0:13 48.6

155-10:33 3-0:25 51.5

J

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 3. Comparison of Bowl-to-Bowl Interconversion Barriers (EBtB), Complexation Energies for the X:C20H10 Complexes (EX, X = H2, 2H2, Cl2, Br2), and Random (ERAND) and Eclipsed (EECL) 1−1 Aggregation Energies (in kcal/mol) as Functions of Wavefunction Type Together with the Def2-TZVP Basis Set complexation energiesa

bowl dynamics

aggregation energiesa

wave function type

EBtB

EH2

E2H2

ECl2

EBr2

ERAND

EECL

HF BP86 BPW91 PBELYP PBE PBEP86 PBEPW91 OLYP B3LYP ωB97 BP86+D2 BP86+D3 B97-D PBE+D2 PBE+D3 B3LYP+D2 B3LYP+D3 DSD-PBEPBE DSD-PBEP86 DSD-BP86 RI-DSD-PBEPBE RI-DSD-PBEP86 RI-DSD-BP86 MP2 SCS-MP2 experiment CCSD(T)/cc-pVQZ//B97-D/cc-pVQZc

8.64 9.97 9.73 10.45 10.21 10.49 10.23 9.52 9.98 9.85 10.91 11.18 10.91 10.88 10.77 10.91 11.06 10.34 10.47 10.46 10.34 10.47 10.46 10.89 9.92 11.1−11.5b 10.48

0.63 0.70 1.07 −0.49 −0.15 −0.57 −0.18 1.19 0.44 −0.96 −0.50 −0.78 −0.77 −1.00 −1.07 −0.76 −0.84 −0.90 −0.96 −0.72 −0.90 −0.96 −0.72 − − −1.3892 −

1.23 1.48 2.27 −1.10 −0.38 −1.29 −0.45 2.53 0.90 −2.24 −1.00 −1.68 −1.62 −2.16 −2.36 −1.58 −1.82 −1.88 −1.88 −1.47 −2.01 −2.01 −1.47 − − −2.7692 −

2.66 −2.40 −1.28 −3.93 −3.61 −4.78 −3.65 0.73 −1.17 −2.91 −5.70 −6.64 −5.13 −5.96 −6.00 −4.47 −4.96 −4.82 −5.14 −4.89 −4.82 −5.15 −4.90 − − − −

2.54 −3.99 −2.73 −4.83 −5.17 −6.44 −5.17 0.26 −2.17 −4.12 −8.50 −9.26 −7.82 −8.38 −8.01 −6.67 −6.89 −7.05 −7.31 −7.22 −7.05 −7.32 −7.22 − − − −

2.82 2.69 4.59 −1.01 −0.24 −2.30 −0.35 6.74 2.14 −6.15 −9.49 −11.48 −10.62 −8.94 −9.14 −10.04 −10.11 −10.48 −10.40 −10.29 −10.48 −10.40 −10.30 −14.11 −10.72 − −

20.00 11.50 16.19 8.33 7.07 2.30 6.89 27.84 13.74 −9.78 −17.58 −19.33 −16.61 −13.70 −11.16 −15.34 −12.75 −21.05 −21.00 −21.27 −21.06 −21.01 −21.79 −35.59 −25.21 −16.193 −

a Interaction energies were determined as the sum of the energy of the separated monomers in the relaxed configuration of the optimized complex. The complexes used for the complexation energies were optimized at the B97-D/Def2-TZVPP level of theory, and the complexes used for the aggregation energies were optimized at the B97-D/6-311G(2d,p) level. bBecause of symmetry, the value for unsubstituted corannulene was estimated from a set of substituted forms.66 cAs shown in the Supporting Information, benchmarking across a variety of functionals revealed a high level of reliability of B97-D for structural data.

given their similar topological characteristics.92 In particular, corannulene derivatives hold considerable interest for use in hydrogen storage devices92,99,100 or as sensor elements in molecular junctions.65,96 Hydrogen is of particular interest as a light-weight and relatively high abundance synthetic fuel with an environmentally benign oxidation product. However, hydrogen storage remains a problem,99,101 making efforts toward the development of new materials for effective and safe storage an active area of research.101 In the complexation of molecular hydrogen, the relatively large permanent dipole of corannulene (∼2.1 D) creates an induced dipole in the absorbed H2, enhancing its physisorption onto the corannulene cap. Figure 10C shows the B97-D/Def2-TZVP-optimized H2:C20H10 complex, with the H2 approaching in a perpendicular fashion onto the convex cap of corannulene at a distance of 2.92 Å. Hydrogen uptake experiments have determined the hydrogen complexation energy to be on the order of 1.4 kcal/mol,92 supporting a weak physisorption process. Complexation of an additional H2 (Figure 10D) was shown to essentially double this value (−2.76 kcal/mol), and further studies confirmed the additive property of hydrogen complexation.92 Studies also indicated a slight preference for complexation with the concave face relative to the convex face at low coverage, while at increased coverage hydrogen complexation with the convex face was preferential.99,101

at the CCSD(T)/cc-pVQZ//B97-D/cc-pVQZ level of theory. The results show the high sensitivity of this property to wave function type, ranging between 8.64 and 11.2 kcal/mol. All of the standard DFs are outside the target accuracy of 0.5 kcal/mol with respect to both the experimental estimate and the CCSD(T) reference value, except for PBELYP and PBEB96. Inclusion of an empirical correction then brings the predictions back within the target accuracy. While the “best” predictions are obtained with B3LYP+D3 and BP86+D3, it is also clear that the D3 empirical correction often overshoots the reference values, so these “best” predictions could be fortuitous. The DSD-DFs also provide reasonable estimates for this property, and the corresponding RI-DSD-DFs are exactly in line with their counterparts, with significant savings in computational time. Application II: H2, Cl2, and Br2 Corannulene Complexation Energies. Nanostructured materials have been shown to be ideally suited for sensor applications with their large surface to volume ratios, making them sensitive to complexation of key target molecules.94,95 Upon functionalization, carbon-based materials of varying dimensionality, although chemically inert in pristine form, can become highly sensitive as detectors for small gas molecules.95,96 Physisorption studies of molecular hydrogen on single- and multiwalled carbon nanotube structures97,98 have motivated the use of curved aromatics as reasonable alternatives K

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 10. Structures of (A) eclipsed-stacked and (B) random-stacked 12 and side views of (C) H2, (D) 2H2, (E) Cl2, and (F) Br2 complexation with the convex side of corannulene.

substituent interactions can also contribute to variations in binding energies. Table 3 compares interaction energy data for the random (ERAND) and eclipsed (EECL) configurations of corannulene (Figure 10A,B) as functions of wave function type in conjunction with the Def2-TZVP basis set. Across all of the standard noncorrected GGA DFs and hybrid DFs in the set, none come close to a reasonable value for either interaction energy, and most do not even predict a stabilizing interaction. Even ωB97, which was seen to be one of the better range-separated DF choices for weak interactions in the performance analysis, significantly underestimates the interaction energies. Dispersion correction with either D2 or D3 brings the energies back in line and shows that EECL is about twice ERAND. This makes sense on the basis of the larger vdW surface benefit with columnar stacking (Figure 10). The degenerate highest occupied molecular orbital (HOMO) of a 1−1 stack is analogous to an antisymmetric filled−filled interaction seen in benzene dimers, but as a result of the bowl shape there is some in-phase overlap of the orbital contributions coming from the hub of one bowl and the rim of the next. As concluded from the performance analysis, DSD-PBEPBE is shown to be on par with CCSD(T) and thus serves as a comparable reference. Comparison of the DSD values with those for the dispersion-corrected DFs for ERAND would suggest only B3LYP+D2(D3) to be within the target accuracy. However, for the eclipsed configuration (EECL), none of the dispersioncorrected DFs are within the target accuracy. The best estimate is provided by BP86+D3, which underestimates EECL by ca. 1.7 kcal/mol with respect to the theoretical references. Although most of the functionals capture the relative interaction energy stability of EECL and ERAND, none predict energies within 0.5 kcal/mol of the reference values. This result will be further elaborated in the next section. Once again, the performance benefit of the RI-DSD-DFT is brought to light. At a fraction of the computational cost (8% of the DSD-DFT cost), the interaction energies of both the random (ERAND) and eclipsed (EECL) configurations are predicted by RI-DSD-PBEPBE and RI-DSD-PBEP86 to be within 0.08 kcal/mol of those for the respective DSD-DFs, consistent with conclusions from the validation sets. Application IV: Pentaindenocorannulene Complexes and Aggregates. Supramolecular order determines to a large extent the performance of curved aromatic building blocks within a device. Considerable effort has been devoted to obtaining robust nanostructures through the aggregation of organic building blocks.78,105−109 The delicate interplay among various intermolecular and interstack interactions, electrostatic interactions (dipole−dipole attractions, repulsion between surfaces), and

Evaluation of the data in Table 3 shows a wide variation in the predicted energies for complexation of H2 with corannulene. Most of the standard uncorrected DFs give values with the wrong sign and well outside the target accuracy. Addition of a dispersion correction brings the values back in line, but only PBE with either the D2 or D3 correction is within the target accuracy. Inclusion of the perturbative correction provides complexation energies within the target accuracy. As seen in the literature studies, across all of the DFs, the complexation energy for two H2 molecules is essentially double that for a single H2 molecule. In the complexation of two H2 molecules, only the PBE+D3 energy is within 0.5 kcal/mol of the literature experimental value. Relatively unexpected is the slightly larger deviation between the DSD-DFs and the respective RI-DSD-DFs for this property. Given the results for H2, the question arises as to the sensitivity of the interaction between corannulene and a target molecule. Here we investigated the interaction energies across the series H2, Cl2 and Br2 with corannulene, which should track roughly with the degree of polarizability, or equivalently, the heat of vaporization. The heats of vaporization for H2, Cl2 and Br2 are 0.904, 20.41, and 29.96, respectively.102 Table 3 summarizes the complexation energies for the three complexes H2:C20H10, Br2:C20H10, and Cl2:C20H10 as functions of wave function type in conjunction with the Def2-TZVP basis set. The optimized geometries for the complexes with Br2 and Cl2 target molecules (Figure 10E,F) are very similar to that of the H2:C20H10 complex discussed above, with interaction distances of 2.85 and 2.84 Å, respectively. As with H2:C20H10, the interaction involves a similar dipole-induced physisorption process. The Cl−Cl and Br−Br distances (2.07 and 2.39 Å, respectively) are still very similar to the distances found in their isolated molecular forms (∼1.99 Å and ∼2.28 Å, respectively), further supporting the physisorption characterization. The interaction energies across the set are predicted to be roughly −1, −5, and −7 kcal/mol for H2:C20H10, Cl2:C20H10, and Br2:C20H10, respectively, indeed correlating with their respective ΔHvap values.102 Application III: Corannulene Crystal Packing and Stacking Features. The curved bowl shape structure of 1 might suggest columnar aggregation (Figure 10A). However, the structure is observed to pack in the crystal phase without any columnar order (Figure 10B). While corannulene does not display any columnar order in the crystal, most functionalized derivatives do show columnar arrangement along a particular stacking axis. Alternatively, assembly of corannulene or associated derivatives onto metallic surfaces reveals stacking tendencies with multiple monolayers.103,104 The preference for stacked bowl-in-bowl packing tends to arise from dispersive interactions. The spatial extent of the anisotropy of the molecular wave functions governs the relative orientation. Dispersive interactions due to L

DOI: 10.1021/acs.jctc.7b00220 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation dispersive interactions governs the arrangement of the molecules in a material. Recent attention has been paid to the indenocorannulene family as attractive targets for next-generation photovoltaics and other optoelectronics applications.110,111 The significant curvature of these components favors the formation of columnar stacking in the crystalline phase, confining the electronic states anisotropically on the nanoscale and imparting important physical characteristics. In this indenocategory of structures, pentaindenocorannulene (C50H20, 2; Figure 11) is a deep-bowl polynuclear aromatic hydrocarbon.

Table 4. Interaction Energies (in kcal/mol) for B97-D/ 6-311G(2d,p)-Optimized Geometries of the Pentaindenocorannulene 1:1 Aggregate (EAGG) and the 1:1 Complex with C60 (ECOMP) (Figure 11 A,B) for Various Wavefunction Types in Conjunction with the Def2-TZVP Basis Set

Figure 11. Pentaindenocorannulene (A) aggregate and (B) complex with C60.

wave function type

EAGG

ECOMP

OLYP HF BPW91 B3LYP BP86 PBELYP PBE PBEPW91 PBEP86 ωB97 B3LYP+D3 B3LYP+D2 PBE+D3 PBE+D2 B97-D BP86+D2 BP86+D3 RI-MP2 RI-SCS-MP2 RI-DSD-PBEPBE RI-DSD-PBEP86 RI-DSD-BP86

69.75 49.32 38.19 33.76 27.10 23.36 17.77 17.45 6.69 −23.14 −36.42 −44.95 −46.45 −47.32 −48.50 −51.60 −53.24 −104.52 −75.11 −61.95 −61.26 −64.83

55.06 40.99 28.67 26.69 20.06 18.94 13.41 13.21 4.87 −17.30 −28.34 −37.07 −30.97 −38.45 −40.84 −43.70 −42.29 −92.97 −67.09 −53.56 −53.19 −56.52

illustrative applications. Typically, compared with traditional evaluation of the integrals, the RI approximation deviates by