Dispersion interactions between semiconducting wires

Alston J. Misquitta,1 James Spencer,2, 3 Anthony J. Stone,2 and Ali Alavi2

arXiv:1005.1332v2 [cond-mat.mes-hall] 20 Jul 2010

Cavendish Laboratory, 19, J J Thomson Avenue Cambridge CB3 0HE, U.K. University Chemical Laboratory, Lenseld Road, Cambridge CB2 1EW, U.K. 3 Department of Physics and Thomas Young Centre, Imperial College London, Exhibition Road, London SW7 2AZ, U.K (Dated: July 21, 2010)

The dispersion energy between extended molecular chains (or equivalently innite wires) with non-zero band gaps is generally assumed to be expressible as a pair-wise sum of atom-atom terms which decay as R6 . Using a model system of two parallel wires with a variable band gap, we show that this is not the case. The dispersion interaction scales as z5 for large interwire separations z, as expected for an insulator, but as the band gap decreases the interaction is greatly enhanced; while at shorter (but non-overlapping) separations it approaches a power-law scaling given by z2 , i.e. the dispersion interaction expected between metallic wires. We demonstrate that these eects can be understood from the increasing length scale of the plasmon modes (charge uctuations), and their increasing contribution to the molecular dipole polarizability and the dispersion interaction, as the band gaps are reduced. This result calls into question methods which invoke locality assumptions in deriving dispersion interactions between extended small-gap systems.
PACS numbers: 34.20.Gj 73.22.-f 31.15.ap

INTRODUCTION

the atomatom R6 interaction.

Conventionally, the dispersion interaction between two systems is formulated in terms of correlated uctuations described by local (frequency-dependent) polarizabilities ab 6 Rab inwhich gives rise to the familiar atomatom C6 teraction at leading order [1]. This form of the interaction has been used with a good measure of success in studies of many systems, from gases to solids, including complexes of biological molecules and organic molecules. These are typically insulators; that is, they have large HOMOLUMO gaps. It has been known for a very long time that the additive atomatom form of the dispersion does not hold very well for metallic systems. The classic case of an atom interacting with a thin metallic surface shows deviations from the additive law: the correct treatment of the metal results in a R3 power law [2], but summing over atomatom terms leads to a R4 interaction. More recently, strong deviations from additivity have been demonstrated in the interactions of extended systems of metals or zero band-gap materials with at least one nano-scale dimension [37]. The semi-classical picture is useful for understanding the source of the dierences in the metallic and insulating cases. In an insulator, electronic perturbations decay exponentially with distance. We can therefore treat electron correlations as being local. At lowest order, these local uctuations give rise to instantaneous local dipoles, and the correlation between these local dipoles gives rise to the atomatom R6 interaction. However, in a zero bandgap material, electronic uctuations are long-ranged, particularly if one or two of the dimensions are nano-scale [3]. It is these long-ranged uctuations that give rise to the non-additivity of the dispersion and deviations from

While the insulating and metallic cases are now well understood, little is known of the intermediate, semiconducting case. The nature of the dispersion interaction between extended molecules with nite but small HOMOLUMO gaps less than about 0.2 a.u. (5 eV) remains an open question. A large number of important nano-molecules fall into this category, such as carbon nano-tubes and the lander-type molecules that are used as organic conductors. The electronic structure of materials made of these molecules depends strongly on structures they assume in the bulk, often via self-assembly. Consequently it is very important to understand exactly how these molecules interact. The dispersion interaction is the dominant source of attraction between these conjugated systems. For want of a clearer understanding, many studies assume the usual insulating case for the atomatom R6 form of this interaction. As we shall show in this article, this is both qualitatively and quantitatively incorrect for semiconducting molecules.

A word of clarication: we use the term non-additive here to describe deviations from the additive atomatom picture of the dispersion interaction between two molecules. It is also commonly used to refer to the deviation from pair additivity seen in the interactions of three or more distinct molecules, but we are not concerned with such eects here, except for a brief comment in the Discussion.

2
INTERACTING WIRES USING HCKEL (TIGHT-BINDING) THEORY

and the set of wavevectors by k= j ; nd n n n j = + 1, + 2, , . 2 2 2 (5)

