Theoretical Surface Science Beyond Gradient-Corrected Density

Subscriber access provided by WEBSTER UNIV

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Theoretical Surface Science Beyond Gradient-Corrected Density Functional Theory: Water at #-AlO(0001) as a Case Study 2

3

Sophia Heiden, Denis Usvyat, and Peter Saalfrank J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b00407 • Publication Date (Web): 25 Feb 2019 Downloaded from http://pubs.acs.org on February 27, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Theoretical Surface Science Beyond Gradient-Corrected Density Functional Theory: Water at α-Al2O3(0001) as a Case Study Sophia Heiden1 , Denis Usvyat2 , and Peter Saalfrank1∗ 1

Institut für Chemie, Universität Potsdam,

Karl-Liebknecht-Straße 24-25, 14476 Potsdam-Golm, Germany 2

Institut für Chemie, Humboldt-Universität zu Berlin Brook-Taylor-Straße 2, 12489 Berlin, Germany ∗

E-mail: [email protected] Phone: +49 331 977-5232 February 23, 2019

Abstract The quantum chemical description of the adsorption, vibrations and reactions of molecules at periodic, solid surfaces is frequently based on a methodological “standard model”: Density functional theory (DFT) in the generalized gradient approximation (GGA), using plane wave bases and three-dimensional supercells. While the computationally efficient GGA functionals can be very successful, cases are known where they don’t perform so well. Most importantly, activation energies for chemical reactions are typically underestimated – with the consequence of computed reaction rates being too large. In this work we consider a well-studied model system: water or water fragments adsorbed on an Al-terminated, α-Al2 O3 (0001) surface, as a test bed for studying the performance of electronic structure methods, both from DFT and wavefunction theory (WFT). On the DFT side, we employ two GGA exchange correlation functionals: PW91 and PBE with and without dispersion corrections, whose results are then compared to those of hybrid functionals B3LYP and HSE06. Further, we follow a periodic wavefunction approach in the form of local second-order Møller-Plesset perturbation theory, LMP2, on a Hartree-Fock (HF) reference. En route, we address issues arising from the choice of the basis set. The key

1 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

findings of our study are: (i) DFT-GGA adsorption energies are in reasonable agreement both with hybrid-DFT and LMP2 values. In particular, the deviations between the relative energies, corresponding to different adsorption structures, are in the range of the error due to missing dispersion corrections or the basis set error. (ii) Harmonic DFT-GGA vibrational frequencies for oxygen-hydrogen stretch modes are by several tens of wavenumbers red-shifted compared to corresponding hybrid-DFT values. The latter are in much better agreement with recent experimental data. (iii) The activation energy for a hydrogen diffusion reaction is grossly underestimated by GGA compared to hybrid-DFT or LMP2, which in turn are quite comparable.

I.

Introduction The adsorption of small molecules on solid surfaces is decisive for many impor-

tant fields, such as heterogeneous catalysis, corrosion or electrochemistry. Besides surface science experiments, the quantum chemical modelling provides important cornerstones for a fundamental understanding of relevant properties and processes. For the latter, the question arises which level of theory provides the best balance between effort (the computational cost) and accuracy. The issue of efficiency is particularly important when treating systems which are complex, e.g., adsorption on reconstructed or defective surfaces. The quest for accuracy depends on the property of interest, usually in the context of experimental determined structures, spectra, and reactivity. In the last two decades or so, an efficient and quite accurate, workhorse or “standard model” emerged for the electronic structure of adsorbate / surface systems: Density functional theory (DFT) in the generalized gradient approximation (GGA),1, 2 using plane wave (PW) bases and three-dimensional supercells.3 On the opposite model / method side, wave function theory (WFT) of various flavours, ranging from Hartree-Fock theory to a multitude of correlated methods, is available, both for periodic and non-periodic (e.g., cluster) models,4–10 typically realized with atomic-orbital (AO) bases. A connection between the two worlds is hybrid-DFT,11, 12 in which a portion of HF-like exchange is used in the exchangecorrelation functional, Exc . Compared to pure (GGA-)DFT, both hybrid-DFT and correlated WFT are “expensive”, notably for periodic systems. The question is then


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


under which circumstances / for which properties the higher effort is worth it. In this paper, we address this question for a concrete example, the adsorption of water on an Al-terminated α-Al2 O3 (0001) (alumina) surface. Of all cuts of alumina, this reconstructed surface is the most stable one under UHV (ultra high vacuum) conditions,13 and widely studied experimentally.14–16 Many experiments also refer to water adsorption on the surface.17, 18 Here we are particularly interested in the low-coverage regime, with less than a monolayer water coverage. This situation was studied for example in Refs.,19, 20 by Vibrational Sum Frequency (VSF) spectroscopy and other methods. Also for theory, water on Al-terminated α-Al2 O3 (0001) was studied by several groups including ours.19, 21–25 These studies were aiming to interpret the experimental findings and to contribute to the understanding of water adsorption on metal oxides in general. Most previous theoretical studies applied GGA based DFT, the “standard model” mentioned above. While DFT-GGA can be quite accurate, it is known to have also its drawbacks. In particular, barrier heights for reactions are often underestimated26 and computed reaction rates are expected to be too high. Hybrid functionals are believed to generally outperform GGAs.27 The actual physical reason why addition of the HF-like exchange should lead to an improved description and, if so, for which properties or systems, is still debatable. Clearly, a fraction of the HF exchange reduces the selfinteraction error and the associated overdelocalization in DFT. From another angle, one can justify inclusion of the HF exchange by means of the adiabatic connection formula for the exchange-correlation energy.11 In contrast to DFT, WFT presents a clear hierarchy of models, that provide a physically justified progressively improved description. MP2 stands in this hierarchy not too high, but its results are definitely substantially better than those of HF (the lowest method in the hierarchy) and, for intermediately polarizable systems, commonly better than those of DFT. The advantage of MP2 is that it treats both long- and short-range correlation. In DFT, among other methods, the latter is frequently treated in some ad hoc (but nevertheless quite successful) way, as a semi-empirical dispersion (“+D”) correction.28, 29 MP2 is known to be especially good for hydrogen-bonded systems and, in particular, for water interactions30–33 including barriers for hydrogen atom transfer.34 3 Environment ACS Paragon Plus


Another issue of interest is the use of atomic orbital vs. plane wave bases, which also relates to computational effort as mentioned. In particular hybrid-DFT/PW calculations are known to be much more expensive compared to DFT-GGA/PW, while, as demonstrated below, using AO bases instead makes the effort of hybridDFT comparable to that of DFT-GGA. Further, AO bases are quite common in the WFT world, both for molecules and periodic solids.35–39 As a clear disadvantage of AO over PW bases in periodic calculations, however, the convergence of results with respect to the basis size is sometimes slow and harder to control. Further, possible basis set linear dependencies and basis set superposition errors (BSSE) have to be taken care of. Further, correlated WFT requires typically larger basis sets for convergence than, e.g., HF or DFT. Specifically for H2 O/α-Al2 O3 (0001), we calculate adsorption energies of molecular and dissociated water, hydroxyl vibrational frequencies and activation energies and rates for a selected reaction: a hydrogen diffusion between two adsorption sites. We, firstly, test the performance of two pure DFT-GGA exchange correlation functionals, with and without dispersion corrections. To go beyond DFT-GGA, we apply periodic hybrid-DFT with a portion of “exact” exchange, as well as WFT in the form of local Møller-Plesset perturbation theory of second order, LMP2.9, 40 For the latter, the periodic Hartree-Fock (HF) solution serves as a reference. En route, we also address the issues of the atomic-orbital bases, such as a systematic improvement of the AO basis and the basis set superposition error. The paper is organized as follows. In section II., surface models and theoretical methods used are described. In Sec. III., we present and discuss the results, focusing on the possible improvements in the description coming from the hybrid-DFT or WFT approaches. A final Sec. IV. summarizes and concludes our findings.


Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


