Binding of polycyclic aromatic hydrocarbons and graphene dimers in density functional theory

An early van der Waals density functional (vdW-DF) described layered systems (such as graphite and graphene dimers) using a layer-averaged electron density in the evaluation of nonlocal correlations. This early vdW-DF version was also adapted to approximate the binding of polycyclic aromatic hydrocarbons (PAHs) (Chakarova S D and Schröder E 2005 J. Chem. Phys. 122 054102). In parallel to that PAH study, a new vdW-DF version (Dion M, Rydberg H, Schröder E, Langreth D C and Lundqvist B I 2004 Phys. Rev. Lett. 92 246401) was developed that provides accounts of nonlocal correlations for systems of general geometry. We apply here the latter vdW-DF version to aromatic dimers of benzene, naphthalene, anthracene and pyrene, stacked in sandwich (AA) structure, and the slipped-parallel (AB) naphthalene dimer. We further compare the results of the two methods as well as other theoretical results obtained by quantum-chemistry methods. We also compare calculations for two interacting graphene sheets in the AA and the AB structures and provide the corresponding graphene-from-graphite exfoliation energies. Finally, we present an overview of the scaling of the molecular–dimer interaction with the number of carbon atoms and with the number of carbon rings.


Introduction
Polycyclic aromatic hydrocarbon (PAH) molecules [1] are formed by the fusion of multiple aromatic rings. Graphene [2,3], i.e. a single isolated sheet of graphite, is a two-dimensional (infinitely) extended system. The PAH molecules can be considered to be small flakes of graphene with edges passivated by hydrogen atoms. Their electronic properties make both PAH molecules and graphene relevant for technological applications. PAH molecules are of environmental interest and concern, as they are highly carcinogenic.
Except for the smallest molecules, PAH molecules are too large to be treated with accurate all-electron wavefunction-based methods. Density functional theory (DFT) presents an efficient way to treat larger systems with many atoms. However, the previous standard implementations of DFT lack the ability to describe the nonlocal dispersive, or van der Waals (vdW), interactions which are important for sparse systems such as dimers of the PAH molecules. This is remedied by the use of DFT with explicit and consistent inclusion of the nonlocal dispersive interactions, as in the functional vdW-DF [4]- [8].
In the study presented here, we report the results from calculations using the vdW functional of [5], and compare them with our earlier results from calculations using the version of the vdW functional in [4]. Like its predecessor, this newer vdW functional version is based on the adiabatic connection formula (ACF) [9,10] and uses a formulation that retains nonlocal correlations expressed in terms of a plasmon model for the electrodynamical interaction. The functional implementation was originally designed for use with finite molecules; in [11] the implementation was augmented to also function for extended systems (e.g. graphene sheets). To facilitate the comparison between the two vdW-DF versions we also include the results of calculations using the older version of the functional that was constructed for translationally invariant layered geometry systems [4,12], which was modified by two of us to be used in an earlier study of benzene [13] and several PAH dimers [14]. These two vdW functional versions will here be termed general geometry (gg) and layered geometry (lg) vdW-DF. Whereas 3 standard DFT implementations fail to predict binding for the PAH and graphene dimers, we obtain values for the binding distances and energies in reasonable or good agreement with quantum-chemistry studies [15]- [18], where available.
The molecular and macromolecular systems investigated in this study are the aromatic dimers of benzene, naphthalene, anthracene and pyrene (C 6 H 6 , C 10 H 8 , C 14 H 10 and C 16 H 10 , respectively) and (infinite) graphene sheets. Our investigation of the graphene systems contrasts the interaction in (co-planar) graphene pairs with that in graphite exfoliation, the process where a graphene sheet is pulled off a graphite surface. For the aromatic dimers and co-planar graphene pairs, we focus our investigation on the sandwich (AA) structure, but also investigate the slipped-parallel (AB) structure for naphthalene and graphene pairs. For graphite exfoliation, we assume a bulk-graphite organization (ABA) and investigate the adhesion of a graphene sheet when aligned in both configurations A and B (as specified by the underlying graphite-surface structure).
In recent years, benzene and PAH dimers have been studied in a number of calculations using quantum-chemistry methods, DFT methods and methods derived from DFT calculations but using additional empirical input. Several methods and sets of results are discussed or mentioned in the comprehensive review in the introduction of [19] and in [15]- [18], [20,21]. We compare the results of vdW-DF calculations to accurate quantum-chemistry calculations for some PAH dimers (mainly in the AA stacking). The vdW-DF method has previously been used to successfully treat benzene dimers in the AA, AB and T-shape structures [22], systems with monosubstituted benzene dimers [23] and PAH molecules on graphite [11].
The rest of the paper is organized as follows. First, the graphene systems and PAH molecules are briefly described, followed by a short review of vdW-DF and how we apply it here in the gg and lg versions. The computational method is presented in section 4. We then present the results in section 5. Finally, in section 6 we discuss the results and compare them with a selection of those from wave function calculations, and with our previous results obtained with the lg version of the vdW-DF.

