Sensitivity Analysis of Cluster Models for Calculating Adsorption

ARTICLE pubs.acs.org/JPCC

Sensitivity Analysis of Cluster Models for Calculating Adsorption Energies for Organic Molecules on Mineral Surfaces M. P. Andersson* and S. L. S. Stipp Nano-Science Center, Department of Chemistry, University of Copenhagen, Denmark ABSTRACT: We calculated the adsorption energy for ethanol on the magnesite {10.4} surface using density functional theory (DFT) and cluster models for the mineral surface, and we quantified the errors introduced by using a finite size cluster, freezing various parts of the mineral cluster during geometry optimization and for altering the edges of the cluster. We also investigated how the adsorption energy changes for increasingly accurate density functionals, PBE, BLYP, B3LYP, and B2PLYP, also when supplemented with empirical dispersion (-D). We concluded that calculations with clusters large enough to include the surface atoms and groups binding to the adsorbate and their nearest-neighbor ions provide accurate adsorption energies, typically for MgCO3
ethanol systems, which is about 60 atoms. Using B3LYP-D and a finite size cluster of 80 atoms, we found that the adsorption energy was underestimated by 0.17 eV for adsorption from vacuum and by 0.10 eV for adsorption from solution, and we estimated a random error of 0.10 eV in adsorption energy when surface atoms were frozen during geometry optimization. The basis set superposition error leads to an overestimation of the adsorption energy of 0.08 eV, which in solution almost cancels the effects of finite size. We also compared adsorption energies for binding to magnesite clusters of ethanol, water, acetic acid, hydroxyl, and acetate from vacuum and from aqueous solution. When adsorbed from solution, using an implicit solvent model, bonding with water is stronger than for the other molecules. No adsorption of any of the other molecules was predicted. The adsorption energy for the acetate ion is considerably more uncertain than for the other molecules. We recommend B3LYP-D as a reliable and computationally effective method for calculating adsorption energies of organic molecules on mineral surfaces using cluster models.

’ INTRODUCTION The interaction between organic molecules and mineral surfaces is central to quite diverse fields, such as enhanced oil recovery, pollutant migration through soil, and biomineralization. Perhaps the most important parameter describing these interactions is the adsorption energy of an organic molecule on a mineral surface. Under normal conditions, this is the dominant term in the free energy of adsorption, which can be used to derive equilibrium constants. By comparing adsorption energy for several molecules, we can determine which species are most likely to adsorb to a mineral surface under a given set of conditions. Adsorbed species can significantly influence chemical and physical properties of a surface. The adsorption energy for a molecule onto a mineral surface can be calculated using several methods. One approach uses periodic boundary conditions combined with density functional theory or empirical force fields describing each of the atomic interactions in the system. Another approach simplifies the system, modeling the surface as a cluster “cut” from it. This nonperiodic approach introduces errors from the treatment of the surface, but it provides a possibility to use more accurate quantum chemistry computational methods. For metals and semiconductors, it is more natural to use periodic boundary conditions, considering that the electronic states themselves are highly delocalized. For insulators, which is the case for most minerals, the cluster approach is quite frequently used. It assumes that the electronic states of the r 2011 American Chemical Society

substrate are localized enough that any edge effects of the cluster do not severely influence the properties of the adsorption site. Recently, several authors have suggested combining periodic density functional theory calculations with highly accurate ab initio calculations for small local clusters to obtain accurate adsorption energies, e.g., refs 1 and 2. Each method has its merits and flaws, and each has its sources of error. For periodic boundary conditions, there are no edge effects, and periodicity in the plane parallel to the surface is modeled accurately. This is true provided no local defects or atomic substitutions are studied because these would be made artificially periodic. Quantum chemistry programs with periodic boundary conditions use either a plane-wave basis set and pure density functionals (e.g., ref 3), or they employ localized basis sets and either pure or hybrid density functionals (e.g., ref 4). The major source of error in these kinds of calculations is thus from the quantum chemical method itself. Perhaps the most widely used density functional for adsorption studies is the PBE (Perdew Burke Ernzerhof) functional.5 It has been estimated from Bayesian statistics that small variations in this functional significantly change the adsorption energy calculated for CO on Cu(100); the error estimate is about 0.3 eV.6 Comparing results Received: December 21, 2010 Revised: March 12, 2011 Published: May 04, 2011 10044

dx.doi.org/10.1021/jp1121189 | J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C derived from a completely different functional, such as BLYP (exchange from Becke7 and correlation from Lee, Yang, and Parr8), showed that this uncertainty is in fact higher. Force-field methods are inherently empirical in their construction. Force fields that take into account the adsorption behavior of molecules on surfaces are parametrized using density functional theory and thus suffer from the same error (e.g., ref 9). Whenever a cluster is used to represent a surface, error arises because the cluster is artificially cut from the solid. Clusters representing covalently bonded minerals, e.g., clays, are not compensated for dangling bonds, and clusters representing ionic compounds such as magnesite lack long-range electrostatic effects. With this simplification, however, it is possible to use more sophisticated quantum chemistry methods, such as hybrid density functionals, double hybrid density functionals, or, in principle, even higher-order ab initio methods, such as coupled cluster methods,1,2 There is a general lack of results quantifying the error that arises from using local cluster models. One recent investigation, for water adsorption on NaCl,1 tested differences for clusters of up to 50 atoms, with and without embedding the cluster in point charges. They found that when electronic embedding is included quite small clusters are required to calculate accurate adsorption energies. The situation is quite different for metals, where quite a few studies have been conducted for adsorption of primarily diatomic molecules such as carbon monoxide on metal clusters.10 For metal clusters, edge and size effects can totally dominate adsorption properties, and chemisorption energies are known to vary significantly as a function of cluster size. For clusters of a certain size, the adsorbate can even be repelled instead of adsorbed, and one might need to choose clusters with particular properties in order for the adsorption energy to be calculated reliably.11 Our aim in this paper is to quantify the errors introduced by a finite sized, nonembedded cluster in a system containing ethanol adsorbing on a magnesite {10.4} surface. Ethanol, CH3CH2OH, is a small organic compound with an OH functional group. Ethanol has been shown to influence the crystallization of minerals.12
14 Magnesite, MgCO3, is a carbonate mineral with a rhombohedral structure, the same as calcite, CaCO3, and several other divalent metal carbonate minerals. Carbonate minerals are relatively common in nature and have been studied widely, particularly calcite and, in particular, their stable {10.4} cleavage face (e.g., refs 15
19). The energy preference and stability of this plane is demonstrated by the clean, sharp, rhombohedral surfaces that appear on crystals grown in nature and in the laboratory (e.g., ref 12, Figure 1a). For modeling, we can construct clusters whose sides consist only of this surface termination. These surfaces are less reactive than clusters created with other surface terminations. An important point in favor of the cluster approach is that periodic boundary conditions make investigation of ions much more difficult because long-range electrostatic interactions need special conditions, or very large unit cells. Thus, the cluster approach is suitable for investigating neutral as well as ionic species, implicit solvation, and empirical dispersion,20 all of which we want to incorporate in our modeling and which we therefore investigate in this paper. The primary goal of this study is to quantify the computational errors that are introduced by using cluster models to describe adsorption of organic molecules on a mineral surface. We have selected the ethanol
magnesite {10.4} system and test for convergence of adsorption energy as a function of the finite size, charge, relaxation, and edge effects of the clusters. Basis set convergence is also tested. Thus, we quantify the convergence of

ARTICLE

Figure 1. Optimized geometry of (a) ethanol and (b) acetate adsorbed on the 9  9  3 Å magnesite cluster with 80 atoms and frozen geometry. The rhombohedral form of the cluster is outlined with blue. The bonds between the adsorbates and the surface are shown in black. Atoms that are removed from the corners for investigating the edge and charge effects are shown by arrows. For the investigations of the effect with higher charge (þ4 and
4), atoms were removed from adjacent corners in the top layer as well. The atoms are colored: Mg, green; C, cyan; O, red; and H, white.

