Toward Accurate Adsorption Energetics on Clay Surfaces.

Clay minerals are ubiquitous in nature, and the manner in which they interact with their surroundings has important industrial and environmental implications. Consequently, a molecular-level understanding of the adsorption of molecules on clay surfaces is crucial. In this regard computer simulations play an important role, yet the accuracy of widely used empirical force fields (FF) and density functional theory (DFT) exchange-correlation functionals is often unclear in adsorption systems dominated by weak interactions. Herein we present results from quantum Monte Carlo (QMC) for water and methanol adsorption on the prototypical clay kaolinite. To the best of our knowledge, this is the first time QMC has been used to investigate adsorption at a complex, natural surface such as a clay. As well as being valuable in their own right, the QMC benchmarks obtained provide reference data against which the performance of cheaper DFT methods can be tested. Indeed using various DFT exchange-correlation functionals yields a very broad range of adsorption energies, and it is unclear a priori which evaluation is better. QMC reveals that in the systems considered here it is essential to account for van der Waals (vdW) dispersion forces since this alters both the absolute and relative adsorption energies of water and methanol. We show, via FF simulations, that incorrect relative energies can lead to significant changes in the interfacial densities of water and methanol solutions at the kaolinite interface. Despite the clear improvements offered by the vdW-corrected and the vdW-inclusive functionals, absolute adsorption energies are often overestimated, suggesting that the treatment of vdW forces in DFT is not yet a solved problem.

Introduction

The accurate treatment of the adsorption of molecules on surfaces is a major challenge of materials modeling, with important applications in nanotechnologies and science: heterogeneous catalysis, sensors, corrosion, lubrication, friction, and coatings, to name but a few. An important case to study is that of clays. Clay minerals are natural aluminosilicates that find use in a wide variety of fields such as medicine, adhesives, paints, and oil drilling.[1−9] They also act as catalysts to ice nucleation in the atmosphere[10] and help cleanse soils and groundwater through adsorption of pollutants. A clear understanding of how molecules interact with the surfaces of clays is of the utmost importance to understand, improve, and control such processes.

Reliable reference data from theory and simulation are of intrinsic value and often important as a complement to experiments.[11,12] Computer simulations of water–surface interactions, at the molecular level, are often based on force fields (FF) and density functional theory (DFT) approaches.[13−15] Although these techniques are incredibly powerful and useful, there are cases where their accuracy is not satisfactory. FF potentials have parameters that have to be fit in order to reproduce experimental results or higher lever theoretical benchmarks, and this is not always straightforward. DFT is traditionally more accurate than FFs but at a larger computational cost. Unfortunately, DFT results are highly sensitive to the choice of the exchange-correlation (XC) functional used, and nowadays there are countless XC functionals to choose from.[16,17] In the field of materials science, the description of weak bonding interactions, and in particular London dispersion forces, is one of the most important challenges. Immense progress has been made in this area recently;[18,19] however, there is no rigorous way to systematically improve XC functionals, and as a result validation with alternative methods is needed. Of the various high level reference methods available,[20−28] quantum Monte Carlo (QMC) is a powerful approach for obtaining benchmark values for solids, surfaces, and large molecular systems. QMC, within the fixed node diffusion Monte Carlo (DMC) approach, has already been used to tackle interesting materials science problems that have been beyond the reach of DFT (see, e.g., refs (29−42)). This has provided reference data which have exposed shortcomings in existing FF models and DFT XC functionals, which in turn aids the development of such approaches.

In this paper we will use two QMC approaches to investigate molecular adsorption on a clay surface, namely, DMC and lattice regularized diffusion Monte Carlo[27,28] (LRDMC). The particular clay we will examine is kaolinite (Al4Si4O10(OH)8), as shown in Figure . Since the first outline of the kaolinite crystal structure by Pauling in 1930,[43] numerous structural studies using X-ray and neutron powder diffraction,[44−47] X-ray single crystal,[48] and electron diffraction methods,[49] as well as theoretical studies[50−53] have been carried out on kaolinite. Consequently it is one of the most suitable aluminosilicate clays to assess the performance of various theoretical methods. In addition, when looking at adsorption processes on kaolinite, cleavage along the (001) basal plane leads to exposure of either aluminate or silicate faces (Figure ). The aluminate (AlOH) face is terminated in hydroxyl groups and as a result is regarded as hydrophilic, whereas the silicate (SiO2) face which exposes saturated Si–O groups is considered to be hydrophobic. The distinct chemical nature of these two surfaces means that adsorbates will interact differently with them, making kaolinite an interesting case study for understanding the role of vdW forces on the adsorption at clay mineral surfaces.

Figure 1

Representation of the kaolinite structure. The hydrogens are sketched in white, the oxygens in red, the silicons by yellow tetrahedra, and the aluminums by pink octahedra. The conventional unit cell is indicated by the blue line. The figure on the left illustrates the layered bulk. The figures on the right are the hydroxyl-terminated face (top) and the silicate-terminated face (bottom). Various adsorption sites on the hydroxyl- and silicate-terminated face are labeled.

