1
Exact exchange and Wilson–Levy correlation: a pragmatic device for studying complex weakly-bonded systems

2
The Wilson–Levy (WL) correlation functional is used together with Hartree–Fock (HF) theory to evaluate interaction energies at intermediate separations (i.e. around equilibrium separation) for several weakly-bonded systems.

3
The HF+WL approach reproduces binding trends for all complexes studied: selected rare-gas dimers, isomers of the methane dimer, benzene dimer and naphthalene dimer, and base-pair stacking structures for pyrimidine, cytosine, uracil and guanine dimers.

4
These HF+WL data are contrasted against results obtained from some popular functionals (including B3LYP and PBE), as well as two newly-developed functionals, X3LYP and xPBE.

5
The utility of HF+WL, with reference to exact-exchange (EXX) density-functional theory, is discussed in terms of a suggested EXXWL exchange–correlation functional.

Introduction

6
The interaction energy between two weakly-interacting fragments comprises many distinct and well-understood (although perhaps not easy to calculate) contributions.

7
Following Stone,1 the interaction energy can be decomposed into first-order, second-order and higher-order contributions.

8
Hartree–Fock (HF) theory adequately captures (at small and large separations) first-order contributions such as the electrostatic interaction and the exchange–repulsion interaction.

9
The performance of HF theory for second-order contributions is less comprehensive.

10
Although HF theory can yield a reasonable estimate of the induction contribution, HF theory cannot recover the dispersion energy.

11
Higher-order terms are recovered (with varying degrees of accuracy) via use of post-HF methods such as Møller–Plesset second-order perturbation theory, MP2, and coupled-cluster singles and doubles with triples corrections, CCSD(T).

12
Both induction and dispersion are present in the regime of large fragment separations (the ‘long-range’ form) where the charge distributions on each fragment do not overlap, and are also present in the small to intermediate separation regime (the ‘short-range’ form) where the charge distributions do overlap.

13
Where charge distributions overlap, second-order contributions are modified by general overlap effects as well as electron exchange, leading to damped induction, exchange-induction, damped dispersion and exchange–dispersion contributions.

14
These short-range terms should behave as a function of overlap between fragments.

15
Therefore, in this article, the general phrase ‘dispersion energy’ does not solely refer to the London dispersion (with r−6 dependence) that characterises the long-range contribution, rather, it refers to the dispersion interaction at all separations.

16
The reader is directed to the body of work concerned with symmetry-adapted perturbation theory (SAPT), where such contributions to the interaction energy are rigorously defined (e.g., the review article of Jeziorski et al2)..

17
Accurate calculation of these second-order contributions, at any separation, for large and/or complex weakly-interacting systems remains a challenging prospect in theoretical chemistry.

18
In the worst case, the use of post-HF methods such as MP2 or CCSD(T) can be limited at best to rather modest basis set sizes.

19
These methods may yield unreliable results for such small basis sets.

20
For example, even with use of large basis sets, MP2 can substantially overestimate the attraction in weakly-bound aromatic systems.3

21
However, recent progress in the development of algorithms for the rapid evaluation of MP2 energies and gradients (e.g. see Saebø et al4)., should in future enable practical geometry optimisations using larger basis sets.

22
Second-order contributions to the interaction energy arise from electronic correlation effects (beyond electron exchange alone).

23
Density-functional theory (DFT) is often proposed as a means of obtaining the correlation contribution to the total electronic energy, in a way that is computationally economical.

24
However, it is now well-known that no mainstream, practical correlation functional can recover the r−6 dependence of London dispersion at large fragment separations (see for example refs. 7 and 25, and references therein).

25
This fact will not be disputed here; lack of London dispersion remains a major limitation of all approaches that involve Kohn–Sham DFT in some way.

26
However, many studies have amply demonstrated that application of DFT to weakly-bound systems can not yield globally reliable estimates of the interaction energy even at intermediate intermonomer separations5–9 where there is overlap between fragments.

27
Several groups have made progress in the development of ‘seamless DFT’ approaches that can describe dispersion contributions at all separations10–12 (including in the overlapping regime).

28
However, all of these approaches are at present too computationally intensive to be applied to large and/or complex weakly-bonded systems.

29
On the basis of Heitler–London calculations, Cybulski and Seversen13 correctly identified that the correlation functionals used in their study only showed binding when it was due to exchange effects from the correlation functional, such that the total intermolecular exchange was negative (for their decomposition scheme).

30
Furthermore, they noted that, as expected, this intermolecular exchange contribution vanishes faster with respect to separation than r−6.

31
This overlap effect has been observed (but not quantified) by others, e.g. see .ref. 7

32
The decomposition analysis used by Cybulski and Seversen, while helpfully demonstrating the overlap effect of binding for some correlation functionals, does not explain why these correlation functionals yield negative contributions to the intermolecular exchange.

33
While it is often the correlation functional which is identified as deficient in recovering second-order contributions in the overlapping regime, several studies have instead identified the exchange functional as problematic,6,8,13–16 especially under the conditions of low-density/high-density-gradient.

34
An accurate description of exchange under such conditions is vital for describing weakly-bonded systems.

35
Empirical scaling of DFT exchange has also been examined by several groups,14,16,17 with a view to correcting this problem.

36
However, the results of these studies do not show transferability of the scaling factor across different systems.