The dispersion energy appears at second-order in intermolecular perturbation theory and is formally expressed in terms of the exact eigenstates and eigenenergies of the non-interacting systems [1]. In a mean-eld theory the dispersion energy between two subsystems (A and B), can be expressed as a sum over the single-electron wavefunctions localised to each subsystem:
(2) Edisp = 1 | i j|r12 |ab |2 , + j a b iA, j B aA,b B i

(1)

where i, j (a, b) are occupied (virtual) single-particle wavefunctions in either subsystem A or B, and i is the eigenvalue of the i-th wavefunction. Assuming two parallel wires, aligned parallel to the x axis and separated by a distance z, the integral which appears above has the form:
1 |ab = i j|r12 i ( x1 ) j ( x2 )a ( x1 )b ( x2 )

= corresponds to a uniform wire with energy eigenvalues given by (k) = 2| cos(kd)|, which in the limit of large n has a vanishing band-gap at half-lling (i.e. is a metal). The opposite limit ( = 0) corresponds to a chain of isolated dimers with energy eigenvalues (k) = , independent of wavevector, and corresponds to the perfect insulator. By varying the ratio / between 0 and 1, we can very conveniently probe the dispersion interaction in the intermediate (semiconducting) regime, with the bandgap given by Eg = 2( ). The required integrals can be evaluated [8] as:
1 |ab = i j|r12

1 nd

K0 (Gz)Yia (G)Y jb (ki + k j ka kb G)

+ zz | | x12 x

dx1 dx2 (2)

and is the Coulomb interaction between the charge density i a due to excitation i a in subsystem A with in B. j b In periodic systems, the wavefunction can be expressed in Bloch form, i.e. j ( x) = eik j x u j ( x), and so i(ka ki ) x the co-density ui ( x)ua ( x) can have i ( x)a ( x) = e a long-wavelength modulation, whose electrostatic eld will not be well described by a multipole expansion at separations comparable to this length scale. Furthermore, in small-gap systems, these excitations have plasmon (2) character and contribute signicantly to Edisp , as we discuss below. Hckel (tight-binding) theory gives us a convenient formalism for evaluating the above expression for interactions between two one-dimensional wires. We considered a two-band model Hamiltonian of the form
n

(6) where {G} are reciprocal lattice vectors of the primitive unit cell and Yia (G) is the (analytic) Fourier transform of u i ua obtained by assuming highly-localised basis functions on each atom. K0 is the zero-th order modied Bessel function of the second kind and is the Fourier transform of the potential generated by a onedimensional lattice at a eld point at z from the lattice[9 11]. K0 decays exponentially with increasingly large arguments and so typically only the 10 smallest reciprocal lattice vectors need to be included in the summation in order to obtain converged results. The calculations presented here use 8401 MonkhorstPack k-points to sample the Brillouin zone, which corresponds in real-space to wires consisting of 16802 sites per crystal cell. Such ne k-point sampling was required to obtained converged results with respect to system size. As shown in g. 1, the dispersion interaction between two uniform wires ( / = 1) varies as z2 across the whole range of z, in good agreement with previously obtained results via RPA[3] and QMC[7] for metallic wires. The interaction between chains of isolated dimers ( = 0) varies as z5 , as would be expected from the conventional picture of dispersion interactions. In the semiconducting regime, there is no unique exponent which characterises the dispersion interaction over all z. However, at small separations z, it tends to the metallic behaviour (z2 ), whereas large z it tends to insulating behaviour z5 .

H=
i

(a 2i a2i1 + a 2i+1 a2i + h.c.),

(3)

where , are alternating bond-strengths between adjacent sites, as would be encountered in a chain of (H2 )n or a -conjugated polyene. We computed the interactions between two such parallel wires, each consisting of 2n identical atoms, equally spaced at intervals of d, such that the unit cell is of length 2d and contains two atoms. Assuming periodic boundary conditions over a crystal cell of length 2dn, then there are two bands per wavevector. The single-particle wavefunctions and energies of sych a system can be found analytically. The band structure is given by (k) = eikd + eikd (4)

INTERACTING H2 CHAINS USING SAPT(DFT)

The principal drawback of the Hckel model is the lack of electron correlation, and hence screening, within