Materials: graphite, graphene, benzene and polycyclic aromatic hydrocarbon (PAH) molecules
Graphite consists of sheets of carbon atoms on a hexagonal lattice. Graphene is an isolated sheet of graphite, produced, for example, by exfoliation from a graphite surface [2]. The sheets can be stacked so that only every second carbon atom is covered with atoms of the sheet on top of it (AB stacking). This is known to be the most stable structure for graphite. The AA stacking, in which the carbon atoms of the different sheets are aligned directly on top of each other, is energetically less stable. In both cases the layers of graphite (or graphene) sheets interact via the vdW interactions, while within each sheet the carbon atoms are covalently bound.
PAHs are a family of molecules with, as the name indicates, multiple aromatic rings. They may be seen as flakes of graphene sheets passivated by hydrogen atoms at the perimeter, except that the C-C bond lengths in the PAHs differ slightly (up to approximately four per cent for the PAHs considered here) from the ones found in graphite. As expected, the C-C bond lengths in large PAH molecules resemble those in graphite more closely than those of small PAHs [24].
The energetically most favored stacking of most of the PAH dimers is slipped-parallel, corresponding roughly to the AB stacking in graphite. However, there are several different AB-stacking types, depending on the direction in which the molecules are slipped with respect 4 to each other. For these systems an exact AB stacking is not possible since the C-C bond lengths are not uniform [15]. Fortunately, small lateral positional changes of the relative dimer structure cause very small energy differences [15]. It is therefore possible to present a simple overview of the PAH dimer interaction by focusing on the cases of AA stacking and by supplementing those results with a few alternative dimer configurations.
PAH molecules are ubiquitous in our environment and in the universe. A common way to inadvertently build PAHs is to burn graphite or other carbon materials at low temperatures, which leads to PAHs in the smoke. This happens in fireplaces, in poorly regulated industrial burning and in cigarette smoke. Furthermore, PAH clusters are believed to self-assemble in interstellar space [25,26] as intermediates between free gas-phase PAHs and amorphous carbon particles, and PAH clusters in interstellar settings have been spectroscopically explored [27].
Since graphite itself does not pose any immediate health danger it is perhaps surprising that PAH molecules are, in general, toxic and carcinogenic. One suggested reason for the latter property is that the aromatic rings of the PAHs interact with the aromatic base-pairs of DNA [28], thereby disturbing its function, although the exact mechanism of action is still under investigation [29]. and [33], it was shown how to extend the original implementation to molecules with infinite extensions in one or more directions, e.g. polymers and sheets of graphene.
The gg version of vdW-DF has been used extensively in the past few years, treating the polyethylene crystal [32], PAH and phenol molecules on graphite and α-Al 2 O 3 (0001) [11,34], K intercalation in graphite [33], stacking interactions of DNA [35], nanotube crystals [36] and a number of other systems [8].
In the lg version of vdW-DF [4], is assumed to be anisotropic in the sense that its components in the directions perpendicular and parallel to the layers differ. A materialsdependent constant q ⊥ is introduced in the model dielectric function . It is determined by requiring that the effective zero-frequency susceptibility (for sheets) or the static image plane position (for surfaces) in the model must correspond to the static response computed by standard (GGA-)DFT [4,6,12]. The lg version requires that the electronic charge density n(r) is averaged laterally, leading to a one-dimensional differential equation for the electric field that can be solved exactly. Thus the lg version assumes nearly translational invariance in the lateral direction. In a previous paper [14], two of us showed how to modify it to also treat finite-but still planar-molecules, like the PAHs, albeit in an approximative way.
In both the lg and gg versions, plasmon-pole-like approximations are used to model the dielectric function, with input only from traditional DFT calculations (implicitly in gg, and explicitly through q ⊥ in lg). In this paper, we present new results using the gg version and we compare them with our previously published results using the lg version [14].