37
Therefore, the fact that mainstream DFT cannot recover short-range second-order contributions must be due, in part, to problems with the exchange functional.

38
Despite these well-known deficiencies of DFT, a number of schemes for empirically correcting DFT calculations by addition of damped dispersion (DFT+D) have been proposed.6,18–20

39
As an alternative to post-HF theories and also DFT+D, another route to gaining economical (and reasonable) estimates of the interaction energy for large weakly-bound systems is to use a Hartree–Fock (HF) plus dispersion approach (HF+D).

40
The use of HF implies that, for a given density, the exchange contribution will be exact (by definition), including under low-density/high-density-gradient conditions.

41
The HF+D approach was first introduced by Hepburn etal.,21 has been successfully applied to small weakly-bound systems,22 and recently Lim and coworkers have applied HF+D to small clusters of polyaromatic heterocycles to good effect.23,24

42
One issue to note with both HF+D and DFT+D is that these methods require good quality dispersion coefficients.

43
Also, the form of the damping function, including anisotropic damping, can be a very important factor in gaining good results from these DFT+D models.

44
While some studies have reported progress in this area,18,19 further studies into damping are probably warranted.

45
Both of these ‘+dispersion’ approaches may succeed for the long-range part of the dispersion interaction, but it is not entirely clear that the DFT+D approach is correct in the overlapping regime.25

46
The problem of calculating the short-range interaction energy can also be approached via the addition of a correlation contribution derived from DFT onto the total HF energy of the system, bearing in mind the deficiencies pointed out by Cybulski and Seversen.13

47
As in the HF+D approach, the exchange contribution should behave correctly at long range.

48
This HF+DFT approach was shown by Pérez-Jordá etal26. to work very well for selected rare-gas dimers in the case where the Wilson–Levy correlation functional (WL)27 was used.

49
Some other correlation functionals also showed binding, but did not recover qualitatively correct equilibrium binding energies.

50
Moreover, this study also showed that in all cases where binding was found, the interaction energy curve decayed to zero faster than r−6, clearly not capturing any long-range dispersion contribution (as expected).

51
In addition, several correlation functionals did not yield any attractive contributions, notably the VWN28 and Perdew–Wang29 correlation functionals.

52
This HF+DFT approach is not the same as Hartree–Fock Kohn–Sham (HFKS) theory, but Pérez-Jordá et al26. reported that the HF+DFT results for rare-gas dimers were almost identical to the HFKS results.

53
In particular, these authors demonstrated that the Wilson–Levy correlation functional (WL)27 worked extremely well in partnership with HF theory, in reproducing the experimental well-depths of the rare-gas dimer systems studied “with remarkable agreement”.

54
Their data, shown in Fig. 1, emphatically reinforces their statement.

55
This approach will herein be referred to as HF+WL.

56
The WL correlation functional is based on the Wigner correlation functional30 (W38), contains gradient corrections, and includes four empirically-determined parameters.

57
The parameters were determined via comparison with the exact (calculated) correlation energy (in the traditional quantum chemistry definition) for several atoms (neutral and cationic), plus the experimental correlation energy of the helium atom.

58
Hartree–Fock densities of these atoms were used in the parametrisation of this functional.

59
As noted by Wilson and Levy,27 since the traditional quantum chemistry definition of correlation energy is not always equivalent to the DFT definition, the WL functional may not be successful if applied to systems where the HF density is substantially different to the true density of the system.

60
Since the WL correlation functional was designed for use with HF densities, it may be an appropriate approximate correlation functional to use with exact exchange methods, at least for studying pure van der Waals interactions.

61
However, the utility of this functional has principally been demonstrated only as a post-HF correction8,31,32 for atoms and molecules, not as a dynamic correlation contribution that is self-consistently evaluated.

62
The performance of WL correlation for covalent bonds has not been as extensively tested.

63
Chermette et al33. reported self-consistent KS-DFT optimisations of the CO and N2O molecules, using the B88 exchange functional34 with WL correlation.

64
Their data indicate that, compared with other correlation functionals used with B88, WL underestimates covalent bond lengths.

65
The HF+WL method should essentially be regarded as an empirical correction to the HF interaction energy, and presumably has inherent deficiencies with the description of intermolecular exchange.

66
It is the applicability of this approach as a pragmatic device that will be further explored in this work.

67
Here, the HF+WL approach is applied to all homoatomic and heteroatomic rare-gas dimers RgKr {Rg = He, Ne, Ar, Kr}, the methane dimer, benzene dimer, naphthalene dimer and also DNA base-pair stacking interactions.

68
Some comments on the choice of example systems studied here are warranted.

69
It is noted that, in general, any perturbative correction to the HF interaction energy should only hold under conditions such that the HF wavefunction will not substantially distort if (some) second-order effects are dynamically included.

70
Therefore, the method has been tested on systems where dynamic inclusion of such effects would not substantially distort the wavefunction, compared with the HF wavefunction.

71
Systems exhibiting pure van der Waals interactions (dispersion and repulsion) are ideal in this case, and the systems chosen here reflect this.

72
In general, all HF + dispersion methods are subject to this condition, which does limit the applicability of this type of approach.

73
However, it appears possible that the HF+WL method can be extended such that all types of non-bonded system could be treated in future.

