Despite a recent flurry of experimental and simulation studies, an accurate estimate of the interaction strength of water molecules with hexagonal boron nitride is lacking. Here, we report quantum Monte Carlo results for the adsorption of a water monomer on a periodic hexagonal boron nitride sheet, which yield a water monomer interaction energy of −84 ± 5 meV. We use the results to evaluate the performance of several widely used density functional theory (DFT) exchange correlation functionals and find that they all deviate substantially. Differences in interaction energies between different adsorption sites are however better reproduced by DFT.

Hexagonal boron nitride (h-BN) has become popular for anyone with an interest in 2-dimensional materials, due to a number of notable properties such as high thermal conductivity, mechanical robustness, and exceptional resistance to oxidation,1 and not least because it is isostructural with graphene. Our interest has been piqued by experimental reports of fascinating behavior of water at h-BN such as superhydrophobicity,2 water cleaning ability,3 or generation of electric current.4 These experiments have already prompted a number of simulation studies of water on h-BN sheets and nanostructures using both density functional theory (DFT) and classical molecular dynamics.5–11 They have been incredibly informative and have helped to, e.g., understand the electrical currents generated in BN nanotubes.4 

However, there is one major unknown at the very heart of any water/BN simulation study: we simply do not know what the interaction strength of a water molecule with h-BN is. DFT calculations yield a range of values for the water monomer adsorption energy depending on the exchange-correlation (xc) functional used5–7 and force fields rely on interaction parameters fitted to particular xc functionals or to experimental data such as contact angles for macroscopic water droplets.8–11 If fitting to experiment, one needs to be certain that the experimental conditions are exactly known; recent lessons learned for water droplets on graphene reveal that contact angle measurements are incredibly sensitive to surface preparation conditions and levels of cleanliness.12–15 

The lack of well-defined reference data for water on h-BN is representative of a much broader problem: there are very few systems for which accurate water monomer adsorption energies have been established. Mainly, this is because even at low temperatures, water molecules cluster into larger aggregates making the determination of monomer adsorption energies with established surface science techniques such as temperature programmed desorption or single crystal adsorption calorimetry highly challenging.16–18 In the absence of experimental data, simulations play an important role, either via explicitly correlated quantum chemistry approaches or quantum Monte Carlo (QMC) (see, e.g., Refs. 19–25). Indeed given recent increases in computational capacity and the fact that it can be applied to periodic systems, QMC has emerged as a powerful technique for obtaining interaction energies of molecules with surfaces20,21 or biomolecules.26–28 

Here, we report results for interaction energy curves for water on a periodic h-BN sheet using fixed node diffusion Monte Carlo (DMC). We obtain from this an estimate of the water/h-BN interaction strength of about −84 ± 5 meV at an equilibrium water-surface distance of ca. 3.4 Å. We have also computed interaction energy curves with a variety of DFT xc functionals and we find that these differ significantly from DMC. Except for local density approximation (LDA), of the functionals considered, those that do not account for van der Waals (vdW) underbind and those that do overbind. DFT based predictions of the equilibrium adsorption height are much better with several functionals coming within 0.2 Å of DMC. In addition, based on DMC and DFT calculations of water on h-BN in a second metastable adsorption structure, we find that many of the xc functionals do reasonably well in predicting the relative energy difference between the stable and metastable adsorption structures.

Two different levels of theory have been used in this study, fixed node DMC and DFT. A standard computational setup has been used for each and so we only discuss the key features here.

QMC calculations were undertaken using the CASINO code,29 with Slater-Jastrow type trial wavefunctions in which the Jastrow factor contains electron-nucleus, electron-electron, and electron-electron-nucleus terms. We used Trail and Needs pseudopotentials30,31 for all atoms, in which the 1s electrons of B, N, and O were treated as core. This setup for the DMC calculations is similar to the one used in Ref. 28 where water adsorption was examined on 1,2-azaborine and agreement between DMC and coupled cluster with single, double and perturbative triple (CCSD(T)) excitations to within 9 meV was obtained. The initial single particle wavefunctions for use in DMC were obtained from DFT plane-wave calculations using the PWSCF package.32 A standard 300 Ry energy cutoff was applied and for efficiency, the resulting wavefunctions were expanded in terms of B-splines33 using a grid multiplicity of 2.0. Trial wavefunctions were generated using the LDA34 which has been validated for weak interactions in previous work.28,35 After optimization of the trial wavefunctions in variational Monte Carlo, we used 6 553 840 walkers across 16 384 cores for each point along the DMC interaction energy curves. The locality approximation was utilized36 with a time step of 0.015 a.u. which we tested against a time step of 0.005 a.u.