II. A

Surface models and computational methods Models

In this study we apply different approaches and models to describe the adsorption, vibrations and reactivity of water and water fragments on an Al-terminated, α-Al2 O3 (0001) surface. In particular, we consider a (2×2) cell with four Coordinatively Unsaturated aluminium Sites (CUS) as shown in Fig.1(a), which readily react with water. We refer to a low-coverage situation with one water molecule per (2×2) cell (coverage=1/4), which has been theoretically studied previously, mostly by DFT-GGA.19, 21, 23, 24 In Ref.,19 this model was used to interpret vibrational (more precisely: VSF) spectra. We will perform two major sets of calculations below, namely (i) DFT-GGA calculations adopting a periodic, three-dimensional supercell model using plane-wave bases and the Vienna Ab initio Simulation Package (VASP), version 5.3/4,3, 41–44 and (ii) DFT-GGA, hybrid-DFT, periodic Hartree-Fock and LMP2 calculations employing a 2D-periodic slab model, AO bases and the CRYSTAL1445–48 (for DFT and HF) or CRYSCOR codes35, 40 (for LMP2), respectively. For comparability, the same structural models / unit cells were chosen for both cases (VASP and CRYSTAL/CRYSCOR), namely the one of Refs.19, 25 The model consists of a 9-layer substrate model (Al4 -O12 -Al4 -Al4 -O12 -Al4 -Al4 -O12 -Al4 ) as shown in Fig.1(b) for the naked (2×2) surface, plus a single water molecule to ensure the 1/4 coverage. Depending on the method and target quantity, geometry optimizations were performed or single point calculations were carried out. In case of geometry optimizations, in addition to adsorbate atoms the five uppermost layers of the surface were allowed to move (see Fig.1(b)), while the four lowest layers were kept fixed to their bulk positions optimized in Ref.19 by PBE+D2/PW. The substrate model has a hexagonal lattice with the lattice constant |a| = 9.66 Å. The plane wave calculations employed 3D-periodic model with the vertical lattice constant |c| = 31.4 Å, corresponding to a vacuum gap of 26.4 Å between neighbouring slabs.



(a)

CUS

Page 6 of 33

(b)

CUS 4 2

CUS

CUS

(c)

Figure 1: (a) Top view of the surface model of the (0001) α-alumina surface, (2×2) supercell. CUS Al atoms are indicated, and O sites labelled 2 and 4 with different distances from a reference CUS site (red label). Colour code: Large grey spheres are Al, pale red spheres are substrate O. (b) The sideview, showing the 9-layer slab Al4 -O12 Al4 -Al4 -O12 -Al4 -Al4 -O12 -Al4 . The horizontal line separates atoms whose positions were movable during geometry optimizations, from atoms kept fixed at positions corresponding to those in the bulk material (optimized in Ref.19 by PBE+D2/PW). (c) The three most stable adsorbed water species, at the 1/4 coverage (one H2 O per (2×2) cell), top view: The left panel shows the molecular adsorption, middle the 1-2 dissociated and the right panel the 1-4 dissociated species. Small red spheres: Water O, small white: H. Optimized geometries were obtained in this case with CRYSTAL, using B3LYP+D3 and AO basis set 3 (see below for the specification).

B

Computational methods In all DFT-GGA plane wave calculations within the 3D supercell approach, we

solve Kohn-Sham equations for the periodic system with VASP. We apply the Projector Augmented Wave (PAW) method,44, 49 using GGA functionals in two flavours: the Perdew-Becke-Ernzerhof (PBE)1, 50 and the Perdew-Wang (PW91)2 exchangecorrelation functionals, respectively. In order to account for dispersion, empirical corrections were introduced, using Grimme’s D228 or D3 schemes.51 As a planewave energy cutoff we use Vc = 400 eV. The total energy was considered to be


Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


converged in the Self Consistent Field (SCF) scheme when differences between two cycles were smaller than 10−5 eV. For geometry optimizations, structural convergence was achieved when forces acting on the ions were below 0.01eV/Å. We used the same converged (3 × 3 × 1) Γ-point centred k-point grid as in prior work,24 with a total of five irreducible k-points in the Brillouin zone. In one example, we also did hybrid-DFT single point calculations with VASP using the HSE06 functional,12 with otherwise unchanged computational parameters. DFT-GGA and hybrid-DFT calculations were also performed with CRYSTAL14,45–48 in this case using 2D slab models and AO bases. Specifically, PBE+D3/AO calculations were compared with corresponding plane wave / 3D supercell calculations. As a hybrid-DFT model, the B3LYP functional11, 52 was chosen, together with Grimme’s D3 dispersion.51 Note that B3LYP is a global hybrid with a fraction of 0.2 of “exact” exchange, while HSE0 is a range-separated hybrid with a portion of 0.25 of “exact” exchange at short range. CRYSTAL and 2D slab / AO models were also adopted for periodic Hartree-Fock calculations.4, 5 By definition, this method is correlation-free. Correlation effects were included by periodic Møller-Plesset perturbation theory,53 in the form of Local Møller-Plesset perturbation theory of second order (LMP2)9, 40 as implemented in the CRYSCOR code. In the LMP2 scheme, the occupied HF orbitals are localized prior to their use for the correlation energy. To span the virtual space, which is truncated for each occupied orbital pair, we employed so called orbital-specific virtuals (OSVs).40 MP2 captures short-range electron correlation effects as well as long-range ones. The latter can be associated with the van der Waals dispersion, which in contrast to the D2 or D3 approach is described in MP2 non-empirically. For CRYSTAL DFT and HF calculations, slightly tighter computational settings were used than for DFT-GGA/PW: Energies were considered converged when differences between the total energies of two subsequent SCF cycles were below 2.7 × 10−6 eV, and sampling in ~k-space was done on a (4 × 4 × 1) ~k-point grid. Various AO basis sets were tested for the DFT and HF calculations using CRYS7 Environment ACS Paragon Plus


Page 8 of 33

TAL, and LMP2 calculations using CRYSCOR. In the end, depending on the properties calculated or method employed the following basis sets were chosen in order to optimize the ratio of accuracy to computational cost (a more detailed specification is given in the supporting information): Basis 1 (or AO1) is a moderately big basis set from the POB-TZVP family of basis sets constructed specifically for periodic systems.54 With this basis set we performed the most computationally demanding transition state optimizations at the B3LYP+D3 level. Basis 2 (or AO2) is an improved version of basis AO1. It was mainly used in the computationally intensive calculations of the vibrational frequencies. Basis 3 (or AO3) is a quite a rich basis of triple zeta quality, constructed by combining the largest periodic basis sets available and high angular momentum orbitals from cc-pVTZ basis sets.55 This basis was used for geometry optimization using PBE or B3LYP and for single point HF / LMP2 calculations. Basis 3-dual (or AO3-dual) is an even bigger basis set, constructed from AO3 by adding diffuse d- and f-orbitals for the oxygen and aluminium atoms. It was employed as the bigger basis within the dual basis set framework for single point LMP2 calculations.56

C

Computed Properties Geometry optimizations of a single water molecule in the (2×2) cell gave both

molecular and dissociated species24 (see below). The (B3LYP+D3/AO1) optimized structures of the molecular and two dissociated species are given in Fig.1(c). Their relative stabilities can be estimated from adsorption energies, calculated as Eads = EH2 O/surf − (EH2 O + Esurf )