74
Specifically, it is conjectured that use of the WL correlation functional, applied self-consistently in conjunction with exact-exchange DFT (e.g. the optimised effective potential, OEP35), denoted here as the EXXWL functional, should have a performance similar to the HF+WL approach.

75
The EXXWL functional would open up possibilities for qualitatively accurate calculations on complex systems not yet accessible by other methodologies, e.g. physisorption onto nanotubes, gas–surface interactions, etc.

76
The example systems selected in this work examine aromatic ring interactions.

77
Such interactions are ubiquitous in Nature, and should be described properly by any theory that claims to model ‘biological systems’.

78
Finally, it must be emphasised that, as demonstrated by Pérez-Jordá et al.,26 the HF+WL method has no pretensions to capturing London dispersion (the long-range contribution).

79
However, this caveat does not necessarily exclude the HF+WL method yielding binding when there is overlap between fragments, crucially, around the equilibrium separation for many weakly-interacting systems.

80
Having accepted the caveats of HF+WL, the utility of this approach as a pragmatic device for locating low-energy structures of large van der Waals systems (some of which have low symmetry) in a simple and quick way will be demonstrated.

81
These HF+WL results will be compared with results obtained using several established mainstream functionals (e.g. B3LYP and PBE), as well as two newly-developed functionals, X3LYP36 and xPBE.37

82
It has been recently suggested that both of these functionals should adequately describe non-bonded interactions for systems involving the first ten elements in the periodic table36,37.

Computational methods

83
All calculations in this work were done using GAUSSIAN.9838

84
The monomers were treated as rigid bodies throughout.

85
Unless otherwise stated, the monomer structures were fully optimised at the MP2/6-31G** level of theory.

86
For a given arrangement of rigid monomers, the HF wavefunction was converged and the total electronic energy (EHF) obtained, for both the complex, and also for each monomer in the complex.

87
The monomer wavefunctions and energies were obtained in the full set of dimer basis functions, such that the interaction energies were corrected for basis-set superposition error (BSSE) via the counterpoise correction39 (CP).

88
The correlation energy was then evaluated using the corresponding HF density.

89
The integration was done over a pruned grid with 99 radial shells and 590 angular points per shell.

90
The quadrature used here ensured that the correlation energy was converged to within 1–2 μEh per atom (depending on atom type).

91
The WL correlation energy (Ecorr) was then extracted from the calculation for both the complex and individual monomers (again using the CP correction).

92
The interaction energy was calculated using the supermolecule method, in the following way:The performance of two newly-developed functionals, X3LYP36 and xPBE,37 was also investigated here, in the case of the molecular dimers.

93
Both of these functionals were coded into GAUSSIAN98, and satisfactorily verified against the published data for He2,36,37 Ne236,37 and Ar237 systems.

94
The basis set dependence of the HF+WL energies was examined.

95
This was done to confirm the dependence of the interaction energy on overlap, for short- to medium-range separations.

96
The interaction energy as calculated here by the supermolecule method is of course not variational, so we do not necessarily expect to see convergence of this energy from above.

97
Furthermore, since the HF+WL method used here is an ad hoc post-Hartree–Fock correction (and not variationally obtained), this will cause further problems when seeking monotonically decreasing behaviour of interaction energies as a function of basis-set size.

98
In all cases, a range of basis sets, chosen from 6-31G*, 6-31G*(0.25),40 6-311G*, 6-31 + G* and the aug-cc-pVXZ41–43 {X = D, T, Q, 5} family were used.

99
The HF+WL interaction energy was calculated as a function of separation in all cases, and the minimum interaction energy was reported in each case.

100
For the methane dimer, three different relative orientations were examined, as shown in Fig. 2.

101
For the benzene dimer, three geometries were considered: a T-shaped structure (C2v symmetry), a slipped-parallel (SP) structure (C2h symmetry) and an undisplaced parallel (F) structure (D6h symmetry), also shown in Fig. 2.

102
For the naphthalene dimer, four low-energy structures and one high-energy structure were examined: two crossed structures (D2d and C2 symmetries) and two slipped-parallel structures (C2h and Ci symmetries), and one undisplaced parallel structure (D2h).

103
The four low-energy structures are shown in Fig. 3.

104
Optimised naphthalene dimer geometries were taken from previous work.44

105
In this case, the HF+WL interaction energies were calculated as a function of separation (along the vector that connects centre of mass points), while the relative orientations of each dimer structure were kept fixed.

106
Finally, for the base-pair stacking dimers, the geometries were taken from Hobza and Šponer.45

107
That is, antiparallel undisplaced structures were used for the pyrimidine and cytosine dimers, and the vertical displacement was optimised using HF+WL.

108
Compared with pyrimidine and cytosine, the procedure for creating the input structures for uracil and guanine was a little different.

109
Following Hobza and Šponer, the uracil dimer geometry was optimised in this work at the MP2/6-31G** level, and agrees closely with the structural data of Hobza and Šponer.46

110
Using this MP2 structure, the separation (along the vector connecting the centre of mass of each molecule) was optimised using HF+WL, while keeping the angle of the monomer normal vectors (normal to the ring plane) constant.

111
The guanine dimer structure was an undisplaced stacked geometry, where one monomer was rotated clockwise about its centre of mass (in the plane of the molecule).

112
The optimum twist angle was identified at the HF+WL/aug-cc-pVDZ level, for a fixed vertical separation of 3.4 Å.