3 ratio of the alternate bond lengths. The HOMO-LUMO gap of a hydrogen chain can be modied by simply introducing a dierence in alternate bond lengths. The distortion parameter gives us a convenient way to control the electronic structure of the chain and allows us to study the interactions between insulating, semiconducting and (near) metallic chains in one framework.

FIG. 1: The second-order dispersion energy of parallel wires modelled using Hckel theory obtained from eq. (1). The bonding interaction is of alternating strength and is controlled by the ratio /. The lines are numerical power-law ts to the small-z and large-z data points. In the metallic ( = ) and perfect insulating ( = 0) cases, a single power-law ts the whole range of z, with the expected exponents 2 and 5 respectively. The intermediate cases show a cross-over between these two limiting regimes as z is varied.
x

FIG. 2: Distorted (H2 )n -chains. The distortion parameter is dened as = y/ x. In all our chains we have chosen x = 1.4487 a.u.

The frequency-dependent density response functions that appear in eq. (7) were evaluated with coupled Kohn Sham perturbation theory, using the PBE0 functional with a hybrid kernel consisting of 75% adiabatic LDA and 25% coupled Hartree-Fock. The accuracy of this approach for the dispersion energy of small molecules surpasses that of Mller-Plesset perturbation theory and rivals that of coupled-cluster methods [13, 14, 16, 17]. Furthermore, there is no qualitative change in our results when the fraction of HF exchange is increased, so the results presented here are unlikely to be artifacts of the shortcomings of coupled KohnSham perturbation theory [18, 19]. All calculations used the cc-pVDZ basis which will result in a signicant underestimation of the strength of the contribution to the dispersion energy that arises from transverse polarization, but should better describe the contribution arising from the longitudinal polarization. Since it is the longitudinal polarization that is most important in these systems, we expect our results to be qualitatively correct.

each subsystem. For more realistic calculations we can use the symmetry-adapted perturbation theory based on density-functional theory [1214] (SAPT(DFT)), where the dispersion energy is evaluated not via a sum-overstates but through a coupled KohnSham formulation based on the density response functions of the interacting molecules[15]:
(2) Edisp = 1 dw dr1 dr 1 d r2 d r2 2 0 A ( r1 , r 1 ; iw) B (r2 , r2 ; iw) |r1 r2 ||r 1 r2 |

(7)

where A/B (r, r ; w) are the frequency-dependent density response functions of the molecules which describe the propagation of a frequency-dependent perturbation in a molecule, from one point to another, within linearresponse theory. The second-order dispersion energy is exact in this formulation if exchange eects are neglected. We have studied the interactions between two parallel nite (H2 )n (g. 2) chains with n = 16 and 32 and distortion parameters = 2.0, 1.5, 1.25 and 1.0, where is the

Dispersion energies are displayed in g. 3 and the effective power laws for the physically important separations between 6 and 20 a.u. are shown in table I together with HOMOLUMO gaps. For insulating chains with an additive dispersion interaction we would expect two regimes determined by chain length L: for z L we should expect the innite chain result, i.e., the eective dispersion interaction should decay as z5 , and for z L we should recover the usual z6 power law. This is what we see with the chains with the largest distortion parameter, = 2.0, which exhibit large HOMOLUMO gaps. As the distortion parameter approaches unity and the HOMOLUMO gap consequently decreases, the deviation from the additive insulating case becomes increasingly apparent, and for the longest of the chains considered here the dispersion interaction is signicantly enhanced. For example at 40 a.u. (roughly half the chain length) an = 1 chain has a dispersion interaction two orders of magnitude larger than an = 2 chain. Additionally, we see increasing nite-size eects: the power laws for the n = 16 and 32 chains, which were very similar for = 2.0, are considerably dierent for = 1.0.

4 (7) we obtain the multipole form for the dispersion energy:

1 ab a b T T 2 tu t u

bb aa tt (iw)uu (iw)dw.

(8)

ab Here T tu is the interaction function between multipole t on site a in subsystem A and multipole u on site b in B, while aa tt are the frequency-dependent non-local polarizabilities for sites a and a , and may be expressed in terms of the frequency-dependent density susceptibility as[20]

