TiO2 Anatase Surface - American Chemical Society

Dec 17, 2015 - Fabrizio Gala,. †. L. Agosta,. ‡ and Giuseppe Zollo*,†. †. Dipartimento di Scienze di Base e Applicate per l'Ingegneria (Sezion...
0 downloads 0 Views 2MB Size
Subscriber access provided by CMU Libraries - http://library.cmich.edu

Article 2

Water Kinetics and Clustering on the (101) TiO Anatase Surface Fabrizio Gala, Lorenzo Agosta, and Giuseppe Zollo J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b10934 • Publication Date (Web): 17 Dec 2015 Downloaded from http://pubs.acs.org on December 20, 2015

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry C 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 20

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

Water Kinetics and Clustering on the (101) TiO2 Anatase Surface Fabrizio Gala,† L. Agosta,‡ and Giuseppe Zollo∗,† †Dipartimento di Scienze di Base e Applicate per l’Ingegneria (Sezione di Fisica), Universit´ a di Roma “La Sapienza”, Via A. Scarpa 14–16, 00161 Rome Italy ‡Department of Materials and Environmental Chemistry, Stockholm university, 2014 Arrhenius Laboratory, Svante Arrhenius v¨ ag 16C, Stockholm Sweden E-mail: [email protected]

1

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 49 50 51 52 53 54 55 56 57 58 59 60

Abstract (101) anatase TiO2 surface in water ambient is an important system for the interaction of biocompatible nano-devices with biological environment. Following the experimental evidences showing that water molecules are mobile at as low temperature as 190 K and tend to form clusters along the [11¯1]/[1¯1¯1] surface directions, a complete theoretical characterization of the dynamical properties of the first water layer on the (101) anatase TiO2 surface is presented. A variety of computational techniques have been employed in the context of the transition state theory in the harmonic regime, ranging from first-principles total energy ground state calculations, density functional perturbation theory, minimum energy path search, kinetic Monte Carlo simulations, in order to explain the experimental results on water kinetics on the (101) anatase TiO2 surface. We have calculated the migration energy barrier of water molecules, the vibrational pre-factor through the phonon density of states, the hopping rate along two principal directions. Lastly, in a kinetic Monte Carlo context, we have simulated and clarified the dynamical processes that are on the basis of the observed experimental behavior.

Keywords Transition state theory, Density Functional Theory; kinetic Monte Carlo; water migration; TiO2 surface hydration.

2

ACS Paragon Plus Environment

Page 2 of 20

Page 3 of 20

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

1

Introduction

Titanium dioxide (TiO2 ) surfaces are employed in many technological area such as medicine and pharmacology, 1 photocatalysis and solar-hydrogen production, 2 surface recognition by organic materials. 3 In particular, while the macroscopic TiO2 ground state phase is rutile, TiO2 nano-particles prefer the anatase phase. As a consequence, because TiO2 nano-devices employed for most of the above applications involve aqueous environment, the interaction between water and TiO2 (101) anatase surface has attracted strong interest and has been extensively studied in the last years. 4–9 However such studies have been mainly focussed on the adsorption properties of water on the titania surface, while little is known on the kinetic properties. Furthermore, it has been shown that water molecules mediate the adsorption of biological molecules on TiO2 (101) surface; 3,10,11 thus selectivity properties of biological matter on TiO2 (101) surface might be strongly affected by the diffusion properties of water on such surfaces. From the experimental point of view, it has been shown 7 that water adsorbed on anatase surface is thermodynamically stable in vacuum up to ∼ 200 K, even if in presence of subsurface defects water adsorption on anatase has been observed up to 400 K. 8 Moreover a tendency to form small water aggregates along the [11¯1]/[1¯1¯1] direction of the (101) anatase surface has been observed by STM 8 that is not fully explained so far. In the following we report on a complete characterization of the kinetic behavior of the first water layer deposited onto the (101) anatase surface giving, for the first time, a complete insight to the above mentioned experimental observation of water clustering along preferential surface directions. The study has been performed in the context of quantum mechanical transition state theory (TST) in the harmonic regime using ground state Density Functional Theory (DFT), Density functional perturbation theory (DFPT), minimum energy path (MEP) search and kinetic Monte Carlo (KMC) simulations.

3

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 49 50 51 52 53 54 55 56 57 58 59 60