113
The vertical separation (for the fixed optimum twist angle) was then optimised for this structure using HF+WL.

114
The geometries used for the base-pair calculations can be found in Fig. 4.

Results and discussion

115
The calculated well-depths and minimum energy separations, obtained using the aug-cc-pV5Z basis, for Kr2, HeKr, NeKr and ArKr are presented in Table 1, along with reference values from Ogilvie and Wang,47,48 and other DFT values from Patton and Pederson.9

116
The interaction energies are also plotted in Fig. 1.

117
In agreement with the previous findings of Pérez-Jordá et al.,26 the trend of the well depths presented here are in very good qualitative agreement with the experimental results, while not being quantitatively accurate.

118
The minimum energy separation is consistently underestimated by not more than 0.5 a0, as is also demonstrated in the results of Pérez-Jordá et al. The HF+WL approach has not captured the slight contraction in separation in going from KrHe to KrNe.

119
Examination of these data plus the data of Pérez-Jordá et al26. show that, with the exception of dimers containing helium, the HF+WL underestimates the interaction energy.

120
The interaction energy for each system was also calculated for separations beyond the minimum-energy point.

121
These calculations confirmed what has been reported previously by Pérez-Jordá et al.,26 namely that the interaction energy decays to zero too rapidly compared with the r−6 decay observed with long-range London dispersion interactions.

122
These results re-emphasise that weak interactions for non-overlapping fragments cannot be captured by the HF+WL approach.

123
However, these data also suggest that failure to recover long-range London dispersion does not preclude a DFT-based methodology from recovering binding at distances where fragments do overlap.

124
As shown in Fig. 1, the qualitative agreement found using HF+WL clearly exceeds the performance of any other mainstream DFT approach (demonstrated, for example, in the results of Patton and Pederson.9)

125
Fig. 1 also re-emphasises the point made by others that the PW91 and PBE functionals do not show great variation in energy between each rare-gas dimer system.

126
The PW91 data fortuitously coincides with the HeKr and NeKr points, while the PBE data fortuitously coincides with the ArKr and Kr2 points.

127
Neither can be said to describe rare-gas equilibrium interaction energies satisfactorily, while the HF+WL equilibrium interaction energies follow the trend of the Ogilvie and Wang data relatively closely.

128
The interaction energies of these rare-gas systems, calculated using HF+WL for a range of basis sets aug-cc-pVXZ {X = D, T, Q, 5} is shown in Fig. 5.

129
For each basis the separation between atoms was optimised, but in most cases did not change appreciably for the different basis sets.

130
It seems clear that the energy has converged by the aug-cc-pVTZ basis, and that the aug-cc-pVDZ basis still follows a reasonable qualitative trend (especially when compared with the PBE and PW91 data), although the trend for HeKr and NeKr is not captured.

131
The basis set dependence was also checked for the smaller 6-31G* basis.

132
For all systems but KrNe, the smaller basis set performed no worse than the aug-cc-pV5Z basis, in terms of capturing the correct trend.

133
However, the disagreement for HeNe is striking – the 6-31G* basis recovering an interaction energy of only −5.5 K. The optimum seperation was increased (to 7.0 a0).

134
As will be confirmed with some of the other systems, the overlap obtained using the 6-31G* basis appears insufficient to ensure confident use of HF+WL here.

135
Clearly the HF+WL approach is not perfect here – but it is outperforming other DFT-based approaches.

136
Three structures of the methane dimer were studied, as shown in Fig. 2; one D3d structure and two D3h structures.

137
The D3d structure has been identified by Tsuzuki and coworkers49 as the global minimum structure, at least at the MP3/aug(df,pd)-6-311G** level of theory.

138
The D3d structure corresponds with isomer H, and the two D3h structures as isomers A and G, in the nomenclature of Tsuzuki et al.49

139
Isomers A and G are other low-energy structures, and therefore these two geometries probe the interaction energy in regions of configuration space away from the global minimum.

140
Separations and interaction energies for these three structures, calculated using the HF+WL method, as well as the X3LYP and xPBE functionals, are shown in Table 2, along with the ab initio data of Tsuzuki et al.49

141
The aug-cc-pVTZ basis was used here for HF+WL, X3LYP and xPBE.

142
These data demonstrate that the HF+WL approach captures the binding trends for these three structures, for all basis sets considered, while clearly not displaying quantitative accuracy.

143
Again, the HF+WL approach underestimates the minimum-energy separation, found to be 4.6, 3.5 and 3.4 Å for the A, G and H structures, respectively.

144
For comparison, the MP3 separations from Tsuzuki et. al. are 4.8, 3.8 and 3.8 Å for the A, G and H structures, respectively.

145
In contrast, the X3LYP and xPBE data seem less reliable here.

146
Although the xPBE data show binding for each structure, and also that the trend of the interaction energies seems correct (A > G > H), the binding is reduced compared with the MP3 results (except in the case of isomer A).

147
The range of xPBE interaction energies is significantly compressedwith respect to the MP3 data.

148
The X3LYP data show binding only for structure A. Out to a separation of 6 Å, the other two structures were not bound using the X3LYP functional.

149
Also in contrast with the rare-gas systems, the isomers studied here showed little sensitivity to changes in basis set using the HF+WL approach.