In what follows, we will provide benchmark values for the adsorption of water and methanol molecules at the pristine hydroxyl- and silicate-terminated (001) faces of kaolinite, by using DMC and LRDMC. Then, we probe DFT XC functionals by considering a range of generalized gradient approximation (GGA) functionals, hybrid functionals, and dispersion-corrected density functionals which account for vdW forces. In the case of molecules adsorbed on kaolinite we find that the bare GGAs and hybrids are quite unreliable: as expected adsorption energies are underestimated, but more importantly, the relative adsorption energies of water versus methanol do not even agree with QMC. Moreover, on the silicate-terminated face the molecules hardly bind at all and move quite far from the surface during geometry optimizations. Accounting for vdW forces improves adsorption energies significantly and stabilizes the structures. However, most of the vdW-corrected and vdW-inclusive functionals predict adsorption energies which are slightly too large compared to QMC. This indicated that there remains room for improvement in terms of how vdW forces are handled in DFT.

The remainder of this paper is organized as follows. In Section 2 we outline the key computational details of our simulations. Since a range of techniques has been used, for brevity the more detailed descriptions of the computational setups are provided in the Supporting Information. QMC and DFT evaluations of the adsorption of water and methanol at both (001) faces of kaolinite are reported and discussed in Section 3. Finally, we summarize our results and draw conclusions in Section 4.

Methods and Computational Setup

Adsorption Energy Evaluation

Adsorption was examined on a single layer of kaolinite, a system with 2D periodicity along the A and B axes as indicated in Figure . The simulated supercell was 1 × 2 the conventional unit cell of bulk kaolinite (ca. 10.38 Å × 9.01 Å) and is comprised of 8 aluminums, 8 silicons, 36 oxygens, and 16 hydrogens for the kaolinite slab plus the atoms of the adsorbed molecule. Note that with ca. 300 electrons the simulations are large for QMC calculations.

The adsorption energy, EadsM@X, of the molecule M at the face X of the kaolinite layer, where M can either be the water (H2O) or the methanol (MeOH) and X can be the hydroxyl-terminated (AlOH) face or the silicate-terminated (SiO2) face, can be evaluated in two ways: the first method, hereafter called complex-minus-fragments, is computed aswhere Eslab+M@X is the total energy for the system with M at the X-face of the kaolinite slab, and Eslab and EM are the total energies of the isolated slab and the isolated molecule M, respectively. The second method, hereafter called complex-minus-far, is computed aswhere Eslab–M is the total energy of a system where the kaolinite slab and the molecule M are far enough apart that their interaction is negligible. The two methods are equivalent if and only if: (i) the size effects due to the periodicity of the system are negligible and (ii) the electronic structure calculations are performed with methods that are exactly size consistent. If these conditions are not satisfied in general we have that Eslab–M ≠ EM + Eslab, meaning that eqs and 2 provide different evaluations of the adsorption energies. In particular, whenever size effects are detected, the complex-minus-far method usually benefits from a larger error cancellation. On the other hand, in cases where size effects are negligible and electronic structure methods are size consistent, there are no residual interactions between the molecule and the periodic partners, and then the complex-minus-fragments method is usually to be preferred. The reason for that is the computational cost: for a system with N electrons the computation is proportional to Nγ, with γ > 1 (e.g., in DFT γ is typically between 2 and 3 and in QMC between 3 and 4), so the cost for calculations of Eslab and EM is cheaper than Eslab–M. Moreover, when several adsorption energies need to be evaluated, Eslab is calculated only once, whereas a different calculation of Eslab–M has to be performed for each molecule M.

QMC Calculations

The two QMC approaches used are DMC and LRDMC. They are both projection Monte Carlo methods: they can access the electronic ground state energy of the system by iteratively projecting an initial trial wave function ψT into the ground state, with the constraint that the projected wave function Φ has the same nodal surface of an appropriately chosen guiding function ψG (fixed node approximation).[20,54] Both the trial and the guiding wave functions are parametrized functions, and they have to be the best approximation of the ground state that we can provide (given the constraint of their ansatz). Thus, usually they are taken such that ψT = ψG = ψVMC, where ψVMC is the best function obtained within a variational Monte Carlo approach, with the variational parameters optimized in order to minimize either the variational energy or the variance. Whenever ψG has the exact nodal surface, the approach is exact; otherwise, it gives the best approximation of the ground state given the fixed node constraint.

In projection Monte Carlo approaches there is a second approximation in how the projection is performed, and it is different in DMC and LRDMC. The projection in DMC comes from the imaginary time Schrödinger equation; it is implemented as an imaginary time evolution, where a time step τ has to be chosen. The chosen τ is a trade-off between accuracy and computational cost: exact results are obtained for τ → 0, but the computational cost is ∝1/τ. The finite time-step error can be controlled by performing several calculations with different values of τ and finally extrapolating to the continuum limit τ → 0. However, in big systems like those considered here, the extrapolation is impractical and sometimes unfeasible or unreliable,[55] but it is sufficient to consider the results for a τ small enough that the expected finite-time error is smaller than the required accuracy. Here, we have chosen τ in order to have an expected time-step error smaller than the stochastic error of the evaluations (see Section I in the Supporting Information). On the other hand, LRDMC is based on the spatial discretization of the molecular Hamiltonian on a lattice of mesh size a, and it resorts to the projection scheme used also in the Green function Monte Carlo algorithm.[56,57] The error induced by the finite mesh size a is analogous to the time-step error appearing in standard DMC calculations. LRDMC preserves the variational principle even when used in combination with nonlocal pseudopotentials[27] (PPs), and it is size consistent for any value of the mesh a, maintaining its efficiency even for systems with a large number of electrons.[28]