FIG. 3: The second-order dispersion energy of pairs of parallel (H2 )n chains. Energies have been normalized by the number of H2 units n. The solid lines are dispersion energies for chains with n = 32 and the dashed lines for chains with n = 16 (only =1.0 and 2.0).

aa tt () =

a a 3 3 Q t (r)(r, r |) Qt (r )d rd r .

(9)

(2) TABLE I: Power-law behaviour of Edisp tted to the form z x in the region from 6 to 20 a.u. Beyond 600 a.u. all chains exhibit dispersion energies with the z6 power law because of their nite length. The HOMOLUMO gaps, Eg , have been calculated from the KohnSham eigenvalues. Energies are in atomic units.

Eg 2.0 1.5 1.25 1.0 0.370 0.280 0.202 0.099

(H2 )16 x 4.89 4.58 4.20 3.52 Eg 0.366 0.270 0.183 0.057

(H2 )32 x 4.84 4.50 4.09 3.17

MULTIPOLE EXPANSION

We can understand these eects in terms of the multipole expansion. The multipole form of the dispersion energy is usually formulated in terms of (local) atomic polarizabilities. This leads to the usual R6 atomatom interaction (at leading order) and would not account for the anomalous dispersion power laws that we have observed in the chains. However the complete distributedpolarizability description is non-local: that is, it describes the change in multipole moments at one atom in response to a change in electrostatic elds at another. We obtain the multipole form of eq. (7) by expanding the Coulomb terms using the distributed form of the ab b a multipole expansion |r r |1 = Q t T tu Qu where a and b denote (atomic) sites and t and u are multipole indices 00 for the charge, 10, 11c and 11 s for the components a of the dipole, and so on. Q t is the multipole moment opab erator for moment t of site a and T tu are the interaction tensors that contain the distance and angular dependence (see Ref.[1] for details). Inserting this expansion in eq.

There are no assumptions made in deriving eq. (8), other than that spheres enclosing the atomic charge densities on dierent molecules do not overlap. We have calculated distributed polarizabilities using the constrained densitytting algorithm [20]. Eq. (8) involves a quadruple sum over sites and is therefore computationally demanding. To reduce this cost, we normally make a simplication by localizing the non-local polarizabilities. That is, the non-local po larizabilities the aa a are transtt () with a formed onto one or other site using the multipole expansion [21, 22]. This transformation is possible only if the non-local terms decay fast enough with inter-site distance. If not, the multipole expansion used in the transformation diverges and the localization is no longer possible. As we shall see, this is precisely what happens in the (H2 )n chains as the distortion parameter decreases. An important feature of the non-local polarizability description is the presence of charge-ow polarizabilities terms with t or u = 00 that describe the ow of charge in a moleculewhich are usually small and die o quickly with distance. The lowest rank charge-ow polarizabil a ity is aa 00,00 : if V is the potential at site a, the change in a = a aa (V a V a ). charge at site a is given by Q 00 00,00 For a system with a large HOMOLUMO gap the chargeow terms are expected to be short-range. The chargeow polarizabilities of the (H2 )32 chain with = 2.0 can be empirically modelled reasonably well with an exponential, that is,
|raa | aa , 00,00 e

(10)

where raa is the intersite distance and 0.5 a.u. In this case, the charge-ow polarizabilities drop by more than an order of magnitude within a few bonds. This is illustrated in g. 4. However, as we reduce the distortion parameter the charge-ow polarizabilities decay more and more slowly with site-site distance, until, at = 1.0, they span the entire length of the chain. Even for the = 1.5 chain, these non-local charge-ow polarizabilities can no

60 50 40

atomic site

30 20 10 0
60 50 40

atomic site

30 20 10 0 0 10 20 30 atomic site

40

50

60

10

20

30 atomic site

40

50

60

FIG. 4: Matrix representation of the charge-ow polarizability matrix aa 00,00 for (H2 )32 chains with = 2, 1.5, 1.25 and 1 (clockwise, from top left). Sites a and a are represented along the x and y axes. Because these terms span a number of orders of magnitude of both signs we have plotted ln(|aa 00,00 |). Colour scheme: Black rectangles correspond to terms of order 1 a.u. and white rectangles to terms of order 103 a.u.