2

Page 4 of 20

Theoretical Method

Ab-initio atomistic modeling based on DFT 12,13 has been already used to model the water adsorption on various TiO2 surfaces both of the rutile and the anatase phases, 5,8,14–16 with a good agreement with the experimental values revealing that water adsorbs molecularly through the formation of a strong bond to a fivefold coordinated surface Ti atom (Ti5c ) atom and weak hydrogen bonds with twofold coordinated surface O atoms (O2c ) . The DFT scheme here employed adopts a generalized gradient approximation (GGA) of the electron exchange and correlation energy using Perdew-Burke-Ernzerhof formula (PBE) 17 as implemented in QUANTUM-ESPRESSO (QE) package. 18 Due to computational reasons, ultra-soft pseudopotentials (US-PP) have been employed being suitable for this system, according also to other authors. 14 US-PPs have allowed the usage of an energy cut-off of 60 Ry for the wave functions and 400 Ry for the electron density in the context of a plane wave expansion basis set. The (101) anatase TiO2 surface has been modeled using a slab geometry laying in the xy plane composed by two surface unit cells along the y and the z directions and one along x. The resulting slab contains 48 TiO2 atoms (being made of two TiO2 molecular layers along z) and the supercell includes two vacuum regions nearly of 16.2 ˚ A thick each. The artificial electric field across the slab induced by the PBCs has been corrected following Ref. 19 and a 4x4x1 Monkhorst-Pack (MP) k-point grid 20 for the Brillouin zone (BZ) sampling has been employed. The ground state configurations have been fully relaxed using the BroydenFletcher-Goldfarb-Shanno (BFGS) algorithm 21 with the Hellmann-Feynman forces acting on the ions, together with empirical long-range corrections. 22 TST predicts the hop rate ν between two local minima as: s kB T Zvib ν= exp h Zvib



∆E − kB T



(1)

s where kB and h are the Boltzmann and the Planck constants, Zvib and Zvib are the partition

functions of the initial and transition (saddle) states and ∆E is the energy difference between 4

ACS Paragon Plus Environment

Page 5 of 20

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

the saddle-point and the minimum energy configurations. ∆E energy barriers and the corresponding diffusion paths, or minimum energy paths, have been studied through the Nudged Elastic Band (NEB) method 23,24 as implemented in QE. On the other hand, quantum statistical mechanics in harmonic approximation allows to rewrite the hop rate defined in Eq.(1) as a function of the Helmholtz free energy F of the system, as: kB T ν= exp h



∆E + ∆Fvib (T ) − kB T



= ν exp ∗



Ea − kB T



(2)

being    Z 3N X kB T ~ωi (q) Fvib (T ) = dq = log 2 sinh Ω BZ i=1 2kB T     Z ∞ ~ω ~ω + kB T log 1 − exp − = dωg(ω) 2 kB T 0

(3)

where Ea = ∆E + ∆Fvib (0) is the activation energy. Due to the slab geometry employed, the integration over the Brillouin zone is restricted to a bi-dimensional space. Hence Ω is the area of the Brillouin Zone (BZ) in such a bi-dimensional space, N is the total number of atoms, and ωi (q) is the frequency of the i-th mode at wave-vector q and g(ω) is the phonon density of states (PHDOS) normalized to the number of phonon branches (i.e. 3N ). The temperature dependent prefactor ν ∗ is then defined as:

ν



kB T = exp h





Z



∞ 0

dω∆g(ω) log 1 − exp



~ω − kB T



In the vibrational partition function of the saddle configuration one of such branches, corresponding to the reaction coordinate, is a pure imaginary number, and must be excluded in the evaluation of the free energy; thus in the following it will be treated as a negative frequency, being excluded from the integral defining Fvib (T ). 5

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 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 20

Phonon frequencies ωi (q), evaluated in the framework of DFPT 25 with QE, are solutions of the secular equation: 1 st 2 Cαβ (q) − ω (q)1 = 0 det √ MM s

(4)

t