our cluster calculations with respect to size, relaxation, and basis set in much the same way that has been done routinely for pseudopotential plane-wave calculations. Our results are intended to provide computational chemists and physicists with an error bar for adsorption studies using cluster models. Once the finite size and basis set errors were quantified, we performed a comparative adsorption study of several small molecules binding through oxygen to a magnesite {10.4} surface. In addition to ethanol, we tested the neutral species: water, H2O, and acetic acid, CH3COOH, as well as their deprotonated counterparts: hydroxyl and acetate. For these five entities, we compared not only gas phase adsorption energies but also adsorption energies in aqueous solution using the implicit solvation model COSMO.21 In an aqueous solution, adsorption energy is defined analogously, but all molecules and clusters are placed in an implicit solvent instead of a vacuum. Implicit dielectric continuum models, while rather crude, give reasonable results for water solvation. 10045

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C Finally, we compared the performance of the PBE functional and the BLYP family of functionals: BLYP, B3LYP,22 and B2PLYP.23 The B2PLYP is a fairly new type of density functional that mixes not only Hartree
Fock and DFT as regular hybrid density functionals but also some ab initio correlation energy from Møller
Plesset perturbation theory (MP2). The B2PLYP functional has been found to give excellent agreement with converged, coupled cluster results for adsorption of hydrocarbons on zeolites.24 It is also in excellent agreement with experimental atomization energies,25 reaction barrier heights,26 and vibrational frequencies.27 It outperforms the B3LYP functional in almost all respects. The BLYP and B3LYP functionals are known to provide good adsorption energies for molecules on solid surfaces, also for metals,28 even though B3LYP is known to fail for properties of the metal substrate.29 We also investigated the effect of empirical dispersion,30
32 which has been quite successful at compensating for known discrepancies in density functional theory. This comparison is relevant in light of the success of empirical dispersion and because the availability of double hybrid functionals is a feature unique to nonperiodic quantum chemistry codes.

’ COMPUTATIONAL DETAILS Software. The magnesite clusters were generated in Materials Studio, and experimental geometries have been used for all frozen clusters. All geometry optimizations have been performed with TeraChem v1.0 and v1.05,33 quantum chemistry programs that run on graphics processor units. This approach allowed us to study clusters of up to 720 atoms using a conventional cubic scaling code. Single-point calculations with triple-ζ basis sets and polarization functions were performed using the ORCA program,34 version 2.7. Implicit solvation calculations using the COSMO model also used the ORCA program. Basis Sets. Both ORCA and TeraChem use atomic orbital type basis sets, whose radial components are expressed as a sum of Gaussian functions. Chemical bonding is mostly through the valence electrons. Therefore, two or more orbitals are used to represent the electron density in the valence region, so-called split-valence basis sets. We have used both Pople type basis sets as well as the more modern Alrich type. The Pople type are denoted as X-YZ(W)G, where G stands for Gaussian functions and X represents the number of Gaussian functions for every orbital in the core region, and Y and Z (and possibly W) are the number of Gaussian functions for the orbitals in the valence region. If only YZ are present, the basis set has two valence orbitals and is called a double-ζ basis set, e.g., 3-21G. If YZW are present, the basis set has three valence orbitals and is referred to as a triple-ζ basis set, e.g., 6-311G. The small basis set 3-21G was used to investigate finite size, charge, cluster termination,and relaxation effects on the adsorption energy for ethanol on the magnesite {10.4} surface. Although not very accurate, it allows for study of much larger clusters. The 6-311G basis set was used to calculate adsorption geometries for the comparative adsorption study and relaxation effects. Polarization functions are not available in TeraChem v1.0 and v1.05, so we used the default basis sets (defbas-X, X = 3
6) in ORCA, which are variants of the TZV basis set.35,36 Defbas-5 and defbas-6 include Pople type diffuse functions,and all of them include polarization functions. Defbas-3 has a single d polarization function for non-hydrogen atoms, and each step up the default basis set size adds more polarization functions for all atoms.

ARTICLE

Quantum Chemistry. All geometry optimizations were performed with the B3LYP functional, and several starting configurations were used for all molecules. B3LYP is known to provide good molecular geometries, even for hydrogen bonded complexes (e.g., ref 37) as well as accurate lattice constants for minerals (e.g., ref 38). Single-point energies on B3LYP/6-311G geometries for the adsorbates were calculated with the pure density functionals PBE and BLYP, the hybrid density functional B3LYP, and the double hybrid functional B2PLYP, with and without empirical dispersion. We have added -D to the functional name for calculations with dispersion: PBE-D, BLYP-D, B3LYPD, and B2PLYP-D in tables and figures. All of these single-point energies use at least a defbas-4 basis set, i.e., TZV(p) for hydrogens and TZV(2d) for main group elements. All geometry optimizations used the following convergence criteria (atomic units): rms gradient < 3e
4; maximum gradient < 4.5e
4; rms step < 1.2e
3; maximum step < 1.8e
3; and maximum change in total energy between two consecutive steps < 1e
5. The ORCA PBE and BLYP single-point energies used the resolution of identity (RI) approximation. B3LYP singlepoint energies used ORCA and the very efficient RIJCOSX39 approximation, which evaluates part of the two-electron integral numerically and part analytically. For situations where all atoms had diffuse basis functions, we were unable to reach wave function convergence with the RIJCOSX approximation. In these cases, we used the RIJONX approximation, which evaluates the coulomb term using the RI approximation. All ORCA B2PLYP single-point energies used the RIJCOSX for the exchange, as well as the RI approximation for the MP2 correlation. Implicit Solvation. We used the COSMO model for water with ORCA default values for C, O, and H. For Mg, we used a cavity radius of 2.02 Å, which is 1.17 times the van der Waals radius 1.73 Å.40 This is analogous to how the default values for the more common elements (such as C, O, and H) have been determined for use in ORCA. No geometry optimization with implicit solvation has been used in this paper, only single-point energies on gas phase geometries. Molecular geometries tend to change only slightly when going from the gas phase to implicit solvent. We have verified this approach for the most crucial part of the geometry, the adsorbate
surface bonding region. For ethanol adsorption on a frozen 80 atom magnesite cluster, we used ORCA to calculate the B3LYP/6-311G equilibrium geometry with and without implicit solvent. The O
Mg bond length connecting the magnesite cluster and the ethanol differed by 0.01 Å and O
Mg
O bond angles, and the direction of the adsorption bond differed by about 1°. These small differences close to a stationary point give only small uncertainties in energies. We found that the energy difference between optimizing in implicit solvent and a single-point energy on the vacuum geometry is only 0.01 eV. For implicit solvation calculations, we also compared the adsorption energies of all adsorbates internally with adsorbed water as reference. The molecules we investigated bind similarly to the magnesite, and most of the errors from using gas phase optimized structures are expected to cancel. All energies from COSMO calculations are corrected for outlying charge. Magnesite Clusters. We created clusters, two molecular layers thick, where lateral size ranged from 60 atoms, with 12 units of MgCO3, to 720 atoms, with 144 units of MgCO3. We also constructed clusters with a 9  9 Å lateral dimension that were three and four atomic layers thick. The 9  9  3 Å, 80 atom cluster is shown in Figure 1, with lines drawn around the 10046

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 1. Optimized Geometry for the Ethanol
Magnesite System for the 3-21G and 6-311G Basis Sets Using Frozen Clusters and with the Binding Surface Atoms and Groups Relaxeda 3-21G frozen cluster

6-311G frozen cluster