Both DMC and LRDMC provide excellent benchmark values for weakly interacting systems, as established in a number of studies.[29−36,39,58,59] We used here a standard setup, described in detail in the Supporting Information. The stochastic error associated with the QMC evaluations of the adsorption energy is ca. 20 meV. The systems under consideration are too large for a QMC-based structural optimization, even at the variational level, so the reference structures were those obtained from the PBE-D3 functional, as described in Section 3.1. As we will see in Section 3, PBE-D3 configurations are in good agreement with those obtained by all the other vdW-inclusive functionals, thus the bias given by the use of PBE-D3 configurations is expected to be small compared to the stochastic error of the DMC evaluation.

There is one aspect of QMC simulations that deserves special care in this specific system, namely, the finite-size errors (FSEs).[60−64] QMC is a many-body method, and in contrast to (effective) one-particle methods such as DFT, QMC cannot simply exploit Bloch’s theorem in calculations for extended periodic systems. FSEs can be taken into account by performing simulations in larger periodic supercells, through the twist-average method,[60] corrections to the Ewald energy,[61] or the Kwee, Zhang, Krakauer (KZK) method.[62] In this work we have used the KZK method (see Section 3 in the Supporting Information).

DFT Calculations

There is by now an almost limitless variety of DFT XC functionals that we could examine.[65] Here we restrict ourselves to the LDA functional;[66] two GGA functionals, PBE[67,68] and RPBE;[69] two hybrid functionals, PBE0[70] and B3LYP;[71−74] three vdW-corrected PBE functionals PBE-D2,[75] PBE-D3,[76] (both from Grimme, D3 correction with “zero-dampling”[77]), and PBE+vdW(TS) from Tkatchenko and Scheffler;[78] two vdW-corrected hybrid functionals, PBE0-D3 and B3LYP-D3[76] (both from Grimme[76]); and four self-consistent nonlocal functionals (often called vdW-inclusive functionals), the original vdW-DF from Dion (also named revPBE-vdW),[79,80] the second generation vdW-DF2,[81] as well as optPBE-vdW[82] and optB86b-vdW[83] from Klimeš et al. We stress that the latter four vdW-inclusive functionals are actually based on GGAs and basically differ from the vdW-corrected GGA functionals (e.g., PBE-D2, PBE-D3, and PBE+vdW(TS)) only in the way the dispersion energy is approximated.[18,19] Other functionals and vdW corrections have been tested, and results obtained using a comprehensive set of approaches is reported in Table S1 of the Supporting Information. Adsorption energies were evaluated using the complex-minus-fragments method (see eq ), but the results are the same as those obtained with the complex-minus-far method, as expected. Further details about the setup of the DFT calculations are reported in Section 4 of the Supporting Information.

Molecular Dynamics Simulations

We also performed a series of molecular dynamics simulations using classical force fields for aqueous water–methanol solutions on kaolinite. The kaolinite slab was modeled as a single sheet of kaolinite (approximately 31 × 36 Å) using the CLAYFF force field,[84] and the OH bond lengths were constrained using the P-LINCS algorithm.[85] Above this slab 538 TIP4P/2005[86] water and 230 OPLS/UA methanol[87] molecules were randomly placed in order to create a liquid film on the kaolinite surface. The standard Lorentz–Berthelot mixing rules were used to compute cross-interactions, except to adjust the adsorption energies as detailed below. Using the GROMACS 4.5 simulation package,[88] constant volume and temperature dynamics were propagated using a leapfrog integrator and a Nosè-Hoover chain thermostat, along with replica-exchange among eight replicas with temperatures ranging from 275 to 310 K in 5 K intervals. Real-space interactions were truncated at 9 Å with corrections to the energy applied, and particle-mesh Ewald was used to account for long-range electrostatics[89,90] with the corrections for the slab geometry of the system.[91] A time step of 2 fs was used and molecular dynamics simulations performed for at least 11 ns, with the first nanosecond being disregarded as equilibration.

Adsorption energies were computed after geometry optimization using the complex-minus-fragments method. This was done in three ways: First, the adsorption energy was computed by applying the standard mixing rules. This yielded adsorption energies of 642 and 640 meV for water and methanol, respectively, and ΔEads = −2 meV. Second, the strength of the Lennard-Jones interaction between the CH3 group of methanol and the oxygen atoms of the kaolinite OH groups was adjusted such that ΔEads matched that of PBE. Finally, the same was done to match ΔEads obtained by DMC.