Computational method
The atomic and electronic structure of the isolated PAH molecules and the graphene sheet is found by performing standard GGA DFT calculations. The bonds within the molecules are covalent and well described by GGA. All atomic positions are allowed to relax self-consistently within GGA according to the Hellmann-Feynman forces. Since we use a plane-wave basis set with periodic boundary conditions, the unit cell chosen must be sufficiently large for the periodic images of the molecule not to interact within GGA. These and the further GGA calculations described below are performed self-consistently within the revPBE flavor of GGA [37], using the plane-wave pseudopotential code Dacapo 5 .
For the graphene-sheet calculations we use a hexagonal unit cell of size (2.46, 2.46 and 26 Å), and for the PAH monomers and dimers, the unit cell size is (17.11, 17.11 and 26 Å), except for pyrene, where the size is (21.39, 21.39 and 20 Å). In all cases, we use ultrasoft pseudopotentials, a plane-wave energy cutoff of 450 eV and a fast-fourier-transform (FFT) grid for the charge density with approximately 0.11-0.13 Å between the grid points. The bond lengths within the PAH molecules were determined in [14]. The lattice constant for the graphene sheets was determined in [11]; for the graphene-on-graphite systems that we use to investigate graphite exfoliation, the computational details were slightly different. 6 5 Open-source plane-wave DFT computer code Dacapo, http://wiki.fysik.dtu.dk/dacapo. 6 The present study of and comparison with graphite exfoliation also represents an extension of the minimal analysis (for configuration ABA-B) that permitted a vdW-DF characterization of potassium absorption and intercalation in graphite, see [33]. That study naturally begins with the bulk graphite bulk structure, for which our vdW-DF characterization yields a 2.476 Å in-plane lattice constant. To be consistent with [33], we choose to characterize graphite exfoliation in both configurations ABA-B and ABA-A, assuming the 2.476 Å in-plane (graphite) lattice constant. 6 We use the vdW-DF as a post-GGA method: the GGA calculations are performed selfconsistently and provide us with a GGA total energy, GGA correlation energy and an electron density n(r) that will be used to evaluate E nl c . A method to evaluate the vdW total energy E vdW-DF self-consistently exists [7], but it has been shown [7] that the effect of using self-consistency is small on E vdW-DF (whereas other quantities, related to the electron distribution, are more clearly affected).
The first step in obtaining the total energy E vdW−DF (given by equation (2)) is to evaluate the GGA-based total energy E GGA and n(r) self-consistently (within GGA) and subtract the correlation energy E GGA c . From n(r) we further evaluate (non-self-consistently) the LDA correlation energy E LDA c and the nonlocal E nl c (as further described below) in order to compute the correlation 7 according to equation (1). The atoms of the two molecules or graphene sheets are kept fixed in their lateral monomer positions while the molecule separation is varied stepwise.
The gg version of vdW-DF does not make use of periodic boundary conditions. This makes it ideal for application to finite molecules [5,6], such as the PAH dimers. The E nl c is determined by equation (3), as described in [5].
The lack of explicit periodic boundary conditions in the gg code for E nl c , in contrast to the underlying GGA calculations, is a minor complication for extended and (in principle) infinite systems like the graphene dimer. For such systems, the original implementation of gg from [5] must therefore be modified. For two parallel interacting graphene sheets we calculate E nl c as the sum of the interaction within the unit cell of the underlying GGA calculations and interactions of spatial points (i.e. FFT grid points) in the unit cell with all spatial points of the (laterally) neighboring periodically repeated images of the unit cell. This is further described in [33] and a similar method is used in [11]. The quantity of grid points needed for convergence of E nl c depends on the system; for graphene we found it sufficient to include the interaction of points within a radius of 10 Å.
As the zero point of the energy scale we choose the dilute gas phase of the molecules. For computational reasons, discussed in [13] for benzene dimers but also more generally valid [33,36], care needs to be taken when choosing how to carry out the reference structure calculations. For the part of the energy not including E nl c , that is, the energy contribution , a dimer calculation where the molecules are separated by a large distance is used. The reason for this is that the technical difficulties associated with adding monomer energies (at the level of GGA) yield a small but noticeable offset in the local part of the energy [13]. For the nonlocal part of the energy (E nl c ), no such problems arise and the reference energy is taken as the energy of the two subsystems of the dimer, that is, two systems each with a single PAH molecule positioned at the same position in the unit cell as the corresponding molecule of the dimer.
The vdW interactions are calculated from charge densities that include sparse regions of space. Because n(r) is described on a discrete mesh in space (namely, on the FFT grid) we enhance the accuracy of the E nl c calculation by a procedure that consistently positions the grid points at the same location relative to the molecule, for the dimer and for the reference systems (isolated molecules). 7 Our calculation of the vdW-DF total energy for the graphite-exfoliation systems proceeds in a slightly different but fully equivalent procedure. To remain consistent with the previous analysis summarized in the study of potassium intercalation and absorption [33], we investigate the exfoliation by performing self-consistent calculations in the PBE [30] flavor of GGA. From the PBE electron densities, we then determine the vdW-DF energy by estimating both the LDA correlation energy and the revPBE-GGA [37] exchange energy. Graphene dimer calculations in the lg version are straightforward and are described in [4,6] and [13]. For the finite PAH molecules, application of lg requires a small modification to account for the finite size of the molecules, as described in [13] and [14].