6-311G 5 cluster atoms relaxed

O
Mg bond length (Å)

2.060

2.070

2.108

2.087

O
H bond length (Å) H
CO3 bond length (Å)

1.004 1.825

1.003 1.796

0.980 1.924

0.979 1.901

Mg
O
H angle (degrees) a

3-21G 5 cluster atoms relaxed

107.6

107.8

109.2

109.9

Data are from the 6  9  3 Å cluster (60 atoms).

Table 2a. Effect of Lateral Cluster Size on Adsorption Energy for Ethanol on Magnesite {10.4}a magnesite cluster size: length along

number of atoms

B3LYP adsorption

maximum component for gradient

B3LYP-D adsorption

rhombohedral edges (Å  Å  Å)

in cluster

energy (eV)

on the ethanol atoms (Ha/Bohr)

energy (eV)

693

60


2.09

1.6e
4

993 15  15  3

80 180


2.10
2.13b

5.4e
5 1.8e
3

21  21  3

320


2.14b

2.4e
3


2.60

24  24  3

405


2.15b

-c


2.60

27  27  3

500


2.14b

-c


2.61

720


2.14

-c


2.61

33  33  3 maximum deviation from smallest cluster

b

0.07d


2.49
2.58

0.13d

a

The clusters with 180 atoms or more use the optimized ethanol adsorption geometry from the 80 atom cluster. Whenever computer memory was sufficient, the gradients were computed to check this approximation. b Single-point energies using the optimized adsorption geometry from the 9  9  3 Å cluster. c Insufficient memory for gradient calculations. d An extra 0.02 eV is added to the results to compensate for incomplete optimization, see text.

edges to mark its rhombohedral form. In addition, the bonds between the adsorbate, ethanol, and the surface are indicated by thick black lines. The atoms that were removed to investigate edge and charge effects (shown by the arrows) were taken from the corners, far away from the adsorption site. Water Boxes. We generated four water boxes with a range of sizes, using the solvate module with the program Visual Molecular Dynamics (VMD).41 The water boxes ranged from 84 atoms to 627 atoms. Adsorption Energy. The adsorption energy, Eads, of a molecule A on cluster C is calculated in the conventional way as Eads ðAÞ ¼ EðA þ CÞ
EðAÞ
EðCÞ where E(A þ C) represents the energy of the cluster with the molecule adsorbed, and E(A) and E(C) are the energies of the free molecule and the bare cluster. This means that greater negative adsorption energy implies stronger adsorption. For the implicit solvent calculations, we assume that the contributions to the solvation energy from the parts of the cluster far away from the adsorption site cancel, much the same as for gas phase adsorption. To make predictions, such as if ethanol or water covers the magnesite surface, we need to consider the reaction EtOHðaqÞ þ H2 OðadsÞ T EtOHðadsÞ þ H2 OðaqÞ where ethanol, EtOH, in solution replaces a water adsorbed on the surface. The reaction energy is ΔE ¼ Eads ðEtOHÞ
Eads ðH2 OÞ where Eads now denotes the adsorption energy of a molecule on the surface in an implicit solvent. Thus, the relevant parameter for obtaining equilibrium constants is the adsorption energy of ethanol (in implicit water) relative to the adsorption energy of water (in implicit water).

’ RESULTS AND DISCUSSION Geometry. We present some key geometric parameters for ethanol adsorption on the 60 atom cluster in Table 1. Bond distances and angles changed considerably when accuracy increased from the 3-21G basis set to 6-311G, as is expected. Slight changes are also observed when the MgCO3 unit in contact with the adsorbing molecule is relaxed. The change is largest in the bonds between the adsorbate and the surface atoms. Lateral Cluster Size. We now turn to testing sensitivity of the calculations to lateral cluster size. In Table 2a, we see that convergence is reached for 320 atom clusters, and the adsorption energy increases slowly and monotonically as the cluster lateral size increases. For clusters of 180 atoms or more, we used the adsorption geometry from the frozen 80 atom cluster. Whenever enough computer memory was available, we calculated the gradients to estimate the uncertainty from using this procedure. The maximum component of the gradient of the 180 and 320 atom clusters was on the order of 1e
3 Ha/Bohr. In the geometry optimizations for ethanol on the frozen 60 and 80 atom clusters, a further 0.02 eV was gained when the maximum component of the gradient was of this order during the optimization. Thus, we add 0.02 eV to the total energy difference between the largest and smallest clusters in the lateral size investigation. The error from using a smaller cluster is quite small, only 0.07 eV (0.05 eV þ 0.02 eV). Therefore, we conclude that the 60 and 80 atom clusters are accurate enough for our purposes because with cubic scaling codes there is a saving of almost 2 orders of magnitude in computational time compared to using the 320 atom cluster. For adsorption in vacuum, various embedding schemes for point charges could reduce this error further,1 but point charge embedding is not compatible with implicit solvation models. A combined quantum mechanics/molecular mechanics approach (QM/MM) would also decrease this error and take into account 10047

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 2b. Comparing the Effect of Lateral Cluster Size on Adsorption Energy for Ethanol, Water, and Acetic Acid on Magnesite {10.4}a magnesite cluster size: length along rhombohedral

number of atoms

B3LYP-D adsorption

B3LYP-D adsorption

B3LYP-D adsorption

edges (Å  Å  Å)

in cluster

energy (eV)

energy (eV)

energy (eV)

ethanol

water

acetic acid

993 15  15  3

80


2.48


2.50


2.49

180


2.58


2.50


2.53

21  21  3

320


2.60


2.50


2.54

27  27  3

500


2.61


2.50


2.55

0.13

0.00

0.06

maximum deviation from smallest cluster a

The clusters with 180 atoms or more use optimized adsorption geometries from the 80 atom clusters.

Table 3. Effect of Vertical Cluster Size on Adsorption Energy for Ethanol on Magnesite {10.4} magnesite cluster size: length along rhombohedral

number of atoms

number of layers perpendicular

B3LYP adsorption

B3LYP-D adsorption

edges (Å  Å  Å)

in a cluster

to adsorption surface

energy (eV)

energy (eV)

993

80

2


2.09


2.48

996

120

3


2.08


2.52

999 maximum deviation from smallest cluster

160

4


2.09 0.01


2.52 0.04

solvent effects but is more computationally demanding and out of the scope of our study. We also performed the same calculations with empirical dispersion included. This raises the energy difference to 0.13 instead of 0.07. Adsorption energies are thus even more underestimated (by 0.07 eV) by neglecting dispersion for large surfaces. Because of the significant difference in results with and without dispersion, we made two error estimates, one with dispersion and one without. The contribution from dispersion is also slightly overestimated. The dispersion parameters used by ORCA for Mg are derived for atomic Mg, not ionic Mg2þ. In ionic compounds the dispersion parameters for Mg are much closer to those of Ne2. To check the transferability of the finite size error to other molecules, we also performed calculations for water and acetic acid on clusters with increasing lateral size, analogous to the ethanol calculations, using the B3LYP-D method. The results are presented in Table 2b. The expansion from an 80 atom cluster up to 500 atoms influences the adsorption energy for the adsorbate molecules quite differently. The adsorption energy for water converges after 80 atoms, and for acetic acid, the change in adsorption energy lies somewhere between water and ethanol. This tells us that the adsorption energy is underestimated with as much as 0.13 eV. Vertical Cluster Size. Compared with the lateral size effect, the number of layers makes a smaller difference. In Table 3, it is clear that there is virtually no effect from including more layers than two in the calculations. The energy difference is only 0.01 eV. We conclude that two layers are quite enough to get converged adsorption energies with respect to cluster size in the direction perpendicular to the surface. If dispersion correction is included, the energy difference is only slightly higher, 0.04 eV, but convergence is reached at three layers. The adsorption energy of water on the similar calcite {10.4} surface was also reached at three layers, within 0.03 eV in a recent study.42 Dispersion. The finite size errors are thus quantified and shown to be small, at least when dispersion is not included. These errors tend to underestimate the adsorption energy. For disper-

