Water on BN doped benzene: A hard test for exchange-correlation functionals and the impact of exact exchange on weak binding

Density functional theory (DFT) studies of weakly interacting complexes have recently focused on the importance of van der Waals dispersion forces whereas, the role of exchange has received far less attention. Here, by exploiting the subtle binding between water and a boron and nitrogen doped benzene derivative (1,2-azaborine) we show how exact exchange can alter the binding conformation within a complex. Benchmark values have been calculated for three orientations of the water monomer on 1,2-azaborine from explicitly correlated quantum chemical methods, and we have also used diffusion quantum Monte Carlo. For a host of popular DFT exchange-correlation functionals we show that the lack of exact exchange leads to the wrong lowest energy orientation of water on 1,2-azaborine. As such, we suggest that a high proportion of exact exchange and the associated improvement in the electronic structure could be needed for the accurate prediction of physisorption sites on doped surfaces and in complex organic molecules. Meanwhile to predict correct absolute interaction energies an accurate description of exchange needs to be augmented by dispersion inclusive functionals, and certain non-local van der Waals functionals (optB88- and optB86b-vdW) perform very well for absolute interaction energies. Through a comparison with water on benzene and borazine (B$_3$N$_3$H$_6$) we show that these results could have implications for the interaction of water with doped graphene surfaces, and suggest a possible way of tuning the interaction energy.

Density functional theory (DFT) studies of weakly interacting complexes have recently focused on the importance of van der Waals dispersion forces whereas, the role of exchange has received far less attention.Here, by exploiting the subtle binding between water and a boron and nitrogen doped benzene derivative (1,2-azaborine)   we show how exact exchange can alter the binding conformation within a complex.
Benchmark values have been calculated for three orientations of the water monomer on 1,2-azaborine from explicitly correlated quantum chemical methods, and we have also used diffusion quantum Monte Carlo.For a host of popular DFT exchangecorrelation functionals we show that the lack of exact exchange leads to the wrong lowest energy orientation of water on 1,2-azaborine.As such, we suggest that a high proportion of exact exchange and the associated improvement in the electronic structure could be needed for the accurate prediction of physisorption sites on doped surfaces and in complex organic molecules.Meanwhile to predict correct absolute interaction energies an accurate description of exchange needs to be augmented by dispersion inclusive functionals, and certain non-local van der Waals functionals (optB88-and optB86b-vdW) perform very well for absolute interaction energies.Through a comparison with water on benzene and borazine (B 3 N 3 H 6 ) we show that these results could have implications for the interaction of water with doped graphene surfaces, and suggest a possible way of tuning the interaction energy.a) Electronic mail: angelos.michaelides@ucl.ac.uk