(1)

from the energy of the adsorbed species (EH2 O/surf ), and the energies of free H2 O (EH2 O , computed in a periodically repeated box coinciding with the supercell) and the naked surface, Esurf . (The AO basis allows for calculating EH2 O also in the non-periodic regime, but that leads to only about 0.01 eV of energy difference to


Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the periodic calculation.) For each minimum structure, with selected methods, normal mode analyses were performed by diagonalizing the force constant matrix at the Γ-point. To compare with experimental VSF spectra which measured O-D frequencies,19 the corresponding hydrogen atoms were replaced by deuterium. In some cases, we also introduced diagonal anharmonic corrections.1 For a selected reaction, i.e., a hydrogen diffusion on the surface between two stable sites (see below), the transition state was optimized using DFT calculations, either with VASP and PWs, or with CRYSTAL and AO basis sets. Since pure density functionals are expected to underestimate barriers (see above), we also performed the transition-state structure optimization using B3LYP/AO1 (with CRYSTAL), as well as single-point calculations of the barriers using HSE06 (with VASP/PW), B3LYP/AO3 (with CRYSTAL) and LMP2/AO3-dual (with CRYSTAL/CRYSCOR). For the reactants and transition states, normal mode analyses were performed and the activation free energies were computed as ∆G‡ (T ) = G‡ (T )−Greactant (T ), where G(X, T ) = E(X)+ZPE+∆Gvib (X, T ) with X = ‡ (transition state) or reactant. Here, E is the electronic energy, ZPE the zero-point energy and ∆Gvib the vibrational contribution to the free energy at a given finite temperature, T , calculated in harmonic approximation. With ∆G‡ obtained in this way, the corresponding reaction rate constant k is calculated using the Eyring equation57 as k(T ) =

kB T exp(−∆G‡ /(kB T )) , h

(2)

with Boltzmann’s constant kB and Planck’s constant h. 1 That is, by varying the O-D distances along O-D bond directions, potential curves V (rO−D ) were obtained and fitted to a polynomial form. Then, one-dimensional Schrödinger equations 2 [− 2m1D d2 /drO−D + V (rO−D )]ψv = Ev ψv were solved, and anharmonic frequencies determined as ωanh = E1 − E0 . This method accounts for “diagonal” anharmonicities, while effects due to the multidimensionality of the potential surface as well as intermode couplings are neglected.



III. A

Results and Discussion Adsorption energies

For a single water molecule adsorbed on a clean (2 × 2) supercell (coverage 1/4), water was shown19, 21 to either adsorb molecularly (denoted “mol” in what follows), or dissociatively. The two most stable structures corresponding to the dissociative adsorption we denote as “1-2” and “1-4” as in Ref.24 The different adsorbed species are depicted in Fig.1(c). The figure shows optimized structures obtained with B3LYP+D3/AO3, however, the corresponding PBE or PW91 structures are very similar.19, 24 The 1-2 species is dissociated in direct neighbourhood, whereas for the 1-4 position the dissociated hydrogen atom is adsorbed to a surface oxygen atom farther away (see Figs.1(a) and (c)). There is agreement in the literature21, 22, 24 that the 1-2 dissociated species is the most stable one. The relative stability between the molecular vs. 1-4 adsorbed structures is yet under debate. It should also be said that at least for higher coverages, of one water per CUS (corresponding to four water per (2×2) cell), dissociation of molecular water on short timescales has recently been questioned on the basis of IR spectra.18 While activation energies for reactions will be considered later, we here already state that this could suggest large activation barriers between molecular and dissociated adsorbates. These were not found by DFT-GGA,21, 24 however, which predicts activation energies in the order of only ∼ 0.2 eV.21, 24 Table 1 presents the adsorption energies for the three species, “mol”, “1-2” and “1-4”, obtained with various methods and basis sets. The GGA/PW data, apart from PBE+D/PW has been taken from our previous work.19, 24 For all the methods, the 1-2 dissociated species is found to be the most stable one, with adsorption energies ranging from -1.59 to -1.81 eV depending on method. The relative stability of the 1-2 structure with respect to the molecularly adsorbed water, ∆1-2/mol = Eads (1-2) − Eads (mol) is always negative. In the PW calculations the spread of ∆1-2/mol for different GGA functionals is only about 0.05 eV around a mean of ∼ −0.35 eV. However, for the 1-4 structure, the sign of the relative stability with respect to the molecular adsorption ∆1-4/mol = Eads (1-4) − Eads (mol) does depend on the method, although in all the calculations ∆1-4/mol is quite close 10 Environment ACS Paragon Plus

Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Adsorption energies for a water molecule on the α-Al2 O3 (0001) (2×2) surface for three adsorption structures: “mol”, “1-2” and “1-4”, as well as the relative stabilities ∆ between them. The PW results were obtained with the VASP code, the AO ones with the CRYSTAL (DFT or HF) and CRYSCOR (LMP2) codes. For the specification of the AO basis sets see Sec.II.B. All geometries were optimized on the respective level, except for HF and LMP2 where single-point calculations at the B3LYP+D3/AO3 geometries were used. “CPC” denotes counterpoise-corrected values. All energies are in eV. method /basis PW91/PW1 PW91+D2/PW1 PBE+D2/PW PBE+D3/PW PBE+D3/AO1 PBE+D3/AO2 PBE+D3/AO3 PBE+D3/AO1 CPC PBE+D3/AO3 CPC B3LYP+D3/AO3 HF/AO32 LMP2/AO32 LMP2/AO3-dual2

mol -1.25 -1.40 -1.31 -1.29 -1.48 -1.47 -1.41 -1.33 -1.32 -1.43 -1.14 -1.34 -1.31

1-2 diss -1.59 -1.81 -1.69 -1.63 -1.73 -1.69 -1.68

1-4 diss -1.25 -1.45 -1.21 -1.30 -1.40 -1.35 -1.32

∆1-2/1-4 -0.34 -0.36 -0.48 -0.33 -0.33 -0.34 -0.36

∆1-2/mol -0.34 -0.41 -0.38 -0.34 -0.25 -0.22 -0.27

∆1-4/mol ± 0.00 -0.05 +0.10 -0.01 +0.08 +0.12 +0.09

-1.81 -1.67 -1.69 -1.61

-1.40 -1.19 -1.26 -1.18

-0.41 -0.48 -0.43 -0.43

-0.38 -0.53 -0.35 -0.30

+0.03 -0.05 +0.08 +0.13

1

From Ref.;24

2

single point calculations at B3LYP+D3/AO3 geometries.

to zero. Comparing PW91 with PW91+D2, dispersion is found to stabilize all adsorbed forms (by up to about 0.2 eV). Comparison of PBE+D2 and PBE+D3 shows that there are also small differences (< 0.1 eV) when different +D corrections are used. Finally we note that the relative energy between the 1-2 and 1-4 structures ∆1-2/1-4 = Eads (1-2) − Eads (1-4) is around -0.4 eV. The LMP2 and B3LYP ∆1-2/1-4 mutually agree here and lie between the PBE-D2 and PBE-D3 values, which deviate by 0.15 eV. The AO-based GGA results are generally comparable to the PW ones, both for the adsorption energies and the relative stabilities. As mentioned above, AO basis set calculations are known to suffer from the basis set incompleteness errors, which are not always easy to eliminate or even estimate. In Table 1, we investigate these issues on the example of the PBE+D3 calculations comparing AO results to the PW ones, for the very same functional. With expansion of the basis (AO1 to AO3) 11 Environment ACS Paragon Plus