sion corrected calculations, the adsorption energy of ethanol is underestimated by 0.17 eV, and for nondispersion corrected calculations, it is 0.08 eV. The effect on water and acetic acid is smaller, and we therefore take the ethanol results as an upper error limit. Because the dispersion corrected functionals behave better for long-range interactions than the uncorrected ones, the larger error estimate is valid for adsorption from the gas phase. However, we are also interested in adsorption from an aqueous phase using implicit solvent, for which no dispersion effects between the adsorbing molecule and the liquid phase are taken into account. We estimate the difference in the two by calculating the adsorption energy of a water molecule on the side of a water box as a function of box size. This procedure is meant to be analogous to estimating finite size effects of the cluster, but now for the solvent part next to the solid/water interface. We have generated water boxes using VMD and used them without geometry optimization and removed one water molecule from the center of a side of the box. We have then performed single-point calculations on these structures with and without dispersion using TeraChem/3-21G, the same method used for cluster finite size error estimation. In Table 4, we present how the inclusion of dispersion affects the adsorption energy as a function of size. The change varies smoothly as the size is increased: from the 15  15  5 Å, 84 atom water box, to the 40  40  5 Å, 627 atom water box, the change is
0.06 eV. This can be compared for the analogous change from adsorption onto the magnesite clusters, comparing the 9  9  3 Å, 80 atom cluster with the 33  33  3 Å, 720 atom cluster, for which the same calculation gives
0.08 eV. Therefore, a large part of the dispersion contribution to the error is canceled. For adsorption from a solvent, we then use the result from the nondispersion calculation plus the resulting difference, i.e., 0.10 eV. All of the dispersion effects do not cancel because there is a higher proportion of heavier elements in minerals than in water. Relaxation Effects. There is a rather large spread in adsorption energies for ethanol on magnesite for the small, 3-21G basis 10048

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 4. Dispersion Effect on Adsorption Energy for a Water Molecule on the Side of a Water Box, Where Box Size Variesa size of water cluster (Å)

a

number of atoms in the

difference between dispersion corrected and uncorrected adsorption energies

water cluster

for water on side of the box, compared with the 84 atom water box [eV]

15  15  5

84

-

15  15  15

297


0.02

20  20  12

402


0.04

40  40  5

627


0.06

The results are from single-point TeraChem/B3LYP/3-21G calculations on VMD generated water boxes.

Table 5a. Relaxation Effects on the Adsorption Energy for Ethanol on Various Magnesite Clusters for the Two Basis sets 3-21G and 6-311G magnesite cluster size: length along rhombohedral edges (Å  Å  Å)

number of atoms number of atoms in the cluster

description of

relaxed

693

60

0

693

60

5

993 696

80 90 150

adsorption energy adsorption energy

relaxed atoms

3-21G (eV)

6-311G (eV)


2.09


1.16

surface groups bound to the adsorbate


1.92


1.09

5 10

surface groups bound to the adsorbate surface groups bound to the adsorbate þ


1.96


1.07
1.16

75

whole surface layer


2.00a


1.16a

whole cluster


1.97

the groups in the layer below 12  15  3 12  15  3

150

150

a

max deviation from unrelaxed surface

0.17

0.09

a

For this size of lateral cluster, we can estimate from Table 2 that the adsorption energy should be lower by about 0.04 eV in response to the larger cluster size. Raising the adsorption energy by 0.04 eV, to compensate for this effect, does not change the maximum deviation from the unrelaxed surface.

Table 5b. Comparison of Relaxation Effects on the Adsorption Energy for Ethanol, Water, Acetic Acid, Acetate, and Hydroxyl on Magnesite Clusters with 80 Atoms for the Basis Set 6-311G number of atoms relaxed

adsorption energy adsorption energy adsorption energy adsorption energy adsorption energy description of relaxed atoms

(eV)

(eV)

(eV)

(eV)

(eV)

ethanol

water

acetic acid

acetate

hydroxyl


2.37


4.31


1.16


1.31


1.20

5

surface groups bound to the adsorbate


1.07


1.22


1.11

10

center atoms/groups in top layer


1.10


1.22


1.13


2.99


4.56

0.09

0.09

0.09

0.62

0.28

0

max. deviation from unrelaxed surface max. dev./adsorption

8%

7%

8%


4.59

26%

6%

energy

set, namely, 0.17 eV. The errors for the more accurate 6-311G basis set are at most only half as large, 0.09 eV (Table 5a). For the comparative adsorption study, which is described later in the paper, we used triple-ζ basis sets, so this is the error to keep in mind. The time required for geometry optimization increases as more geometric parameters are included. This means that freezing more atoms saves more computing time, and our results show that relaxing more atoms does not provide significantly more accurate results for the adsorbates binding to the surface through a single bond. The uncertainty from various degrees of relaxation is smaller than effects from finite size, and the extra computational effort can better be spent elsewhere. Thus, we conclude that the error from the various degrees of surface relaxation for accurate basis sets is 0.09 eV for ethanol. In Table 5b, we compare relaxation effects for adsorption of ethanol, water, acetic acid, and the ions acetate and hydroxyl on clusters with 80 atoms. Contrary to the finite size effect, the

change in adsorption energy for the three neutral molecules is very similar when the surface atoms are relaxed. This means that later in the paper, when we use adsorption energies for the other adsorbing molecules relative to adsorbed water, the relaxational effects cancel. The charged species behave differently. The change in adsorption energy for the hydroxyl is opposite in sign to that of the neutral species, but still rather small. The relative error compared to the adsorption energy is very similar for all species binding to a single surface Mg atom. The acetate ion binds in a different manner, each of its oxygen atoms binding to a surface Mg atom (Figure 1b). Perhaps as a consequence of this, the adsorption energy is much more sensitive to relaxation effects. Using frozen cluster geometry thus considerably reduces the accuracy of acetate bonding. It might be a general phenomenon that molecules bond to a surface with two or more bonds instead of one, but such a generalization is outside the scope of this paper. 10049

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 6. Effect of Cluster Charge and Edge Geometry on the Adsorption Energy for Ethanol, on a 9  9  3 Å (80 Atom) Magnesite Clustera ions removed from

resulting cluster

adsorption energy

corners of cluster

charge (e)

(eV)


4


1.98


2


2.04

0


2.09

Mg2þ, CO32
CO32


0 þ2


2.06
2.15

2CO32


þ4


2.45

2Mg

2þ

Mg2þ none

a

All Mg2þ and CO32
ions were removed as far as possible from the adsorption site; i.e., only next-nearest neighbors were removed (shown in Figure 1).

Edge Effects. We have also constructed a series of 9  9  3 Å