150
For each basis set used (6-31G*, 6-31G*(0.25), aug-cc-pVXZ {X = D, T, Q}) the trend of interaction energies (A > G > H) was always followed.

151
The minimum energy separations and interaction energies for the three benzene dimer structures, calculated at the HF+WL/aug-cc-pVDZ level, are shown in Tables 3 and 4, respectively.

152
Again, the D6h (F) structure provides a probe of the HF+WL method away from the two low-energy structures.

153
For comparison, the data of Tsuzuki et al50. are shown.

154
These data were obtained using their extrapolation method where the interaction energy is calculated at the MP2 basis set limit plus a CCSD(T) correction evaluated with a large basis set.

155
Again, it can be seen that the HF+WL separations are underestimated compared with the best estimate data.

156
The interaction energies, although overestimated for the T and SP structures, still show a reasonable trend compared with the F structure, although the T and SP structures are not isoenergetic as suggested in the results of the ab initio calculations.

157
The F structure is bound in the HF+WL model, but underestimates the attraction, again emphasising that the HF+WL method cannot compete with post-HF approaches.

158
Also shown in Table 4 for comparison are results from full DFT calculations, as reported by Tsuzuki and Lüthi51 and also from current work, that confirm the unsuitability of these functionals as applied to weakly-bound systems.

159
The HF+WL approach provides a reasonable estimate of the energetic and structural trends of this system (while not being quantitatively accurate), where full DFT can not.

160
The similarity between the PBE and PW91 results is also noted.

161
Table 4 shows the interaction energy of the T, F and SP structures calculated at the X3LYP/aug-cc-pVDZ and xPBE/aug-cc-pVDZ levels of theory, at the geometries taken from Tsuzuki et al.,50 revealing a similar performance of these two functionals compared with existing functionals.

162
The interaction energy as a function of centre-of-mass separation was also investigated for these two functionals.

163
The F structure was unbound at all separations up to and including 6.0 Å at both the xPBE/aug-cc-pVDZ and X3LYP/aug-cc-pVDZ levels.

164
The T structure was bound at a separation of 5.7 Å at the xPBE/aug-cc-pVDZ level, with an optimum interaction energy of −0.48 kcal mol−1.

165
Similarly, at the X3LYP/aug-cc-pVDZ level, the T structure was bound at an optimum distance of 5.6 Å with an interaction energy of −0.37 kcal mol−1.

166
The SP structure was investigated in the neighbourhood of the structure of Tsuzuki et al., such that the vertical and horizontal separations, as well as the distance along the vector connecting the two centre of mass points, was varied by displacements of ±0.2 Å.

167
None of these locally-perturbed structures showed binding at the xPBE/aug-cc-pVDZ level.

168
The performance was very similar at the X3LYP/aug-cc-pVDZ level.

169
The only structure that yielded an interaction energy of less than +1 (but greater than zero) kcal mol−1, for either xPBE or X3LYP, was for the case when the distance along the vector connecting the two centre of mass points was extended by 0.2 Å.

170
These results, as well as those summarised in Table 4, demonstrate that the recovery of binding for molecular complexes, even at separations where the fragments overlap, remains an elusive goal for mainstream DFT exchange–correlation functionals.

171
The performance of all of the functionals examined here is poor in comparison with the HF+WL approach.

172
The basis set dependence of the the HF+WL approach was also gauged for the benzene dimer structures.

173
At the minimum-energy separations found at the HF+WL/aug-cc-pVDZ level, the interaction energies were recalculated using the 6-31G* basis set.

174
While the T and SP structures remained bound (interaction energies of −2.83 kcal mol−1 and −1.85 kcal mol−1, respectively), the F structure is unbound at this level of theory, at +0.78 kcal mol−1.

175
Exploration of the interaction energy of the F structure at a range of separations confirmed that this structure remained unbound.

176
Increasing flexibility by using a triple split valence basis, in this case using 6-311G*, leads to binding of the F structure (−0.53 kcal mol−1).

177
At the HF+WL/6-311G* level, the T and SP structures also remained with reasonable interaction energies (−3.16 and −2.84 kcal mol−1, respectively).

178
Similarly, diffuse functions were added to the 6-31G* basis, such that the 6-31+G* basis set was considered.

179
At this level, all structures are bound, with F, T and SP structures having interaction energies −0.20, −3.08 and −2.58 kcal mol−1, respectively.

180
These data suggest that a basis with at least split valence plus diffuse and polarization functions, or triple split valence plus polarization should be used to obtain qualitatively correct HF+WL interaction energies.

181
Finally, by way of comparison with other empirically-corrected methods, Wu et al6. report a DFT + disp interaction energy for the SP structure (at the Tsuzuki et al50. geometry) of −3.99 kcal mol−1.

182
Gonzalez and Lim24 report a HF + disp interaction energy of the T structure (at the geometry used by Tsuzuki et al50). of −2.80 kcal mol−1, albeit with a rather modest basis for the HF component.

183
The naphthalene dimer presents a significant challenge in the calculation of competitive low-energy structures.

184
Calculations at the MP2 level tend to overestimate the attraction in the naphthalene dimer, most significantly for the structures lowest in energy.52

185
Application of DFT to this system is also particularly unsuitable: the high-energy structures have attractive electrostatic energy contributions and relatively small dispersion contributions, and the low-energy structures have electrostatic contributions that are less favourable, while having relatively large dispersion contributions.44,52