longer be localized without incurring a signicant error, but for the chains with = 1.25 and 1.0 localization results in a qualitatively incorrect physical picture. Moreover the charge-ow polarizabilities contribute to the dipoledipole polarizability in the direction of the chain, and this contribution increases dramatically as the band-gap decreases. For a 64-atom H2 chain with = 2, the static dipole polarizability xx = 11c,11c parallel to the chain, calculated as described above, is about 410 a.u., of which 180 a.u. is contributed by charge-ows. When = 1, so that the H atoms are equally spaced, xx is about 11350 a.u., larger by about two orders of magnitude, and all of this increase is attributable to the charge-ow eects. We have calculated the dispersion energy using eq. (8) with distributed polarizabilities, including terms to rank 1. Higher ranking terms can be included, but these are not needed for the hydrogen chains. These results are presented in g. 5 for the (H2 )32 chains with = 2.0 and 1.0. First of all, consider the insulating chain ( = 2.0): the charge-ow polarizabilities alone severely underestimate the dispersion energy, but when terms up to rank 1 are included the agreement of eq. (8) with the non(2) expanded SAPT(DFT) value for Edisp is excellent for all inter-chain separations shown in the gure. However, at

FIG. 5: The second-order dispersion energy calculated using the polarizabilities of (H2 )32 chains with distortion parameters = 1 and 2. The contribution from the charge-ow polarizabil ities, aa 00,00 , is displayed separately from the sum of contributions up to the dipole-dipole polarizabilities.

= 1, the dispersion interaction is entirely dominated aa by the pure charge-ow (i.e. 00 ,00 ) terms, the higher order terms involving the dipole polarisabilities making a negligible contribution. Since the overall dispersion interaction is greatly enhanced as is reduced, this implies that these pure charge-ow terms enhance the dispersion

6 interaction in this limit. Why do the charge-ow terms result in the anomalous power laws? The lowest rank charge-ow terms, aa 00,00 , that appear in eq. (8) are associated with T -functions for ab 1 the charge-charge interaction: T 00 ,00 = Rab , where Rab is the distance between site a on one wire and site b on the other. If the wires are separated by distance vector R = (0, 0, z) and xa and xb are distance vectors for sites a and b along the chains then Rab = R (xa xb ). The dispersion energy arising from just the charge-ow terms is
(2) Edisp (00, 00) =

L z: In this limit both Rab and Ra b can be expanded in a multipole expansion about R = 1 1 (0, 0, z): R = |R xab |1 ab = |R (xa xb )| 1 2 3 1 1 z 2 xab z , and likewise for Ra b . Once again, the leading terms vanish because of the sum rule, leaving only the eective dipole-dipole contribution:
(2) Edisp (00, 00)

1 1 2 4z6

2 x2 ab xa b aa bb

bb aa 00,00 (iw)00,00 (iw)dw

1 2

aa bb

1 1 Rab Ra b

(11)

C6 (00, 00) . z6

(13)

bb aa 00,00 (iw)00,00 (iw)dw.

In contrast to the dipole-dipole polarizabilities, the charge-ow polarizabilities satisfy the sum-rule aa a t,00 (w) = 0, which is a direct consequence of the charge conservation requirement: A (r, r ; w)d3 r = 0. This leads to cancellation between the charge-ow contributions to the total dispersion energy, and to an R6 distance dependence at long range, but the cancellation is incomplete at short range, and terms in Rn , 2 n 5 also occur. To see how this arises, consider the following length scales: (1) Lc = 1/, which is a measure of the extent of the charge uctuations determined by the exponential decay of the charge-ow polarizabilities assumed in eq. (10), and (2) L, the chain length. z Lc : Here the charge-ow terms dominate and contribute R2 terms to the dispersion energy. From g. 4 we see that Lc is largest for the nearmetallic wire for which we see the eects of these terms (g. 5) to z 60 a.u. Lc z < L: In this region the extent of the charge uctuations is small compared with Rab and only those Ra b close to Rab are important. We can therefore expand Ra b about Rab using the multipole expansion. The leading term in this expansion is
(2) Edisp (00, 00)

This explains the large-z charge-ow contribution to the dispersion energy shown in g. 5. Notice that for the near-metallic wires, this contribution to the total molecular C6 coecient dominates that from the dipole-dipole polarizabilities, but for the largegap wires the opposite is true. In short, the charge-ow polarizabilities give rise to the changes in power-law of the dispersion energy and also contribute to an enhancement of the eective C6 coecient, applicable at long range.
DISCUSSION