st where Cαβ (q) are the interatomic force constants (IFCs) between atoms s and t (of mass

Ms and Mt respectively) evaluated in reciprocal space, 25 and α, β run over the cartesian indexes {x, y, z}; a 2x2x1 q point uniform grid has been employed to calculate the entire phonon spectra: indeed such mesh has been tested to be sufficient for convergence. All the atoms have been involved in the phonon calculations, none of them being kept fixed. The semi-empirical term accounting for long-range interactions has been neglected in the calculation of the phonon spectra because its contribution to the IFCs, evaluated analytically, resulted negligible for all the configurations studied. Lastly, Kinetic Monte Carlo (KMC) simulations have been performed to follow the diffusion and clustering events of water molecules on the (101) anatase surface at various temperature and coverage values. The rates used in the KMC simulations have been defined on the basis of the migration energy and hopping rates calculated as specified above and corrected for local effects as discussed in the following section.

3

Results and Discussion

As mentioned above, the adsorption of a single water molecule on a (101) TiO2 surface occurs through three bonds: one Ow − T i5c and two Hw − O2c hydrogen bonds (H-bonds). 3 The adsorption energy of such ground state configuration (”stable configuration”) is equal to Eads = -0.73 eV (see Fig.1(a)). This fact suggests to look for diffusion paths connecting Ti5c equivalent surface atoms, i.e. the adsorption sites where a water molecule can be adsorbed; following this criterion, two migration paths have been chosen, namely the ones along the [010] and [11¯1]/[1¯1¯1] directions parallel to the surface; they basically connect two adjacent 6

ACS Paragon Plus Environment

Page 7 of 20

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

adsorption sites where the water can stay in its ground state configuration. In Figs. 1 and Fig. 2 the energy barriers and the corresponding configurations are reported for respectively the [010] and the [11¯1]/[1¯1¯1] diffusion paths with the relevant error bars.

Figure 1: (Color online) MEP and total energy profile for water diffusion along the [010] direction on the TiO2 surface, together with the corresponding configurations.

Figure 2: (Color online) MEP and total energy profile for water diffusion along the [11¯1]/[1¯1¯1] direction on the TiO2 surface, together with the corresponding configurations. The total energy profiles here calculated, exhibit a symmetric and an asymmetric barrier for [010] and [11¯1]/[1¯1¯1] respectively, similarly to previous ones obtained with a semi7

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 49 50 51 52 53 54 55 56 57 58 59 60

empirical model; 26 the energy barriers here measured for migration along the [010] and the [11¯1]/[1¯1¯1] paths are respectively ∆E[010] ≃ 0.52 eV and ∆E[11¯1]/[1¯1¯1] ≃ 0.57 eV, indicating that the corresponding ones calculated by semi-empirical methods 26 are severely overestimated by 130 and 150 meV respectively.

Figure 3: (Color Online) Calculated PHDOS (convoluted with gaussian functions) of the anatase (101) surface (grey area), of the stable configuration (black line), and of the two saddle configurations studied (red and orange lines respectively) (a); variation of the vibrational free energy as a function of the temperature for both the saddle configurations (b); frequency pre-factor for the two directions of migrating water (c). In both cases TiO2 surface atoms are little affected by relaxation and the main contribution to the calculated energy barrier is due to bond breaking and forming between the diffusion water and the surface. Along [010], water shifts almost rigidly and the largest contribution to the energy barrier is due to the breaking of the Ow − T i5c bond, while at least one hydrogen bond is maintained along the entire migration path. The saddle point configuration of water migration along [010] (”[010] saddle configuration”) is such that the water oxygen stays in the middle point between the initial and the final Ti5c adsorption sites of the diffusion step. Along [11¯1]/[1¯1¯1], instead, a sort of ”rolling stone” behavior of the water molecule occurs; 27 again the covalent bond breaking contributes most to the energy barrier, but in this case we can notice that, during the first half, the molecule rotates about the axis that contains the two hydrogens (bonded to the surface through the Hw − O2c H-bonds), while in the 8

ACS Paragon Plus Environment

Page 8 of 20

Page 9 of 20

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