I. INTRODUCTION
An accurate description of the structures and energies of weakly interacting systems is important in materials science and biology, but it is often difficult to obtain reference data either theoretically or experimentally.A key challenge lies in the ability to capture small energy differences -on the order of a few meV -that can have drastic effects on the structure.2][3][4][5][6] Likewise for water clusters, most notably the water hexamer, there are several isomers that have energies within just 5 meV/H 2 O. 5,7,8 Furthermore, in biological applications there can be numerous shallow energy minima with different conformations.Particularly for complex organic systems, predicting the exact lowest energy conformation is crucial to determine the crystal structure of drugs 9,10 and the mechanisms by which proteins function. 11e particularly interesting weakly interacting system of relevance to potential applications in clean energy, water purification, hydrogen storage, and bio-sensing, [12][13][14][15][16] is the interaction of water with layered materials.Notably interfacial water on graphene, hexagonal boron nitride (h-BN) and hybrid composites of these materials.[29] One can use smaller model systems for graphene 25,[29][30][31][32] and h-BN, such as benzene and the inorganic counterpart borazine (B 3 N 3 H 6 ), to help understand the interaction with water.
With these small molecules, it is possible to use high accuracy methods to calculate benchmark interaction energies and binding conformations, [33][34][35][36][37] that would otherwise be infeasible for the extended surfaces.Given the shortage of reference data across these systems, there is a strong incentive to deliver accurate benchmark calculations, and for our purposes, we require a model system that is a hybrid of benzene and borazine.With that in mind, we study the weak binding between water and an aromatic molecule known as 1,2-dihydro-1,2azaborine (or 1,2-azaborine for short) as a reference system.The 1,2-azaborine molecule shares many similarities to benzene, with the clear distinction being the asymmetry of the molecule due to boron and nitrogen substitution (see Fig. 1).The mixture of boron, nitrogen and carbon atoms in this molecule makes it a suitable model for benchmarking in relation to extended surfaces that are hybrids of graphene and h-BN. 38Moreover, the asymmetry makes 1,2-azaborine an ideal system for testing the performance of computational methods because each atom in the ring has a distinct chemical environment that serves as a tag, in contrast to benzene in which carbon atoms are obviously indistinguishable.
5][46] Here, we have examined the water/1,2-azaborine system with coupled cluster with single, double and perturbative triple excitations (CCSD(T)), and a range of density functional theory 47,48 (DFT) exchange-correlation (xc) functionals.0][51][52] As part of this study we have also tested diffusion quantum Monte Carlo (DMC).For selected systems DMC has been shown to produce very accurate results 27,36,53,54 including "subchemical" accuracy (< 1 kcal mol −1 ) in dispersion dominated systems. 55Here we show that DMC does well with respect to the CCSD(T) benchmarks, again achieving subchemical accuracy.This short article will start with a brief summary of the employed methods, followed by results from a set of calculations that allow us to directly compare the performance of different xc functionals with other explicitly correlated methods.After analysis and discussion of the various DFT results on the intermediate 1,2-azaborine system, we investigate the effect of the boron nitride doping, by also studying the interaction of water with the pure systems of benzene and borazine.We close with a discussion and some general conclusions.

A. Computational Methods
7][58] Second order Møller-Plesset perturbation theory (MP2) with up to aug-cc-pV5Z basis sets 59 along with CCSD(T) calculations at the aug-cc-pVTZ level have been conducted.Due to the unfavorable scaling of CCSD(T), it is more feasible to conduct MP2 calculations with larger basis sets, deducing the complete basis set (CBS) limit, and subsequently calculating the ∆CCSD(T) value of absolute interaction energy at the CBS limit.For a description of this procedure along with analysis of errors, the reader is referred to the recent work of Sherrill and coworkers. 60garding the CBS limit, various extrapolation schemes have been discussed [61][62][63][64][65] and we have chosen to use the one proposed by Halkier et al. [62][63][64] Gaussian03 66 was used for the Hartree-Fock (HF) and post-HF calculations.
The initial single particle wavefunctions for use in DMC were obtained from DFT planewave (PW) calculations using the PWSCF package 67 and Trail and Needs pseudopotentials (PPs) were used for all atoms in the system, 68,69 warranting a standard 300 Ry energy cut-off.
Previous work by Ma et al. 36 indicates that weak binding energies are not overly sensitive to the trial wavefunctions (TWs), having tested a few xc functionals (including hybrids) and also HF.We have generated TWs using the local density approximation (LDA) 70 and also the Perdew-Burke-Ernzerhof (PBE) 71 xc functional.The resulting wavefunctions were expanded in terms of B-splines 72 for efficiency.DMC calculations have been performed using the CASINO code, 73 and we have used Slater-Jastrow type TWs, in which the Jastrow factor contains electron-nucleus, electron-electron, and electron-electron-nucleus terms.We used a combination of DMC calculations using 16000 walkers across 160 cores and 64000 walkers across 640 cores.Final DMC results have been derived by the weighted averaging of the results and errors.Time steps of 0.0025, 0.005 and 0.01 a.u.have been tested and the locality approximation was utilized. 74We obtained statistical error bars for our interaction energies of ±3 meV, which corresponds to 1σ.6][77][78] was used for all the DFT calculations.VASP employs plane-wave basis sets and uses projector augmented wave (PAW) potentials 79,80 to model the core region of atoms.After a series of convergence tests for the plane-wave cut-off energy and unit cell, we chose to use a 500 eV cut-off energy and a 15 Å length cubic unit cell, along with Γ-point sampling of reciprocal space.
There is, of course, an almost endless list of xc functionals that could be considered, and here we benchmark a selection of fairly widely used functionals.The functionals tested include PBE 71 which is a generalized gradient approximation (GGA) functional that does not contain long-range correlation.4][85][86] There have been many developments to include van der Waals (vdW) dispersion in xc functionals, as discussed in the perspective of Klimeš and Michaelides 49 and reviewed by Grimme, 87 and here we have tested several of these vdW-inclusive DFT approaches.Specifically, PBE-D2, 88 a semiempirical functional that contains Grimme's D2 correction, and also two correction schemes from Tkatchenko and Scheffler, namely vdW-TS 89 and vdW-TS+SCS, 90 referred to here as TS and TS+SCS, respectively.Using the TS and TS+SCS schemes, C 6 coefficients and vdW radii are determined from ground state electron densities, 89 whilst TS+SCS also includes long-range screening effects. 90The TS and TS+SCS corrections will be applied to PBE, PBE0 and B3LYP 91 .][94][95][96][97][98][99][100][101][102][103] The vdW-DFs considered include the original vdW-DF which we refer to as revPBE-vdW, 92 several optimized vdW functionals (optPBE-vdW, 93 optB88-vdW, 93 and optB86b-vdW 94 ), and also vdW-DF2. 95The exchange part of revPBE-vdW and the optimized vdW functionals are different but they all share the same non-local correlation part.In vdW-DF2, both the exchange and correlation components have been modified.