The physical picture which emerges from both the Hckel approach and the ab initio calculations is the following. In systems with a nite but small gap, spontaneous charge-uctuations (plasmon modes in the innite case) introduce a secondary length-scale intermediate in size to the interatomic distance and the system size. This length scale grows as the gap gets smaller. For separations z small compared with this length scale, the dispersion energy arising from the correlated uctuations has metallic character. At z that are large compared with the length scale of the uctuations, the dispersion can be described using Londons dipole approximation, giving z5 behaviour, but the magnitude of the uctuations now depends strongly on the band gap, leading to orders of magnitude enhancement over the insulator case. In small-gap systems these uctuations give rise to a strong non-additivity in the polarizability. For example, the ratio of the longitudinal static polarizabilities for the near-metallic and insulating (H2 )32

=0

1 1 R ab Rab ab aa bb 00,00 (iw) 00,00 (iw) dw 1 2

(12)

where we have used the charge-ow sum rule. So we see that the R2 contributions sum to zero in this region. However, the higher-order contributions are non-zero.

7 chains is 28. The dispersion energy is proportional to the square of the polarizability, that is, 784. This is roughly the ratio of the dispersion energies for these two cases but only at separations z much greater than the chain length. Any attempt to extrapolate this result to shorter distances results in a severe overestimation of the dispersion energy. This eect is not a consequence of retardation (we use the non-relativistic Hamiltonian) or damping (charge-density overlap is negligible). It originates from the complex behaviour of the nonlocal charge-ow polarizabilities. These are terms of rank zero that describe charge uctuations in the system and are associated with a delocalized exchange-correlation hole [23, 24]. This delocalization can be quantied using the localization tensor [24, 25] and may give us a quantitative method for dening the charge uctuation length-scale, Lc . We are currently investigating this possibility. We have demonstrated that both the change in power-law with distance and the enhancement of the dispersion energy can be understood using nonlocal polarizability models containing charge-ow polarizabilities. It should come as no surprise that these eects are also strongly dependent on system size (see g.3). These results call into question theoretical methods that impose locality so as to scale linearly with system size (LCCSD(T), LMP2), or that approximate the dispersion energy using a pair-wise C6 R6 interaction, as is done with dispersion-corrected DFT methods and empirical potentials. In the former, these eects can be included by extending the region of locality, though at the cost of losing linearity in scaling, but the latter methods should not be applied to systems such as these. Even DFT functionals with a non-local dispersion correction, such as the van der Waals functional of Dion et al.[26] and the more recent functional of Vydrov and Van Voorhis [27] are unlikely to contain the correct physics because of an implicit assumption of locality in the polarizability. These functionals will include many-body non-additive eects between non-overlapping systems, but not the non-additive eects within each system such as those described here. On the other hand, methods based on the random phase approximation and quantum Monte Carlo should be able to describe the non-additive eects described in this paper if nite-size eects are kept under control. In systems containing carbon atoms, such as conjugated chains, the contributions from the core electrons will at least partially mask those from the more mobile electrons, so it is possible that the changes in power-law of the dispersion interaction will not be as dramatic as those we see in the (H2 )n chains. But we should nevertheless expect a high degree of non-additivity arising from the charge-ow terms. One indicator of this nonadditivity is the non-linear dependence of the molecular polarizability on system size. This has been already demonstrated above and is further supported by the experimental work of Compagnon et al. [28] on the polarizabilities of the fullerenes C70 and C60 which are found to be in the ratio 1.33 rather than the ratio 70/60 = 1.17 that would be expected if the systems were additive. Therefore, if reliable C6 -models are to be constructed for extended -conjugated systems, these models will need to absorb the eects of the non-additivity of the chargeow polarizabilities. That is, they should be tailored to suit the electronic structure of the system rather than be transferred from calculations on smaller systems. The WilliamsStoneMisquitta (WSM) procedure [2931] is one such method that is capable of constructing eective local polarizability and dispersion models that account for the electronic structure of the system. In fact, the non-local models presented in this paper are derived in the rst step of the WSM procedure. We are currently investigating the behaviour of WSM dispersion models for a variety of carbon systems. Finally, strongly delocalized systems will also have important contributions to the dispersion energy from terms of third order in the interaction operator, that is, from the hyperpolarizabilities. We have not considered such terms in this paper. Furthermore, in ensembles of such systems there will be strong non-additive eects between molecules. We are currently investigating the nature and importance of this type of non-additivity.