second half of the path the water molecule rotates about the axis containing the oxygen atom that is bonded to the Ti5c atom. The saddle point (the ”[11¯1]/[1¯1¯1] saddle configuration” shown in Fig.2(c) ) is such that the H2 O molecule lies in the plane orthogonal to the surface. Following TST, the hopping rate of the water molecules along the two diffusion paths depends on the pre-factor ν ∗ that is related to the vibrational properties of both the saddle and the stable configurations. The PHDOS of the stable configuration, reported in Fig.3(a) shows that the adsorption of the water molecule affects only slightly the TiO2 surface g(ω) (that is in the whole range [0 ≤ ω ≤ 27 T Hz]) in the range [15 T Hz ≤ ω ≤ 25 T Hz]; moreover three local water vibrational peaks are visible outside the surface region. Those peaks (corresponding to the modes 151, 152, and 153 with frequency values reported in Table 1) are respectively a bending mode and a symmetrical and an asymmetrical stretching modes of the O-H bonds of the H2 O molecule. In particular the 153 peak is down-shifted by ∼ 1.53 THz with respect to the isolated water molecule. In the [010] saddle configuration, a new imaginary eigenvalue of the secular equation (an imaginary frequency peak) appears at 2.08i THz; such imaginary solution corresponds to the reaction coordinate of the diffusion process and shows no dispersion in the BZ. The other localized peaks are still present, even if they are shifted (see Table 1), and the 152 mode is an asymmetrical stretching mode instead of a symmetrical one. This circumstance is due to the breaking of the local symmetry of the system at the saddle point where the water molecule is rotated in the (101) plane (see Fig.1(d)). Similarly, the PHDOS of the [11¯1]/[1¯1¯1] saddle configuration has an imaginary eigenvalue at 1.60i THz (along the reaction coordinate), and three localized modes outside the frequency region of the surface; although the molecule has a symmetry axis orthogonal to the surface (see Fig.2)(c)), the characters of the modes are the same as the stable configuration ones while the frequencies undergo some shift with respect to the stable configuration. On the basis of the above PHDOS, the change of the vibrational free energy of the saddle

9

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 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 20

point with respect to the stable configuration ∆Fvib (T ) have been calculated and are reported in Fig.3(b); the zero point energy (ZPE) differences (i.e. ∆Fvib (0) ) are quite small and negative (see Table 2) being almost negligible with respect to the energy barriers ∆E (first column of Table 2). The resulting vibrational prefactors ν ∗ shown in Fig.3(c) are rather different: while in the [010] diffusion direction ν ∗ is already nearly saturated at T=150K (being withing 5% of the saturation value that is obtained at T=200 K), in the [11¯1]/[1¯1¯1] case ν ∗ increases monotonically in a much wider range containing the room temperature (RT); hence the vibrational prefactor at RT results four times larger for [11¯1]/[1¯1¯1] than for [010] (see Table 2). This circumstance makes the rate along [010] to behave nearly as a simple Arrhenius function of the temperature (with a constant prefactor) in the range 150K ≤ T ≤ 300K where the harmonic approximation is likely to hold, with the saturated pre-factor value being ν∗ = 3.64T Hz (see Table 2). In any case the two pre-factors are of the same order of magnitude. This is the reason why the rates are mainly dependent on the activation energies along the two directions considered. Table 1: Frequencies and classification of normal modes for isolated frequencies for the configurations studied. Modes are classified as: bending(B), asymmetrical stretching (aS), and symmetrical √ stretching (sS). The symbol i refers to the imaginary unit and it is equal to −1. conf. isolated water stable [010] saddle [11¯1]/[1¯1¯1] saddle

Imaginary Mode [THz] 2.08i 1.60i

Mode 151 [THz] 47.93 (B) 47.34 (B) 47.24 (B) 47.76 (B)

Mode 152 [THz ] 108.72 (sS) 108.27 (sS) 107.60 (aS) 110.27 (sS)

Mode 153 [THz] 112.031 (sS) 110.50 (aS) 114.12 (aS) 113.43 (aS)

Indeed the final calculated jump frequencies at room temperature, that are of the same order of magnitude, are such that ν[010] = 7.82 · 103 Hz is greater than ν[11¯1]/[1¯1¯1] = 3.94 · 103 Hz, due to the larger diffusion energy barriers along [11¯1]/[1¯1¯1] (with respect to the [010] one) affecting the exponential factor in Eq. (2). On the basis of the calculated values of the pre-factor ν ∗ and of the activation energy Ea , we have implemented a kinetic Monte Carlo model in order to explain the preferential 10

ACS Paragon Plus Environment

Page 11 of 20

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