B. Water and 1,2-azaborine Setup
The absolute interaction energy, E int , that we refer to throughout this study is defined as, where E tot com is the total energy of the bound complex between water and the substrate (1,2azaborine benzene or borazine), and E tot sub and E tot wat are the total energies of relaxed substrate molecule and water, respectively.
In 1,2-azaborine the electronic environment of individual carbon atoms differs due to the asymmetry introduced through the boron and nitrogen atoms.We use the numbering scheme shown in Fig. 1(a) to make the distinction between the atoms.
We obtained three distinct orientations of water over 1,2-azaborine using PBE-D2 104 in order to carry out the benchmarking study, namely C3, C5 and C5C3, also depicted in Fig. 1.The nomenclature of the complexes refers to the specific carbon atoms involved in hydrogen bonding, keeping in line with the numbering scheme in Fig. 1(a).The hydrogen atoms of water point toward carbon-3 and carbon-5 in all three complexes; most likely as a result of the higher electron charge around these carbon atoms due to the conjugation of the localized lone pair of electrons from the nitrogen atom.Indeed, Bader analysis shows that carbon-3 and carbon-5 have a larger atomic volume compared to carbon-2 and carbon-4, as a result of increased electron charge around them.Some characteristic structural parameters are listed in Table I for the three different configurations.Although the PBE-D2 geometries are not benchmark accuracy, there have been various studies on polymer crystals and layered materials, [105][106][107] and even for water on graphene, 28 indicating that PBE-D2 can provide reasonable structures for weakly interacting systems.Interested readers may also refer to the Supporting Information 108 (SI) for bond lengths as well as the structures of fully relaxed complexes with MP2 and several of the xc functionals that we test.We computed benchmark absolute interaction energies for the three water adsorption complexes (C3, C5, and C5C3) mentioned in Section II B, using ∆CCSD(T).The results for ∆CCSD(T) have been extrapolated to the CBS limit and the computed interaction energies reveal that C3 is the most stable binding configuration with an interaction energy of −155 meV, followed by C5 (−146 meV) and C5C3 (−143 meV), as listed in Table II.
The interaction energies we computed with DMC are in reasonable agreement with ∆CCSD(T); within 5-14 meV depending on the adsorption structure considered and the TW used.This level of agreement is in line with several recent studies in which DMC has been compared to coupled cluster. 27,36,55,109The small difference between DMC and ∆CCSD(T) interaction energies could be due to issues such as the use of PPs in DMC, the fixed node approximation in DMC, or the approximations used to obtain the ∆CCSD(T) values (including CBS extrapolation). 60With regard to the relative stabilities of the complexes, the DMC results suggest the same trend as ∆CCSD(T): C3 is more stable than C5 and C5 is more stable than C5C3.Considering the statistical error bars for each DMC interaction energy however, C3 and C5 could be degenerate according to DMC with LDA TWs.Whilst with PBE TWs, C3 is 9 meV more stable than C5, indicating a better trend prediction with the latter.Note the total energies using LDA TWs are slightly lower than those obtained with PBE TWs 110 and since DMC is a variational method, we consider the interaction energies with LDA TWs to be slightly more reliable in this particular system.
Therefore, once again, it appears that DMC is useful for obtaining reliable interaction energies, but there is also an inherent difficulty in using a stochastic method like DMC, to clearly distinguish between complexes with very small energy differences.The TS and TS+SCS corrections perform significantly better than the D2 correction, with only 10% and 5% errors, respectively.
One can see from Fig. 2 that the percentage error lines across C3, C5, and C5C3 have a very similar form, but they shift with regard to the reference ∆CCSD(T) binding energy.This means that for the C5 complex vdW-DF2 is providing the best agreement with ∆CCSD(T) (2% error) and the optimized vdW functionals are overbinding by 10-15%.
Whereas for C5C3, the xc functionals perform in a similar manner as for C3, with the optimized vdW functionals performing the best once again (< 5% error).Of PBE and its dispersion corrected forms, PBE-TS+SCS performs the best for all three complexes (under-binding by 5-15%).Note that MP2 consistently overbinds all three structures by ∼ 9%.
Regardless of the absolute interaction energies of PBE, dispersion corrected PBE and the non-local vdW functionals, they all fail to predict C3 as the most stable complex.In addition, the TS and TS+SCS corrections are not satisfactory as they stabilize the C5C3 complex such that it becomes degenerate with either C3 or C5.Clearly with a fairly flat potential energy surface, the difference between C3 and C5 is a considerable challenge for the xc functionals.According to the benchmark ∆CCSD(T) values, C5 is only 9 meV less stable than C3, making it difficult to assign the source of error that leads to so many different xc functionals predicting the wrong trend.
One possible source of error is an inadequate description of exchange and to address this we initially performed HF calculations.We find that with HF the trend is correctly predicted with C3 as the most stable configuration, despite the lack of correlation and highly underestimated interaction energies.The HF results suggest that the lack of exact exchange is perhaps the main reason for many of the xc functionals predicting C5 instead of C3.
For further insight, HF symmetry adapted perturbation theory (HF-SAPT) calculations 111 revealed that it is the exchange-repulsion energy (mostly electrostatic) that puts the binding energy in favor of C3 and not C5.
Guided by this insight, we computed binding energies using the hybrid PBE0 and B3LYP functionals with varying amounts of exact exchange.The results are listed in Table II with the proportion of exact exchange indicated by the superscripts, along with hybrid functional results corrected with the TS and TS+SCS vdW schemes.In addition, Fig. 3 shows how the interaction energies of the C3 and C5 complexes vary with the proportion of exact exchange with PBE0-like functionals.Using standard PBE0 0.25 the lowest energy complex is C5, but as the percentage of exact exchange is increased to 50% and above, C3 becomes the most stable complex.Similarly with B3LYP, 60% exact exchange is needed in order to switch the site preference to C3. However PBE0 0.75 and B3LYP 0.60 still underbind by 20% with respect to ∆CCSD(T), and including the TS+SCS corrections leads to overbinding by the same amount.The combination of the hybrid functionals with TS and TS+SCS corrections tends to decrease the energy differences between the different configurations.Depending on the particular combination of exact exchange and TS-like dispersion, both C3 and C5 can be degenerate or almost degenerate or even all three structures can be almost degenerate.
Overall we learn from these various calculations that a high proportion of exact exchange improves the relative energies of the various structures.However, still, in the future a more refined description of long-range correlation and exchange is needed in order to predict the correct trend and absolute interaction energies.
A change in the stability trend of weakly interacting complexes due to the amount of exact exchange is not limited to the systems studied here; Thonhauser et al. 44 observed a similar distortion in the ordering of conformers in their study of monosubstituted benzene dimers, which they corrected by using HF exchange.The need for a very high fraction of exact exchange in xc functionals to give correct predictions warrants further discussion.We have very carefully looked at the electronic structures obtained from the various functionals for the different complexes and a thorough inspection of the individual energy contributions to the absolute interaction energies was particularly informative.By decomposing the interaction energies into the kinetic, potential, Hartree, exchange and correlation energies, we have analyzed the effect of increased exact exchange on the individual energy contributions for C3 and C5.We find that the main distinguishing feature between C3 and C5, is that the kinetic energy and the repulsive terms (Hartree and exchange energies) change to different extents as the fraction of exact exchange is increased.Kinetic energy favors C5 but with a high fraction of exact exchange, the C3 complex is stabilized by less repulsion (Hartree and exchange) than in C5.This observation is in accord with the HF-SAPT results mentioned earlier.
It is known that exact exchange localizes the electrons and alleviates the surplus delocalization which is otherwise present in standard GGA calculations.The localization of electrons determines the electron density distribution over a system and therefore determines the interaction sites within complexes.In addition to orbital overlap, hydrogen bonding and dispersion interactions are also both affected by a change in the electron density distribution.The interactions in our water/1,2-azaborine complexes have contributions from hydrogen bonding, dispersion interactions and weak orbital mixing.We have seen that the delicate balance between such interactions is sensitive to the relative contributions from individual energy terms, prompting the need for both exact exchange and a good description of the correlation energy for the reliable prediction of configurations in weakly interacting systems, such as ours.