VASP 5.3.537–40 was used for the DFT calculations, making use of projector augmented wave (PAW) potentials41,42 to model the core regions of atoms (again, the 1s electrons of B, N, and O were treated as core). Following careful tests, we chose a 500 eV plane-wave cutoff and a (4 × 4) unit cell of h-BN with 16 Å between sheets, along with Γ-point sampling of reciprocal space.43 The proliferation of DFT xc functionals over the last decade44 means that we cannot possibly consider all xc functionals or even all modern xc functionals designed to capture weak interactions. Rather, we consider a small selection that have been widely used in adsorption studies. This includes the LDA, the PBE45 generalised gradient approximation (GGA), two hybrid functionals (PBE046,47 and B3LYP48–51), and several vdW inclusive functionals (PBE + D2,52 PBE + D3,53,54 DFT + vdW,55 optB86b-vdW,56–58 and vdW-DF256,59). The DFT + vdW correction (from Tkatchenko and Scheffler) was applied to three xc functionals (PBE, PBE0, and B3LYP).

Results for interaction energy curves of a water monomer on h-BN in two different adsorption modes (Figure 1), obtained from DMC and a range of DFT xc functionals, are shown in Figure 2. The interaction energy between the adsorbate and substrate is plotted as a function of the perpendicular distance between the oxygen atom of the water molecule and the h-BN sheet. The absolute interaction energy between water and the substrate was calculated as follows:

E int = E d tot E far tot ,
(1)

where E d tot is the total energy of water and h-BN at a given oxygen-surface separation distance, d, and E far tot is the total energy of water and h-BN at 8 Å oxygen-surface distance. This definition allows the same Jastrow factor to be used for all configurations, including the reference structure in DMC. Adsorption structures were obtained from optB86b-vdW optimizations of water on a fixed h-BN sheet, see supplementary material for coordinates of the adsorption and reference configurations.60 We chose to use the same structures for DMC and all xc functionals because this makes for a cleaner comparison. For many of the xc functionals, we have computed interaction energy curves with fully relaxed structures and the differences between the relaxed and the fixed structures are <5 meV, except in the repulsive wall at short oxygen-surface separations. The first adsorption structure considered has the oxygen of the water molecule above an N site with one of the OH bonds directed at that N atom (Figures 1(a) and 1(c)). This is the most stable adsorption structure according to previous DFT studies.5 The second structure has the oxygen atom of the water molecule above a B site with the plane of the molecule tilted away from the substrate by 128° (Figures 1(b) and 1(d)). According to our DFT calculations, this is the most stable structure for water at the B site but ∼20 meV less stable than the N site adsorption structure. We consider water adsorption at the B site to establish if the DFT site preference for this system is correct.61 

FIG. 1.

Structures of the two adsorption modes of water on h-BN considered in this study. (a) and (c) Top and side view of water above an N site of h-BN. (b) and (d) Top and side view of water above a B site of h-BN. Boron is pink, nitrogen is blue, oxygen is red and hydrogen is white. All calculations have been performed on periodic unit cells, with the periodic unit cell in the x,y plane indicated by the blue frames in (a) and (b).

FIG. 1.

Structures of the two adsorption modes of water on h-BN considered in this study. (a) and (c) Top and side view of water above an N site of h-BN. (b) and (d) Top and side view of water above a B site of h-BN. Boron is pink, nitrogen is blue, oxygen is red and hydrogen is white. All calculations have been performed on periodic unit cells, with the periodic unit cell in the x,y plane indicated by the blue frames in (a) and (b).

Close modal
FIG. 2.

(a) Interaction energy curves for water situated above the N site in h-BN as shown in the inset. (b) Interaction energy curves for water situated above the B site in h-BN as shown in the inset. The lines connecting the data points are merely there to guide the eye.

FIG. 2.