clusters from which we removed atoms from the corners of the cluster, to investigate the effect of cluster charge and termination on the adsorption energy. In Table 6, we compare the central two rows for neutral clusters to find that the adsorption energy differs only by 0.03 eV, when one Mg2þ and one CO32
are removed from two different corners, neither of which is adjacent to the adsorption size. The atoms removed are indicated by arrows in Figure 1a. Thus, the details of the surface termination are secondary to relaxation effects and the lateral size of the cluster. This is a very limited test, and larger edge modifications might result in larger variations. We anticipate that as long as all nearest neighbors of the surface atoms or groups involved in bonding to the adsorbate are included the effects of changes beyond that are small. Cluster Charge. For charged clusters, when the total charge is small (þ2 e or
2 e) the adsorption energy changes by about 0.05 eV, compared with the change for the neutral cluster (Table 6). This is again smaller than the lateral size and relaxation effects. When the total charge of the cluster goes beyond þ2 e or
2 e, larger deviations in the adsorption energy are observed. On the other hand, it is up to the computational researcher to choose the charge of the cluster. It is sometimes easier to add or remove an electron from a cluster so as to work with filled electron shells, thus saving computational time and effort because the calculation is restricted rather than unrestricted. Our results indicate that such a procedure has a very small impact on adsorption energies for clusters of about 80 atoms. Some mineral surfaces do have a net charge under certain conditions. One example is clay mineral surfaces, where Al substitutes for Si in the layers and the overall surface charge varies with pH and the ions present in solution. For similar systems, this change in adsorption energy as a result of cluster charge could be important, but it is beyond the scope of this work. Estimating Total Error. We now assume that the errors from edge effects and relaxation are independent of one another and exclude any effects of charge. The standard error is then the root of the sum of square errors and amounts to 0.10 eV (10 kJ/mol), which is quite low. This error applies to all the neutral molecules in our study and can for example be compared with the Bayesian error estimate, (0.3 eV, for the adsorption energy of CO on Cu(100) using a similar class of density functionals6 and an energy difference larger than 0.3 eV for different density functionals used for water adsorption on NaCl.1 The structure of the

Figure 2. Adsorption energy for ethanol on magnesite clusters as a function of the change upon adsorption in Mulliken charge of the Mg atom binding to ethanol. The adsorption energies are from Tables 2a and 3.

various rhombohedral carbonates (magnesite, MgCO3; calcite, CaCO3; siderite, FeCO3; smithsonite, ZnCO3; rhodochrosite, MnCO3; and others) is the same, so it is likely that this error can be transferred for similar studies on these minerals. Transferability to other, noncarbonate minerals remains to be investigated. In addition to this error, there is also the underestimation for adsorption energy, i.e., as much as 0.10 eV, for adsorption from water (implicit solvent), and as much as 0.17 eV, for adsorption at a gas/solid interface. Relative Importance of Cluster Size, Edge Geometry, and Charge. To elucidate why the adsorption energy changes with increasing cluster size, we have looked at the Mulliken charges of the O atom of ethanol and the Mg atom of the surface in what forms the adsorbate
surface bond. Their sum should be a good indicator of the electron density in the bond region. It turns out that almost all of the change in Mulliken charge upon adsorption takes place on the Mg atom. Figure 2 is a plot of the adsorption energies from Tables 2a and 3 as a function of the change in Mulliken charge for that Mg atom. The good correlation demonstrates that the change in adsorption energy comes at least in part from a change in electron density in the adsorbate
surface bond. It thus shows that the bond is altered when cluster size is increased and that the change in adsorption energy for ethanol does not result from direct electrostatic interaction. It is possible that the additional atoms present in the larger clusters provide a more accurate description of the electrostatic environment in the bond region and that the self-consistent electron density is affected. One must also take into account the fact that the clusters are always expanded in a systematic way so the total charge of the added atoms is zero. While Mg is charged positively and CO3 negatively, they are added at positions with similar distances to the adsorbate region so the electrostatic interactions are canceled to a large degree. Part of the long-range electrostatic effects is likely screened out as well because of the regular arrangement of ions of alternating charge. If there is too much positive or negative charge at the cluster edges, the electrostatic effects start to dominate the interactions (Table 6). We therefore advocate that the clusters should be constructed in such a way as to allow the adsorbate to be placed close to the center of the clusters. 10050

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 7a. Effect of Basis Set on Adsorption Energy for Ethanol on the 6  9  3 Å Cluster with Contact Atoms MgCO3 Relaxeda basis set

program used

3-21G


1.92

TeraChem

6-311G/3-21G mixedb


1.47

TeraChem

6-311G


1.09

TeraChem

6-31þþG


0.94

TeraChem

default basis 3c default basis 4c


0.75
0.65

ORCAd ORCAd

default basis 5c


0.63

ORCAd

c


0.63

ORCAd

c


0.56

default basis 6

default basis 5 BSSE corrected a

adsorption energy (eV)

ORCAd b

The B3LYP functional was used in all cases. 6-311G on ethanol and the Mg and CO3 binding to it. 3-21G on all surrounding atoms. c These are ORCA standard basis sets; more details are available in the text and in the ORCA manual. d Single-point energy on the optimized 6-311G geometry from TeraChem

A most important point is that both edge and charge effects decrease for larger clusters, making the cluster size a much more important factor than either edge or charge effects. In an insulator, changes in electronic structure are localized by the “nearsightedness” of density functional theory,43,44 and for large clusters, the exponential decay44 of changes in electronic density would eventually make the exact nature of the edges irrelevant. Even for clusters as small as 80 atoms, the edge effects are small (Table 6). It is therefore quite safe to saturate dangling bonds of large clusters with added atoms or molecular entities to allow for more reliable SCF convergence or to obtain a neutral cluster. In addition, the edges of large enough clusters can be frozen to retain bulk geometry, while at the same time, enough cluster atoms can be allowed to relax to provide a more realistic environment at the adsorption site. For multiple bonds to the surface, our results for the acetate ion indicate that this can be crucial for obtaining accurate adsorption energies. The most important aspect for constructing the cluster is that it is made large enough. Our aim for this work was to provide data that would help researchers make educated decisions for constructing clusters as surface models. Basis Set. Results for the effects of the basis set are presented in Table 7a. As expected, the influence is quite large. The 3-21G basis set is inadequate for obtaining quantitatively relevant adsorption energies, but its computational efficiency enabled us to perform convergence studies for up to 720 atom clusters. To get convergence with respect to polarization functions and diffuse functions for the TZV basis set, we needed at least two sets of polarization functions for heavy atoms and one for hydrogens (the defbas4 basis set in ORCA). This gives an error of 0.02 eV compared with the largest basis set used (defbas6). This triple-ζ basis set has three sets of polarization functions on heavy atoms (2df) and two sets on hydrogens (2p) as well as two sets of diffuse functions on heavy atoms (sp) and one on hydrogens (s). Furthermore, the nearest neighbors of the atoms bound to the adsorbate need an accurate basis set. The energy difference between using 6-311G on all atoms and using 3-21G on the nearest neighbors is quite large, 0.38 eV. For the charged species, we needed to add diffuse functions at least to the oxygen atoms to reach a similar level of basis set convergence. The results are shown in Table 7b. There is a very small difference between using

diffuse functions on all atoms and only on oxygen atoms; it is less than 0.03 eV. There was a very large impact on computational time, however. When only oxygen atoms have diffuse functions, the RIJCOSX approximation could be used, which saves about a factor of 30 in computational time compared to the RIJONX approximation that we had to resort to, when all atoms had diffuse functions. In any study with localized Gaussian basis sets, interaction energies between molecules suffer from basis set superposition error (BSSE). This means that when two molecules or fragments are bound to each other each can use the basis set for the other molecule or fragment for its own electron density and gain extra (unphysical) energy. We performed counterpoise corrections45 for all molecules adsorbed on frozen 80 atom magnesite clusters using the defbas4 basis set with diffuse functions on oxygen atoms. We found BSSE in the range of 0.07
0.18 eV (Table 7b). The BSSE was about a factor of 2 higher than for the neutral species. This is quite reasonable compared with energies determined from other studies of adsorption on insulator surfaces. One example is aromatic molecules adsorbed on silica using B3LYP and a triple-ζ basis set with polarization functions which gives BSSE of 0.2
0.4 eV.46 A second example is CO and NO adsorption on NiO and MgO where energy was calculated using various methods using triple-ζ basis sets with polarization functions, for which BSSE was in the range 0.1
0.4 eV.47 The last and most relevant example is the adsorption of hydrocarbons, water, and alkoxides on zeolites with the B3LYP and B2PLYP functionals using triple-ζ basis sets and polarization functions.24 For water, the BSSE for B3LYP was 0.04 eV. We have chosen to neglect the effects of the BSSE for the rest of our study, considering that the error is small and that the limited cluster size underestimates the adsorption energy by a similar amount, thus making a fortuitous error cancellation. We also compare all adsorption energies to the adsorption energy of water, so some of the effects for both BSSE and cluster size cancel as well. Comparative Adsorption Study. We now turn to the geochemically more interesting results from the comparative adsorption study of ethanol, water, acetic acid, hydroxyl, and acetate. Because the error from edge effects and relaxation is only 0.10 eV for the neutral molecules, we have performed adsorption studies on either completely frozen 80 atom clusters or 60 atom clusters with 5 atoms relaxed (the MgCO3 unit in the center of the top layer). Not surprisingly, for adsorption in a vacuum, the ions, OH
and CH3COO
, interact much more strongly with the magnesite surface than the protonated neutral species (Table 8). Hydroxyl adsorbs more strongly than acetate, in spite of acetate bridging over two Mg surface atoms. If we now add implicit water by using the COSMO model, the situation changes dramatically. The ions are stabilized in the polar aqueous solvent much more than the neutral species. After introducing implicit solvent, water now adsorbs more strongly than the other species, even though the energy differences for water, ethanol, and acetic acid are within the determined error of 0.10 eV. The acetate is the least stable species, less stable by 0.76 eV, which is even larger than the error from relaxation, 0.62 eV. In the last column of Table 8, we have also presented all adsorption energies relative to adsorbed water. All molecules bind to one Mg except the acetate ion. Thus, acetate replaces two adsorbed water molecules, but all of the other molecules replace one water. All solvated molecules are thus compared to one (or two for the acetate ion) explicit water molecules, as opposed to using only implicit water. We adopted this procedure because it is likely that 10051

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 7b. Effect of Basis Set on Adsorption Energy of Ethanol, Water, Acetic Acid, Hydroxyl, and Acetate on the Frozen 9  9  3 Å Clustera adsorption energy (eV) basis set

ethanol

water

acetic acid

hydroxyl

acetate

default basis 4


0.62


0.70


0.80


3.27


1.88

default basis 5


0.62


0.75


0.78


2.78


1.82

default basis 4 þ diffuse oxygen atoms BSSE default basis 4 þ diffuse oxygen atoms


0.61 0.07


0.72 0.10


0.79 0.10


2.77 0.18


1.81 0.13

a The basis set superposition errors (BSSE) for all species calculated with the default basis 4 þ diffuse oxygen atoms are also included. The B3LYP functional was used in all cases.

Table 8. Adsorption Energy (Eads) on Magnesite {10.4} for Some Small Molecules with Functional Groups Containing Oxygena 6  9  3 Å cluster, relaxed contact atoms molecule

Eads (vac.)

Eads (aq)

Eads (aq)
Eads (vac.)

9  9  3 Å cluster, all atoms frozen Eads (vac.)

Eads (aq)

Eads (aq)
Eads (vac.)

Eads (aq) replacing water at surface

ethanol


0.63


0.22

0.41


0.62


0.20

0.42

0.16

water acetic acid


0.68
0.75


0.36
0.29

0.32 0.46


0.75
0.78


0.36
0.31

0.39 0.47

0.00 0.05

acetate

-

-


1.82

0.04

1.86

0.76

hydroxyl

-

-


2.78


0.29

2.49

0.07

a

Solvation effect was calculated using implicit solvent (water), using the COSMO method and the ORCA program. All energies were calculated as B3LYP/defbas-5 single-point energies using ORCA on geometries obtained from TeraChem using the 6-311G basis and the B3LYP functional. The molecules are adsorbed on the 6  9  3 or the 9  9  3 Å magnesite cluster. All energies are in eV. The lowest adsorption energies for the 9  9  3 Å cluster are shown in bold.

it cancels errors from the BSSE, relaxation, and finite size effects and implicit solvation. Also from the data presented in Table 8, it is clear that the differences in adsorption energies determined assuming vacuum for ethanol, water, and acetic acid for the 60 and 80 atom clusters are 0.01, 0.07, and 0.03, all of which are lower than our error estimate. Our results thus predict that in an aqueous solution neither of the molecules, ethanol or acetic acid, or the acetate or hydroxyl ion are present at the magnesite{10.4} surface. Our error estimate of 0.10 eV makes the prediction uncertain because the relative adsorption energies of acetic acid and hydroxyl ion are within this error, and the relative adsorption energy of ethanol is only slightly higher. It is interesting to compare these predictions with results from a similar system, on calcite {10.4}.14,48 Experiments show that in a mixture of ethanol and water ethanol binds more strongly than water. Molecular dynamics simulations show that for a random 50/50 mixture ethanol replaces adsorbed water. There are at least two explanations for the differences between the results of this study and those published by Sand, Cooke, and colleagues. First, our calculations have been made on the magnesite surface; the smaller hydration radius for Mg compared with Ca could result in much stronger attraction for water than ethanol. Second, the ordered structure of ethanol allows the aliphatic CH3 groups to be in close contact with each other and interact through van der Waals forces. This type of collective behavior is not included in our model and would require an extension to higher coverages. Studying the interaction of nearestneighbor adsorbates is beyond the scope of this paper. In an experimental study on magnesite powders, fatty acids were found to adsorb onto the mineral from an aqueous solution.49 Again this could be an effect of van der Waals interaction between the fatty ends of the acids, which is not included in our model. The

adsorption was furthermore measured using thermogravimetric analysis, and results for the {10.4} surface cannot be singled out. Implicit solvation models have improved in the past few years in parallel with advances in the rest of computational chemistry (e.g., ref 50 and references therein). The COSMO model used in this paper is quite a crude one in comparison with more recent implementations, such as the -RS variant of COSMO,51
53 which has also been successfully applied to adsorption equilibria.54 Many newer models take into account entropic effects as well, providing free energies of adsorption as opposed to differences in total electronic energies that we calculate. Here we try to estimate the errors for solvation energies as follows: We compare adsorption energies using the COSMO model for ethanol, water, and acetic acid adsorbed on both the 60 atom and the 80 atom magnesite cluster. From the data of Table 8, we can calculate the solvation energy contribution to the adsorption energy and compare for the various species. The absolute solvation energy differences for ethanol, water, and acetic acid are 0.01, 0.07, and 0.01 eV, similar in size to errors arising from using the cluster approximation. This means that while these errors are small they are still one of the limiting factors for improving the accuracy of energies for molecules adsorbing from solution onto mineral surfaces. It is also quite likely that some finite size errors are smaller than we have calculated them to be, when implicit solvent is present, because of the screening of long-range electrostatic effects. Quantum Chemical Method. We discussed earlier that using a finite size cluster underestimates the adsorption energies by as much as 0.18 eV, and freezing most or all of the cluster atoms yields an error of 0.09 eV, so we concluded that the computational effort might be better spent elsewhere. We therefore explored the increasingly sophisticated quantum chemical methods of the BLYP family of functionals, BLYP, B3LYP, and 10052

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 9a. Deviations from the B2PLYP-D Method for Adsorption Energies of Ethanol, Acetic Acid, Acetate, and Hydroxyl from the Gas Phase onto Magnesite {10.4} of the Unrelaxed 9  9  3 Å, 80 Atom Clustera PBE

BLYP B3LYP B2PLYP PBE-D BLYP-D B3LYP-D

MUE

0.20

0.38

0.31

0.16


0.02

0.01

0.00

MAE

0.28

0.40

0.31

0.16

0.12

0.12

0.03

max. dev.

0.34

0.57

0.43

0.21

0.09

0.09

0.04

min. dev.
0.18
0.07

0.12

0.10


0.32


0.29


0.07

a

Mean unsigned errors (MUE), mean absolute errors (MAE), as well as maximum and minimum deviations are presented.

Figure 3. Comparison of the PBE, BLYP, B3LYP, and B2PLYP functionals with and without correction for dispersion (-D) for adsorption onto the magnesite {10.4} surface using the 80 atom, 9  9  3 Å cluster. Calculations are for: (a) adsorption from the gas phase, (b) solvation energies for implicit water, (c) adsorption from implicit water, and (d) adsorption from implicit water relative to the adsorption energy of (explicit) water. Error bars of 0.10 eV are added to the most accurate functional, B2PLYP-D, to visualize the difference in quantum chemical method compared with the error from relaxation and edge effects.

B2PLYP, as well as the effect of empirical dispersion. In addition to these functionals, we also include the PBE functional, with and without dispersion, because it is probably the most used density functional for adsorption studies. The results are presented in Figure 3a
d. All calculations were performed as single-point energies using the ORCA defbas-4 basis set on the TeraChem 6-311G geometries. The first point, by comparison with Table 8, is that the lack of diffuse functions in defbas-4 significantly changes the adsorption energy of OH
. If these results are to be used for anything other than comparison between functionals, a recalculation with added diffuse functions to oxygen atoms or an even more accurate basis set is required, so convergence can be reached with respect to the basis set for all species. For B3LYP, we checked the adsorption energy of OH
for the defbas5 and defbas6 bases as well; convergence was reached to within 0.03 eV at the defbas4 level with diffuse functions on oxygen atoms. Table 7b presents data for convergence with respect to diffuse functions. Figure 3 and Table 9a show that adsorption energies from the gas phase are similar for all dispersion corrected functionals. B2PLYP-D results on zeolites are very close to highly accurate coupled cluster results for adsorbed molecules.24 Tongying et al.24 further showed that BLYP-D and B3LYP-D give adsorption energies very close to B2PLYP-D results, and we can confirm this for adsorption on magnesite {10.4}. Larger discrepancies were observed for the uncorrected functionals. The solvation energy, which is the difference between the adsorption energy in solution and in vacuum, is shown for the series in Figure 3b. It is very similar for the three functionals; there is no effect from dispersion because this empirical add-on does not change the electron density. The only deviating point is BLYP for the hydroxyl ion. In Figure 3c, the adsorption energies in implicit water are shown. Here we find a result very similar to adsorption from the gas phase because the solvation energies are so similar. In Figure 3d, all adsorption energies are compared with adsorbed water. Most of the differences between functionals cancel. We see in Table 9b that only small differences between the functionals remain and that the choice of functional is not as crucial as for adsorption from the gas phase. Our results show that PBE-D, BLYP-D, and B3LYP-D provide very reliable adsorption energies regardless of whether the adsorbate originates as a gas or from implicit solvent. The inclusion of empirical dispersion greatly improves the calculated adsorption energies of DFT, particularly the pure density functionals. We suggest that it should be included in future adsorption studies using DFT. BLYP-D was a little more accurate than B3LYP-D for adsorption from a solvent, but this could result from an error cancellation effect. The BLYP-D adsorption energy of OH
10053

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C

ARTICLE

Table 9b. Deviations from the B2PLYP-D Method for Adsorption Energies of Ethanol, Acetic Acid, Acetate, and Hydroxyl Compared with the Adsorption Energy of Water from Implicit Water onto the Magnesite {10.4} Surface of the Unrelaxed 9  9  3 Å, 80 Atom Clustera PBE

0.02
0.06

0.00


0.06


0.01


0.05

0.09 0.16

0.07 0.03

0.04 0.04

0.01 0.03

0.01 0.01

0.05 0.00

min. dev.
0.30
0.12
0.17


0.07


0.26


0.04


0.11

MUE MAE max. dev.


0.06

BLYP B3LYP B2PLYP PBE-D BLYP-D B3LYP-D

0.10 0.09

a

Mean unsigned errors (MUE), mean absolute errors (MAE), as well as maximum and minimum deviations are presented.

in vacuum is 0.29 eV lower than B2PLYP-D, and at the same time, the solvation energy for OH
is 0.35 eV higher than for B2PLYP-D. This less consistent behavior for BLYP-D leads us to recommend the use of B3LYP-D for adsorption energies of organic molecules on mineral surfaces. It is almost as accurate as B2PLYP-D but computationally much cheaper for systems of about 90 atoms, where the poorer scaling of the MP2 part of B2PLYP starts to dominate computational time. On the other hand, if reaction barriers and transition states are of interest, we recommend the use of B2PLYP-D because of its much higher accuracy than B3LYP-D.26 If more sophisticated approximations are used, that include entropic and nearest-neighbor effects on adsorption, accurate phase diagrams for surface coverage as a function of temperature and pH can be drawn similarly to what has been done for gas/solid interfaces, such as refs 55 and 56. This would allow quite accurate predictions of surface composition under given external conditions that influence chemical properties, such as reactivity, and physical properties, such as hydrophobicity.

’ CONCLUSION From our hybrid density functional theory calculations (B3LYP) of the ethanol/magnesite {10.4} system, we can quantify errors resulting from using a finite size mineral cluster, frozen cluster, and varying edge termination. As a result of the finite size of a cluster, there is a systematic underestimation of the binding energy of as much as 0.17 eV for gas phase adsorption and as much as 0.10 eV for adsorption from an aqueous solution. This difference comes from the fact that dispersion effects, resulting from size effects, partly cancel for systems in solution. In addition, there is a random error from relaxation and edge effects that is estimated to be 0.09 eV. This error is similar for molecules binding in a similar fashion to the surface. For the acetate ion, which binds to the surface through two Mg atoms, the relaxational effect was much higher, 0.62 eV. The basis set superposition error (BSSE) for ethanol, water, and acetic acid on magnesite{10.4} is estimated to be 0.05
0.10 eV. The BSSE for the ions, acetate and hydroxyl, is slightly higher, about 0.15 eV. This means that if BSSE is neglected the errors from finite size effects and BSSE partly offset each other, a fortuitous error cancellation. In a comparative study of adsorption from vacuum and aqueous solution, using an implicit solvent model for ethanol, water, acetic acid, hydroxyl, and acetate, we predicted that water adsorbs more strongly than the other species, and none of the other molecules would produce any significant coverage on magnesite {10.4}. We compared several levels of theory: the pure density functionals PBE and BLYP, the hybrid functional B3LYP,

and the double hybrid functional B2PLYP, both with and without dispersion for all functionals. We found for adsorption from the gas phase that the addition of dispersion is crucial to obtain quantitative convergence with respect to the quantum chemical method. The effect is more or less the same for absolute adsorption energies from water modeled using an implicit solvent because the implicit solvation energies are almost identical for all quantum chemical methods studied here. On the other hand, when investigating equilibrium behavior for these kinds of systems, it is only the relative adsorption energies between the competing species that matter. Here the differences between the functionals are much smaller, with mean absolute errors on the order of 0.1 eV or less. We conclude that B3LYP-D is a good method for calculating adsorption energies, both from vacuum and solution. It is more reliable than the pure density functionals with fewer large errors and at the same time more than an order of magnitude faster than the B2PLYP(-D) methods for systems of this size.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]; [email protected].