ACKNOWLEDGEMENTS

We thank University College London for computational resources. AJM and JS thank EPSRC for funding.

[1] A. J. Stone, The Theory of Intermolecular Forces (Clarendon Press, Oxford, 1996). [2] H. B. G. Casimir and D. Polder, Phys. Rev. 73, 360 (1948). [3] J. F. Dobson, Surface Science 601, 5667 (2007). [4] J. F. Dobson, A. White, and A. Rubio, Phys. Rev. Lett. 96, 073201(4) (2006). [5] J. F. Dobson, K. McLennan, A. Rubio, J. Wang, T. Gould, H. M. Le, and B. P. Dinte, Aust. J. Chem. 54, 513 (2001). [6] J. F. Dobson, J. Wang, B. P. Dinte, K. McLennan, and H. M. Le, Int. J. Quantum Chem. 101, 579 (2004).

8
[7] N. D. Drummond and R. J. Needs, Phys. Rev. Lett. 99, 166401(4) (2007). [8] J. Spencer, Ph.D. thesis, University of Cambridge (2009). [9] P. Epstein, Mathematische Annalen 56, 615 (1903). [10] P. Epstein, Mathematische Annalen 63, 205 (1906). [11] M. P. Tosi, Cohesion of ionic solids in the Born model (Academic Press, 1964), vol. XVI of Solid State Physics. [12] A. J. Misquitta, B. Jeziorski, and K. Szalewicz, Phys. Rev. Lett. 91, 33201 (2003). [13] A. Hesselmann, G. Jansen, and M. Schutz, J. Chem. Phys. 122, 014103 (2005). [14] A. J. Misquitta, R. Podeszwa, B. Jeziorski, and K. Szalewicz, J. Chem. Phys. 123, 214103 (2005). [15] H. C. Longuet-Higgins, Disc. Farad. Soc. 40, 7 (1965). [16] A. Hesselmann and G. Jansen, Phys. Chem. Chem. Phys. 5, 5010 (2003). [17] R. Podeszwa and K. Szalewicz, Chem. Phys. Lett. 412, 488 (2005). [18] B. Champagne, E. A. Perpete, S. J. A. van Gisbergen, E.-J. Baerends, J. G. Snijders, C. Soubra-Ghaoui, K. A. Robins, and B. Kirtman, J. Chem. Phys. 109, 10489 (1998). [19] H. Sekino, Y. Maeda, M. Kamiya, and K. Hirao, J. Chem. Phys. 126, 014107 (2007). [20] A. J. Misquitta and A. J. Stone, J. Chem. Phys. 124, 024111 (2006). [21] C. R. Le Sueur and A. J. Stone, Mol. Phys. 83, 293 (1994). [22] T. C. Lillestolen and R. J. Wheatley, J. Phys. Chem. A 111, 11141 (2007). [23] J. G. Angyan, J. Chem. Phys. 127, 024108 (2007). [24] J. G. Angyan, Int. J. Quantum Chem. 109, 2340 (2009). [25] R. Resta and S. Sorella, Phys. Rev. Lett. 82, 370 (1999). [26] M. Dion, H. Rydberg, E. Schroder, D. C. Langreth, , and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004). [27] O. A. Vydrov and T. V. Voorhis, Phys. Rev. Lett. 103, 063004 (2009). [28] I. Compagnon, R. Antoine, M. Broyer, P. Dugourd, J. Lerme, and D. Rayane, Phys. Rev. A 64, 025201 (2001). [29] A. J. Misquitta and A. J. Stone, J. Chem. Theory Comput. 4, 7 (2008). [30] A. J. Misquitta, A. J. Stone, and S. L. Price, J. Chem. Theory Comput. 4, 19 (2008). [31] A. J. Misquitta and A. J. Stone, Mol. Phys. 106, 1631 (2008).