Results
The main aim of this paper is to report the results on the calculated vdW binding in small PAH dimers and in graphene systems using the gg version of vdW-DF and to compare these results with our previous results from the lg version. The lg results for the PAH dimers were obtained mostly with AA stacking; therefore our present gg results for PAH dimers are mainly for AA-stacked molecules. The lg results for the PAH dimers included a modification of the original lg version, in order to treat the finite (lateral) size of the molecules, so in order to compare the two vdW-DF versions more directly we also include a comparison of graphene dimer results in AA and AB stacking. Figure 1 shows the interaction energies of two graphene sheets, arranged in AA and AB stacking, obtained with the two versions of vdW-DF. Both show a clear binding minimum. This is in contrast to GGA, which for such systems gives repulsive or very weak binding at unrealistically large binding distances. For a detailed discussion on the differences between the GGA and the vdW functional, see, for example, [5,6,36]. Table 1 presents a summary of the binding distances (d b ) and binding energies per surface C atom (E b /N ) for graphene pairs, graphite exfoliation and bulk graphite. With gg we obtain for the graphene dimer a binding energy of 42 meV atom −1 at the distance of 3.75 Å for the AA Table 1. Comparison of carbon-sheet interaction in graphene pairs (AB and AA), in graphite exfoliation (ABA-B and ABA-A) and in graphite bulk (AB n→∞ ). The table lists binding energies per surface C atom, E b /N , and binding distances, d b , as calculated in the lg version (when available) and gg version of vdW-DF.  stacking, and the slightly larger 47 meV atom −1 at the distance of 3.61 Å for the AB stacking. This is in agreement with the experimental fact that graphite organizes in AB stacking. From the lg calculations we obtained a somewhat weaker binding for both the graphene dimer and the bulk graphite. This is further discussed in section 6. Figure 2 shows the graphite-sheet interaction for graphite exfoliation, investigated with the gg version of vdW-DF and comparing two possible alignments of the graphene overlayer on the graphite surface, ABA-A and ABA-B. For both alignments, we find a clear minimum. Table 1 presents the values of d b and E b /N obtained from a fit to the calculated vdW-DF energy  The four molecules are illustrated with dark spheres for the C atoms and with small white spheres for the H atoms. The insets show the fits to fourth-order polynomial curves of data points that are close to the minimum.