Results

Reference Structures for Water and Methanol Adsorbed on Kaolinite

Water adsorption on the hydroxyl-terminated face of kaolinite has been studied experimentally[92−94] and theoretically,[95−101] whereas adsorption on the silicate-terminated face is less well studied. Very little is known about methanol adsorption on either face. In the following, the most stable adsorption structures identified for water and methanol on the two kaolinite surfaces are presented.

On each surface a range of adsorption sites were considered, as indicated in Figure . According to the number of H-bonds formed between the adsorbate and the surface, the adsorption sites can be classified into three categories: 3-fold, 2-fold, and 1-fold sites. The most stable configurations obtained, using DFT with the PBE-D3 functional, are shown in Figures A–2H. These structures have been taken as the reference for DMC, LRDMC, and the other DFT calculations. In addition, starting from the reference PBE-D3 structures, we have relaxed the geometries for each of the different functionals considered, as shown in Figures A–3D. Figure E compares the distance of the molecules from the slab as obtained with different functionals.

Figure 2

Adsorption of water and methanol on the hydroxyl-terminated and the silicate-terminated faces of kaolinite (side view in first row, top view in second row). Geometries are relaxed using the PBE-D3 functional and have been taken as reference for the other calculations. The adsorbed molecule on kaolinite is depicted in cyan and gray, and the H-bonds are represented by the blue dashed lines.

Figure 3

Panels A, B, C, and D show the most stable DFT structures for the adsorption of water and methanol on the hydroxyl- and silicate-terminated faces of kaolinite, provided by the different XC functionals considered. The color scheme for the various functionals is blue for PBE, cyan for RPBE, white for PBE-D2, black for PBE-D3 (that is also the reference for QMC calculations), pink for PBE+vdW(TS), violet for vdW-DF, green for vdW-DF2, gray for optPBE-vdW, and yellow for optB86b-vdW. Panel E shows the height of the center-of-mass of the adsorbed molecules from the average surface plane defined by the surface oxygens for the different XC functionals. The four dashed horizontal lines correspond to the values for the reference PBE-D3 structures.

Concerning water at the hydroxyl-terminated face (H2O on AlOH), all structures initially put in 2-fold and 1-fold sites (A5–A8) moved to the 3-fold site A1. This preference for the A1 site agrees with previous DFT studies with local[99] and semilocal[96] XC functionals. In the most stable configuration, shown in Figures A and 2E, the C2 axis of the water molecule lies almost parallel to the plane of the surface. The water molecule donates one H-bond (OH distance of 1.69 Å) to and accepts two H-bonds (2.01 and 2.04 Å) from the surface (PBE-D3 values).

Let us now consider the adsorption of methanol at the hydroxyl-terminated face (MeOH on AlOH). One way of viewing methanol is as a water molecule with one of its hydrogen atoms replaced by a methyl group. This leads to two possible types of interaction with the surface: (i) hydrogen bond formation with the hydroxyl functional group and (ii) dispersion interactions arising from the −CH3 group. All calculations in which the methanol began parallel to the surface ended with the methanol perpendicular to the surface, maximizing the distance between the −CH3 group and the kaolinite. The −CH3 group can therefore be considered a “spectator group” that does not participate directly in the adsorption on the surface. The adsorption of methanol is therefore very similar to that of water, and indeed, we find that A1 is the most favorable site, with the methanol donating one H-bond to and accepting two from the surface. As was the case for the water structure, the H-bond donated by the methanol is much stronger than the two H-bonds it accepts: 1.68 vs 1.97 and 2.03 Å, respectively, with the PBE-D3 functional. The most stable configuration is shown in Figures B and 2F.

As noted above, adsorption of water at the silicate-terminated face (H2O on SiO2) is less well studied than adsorption at the hydroxyl-terminated face.[98−101] Of the six adsorption sites (S1–S6) considered here, the 1-fold S5 site turned out to be the most stable at the GGA level, and the 2-fold S1 generally is the most stable for the vdW-corrected and vdW-inclusive functionals. The PBE-D3 structure is depicted in Figures C and 2G.

The most stable structure found for methanol at the silicate-terminated face (MeOH on SiO2) using the PBE-D3 functional is shown in Figures D and 2H. The leading interaction here is dispersion; there is no H-bond-like interaction because the OH group of the methanol is parallel to the surface of the slab.

Benchmark Results from DMC and LRDMC