186
This is opposed to the benzene dimer where the low-energy structures have favourable electrostatic interactions.

187
Accordingly, with the exception of LDA, most functionals will incorrectly predict that the T structures (T1 and T2 as shown in ref. 44) are the lowest in energy (no isomers are bound using PW91 and B3LYP).53

188
In other words, for the naphthalene dimer, fortuitous location of the low-energy structures, by virtue of attractive electrostatic interactions, is not possible, making this system an excellent test case for new functionals.

189
Most published theoretical studies of the naphthalene dimer have made use of MP2 calculations.44,52,54–56

190
Full CCSD(T) calculations with any reasonable basis set are currently out of reach, at least for the low symmetry structures.

191
Several studies have reported low-energy structures of the naphthalene dimer using MP2 calculations with a range of basis sets,44,52,54,55 including correlation corrections to MP2 interaction energies, using CCSD(T) calculations.52

192
It appears that, at present, the lowest-energy structure corresponds to a crossed C2 geometry, with a slipped Ci structure also energetically competitive.44

193
This result has also been reported by Gonzalez and Lim, using a HF+D approach.24

194
Here, four low-energy structures (Fig. 3) and one high-energy structure are considered.

195
The optimum interaction energies and centre-of-mass separations are summarised in Table 5, with the MP2 (CCSD(T)-corrected) data of Tsuzuki et al52. presented for comparison, along with the HF+D data of Gonzalez and Lim.24

196
For each low-energy structure, the separations were varied with displacements of ±0.2 Å away (along the vector connecting the centre ofmass points) from the optimum separation (as was determined by MP2/6-31G* geometry optimizations.44)

197
For the high-energy D2h structure, the optimum separation was taken from Tsuzuki et al.52

198
This separation was varied along the vector connecting the centre of mass points in the same way.

199
It is seen in Table 5 that the HF+WL separations all seem reasonable when compared with the ab initio data, although are consistently too short.

200
The trend in energetics for the D2d, Ci and C2h structures agrees with the best estimates, although there is not quantitative agreement.

201
In particular, the best estimates indicate that the D2d and C2h structures are almost isoenergetic, while the HF+WL data support a difference of over 0.5 kcal mol−1.

202
Tsuzuki et al. did not examine the C2 structure, so we have made comparisons with the trend from the HF+D results of Gonzalez and Lim.24

203
In agreement with their findings, and with previous MP2 optimisations,44 the HF+WL approach also supports the C2 structure as the putative global minimum for the naphthalene dimer system.

204
The absolute HF+WL interaction energy for the D2h structure deviates significantly from the ab initio value, although it is bound, and exhibits the smallest interaction energy of the set, in accord with Tsuzuki et al. Again, these data emphasise that HF+WL cannot outperform high-quality post-HF calculations.

205
The X3LYP and xPBE functionals were also tested for all five naphthalene dimer structures.

206
For both of these functionals, none of the five structures were bound, out to a centre-of-mass separation of 6.0 Å.

207
For the xPBE functional, at a 6.0 Å separation, interaction energies were 0.19, 0.12, 0.72, 0.06 and 0.23 kcal mol−1 for the D2d, C2, C2h, Ci and D2h structures, respectively.

208
The X3LYP values showed a similar trend, with absolute interaction energies slightly higher in all cases.

209
This behaviour is similar to that noted for other popular functionals.53

210
It is clear that neither of these functionals can reliably yield binding for aromatic ring interactions at a qualitative level.

211
The basis set dependence of the HF+WL approach was also tested for the naphthalene dimer.

212
At first the 6-31G* basis was used.

213
This appeared to give satisfactory results for the four low-energy structures–yielding −4.77, −3.10, −4.70 and −3.93 kcal mol−1 for the C2, D2d, Ci and C2h structures, respectively.

214
At this level of theory, the D2h was not bound for any separation.

215
This result confirms the behaviour of the HF+WL/6-31G* level as noted for the benzene dimer calculations.

216
The D2h was checked instead using both the 6-31+G* and 6-311G* basis sets, with both yielding binding, optimised at −1.07 kcal mol−1 and −1.80 kcal mol−1, respectively.

217
The HF+WL/6-31+G* also yields reasonable results for the low-energy structures, giving interaction energies of −5.83, −4.14, −5.75 and −4.86 kcal mol−1 for the C2, D2d, Ci and C2h structures respectively.

218
This confirms the necessity for using at least split-valence plus diffuse functions in the basis set when using HF+WL.

219
The final test system comprises the base-stacking interaction for dimers of pyrimidine, cytosine, uracil and guanine.

220
The stabilisation of stacked base-pairs is thought to arise from van der Waals interactions.57

221
The geometries of these systems are shown in Fig. 4.

222
These systems have been the subject of studies using MP2,45 CCSD(T),45 empirically-corrected DFT16,18,58 and HF+D calculations.58

223
The pyrimidine system has been best characterised, at the CCSD(T)/aug-cc-pVDZ(+f functions and bond functions) level of theory,45 but only at one geometry (see Fig. 4) with the vertical separation fixed to 3.3 Å.

224
For the remaining systems (fixed at the geometries in Fig. 4) examined by Hobza and Šponer, uracil was characterised up to the MP2/aug-cc-pVTZ(+ bond functions) level, cytosine up to the MP2/aug-cc-pVTZ leveland guanine up to the MP2/aug-cc-pVDZ level.