Graphene-sheet interaction
variation. For the ABA-A configuration, the graphene sheet is pushed further out and couples less to the underlying graphite surface. Contrasting the binding with that calculated for a pair of graphene sheets, we find that the subsurface graphite sheets on their own cause only a very small strengthening of the adhesion.

Benzene and PAH dimers
The interaction energy curves calculated with both lg and gg for the benzene, naphthalene, anthracene and pyrene dimers in AA stacking are shown in figure 3. Both the lg and gg versions give realistic energy curves for all the PAH dimers, i.e. with significant binding. For the benzene dimer the two versions perform similarly, while for the larger PAHs the quantitative difference in binding energy successively increases with increasing molecule size.
In table 2, we report the binding distances (d b ) and the binding energies (E b ) for the benzene and PAH dimers. The minima of the curves are found by fitting a fourth-order Table 2. Dimer binding distances and energies obtained with the lg and gg versions of vdW-DF. polynomial to the data points close to the minimum (the points shown in the insets of figure 3) to compensate for the small noise present in the calculations, originating mainly from the exchange part of the GGA energy (see, for example, the revPBE interaction energy curves in [14]).
Our methods are not limited to treating PAH dimers in the AA stacking only, although this has been the focus here. Figure 4 shows the performance for the naphthalene dimer in one possible AB stacking. The figure also illustrates a top view of the dimer structure with the lateral displacements R 1 = 1.2 Å and R 2 = 0.7 Å. In this case, the top molecule is placed in such a way that the position of its carbon atoms resembles as much as possible the positions of the atoms in graphene in the AB stacking. The binding energy is 0.29 eV at the optimal vertical separation of 3.8 Å (see table 2). Hence, as expected, a naphthalene dimer in AB stacking is more energetically favorable than a naphthalene dimer in AA stacking. The stronger bond between the naphthalene molecules results in a smaller equilibrium distance for the AB stacking compared with that for the AA stacking.
In summary, we find that including the vdW interaction in the DFT calculations of PAH and graphene dimers is necessary for a thorough description of their binding. As expected from experiments we find binding of the dimers, and a stronger binding for the AB-stacked molecules (where calculated) than for the AA stacked. The two different versions of the vdW-DF, the layered geometry [4,6,12] and the general geometry [5] versions, give similar results for the benzene dimer in the AA stacking. The general tendency is that the gg version gives somewhat larger binding energies than lg, with the differences increasing with the size of the PAH. This is especially clear from figure 3 for the PAH dimers in AA stacking. The quantitative differences in the results from the two versions of the functional are further discussed in section 6.

Discussion
First we discuss the graphene-dimer and exfoliation results in comparison to experiments and previous calculations and then we discuss the benzene and PAH dimer results compared to quantum-chemical calculations. We also discuss the scaling of the dimer binding energy with the size of the molecule.