B. From benzene to borazine
Thus far we have studied the interaction of a water monomer with the intermediate 1,2-azaborine molecule as a model for boron nitride doped graphene.We have found that an inadequate description of exchange in xc functionals can alter the binding orientation, despite reasonable interaction energies being predicted by some of the vdW xc functionals.
Here we establish whether xc functionals without exact exchange can correctly describe the interaction between the water monomer and the pure counterparts of 1,2-azaborine: benzene and borazine (B 3 N 3 H 6 ).It is also useful, from a materials design perspective, to understand how the boron nitride doping affects the interaction with water, compared to graphene and h-BN; which can be mimicked to some extent by benzene and borazine, respectively.
Adhering to the use of model systems allows us to compare interaction energies from DFT and benchmark CCSD(T) calculations.To this end, we have examined water on benzene and borazine and we find that in contrast to 1,2-azaborine, these pure systems are much less challenging, with different functionals able to predict the same orientation as obtained in previous benchmark calculations. 33,34,37Previous work has shown that water tilts towards a carbon atom in benzene, and the complex has an absolute interaction energy of −145 meV calculated with CCSD(T). 33,34Similarly, in the water-borazine complex, water tilts towards a nitrogen atom and the CCSD(T) interaction energy reported by Wu et al. 37 is −92 meV.We have relaxed several starting geometries of the water-benzene and water-borazine complexes with PBE and optB86b-vdW and the absolute interaction energies are shown in Table III.For the water-benzene complex, PBE underbinds as previously shown, 35 whilst optB86b-vdW performs very well: only 2% error compared to the benchmark.Moreover with PBE, there is almost no distinction between water binding to benzene or 1,2-azaborine, but using the dispersion inclusive optB86b-vdW, the binding on 1,2-azaborine is almost 20 meV stronger.Similarly using CCSD(T), our benchmark absolute interaction energy of the C3 complex (−155 meV) is 10 meV stronger than the absolute interaction energy of water to benzene (−145 meV calculated by Min et al. 34 ).On borazine the absolute interaction energy is weaker according to all methods: PBE (−85 meV), optB86b-vdW (−122 meV) and also with CCSD(T) (−92 meV). 37In all cases, the water tilts towards nitrogen in the borazine ring 37 and this tilting, that also occurs on 1,2-azaborine and benzene, is indicative of a weak hydrogen bond forming in all three complexes.Not only do PBE and optB86b-vdW obtain the correct binding orientations for water on benzene and borazine, but also optB86b-vdW yields only a 2-3% error for the water-benzene and C3 complexes.For completeness, we also include the results obtained with vdW-DF2 in Table III, which are similar to those obtained with optB86b-vdW.Using vdW-DF2 the error for water-benzene is 5% and for the C3 complex the binding energy agrees with the benchmark.Even so, the wrong orientation is predicted by PBE, optB86b-vdW and vdW-DF2 for water binding to 1,2-azaborine as discussed in Section III A. Therefore, the prediction of binding sites with weak interactions on doped surfaces and between complex organic molecules could be compromised by an inadequate description of exchange, whilst being correct in pure systems.
Finally, it is interesting to note that water interacts more strongly with the intermediate 1,2-azaborine molecule than either benzene or borazine.If upon doping graphene with boron and nitrogen a similar increase in the interaction with water was found then this could provide a means of tuning the strength of the water substrate interaction to have specific, targeted wetting properties.