’ ACKNOWLEDGMENT We are grateful to Andreas Klamt, COSMOlogic GmbH & Co. KG, for valuable discussions on cavity radii in COSMO calculations. Funding was provided by the Danish Center for Scientific Computing (DCSC) and the Nano-Chalk Venture, supported by the Danish National Advanced Technology Foundation (HTF) and Maersk Oil and Gas A/S. ’ REFERENCES (1) Li, B.; Michaelides, A.; Scheffler, M. Surf. Sci. 2008, 602 (23), L135–L138. (2) Tosoni, S.; Sauer, J. Phys. Chem. Chem. Phys. 2010, 12 (42), 14330–14340. (3) Giannozzi, P.; J. Phys.: Condens. Matter 2009, 21 (39), 395502. (4) Dovesi, R. S., V. R.; Roetti, R.; Orlando, R.; Zicovich-Wilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL06; University of Torino: Torino, 2006. (5) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77 (18), 3865–3868. (6) Mortensen, J. J.; Kaasbjerg, K.; Frederiksen, S. L.; Norskov, J. K.; Sethna, J. P.; Jacobsen, K. W. Phys. Rev. Lett. 2005, 95 (21), 4. (7) Becke, A. D. Phys. Rev. A 1988, 38 (6), 3098. (8) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. Rev. B 1988, 37 (2), 785–789. (9) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliott, J. A.; Lardge, J. S.; Duffy, D. M. J. Phys. Chem. C 2007, 111 (32), 11943–11951. (10) Hermann, K.; Bagus, P. S.; Nelin, C. J. Phys. Rev. B 1987, 35 (18), 9467. (11) Panas, I.; Schule, J.; Siegbahn, P.; Wahlgren, U. Chem. Phys. Lett. 1988, 149 (3), 265–272. (12) Dickinson, S. R.; McGrath, K. M. J. Mater. Chem. 2003, 13 (4), 928–933. (13) Chen, S. F.; Yu, S. H.; Jiang, J.; Li, F. Q.; Liu, Y. K. Chem. Mater. 2006, 18 (1), 115–122. (14) Sand, K. K.; Yang, M.; Makovicky, E.; Cooke, D. J.; Hassenkam, T.; Bechgaard, K.; Stipp, S. L. S. Langmuir 2010, 26 (19), 15239–15247. 10054

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055