The DMC and LRDMC results for water and methanol adsorption on the two faces of kaolinite are reported in Table . As mentioned in the previous section, “bare” DMC and LRDMC evaluations have to be corrected for finite-size effects, and in our LRDMC simulations there is also an unphysical dipole interaction between slabs due to the 3D periodicity employed. The DMC calculations have been performed with 2D periodicity and so do not suffer from the latter problem. The bare and corrected results are reported in Table . From this it can be seen that our best estimates of the adsorption energy of water on the hydroxyl-terminated face are −648 ± 18 meV with DMC and −675 ± 14 meV with LRDMC. For methanol our best estimates of the adsorption energy are −694 ± 18 meV with DMC and −701 ± 13 meV with LRDMC. We notice that we are in the chemisorption regime both for water and for methanol, although the adsorption energy of methanol is slightly larger. Note that for both molecules the DMC and LRDMC evaluations are in good agreement, with the differences falling within the stochastic error of the evaluations. This shows that fixed-node projection QMC schemes are robust approaches: they are only slightly affected by the actual computational setup and implementation. Nonetheless, two slightly different adsorption energies for each case are obtained, and we should choose only one of them to use as our benchmark. We feel that in the specific case considered here the DMC values are likely to be more reliable since they have been obtained in 2D, as opposed to the LRDMC results which have been corrected for the dipole in the 3D cell. Moreover, the reported LRDMC evaluations use eq , which has larger FSE than the reported DMC evaluations, which use eq .

Table 1

DMC and LRDMC Evaluations (in meV) of the Adsorption Energy of Water and Methanol Molecules on the Hydroxyl- and Silicate-Terminated Faces of Kaolinite and the Water Minus Methanol Difference, ΔEads = EadsH – EadsMeOH, for Each Face of Kaolinite (ΔEads is Positive When Methanol Is More Strongly Adsorbed, Negative Otherwise)a

	hydroxyl-terminated face			silicate-terminated face
	H2O	MeOH	ΔEads	H2O	MeOH	ΔEads
bare DMC (eq 2)	–632 ± 18	–677 ± 18	45 ± 25	–172 ± 23	–236 ± 18	64 ± 29
FSE correction	–16	–17	+1	–12	–14	+2
corrected DMC	–648 ± 18	–694 ± 18	46 ± 25	–184 ± 23	–250 ± 18	66 ± 29
bare LRDMC (eq 1)	–674 ± 14	–736 ± 13	64 ± 13
FSE correction	+35	+73	–38
dipole correction	–36	–38	+2
corrected LRDMC	–675 ± 14	–701 ± 13	26 ± 13

As discussed in the text, bare DMC and LRDMC results are affected by finite-size errors (see Section 3 and Table S2 in the Supporting Information) that we have estimated and corrected the adsorption energies for accordingly. In addition, bare LRDMC evaluations are affected by an unphysical dipole–dipole interaction between the periodic slabs (because in this case 2D periodicity was not available, and we had to use 3D periodicity), thus we have included a dipole interaction correction. LRDMC simulations have not been performed for adsorption on the SiO2 face.

Having compared the results of the two QMC approaches on the hydroxyl-terminated face, we have only performed a DMC evaluation on the silicate-terminated face. The DMC adsorption energy at the silicate-terminated face is −184 ± 23 meV for water and −250 ± 18 meV for methanol. The methanol adsorbs more strongly than water, as for the hydroxyl-terminated face; however, in this case the adsorption is weaker, and we are in the physisorption regime.

Evaluation of DFT XC Functionals: Adsorption Energies and Structures

We now examine how the various DFT XC functionals considered in this study perform for water and methanol adsorption on the two faces of kaolinite.

Water Adsorption on the Hydroxyl-Terminated Face of Kaolinite

In Table and Figure we summarize the adsorption energies obtained with the different density functionals. At the GGA level PBE and RPBE give significantly different adsorption energies, with RPBE providing a value that is roughly 50% that of PBE. In line with the smaller adsorption energy we also see that the bonds the molecule makes with the surface with the RPBE functional are considerably longer than what is obtained with PBE. Specifically, with RPBE the two H bonds accepted from the surface are 2.30 and 2.36 Å, and the one donated is 1.81 Å, versus 2.03, 2.06, and 1.70 Å with PBE. Including dispersion interactions does not drastically change the geometry of the adsorbed water monomer at the hydroxyl-terminated face: the bond lengths at the PBE-D2, PBE-D3, PBE+vdW(TS), and opt-B86b-vdW level are slightly shortened, but they remain within 0.05 Å of the PBE structure. PBE-D2 predicts the shortest distance from the surface and the shortest H-bonds. The other functionals give H-bond lengths between the values provided by PBE and RPBE. From the shortest to the longest interaction distance, the functionals are ranked in the following order: PBE-D2 < PBE-D3 ∼ optb86b-vdW ∼ PBE+vdW(TS) < PBE < optPBE-vdW < vdW-DF2 < vdW-DF < RPBE. This trend roughly follows the sequence of adsorption energy predicted by the functionals and is also consistent with previous studies of DFT XC functionals for hydrogen-bonded systems.[58,102−104] Relaxation from the PBE-D3 geometry (performed for all the GGA and GGA+vdW functionals) results in rather small increases in adsorption energies. The maximum difference is observed for vdW-DF, with an increase of 36 meV upon relaxation.

Table 2

Adsorption Energy of Water, EH, and of Methanol, EMeOH, on the Hydroxyl-Terminated Face of Kaolinite and Adsorption Energy Difference, ΔEads = EadsH – EadsMeOH, between Water and Methanol, Obtained with DMC and Several DFT XC Functionalsa

	hydroxyl-terminated face