(a) Interaction energy curves for water situated above the N site in h-BN as shown in the inset. (b) Interaction energy curves for water situated above the B site in h-BN as shown in the inset. The lines connecting the data points are merely there to guide the eye.

Close modal

Let us now focus on the DMC interaction energy curves for water on h-BN. Because of the enormous computational cost of DMC, we can only compute a small number of points for each energy curve, which limits the resolution of the curves. Nonetheless, they are sufficiently well defined to yield an adsorption energy of −84 ± 5 meV at a height of ∼3.4 Å at the N site and an adsorption energy of −63 ± 5 at a height of ∼3.2 Å at the B site. The O atom sits slightly further away from the substrate at the N site because of the orientation of the molecule at this site, wherein the H atoms points to the N, forming a weak hydrogen bond like interaction. The relative energies of the two sites confirm the DFT site preference but more importantly provide an estimate of the water monomer interaction energy that is free of any arbitrary choices of DFT xc functional.

Obtaining an accurate estimate of the interaction strength between water and h-BN is important in its own right, however, it also provides a valuable benchmark which we now exploit. Here, we use our DMC derived interaction energy curves to evaluate how various DFT xc functionals perform for this system. Interaction energy curves from several functionals are included in Figure 2 and in some respects, these reveal a familiar story. Looking at the most stable site first, LDA overbinds by predicting an adsorption energy of −183 meV with the molecule 0.4 Å closer to the surface than DMC. In contrast, the GGA and the hybrid functionals underbind: PBE is ∼−45 meV, PBE0 is ∼−40 meV, and B3LYP is ∼−15 meV. The PBE and PBE0 adsorption heights are fairly reliable at 3.40 Å whereas the shallow B3LYP minimum is located at 3.55 Å. More interesting are the results from the vdW inclusive functionals since these are, in principle, designed to accurately describe weak interaction systems. Surprisingly, we find that all vdW inclusive functionals considered significantly overbind this adsorption system. Specifically, the adsorption energies are in the −140 to −170 meV range, with vdW-DF2 predicting the smallest adsorption energy and optB86b-vdW the largest. This overbinding also persists at large adsorbate-substrate distances; compare, for example, the DMC and vdW-inclusive DFT results at 4 − 5 Å from the surface. The predicted height above the surface is in reasonably good agreement with DMC, only around 0.1 Å closer to the surface for all vdW-inclusive functionals.

Moving to the B site adsorption structure, we find that systematically, with the exception of PBE + D2, the interaction strength is reduced by ∼20 − 30 meV. This is in very good agreement with the DMC energy difference between these two sites (PBE-D2 predicts that the B site is ∼60 meV less stable than the N site). Thus, although none of the xc functionals considered come within 40% of the DMC interaction energy, the change in interaction energies between adsorption sites is in most cases described fairly accurately. We note that we have considered just two adsorption structures and considerably more work is needed to fully substantiate this conclusion. Moreover, at this stage, we do not fully understand the poor performance of the vdW functionals and defer a more detailed analysis to a future publication in which results from on-going water adsorption calculations on BN clusters will also be presented.

In summary, we have obtained DMC interaction energy curves for water on a periodic hexagonal sheet of BN and used these to evaluate the performance of a number of xc functionals. The interaction energy obtained is −84 ± 5 meV. This is clearly a small number, corresponding to the physisorption regime.62 It is, however, about 15 meV larger than the value predicted by DMC for water on graphene. Interestingly, many of the van der Waals inclusive functionals also predict a similar 15–20 meV increase upon going from graphene to h-BN.5,20 This suggests that although interaction energies are overestimated with these functionals, the relative interaction energies between the two materials are fairly well described.

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 No. 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. PP00P2 138932). This research used resources as part of an INCITE 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 No. DE-AC02-06CH11357. In addition, we are grateful for computing resources provided by the London Centre for Nanotechnology and University College London.