225
In this by work of Hobza and Šponer,45 the vertical separation of cytosine and guanine dimers were fixed at 3.3 Å and 3.4 Å, respectively, and the uracil dimer geometry came from geometry optimisation at the MP2/6-31G** level.

226
The minimum energy separations and interaction energies for these four systems, calculated at the HF+WL/aug-cc-pVDZ level, are reported in Table 6, as are ab initio data obtained by Hobza and Šponer.45

227
The ab initio data of Hobza and Šponer are the best estimates (obtained using an extrapolation procedure involving MP2 complete basis set data with CCSD(T) corrections) available at present.

228
The agreement with the best energetic estimates from Hobza and Šponer is good for all cases.

229
The trend in the interaction energies is reproduced, although the absolute values appear fortuitously close to the best estimates.

230
It is interesting to note that the lowest-energy separation for pyrimidine is larger than the separation estimated by Hobza and Šponer.

231
All other separations obtained using the HF+WL approach are, as expected, underestimated.

232
These systems were also investigated using the xPBE and X3LYP functionals, both with the aug-cc-pVDZ basis set.

233
The results are summarised in Table 6.

234
Every structure was bound for both functionals, although the attraction is decreased in comparison with the data of Hobza and Šponer.

235
Furthermore, the separation corresponding to the optimum interaction energy is increased relative to the results from Hobza and Šponer. xPBE manages to capture the trend of the interaction energies, although as already observed for the methane dimer, the range of interaction energies is considerably smaller than is seen in the ab initio results.

236
X3LYP is seen to perform similarly to xPBE for these systems.

237
The basis set dependency was again explored by repeating these calculations at the HF+WL/6-31G* level.

238
Again, the vertical separation between monomers was optimised to yield the optimum interaction energy in each case.

239
The results with the 6-31G* basis do not show much difference with the data obtained with the aug-cc-pVDZ basis.

240
The interaction energies of pyrimidine, uracil, cytosine and guanine were −3.1, −11.3, −10.3 and −14.0 kcal mol−1, respectively.

241
The corresponding optimised separations at the HF+WL/6-31G* level were 3.4, 2.97, 3.1 and 3.0 Å, respectively.

242
This behaviour in stark contrast with the behaviour noted for the D2h isomer of the naphthalene dimer, and the F structure of the benzene dimer, both of which were unbound at the HF+WL/6-31G* level.

243
It must be re-emphasised that, as a general rule, it appears that HF+WL calculations must be used with at least split-valence plus diffuse functions in order to capture qualitative behaviour.

244
It is clear from all of the test cases examined that the HF+WL approach is capable of recovering qualitative trends in interaction energies for pure van der Waals (dispersion and repulsion) systems, albeit with shorter equilibrium distances (most of the time).

245
Relative orientations seemed to be captured well for the HF+WL approach.

246
However, no claim is made to suggest that the HF+WL approach can describe a pure van der Waals system with quantitative accuracy, and therefore is not competitive with post-Hartree–Fock methodologies.

247
To quantify the conditions under which this approach appears to be of value it is noted that, for example, in the case of the rare gas dimers, the HF+WL approach yields reasonable results for systems with overlaps no smaller than around 10−5 (with the exception of the helium dimer).

248
A HF+WL-type of approach may be extremely useful for optimising geometries of low-symmetry structures.

249
Low-symmetry structures can be very time-consuming to optimise at any post-HF level (with a sufficiently large basis set), although advances are being made in this area (e.g. see .ref. 4)

250
The results of HF+WL optimisations would at least indicate which of the structures should be considered for further refinement by more expensive post-HF methods.

251
For example, the location of the C2 structure of the naphthalene dimer was missed by several MP2-based geometry optimisation studies, presumably because of the low symmetry and the computational demands this creates.

252
Another example is the geometry optimisation of small molecules adsorbed onto periodic substrates.

253
Periodic EXX calculations are now possible via some mainstream plane-wave codes, making EXXWL presumably more efficient than currently available periodic MP2 implementations (as found in GAUSSIAN03).

254
It is also noted that the HF+WL approach appeared to work reliably only when a basis set featuring at least split-valence plus diffuse functions was used.

255
The HF+WL/6-31G* level was demonstrated to be insufficient in several instances.

256
Furthermore, these conclusions only hold at separations where the charge distribution on each fragment overlaps.

257
The HF+WL approach cannot be successfully applied in the large separation regime.

258
It is also clear that the current group of exchange–correlation functionals, B3LYP, PBE, X3LYP and xPBE, cannot be reliably used even for qualitative work.

259
For these functionals, of those structures that showed any binding at all, the binding was substantially underestimated compared with best estimates from ab initio calculations.

260
Two clear limitations of these current calculations are acknowledged.

261
The first is that for all the systems studied here, the monomers were treated as rigid body fragments.

262
However, the intramolecular degrees of freedom could in principle be optimised using the HF+WL approach.

263
The second limitation is on the type of system that can be studied with this method (the same restrictions apply to HF+D methods), and is a consequence of how the correlation contribution forms a perturbative correction to the total Hartree–Fock energy.

264
This limitation could be lifted if the WL correlation could be included dynamically, within Kohn–Sham DFT, via the implementation of exact-exchange methods (denoted here as the proposed EXXWL functional).