method	EH2Oads	EMeOHads	ΔEads
DMC	–648 ± 18	–694 ± 18	46 ± 25
LDA	–1102	–1138	36
GGA functionals
PBE	–607(−608)	–614(−616)	7(8)
RPBE	–360(−381)	–354(−380)	–6(−1)
hybrid functionals
PBE0	–599	–615	16
B3LYP	–524	–528	4
GGA+vdW functionals
PBE-D2	–822(−826)	–879(−882)	57(56)
PBE-D3	–767	–829	62
PBE+vdW(TS)	–769(−764)	–841(−833)	72(69)
vdW-DF	–530(−566)	–597(−635)	66(69)
vdW-DF2	–616(−641)	–658(−686)	42(44)
optPBE-vdW	–689(−699)	–767(−779)	78(80)
optB86b-vdW	–751(−752)	–835(−836)	84(85)
hybrid+vdW functionals
PBE0-D3	–768	–840	72
B3LYP-D3	–776	–849	72

The best performing functional is indicated in bold. All energy values are in meV. Energies have been obtained on PBE-D3 optimized structures, but in parentheses we also report the adsorption energies when geometries are fully relaxed consistently for each GGA and GGA+vdW functional.

Figure 4

Adsorption energies on kaolinite obtained by various XC functionals and DMC, for water on the hydroxyl-face (H2O on AlOH, blue upper triangles), methanol on the hydroxyl-face (MeOH on AlOH, red squares), water on the silicate face (H2O on SiO2, cyan lower triangles), and methanol on the silicate face (MeOH on SiO2, orange diamonds). Filled points represent the values with the reference structures (obtained using PBE-D3), and empty points (reported only for GGA and GGA+vdW functionals) correspond to relaxed structures for the specific functional. The solid lines are the reference DMC adsorption energies, and the shaded areas show the stochastic error.

A comparison with the QMC adsorption energies shows that vdW-DF2 and optPBE-vdW yield the best agreement, with the former providing a slightly underestimated adsorption energy (by −32 ± 18 meV) and the latter a slightly overestimated one (by 41 ± 18 meV). It also appears that the two GGA functionals (PBE and RPBE), the two hybrid functionals (PBE0 and B3LYP), and vdW-DF underestimate the interaction energy, whereas all the other functionals (PBE-D2, PBE-D3, PBE+vdW(TS), optB86b-vdW, PBE0-D3, and B3LYP-D3) overestimate the interaction. In particular, evaluations of the adsorption using vdW-corrected hybrid functionals do not seem to improve significantly compared to the GGA+vdW approaches.

Methanol Adsorption on the Hydroxyl-Terminated Face of Kaolinite

A careful investigation shows that the ordering of the functionals according to the H-bond lengths and to the adsorption energy is almost the same as that observed for water. The only exception is vdW-DF, which for methanol gives a larger adsorption energy than that obtained with PBE. The comparison with DMC confirms, as for water, that the best performing functionals are vdW-DF2 and optPBE-vdW. Again the GGA functionals and vdW-DF underestimate the adsorption energy, while all the other functionals (PBE-D2, PBE-D3, PBE+vdW(TS), optB86b-vdW, PBE0-D3, and B3LYP-D3) overestimate the interaction. As for the previous system, vdW-corrected hybrid functionals do not seem to improve significantly with respect to the GGA+vdW approaches.

Before discussing adsorption on the silicate face of kaolinite, we briefly compare water and methanol adsorption. An important finding from the results presented in Table is that with the GGA functionals the adsorption energies of water and methanol are similar (e.g., Eads,PBEH = −608 meV and Eads,PBEMeOH = −616 meV), but upon inclusion of dispersion interactions methanol is stabilized to a greater extent (e.g., Eads, optPBE–vdWH = −699 meV and Eads,optPBE–vdWMeOH = −779 meV). This is apparent in Figure , where the difference in adsorption between water and methanol is plotted. Therefore, even though the methyl group is considered a spectator, its vdW interaction with the surface is non-negligible, and it is clearly desirable to properly account for dispersion interactions in these systems. DMC confirms the reliability of the vdW-inclusive functionals on this issue as DMC also finds that methanol binds more strongly than water.

Figure 5

Difference in adsorption energy, between water and methanol on the hydroxyl-terminated (left panel) and silicate-terminated (right panel) faces of kaolinite, as obtained from various XC functionals and DMC. Positive values mean that methanol binds more strongly than water. Filled points represent the values for the configurations optimized using PBE-D3, and empty points (reported only for GGA and GGA+vdW functionals) correspond to relaxed structures for the specific functional. The solid lines are the reference DMC adsorption energies, and the shaded areas show the stochastic error.

Water Adsorption on the Silicate-Terminated Face of Kaolinite

In Table we present the results from all the functionals for adsorption on the silicate-terminated face. Irrespective of which functional is used the adsorption energies obtained are in the physisorption regime. Consequently, the inclusion of vdW forces is expected to have a more obvious impact than on the hydroxyl-terminated face.

Table 3