Table 2: Diffusion barriers, variations in the ZPE, vibrational prefactors and jump frequencies for the [010] and [11¯1]/[1¯1¯1] diffusion directions. Prefactors and jump frequencies have been calculated at room temperature (i.e. T= 300 K). direction

∆E (eV)

∆Fvib (0) (eV)

∗ ν(300) (THz)

ν(300) (kHz)

[010] [11¯1]/[1¯1¯1]

0.52 0.57

-4.04 ·10−3 -7.20 ·10−3

3.64 11.23

7.82 3.94

water clustering along the [11¯1]/[1¯1¯1] directions observed by STM. 8 We have considered six possible jumping events for each isolated water molecule: two towards the two nearest neighbours Ti adsorption sites (distance d = 3.47 ˚ A) located along the [010] surface direction (occurring with the relevant rate along [010]) and four towards the next nearest neighbours Ti surface adsorption sites (d = 5.491 ˚ A) located along the two sets of [11¯1]/[1¯1¯1] type directions (occurring with the relevant rate for the [11¯1]/[1¯1¯1] direction). In case a given nearest or the next nearest neighbour was occupied, the relevant jumping event would be inhibited so that the total number of events would be reduced by one. The activation energy must be modified depending on the local environment of the migrating water molecule; an exact treatment would require the explicit calculation of all the possible migration barriers depending on the local environment; however this is unattainable with reasonable computational resources because the scenario under study is rather complicated. Therefore we have decided to treat the local environment using the ground state data reported in the same article by He and co-workers. 8 Indeed it can be easily shown that, adopting the isolated water molecule as a reference for the total energy, the migration energy barrier can be corrected for local effects using the total energy values reported in that article. Hence, consider, for instance, the case of a migrating water that stays close (at the nearest neighbour Ti site) to another water along the [010] row; then, according to He and coworkers, the migration barrier obtained for an isolated water molecule should be reduced by 66 meV (i.e. 33 meV per water molecule). Following this choice, we have corrected the relevant energy barrier (i.e. the [010] or the [11¯1]/[1¯1¯1] ones) by an additive term depending on the local environment: all the possible local environments have been

11

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 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 12 of 20

Page 13 of 20

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

In Fig.4(a) it is shown the number of water pairs, i.e. the pairs of waters staying at nearest (along the [010] direction) or next nearest (along the two [11¯1]/[1¯1¯1]-type directions) neighbours during time. It can be seen that, after the initial transient of about 100 s, the average number of waters pairs aligned along the [11¯1]/[1¯1¯1] directions is about 60, value that is kept nearly constant during the simulation. On the contrary the number of waters aligned along the [010] direction is negligible. We have also measured the average length of the clusters aligned along the [11¯1]/[1¯1¯1] direction that is shown in Fig. 4(b). It is seen that the average length (the value shown has been further averaged along a time interval of 1000 s) is about ¯l[11¯1]/[1¯1¯1] ≈ 6.5 ˚ A which means that about 18% of all the alignments measured along [11¯1]/[1¯1¯1] include three waters while about 82% of water pairs are isolated. These results are consistent with the experimental results by He and co-workers that have observed the occurrence of alignments along the [11¯1]/[1¯1¯1] directions by sequential STM images; in our simulations, we have no evidence of water clusters with more than three waters aligned along the same [111]-type direction at T = 190K and 0.11 ML of coverage. The STM images reported by He and co-workers were taken subsequently but it was not specified the exact time-scale of the STM pictures. Then we have calculated the average time latency of the water pairs aligned either along the [11¯1]/[1¯1¯1] or the [010] directions, as shown in Fig 4(c). The time latency at t is the average life time of water pairs alignments calculated in the time interval [0, t]. After an initial transient, the average latency of water pairs along the [11¯1]/[1¯1¯1] saturates at the value of tlat = 10s while the one measured for the alignment along the [010] is two orders of magnitude smaller. This means that the time needed to observe water migrating on the surface should be at least 10 s in these conditions. A typical snapshot obtained during the KMC simulation is reported in Fig. 4(d) where two clusters with three waters are indicated. In the same article, He and co-workers reported on the formation of a water superstructure at lower temperature (150 K) and larger coverage (0.24 ML). We have replicated this experiment in the context of KMC showing the frequent formation of water superstruc-

13

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 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 14 of 20

Page 15 of 20

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