1.
A.
Pakdel
,
Y.
Bando
, and
D.
Golberg
,
Chem. Soc. Rev.
43
,
934
(
2014
).
2.
A.
Pakdel
,
C.
Zhi
,
Y.
Bando
,
T.
Nakayama
, and
D.
Golberg
,
ACS Nano
5
,
6507
(
2011
).
3.
W.
Lei
,
D.
Portehault
,
D.
Liu
,
S.
Qin
, and
Y.
Chen
,
Nat. Commun.
4
,
1777
(
2013
).
4.
A.
Siria
,
P.
Poncharal
,
A.-L.
Biance
,
R.
Fulcrand
,
X.
Blase
,
S. T.
Purcell
, and
L.
Bocquet
,
Nature
494
,
455
(
2013
).
5.
G.
Tocci
,
L.
Joly
, and
A.
Michaelides
,
Nano Lett.
14
,
6872
(
2014
).
6.
Y.
Ding
,
M.
Iannuzzi
, and
J.
Hutter
,
J. Phys. Chem. C
115
,
13685
(
2011
).
7.
H.
Li
and
X. C.
Zeng
,
ACS Nano
6
,
2401
(
2012
).
8.
C. Y.
Won
and
N. R.
Aluru
,
J. Am. Chem. Soc.
129
,
2748
(
2007
).
9.
C.
Won
and
N.
Aluru
,
J. Phys. Chem. C
112
,
1812
(
2008
).
10.
M. C.
Gordillo
and
J.
Martí
,
Phys. Rev. E
84
,
011602
(
2011
).
11.
R. C.
Dutta
,
S.
Khan
, and
J. K.
Singh
,
Fluid Phase Equilib.
302
,
310
(
2011
).
12.
F.
Taherian
,
V.
Marcon
,
N. F. A.
van der Vegt
, and
F.
Leroy
,
Langmuir
29
,
1457
(
2013
).
13.
J.
Rafiee
,
X.
Mi
,
H.
Gullapalli
,
A. V.
Thomas
,
F.
Yavari
,
Y.
Shi
,
P. M.
Ajayan
, and
N. A.
Koratkar
,
Nat. Mater.
11
,
217
(
2012
).
14.
C.-J.
Shih
,
Q. H.
Wang
,
S.
Lin
,
K.-C.
Park
,
Z.
Jin
,
M. S.
Strano
, and
D.
Blankschtein
,
Phys. Rev. Lett.
109
,
176101
(
2012
).
15.
Z.
Li
,
Y.
Wang
,
A.
Kozbial
,
G.
Shenoy
,
F.
Zhou
,
R.
McGinley
,
P.
Ireland
,
B.
Morganstein
,
A.
Kunkel
,
S. P.
Surwade
 et al,
Nat. Mater.
12
,
925
(
2013
).
16.
A.
Hodgson
and
S.
Haq
,
Surf. Sci. Rep.
64
,
381
(
2009
).
17.
C. T.
Campbell
and
J. R. V.
Sellers
,
Chem. Rev.
113
,
4106
(
2013
).
18.
J.
Carrasco
,
A.
Hodgson
, and
A.
Michaelides
,
Nat. Mater.
11
,
667
(
2012
).
19.
B.
Li
,
A.
Michaelides
, and
M.
Scheffler
,
Surf. Sci.
602
,
L135
(
2008
).
20.
J.
Ma
,
A.
Michaelides
, and
D.
Alfè
,
J. Chem. Phys.
134
,
134701
(
2011
).
21.
J.
Ma
,
A.
Michaelides
,
D.
Alfè
,
L.
Schimka
,
G.
Kresse
, and
E.
Wang
,
Phys. Rev. B
84
,
033402
(
2011
).
22.
E.
Voloshina
,
D.
Usvyat
,
M.
Schütz
,
Y.
Dedkov
, and
B.
Paulus
,
Phys. Chem. Chem. Phys.
13
,
12041
(
2011
).
23.
G. R.
Jenness
,
O.
Karalti
, and
K. D.
Jordan
,
Phys. Chem. Chem. Phys.
12
,
6375
(
2010
).
24.
L.
Shulenburger
and
T. R.
Mattsson
,
Phys. Rev. B
88
,
245117
(
2013
).
25.
D.
Lu
,
Y.
Li
,
D.
Rocca
, and
G.
Galli
,
Phys. Rev. Lett.
102
,
206411
(
2009
).
26.
M.
Dubecký
,
P.
Jurečka
,
R.
Derian
,
P.
Hobza
,
M.
Otyepka
, and
L.
Mitas
,
J. Chem. Theory Comput.
9
,
4287
(
2013
).
27.
A.
Benali
,
L.
Shulenburger
,
N. A.
Romero
,
J.
Kim
, and
O. A.
von Lilienfeld
,
J. Chem. Theory Comput.
10
,
3417
(
2014
).
28.
Y. S.
Al-Hamdani
,
D.
Alfè
,
O. A.
von Lilienfeld
, and
A.
Michaelides
,
J. Chem. Phys.
141
,
18C530
(
2014
).
29.
R. J.
Needs
,
M. D.
Towler
,
N. D.
Drummond
, and
P.
López Ríos
,
J. Phys.: Condensed Matter
22
,
023201
(
2010
).
30.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
122
,
174109
(
2005
).
31.
J. R.
Trail
and
R. J.
Needs
,
J. Chem. Phys.
122
,
014112
(
2005
).
32.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
 et al,
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
33.
D.
Alfè
and
M. J.
Gillan
,
Phys. Rev. B
70
,
161101(R)
(
2004
).
34.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
35.
J.
Ma
,
D.
Alfè
,
A.
Michaelides
, and
E.
Wang
,
J. Chem. Phys.
130
,
154303
(
2009
).
36.
L.
Mitas
,
E. L.
Shirley
, and
D.
Ceperley
,
J. Chem. Phys.
95
,
3467
(
1991
).
37.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
38.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
49
,
14251
(
1994
).
39.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
40.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
41.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
42.
G.
Kresse
and
D.
Joubert
,
Phys. Rev. B
59
,
1758
(
1999
).
43.

