Structure and binding in crystals of cage-like molecules: hexamine and platonic hydrocarbons Kristian Berland and Per Hyldgaard arXiv:1010.1487v1 [physics.chem-ph] 7 Oct 2010 Department of Microtechnology and Nanoscience, MC2, Chalmers University of Technology, SE-41296 Göteborg, Sweden (Dated: October 8, 2010) In this paper, we show that first-principle calculations using a van der Waals density functional (vdW-DF), [Phys. Rev. Lett. 92, 246401 (2004)] permits determination of molecular crystal structure. We study the crystal structures of hexamine and the platonic hydrocarbons (cubane and dodecahedrane). The calculated lattice parameters and cohesion energy agree well with experiments. Further, we examine the asymptotic accounts of the van der Waals forces by comparing full vdW-DF with asymptotic atom-based pair potentials extracted from vdW-DF. The character of the binding differ in the two cases, with vdW-DF giving a significant enhancement at intermediate and relevant binding separations. We analyze consequences of this result for methods such as DFT-D, and question DFT-D’s transferability over the full range of separations. PACS numbers: 61.50.Lt, 61.66.Hq, 71.15.Mb I. INTRODUCTION Understanding supramolecular structure and interactions is essential for understanding many biological processes.1 Biological and other supramolecular complexes as polymers and overlayers are sparse matter, that is, they contain low electronic density in essential regions. The general lack of order in these systems prohibits precise measurements of atomic structure and therefore challenges development of theoretical methods. In turn, this makes transferable first principles schemes attractive; an accurate account of simple periodic structures permits accurate characterization to be made and reliable conclusions to be drawn for more complex (non-periodic) systems. Molecular crystals and simple polymer crystals2 are ideal testing grounds for applications of first-principle descriptions of sparse matter. Unlike most sparse matter, they constitute ordered systems and therefore leads to unambiguous comparison between theory and experiments. The dispersion forces underpin the cohesion of sparse matter. Modeling of sparse matter at the electronic level, therefore, requires that we take these effects into account. Traditional implementations of Density-Functional theory (DFT), which parameterizes the density functional within the local density approximation (LDA)3 or the generalized gradient approximation (GGA)4,5 have allowed for routine modeling of matter with dense electron distributions. However, the lack of non-local correlation and hence dispersion forces in these implementations, has effectively inhibited widespread use for general sparse matter systems. Computational methods able to predict molecular crystal structure and stability are also of value in the pharmaceutical industry, for example in mapping out competing crystalline phases for the fabrication of pharmaceuticals.6 The electronic structure and response are often required to understand and compute properties of technological significance. Consequentially, a range of developments has aimed to extend density functional approximations with an account of dispersive interactions and thus capability to address sparse matter challenges. Such approaches include non-local functionals7–12 , intermolecular perturbation theory13,14 , and the addition of a semi-emperical atom-based account of the dispersion forces (DFT-D).15–18 DFT-D has also been used to characterize molecular crystals.19,20 We study molecular crystal using a van der Waals density functional (vdW-DF).10 The correlation part of this functional depends non-locally on the density and accounts for dispersion forces. vdW-DF has been used in a series of first-principle studies of sparse matter,21 yielding for example new insight to the twist of DNA,22 nanotube bundles,23 water hexamers,24 metal-organic frameworks,25 and the binding mechanisms at organicmetal interfaces.26–29 In this paper, we first demonstrate that vdW-DF permits structural determination of molecular crystals and second asses the character of the attractive dispersion interactions. For the first, we calculate the lattice parameters for crystals of the cage-like hexamine, dodecahedrane, and cubane molecules. For the second, we compare an atom-centered asymptotic-1/r6 approximation, frequently encountered in force-field methods6,30,31 and in DFT-D, with the correlation energy provided by vdWDF. We also analyze consequences for DFT-D and question its performance at intermediate separations. The particular choice of molecules was motivated in part by their aesthetic appeal and in part for the property that their symmetric geometries permits a simple assignment and analysis of the strength of the asymptotic interactions even when expressed in terms of an atom-pair basis. This paper has the following plan: The next section presents the molecules hexamine, cubane, and dodecahedrane and their experimental crystal structures. The third section deals with the computational details of vdW-DF. The fourth gives the results for lattice parameters and bulk modulus. In the fifth section, we compare the asymptotic 1/r6 -form of van der Waals interactions with the full non-local correlation for different molecules, and the vdW-DF potential energy curve with DFT-D cal- 2 FIG. 1: The molecular and crystal structure of hexamine, dodecahedrane, and cubane. The upper panels show the respective molecular structures. The mid panels give schematics of the crystal structures, where the molecular symmetry and orientation are highlighted by the use of a tetrahedron, dodecahedron, and a cube in place of the molecules. The lower panel shows, for hexamine and dodecahedrane, facets of the crystal structure and for cubane the orientation of the cubane molecule in the unit cell. Hexamine has a bcc unit cell, dodecahedron an fcc unit cell, and cubane a rhombohedral unit cell, with equal angles α between the lattice vectors. culations for a cubane dimer. The final section holds our conclusions. II. PLATONIC MATTER Plato asserted that the basic building blocks of Nature were five geometrical structures, today known as platonic solids.32 These are the five convex polyhedra with all faces, edges, and angles congruent. On the molecular level, two of these solids have synthetic hydrocarbon (compounds of only carbon and hydrogen) analogues: cubane (C6 H6 ), which corresponds to the cube, and dodecahedrane (C20 H20 ), which corresponds to the twelve-faced dodecahedron. The fascination mathematicians held for these geometrical structures since antiquity echoed in the more recent struggle to synthesize their hydrocarbon representatives. The highly strained bonds posed the main challenge. Obstacles were eventually overcome and in 1964 cubane was synthesized33 followed by dodecahedrane in 1978.34 A third, tetrahedrane exists only within a larger chemical structure35 and no crystallographic characterization exists. Tetrahedrane does therefore not represent a good testing ground for vdWDF and we instead study the molecular crystal of hexamine (C6 H12 N4 ). Although not a hydrocarbon, it shares several features with the platonic hydrocarbons: it is organic, its non-hydrogen atoms form a cage, and this cage has the symmetry of platonic solid tetrahedron. The top row of Fig. 1 shows the molecular structures (from left to right) of hexamine, cubane and dodecahedrane, ordered according to the ascendancy of their platonic analogues. The second and third row of Fig. 1 show the crystal structures. The first column shows hexamine, the mid, cubane, and the final, dodecahedrane. The crystal structure of hexamine forms a bcc with a single molecule in each unit cell and I4̄3m symmetry. Hexamine has been used as a model system in numerous studies36–38 and the crystal structure was determined as early as in 1923.39 The crystal structure of cubane is rhombohedral with equal external angles and a single molecule in each cell. The molecule is oriented according to the R3̄ space group; 3 a rotation about the [111] axis specifies the configuration of the molecule.31 The crystal structure of dodecahedrane is fcc, also with a single molecule in each cell and Fm3̄space group.40 Hexamine finds numerous industrial uses.41 For instance, it serves as a component in fuel tablets and as an antibiotic.42 The platonic hydrocarbons find only hypothetical applications; cubane has been identified as a potential high-energy fuel and explosive.43 In addition to being well-studied crystals, especially hexamine and cubane, these crystals make attractive testing grounds for vdW-DF for two more reasons. First, their simple structures allow for brute-force determination of lattice parameters. Within the crystal symmetry, this determination corresponds to mapping out the potential energy of a single parameter for hexamine and dodecahedrane, and three for cubane. The brute-force approach facilitates post-processing analysis, such as computation of bulk modulus,44 and makes it easier to evaluate how the choice of exchange in DFT affects crystal structure and cohesion energy. Second, the high symmetry of the molecules reduces the large set of atom-to-atom C6 coefficients to only a few equivalent values. This reduction simplifies the comparison of full vdW-DF calculations with atom-based asymptotic accounts of the van der Waals forces. III. COMPUTATIONAL DETAILS The crystal structure and bulk modulus of hexamine, cubane, and dodecahedrane are determined with DFT using a nonlocal density functional called vdW-DF. The traditional semi-local GGA for exchange-correlation provides accurate bond lengths and charge density n(r), but fails to capture correlated motion of separated electrons: the long-range dispersion forces. vdW-DF includes these correlations and can therefore account for the structure and cohesion of sparse matter. Since details of the functional and its implementation are given elsewhere,10,21,45–47 we focus mostly on computational steps specific for determination of crystal structures. The non-local correlation of vdW-DF takes the form of a double-space integral, Z Z 1 nl Ec [n] = dr′ n(r) φ(r, r′ )n(r′ ) , (1) dr 2 V0 V over an interaction kernel φ(r, r′ ). Here V0 denotes the central unit cell and V (formally) the entire space. The kernel can be tabulated in terms of two parameters d and d′ related to the local response q0 (r) and spatial separation |r−r′ | by d = q0 (r)|r−r′ | and d′ = q0 (r′ )|r−r′ |. The remaining part of the exchange-correlation functional of vdW-DF consists of the exchange part of revPBE48 and the correlation of LDA: vdWDF Exc = EcLDA + ExrevPBE + Ecnl [n] . (2) The total energy functional of vdW-DF, E vdWDF [n], also includes the standard elecrostatic and kinetic-energy terms within the Kohn-Sham scheme.49 It is convenient to write this energy in terms of a semi-local part E0 [n] containing all but the non-local correlation, so that E vdWDF [n] = E0 [n] + Ecnl [n]. For input charge density n(r), we use the result of semi-local calculations with the PBE5 flavor of GGA. We will refer to calculations with the PBE flavor of GGA as DFT-GGA. The charge density could also have been obtained within fully selfconsistent vdW-DF.45,50 However, the two-step non-self consistent procedure introduces only a slight approximation. Previous studies have documented that for systems with small charge transfer, binding energies of non-self consistent vdW-DF only differ by tiny amounts from fully self-consistent energies.45,51 To speed up evaluation of the non-local correlation, we introduce a radius cutoff based on the decay of van der Waals forces at large separations. With this cutoff, the kernel takes form Z Z 1 nl Ec [n] ≈ dr′ φ(r, r′ )n(r′ ) . (3) dr n(r) 2 V0 |r−r′ |<R In the above expression, we see that the CPU-cost for evaluating the non-local part of vdW-DF goes as R3 O(V0 ). Thus, for large or periodic systems the computational costs increase linearly with system size. To cut computational costs further, we introduce an extra radius cutoff corresponding to the separation between dense and sparse sampling of the charge-density grid. We note that Ref. 50 reports a more elaborate scheme which considerably reduces CPU-costs, yet our simple measure was sufficient for our non-self consistent calculations as the underlying DFT-GGA calculations of electronic density dominated time consumption. The use of revPBE for exchange in vdW-DF was motivated by the fact that this exchange functional excludes unphysical binding effects at large distances.46,52 For a range of systems, vdW-DF overestimates binding separations.21,46 Several studies indicate that this discrepancy can be attributed to the details of the exchange functional.21,51,53–55 Puzder et al53 have demonstrated that replacing Hartree-Fock exchange with revPBE improves binding separations for benzene dimers. Gulans et al54 have shown that for a selected range of molecular complexes the PBE exchange functional improves binding energies. We furthermore illustrate the sensitivity to exchange by including results based on use of an alternative vdW-DF(PBE), where revPBE exchange has been replaced by that of PBE. We do not argue for replacing vdW-DF with vdW-DF(PBE), instead we simply explore consequences of a different account of exchange.55 To calculate crystal parameters and cohesion energy, we minimize the potential energy. In many respects this vdW-DF structure determination is similar to those in Refs. 2,23,47,56. The potential energy is given by the difference between the total energy of the full crystal and 4 a reference energy for a system of isolated molecules TABLE II: The vdW-DF prediction of lattice parameters, cohesion energy and bulk modulus for the crystals of hexamEcoh (a, {αθ}) = E (a, {αθ})−E (a → ∞, {αθ}) . ine, cubane and dodecahedrane compared with an alternative (4) vdW-DF(PBE) based on PBE-exchange and with experimental values. The experimental lattice parameters are based on In the above equation the curly brackets are specific for low temperature measurements, except for dodecahedrane. the cubane crystal as it depends on three rather than one vdWDF vdWDF parameter. In the reference calculation corresponding to the reference energy, E vdWDF (a → ∞, {αθ}), the semilocal part is obtained somewhat differently from the part containing the non-local correlation. For the former, we effectively isolate the molecules in our periodic boundary calculations, by using a unit cell of doubled size in all directions. This measure secures negligible charge overlap between the molecules in the supercells. For the latter (nonlocal) part we restrict the integral of Eq. 3 to the central supercell to avoid coupling between the enlarged unit cells. Hence, only non-local correlations within the molecule contribute to the reference energy. To enhance accuracy in the evaluation of the non-local part of the potential energy, we systematically cancel a small, but noticeable, grid dependence in the evaluation of Eq. 3. This cancellation is performed by making sure to use the same FFT grid spacing in the reference calculation as in the main calculation. Furthermore, we make sure to place the isolated molecules in the same relative configuration to the underlying grid in the reference calculation as in the main calculation. We thus perform an additional reference calculation for every molecular configuration investigated. These measures have been used to secure a high accuracy of the non-local part of vdWDF in several earlier studies.2,46,47,56 We map the potential energy landscape by varying the lattice parameters of the molecular crystals within the experimental crystal symmetry. The stiff cage-molecules allow us to keep the internal coordinates of the molecules frozen for all configurations. The molecular structures are determined in isolation using the PBE flavor of GGA. The resulting bond lengths will be compared with experimental data in the next section to verify the utility of conventional DFT-GGA for the internal structure of strained molecules. TABLE I: Experimental and calculated bond lengths of hexamine, cubane, and dodecahedrane. The calculations were done with the PBE flavor of GGA. l denotes the C-C bond length for cubane and dodecahedrane and the C-N bond length for hexamine. lCH denotes the carbon-hydrogen bond length. Parameter l[Å] lexp [Å] lCH [Å] exp lCH [Å] a Ref. 58. b Ref. 59. c Ref 60. Hexamine 1.472 1.476a 1.101 1.088a Cubane 1.566 1.562b 1.095 1.097b Dodecahedrane 1.549 1.544c 1.100 - Parameter vdW-DF vdW-DF (PBE) Exp. a (Å) Ecoh (eV) B0 (GPa) 7.14 -1.0127 10.0 6.93 -1.427 14.0 6.910a -0.827b 7.0c a (Å) α θ Ecoh (eV) B0 (GPa) 5.45 73 47.5 -0.77 7.2 5.25 72.5 46.5 -1.15 14.8 5.20d 72.7d 46 d -0.857e - a (Å) Ecoh (eV ) B0 (GPa) 10.92 -1.46 12.2 10.56 -2.06 18.6 10.60f Hex. Cub. Dod. a Ref. - 58. b Ref. 61. c Ref. 62. d Ref. 58. e Ref. 63. f Ref. 64. The electronic-structure calculations rely on the planewave code DACAPO 57 using ultra-soft pseudopotentials. In combination with a separate reference calculation as previously discussed, we secure the convergence of the non-local correlation by specifying an FFT-grid spacing less than 0.13 Å. This spacing leads to an effective planewave energy cutoff of at least 500 eV. For all crystals, we set the Monkhorst-Pack k-sampling to 4 × 4 × 4. IV. STRUCTURE DETERMINATION A. Molecular structure Table I shows the calculated molecular structures. It demonstrates DFT-GGA can account for the intramolecular bonding even in the highly strained cubane molecules. The resulting bond lengths differ by less than 1 % from the experimental values. In first principle studies of molecular crystals, accurate determination of lattice parameters require accurate account of constituent molecules. Unlike many calculations, where determination of crystal structure often starts from the structure of the individual molecules, most experiments resolve the molecular structure by looking at the diffraction pattern of a full molecular crystal. For our purposes of testing the first-principle vdW-DF method, this is fortunate, as efforts to characterize molecular structures also generate an abundance of experimental data on molecular crystal structures. 5 FIG. 2: The potential energy curves for the crystal of hexamine and dodecahedrane (left panel) and corresponding contour plots for the crystal of cubane (right panel). The left panel shows potential energy where the curves are normalized separately so that the experimental lattice parameter and the calculated cohesion energy equals unity. The solid curve represents dodecahedrane and the dashed hexamine. The difference in the optimal a/aexp value can be attributed to the experimental lattice parameter being measured at low temperature for hexamine, but not for dodecahedrane. The right panel shows spline-interpolated contour plots of the two-dimensional intersection of the three-dimensional potential energy landscape. The main figure corresponds to the optimal value of the internal angle θ, while the insert corresponds to the optimal value of the unit cell length a. The pronounced asymmetry of the curves and contours reflect the hard wall provided by Pauli repulsion. B. Crystal structure Figure 2. shows the binding curves and contours found by varying the molecular crystals’ lattice parameters identified in Fig. 1. The simple crystal structures of hexamine and dodecahedrane give a one-dimensional potential energy landscape, while the cubane crystal has a three dimensional one. In the right panel of Fig. 2, we display the αa and the αθ intersections, for the optimal values of θ and a. The curves exhibit a pronounced asymmetry around their minimum. This asymmetry arise because the van-der Waals attraction is much softer than the kinetic-energy repulsion. Table. II contains the calculated results and experiment values. Standard vdW-DF performs well both for lattice parameters and cohesion energy. If not directly available, we obtain experimental cohesion energies by correcting for gas phase and vibrational contributions to the enthalpy of sublimation, using the method described in Refs. 65,66. The bulk modulus is obtained with use of polynomial interpolation according to the scheme of Ziambaras and Schröder.44 The required polynomial were constructed using data from selected one-dimensional deformations.67 For hexamine, where the experimental bulk modulus is available, the computed value show fair agreement with the experimental value. There is also a trend for bigger cage molecules to have a larger bulk modulus. We attribute this trend to the fact that for bigger molecules a smaller relative part of the unit cell consists of soft intramolecular regions. Therefore, as the relative unit cell dimensions change, the distance between (stiff) molecules changes more for big molecules than for small molecules. For all three crystals, we find unit-cell volumes somewhat larger than the experimental ones. Similar overestimations have also been encountered in previous studies.21,46,55 The alternative choice of PBE as exchange functional, vdW-DF(PBE), influences results substantially. On one hand, it improves lattice parameters, almost to level of standard DFT for intra-molecular bonds. On the other hand, the value of cohesive energy and bulk modulus worsens. These results signal that the main discrepancy between experiments and vdW-DF stems from the specific form of semi-local exchange.21,53–55 We note that an ab initio study of cubane has previously been performed at the LDA level.68 Being a local functional, LDA has no physical basis for the van der Waals binding that provides the cohesion of this molecular crystal. The spurious LDA binding arises from an unphysical accounts of exchange.46,52 Once the molecular crystals are investigated with DFT-GGA, which has an improved account of exchange, the binding essentially vanishes.68 V. ASYMPTOTIC PAIR POTENTIALS VERSUS NON-LOCAL CORRELATION The asymptotic van der Waals interactions between two atoms or molecules goes as the simple power law C6 /r6 , where C6 gives the strength of the interaction. This familiar result can be derived from second-order perturbation theory,69 or from an analysis of the shifts in the zero-point motions of the electron.70 A common strategy in force-field methods and empirical extensions of DFT 6 than to the surface atoms. Within vdW-DF, Kleis et al23 demonstrated that for interactions in nanotube bundles, the force stems primarily from the electron tail around the nanotube, and that as the tubes get closer, higher order moments dominate over the asymptotic interaction. Part of this enhancement can be interpreted as an imageplane effect. Here we investigate whether an atom-centered 1/r6 form is a good approximation for the non-local correlation of molecular dimers. As our argument is based on the asymptotic vdW-DF account of van der Waals forces, we also discuss for these molecules the accuracy of the asymptotic account. In the first subsection, we compare the non-local correlation of vdW-DF for dimers of hexamine, cubane, and dodecahedrane with its corresponding atom-centered pair potentials. We document a significant enhancement of the non-local correlation at short (binding) separations and at intermediate separations, one to three Ångström beyond typical binding separations. In the second subsection, we discuss the accuracy of our C6 coefficients. In the third, we investigate consequences of this result, in particular for the use of DFT-D. We will argue that although standard DFT-D methods can provide good descriptions of both short and asymptotic separations, their asymptotic atom-centered form does not describe the enhancement of correlation at intermediate separations exhibited by vdW-DF. FIG. 3: Comparison between the non-local correlation of vdW-DF and atomic pairs potentials generated with asymptotic vdW-DF for a dimer of hexamine, dodecahedrane, and cubane (from top to bottom) in configurations corresponding to nearest neighbors in their respective crystal. The dashed curve gives the APP-vdWDF result, while the solid curve gives non-local correlation of vdW-DF: ∆Ecnl (d) = Ecnl (d) − Ecnl (d → ∞). The horizontal axis gives the difference from the vdW-DF crystal binding separation d0 , which respectively takes the values 6.18, 5.45 and 7.72 Å. The difference between the two curves demonstrates the enhancement of non-local correlations over the vdW-DF asymptotic atombased account at relevant binding separations and intermediate separations. A. Asymptotic vdW-DF For the the asymptotic van der Waals forces, the C6 coefficient between two fragments, A and B, can be computed from the general formula Z 3 ∞ C6AB = du αA (iu)αB (iu) , (5) π 0 where α(ω) is the polarizability of the fragment. In order to calculate the C6 coefficients we approximate α(ω) with the local external-field susceptibility of vdW-DF, χvdW−DF (ω, r) = A is to adopt such an asymptotic form at all separations in terms of atom-centered pair-potentials (APP). However, we can not take for granted that the asymptotic behavior should hold for separations closer to that of intramolecular binding. On the contrary, because van der Waals forces arise from correlated motion of electrons and not from the atomic nuclei, several mechanisms affect and enhance the interaction at short and intermediate separations: higher order moments contribute, polarizability changes as charge is distorted and finite-k dispersion of the electronic response becomes important. Zaremba and Kohn71 considered adsorption of noble atoms on surfaces and documented a significant enhancement of dispersion energies over an atom-centered account; their asymptotic 1/d3 form use the distance to an image plane d, rather nA (r) 2 [9q0 (r)2 /8π] − ω 2 , and a polarizability given by Z vdW−DF αA (ω) = d3 r χvdW−DF (ω, r) . A (6) (7) To generate vdW-DF based atom-centered pair potentials (APP-vdWDF), P we first partition the full charge density n(r) = i ni (r) among the atoms of the molecules with aid of Bader analysis.72,73 Based on this charge partition, we calculate the atom-to-atom C6 coefficients using Eq. (5). Initially, this procedure generates N 2 different C6 coefficients for a molecular dimer of N atoms per molecule. For a dodecahedrane dimer, we get as much as 40 × 40 = 1600 coefficients. Fortunately, because of the high symmetry of the isolated molecules, 7 this number reduces to only three equivalent values for cubane and dodecahedrane: C6C−C , C6C−H , and C6H−H . For hexamine, the extra nitrogen atoms lead to six coefficients. As noise in the electronic density affects the value of the coefficients, we average over a large set of equivalent values to obtain the final values. TABLE III: Computed values of the C6 coefficients (Hartree atomic units) for different pairs of atoms within respective molecules, using the asympotic form of vdW-DF with charge density as partitioned with a Bader analysis. The figure also shows the molecule-molecule C6 coefficients obtained with the Andersson-Langreth-Lundqvist (ALL) scheme,74 and the molecule-molecule C6 coefficents obtained with the use of parameters given in DFT-D schemes.16–18 C6vdWDF C-C C-H H-H N-N N-C N-H mol-mol mol-mol mol-mol mol-mol mol-mol (vdW-DF) (ALL) (Wu) (Grimme) (Jurečka) hex 4.44 4.04 3.98 32.6 12.0 11.0 3470 3270 4340 4000 4790 cub 13.5 6.87 3.72 1990 1940 2600 2630 3120 dod 11.2 6.45 4.07 11300 10200 16300 16400 19500 The upper part of Table III shows the C6 coefficients calculated within asymptotic vdW-DF, both as partitioned according to the Bader analysis and as evaluated for the entire molecule. The calculated coefficients per atomic pair deviates much from a naive assignment of the full coefficient of the molecule according to the number of valence electrons of the underlying atom. In such a scheme the C-C coefficient would be 16 times larger than the H-H coefficient, while in fact, for cubane and dodecahedrane, it is merely three-four times stronger. This result can be attributed to the relative stronger response of the low-density regions surrounding the hydrogen atoms as q0 (r) ∝ [n(r)]1/3 (in the homogeneous limit), and these areas dominate Ecnl . The somewhat anomalous values for hexamine can be attributed to our Bader analysis scheme partitioning a significant portion of the charge density near the carbon atoms to the centrally located nitrogen atoms.75 For cubane and dodecahedrane the partitioning was similar, and the differences in atomistic C6 values show that they are influenced by their local enviroment. Having generating C6 coefficients appropriate for a comparison between the asymptotic account and the full correlation of vdW-DF, we study dimers of hexamine, cubane and dodecahedrane at different separations. We choose orientations that are given by the nearestneighbor configurations in the crystals. The total asymp- totic non-local correlation of APP-vdWDF reads vdWDF = Eapp XX i j C6ij , |ri − r′j | (8) where i and j label atoms in separate molecules of the dimer. Figure 3 shows results for the full non-local correlation of vdW-DF (solid curve) and that of APP-vdWDF (dashed curve). At large separations these two curves converge. In contrast, they differ significantly at relevant binding (short) separations and at intermediate separations. At these separations the non-local correlation is almost double as large as that of APP-vdWDF, for both hexamine and cubane, while for the biggest molecule, dodecahedrane, the difference is smaller. For all dimers, we also need to go to relatively large separations to recover asymptotic values for the non-local correlation. As vdW-DF is an ab initio functional, based on a set of exact sum rules,10 the strong enhancement of non-local correlations at short and intermediate separations indicates that the asymptotic form neglects important contributions. It also highlights that in constructing modified semi-empirical van der Waals functionals,11,12 a fitting of the asymptotic functional to C6 coefficients does not guarantee an accurate description at short binding separations. B. Comparison of C6 coefficients The lower part of of Table III, gives the molecular C6 coefficients calculated with asymptotic vdW-DF, the Andersson-Langreth-Lundqvist (ALL) scheme,74 and computed using the atomic coefficients of Wu et al,16 Grimme,17 and Jurečka et al.18 The ALL and vdWDF give coefficients smaller than that used in DFT-D schemes. We expect that Grimme and Wu provide good values for molecular C6 coefficients, because their underlying atomistic coefficients were fitted to reproduce a range of accurate molecular C6 coefficients calculated from experimental molecular polarizabilities (Ref. 16 and references therein). The coefficients of Jurečka,18 give somewhat larger molecular C6 coefficients. This comes from the use of Slater-Kirkwood average76 for C6 coefficients between different atomic species, while keeping those of Grimme for indentical atomi species (the C6 coefficients of Grimme are optimized for a different average. Asymptotic vdW-DF and ALL likely underestimates the C6 coefficients for these molecules; they are 20-30 % smaller than that used in DFT-D methods. For the similar ALL scheme, Ref. 77 reports an underestimation of C6 coefficients for larger molecules, in particular for benzene and C60. A difference between the ALL scheme and asymptotic vdW-DF is that for the former a hard cutoff accounts for plasmon damping, while for the latter, the local response q0 (r) provides a smooth cutoff. There is 8 good consistency between the two methods. Both methods also assume a local, scalar, relationship between the applied and the full electric field, which is an approximation for finite-sized objects.74 We speculate that this approximation contributes to the underestimation of C6 coefficient for the investigated, relatively large, molecules. TABLE IV: Atomistic C6 coefficients for cubane used in APPvdWDF (calculated with asymptotic vdW-DF and charge density partitioned with Bader analysis), and coefficients used in DFT-D methods. C6vdWDF C-C C-H H-H a Ref. b Ref. c Ref. vdW-DF 13.5 6.8 3.7 Wua 22.06 7.89 2.83 . . . Grimmeb 28.3 5.01 2.75 Jurečkac 28.3 8.82 2.75 16 17 18 Table IV shows atomistic C6 coefficients for cubane as calculated with asymptotic vdW-DF, and given by Wu,16 Grimme,17 and Jurečka,18 for use in DFT-D schemes. vdW-DF weights the relative response of the carbon less than that of the hydrogen, compared with the values of DFT-D.75 This property could relate to the above mentioned approximate treatment of electrodynamics. It could also relate to the carbon atoms being located somewhat inside the molecule, having a different local charge density and responding less to external fields than an atom on the exterior would; in contrast, DFT-D does not discriminate between atoms at different locations. C. Consequences for atom-based pair potentials: the binding curve of cubane The vdW-DF results presented in the first subsection shows that a simple asymptotic account only partially captures correlation effects at short and intermediate separations. As DFT-D use such an asymptotic form to describe non-local correlations, this result stands in apparent contrast to the many successful applications of DFT-D.16–18,78–80 To understand consequences of our result for methods such as DFT-D, we must first consider other effects that could contribute to the difference between APPvdWDF and the non-local correlation of vdW-DF. In vdW-DF, correlations are described by EcLDA + Ecnl and hence Ecnl also accounts for semi-local correlations.10,46 Second, we must consider the specific designs of actual DFT-D schemes, because these could counteract the lack of enhancement of non-local correlations. To describe exchange-correlation, DFT-D combines the asymptotic atom-centered form with a semi-local GGA account. It also introduces fitting parameters to be used in combination with a specific GGA flavors. The top panel of Fig. 4 shows that gradient corrections does not account for the difference between the full vdWDF and the APP-vdWDF results. For APP-vdWDF, we can combine the purely semi-local correlation of PBE with APP-vdWDF (in a new description APP-mod), vdWDF vdWDF EAPP−mod (d) = EAPP (d) + ∆EcPBE (d) − ∆EcLDA (d) , (9) to assess the magnitude of purely semi-local corrections relative to the difference between vdW-DF and APPvdWDF. In this APP account, we have introduced the LDA and PBE terms: ∆Ec (d) = Ec (d) − Ec (d → ∞). We focus our discussion on the cubane dimer. The lower thin solid curve gives APP-mod; the thick solid curve gives the non-local correlation of vdW-DF. The curves are shown as a function of the intermolecular distance d, with d0 indicating the vdW-DF binding distance in the crystal. At short separations, the thin curve lies closer to the thick curve than the corresponding APPvdWDF result (thick dashed curve). Thus, some of the difference between APP-vdWDF and the non-local correlation arises from a lack of semi-local correlation contributions in APP-vdWDF.10 However, even with this inclusion the difference is still significant, and at intermediate separations it remains undiminished. The middle panel of Fig. 4 details the effects of using the semi-empirical fitting of DFT-D and the larger C6 coefficients. We compare the binding curve for cubane obtained by use of vdW-DF with the binding curve obtained with DFT-D calculations. We select the schemes of Grimme17 and Jurečka,18 as these provide parameters for the PBE flavor of exchange-correlation, which is available to us.81 The DFT-D scheme of Grimme scales the strength of the dispersive interaction to a particular semi-local exchange-correlation to achieve good performance at short separations. For the PBE flavor of GGA, the cohesion energy of a dimer reads DFT−D,G PBE Ecoh = Ecoh +sPBE 6 X ij fPBE,G (|ri −rj |) ij C6ij , |ri − rj |6 (10) where s6 = 0.7. This scheme therefore sacrifices the asymptotic description in favor of the a good description of binding separations. The scheme of Jurečka instead adjust parameters of the damping function to the flavor of semi-local exchange-correlation, DFT−D,J PBE Ecoh = Ecoh + X ij ij fPBE,J (|ri − rj |) C6ij , |ri − rj |6 (11) and ensure that → 1 for large separations. The DFT-D scheme can therefore, in principle, describe both asymptotic and binding separations. However, the question remains on how it performs for systems where characteristic separations lies between these two limits. The binding curves of vdW-DF(revPBE) (the solid curve) and for the DFT-D scheme of Grimme(PBE) (upij fPBE,J 9 responding DFT-D scheme binds stronger than vdW-DF. The lower panel of Fig. 4 shows that the difference at intermediate separations comes primarily from the varying accounts of correlation. The solid curve shows the non-local correlation of vdW-DF. The dash-dotted upper (lower) curve shows DFT−D,G(J) Eapp = X ij fPBE,G(J) (|ri − rj |) ij C6ij |ri − rj |6 + ∆EcPBE (d) − ∆EcLDA (d) , FIG. 4: Comparison between different accounts of non-local correlation for a dimer of cubane. In the upper panel, the dashed curve gives APP-vdWDF. The thin solid curve gives the sum of APP-vdWDF and the correlation of PBE. The thick solid curve gives the non-local correlation of vdW-DF. In the middle panel, the thick curve gives the cohesion envdWDF ergy using vdW-DF, Ecoh (d), the upper (lower) dashdotted curve gives the DFT-D binding curve as given by DFT−D,G DFT−D,J Grimme, Ecoh (d) (Jurečka, Ecoh (d)). The lower panel gives the corresponding non-local correlation of vdWDF and Grimme (Jurečka). The thin solid curve gives the difference between exchange of revPBE and PBE. The horizontal axis gives the difference from the vdW-DF crystal binding separation d0 for cubane. per dash-dotted curve) and Jurečka(PBE) (lower dashdotted curve) indicates an underestimation of DFT-D at intermediate separations. The minimum at negative d − d0 shows that use of DFT-D would improve lattice constants over vdW-DF.55 The two DFT-D schemes give quite differing binding energies. Both results show that a binding energy at the same magnitude as vdW-DF can be achieved, even with an atom-based asymptotic form of the attractive potential. The energy of vdW-DF is significantly larger at intermediate separations than that of Jurečka despite that C6 coefficients of asymptotic vdWDF are underestimated (while those of Jurečka are likely to be somewhat overestimated) and despite that the cor- (12) which for DFT-D corresponds best to the non-local correlations provided by ∆Ecnl in vdW-DF. The figure also shows that the difference partly cancels, at short but not at intermediate separations, with the energy difference between the exchange flavors of revPBE and PBE, δEx = ∆ExrevPBE − ExPBE . Thus, for certain exchange functionals, adding an asymptotic atom-based account of non-local correlation can generate good binding values, yet our vdW-DF results indicate that this framework is not optimal for describing interactions at intermediate separations. Our results suggest that APPs and DFT-Ds could be improved at intermediate separations. A possible strategy is to replace the atomic separations r by an effective separation r−r0 , where r0 reflects the image-planes found for nanotubes and surfaces in Refs. 23,71. Keeping this (surface-physics) effect would increase the strength of the dispersion interactions at shorter separations. In summary, the results of this section suggests that an asymptotic atom-based pair potentials has a limited transferability over the full range of separations. Thus, for schemes using such a form, our results raises questions on their ability to generate accurate results under broad condition (having multiple characteristic separations), for instance involving phase transitions or processes that drive the system out off equilibrium, in protein unfolding, in phase transitions, or simply for systems which have competing interactions.1 We note that there is no guarantee that vdW-DF, in its current form, can provide an accurate account under such broad conditions. A vdW-DF limitation is here exemplified by the likely underestimation of C6 coefficients for the platonic molecules. Nevertheless, we argue that non-local functionals, like vdW-DF, hold the most promise for dealing with molecular configurations under broad conditions. This is because an electron-based approach provides a framework which naturally includes image-plane and multipole effects. It therefore holds the key to an account which describe the variation in dispersive response over the full range of separations. VI. CONCLUSIONS For the three cage-like organic molecular crystals, hexamine, cubane and dodecahedrane, vdW-DF gives lattice 10 parameters and cohesion energy that agree well with experiments, although for all three crystals the unit cell volumes are overestimated. A substantial sensitivity of lattice parameters and cohesion energy to the flavor of semi-local exchange signals that this overestimation stems mostly from the chosen form of exchange functional. We have also shown that, at short and intermediate separations, the full non-local correlation of vdW-DF is considerably larger than its corresponding atom-based asymptotic account. Notwithstanding that the asymptotic account of vdW-DF likely needs improvement, this enhancement indicates that the asymptotic 1/r6 form of atomic pair potentials, by construction, can not give a transferable account over a large range of separations. This paper underlines the usefulness of studying simple model sparse systems as molecular crystals to gain insight into methods intended for the study of sparse and supramolecular systems.1 Both DFT-D and vdWDF benefit from such testing, because they are designed to be parameter-free, and an accurate account of molec- 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 J.M. Lehn, Science 260, 1762 (1993). J. Kleis, B. I. Lundqvist, D. C. Langreth, E. Schröder, Phys. Rev. B 76, 100201 (2007). S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980). J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). Day, G. M., Chisholm, J., Shan, N., Motherwell, W. D. S. and Jones, W, Cryst. Growth Des.4, 1327 (2004). W. Kohn, Y. Meir, and D. E. Makarov, Phys. Rev. Lett. 80, 4153 (1998). J.F. Dobson, and J. Wang, Phys. Rev. Lett. 82, 2123 (1999). H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 91, 126402 (2003). M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92 (2004); Phys. Rev. Lett. 95, 109902(E) (2005). O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 130, 104105 (2009). O. A. Vydrov and T. Van Voorhis, Phys. Rev. Lett. 103, 063004 (2009). M. A. Basanta, Y. J. Dappe, J. Ortega, and F. Flores, Europhys. Lett. 70, 355 (2005). Y. J. Dappe, M. A. Basanta, F. Flores, and J. Ortega, Phys. Rev. B 74, 205434 (2006). X. Wu, M. C. Varga, S. Nayak, V. Lotrich, G. Scoles, J. Chem. Phys. 115, 8748 (2001). Q. Wu and W. Yang, J. Chem. Phys. 116, 515 (2002). S. Grimme, J. Comput. Chem. 25, 1463 (2004). P. Jurečka, J. Černy, P. Hobza, D. R. Salahub, J. Comput. Chem. 28, 555 (2007). V. Barone, M. Casarin, D. Forrer, M. Pavone, M. Sambi, A. Vittadini, J. Comp. Chem. 30, 934 (2008). ular crystals would suggest an accurate account of more complex assemblies of similar molecules. The presented molecular crystals provide particularly accessible cases of extended systems and can therefore be used in conjunction with future development of pair-potential methods and exchange-correlation functionals. VII. ACKNOWLEDGEMENTS We acknowledge E. Schröder for useful comments and discussions. We thank I. Sinno for help in designing the molecular-crystal schematics and Ø. Borck for access to his code for Bader analysis. The Swedish National Infrastructure for Computing (SNIC) is acknowledged for computer allocation and for KBs participation in the graduate school NGSSC. The work was supported by the Swedish research Council (Vetenskapsrådet VR) under 621-2008-4346. 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 M. A. Neumann and M. Perrin. J. Phys. Chem. B 109, 15531 (2005). D. C. Langreth, B. I. Lundqvist, S. D. Chakarova-Kack, V. R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, J. Kleis, L. Kong, S. Li, P. G. Moses, E. Murray, A. Puzder, H. Rydberg, E. Schrder, and T. Thonhauser, J. Phys.: Condens. Matter 21, 084203 (2009). V.R. Cooper, T. Thonhauser, A. Puzder, E. Schröder, B.I. Lundqvist, and D.C. Langreth, J. Am. Chem. Soc. 130, 1304 (2008). J. Kleis, E. Schröder, and P. Hyldgaard, Phys. Rev. B 77 205422 (2008). A. Kelkkanen, B. I. Lundqvist, and J. K. Nørskov, J. Chem. Phys. 131, 046102 (2009). L. Kong, G. Román-Pérez, J. M. Soler, and D. C. Langreth, Phys. Rev. Lett. 103, 096103 (2009). P. Sony, P. Puschnig, D. Nabok, and C. Ambrosch-Draxl, Phys. Rev. Lett. 99, 176401 (2007). N. Atodiresei, V. Caciuc, P. Lazic, S. Blügel, Phys. Rev. Lett. 102, 136809 (2009). K. Berland, T. L. Einstein, and P. Hyldgaard, Phys. Rev. B 80, 155431 (2009). K. Toyodaa, Y. Nakanoa, I. Hamadaa, K. Leec, S. Yanagisawaa, and Y. Morikawaa, Surf. Sci. 603, 2912 (2009). A. Gavezzotti, Modell. Simul. Mater. Sci. Eng. 10, R1 (2002). T. Yildirim, P. M. Gehring, D. A. Neumann, P. E. Eaton, and T. Emrick, Phys. Rev. Lett. 78, 4938 (1997). Plato, Timaeus, http://classics.mit.edu/Plato/timaeus.html. P. E. Eaton, and T. W. Cole, J. Am. Chem. Soc., 86, 3157 (1964). L. A. Paquette, R. J. Ternansky, D. W Balogh, G. J. Kentgen, J. Am. Chem. Soc. 105, 5446 (1983). M. N. Glukhovtsev, S. Laiter, A. Pross, J. Phys. Chem. 99, 6828 (1995) M. S. Ahmed, Acta Cryst. 5, 587 (1952). 11 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 R. Rao, T. Sakuntala, S.K. Deb, A.P. Roy, V. Vijayakumar, B.K. Godwal, and S.K. Sikka, Chem. Phys. Lett. 313, 749 (1999). G. H. Heilmeier, App. Opt. 3, 1281 (1964). R. G. Dickison and A. L. Raymond, J. Am. Chem. Soc. 45, 22 (1923). B. S. Hudson, D. G. Allis, S. F. Parker, A. J. RamirezCuesta, H. Herman, and H. Prinzbach, J. Phys. Chem. A 109, 3418 (2005). A. Almdari and F. Tabkhi, Chem. Eng. and Proc. 43, 803, (2003). Ullmann’s Encyclopedia of Industrial Chemistry 2008. S. Borman, Chem. Sci. News 72, 34 (1994). E. Ziambaras and E. Schröder, Phys. Rev. B 68, 064112 (2003). T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard, and D. C. Langreth, Phys. Rev. B 76, 125112 (2007). D. C. Langreth, M. Dion, H. Rydberg, E. Schröder, P. Hyldgaard, B. I. Lundqvist, Int. J. Quantum Chem. 101, 599 (2005). E. Ziambaras, J. Kleis, E. Schröder, and P. Hyldgaard, Phys. Rev. B 76, 155425 (2007). Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965). G. Román-Párez, and J. M. Soler Phys. Rev. Lett. 103, 096102 (2009). S. Li, V. R. Cooper, T. Thonhauser, A. Puzder, and D.C. Langreth, J. Phys. Chem. A 112, 9031 (2008). X. Wu, M. C. Vargas, S. Nayak, V. Lotrich, and G. Scoles, J. Chem. Phys. 115, 8748 (2001). A. Puzder, M.Dion, and D. C. Langreth, J. Chem. Phys. 124, 164105 (2006). A. Gulans, M. J. Puska, and R. M. Nieminen, Phys. Rev. B 79, 201105 (2009). Several new exchange flavors designed for use within the vdW-DF framework promises to improve the vdW-DF account of binding separations: V. R. Cooper, arXiv: 0910.1250v1; J. Klimeš, D. R. Bowler, and A. Michaelides J. Phys.: Condens. Matter 22, 022201 (2010). S. D. Chakarova-Käck, E. Schröder, B. I. Lundqvist, and D. C. Langreth, Phys. Rev. Lett. 96, 146107 (2006). Opensource code DACAPO, http://www.fysik.dtu.dk/CAMPOS/ L.N. Becka and D.W.J. Cruickshank, Proc. Roy. Soc. A 273 (1963). L. Hedberg, K. Hedberg, P.E. Eaton, N. Nodari and A.G. Robiette. J. Am. Chem. Soc. 113, 1514 (1991). K. Schmidt, M. Springborg, M. Bertau, F. Wahl, A. 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 Weiler, K. Scheumann, J. Worth, M. Keller M, H. Prinzbach, Tetrahedron 53, 10029 (1997). E. A. Arnautova, M. V. Zakharova, T. S. Pivina, E. A. Smolenskii, D. V. Sukhachev, and V. V. Shcherbukhin, Rus. Chem. Bul. 45, 1066 (1996). G. N. Ramachandran and W. A. Wooster, Acta Cryst. 4, 431 (1951). J. S. Chickos and W. E. Acree, J. Phys. Chem. Ref. Data 31, 537 (2002). J. C. Galluci, C. W. Doecke, and L. A. Paquette, J. Chem. Phys. 108, 1343 (1986). D. Nabok, P. Puschnig, and C. Ambrosch-Draxl, Phys. Rev. B 77, 245316 (2008). A. T. Hagler, E. Huler, and S. Lifson, J. Am. Chem. Soc. 96, 5319 (1974). W.F. Perger, J. Criswell, B. Civalleri, and R. Dovesi, Comp. Phys. Comm. 180, 1753 (2009). T. Yildirim, S. Ciraci, Ç. Kiliç, and A. Buldum, Phys. Rev. B 62, 7625 (2000). H. Margenau and N. R. Kestner, Theory of Interatomic Forces (Pergamon Press, Oxford, 1969), 2nd ed. G.Mahan, J. Chem. Phys. 43, 1569 (1965). E. Zaremba and W. Kohn, Phys. Rev. B 13, 2270 (1976). R. F. W. Bader, Atoms in Molecules: A Quantum Theory (Oxford University Press, Oxford, 1990). G. Henkelman, A. Arnaldsson, and H. Jónsson, Comput. Mater. Sci. 36, 254 (2006). Y. Andersson, D.C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 76, 102 (1996). We have explicitly verified that the precise assignment of the total molecular response onto the atoms does not alter the APP-vdWDF curves significantly. An assignment of the total asymptotic vdW-DF response according to the relative strength of Grimme yield a weak increase in the difference between APP-vdWDF and the non-local correlation. J. C. Slater and J. G. Kirkwood, Phys. Rev. 37, 682 (1931). Y. Andersson and H. Rydberg, Phys. Scr. 60, 211 (1999). M. Pavone, N. Rega and V. Barone, Chem. Phys. Lett. 452, 333 (2008). Claudio Morgado, Mark A. Vincent, Ian H. Hillier and Xiao Shan, Phys. Chem. Chem. Phys. 9, 448 (2007). S. Grimme, J. Antony, T. Schwabe and C. MckLichtenfeld, Org. Biomol. Chem. 5, 741, (2007). In our comparison, we neglect that these schemes are optimized for a particular basis set, arguing that our comparison does not require a fine-tuned version of DFT-D.