along [11¯1]/[1¯1¯1] (the average length is measured as about 7.5 ˚ A); moreover, these clusters are likely to stay at [11¯1]/[1¯1¯1] next nearest rows because of the large surface density and because the clustering at adjacent rows is inhibited for energetic reasons (the energy penalty for the alignment along [010] is quite large). Then, while the clustering along [11¯1]/[1¯1¯1] is explained on the basis of the different diffusion properties and energetic penalty values for the alignments along [11¯1]/[1¯1¯1] with respect to [010], the formation of the super-stuctured water patterns observed at T=150 K and 0.24 ML of coverage is also strictly related to the large coverage ratio used in the experiments and to the large latency that makes such structures stable.

4

Conclusions

To summarize, in the context of quantum-mechanical TST and in connection with harmonic approximation, we have employed ab-initio density-functional calculations to evaluate the prefactors and the activation energies for water diffusion on a defect free (101) TiO2 surface along the [010] and the [11¯1]/[1¯1¯1] directions. We have found that the pre factor of the diffusion rate along the [010] direction saturates at lower temperature (nearly T=150 K) than the one along the [11¯1]/[1¯1¯1] that, instead, increases monotonically in the same temperature range. The two pre-factors, are, however, of the same order of magnitude. We have also found that the activation energy for migration along [010] is nearly 50 meV lower than the one along [11¯1]/[1¯1¯1] that makes the diffusion along [010] more favorable. On the basis of the water molecule jump rates so calculated, the migration/aggregation processes of water belonging to the first layer on the (101) anatase TiO2 surface have been simulated in the context of kinetic Monte Carlo method. Our results clarify the experimental results reported in the recent literature showing preferential aggregation along the [11¯1]/[1¯1¯1] direction for low coverage (0.11 ML) at T=190 K. Moreover the formation of water superstructures previously evidenced by STM at T = 150K and 0.24 ML of coverage is definitely attributed to the

15

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 49 50 51 52 53 54 55 56 57 58 59 60

formation of parallel alignments along the [11¯1]/[1¯1¯1] with one empty row in between. Such super structures are quite stable at the chosen temperature and coverage which means that stable water patterning of the (101) TiO2 anatase surface can be obtained in these conditions.

5

Acknowledgments

The computing resources for this work have been provided by CINECA under the HPC Grant Project HP10CPAIB5 and by CRESCO/ENEAGRID High Performance Computing infrastructure and its staff. 28 CRESCO/ENEAGRID High Performance Computing infrastructure is funded by ENEA, the Italian National Agency for New Technologies, Energy and Sustainable Economic Development and by Italian and European research programs, see http://www.cresco.enea.it/english for information. Part of this research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The authors are also grateful to Dr. Luigi Bagolini for the many useful discussions about the topics of the present work.

References (1) Shiba, K. Exploitation of Peptide Motif Sequences and their Use in Nanobiotechnology. Curr. Opin. Biotechnol. 2010, 21, 412–25. (2) Fujishima, A.; Zhang, X.; Trykc, D. A. TiO2 Photocatalysis and Related Surface Phenomena. Surf. Sci. Rep. 2008, 63, 515–582. (3) Agosta, L.; Zollo, G.; Arcangeli, C.; Buonocore, F.; Gala, F.; Celino, M. Water Driven Adsorption of Amino Acids on the (101) Anatase TiO2 Surface: an Ab Initio Study. PCCP 2015, 17, 1556–1561. 16

ACS Paragon Plus Environment

Page 16 of 20

Page 17 of 20

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