For example, tests with a higher plane-wave cutoff (600 eV) and denser k-point mesh (5 × 5 × 1) performed for the PBE functional yielded an interaction energy that differed from the reported one by <3 meV. Similarly when we tested the current setup for water adsorption on the 1,2-azaborine system considered in Ref. 28 against all-electron PBE calculations with an aug-cc-PV5Z basis set, we found that the results with the two approaches differed by only 1 meV, the PBE adsorption energy for that system being 109-110 meV.

44.
K.
Burke
,
J. Chem. Phys.
136
,
150901
(
2012
).
45.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
46.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
47.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
,
J. Chem. Phys.
105
,
9982
(
1996
).
48.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
49.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
50.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
,
Can. J. Phys.
58
,
1200
(
1980
).
51.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
,
J. Phys. Chem.
98
,
11623
(
1994
).
52.
S.
Grimme
,
J. Comput. Chem.
27
,
1787
(
2006
).
53.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
54.
S.
Grimme
,
S.
Ehrlich
, and
L.
Goerigk
,
J. Comput. Chem.
32
,
1456
(
2011
).
55.
A.
Tkatchenko
,
R. A.
DiStasio
,
R.
Car
, and
M.
Scheffler
,
Phys. Rev. Lett.
108
,
236402
(
2012
).
56.
M.
Dion
,
H.
Rydberg
,
E.
Schröder
,
D. C.
Langreth
, and
B. I.
Lundqvist
,
Phys. Rev. Lett.
92
,
246401
(
2004
).
57.
A. D.
Becke
,
J. Chem. Phys.
84
,
4524
(
1986
).
58.
J.
Klimeš
,
D. R.
Bowler
, and
A.
Michaelides
,
Phys. Rev. B
83
,
195131
(
2011
).
59.
K.
Lee
,
É. D.
Murray
,
L.
Kong
,
B. I.
Lundqvist
, and
D. C.
Langreth
,
Phys. Rev. B
82
,
081101
(
2010
).
60.
See supplementary material at http://dx.doi.org/10.1063/1.4921106 for structural details of the configurations used to obtain interaction energies.
61.
P. J.
Feibelman
,
B.
Hammer
,
J. K.
Nørskov
,
F.
Wagner
,
M.
Scheffler
,
R.
Stumpf
,
R.
Watwe
, and
J.
Dumesic
,
J. Phys. Chem. B
105
,
4018
(
2001
).
62.

Zero point energy contributions (computed within the harmonic approximation) weaken the optB86b-vdW interaction strength by ∼30 meV. Since this is the most strongly binding xc functional, others are likely to show a smaller reduction than this.

All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.

Supplementary Material