265
Methods based on exact exchange, such as the optimised effective potential (OEP),35 the KLI approximation,59 and more recently, the localised HF approach60 and the direct optimisation method,61 will also ensure that the intermolecular exchange behaves as HF exchange.

266
The fully self-consistent EXXWL functional would allow treatment of systems where a general HF+D approach currently fails (such as the NaAr complex, molecular physisorption on metal surfaces, etc)., where the HF wavefunction cannot distort if dynamic correlation effects are not accounted for.

267
Work is currently underway to enable calculations with the EXXWL functional in future.

268
Compared with the HF+WL approach, a much wider range of systems could be treated with the EXXWL functional.

269
Furthermore, promising developments of inclusion of ‘left–right’ correlation with exact exchange within the KS-DFT framework62 suggest that an approach that combines non-dynamical correlation with the EXXWL functional may in future permit reasonable modelling of bond-breaking in weakly-bonded systems.

270
Such an EXXWL functional is not intended to act as a replacement for more accurate post-HF methods, but rather as a device for optimising the structure and energetics of weakly-bonded systems that are not amenable even to MP2 optimisations.

271
Compelling reasons for why the WL functional works so well in the cases studied here remain unclear.

272
Of course, there is no way that the WL correlation functional can provide a realistic description of London dispersion interactions at long range.

273
WL was parametrised from high-quality quantum chemical data for a small number of atoms.

274
It also, by design, should work well with Hartree–Fock densities.

275
Wilson and Levy point out that, due to the way in which the functional was parametrised, it may not perform well for systems where the Hartree–Fock density differs substantially from the true density.

276
Because of the constraints imposed by the HF+D-like approach used here, the systems chosen for study should have HF densities that are reasonably close to the true density.

277
Therefore, an investigation of systems where there is less similarity between HF and real densities is, at present, desirable.

278
However, these systems should only be studied within a fully self-consistent implementation of the EXXWL functional.

279
The WL correlation potential has been examined, and appears to have unremarkable properties in comparison with other correlation potentials.

280
It is similar to other correlation potentials, e.g. the Perdew(86)63 correlation functional, in the sense that the potential decays to the asymptotic limit exponentially.

281
This behaviour is reasonable in the sense that the correlation potential should decay faster than the exchange potential (in the case of the HF potential, the asymptotic decay goes like −1/r).

282
However, it has been remarked by others64 that perhaps a more realistic asymptotic behaviour of the correlation potential would be −α/r4, where α is the polarizability.

283
The WL correlation functional is based on the local Wigner correlation functional30 (W38).

284
The LYP correlation functional,65 due to its basis on the Colle–Salvetti formula, is also based on the form of the W38 functional.

285
The LYP functional comprises a (local) W38 component multiplied by another functional that contains gradient correction terms.

286
This is in common with many exchange functionals that analogously comprise the Dirac (local) exchange, multiplied by another functional that contains gradient correction terms.

287
Functionals of a form such as this lend themselves easily to an analysis in the low-density/high-density-gradient regime, since one can construct an ‘enhancement factor’ (the function that is multiplied by the local component of the functional).

288
Unfortunately, the WL functional is not multiplicatively separable in this sense from the W38 functional.

289
Since it is difficult to disentangle the W38 functional from the WL functional, the WL functional is therefore not amenable to an analogous analysis of the ‘enhancement factor’.

290
It is worth mentioning that even W38 (without gradient corrections) in partnership with Becke8834 exchange (the B-W functional)66 has been shown to perform comparably with BLYP, although the reasons for this also remain unclear.

291
Finally, it is plausible that the WL functional yields unphysical attractive contributions for some of the short-ranged contributions to the interaction energy.

292
A decomposition scheme based on Heitler–London type calculations, as reported by Cybulski and Seversen,13 may be helpful in this instance.

293
The fact that the HF+WL almost always underestimates monomer separations seems consistent with this idea.

294
Even if the origins of this unphysical behaviour can be revealed, further work is required to understand why these (potentially) unphysical terms combine such that very good interaction energies are ultimately recovered, albeit fortuitously.

295
Thus, while admittedly not having a satisfying and physically rigorous theoretical justification, the HF+WL approach (or in future, the EXXWL functional) may become a useful pragmatic device for optimising the geometries of large and/or complicated van der Waals complexes that are currently beyond the reach of optimisation by more sophisticated ab initio methods.

296
Other functionals, including those that are claimed to describe weak interactions such as xPBE and X3LYP, while also not capable of recovering London dispersion, cannot reliably perform the role of a pragmatic device, even when monomer charge distributions overlap.

Conclusions

297
The Wilson–Levy (WL) correlation functional was used as a post-Hartree–Fock (HF) perturbative correction in the evaluation of the interaction energy for a range of weakly-bound dimers, in the regime of small to intermediate fragment separations.

298
It was found that this HF+WL approach works well in a pragmatic and qualitative sense for the energetics of closed-shell van der Waals systems, but without a satisfying physical justification, and without quantitative accuracy.

299
Two newly-developed exchange–correlation functionals, designed for applications to biological systems (X3LYP and xPBE), were also investigated and found to perform not as well as the HF+WL approach.

300
The results of this work suggest an apparently promising device in the form of the EXXWL functional, for determining the energetic trends of molecules physisorbed on surfaces, nanotubes and other periodic structures.