(4) Zhao, Z.; Li, Z.; Zou, Z. Understanding the Interaction of Water with Anatase TiO2 (101) Surface from Density Functional Theory Calculations. Physics Letters A 2011, 375, 2939–2945. (5) Sun, C.; Liu, L.-M.; Selloni, A.; Lu, G. Q. M.; Smith, S. C. Titania-Water Interactions: a Review of Theoretical Studies. J. Mater. Chem. 2010, 20, 10319–10334. (6) Vittadini, A.; Casarin, M.; Selloni, A. Chemistry of and on TiO2 -Anatase Surfaces by DFT Calculations: a Partial Review. Theor. Chem. Acc. 2007, 117, 663–671. (7) Herman, G. S.; Dohnhalek, Z.; Ruzycki, N.; Diebold, U. Experimental Investigation of the Interaction of Water and Methanol with Anatase-TiO2 (101). J. Phys. Chem. B 2003, 107, 2788. (8) He, Y.; Tilocca, A.; Dulub, O.; Selloni, A.; Diebold, U. Local Ordering and Electronic Signatures of Submonolayer Water on Anatase TiO2 (101). Nat. Mater. 2009, 8, 585– 589. (9) Aschauer, U.; He, Y.; Cheng, H.; Li, S.-C.; Diebold, U.; Selloni, A. Influence of Subsurface Defects on the Surface Reactivity of TiO2 : Water on Anatase (101). J. Phys. Chem. C 2010, 114, 1278–1284. (10) K¨oppen, S.; Bronkalla, O.; Langel, W. Adsorption Configurations and Energies of Amino Acids on Anatase and Rutile Surfaces. J. Phys. Chem. C 2008, 112, 13600– 13606. (11) Schneider, J.; Ciacchi, L. C. Specific Material Recognition by Small Peptides Mediated by the Interfacial Solvent Structure. J. Am. Chem. Soc. 2012, 134, 2407–13. (12) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864– B871.

17

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 49 50 51 52 53 54 55 56 57 58 59 60

(13) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (14) Lazzeri, M.; Vittadini, A.; Selloni, A. Structure and Energetics of Stoichiometric TiO2 Anatase Surfaces. Phys. Rev. B 2001, 63, 155409. (15) Vittadini, A.; Selloni, A.; Rotzinger, F. P.; Gr¨atzel, M. Structure and Energetics of Water Adsorbed at TiO2 Anatase (101) and (001) Surfaces. Phys. Rev. Lett. 1998, 81, 2954–2957. (16) Bandura, A. V.; Kubicki, J. D. Derivation of Force Field Parameters for TiO2 -H2 O Systems from Ab Initio Calculations. The Journal of Physical Chemistry B 2003, 107, 11072–11081. (17) Perdew, J.; Burke, K.; M.Ernzerhof, Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77(4), 3865–68. (18) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I. et al. QUANTUM ESPRESSO: a Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502–1–19. (19) Bengtsson, L. Dipole Correction For Surface Supercell Calculations. Phys. Rev. B 1999, 59, 12301. (20) Monkhorst, H.; Pack, J. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1973, 13(5), 5188–92. (21) Fletcher, R. A New Approach to Variable Metric Algorithms. The Computer Journal 1970, 13(6), 317–22. (22) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a LongRange Dispersion Correction. J. Comput. Chem. 2006, 27(15), 1787–99. 18

ACS Paragon Plus Environment

Page 18 of 20

Page 19 of 20

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

(23) Henkelman, G.; Uberuaga, B. P.; J´onsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Point and Minimum Energy Paths. J. Chem. Phys. 2000, 113(22), 9901–9904. (24) Sheppard, D.; Terrel, R.; Henkelman, G. Optimization Methods for Finding Minimum Energy Paths. J. Chem. Phys. 2008, 128, 134106–1/134106–10. (25) Baroni, S.; de Gironcoli, S.; Corso, A. D.; Giannozzi, P. Phonons and Related Crystal Properties from Density-Functional Perturbation Theory. Rev. Mod. Phys. 2001, 73, 515–562. (26) Zhao, Z.; Li, Z.; Zou, Z. Understanding the Interaction of Water with Anatase TiO2 (101) Surface from Density Functional Theory Calculations. Phys. Lett. A 2011, 375, 2939–2945. (27) Agosta, L.; Gala, F.; Zollo, G. Water Diffusion on TiO2 Anatase Surface. AIP Conf. Proc. 2015, 1667, 020006–1–020006–9. (28) Ponti, G.; Palombi, F.; Abate, D.; Ambrosino, F.; Aprea, G.; Bastianelli, T.; Beone, F.; Bertini, R.; Bracco, G.; Caporicci, M. et al. The Role of Medium Size Facilities in the HPC Ecosystem: the Case of the New CRESCO4 Cluster Integrated in the ENEAGRID Infrastructure. IEEE HPCS 2014, 6903807, 1030–1033.

19

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 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6: Graphical table of contents

20

ACS Paragon Plus Environment

Page 20 of 20