Adsorption Energy of Water, EH, and of Methanol, EMeOH, on the Silicate-Terminated Face of Kaolinite and Adsorption Energy Difference, ΔEads = EadsH – EadsMeOH, between Water and Methanol, Obtained with DMC and Several DFT XC Functionalsa

	silicate-terminated face
Method	EH2Oads	EMeOHads	ΔEads
DMC	–184 ± 23	–250 ± 18	66 ± 29
LDA	–295	–315	20
GGA functionals
PBE	–85(−92)	–93(−118)	8(27)
RPBE	25(−41)	21(−185)	4(144)
hybrid functionals
PBE0	–75	–86	11
B3LYP	–24	–19	–5
GGA+vdW functionals
PBE-D2	–262(−276)	–331(−340)	69(64)
PBE-D3	–248	–314	66
PBE+vdW(TS)	–271(−273)	–356(−362)	85(89)
vdW-DF	-183(−193)	–297(−303)	115(110)
vdW-DF2	–210(−214)	-292(−295)	82(81)
optPBE-vdW	–248(−252)	–369(−372)	121(121)
optB86b-vdW	–236(−238)	–350(−358)	115(119)
hybrid+vdW functionals
PBE0-D3	–241	–311	69
B3LYP-D3	–273	–345	71

The best performing functionals are indicated in bold. All energy values are in meV. Energies have been obtained on PBE-D3 optimized structures, but in parentheses we also report the adsorption energies when geometries are fully relaxed consistently for each GGA and GGA+vdW functional.

As on the hydroxyl-terminated face, RPBE gives an adsorption energy that is noticeably less exothermic than PBE. In the case of the vdW-DFs we see an across-the-board stabilization relative to the GGA functionals. Like at the hydroxyl-terminated face, vdW-DF and vdW-DF2 give the weakest adsorption energy of the vdW-DFs. The PBE-D2 and PBE+vdW(TS) functionals give the strongest overall adsorption energy, with 276 and 273 meV, respectively, followed by optPBE-vdW, PBE-D3, and optB86b-vdW with values close to 250 meV. Overall, the spread of the vdW-based evaluations is much smaller than for the hydroxyl-terminated face. Whereas at the hydroxyl-terminated face the adsorption structure was altered only moderately upon inclusion of vdW, at the silicate-terminated face more significant changes are observed. Specifically, the GGA functionals predict the molecule to be much further away from the surface than the vdW inclusive functionals do. This difference is also reflected in Figure , if we consider the difference between the adsorption energies of the GGA functionals at the PBE-D3 geometry and when the structures are relaxed. On the other hand, the geometries provided by the vdW-inclusive approaches are in very good agreement with PBE-D3, so Eads evaluated on either the PBE-D3 geometry or on the relaxed structures are similar.

Comparison with DMC supports the general reliability of the vdW-corrected and vdW-inclusive approaches over the bare GGA and hybrid functionals. However, it also shows that almost all the GGA+vdW and the two hybrid+vdW functionals overestimate the adsorption energy. Similar overestimates have been seen recently for physisorbed water on hexagonal boron-nitride.[33] On this surface the best performance is seen for the vdW-DF and the vdW-DF2 functionals, both of them being in agreement with DMC, given the DMC stochastic error. It is also interesting to note that even though water still binds preferentially to the hydroxyl-terminated face the relative adsorption strengths are significantly altered: the ratio EadsH/EadsH is 6.6 for PBE, 9.3 for RPBE, ca. 3 for the vdW-inclusive functionals, and 3.5 ± 0.5 at the DMC level.

Methanol Adsorption on the Silicate-Terminated Face of Kaolinite

Similar to water, we again expect dispersion interactions to play more of a role at the silicate-terminated than at the hydroxyl-terminated face. Indeed, we again observe big differences between functionals including or not the vdW interaction. In the most stable structure found using the PBE-D3 functional there is no hydrogen-bond-like interaction, and methanol is parallel to the surface of the slab (see Figure D). The geometry at the GGA level has the methanol molecule found at a much larger distance from the surface, as depicted in Figure D.

As on the hydroxyl-terminated face, the degree of stabilization due to dispersion interactions is greater for methanol than it is for water. At the GGA level, water and methanol bind with similar interaction strengths (methanol binds more strongly by 27 meV at the PBE level), but when vdW is accounted for in the functional we observe that methanol binds more strongly by 64 meV (for PBE-D2) to 121 meV (for optPBE-vdW), as shown in Figure .

The comparison with DMC shows that in this case the GGA functionals underestimate the adsorption energy and that all the GGA+vdW and hybrid+vdW functionals overestimate Eads. The best agreement is again obtained for the vdW-DF and vdW-DF2 functionals, the former overestimating the interaction by 47 ± 18 meV and the latter by 42 ± 18 meV. The ratio EadsMeOH@AlOH/EadsMeOH@SiO is 5.2 for PBE, 2.1 for RPBE, between 2.1 and 2.6 for the vdW-corrected and vdW-inclusive functionals, and 2.8 ± 0.2 at the DMC level.

Discussion and Conclusions