IV. DISCUSSION AND CONCLUSIONS
We have calculated three benchmark values of the absolute interaction energy between a water monomer and 1,2-azaborine using ∆CCSD(T) extrapolated to the CBS limit.
The lowest energy complex according to the explicitly correlated, exact exchange methods (∆CCSD(T), DMC and MP2) is C3 with an absolute interaction energy of −155 meV from ∆CCSD(T).The DMC energies are, on average, within 6-8% of ∆CCSD(T).Meanwhile, xc functionals including PBE, dispersion corrected PBE and an array of vdW-DFs, fail to predict the same binding configuration as ∆CCSD(T), MP2 or DMC.Instead, these xc functionals indicate that C5 is the lowest energy complex.Although interatomic many-body dispersion (MBD) forces have previously been shown to be neglected or wrongly described in DFT xc functionals, 112 HF and HF-SAPT calculations demonstrate that in this case, exact exchange corrects the trend without including MBD forces; thus we can deduce that MBD is not a source of error for the stability trend.Previous work has shown that for strongly interacting systems, like CO on Pt (111), exact exchange improves the alignment of electronic states between the substrate and the adsorbate and to recover the experimental chemisorption site. 113Also, for defects in semiconductors, exact exchange is a crucial factor in finding the correct defect states (see e.g.ref 114 ).Here, it is demonstrated that exact exchange and the associated changes to the electronic structure plays a decisive role in the prediction of the lowest energy configuration for weakly interacting (asymmetric) complexes.As such, it is imperative that more emphasis is placed on the accurate treatment of exchange when using DFT to examine weakly bound complexes and adsorption.The wrong prediction of the configuration for our relatively small system of water on 1,2-azaborine suggests that the delocalization error is likely to be even more pervasive for larger systems with more shallow energy minima, such as for physisorption on doped surfaces and in crystal polymorph prediction.
In terms of the absolute interaction energies, however, the inclusion of dispersion interactions is essential for weakly bound systems, and in the systems studied here this is more accurately achieved by non-local vdW functionals than dispersion corrected PBE, PBE0 and B3LYP.By comparing absolute interaction energies, we determined that optB88-vdW and optB86b-vdW are generally the best performing xc functionals from those tested.Our findings imply that in order to predict the correct binding configuration as well as the energy, DFT xc functionals must simultaneously contain exact exchange (so as to avoid delocalization error) and account for dispersion interactions.However, even with dispersion corrected hybrid functionals we were not able to obtain a perfect trend or accurate interaction energies; echoing the need to develop xc functionals that contain a better description as well as balance, of both exchange and non-local correlation. 51nally, we have found that water binds more strongly to the intermediate system, 1,2azaborine than to either benzene or borazine.Should equivalent behavior be observed on doped graphene surfaces, then this could represent a means of tuning the binding and wetting of interfacial water to graphene and h-BN.In future work we will focus on exploiting this possibility.

