Dispersion interactions between semiconducting wires

Conventionally, the dispersion interaction between two systems is formulated in terms of correlated fluctuations described by local (frequency-dependent) polarizabilities which gives rise to the familiar atom–atom $-C_{6}^{ab}R_{ab}^{-6}$ interaction at leading order Stone (1996). 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 HOMO–LUMO gaps.

It has been known for a very long time that the additive atom–atom 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 $R^{-3}$ power law Casimir and Polder (1948), but summing over atom–atom terms leads to a $R^{-4}$ 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 Dobson (2007); Dobson et al. (2006, 2001, 2004); Drummond and Needs (2007).

The semi-classical picture is useful for understanding the source of the differences 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 fluctuations give rise to instantaneous local dipoles, and the correlation between these local dipoles gives rise to the atom–atom $R^{-6}$ interaction. However, in a zero band-gap material, electronic fluctuations are long-ranged, particularly if one or two of the dimensions are nano-scale Dobson (2007). It is these long-ranged fluctuations that give rise to the non-additivity of the dispersion and deviations from the atom–atom $R^{-6}$ interaction.

While the insulating and metallic cases are now well understood, little is known of the intermediate, semi-conducting case. The nature of the dispersion interaction between extended molecules with finite but small HOMO–LUMO 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 $\pi$-conjugated systems. For want of a clearer understanding, many studies assume the usual insulating case for the atom–atom $R^{-6}$ form of this interaction. As we shall show in this article, this is both qualitatively and quantitatively incorrect for semiconducting molecules.

A word of clarification: we use the term ‘non-additive’ here to describe deviations from the additive atom-atom 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 effects here, except for a brief comment in the Discussion.

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 Stone (1996). In a mean-field 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:

where $i,j$ ($a,b$) are occupied (virtual) single-particle wavefunctions in either subsystem $A$ or $B$, and $\varepsilon_{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:

and is the Coulomb interaction between the charge density $\psi_{i}^{*}\psi_{a}$ due to excitation $i\rightarrow a$ in subsystem $A$ with $\psi_{j}^{*}\psi_{b}$ in $B$.

In periodic systems, the wavefunction can be expressed in Bloch form, i.e. $\psi_{j}(x)=e^{ik_{j}x}u_{j}(x)$, and so the co-density $\psi_{i}^{*}(x)\psi_{a}(x)=e^{i(k_{a}-k_{i})x}u_{i}^{*}(x)u_{a}(x)$ can have a long-wavelength modulation, whose electrostatic field 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 character and contribute significantly to $E^{(2)}_{\rm disp}$, as we discuss below.

Hückel (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

where $\beta,\beta^{\prime}$ are alternating bond-strengths between adjacent sites, as would be encountered in a chain of (H₂)_n or a $\pi$-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

and the set of wavevectors by

$\beta=\beta^{\prime}$ corresponds to a uniform wire with energy eigenvalues given by $\epsilon(k)=\pm 2\beta|\cos(kd)|$, which in the limit of large $n$ has a vanishing band-gap at half-filling (i.e. is a metal). The opposite limit ($\beta^{\prime}=0$) corresponds to a chain of isolated dimers with energy eigenvalues $\epsilon(k)=\pm\beta$, independent of wavevector, and corresponds to the perfect insulator. By varying the ratio $\beta^{\prime}/\beta$ between 0 and 1, we can very conveniently probe the dispersion interaction in the intermediate (semiconducting) regime, with the band-gap given by $\Delta E_{g}=2(\beta-\beta^{\prime})$. The required integrals can be evaluated Spencer (2009) as:

where $\{G\}$ are reciprocal lattice vectors of the primitive unit cell and $Y_{ia}(G)$ is the (analytic) Fourier transform of $u_{i}^{*}u_{a}$ obtained by assuming highly-localised basis functions on each atom. $K_{0}$ is the zero-th order modified Bessel function of the second kind and is the Fourier transform of the potential generated by a one-dimensional lattice at a field point at $z$ from the latticeEpstein (1903, 1906); Tosi (1964). $K_{0}$ 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 Monkhorst-Pack $k$-points to sample the Brillouin zone, which corresponds in real-space to wires consisting of 16802 sites per crystal cell. Such fine k-point sampling was required to obtained converged results with respect to system size.

As shown in fig. 1, the dispersion interaction between two uniform wires ($\beta^{\prime}/\beta=1$) varies as $\sim z^{-2}$ across the whole range of $z$, in good agreement with previously obtained results via RPADobson (2007) and QMCDrummond and Needs (2007) for metallic wires. The interaction between chains of isolated dimers ($\beta^{\prime}=0$) varies as $z^{-5}$, 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 ($z^{-2}$), whereas large $z$ it tends to insulating behaviour $z^{-5}$.

The principal drawback of the Hückel model is the lack of electron correlation, and hence screening, within each subsystem. For more realistic calculations we can use the symmetry-adapted perturbation theory based on density-functional theory Misquitta et al. (2003); Hesselmann et al. (2005); Misquitta et al. (2005) (SAPT(DFT)), where the dispersion energy is evaluated not via a sum-over-states but through a coupled Kohn–Sham formulation based on the density response functions of the interacting moleculesLonguet-Higgins (1965):

where $\alpha_{A/B}({\mathbf{r}},{\mathbf{r}}^{\prime};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 linear-response theory. The second-order dispersion energy is exact in this formulation if exchange effects are neglected.

We have studied the interactions between two parallel finite (H₂)_n (fig. 2) chains with $n=16$ and $32$ and distortion parameters $\eta=2.0,1.5,1.25$ and $1.0$, where $\eta$ is the ratio of the alternate bond lengths. The HOMO-LUMO gap of a hydrogen chain can be modified by simply introducing a difference in alternate bond lengths. The distortion parameter $\eta$ 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.

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 Møller-Plesset perturbation theory and rivals that of coupled-cluster methods Misquitta et al. (2005); Hesselmann and Jansen (2003); Hesselmann et al. (2005); Podeszwa and Szalewicz (2005). 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 Kohn–Sham perturbation theory Champagne et al. (1998); Sekino et al. (2007). All calculations used the cc-pVDZ basis which will result in a significant 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.

Dispersion energies are displayed in fig. 3 and the effective power laws for the physically important separations between 6 and 20 a.u. are shown in table 1 together with HOMO–LUMO gaps. For insulating chains with an additive dispersion interaction we would expect two regimes determined by chain length $L$: for $z\ll L$ we should expect the infinite chain result, i.e., the effective dispersion interaction should decay as $z^{-5}$, and for $z\gg L$ we should recover the usual $z^{-6}$ power law. This is what we see with the chains with the largest distortion parameter, $\eta=2.0$, which exhibit large HOMO–LUMO gaps. As the distortion parameter $\eta$ approaches unity and the HOMO–LUMO 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 significantly enhanced. For example at $40$ a.u. (roughly half the chain length) an $\eta=1$ chain has a dispersion interaction two orders of magnitude larger than an $\eta=2$ chain. Additionally, we see increasing finite-size effects: the power laws for the $n=16$ and $32$ chains, which were very similar for $\eta=2.0$, are considerably different for $\eta=1.0$.

We can understand these effects 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 $R^{-6}$ atom–atom interaction (at leading order) and would not account for the anomalous dispersion power laws that we have observed in the chains. However the complete distributed-polarizability description is non-local: that is, it describes the change in multipole moments at one atom in response to a change in electrostatic fields at another.

We obtain the multipole form of eq. (7) by expanding the Coulomb terms using the distributed form of the multipole expansion $|{\mathbf{r}}-{\mathbf{r}^{\prime}}|^{-1}=\hat{Q}^{a}_{t}T^{ab}_{tu}\hat{Q}^{b}_{u}$ where $a$ and $b$ denote (atomic) sites and $t$ and $u$ are multipole indices — 00 for the charge, 10, $11c$ and $11s$ for the components of the dipole, and so on. $\hat{Q}^{a}_{t}$ is the multipole moment operator for moment $t$ of site $a$ and $T^{ab}_{tu}$ are the interaction tensors that contain the distance and angular dependence (see Ref.Stone (1996) for details). Inserting this expansion in eq. (7) we obtain the multipole form for the dispersion energy:

Here $T^{ab}_{tu}$ is the interaction function between multipole $t$ on site $a$ in subsystem $A$ and multipole $u$ on site $b$ in $B$, while $\alpha^{aa^{\prime}}_{tt^{\prime}}$ are the frequency-dependent non-local polarizabilities for sites $a$ and $a^{\prime}$, and may be expressed in terms of the frequency-dependent density susceptibility asMisquitta and Stone (2006)

There are no assumptions made in deriving eq. (8), other than that spheres enclosing the atomic charge densities on different molecules do not overlap. We have calculated distributed polarizabilities using the constrained density-fitting algorithm Misquitta and Stone (2006).

Eq. (8) involves a quadruple sum over sites and is therefore computationally demanding. To reduce this cost, we normally make a simplification by localizing the non-local polarizabilities. That is, the non-local polarizabilities — the $\alpha^{aa^{\prime}}_{tt^{\prime}}(\omega)$ with $a\neq a^{\prime}$ — are transformed onto one or other site using the multipole expansion Le Sueur and Stone (1994); Lillestolen and Wheatley (2007). 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 (H₂)_n chains as the distortion parameter $\eta$ decreases.

An important feature of the non-local polarizability description is the presence of charge-flow polarizabilities—terms with $t$ or $u=00$ that describe the flow of charge in a molecule—which are usually small and die off quickly with distance. The lowest rank charge-flow polarizability is $\alpha^{aa^{\prime}}_{00,00}$: if $V^{a}$ is the potential at site $a$, the change in charge at site $a$ is given by $\Delta\hat{Q}^{a}_{00}=-\sum_{a^{\prime}}\alpha^{aa^{\prime}}_{00,00}(V^{a^{\prime}}-V^{a})$. For a system with a large HOMO–LUMO gap the charge-flow terms are expected to be short-range. The charge-flow polarizabilities of the (H₂)₃₂ chain with $\eta=2.0$ can be empirically modelled reasonably well with an exponential, that is,

where $r_{aa^{\prime}}$ is the intersite distance and $\gamma\approx 0.5$ a.u. In this case, the charge-flow polarizabilities drop by more than an order of magnitude within a few bonds. This is illustrated in fig. 4. However, as we reduce the distortion parameter $\eta$ the charge-flow polarizabilities decay more and more slowly with site-site distance, until, at $\eta=1.0$, they span the entire length of the chain. Even for the $\eta=1.5$ chain, these non-local charge-flow polarizabilities can no longer be localized without incurring a significant error, but for the chains with $\eta=1.25$ and $1.0$ localization results in a qualitatively incorrect physical picture.

Moreover the charge-flow polarizabilities contribute to the dipole–dipole polarizability in the direction of the chain, and this contribution increases dramatically as the band-gap decreases. For a 64-atom H$2$ chain with $\eta=2$, the static dipole polarizability $\alpha_{xx}=\alpha_{11c,11c}$ parallel to the chain, calculated as described above, is about 410 a.u., of which 180 a.u. is contributed by charge-flows. When $\eta=1$, so that the H atoms are equally spaced, $\alpha_{xx}$ is about 11350 a.u., larger by about two orders of magnitude, and all of this increase is attributable to the charge-flow effects.

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 fig. 5 for the (H₂)₃₂ chains with $\eta=2.0$ and $1.0$. First of all, consider the insulating chain ($\eta=2.0$): the charge-flow polarizabilities alone severely underestimate the dispersion energy, but when terms up to rank 1 are included the agreement of eq. (8) with the non-expanded SAPT(DFT) value for $E^{(2)}_{\rm disp}$ is excellent for all inter-chain separations shown in the figure. However, at $\eta=1$, the dispersion interaction is entirely dominated by the pure charge-flow (i.e. $\alpha^{aa^{\prime}}_{00,00}$) terms, the higher order terms involving the dipole polarisabilities making a negligible contribution. Since the overall dispersion interaction is greatly enhanced as $\eta$ is reduced, this implies that these pure charge-flow terms enhance the dispersion interaction in this limit.

Why do the charge-flow terms result in the anomalous power laws? The lowest rank charge-flow terms, $\alpha^{aa^{\prime}}_{00,00}$, that appear in eq. (8) are associated with $T$-functions for the charge-charge interaction: $T^{ab}_{00,00}=R_{ab}^{-1}$, where $R_{ab}$ is the distance between site $a$ on one wire and site $b$ on the other. If the wires are separated by distance vector ${\mathbf{R}}=(0,0,z)$ and ${\mathbf{x}}_{a}$ and ${\mathbf{x}}_{b}$ are distance vectors for sites $a$ and $b$ along the chains then ${\mathbf{R}}_{ab}={\mathbf{R}}-({\mathbf{x}}_{a}-{\mathbf{x}}_{b})$. The dispersion energy arising from just the charge-flow terms is

In contrast to the dipole-dipole polarizabilities, the charge-flow polarizabilities satisfy the sum-rule $\sum_{a^{\prime}}\alpha^{aa^{\prime}}_{t,00}(w)=0$, which is a direct consequence of the charge conservation requirement: $\int\alpha_{A}({\mathbf{r}},{\mathbf{r}}^{\prime};w){\mathrm{d}}^{3}{\mathbf{r}^{\prime}}=0$. This leads to cancellation between the charge-flow contributions to the total dispersion energy, and to an $R^{-6}$ distance dependence at long range, but the cancellation is incomplete at short range, and terms in $R^{-n},2\leq n\leq 5$ also occur. To see how this arises, consider the following length scales: (1) $L_{c}=1/\gamma$, which is a measure of the extent of the charge fluctuations determined by the exponential decay of the charge-flow polarizabilities assumed in eq. (10), and (2) $L$, the chain length.

• •

  $z\leq L_{c}$: Here the charge-flow terms dominate and contribute $R^{-2}$ terms to the dispersion energy. From fig. 4 we see that $L_{c}$ is largest for the near-metallic wire for which we see the effects of these terms (fig. 5) to $z\sim 60$ a.u.

• •

  $L_{c}\ll z\mathchar 12604\relax L$: In this region the extent of the charge fluctuations is small compared with $R_{ab}$ and only those $R_{a^{\prime}b^{\prime}}$ close to $R_{ab}$ are important. We can therefore expand $R_{a^{\prime}b^{\prime}}$ about $R_{ab}$ using the multipole expansion. The leading term in this expansion is

  	$\displaystyle E^{(2)}_{\rm disp}(00,00)\approx-\frac{1}{2\pi}\sum_{ab}\frac{1}{R_{ab}}\frac{1}{R_{ab}}$
  	$\displaystyle\qquad\times\int_{0}^{\infty}\left(\sum_{a^{\prime}}\alpha^{aa^{\prime}}_{00,00}({\rm i}w)\right)\left(\sum_{b^{\prime}}\alpha^{bb^{\prime}}_{00,00}({\rm i}w)\right)dw$
  	$\displaystyle=0$		(12)

  where we have used the charge-flow sum rule. So we see that the $R^{-2}$ contributions sum to zero in this region. However, the higher-order contributions are non-zero.

• •

  $L\ll z$: In this limit both $R_{ab}$ and $R_{a^{\prime}b^{\prime}}$ can be expanded in a multipole expansion about ${\mathbf{R}}=(0,0,z)$: $R_{ab}^{-1}=|{\mathbf{R}}-({\mathbf{x}}_{a}-{\mathbf{x}}_{b})|^{-1}=|{\mathbf{R}}-{\mathbf{x}}_{ab}|^{-1}\approx z^{-1}-\frac{1}{2}x_{ab}^{2}z^{-3}$, and likewise for $R_{a^{\prime}b^{\prime}}^{-1}$. Once again, the leading terms vanish because of the sum rule, leaving only the effective dipole-dipole contribution:

  	$\displaystyle E^{(2)}_{\rm disp}(00,00)$	$\displaystyle\approx-\frac{1}{2\pi}\frac{1}{4z^{6}}\sum_{aa^{\prime}}\sum_{bb^{\prime}}x_{ab}^{2}x_{a^{\prime}b^{\prime}}^{2}$
  		$\displaystyle\times\int_{0}^{\infty}\alpha^{aa^{\prime}}_{00,00}({\rm i}w)\alpha^{bb^{\prime}}_{00,00}({\rm i}w)dw$
  		$\displaystyle\equiv-\frac{C_{6}(00,00)}{z^{6}}.$		(13)

  This explains the large-$z$ charge-flow contribution to the dispersion energy shown in fig. 5. Notice that for the near-metallic wires, this contribution to the total molecular $C_{6}$ coefficient dominates that from the dipole-dipole polarizabilities, but for the large-gap wires the opposite is true.

In short, the charge-flow polarizabilities give rise to the changes in power-law of the dispersion energy and also contribute to an enhancement of the effective $C_{6}$ coefficient, applicable at long range.

The physical picture which emerges from both the Hückel approach and the ab initio calculations is the following.

• •

  In systems with a finite but small gap, spontaneous charge-fluctuations (plasmon modes in the infinite 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 fluctuations has metallic character.

• •

  At $z$ that are large compared with the length scale of the fluctuations, the dispersion can be described using London’s dipole approximation, giving $z^{-5}$ behaviour, but the magnitude of the fluctuations now depends strongly on the band gap, leading to orders of magnitude enhancement over the insulator case.

• •

  In small-gap systems these fluctuations 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 (H₂)₃₂ 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 effect 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 non-local charge-flow polarizabilities. These are terms of rank zero that describe charge fluctuations in the system and are associated with a delocalized exchange-correlation hole Angyan (2007, 2009). This delocalization can be quantified using the localization tensor Resta and Sorella (1999); Angyan (2009) and may give us a quantitative method for defining the charge fluctuation length-scale, $L_{c}$. 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 non-local polarizability models containing charge-flow polarizabilities.

• •

  It should come as no surprise that these effects are also strongly dependent on system size (see fig.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 $-C_{6}R^{-6}$ interaction, as is done with dispersion-corrected DFT methods and empirical potentials. In the former, these effects 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.Dion et al. (2004) and the more recent functional of Vydrov and Van Voorhis Vydrov and Voorhis (2009) 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 effects between non-overlapping systems, but not the non-additive effects 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 effects described in this paper if finite-size effects are kept under control.

In systems containing carbon atoms, such as $\pi$-conjugated chains, the contributions from the core electrons will at least partially mask those from the more mobile $\pi$ 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 (H₂)_n chains. But we should nevertheless expect a high degree of non-additivity arising from the charge-flow terms. One indicator of this non-additivity 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. Compagnon et al. (2001) 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 $C_{6}$-models are to be constructed for extended $\pi$-conjugated systems, these models will need to absorb the effects of the non-additivity of the charge-flow 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 Williams–Stone–Misquitta (WSM) procedure Misquitta and Stone (2008a); Misquitta et al. (2008); Misquitta and Stone (2008b) is one such method that is capable of constructing effective 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 first 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 effects between molecules. We are currently investigating the nature and importance of this type of non-additivity.

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