In this paper we have used QMC to examine the adsorption of water and methanol on the hydroxyl- and silicate-terminated (001) faces of kaolinite. The QMC results on the hydroxyl-terminated face have been obtained independently with two different fixed-node projection QMC methods: DMC and LRDMC. The two methods differ in terms of algorithms (DMC is based on a time-discretization approximation, LRDMC on a space-discretization approximation), implementation (DMC calculations have been performed using the CASINO code; LRDMC using the TurboRVB code), and setup (for instance, different PPs, basis sets, and Jastrow terms). Nonetheless, both approaches produce results in good agreement with the small differences between the approaches coming within the stochastic error of the evaluations.

QMC results indicate that both water and methanol adsorb on the hydroxyl-terminated face, forming three H-bonds, with an interaction energy larger than 0.6 eV. The adsorption on the silicate-terminated face is much weaker, smaller than 0.3 eV. In both cases the methanol binds slightly more strongly than water.

As discussed, the QMC results provide a benchmark that can help in further understanding of how other computationally cheaper methods perform for adsorption. Specifically, we have compared them with the results provided by a selection of commonly used XC functionals in DFT (covering GGA, hybrid, vdW-corrected GGA, vdW-corrected hybrid, and vdW-inclusive functionals). This shows that the vdW-corrected and vdW-inclusive functionals predict adsorption energies that are considerably larger than those calculated using the bare GGA or hybrid functionals, but the degree of stabilization is system dependent. As discussed, in the systems under consideration in this work the QMC references indicate that bare GGA and hybrid-based predictions are often underestimated, whereas approaches that account for the vdW interaction yield results in qualitative agreement with QMC, although the absolute value of the adsorption energy can be overestimated, particularly on the silicate-terminated face. Overall, the best results are provided by vdW-DF2, and among the vdW-corrected approaches we notice good performance from PBE-D3. Inclusion of exact exchange does not appear to lead to any improvement for the systems considered here; for instance, results from PBE-D3 and PBE0-D3 are almost identical. The GGA+vdW functionals, although based on GGA, perform better than the bare GGAs also in terms of geometries. Indeed, on the silicate-terminated face (where the interaction is weaker) structure relaxation performed with the vdW-corrected and vdW-inclusive functionals leads to very similar configurations, whereas with the GGAs the adsorbates sometimes strayed away from the surface. Looking forward there certainly still seems to be scope for further improvements in the treatment of these systems with DFT. Of the functional considered vdW-DF2 offers the best performance, but it does not convincingly deliver chemical accuracy for all four adsorption scenarios considered. Approaches such as Hamada’s revised vdW-DF2 functional[105] or Tkatchenko’s many-body dispersion[106] would be interesting to explore.

The comparison between the adsorption of water and methanol is also interesting. At the GGA level there is very little difference in adsorption energies, whereas methanol becomes more strongly bound when vdW interactions are accounted for. As clay minerals can cleanse groundwater through the uptake of pollutants, the relative adsorption energies with respect to water is a highly important quantity. Even for methanol, which is one of the simplest organic molecules able to form a hydrogen bond, we see that including vdW interactions can significantly alter the adsorption energy relative to water; on the hydroxyl-terminated face, water and methanol bind with similar energies, but inclusion of dispersion forces tips the balance in favor of methanol.

Before closing we note that we have examined the consequences of altering the relative interaction strength of water and methanol with kaolinite through a series of classical molecular dynamics simulations of liquid water–methanol solutions on kaolinite. The results of these simulations are shown in Figures A–6B. Specifically in Figure B we show results obtained with water and methanol interaction parameters that use standard Lorentz–Berthelot mixing rules, values matching PBE, or values matching DMC. As can be seen Figure B, with the DMC value of ΔEads the adsorption of methanol yields a density profile with a much more pronounced first peak compared to the ΔEads, corresponding to PBE or standard Lorentz–Berthelot mixing rules. Thus, we see that the standard approach for exploring aqueous solutions at a clay surface with force fields leads to a rather poor description of the interface. This effect is likely to become even more significant as the size of the organic tail of the adsorbate increases and demonstrates the importance of an accurate modeling of dispersion interactions when exploring wet interfaces of environmental relevance. The ability to accurately incorporate nonlocal dispersion interactions is therefore extremely important if one aims to model environmentally relevant adsorption processes on kaolinite and other clays.

Figure 6

Molecular dynamics simulations of a water–methanol mixture on the hydroxyl-terminated face of kaolinite. Panel A: Snapshot of the simulation system, where the methanol molecules are shown in red and the water molecules as glassy blue circles, and the kaolinite slab is the same as in Figure . The black lines show part of the periodic simulation cell boundaries. Panel B: Density profiles of methanol above the kaolinite surface, for the eight replicas with temperatures ranging from 275 to 310 K, obtained for values of ΔEads (namely, the difference between adsorption energy of water and methanol (see text)) corresponding to Lorentz–Bethelot, PBE, and DMC. “z” is the distance from the average position of the top layer of oxygen atoms in the kaolinite surface. It is clear that an appropriate choice of ΔEads can significantly affect the density of the water–methanol solution at the interface.