FIG. 2 .
FIG. 2. The percentage difference between ∆CCSD(T) interaction energy and that from DMC, MP2 and the various DFT xc functionals, of water to 1,2-azaborine.The solid black lines at zero represent the ∆CCSD(T) reference.Note that DMC results with both LDA (Ψ LDA ) and PBE (Ψ PBE ) TWs are reported and that the HF results are off the chart.The superscripts for the hybrid PBE0 and B3LYP functionals indicate different proportions of exact exchange.

FIG. 3 .
FIG.3.The interaction energy of the C3 and C5 complexes are plotted against the proportion of exact exchange in PBE0-style calculations, with zero exact exchange corresponding to PBE.
We are grateful for support from University College London and Argonne National Laboratory (ANL) through the Thomas Young Centre-ANL initiative.Some of the research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 616121 (HeteroIce project).A.M. is supported by the Royal Society through a Wolfson Research Merit Award.O.A.v.L. acknowledges funding from the Swiss National Science foundation (No.PPOOP2 138932).This research used resources as part of an IN-CITE project (awarded to D.A.) at the Oak Ridge National Laboratory, (Titan) which is supported by the Office of Science of the U.S. Department of Energy (DOE) under Contract No. DEAC05-00OR22725.This research also used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. DOE under contract DE-AC02-06CH11357.In addition, we are grateful for computing resources provided by the London Centre for Nanotechnology and University College London.We would like to thank G. Tocci, C. Gattinoni, and R. Ramakrishnan for useful discussions.