The Journal of Physical Chemistry C (15) Wright, K.; Cygan, R. T.; Slater, B. Phys. Chem. Chem. Phys. 2001, 3 (5), 839–844. (16) Villegas-Jimenez, A.; Mucci, A.; Whitehead, M. A. Langmuir 2009, 25 (12), 6813–6824. (17) Freeman, C. L.; Asteriadis, I.; Yang, M. J.; Harding, J. H. J. Phys. Chem. C 2009, 113 (9), 3666–3673. (18) Villegas-Jimenez, A.; Mucci, A.; Pokrovsky, O. S.; Schott, J. Geochim. Cosmochim. Acta 2009, 73 (15), 4326–4345. (19) Lardge, J. S.; Duffy, D. M.; Gillan, M. J.; Watkins, M. J. Phys. Chem. C 2010, 114 (6), 2664–2668. (20) Grimme, S. J. Comput. Chem. 2006, 27 (15), 1787–1799. (21) Klamt, A.; Schuurmann, G. J. Chem. Soc., Perkin Trans. 2 1993, No. 5, 799–805. (22) Becke, A. D. J. Chem. Phys. 1993, 98 (7), 5648–5652. (23) Grimme, S. J. Chem. Phys. 2006, 124 (3), 16. (24) Tongying, P.; Tantirungrotechai, Y. J. Mol. Struct.-Theochem 2010, 945 (1
3), 85–88. (25) Schwabe, T.; Grimme, S. Phys. Chem. Chem. Phys. 2006, 8 (38), 4398–4401. (26) Zheng, J. J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2009, 5 (4), 808–821. (27) Biczysko, M.; Panek, P.; Scalmani, G.; Bloino, J.; Barone, V. J. Chem. Theory Comput. 2010, 6 (7), 2115–2125. (28) Stroppa, A.; Kresse, G. New J. Phys. 2008, 10, 063020. (29) Paier, J.; Marsman, M.; Kresse, G. J. Chem. Phys. 2007, 127 (2), 024103. (30) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114 (12), 5149–5155. (31) Grimme, S. J. Comput. Chem. 2004, 25 (12), 1463–1473.  erny , J.; Hobza, P.; Salahub, D. R. J. Comput. (32) Jurecka, P.; C Chem. 2007, 28 (2), 555–569. (33) Ufimtsev, I. S.; Martinez, T. J. J. Chem. Theory Comput. 2009, 5 (10), 2619–2628. (34) Neese, F. ORCA, 2.6; University of Bonn: Germany, 2008. (35) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7 (18), 3297–3305. (36) Schafer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97 (4), 2571–2577. (37) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. J. Chem. Phys. 2003, 119 (23), 12129–12137. (38) Prencipe, M.; Pascale, F.; Zicovich-Wilson, C. M.; Saunders, V. R.; Orlando, R.; Dovesi, R. Phys. Chem. Miner. 2004, 31 (8), 559– 564. (39) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Chem. Phys. 2009, 356 (1
3), 98–109. (40) Mantina, M.; Chamberlin, A. C.; Valero, R.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2009, 113 (19), 5806–5812. (41) VMD/NAMD/BioCoRE/JMV/ is developed with NIH support by the Theoretical and Computational Biophysics group at the Beckman Institute, University of Illinois at Urbana
Champaign. (42) Lardge, J. S.; Duffy, D. M.; Gillan, M. J. J. Phys. Chem. C 2009, 113 (17), 7207–7212. (43) Kohn, W. Phys. Rev. Lett. 1996, 76 (17), 3168–3171. (44) Prodan, E.; Kohn, W. Proc. Natl. Acad. Sci. U.S.A. 2005, 102 (33), 11635–11638. (45) Boys, S. F.; Bernardi, F. Mol. Phys. 2002, 100 (1), 65–73. (46) Rimola, A.; Civalleri, B.; Ugliengo, P. Phys. Chem. Chem. Phys. 2010, 12 (24), 6357–6366. (47) Valero, R.; Gomes, J. R. B.; Truhlar, D. G.; Illas, F. J. Chem. Phys. 2010, 132 (10), 104701. (48) Cooke, D. J.; Gray, R. J.; Sand, K. K.; Stipp, S. L. S.; Elliott, J. A. Langmuir 2010, 26 (18), 14520–14529. (49) Thomas, M. M.; Clouse, J. A.; Longo, J. M. Chem. Geol. 1993, 109 (1
4), 201–213. (50) Cramer, C. J.; Truhlar, D. G. Acc. Chem. Res. 2009, 42 (4), 493–497. (51) Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. C. W. J. Phys. Chem. A 1998, 102 (26), 5074–5085.

ARTICLE

(52) Klamt, A.; Diedenhofen, M. J. Comput.-Aided Mol. Des. 2010, 24 (4), 357–360. (53) Klamt, A.; Eckert, F.; Arlt, W. COSMO-RS: An Alternative to Simulation for Calculating Thermodynamic Properties of Liquid Mixtures. In Annual Review of Chemical and Biomolecular Engineering; 2010; Vol. 1, pp 101
122. (54) Mehler, C.; Klamt, A.; Peukert, W. AIChE J. 2002, 48 (5), 1093–1099. (55) Andersson, M. P.; Abild-Pedersen, E.; Remediakis, I. N.; Bligaard, T.; Jones, G.; Engbwk, J.; Lytken, O.; Horch, S.; Nielsen, J. H.; Sehested, J.; Rostrup-Nielsen, J. R.; Norskov, J. K.; Chorkendorff, I. J. Catal. 2008, 255 (1), 6–19. (56) Reuter, K.; Scheffler, M. Phys. Rev. Lett. 2003, 90 (4), 046103.

10055

dx.doi.org/10.1021/jp1121189 |J. Phys. Chem. C 2011, 115, 10044–10055