Character of van der Waals (vdW) binding
Contrasting the graphene-dimer and the graphene-from-graphite exfoliation energies provides an illustration of the limited effective range that the dispersive or vdW interactions have. The relatively short-range character of the vdW binding is also documented and discussed in a range of other vdW-DF studies, for example [8,11,36]. The findings suggest the additional possibility of using vdW-DF calculations to model (nonperiodic) supramolecular systems. For dense matter, the DFT community often uses accurate calculations for system segments to also provide insight into larger-scale system behavior. This will be possible if the segment includes all relevant interaction effects. The fact that the main contribution to the vdW binding arises in a limited region makes it possible to adapt the strategy for the investigation of at least some supramolecular problems.
For the graphene dimer, we find in section 5 that the gg approach gives stronger binding than the lg approach. In figure 1, two experimentally obtained estimates of the interaction energy in (AB-stacked) graphite are also marked. The first experimental value, 35 +15 −10 meV atom −1 , is deduced from the interactions in collapsed carbon nanotubes [38], while the second, 52 ± 5 meV atom −1 , results from an analysis of thermal desorption experiments of PAHs on graphite [39]. 8 None of these experiments give information about the binding separation; in figure 1, we therefore set the experimental data points at the bulk sheet separation of graphite (3.34 Å) [40]. We expect the true separation of AB graphene to be similar to or slightly larger than this bulk graphite interlayer distance.
Thus, from figure 1 and table 1, we find that for the graphene dimer, gg gives results that are numerically more close to experiments than the results of lg. This may seem somewhat unexpected as the lg method is designed for (translationally invariant and extended) layered systems, while the gg method is designed more generally with also other geometries in mind. The lateral averaging of the electronic charge density introduced in lg cannot be held responsible, as gg applied with the same averaging gives an even stronger binding [41]. The difference in performance for the graphene sheet systems must rather be attributed to the different sets of approximations used in the two versions [5,12].
Without going into the details, we stress the difference in nature and applicability between planar-geometry and general-geometry versions of vdW-DF. Whether the graphene dimer system is to be considered planar or point-like depends on the scale. On a large scale, the graphene sheets look very planar, but on a small (atomic) scale each atom in one of the graphene sheets only experiences the presence of a handful of other atoms in the other graphene sheet; a very general geometry. The integration of the six-dimensional integral for the nonlocal correlation energy (defined in equation (3)) of the graphene sheets hints at the small-scale picture being the most appropriate one. In practice, we find that a 10 Å radius circle of graphene in the extended graphene unit cell is sufficient to already include 99% of the interaction strength, and to obtain a slightly smaller part of the interaction strength we need to include a circle of a radius comprising just a couple of atomic distances. Reasoning along these lines leads us to also favor the gg version for graphene dimers.
Since lg gives weaker binding than both gg and the true binding for the graphene sheets, it is to be expected that this will also be the case for the PAH dimers. This is indeed what we find in section 5; namely, the larger the PAH system, the larger the binding energy difference between lg and gg. A class of systems for which experimental data were obtained recently is that of PAHs adsorbed on graphite [39]. Use of the gg version of the vdW-DF on benzene and naphthalene adsorbed on graphite yields results in very promising agreement with these data [11].