the AO-based adsorption energies and the relative stabilities appear to approach the PW values. However, even for the AO3 basis set, the adsorption energy for the molecular adsorbate still deviates from the PW result by about 0.1 eV. The PBE+D3/AO adsorption energies approach the PW ones from below, suggesting that the deviation is to a large extent due to the BSSE in the AO calculations. Indeed, counterpoise corrected (CPC) “mol” interaction energies for both AO1 and AO3 are very close to the PW results. Unfortunately, for the dissociated adsorption cases “1-2” and “1-4”, the BSSE could not be corrected successfully using the counterpoise scheme, as it would involve isolated H and OH radicals rather than stable species studied here. Furthermore, since the discrepancy between AO3 and PW adsorption energy is noticeable larger for “mol” than for the “1-2” case, the BSSE is likely to be not uniform for the three structures and thus does not cancel in the relative stabilities. This suggests that the tendency of the AO-based calculations to stabilize the “mol” adsorption with respect to the “1-2” one might be non-physical and driven by the BSSE. To summarize the DFT-GGA results of Table 1 so far, the “1-2” structure is the global minimum for the water adsorbed on α-Al2 O3 (0001). The “1-4” and “mol” structures are energetically very close to each other with the adsorption energy being around 0.3-0.4 eV higher than for “1-2”. Next, from Table 1 we see that neither hybrid DFT nor MP2 provide essentially new results for the adsorption energies compared to DFT-GGA. Indeed, as follows from the table the LMP2 energies are very close for example to those of PBE+D2, while the B3LYP+D3 ones to those of PW91+D2. HF seems to underbind, in particular the less stable species. Without an accurate experimental reference, however, it is difficult to judge which of these methods are the most accurate ones. From this we conclude that for the adsorption energies of water on α-Al2 O3 (0001), AO-based hybrid DFT or MP2 calculations with a triple-zeta basis hardly provide much improvement over the standard PWbased GGA. At the same time we note that the WFT approach can potentially deliver a much better description than DFT, if one further corrects both the basis set and method errors of LMP2 by evaluating a ∆CCSD(T) (coupled cluster singles doubles with perturbative triples) correction and extrapolating the interaction energies to the basis set limit.39, 58–60 This can be done both by finite, embedded 12 Environment ACS Paragon Plus

Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cluster models or periodically. For the system size at hand, however, a CCSD(T) calculation – as an alternative to an experimental reference – is currently a challenge. We finally note that adsorption geometries differ only very slightly between AO-based DFT-GGA and B3LYP, and also in comparison with the PW-based DFT-GGA calculations. This is demonstrated in Table 2 for bond lengths and angles of the OH bond with respect to the surface normal, for PBE+D3/AO3, B3LYP+D3/AO3 and PBE+D3/PW, respectively. Table 2: Optimized geometric parameters of various adsorbed species: OH bond lengths (BL) and OH angles (θ) with respect to the surface normal. Bond lengths are given in Å and angles in °. For the dissociated species, the first number refers to the adsorbed OH group and the second to the surface OH group. The first block gives results for PBE+D3/PW, the second for PBE+D3/AO3 and the last one for B3LYP+D3/AO3. method / basis PBE+D3/PW PBE+D3/AO3 B3LYP+D3/AO3

B

parameter BL θ BL θ BL θ

mol 0.98, 0.99 87.9, 99.1 0.98, 0.98 87.4, 97.9 0.97, 0.97 85.0, 95.8

1-2 diss 0.97, 0.98 51.7, 34.9 0.96, 0.98 52.0, 34.2 0.95, 0.97 49.8, 33.0

1-4 diss 0.97, 0.98 51.2, 27.2 0.96, 0.98 52.4, 26.5 0.96, 0.97 49.9, 26.7

Vibrational frequencies As a second set of properties, we looked at the vibrational frequencies of adsorbed

species. To relate the calculated frequencies to the experimental VSF spectra of OD vibrations, corresponding to the D2 O adsorbed in low coverage on Al-terminated α-Al2 O3 (0001),19 the hydrogen atoms were replaced by deuterium. In Ref.,19 an assignment of experimental peaks to the computed O-D vibrations of the dissociated species (molecularly adsorbed species were not observed in that experiment) was made based on frequency differences between the individual peaks. For each dissociated water, two O-D vibrations are expected, one from the water O-D fragment (denoted as “ads”), another one from the O-D formed when D (or D+ ) of the dissociated water attaches to a nearby surface O atom (denoted as “surf”).



Vibrational frequencies computed for various species with various methods are shown in Table 3. The PBE+D2/PW frequency calculations were performed in Ref.19 We include these results in Table 3 for comparison, along with the experimental values. While the PBE+D2 and the experimental vibrational energy differences match quite well, the absolute frequencies noticeably deviate, with computed frequencies being considerably red-shifted with respect to experiment. In this work to address the differences due to the basis (plane wave vs. AO) and, more importantly, due to inclusion of the “exact” exchange, we computed the same quantities at the PBE+D3/AO2 and B3LYP+D3/AO2 level using CRYSTAL. The corresponding data are also given in the table. Besides the pure harmonic frequencies ν˜, we also provide the scaled ones in the brackets, which are obtained by applying scaling factors of 1.031 for PBE and 1.004 for B3LYP, respectively. These factors are rounded averages over scaling factors recommended in Ref.61 for molecules, to compare with measured frequencies. We also calculated anharmonically corrected wavenumbers νãnh (without scaling), as described in Sec.II.C.2 In the CRYSTAL code this procedure is automatized. We observe that for PWs and AOs both harmonic and anharmonic PBE frequencies (+D2/PW, +D3/AO2) are in reasonably good agreement. This suggests that the vibrational frequencies are appreciably insensitive either to the basis set issues in the AO-based calculations, or to the choice of dispersion corrections. We find that PBE+D2 unscaled harmonic frequencies are redshifted with respect to experiment by more than 100 cm−1 . However, frequency differences ∆1 to ∆3 are in good agreement with experiment as mentioned above and in Ref.19 In contrast to the adsorption energies, for the vibrational frequencies B3LYP shows a profound improvement compared to PBE, at least when the harmonic approximation is invoked. Unscaled harmonic B3LYP+D3 frequencies deviate from experimental ones much less than PBE, by at most 50 cm−1 and by around 30 cm−1 on average. After scaling, however, both PBE and B3LYP values become close to 2