TABLE I .
Perpendicular separation distances (in Ångstrom) of the oxygen atom to the plane of the 1,2-azaborine ring (R O-plane ), and hydrogen atoms of water to C3 and C5 (R HW-C3 and R HW-C5 ).The shortest C-H distance is reported in each case.

TABLE II .
Absolute interaction energies in (meV) of water on 1,2-azaborine, using the structures shown in Fig.1.The benchmark values from ∆CCSD(T) are presented in addition to DFT, HF, MP2 and DMC results.Lowest energies for each method are highlighted in bold.
and optPBE-vdW: remarkably less than 3% errors.The other vdW functionals do not perform as well, with vdW-DF2 underbinding by 15% and revPBE-vdW underbinding by 25%.As anticipated, PBE is strongly underbinding by almost 35% due to the lack of longrange correlation, whilst dispersion corrected PBE-D2 overestimates the binding by 20%.

TABLE III .
34solute interaction energies of water to benzene, 1,2-azaborine, and borazine using PBE, optB86b-vdW, vdW-DF2 and CCSD(T).The absolute interaction energies with DFT correspond to optimized structures, and the values for water-1,2-azaborine correspond to the C3 complex.CCSD(T) interaction energy by Min et al. with CBS extrapolation.34b ∆CCSD(T) interaction energy calculated in Section III A for the C3 complex.cCCSD(T)interaction energy without CBS extrapolation by Wu et al.37  dAs part of this study our absolute interaction energy is −84 ± 4 meV for water-borazine using DMC.