Comparison with quantum-chemistry calculations: benzene and PAH
For the benzene and PAH dimer interactions we are not aware of any existing experimental results that offer a direct comparison to our vdW-DF calculations. The fact that the AA structures, on which we mainly focus here, are energetically less favorable than the AB structures makes such a comparison difficult. We can, however, compare our results with the ones available from quantum-chemistry calculations.
A number of quantum-chemistry calculations are available for the benzene dimer, as it is the smallest aromatic system and thus the most tractable one. A coupled cluster method with single, double and noniterative triple excitations [CCSD(T)] study [16] found the binding energy in the AA stacking to be 1.63 kcal mol −1 (0.071 eV dimer −1 ) and older calculations, with various methods, give binding energies in the range 0.06-0.16 eV for separations around 3.7-4.1 Å (references given within [13], [17] and [22]). We find these separations to be similar to the 4.0 Å obtained with our gg version of vdW-DF. Also the binding energy agrees reasonably well with our 0.12 eV. The same gg method applied by our Rutgers collaborators using a different plane-wave code and with norm-conserving pseudopotentials gave a binding of 0.10 eV dimer −1 at 4.1 Å [22]. The interaction energy curve is shallow around the binding distance, as shown in figure 3; thus there is some uncertainty and sensitivity to details in the calculations associated with the resulting binding distance. The lg and gg calculations are expected to give large equilibrium separations due to our use of the revPBE flavor of GGA for the exchange contribution [5]. Using exact Hartree-Fock (HF) exchange instead of revPBE exchange yields the values 3.8 Å and 0.12 eV for the benzene dimer in the AA stacking [22]. As this effect is quite general [22], we expect to find in this work slightly larger binding distances for all PAH dimers. Also, the AB-stacked and the T-shaped benzene dimers were considered in some of the above studies, but those structures are not investigated here.
In general, accurate quantum-chemistry calculations for the PAH dimers are much more scarce than for the benzene dimers, as the size of these systems pushes the methods toward their limits. The few calculations that exist are mostly on different stackings than the AA stacking that is the focus here.
For the naphthalene dimer in the AB structure, CCSD(T)-based calculations of a similar parallel-displaced structure with slightly different in-plane displacements, R 1 = 1.4 Å and R 2 = 1.0 Å [18], give an optimal vertical separation of 3.5 Å. This is slightly smaller than the separation of 3.8 Å found for gg and 3.7 Å found for lg. The binding energy of 0.247 eV dimer −1 [18] agrees better with the lg value (0.25 eV dimer −1 ) but is also close to the gg value (0.29 eV dimer −1 ). Other calculations referred to in [18] give results in the range 0.16-0.46 eV dimer −1 .
A recent MP2 study deals with the naphthalene AB structure similar to the structure used here [15]. An equilibrium separation of 3.35 Å and an AB binding energy of 0.256 eV is found for the naphthalene dimer.
Comparing the vdW-DF results for naphthalene dimers in AA and AB stacking (figures 3 and 4 as well as table 2) shows that the influence on E b of the in-plane position of the molecules is several orders of magnitude smaller than that of the interplane distance. Similar results have been shown in [15]. Figure 5 presents the AA stacking binding energy per C atom (E b /N ) and the corresponding equilibrium separation d b as a function of the size of the PAH molecule, measured as the number of C atoms in the molecule. Also shown are the corresponding values for the AA-stacked graphene dimer. The binding becomes stronger (the value of E b /N increases) as the molecules grow in size and the d b value decreases. This is in agreement with even the simplest effective medium theory picture. However, in the PAH systems that we consider here, the binding is purely due to vdW interactions and has no covalent contribution; therefore the intuitive picture is not a priori valid. It is clear that E b /N is not constant, and neither E b /N nor d b scales linearly (which would also be counter-intuitive, as the graphene dimer has finite values of E b /N and d b and can be seen as the asymptote of growing PAH molecules). This illustrates that the intriguing vdW interaction in these systems cannot be accounted for in a simple manner. The vdW interaction is not only present between C atoms that are on top of each other (for the AA stacking) but also between atoms in the two molecules that are not directly above each other. Figure 6 shows the vdW-DF results for the PAH-dimer binding energy E b versus the number of cyclic rings and shows an approximately linear scaling. However, the group of linear PAHs, such as benzene, naphthalene and anthracene, stands out with an almost exact linear scaling as a function of the number of carbon rings. The linear PAHs, also called acenes, are linearly fused benzene rings. Pyrene is not an acene, and the E b of pyrene falls somewhat off the linear curve in figure 6. The linear scaling of the AA-stacked acenes predicted by our vdW-DF calculations is in agreement with the results of symmetry-adapted perturbation theory calculations in [43] concerning other stackings than AA. It is also consistent with the experimental observation of a linear scaling in adsorption energies on Cu(111) for acenes [43], and in the adsorption energy of a number of cyclic molecules on the basal planes of MoS 2 [44]. The linear scaling suggests the possibility of seeking approximation schemes for vdW binding effects in molecular clusters (for example [45]).

Conclusion
We calculated the binding energies for a group of PAH dimers and dimers of graphene sheets, as well as the exfoliation energy of a graphene layer off the graphite surface, all within firstprinciples DFT. We applied a recent functional that included the nonlocal dispersion (van der Waals) interactions, which are essential for binding in graphene and PAH dimers. Our results agree reasonably well with the few available experimental results and with highly accurate quantum-chemistry calculations. A linear scaling is found for the binding energy of the dimers of acenes (linear PAHs). However, the dimer of the one PAH in this study that is not linear deviates from this trend, revealing the intriguing complexity of the vdW interaction even in such simple systems as the PAH dimers.