The PBE+D2/PW anharmonic frequencies shown in Table 3 were evaluated in Ref.19 by fitting the potential curve along the normal mode (rather than O-D bond, as here), to a Morse potential (rather than to a polynomial as here).


Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: The frequencies of the O-D stretch vibrations of H2 O/α-Al2 O3 (0001) (2×2). The harmonic ν˜ and anharmonic νãnh frequencies were calculated with PBE+D3/AO2, B3LYP+D3/AO2 (both using CRYSTAL) and PBE+D2/PW (using VASP). The scaled harmonic frequencies with the scaling factors of 1.031 and 1.004 for PBE and B3LYP, respectively, are given in the parentheses. The available experimental values assigned to the O-D vibrations are also included. We also provide the frequency differences between 1-2 ODads and 1-2 ODsurf , 1-4 ODsurf , and (1-4 ODads : ∆1 = ν˜(1-2 ODads )-˜ ν (1-2 ODsurf ), ∆2 = ν˜(1-2 ODads )-˜ ν (1-4 ODsurf ), −1 ∆3 = ν˜(1-2 ODads )-˜ ν (1-4 ODads ). All values in cm .

stretch mode mol: OD1 mol: OD2 1-2: ODsurf 1-2: ODads 1-4: ODsurf 1-4: ODads ∆1 ∆2 ∆3 1

PBE+D3/AO2 ν˜ νãnh 2657 (2739) 2562 2539 (2618) 2490 2596 (2676) 2521 2808 (2895) 2727 2621 (2702) 2531 2795 (2882) −1 212 (219) 206 187 (193) 196 13 (13)

B3LYP+D3/AO2 ν˜ νãnh 2747 (2758) 2650 2627 (2638) 2584 2697 (2708) 2623 2883 (2895) 2805 2715 (2726) 2631 2873 (2884) 2788 186 (187) 182 168 (169) 174 10 (11) 17

PBE+D2/PW19 ν˜ νãnh 2664 (2747) –1 2550 (2629) 2486 2629 (2710) 2523 2810 (2897) 2713 2647 (2729) 2541 2795 (2882) 2696 181 (187) 190 163 (168) 172 15 (15) 17

The SCF did not converge for some of the points on O-D potential energy curve, or the fit was

unreliable, so the corresponding anharmonic frequency could not be computed.

each other and agree well with the experimental values. Interestingly, correcting the unscaled harmonic frequencies anharmonically, noticeably worsens the agreement with experiment for B3LYP and even more so for PBE. Indeed, the computed anharmonic corrections are relatively large (up to 100 cm−1 ) and shift the frequencies (further) to the red. The reason of this worsening is not entirely clear. It could be related to the idealized structure model (ideal surface, exact 1/4 coverage). In the real experimental structure the potential might be less anharmonic due the presence of the environment molecules. Another possible reason could be an incorrect interpretation or assignment of experimental peaks. In passing we also provide the frequencies of the O-D stretch vibrations for D2 O adsorbed on the (11¯20) surface of alumina (one molecule per (2 × 2) cell) in Table 4. The experimental VSF results as well as the calculated one are taken from a recent publication Ref.25 For this surface, there is a larger variety of possible (dis-


Exp.19 ν˜ – – 2729 2910 2764 2900 191 146 10


Page 16 of 33

Table 4: Harmonic and anharmonic O-D stretch vibrational frequencies for water adsorbed on of H2 O/α-Al2 O3 (11¯20), one water per (2 × 2) cell. The PBE+D2/PW and the experimental results are taken from Ref.25 The scaled harmonic frequencies are given in brackets. The abbreviations stand for different adsorption modes as described in Ref.:25 iCa2=inter-CUSakO-µ2 , Cb2=CUSbkO-µ2 and iCb2=interCUSbkO-µ2 . All values are in cm−1 .

stretch mode iCa2: ODsurf iCa2: ODads Cb2: ODsurf Cb2: ODads iCb2: ODsurf iCb2: ODads

PBE+D3/AO2 ν˜ νãnh 2682 (2765) 2599 2728 (2813) 2644 1658 (1709) 1310 2769 (2855) 2687 2657 (2739) 2595 2688 (2771) 2569

B3LYP+D3/AO2 ν˜ νãnh 2773 (2784) 2695 2811 (2822) 2729 2019 (2027) 1722 2843 (2854) 2765 2778 (2799) 2694 2777 (2788) 2689

PBE+D2/PW25 ν˜ 2694 (2778) 2731 (2816) 1711 (1764) 2785 (2871) 2689 (2772) 2692 (2775)

Exp.25 ν˜ 2762 2812 2839 -

sociative) adsorption modes. The AO-based calculations use the same structural models for dissociated species as in Ref.,25 reoptimized on the respective level. The general trends here are similar to those observed for the (0001) surface: (i) the harmonic B3LYP+D3 frequencies agree well with assigned experimental frequencies, (ii) PBE+D results are not as accurate unless scaled, (iii) inclusion of anharmonicities noticeably worsens the agreement between theory and experiment.

C

Activation energies and reaction rate constants In Ref.,24 several surface reactions for H2 O/α-Al2 O3 (0001) (one molecule per

(2 × 2) cell) were studied in the framework of the transition state theory. For the current case study we have selected one reaction, namely the hydrogen transfer from the 1-4 species to 1-2, which we denote as Df-H-4-2 (“Df”=“diffusion”). The reactant and product states for this reaction as well as the transition state, optimized at the B3LYP+D3 level, are shown in Fig.2. (Geometries obtained with GGA/PW look similar.24 ) In Ref.,24 the Df-H-4-2 reaction was studied with PW91/PW using VASP, and the transition state determined with the Nudged Elastic Band (NEB) method.62 Also the reactant (1-4) was optimized on this level, and normal mode analyses for reactant and TS were performed to obtain ∆G‡ (300 K), for use in the Eyring ex-


Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: The Df-H-4-2 reaction. Left panel shows the reactant, 1-4, the right panel the product, 1-2, the middle panel the transition state of the reaction. The minima were optimized at the B3LYP+D3/AO3 level of theory, the transition state with B3LYP+D3/AO1. pression of the rate constant k(300 K). Here we also calculated in the PW format the same quantities using another GGA, namely, PBE+D2, and a hybrid functional HSE06. For the latter, expensive method the PBE+D2/PW TS and reactant structures as well as the vibrational contributions ∆ZPE +∆Gvib were used. Further, in the AO format, we applied the B3LYP+D3 and LMP2 methods. The TS structure was obtained by the TS optimization procedure as implemented in the CRYSTAL code,63 carried out at the B3LYP+D3/AO1 level. The reactant structure was optimized with the AO3 basis. The barrier was evaluated from single point calculations using the AO3 basis for B3LYP+D3 and HF, and AO3-dual for LMP2. To obtain ∆G‡ for HF and LMP2, the B3LYP+D3 corrections ∆ ZPE+∆Gvib were added to the respective ∆E ‡ . The electronic barriers, activation free energies at 300 K and the corresponding rates are compiled in Table 5. For the barriers, as expected, the use of hybrid functionals or LMP2 makes substantial differences compared to GGA. The electronic LMP2 and HSE06 barriers are mutually similar and higher than the GGA ones by about 0.15 eV. The B3LYP barrier is higher than the former by further 0.1 eV. Finally Hartree-Fock predicts a very high barrier of more than 1 eV. The activation free energies are lower than ∆E ‡ : by 0.15 eV for PBE+D3 or PW91 and by 0.11 eV for B3LYP+D3. This leads to reaction rates of the order of 108 s−1 predicted by GGA and substantially smaller rates by other methods: the hybrid functionals predict values



Page 18 of 33

Table 5: ∆E ‡ , ∆G‡ (300 K) in eV and k(300K) in s−1 for the reaction Df-H-4-2, obtained with CRYSTAL for B3LYP and LMP2 using AO basis sets and with VASP for GGA and HSE06 using plane wave bases.

∆E ‡ ∆G‡ k

B3LYP+D3/ AO31 0.69 0.58 1.2 × 103

HF/ AO32 1.17 1.06 1.0 × 10−5

LMP2/ AO3-dual2 0.60 0.49 3.7 × 104

PW91/ PW 0.42 0.27 2.0 × 108

PBE+D2/ PW 0.44 0.29 9.1 × 107

HSE06/ PW3 0.56 0.41 8.1 × 105

∆E ‡ calculated from a TS energy taken from the single-point B3LYP+D3/AO3 calculation at the B3LYP+D3/AO1 optimized geometry, and the reactant energy at the B3LYP+D3/AO3 level and the B3LYP+D3/AO3 optimized geometry. ∆ZPE and ∆Gvib contributions to ∆G‡ were calculated with B3LYP+D3/AO3. 2 ∆E ‡ was obtained from single-point calculations with the respective method (HF/AO3 or LMP2 /AO3-dual) at the B3LYP+D3/AO1-optimized TS and B3LYP+D3/AO3-optimized reactant structures with the respective method (HF/AO3 or LMP2 /AO3-dual); the ∆ZPE and ∆Gvib contributions to ∆G‡ were computed at the B3LYP+D3/AO3 level. 3 ∆E ‡ was obtained from single-point calculations at PBE+D2/PW-optimized TS and reactant geometries; the ∆ZPE and ∆Gvib contributions to ∆G‡ were computed at PBE+D2/PW level. 1

of 106 s−1 by HSE06 or 103 s−1 by B3LYP+D3, while LMP2 is in between with 104 s−1 . Hartree-Fock with a rate of 10−5 s−1 is obviously an outlier, that predicts no diffusion is taking place. The very fast GGA rates, on the other hand, are also hard to bring in line with the experimental observation of 1-4 vibrations. They may also be a reason why DFT-GGA predicts fast dissociation from the molecular state, as challenged by recent experiments.18 In this line, the hybrid DFT or LMP2 results are more reliable. In the absence of a high-level benchmark of the barrier or a measured reaction rate, it is however difficult to make a quantitative statement.

IV.

Summary and Conclusions

GGA functionals like PBE in connection with plane wave basis sets are commonly used in surface science and solid state theory in general, since they deliver qualitatively good results at comparatively low computational cost. However, certain weaknesses of this approach have been shown to be present even for “ground state properties”. In this work we examined, for one and the same system, namely water (fragments) adsorbed at an ideal α-Al2 O3 (0001) surface, the performance of GGA/plane wave as the “standard model”, for three different ground state properties of general


Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interest to surface science: Adsorption energies, adsorbate vibrational frequencies, and reaction barriers including resulting rates. The performance test was done in comparison to methods going beyond GGA and even beyond DFT: Hybrid-DFT with some portion of exact exchange, as well as WFT such as Hartree-Fock and Møller-Plesset perturbation theory. We also compared plane wave and atom-centred bases, and the importance of dispersion corrections (within DFT) or post-HF correlation energy (within WFT). Our main findings can be summarized as follows: • For the adsorption energies it was not possible to detect, whether B3LYP+D or LMP2 were providing more accurate energies than GGAs+D, as the former results were in the same ballpark and scattered by a similar magnitude as the GGA ones. Furthermore, the BSSE, which is not easy to correct for the dissociative adsorption, was also of the magnitude of this scatter. We note, however, that a hierarchical application of the WFT approach (rather than LMP2 alone) would be advantageous to DFT. For that the basis set and the method errors should be systematically corrected,32, 38, 39, 58–60 e.g., by evaluating the corresponding basis-set-limit and CCSD(T) corrections using embedded, finite clusters. • Concerning vibrations, the use of hybrid DFT leads to a noticeable improvement over GGA. For the present system, the harmonic GGA vibrational frequencies (both for plane wave and AO bases) are typically red-shifted by about 100 cm−1 compared to the corresponding (assigned) experimental values, while the hybrid DFT (B3LYP+D3 in this case) results agree with experiment much better. The GGA frequencies can however be also corrected by scaling them with empirical scaling factors.61 Interestingly, for both GGA and hybrid DFT, the inclusion of anharmonicity of the potential along vibrational coordinates worsens the agreement with experiment. The origin of this effect is not entirely clear. It may result from the environment molecules, absent in the theoretical model, whose influence could damp the anharmonicity of the potential. • The need to go beyond GGA is most apparent when it comes to reaction barriers and respective reaction rates. GGA noticeably underestimates the 19 Environment ACS Paragon Plus


barriers and overestimates rates. Hartree-Fock is another extreme in this respect: it grossly overestimates barriers virtually prohibiting the reaction. The use of hybrid functionals (in the DFT framework) or correlated wavefunction methods (in the WFT world), improves these quantities considerably. The applicability of more accurate electronic structure methods relies also on the question of gain in accuracy vs. computational effort. GGA DFT is known to be computationally very efficient, especially so in the combination with plane waves and pseudopotential approximation. At the same time evaluation of the exchange contribution and thus the hybrid DFT is computationally very costly with plane waves. For example, a VASP single-point PBE calculation for the transition state geometry (Fig.2, middle) with 400 eV energy cutoff on a 24-processor Intel-Xeon machine with 2.3 GHz clock rate took just slightly more than 7 minutes, while the corresponding HSE06 calculation was more than 36 hours. With AOs the GGA DFT calculations become somewhat more expensive, but, at the same time, evaluation of the exchange there is not much more difficult than the Coulomb part, at least for moderately sized basis sets. For example, a calculation for the same structure, carried out on the same machine, but performed with CRYSTAL with basis set AO1, took 31 minutes for PBE and 36 minutes for B3LYP. The AO1 or AO2 (triple-zeta-quality) basis sets are standard for AO-based DFT applications with CRYSTAL. For the much larger basis set AO3, the calculations become noticeably more expensive: 107 minutes for PBE and 8 hours for B3LYP, but still much faster than a hybrid-DFT calculation with PWs. Although MP2 is generally considered more expensive than DFT, current efficient implementations, which use density fitting and local correlation techniques, reduce the cost of MP2 calculations to that of (hybrid) DFT, at least for basis sets of comparable size. Unfortunately, MP2 is more sensitive to the basis set quality, so moderate basis sets, that are standardly applied in conjunction with DFT, might not be sufficient for MP2. Therefore in this work we performed the LMP2 calculations with quite extended basis sets AO3 and AO3-dual. The LMP2 calculations were conducted using the serial compilation of the CRYSCOR code, as the parallelization of OSV-based LMP2 is yet in the test-phase. A serial single-point LMP2/AO3 20 Environment ACS Paragon Plus

Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculation for the same geometry took 31 hours. The corresponding parallel calculation on 24 processors would take only 105 minutes, assuming a 75% efficiency of the parallelization, observed on the PAO-based LMP2 for 2D periodic systems.64 This is comparable to an AO-based GGA calculation. The actual bottleneck of such a calculation is clearly the HF part which was here 15 hours. In summary, hybrid functionals with moderate AO basis set when applied to adsorption systems, are only moderately more expensive than GGA calculations with plane waves and can deliver more accurate results than the latter for example for vibrational frequencies or barriers and reaction rates. MP2 is noticeably more expensive than GGA, but mainly due to the HF step, which becomes difficult with the extended (and more diffuse) basis sets, required for MP2. MP2 itself provides also a better description than GGA. However, the main advantage of the WFT / MP2 methodology is the possibility of the systematic elimination of errors, for instance, by applying ∆CCSD(T) corrections. Work along these lines can open the way to highly accurate electronic structure calculations for complex adsorbate systems.

V.

Supporting Information Supporting information contains the detailed specification of the AO basis sets

employed in this study.

VI.

Acknowledgments

We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for support of this work through Collaborative Research Center 1109 “Understanding of Metal Oxide/Water Systems at the Molecular Scale: Structural Evolution, Interfaces and Dissolution”, and under Germany’s Excellence Strategy – EXC 2008/1 (UniSysCat) – 390540038. We also thank the North-German Supercomputing Alliance (HLRN) for providing HPC resources that have contributed to the research results reported in this paper.



References 1 Perdew,

J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation

Made Simple, Phys. Rev. Lett. 1996, 77, 3865–3868. 2 Perdew,

J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.;

Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation, Phys. Rev. B 1992, 46, 6671-6687. 3 Kresse,

G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for

Metals and Semiconductors Using a Plane-Wave Basis Set, Comput. Mat. Sci. 1996, 6, 15-50. 4 Del

Re, G.; Ladik, J.; Biczø, G. Self-Consistent-Field Tight-Binding Treatment

of Polymers. I. Infinite Three-Dimensional Case, Phys. Rev. 1967, 155, 997–1003. 5 Pisani,

C.; Dovesi, R.; Roetti, C. Hartree-Fock Ab Initio Treatment of Crystalline

Solids. In Lecture Notes in Chemistry Series, Vol. 48; Springer Verlag, Berlin: 1988. 6 Stoll,

H. Correlation Energy of Diamond, Phys. Rev. B 1992, 46, 6700–6704.

7 Knab,

R.; Förner, W.; Čížek, J.; Ladik, J. Numerical Application of the Coupled

Cluster Theory with Localized Orbitals to Polymers II. Optimal Localization of Wannier Functions and the Correlation Energy in Different Approximations, J. Mol. Struct.: THEOCHEM 1996, 366, 11–33. 8 Tuma,

C.; Sauer, J. Treating Dispersion Effects in Extended Systems by Hybrid

MP2:DFT Calculations – Protonation of Isobutene in Zeolite Ferrierite, Phys. Chem. Chem. Phys. 2006, 8, 3955–3965. 9 Maschio,

L.; Usvyat, D.; Manby, F. R.; Casassa, S.; Pisani, C.; Schütz, M. Fast

Local-MP2 Method with Density-Fitting for Crystals. I. Theory and Algorithms, Phys. Rev. B 2007, 76, 075101-1–9. 10 Booth,

G.; Grüneis, A.; Kresse, G.; Alavi, A. Towards an Exact Description of

Electronic Wavefunctions in Real Solids, Nature 2013, 493, 365–370.


Page 22 of 33

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11 Becke,

A. D. A New Mixing of Hartree-Fock and Local Density-Functional The-

ories, J. Chem. Phys. 1993, 98, 1372–1377. 12 Krukau,

A.; Vydrov, O.; Izmaylov, A.; Scuseria, G. Influence of the Exchange

Screening Parameter on the Performance of Screened Hybrid Functionals, J. Chem. Phys. 2006, 125, 224106-1–5. 13 Kurita,

T.; Uchida, K.; Oshiyama, A. Atomic and Electronic Structures of

α-Al2 O3 Surfaces, Phys. Rev. B 2010, 82, 1553191-1–14. 14 Chang,

C. C. Silicon-on-Sapphire Epitaxy by Vacuum Sublimation: LEED-Auger

Studies and Electronic Properties of the Films, J. Vac. Sci. Technol. 1971, 8, 500–511. 15 Tsyganenko,

A. A.; Mardilovich, P. P. Structure of Alumina Surfaces, J. Chem.

Soc., Faraday Trans. 1996, 92, 4843–4852. 16 Lee,

W. E.; Lagerlof, K. P. D. Structural and Electron Diffraction Data for

Sapphire (α-Al2 O3 ), J. Electron Microsc. Tech. 1985, 2, 247–258. 17 Kelber,

J. Alumina Surfaces and Interfaces Under Non-Ultrahigh Vacuum Con-

ditions, Surf. Sci. Repts. 2007, 62, 271–303. 18 Petrik,

N.; Huestis, P.; LaVerne, J.; Aleksandrov, A.; Orlando, T.; Kimmel, G.

Molecular Water Adsorption and Reactions on α-Al2 O3 (0001) and α-Alumina Particles, J. Phys. Chem. C 2018, 122, 9540–9551. 19 Kirsch,

H.; Wirth, J.; Tong, Y.; Wolf, M.; Saalfrank, P.; Campen, R. K.

Experimental Characterization of Unimolecular Water Dissociative Adsorption on α-Alumina, J. Phys. Chem. C 2014, 118, 13623-13630. 20 Huang,

P.; Pham, T. A.; Galli, G.; Schwelger, E. Alumina(0001)/Water Inter-

face: Structural Properties and Infrared Spectra from First-Principles Molecular Dynamics Simulations, J. Phys. Chem. C 2014, 118, 8944–8951. 21 Hass,

K. C.; Schneider, W. F.; Curioni, A.; Andreoni, W. First-Principles

Molecular Dynamics Simulations of H2 O on α-Al2 O3 (0001), J. Phys. Chem. B 2000, 104, 5527–5540.



22 Ranea,

V. A.; Carmichael, I.; Schneider, W. F. DFT Investigation of Intermediate

Steps in the Hydrolysis of α-Al2 O3 (0001), J. Phys. Chem. C 2009, 113, 2149– 2158. 23 Thissen,

P.; Grundmeier, G.; Wippermann, G.; Schmidt, S. Water Adsorption

on the α-Al2 O3 (0001) Surface, Phys. Rev. B 2009, 80, 245403. 24 Wirth,

J.; Saalfrank, P. The Chemistry of Water on α-Alumina: Kinetics and

Nuclear Quantum Effects from First Principles, J. Phys. Chem. C 2012, 116, 26829–26840. 25 Heiden,

S.; Yue, Y.; Kirsch, H.; Wirth, J.; Saalfrank, P.; Campen, R. K.

Water Dissociative Adsorption on α-Al2 O3 (11¯20) is Controlled by Surface Site Undercoordination, Density, and Topology, J. Phys. Chem. C 2018, 122, 6573– 6584. 26 Zhao,

Y.; González-García, N.; Truhlar, D. Benchmark Database of Barrier

Heights for Heavy Atom Transfer, Nucleophilic Substitution, Association, and Unimolecular Reactions and its Use to Test Theoretical Methods, J. Phys. Chem. A 2005, 109, 2012-2018. 27 Cora,

F.; Alfredsson, M.; Mallia, G.; Middlemiss, D. S.; Mackrodt, W. C.;

Dovesi, R.; Orlando, R. The Performance of Hybrid Density Functionals in Solid State Chemistry. In Principles and Applications of Density Functional Theory in Inorganic Chemistry II, Vol. 113; Springer, Berlin: 2004. 28 Grimme,

S. Semiempirical GGA-Type Density Functional Constructed with a

Long-Range Dispersion Correction, J. Comput. Chem. 2006, 27, 1787–1799. 29 Tkatchenko,

A.; Scheffler, M. Accurate Molecular van der Waals Interactions

from Ground-State Electron Density and Free Atom Reference Data, Phys. Rev. Lett. 2009, 102, 0730051-1–4. 30 Hermann,

A.; Schwerdtfeger, P. Ground-State Properties of Crystalline Ice from

Periodic Hartree-Fock Calculations and a Coupled-Cluster-Based Many-Body Decomposition of the Correlation Energy, Phys. Rev. Lett. 2008, 101, 183005-1–4.


Page 24 of 33

Page 25 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


31 Heßelmann,

A. Improved Supermolecular Second Order Møller-Plesset Inter-

molecular Interaction Energies Using Time-Dependent Density Functional Response Theory, J. Chem. Phys. 2008, 128, 144112-1–9. 32 Tsatsoulis,

T.; Hummel, F.; Usvyat, D.; Schütz, M.; Booth, G. H.; Binnie, S. S.;

Gillan, M. J.; Alfe, D.; Michaelides, A.; Grüneis, A. A Comparison Between Quantum Chemistry and Quantum Monte Carlo Techniques for the Adsorption of Water on the (001) LiH Surface, J. Chem. Phys. 2017, 146, 204108-1–9. 33 Kosata,

J.; Merkl, P.; Teeratchanan, P.; Hermann, A. Stability of Hydrogen

Hydrates from Second-Order Møller-Plesset Perturbation Theory, J. Phys. Chem. Lett. 2018, 13, 5624–5629. 34 Pisani,

C.; Maschio, L.; Casassa, S.; Halo, M.; Schütz, M.; Usvyat, D. Periodic

Local MP2 Method for the Study of Electronic Correlation in Crystals: Theory and Preliminary Applications, J. Comp. Chem. 2008, 29, 2113-2124. 35 Pisani,

C.; Schütz, M.; Casassa, S.; Usvyat, D.; Maschio, L.; Lorenz, M.;

Erba, A. Cryscor: A Program for the Post-Hartree-Fock Treatment of Periodic Systems, Phys. Chem. Chem. Phys. 2012, 14, 7615–7628. 36 Ben,

M. D.; Hutter, J.; VandeVondele, J. Forces and Stress in Second Or-

der Møller-Plesset Perturbation Theory for Condensed Phase Systems within the Resolution-of-Identity Gaussian and Plane Waves Approach, J. Chem. Phys. 2015, 143, 102803-1–22. 37 Booth,

G.; Tsatsoulis, T.; Chan, G. K.-L.; Grüneis, A. From Plane Waves to

Local Gaussians for the Simulation of Correlated Periodic Systems, J. Chem. Phys. 2016, 145, 084111-1–10. 38 McClain,

J.; Sun, Q.; Chan, G. K.-L.; Berkelbach, T. C. Gaussian-Based

Coupled-Cluster Theory for the Ground-State and Band Structure of Solids, J. Chem. Theory Comput. 2017, 13, 1209–1218. 39 Gruber,

T.; Liao, K.; Tsatsoulis, T.; Hummel, F.; Grüneis, A. Applying

the Coupled-Cluster Ansatz to Solids and Surfaces in the Thermodynamic Limit, Phys. Rev. X 2018, 8, 021043-1–6.



40 Usvyat,

Page 26 of 33

D.; Maschio, L.; Schütz, M. Periodic Local MP2 Method Employing

Orbital Specific Virtuals, J. Chem. Phys. 2015, 143, 102805-1–12. 41 Kresse,

G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals, Phys.

Rev. B 1993, 47, 558-561. 42 Kresse,

G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-

Metal – Amorphous-Semiconductor Transition in Germanium,

Phys. Rev. B

1994, 49, 14251–14269. 43 Kresse,

G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy

Calculations using a Plane-Wave Basis Set, Phys. Rev. B 1996, 54, 11169–11186. 44 Kresse,

G.;

Joubert, D.

From Ultrasoft Pseudopotentials to the Projector

Augmented-Wave Method, Phys. Rev. B 1999, 59, 1758–1775. 45 Dovesi,

R.;

Orlando, R.;

Erba, A.;

Zicovich-Wilson, C. M.;

Civalleri, B.;

Casassa, S.; Maschio, L.; Ferrabone, M.; Pierre, M. D. L.; D’Arco, P.; Noël, Y.; Causà, M.; Rerat, M.; Kirtman, B. CRYSTAL14: A Program for the Ab Initio Investigation of Crystalline Solids, Int. J. Quantum Chem. 2014, 114, 1287–1317. 46 Pascale,

F.; Zicovich-Wilson, C. M.; Lopez, F.; Civalleri, B.; Orlando, R.;

Dovesi, R. The Calculation of the Vibration Frequencies of Crystalline Compounds and its Implementation in the CRYSTAL Code, J. Comput. Chem. 2004, 25, 888–897. 47 Zicovich-Wilson,

C.;

Pascale, F.;

Roetti, C.;

Saunders, V.;

Orlando, R.;

Dovesi, R. The Calculation of the Vibration Frequencies of α-Quartz: The Effect of Hamiltonian and Basis Set, J. Comput. Chem. 2004, 25, 1873–1881. 48 Tosoni,

S.; Pascale, F.; Ugliengo, P.; Orlando, R.; Saunders, V. R.; Dovesi, R.

Vibrational Spectrum of Brucite, Mg(OH)2 : A Periodic Ab Initio Quantum Mechanical Calculation Including OH Anharmonicity, Chem. Phys. Lett. 2004, 396, 308–315. 49 Blöchl,

P. E. Projector Augmented-Wave Method, Phys. Rev. B 1994, 50, 17953–

17979.


Page 27 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50 Perdew,

J. P.; Burke, K.; Ernzerhof, M. Erratum: Generalized Gradient Ap-

proximation Made Simple, Phys. Rev. Lett. 1997, 78, 1396. 51 Grimme,

S.; Antony, J.; Ehrlich, S.; Goerigk, L. A Consistent and Accurate Ab

Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu, J. Chem. Phys. 2010, 132, 154104-1–19. 52 Stephens,

P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Cal-

culation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Theory, J. Phys. Chem. 1994, 98, 11623–11627. 53 Suhai,

S. Quasiparticle Energy-Band Structures in Semiconducting Polymers:

Correlation Effects on the Band Gap in Polyacetylene, Phys. Rev. B 1983, 27, 3506–3518. 54 Peintinger,

M. F.; Oliveira, D. V.; Bredow, T. Consistent Gaussian Basis Sets

of Triple-Zeta Valence with Polarization Quality for Solid-State Calculations, J. Comput. Chem. 2012, 34, 451–459. 55 Dunning,

T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations.

I. The Atoms Boron Through Neon and Hydrogen, J. Chem. Phys. 1989, 90, 1007–1023. 56 Usvyat,

D.; Maschio, L.; Pisani, C.; Schütz, M. Second Order Local Møller-

Plesset Perturbation Theory for Periodic Systems: The CRYSCOR Code, Z. Phys. Chem. 2010, 224, 441–454. 57 Eyring,

H. The Activated Complex and the Absolute Rate of Chemical Reactions,

Chem. Rev. 1935, 17, 65-77. 58 Müller,

C.; Usvyat, D. Incrementally Corrected Periodic Local MP2 Calculations:

I. The Cohesive Energy of Molecular Crystals, J. Chem. Theory Comput. 2013, 9, 5590–5598. 59 Boese,

A.; Saalfrank, P. CO Molecules on a NaCl(100) Surface: Structures,

Energetics and Vibrational Davydov Splittings, J. Phys. Chem. C 2016, 120, 12637–12653.



60 Alessio,

M.; Bischoff, F. A.; Sauer, J. Chemically Accurate Adsorption Energies

for Methane and Ethane Monolayers on the MgO(001) Surface, Phys. Chem. Chem. Phys. 2018, 20, 9760–9769. 61 Kershawani,

M.; Brauer, B.; Martin, J. Frequency and Zero-Point Vibrational

Energy Scale Factors for Double-Hybrid Density Functionals (and Other Selected Methods): Can Anharmonic Force Fields be Avoided?, J. Chem. Phys. A 2015, 119, 1701-1714. 62 Jónsson,

H.; Mills, G.; Jacobsen, K. W. Classical and Quantum Dynamics in

Condensed Phase Simulations; World Scientific: 1995. 63 Zicovich-Wilson,

C. M.; San-Roman, M. L.; Ramirez-Solis, A. Mechanism of F−

Elimination from Zeolitic D4R Units: A Periodic B3LYP Study on the Octadecasil Zeolite, J. Phys. Chem. C 2010, 114, 2989–2995. 64 Maschio,

L. Local MP2 with Density Fitting for Periodic Systems: A Parallel

Implementation, J. Chem. Theory Comput. 2011, 7, 2818–2830.

TOC Graphic

1−4

1−2 TS MP2 GGA


Page 28 of 33

Page 29 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34


1−4

1−2 TS MP2 GGA

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

CUS

Page 30 of 33

CUS 4 2

CUS

CUS


Page 31 of 33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18


Page 32 of 33

Page 33 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18



Theoretical Surface Science Beyond Gradient-Corrected Density

Recommend Documents