7678
J. Phys. Chem. C 2008, 112, 7678–7688
A Shell Model for Atomistic Simulation of Charge Transfer in Titania Sebastien Kerisit,* N. Aaron Deskins, Kevin M. Rosso, and Michel Dupuis Chemical and Materials Sciences DiVision, Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: January 26, 2008; ReVised Manuscript ReceiVed: February 23, 2008
The derivation of atomistic potential parameters, based on electronic structure calculations, for modeling electron and hole polarons in titania polymorphs is presented. The potential model is a polarizable version of the Matsui and Akaogi model (Matsui, M.; Akaogi, M. Mol. Simul. 1991, 6, 239) that makes use of a shell model to account for the polarizability of oxygen ions. The -1 and +1 formal charges of the electron and hole polarons, respectively, are modeled by delocalizing the polaron’s charge over a titanium or oxygen ion, respectively, and its first nearest-neighbors. The charge distributions and the oxygen polarizability were fitted to the reorganization energies of a series of electron and hole polaron transfers in rutile and anatase obtained from electronic structure calculations at zero Kelvin and validated against lattice deformation due to both types of polaron. Good agreement was achieved for both properties. In addition, the potential model yields an accurate representation of a range of bulk properties of several TiO2 polymorphs as well as Ti2O3. The model thus derived enables us to consider systems large enough to investigate how the charge transfer properties at titania surfaces and interfaces differ from those in the bulk. For example, reorganization energies and free energies of charge transfer were computed as a function of depth below vacuum-terminated rutile (110) and anatase (001) surfaces using a mapping approach first introduced by Warshel (Warshel, A. J. Phys. Chem. 1982, 86, 2218). These calculations indicate that deviations from bulk values at the surface are substantial but limited to the first couple of surface atomic layers and that polarons are generally repelled from the surface. Moreover, attractive subsurface sites may be found as is predicted for hole polarons at the rutile (110) surface. Finally, several charge transfers from under-coordinated surface sites were found to be in the so-called Marcus inverted-region. Introduction Titania (TiO2) has been the subject of extensive research due in part to its photocatalytic properties1 with potential applications such as the degradation of organic pollutants2 and water splitting for hydrogen production.3 For example, some fundamental mechanisms of photocatalysis have been investigated in a series of papers on the photoinduced oxidation of trimethylacetic acid (TMA) adsorbed on the rutile (110) surface.4–9 In these experiments, TMA is initially deposited as a monolayer on the surface and subsequently, under UV irradiation, photons are absorbed in the crystal to create electron-hole pairs. The charges migrate to the crystal surface where holes react with adsorbed TMA molecules, which ultimately decompose into isobutene and CO2. Alternatively, electrons and holes may also recombine within the crystal or at its surface and thus lower the quantum yield of the photochemical process. Therefore, a better understanding of the transport properties of electrons and holes in these materials, including at surfaces and interfaces, is crucial to the improvement of photoactive materials. The transport of electrons and holes in titania has been described by several researchers using a small polaron model.10–14 Polarons can be considered as quasiparticles that consist of an electron or hole, self-trapped at a lattice site, and its accompanying lattice polarization. A polaron is termed small or large depending on the strength of the electron-phonon coupling. In titania, an electron small polaron is localized on a titanium site thereby reducing it to Ti3+, whereas a hole small polaron is * To whom correspondence should be addressed. E-mail: sebastien.kerisit@ pnl.gov.
centered on a oxygen site thus formally forming a O- species. Small polarons have a high effective mass and the spatial extent of their accompanying distortion is smaller than a lattice constant. Small polaron motion in crystal lattices can follow two mechanisms:15–17 a band process, which occurs at low temperature, and an activated hopping process, which dominates at high temperature. The temperature above which the hopping mechanism dominates is often estimated as being approximately half-the Debye temperature.16 The Debye temperature of TiO2 was derived experimentally to be 530 K,18 suggesting that polaron motion occurs via a hopping mechanism in TiO2 at room temperature. Recently, we computed the rate of elementary electron polaron hops in stoichiometric rutile and anatase19 using electronic structure calculations and a small polaron hopping model based on the seminal work of Holstein,15,16,20,21 Austin and Mott,22 and Marcus.23,24 A similar approach was used to model electron and hole polaron transfers in chromia (RCr2O3),25 hematite (R-Fe2O3),26,27 and wu¨stite (FeO).28 Good quantitative agreement with experiment was obtained with this approach. In particular, the calculated parameters for electron transfers in the basal plane and c direction of hematite successfully reproduced the known anisotropy of the electrical conductivity.26,27 In addition, the temperature dependence of the electrical conductivity in hematite obtained from a kinetic Monte Carlo model based on the calculated elementary rates was also well reproduced.29 Our recent density functional theory (DFT) study19 of titania was focused on the intrinsic charge transport rates. For each polymorph two directions for electron transfer were considered. The diabatic energy, calculated using the
10.1021/jp8007865 CCC: $40.75 2008 American Chemical Society Published on Web 04/25/2008
Atomistic Simulation of Charge Transfer in Titania
J. Phys. Chem. C, Vol. 112, No. 20, 2008 7679
Figure 1. Free energy diagram of a charge transfer reaction as a function of the energy gap.35–37
TABLE 1: Core and Shell Charges Used in This Work, where (l) Indicates Lattice Ions and (e) and (h) Indicate the Ions That Belong to an Electron or Hole Polaron, Respectivelya species
core (e)
species
core (e)
shell (e)
Ti(l) Ti(e) Ti(e)u Ti(h) Ti(h)u
2.196000 1.647000 1.571833 2.374000 2.412000
O(l) O(e) O(e)u O(h) O(h)u
0.500 0.500 0.500 0.500 0.500
-1.598000 -1.673167 -1.673167 -1.132000 -1.030000
a The subscript “u” indicates under-coordinated surface sites. All oxygen core-shell spring constants, k, are set to 44.3 eV.Å-2. 2 Core-shell interaction potential form: V ) k × rc-s .
TABLE 2: Buckingham Potential Parameters Used in This Worka ion pair (ij)
Aij (eV)
Fij (Å)
Cij (eV. Å6)
Ti(l,e,h)-Ti(l,e,h) Ti(l,e,h)-O(l,e,h) O(l,e,h)-O(l,e,h)
31120.53 16957.71 11782.88
0.154 0.194 0.234
5.25 12.59 30.22
a
These potential parameters also apply to under-coordinated surface sites. Buckingham potential form: Vij ) Aij × exp(-rij/Fij) Cijrij-6.
DFT+U approach,30 was found to be approximately 0.3 eV for all four transfers. The electronic coupling matrix element, however, was found to be much larger for the [001] transfer in rutile than for the other three transfers due to a greater overlap of the Ti d orbitals in that direction. Using this work as a baseline, our objective is to extend the computational model of charge transport in titania to include the effect of defects, surfaces, and interfaces. Investigating surfaces and interfaces will require us to consider much larger systems that may contain up to several thousands of particles and thus go beyond what is currently achievable with electronic structure calculations. Therefore, in this paper, we calibrate a potential model for simulating electron and hole polarons in titania with the DFT data. It was shown previously for hematite31 and wu¨stite28 that potential-based models are able to reproduce ab initio reorganization energies. The model reported here is an improved version of that of Matsui and Akaogi (MA).32 In previous work,
Figure 2. Electron (blue lines) and hole (green lines) polaron transfer directions in rutile (top) and anatase (bottom). Electron transfer directions are labeled e and c for the edge and corner transfer, respectively. Hole polaron directions are labeled A- E in rutile and A- F in anatase. Only nearest-neighbor hops are considered. Titanium ions are shown in white and oxygen ions are shown in red.
TABLE 3: Comparison of the Experimental and Calculated Bulk Properties of Ti2O3 property lattice parameter - a (Å) lattice parameter - c (Å) volume - V (cm3.mol-1) Ti-O distance (Å) bulk modulus (GPa) elastic constant - C11 elastic constant - C12 elastic constant - C13 elastic constant - C14 elastic constant - C33 elastic constant - C44 a
(GPa) (GPa) (GPa) (GPa) (GPa) (GPa)
expt.
calc.
5.157a 13.610a 15.73a 2.067 (3)a 2.026 (3)a 207b 309.7c 131.6c 162.3c -2.4c 296.9c 106.9c
5.111 13.642 15.49 2.074 (3) 1.989 (3) 188 298.9 147.7 139.2 -3.0 249.5 99.1
Reference 60. b Reference 61. c Reference 62.
we used the original MA model to characterize a series of rutileanatase interfaces33 constructed using the near-coincidence-site lattice theory.34 However, the MA model is a rigid-ion model and such models have been shown to overestimate reorganization energies to a large extent.28,31 Therefore, in the first part of this work, we introduce a polarizable version of the MA model that correctly reproduces the DFT parameters for electron and
7680 J. Phys. Chem. C, Vol. 112, No. 20, 2008
Kerisit et al.
TABLE 4: Reorganization Energies (λ) of a Series of Electron and Hole Polaron Transfers in Rutile and Anatase as Obtained from the Rigid-Ion and Shell Models and DFT Calculationsa λ - rigid-ion
λ - shell model
λ - DFT
polymorph
direction
rutile rutile anatase anatase
e c e c
electron polaron 2.10 2.35 2.12 2.41
1.12 1.21 1.15 1.24
1.15 1.23 1.19 1.21
rutile rutile rutile rutile rutile anatase anatase anatase anatase anatase anatase
A B C D E A B C D E F
hole polaron 2.44 2.13 2.33 3.14 3.18 2.63 2.75 3.15 2.13 3.30 2.51
1.89 1.80 1.84 2.36 2.41 1.74 2.19 2.30 1.79 2.51 1.97
2.21 2.16 2.00 2.35 2.49 1.96 2.04 2.09 1.68 2.35 2.04
a
a
All reorganization energies are in eV.
TABLE 5: Comparison of the Changes in Titanium-Oxygen Distances (∆r) Obtained with the Shell Model and DFT Calculations upon Formation of Electron and Hole Polarons in Rutile and Anatasea polymorph distance change (∆r) (Å) ∆Ti-O (1st shell) (×2) ∆Ti-O (1st shell) (×2) ∆Ti-O (1st shell) (×2) ∆O-Ti ∆O-Ti ∆O-Ti ∆O-Ti ∆O-Ti ∆O-Ti ∆O-Ti a
rutile ∆r - shell model
anatase ∆r DFT
∆r - shell model
∆r DFT
electron polaron 0.042 0.037 0.087 0.095 0.087 0.095
0.146 0.047 0.047
0.154 0.068 0.035
hole polaron 0.157 0.194 0.157 0.199 0.136 0.156 0.026 0.015 0.025 0.022 -0.006 -0.145 -0.007 0.036
0.135 0.135 0.282 -0.103 -0.104 -0.103 -0.103
0.120 0.120 0.290 -0.114 -0.114 -0.112 -0.112
(1st shell) (1st shell) (1st shell) (2nd shell) (2nd shell) (2nd shell) (2nd shell)
All distance changes are in Å.
TABLE 6: Energy Difference (∆Ef) between the Polaron Formation Energies in Rutile (EfR) and Anatase (EfA) Obtained with the Shell Model and DFT Calculations (∆Ef ) EfR - EfA)a
a
TABLE 7: Electron and Hole Polaron Binding Energies in Rutile and Anatasea
polaron
shell model
DFT
electron polaron hole polaron
0.273 -0.436
0.709 -1.300
All energies are in eV.
hole transfers in both rutile and anatase as well as a series of bulk properties of several TiO2 polymorphs and Ti2O3. In addition, in the second part of the paper, we make use of this model to determine the influence of a free surface on the energetics of polaron transfer, using the rutile (110) and anatase (001) surfaces as examples. Computational Methods Charge Transfer Model. In this section, we give a brief description of the charge transfer model used in this work. For a complete discussion of the charge transfer model, see refs 26
polaron
rutile
anatase
electron polaron hole polaron
1.02 1.49
1.01 1.53
All energies are in eV.
and 27 and references therein. In this model, charge transfer is modeled as the activated thermal hopping of electrons and holes localized as small polarons. In the polaron model, electrons and holes introduced in a crystal lattice interact strongly with phonons and lower their energy by localizing at a lattice site thereby creating a local distortion of the lattice. A small polaron refers to the unique case where the distortion is localized to a region smaller than a lattice constant. Polaron transfer is described by a two-state model whereby the reactants (before charge transfer) and products (after charge transfer) potentialenergy surfaces are assumed to be parabolic with respect to the reaction coordinate, which is a function of the nuclear coordinates. Charge transfer occurs when thermal fluctuations cause the reactants and products states to be energetically degenerate. This model provides a basis for several theories, such as those of Holstein,15,16,20,21 Austin and Mott,22 and Marcus,23,24 which describe the transfer of polarons in various systems. In what follows, for simplicity, we use the terminology established by Marcus. The charge transfer between two sites is described by the following charge transfer parameters: the reorganization energy, the electronic coupling matrix element, the diabatic free energy of activation, and the free energy of reaction. The reorganization energy, λ, is the energy to distort the configuration of the reactants into that of the products, and vice versa, without changing the electronic distribution; the electronic coupling matrix element, VAB, is the amount of electronic interaction between the reactants and products states at the transition state; the diabatic free energy of activation, ∆G*′, is the energy required to thermally excite the system to the transition state configuration; and the free energy of reaction, ∆G°, is the change in free energy upon charge transfer. Warshel and co-workers35–37 introduced an approach for evaluating, using computer simulations, the potential free energy surfaces as a function of the reaction coordinate, which is defined as the energy gap, ∆E, that is, the energy difference between the reactants and products charge distributions for a particular configuration of the system. The electron transfer parameters are illustrated in Figure 1. As the potential energy surfaces are parabolic, ∆G*′ is directly related to λ and ∆G° and therefore the number of independent parameters needed to describe a charge transfer reaction is reduced to three. Charge transfer can take place in two different regimes: adiabatic and nonadiabatic, depending mainly on the magnitude of VAB. Different rate expressions are used to calculate the rate of charge transfer depending on whether the charge transfer is adiabatic or nonadiabatic. The primary goal of this work is the calibration of a potential model to a DFT data set; therefore, we do not calculate VAB, which cannot be obtained from potential-based calculations, for the charge transfers of interest and instead focus on the calculation of λ and ∆G°. The DFT data set consists of symmetrical polaron transfers in bulk rutile and anatase. Therefore, ∆G° is zero for all these polaron transfers. In addition, the polaron transfers were modeled by carrying out calculations along the reaction coordinate, which was determined from linear interpolation of the energy mini-
Atomistic Simulation of Charge Transfer in Titania
J. Phys. Chem. C, Vol. 112, No. 20, 2008 7681
TABLE 8: Calculated and Experimental Bulk Properties of a Series of Titania Polymorphs polymorph property a (Å) b (Å) c (Å) V (cm3.mol-1) Ti-O (Å) a
rutile
anatase
expt. a
4.594 4.594a 2.959a 18.80 1.980 (2)a 1.949 (4)a
calc.
expt.
calc.
3.785b
4.506 4.506 3.007 18.39 1.925 (2) 1.955 (4)
brookite
3.788 3.788 9.483 20.49 1.995 (2) 1.925 (4)
3.785b 9.512b 20.51 1.980 (2)b 1.934 (4)b
9.174c 5.449c 5.138c 19.33 1.863-2.052
calc.
expt. 4.515d
9.153 5.408 5.140 19.16 1.916-1.999
5.497d 4.939d 18.45 2.050 (2)d 1.919 (4)d
calc. 4.524 5.395 4.939 18.15 1.986 (2) 1.934 (4)
Reference 63. b Reference 64. c Reference 65. d Reference 66.
TABLE 9: Comparison of the Calculated and Experimental Elastic Moduli and Elastic Constants of Rutile and Anatase polymorph property
rutile expt.
anatase calc.
expt./ ab inito
bulk modulus shear modulus
elastic moduli (GPa) 236 212a 113a 123
C11 C33 C44 C66 C23 C12
elastic constants (GPa) 268.0a 316.2 320c 484.2a 438.3 190c 123.8a 117.5 54c 190.2a 222.1 60c 147.4a 151.7 143c 174.9a 226.6 151c
a
expt.
TiO2 II
178b 58c
calc. 173 61 442.3 199.1 27.1 49.5 102.2 81.4
Reference 67. b Reference 68. c Reference 69.
[
P(X) P(〈∆E〉)
]
Vθ ) (1 - θ)VA + θVB
(1)
where 〈∆E〉 is the ensemble average of the energy gap. In theory, each free energy curve could be obtained from a single molecular dynamics simulation. The diabatic activation free energy would be the value of ∆G for ∆E ) 0. However, as it is often the case, the activation free energy is too high for the region of phase space in the vicinity ∆E ) 0 to be significantly sampled with a single molecular dynamics calculation. To overcome this problem, we employed the umbrella sampling technique as first introduced, for an electron transfer reaction, in seminal papers by Hwang and Warshel,36 Kuharski et al.,39 and King and Warshel.37 In this method, the atomic charges of the reactants are gradually changed into those of the products
(2)
θ is varied from 0 to 1 in a series of discrete values and a simulation is carried out for each value. This biases the system to generate configurations that have values of ∆E significantly different from 〈∆E〉. Therefore, a thorough sampling of the phase space along the reaction coordinate is achieved. The bias thus introduced is removed, once the simulation has been run, by adjusting the free energy as follows (see Kuharski et al.39 for detail of the derivation of eq 3):
[
∆G(X) ) -RT ln
mized initial and final polaron configurations, that is, temperature was not taken into account.19 Consequently, in the first part of the paper, which is concerned with the derivation of a potential model based on the DFT data set, the reorganization energies were calculated using the so-called “direct” method38 at zero Kelvin to enable a direct comparison with the DFT data. In the direct method, the reorganization energy is calculated as the difference between the energy of the products state (∆GP curve in Figure 1) in the reactants (A) and products (B) configurations. Alternatively, the reorganization energy can be computed from the difference between the energy of the reactants state (∆GR curve) in the products and reactants configurations. In the second part of this work, we focus on charge transfer reactions at titania surfaces at room temperature. Therefore, to determine λ and ∆G°, we combine molecular dynamics simulations and the umbrella sampling technique and calculate the free energy surfaces of the reactants and products states as a function of the energy gap ∆E. Following Warshel’s work,35 the free energy curves are calculated from the probability of finding the system in a state for which ∆E ) X versus the probability of being in a state for which ∆E ) 〈∆E〉:
∆G(X) ) -RT ln
and vice-versa. The new potential of the system is calculated from a linear combination of the reactants and products state potentials (VA and VB respectively):
P(X)θ P(〈∆E〉)
]
- θX + C
(3)
where C is a normalization constant and P(X)θ is the probability of finding the system at ∆E ) X when its potential is Vθ. Finally, charge transfer is a Franck-Condon process in that the electronic transfer is instantaneous relative to the nuclei motion and hence the system polarizability needs to be calculated for both the reactants and products states. As our model accounts for the polarizability of anions by means of a shell model, the position of the shells must be relaxed for each reactants and products state configuration used in eq 3. The energy minimization of the shells is carried out using the steepest descent method. Potential Model. The calculations reported in this work are based on the Born model of solids,40 in which the atoms of a system are represented as point-charge particles that interact via long-range Coulombic forces and short-range interactions. The latter are described by parametrized functions and represent the repulsion between electron-charge clouds and the van der Waals attraction forces. In this work, the short-range interactions are described using a Buckingham potential and therefore the pairwise interaction energy takes the following form:
Uij )
( )
rij Cij 1 qiqj + Aij exp 4πε0 rij Fij rij
(4)
The parameters (q, A, F, and C) used in this study are those optimized by Matsui and Akaogi32 (MA) for the TiO2 polymorphs rutile, anatase, brookite, and TiO2-II (see Tables 1 and 2). It should be noted that this model uses partial atomic charges rather than formal charges. This potential model has proved to be an exceptionally reliable and transferable model. For example, Swamy et al.41 showed that the MA potential was able to reproduce the crystal structures of a range of TiO2 polymorphs, many of which were not included in the original derivation set and also had titanium coordination numbers ranging from 6 to 9 whereas the polymorphs used for the derivation were only examples of 6-fold coordinated titanium. In addition to structural aspects, the MA potential provides an accurate description of the relative stability of TiO2 polymorphs
7682 J. Phys. Chem. C, Vol. 112, No. 20, 2008
Kerisit et al.
Figure 3. Electron (left) and hole (right) polaron transfer directions at the rutile (110) (top) and anatase (001) (bottom) surfaces. The transfer direction labels are the same as used in Figure 2, with additional numbering when more than one transfer of a particular type is possible due to the lowering of the symmetry by the surface. In all cases, only the first two surface layers are shown.
TABLE 10: Ti-O Bond Lengths and Reorganization Energies for the e Electron and A Hole Polaron Transfers in Bulk Rutile and e2 Electron and A1 Hole Polaron Transfers at the Rutile (110) Surface electron polaron
a
hole polaron
model
Ti-Oax (Å)
Ti-Oeq (Å)
Ti-Oeq (Å)
λ (eV)
Ti-O (Å)
λ (eV)
DFT - bulk w/o polaron DFT - surface w/o polarona DFT - surface w/ polarona shell model - bulk w/o polaron shell model - surface w/o polarona shell model - surface w/ polarona
2.01
1.96
1.96
1.15
1.96
2.20
1.93
1.96
1.96
-
1.82
-
2.13
2.03
2.04
1.59
2.04
2.48
1.92
1.96
1.96
1.12
1.96
1.89
1.89
1.95
1.95
-
1.90
-
1.95
2.02
2.02
1.23
2.14
2.12
Five-fold coordinated surface Ti site for electron polaron column and topmost (i.e., surface bridging) oxygen site for hole polaron column.
and, as shown by Dubrovinskaia et al.,42 was able to predict the existence of an intermediate phase at high pressure, which was confirmed experimentally. Finally, Bandura and Kubicki43 showed that the MA potential gave good agreement with the surface bond lengths obtained from periodic density functional theory calculations at the rutile (110) surface thus indicating that, although originally derived for bulk systems, the model is transferable to surfaces. However, the MA potential is a rigidion model and we have shown in the past28 that such models, in which there is no explicit account of the ions’ polarizability, are not able to reproduce the reorganization energy associated with polaron transfers in bulk oxides. Therefore, we derived a shell model version of the MA potential that improves the agreement with reorganization energies obtained from DFT calculations of several polaron transfers in titania polymorphs. In the shell model,44 a polarizable ion is composed of two
particles, a core and a shell, which share the ion’s charge and are linked by a harmonic spring, k:
Uc-s ) k × r2c-s
(5)
where rc-s is the core-shell separation distance. The free-ion polarizability, R, is given by the expression:
R)
Y2 k
(6)
where Y is the shell charge. In addition, we also introduce in this paper potential parameters for modeling electron and hole small polarons. A full description of the approach taken to derive these potential parameters is given in the Results and Discussion section. All the potential parameters used in this work are summarized in Tables 1 and 2.
Atomistic Simulation of Charge Transfer in Titania
J. Phys. Chem. C, Vol. 112, No. 20, 2008 7683 ture), whereas the surface simulations were carried out in the NVT ensemble (constant number of particles, constant volume, and constant temperature). In these simulations, the temperature and pressure were kept constant by use of the Nose´-Hoover thermostat46 and the Hoover barostat,47 respectively. The electrostatic forces were calculated by means of the Ewald summation method.48 A 9 Å cutoff was used for the short-range interactions and the real part of the Ewald sum. The Verlet leapfrog integration algorithm was used to integrate the equations of motion with a time step of 0.2 fs. The shells were given a mass of 0.2 au, and their motion was treated as that of the cores following the adiabatic shell model first introduced by Mitchell and Fincham.49 Results and Discussion
Figure 4. Examples of free energy profiles obtained for the c2 electron transfer (“layer 1 to 2”) and equivalent transfers in layers below (“layer 2 to 3” and “layer 3 to 4”) (see Figure 3). In all three profiles, the free energy surface on the left-hand side is that for the case where the electron polaron is in the upper layer (i.e., the direction indicated in the key is from the parabola on the left to the parabola on the right).
TABLE 11: Reorganization Energies (λ) Obtained for Electron Transfers at the Rutile (110) Surface from Umbrella Sampling Calculations at 300 Ka electron transfer surface layer
λ - e1
λ - e2
λ - c1
λ - c2
λ - c3
1 2 3 4 Bulk
1.02 0.95 0.93 0.90 0.93
1.01 0.94 0.89 0.94 0.93
1.11 1.03 1.01 1.07 1.03
1.21 1.09 1.05 1.03
1.10 1.06 1.03 1.03
a All reorganization energies are in eV. Transfer directions (e1, e2, c1, c2, and c3) are illustrated in Figure 3. Reorganization energies in bold indicate polaron transfers in the inverted region. Reorganization energies of equivalent transfers in the bulk are given as a comparison.
TABLE 12: Reorganization Energies (λ) Obtained for Hole Polaron Transfers at the Rutile (110) Surfacea hole transfer surface layer λ - A1 λ - A2 λ - B1 λ - B2 λ - C1 λ - C2 1 1 2 2 3 3 4 4 Bulk
2.08 1.74 1.96 1.71 1.71 1.71 1.69 1.69 1.66
1.89 1.72 1.67 1.64 1.66
1.91 1.77 1.63 1.58 1.64
1.86 1.57 1.53 1.63 1.64
1.91 1.82 1.62 1.66 1.62 1.64 1.66
1.96 1.64 1.63 1.61 1.63 1.68 1.66
a All reorganization energies are in eV. Transfer directions (A1, A2, B1, B2, C1, and C2) are illustrated in Figure 3. Reorganization energies in bold indicate polaron transfers in the inverted region. Reorganization energies of equivalent transfers in the bulk are given as a comparison.
Atomistic Simulations. Unless otherwise noted, all the atomistic simulations in this work where carried out with the computer program DL_POLY.45 Calculations that followed the “direct” method were performed at zero Kelvin with the steepest descent approach. All the umbrella sampling calculations were run at 300 K and zero-applied pressure. The bulk simulations were performed in the NPT ensemble (constant number of particles, constant pressure, and constant tempera-
Potential Model Derivation. We first set out to derive parameters for modeling electron polarons in rutile and anatase. For each polymorph, the two electron transfers to nearestneighbor sites studied previously19 were considered and are illustrated in Figure 2. The first transfer takes place along the [001] direction in rutile and the [201] direction in anatase, that is, between two edge-sharing TiO6 octahedra in both polymorphs. The second transfer occurs between two corner-sharing TiO6 octahedra along the [111] direction in rutile and the [100] direction in anatase. The two transfers will be referred to as the edge and corner transfers hereafter, respectively. Electron polarons localize on titanium sites thereby formally reducing Ti4+ to Ti3+. In Ti2O3, all Ti atoms are formally Ti3+ ions. Therefore, given the oxygen partial charge in the MA model (qO(l)), the Ti3+ partial charge (qTi(e)) was determined from the Ti2O3 stoichiometry, that is, qTi(e) ) -3 × qO(l)/2. We evaluated whether the derived partial charge was adequate to model Ti3+ ions by comparing a range of calculated bulk properties of Ti2O3 with experimental data, as shown in Table 3. The agreement is excellent with lattice parameters within 1% and elastic constants in close agreement, which is a testament to the quality of the model derived by Matsui and Akaogi. Because the MA potential employs partial charges, the charge difference for reducing Ti4+ to Ti3+ is less than -1.0 e. Consequently, the remaining charge difference required to introduce an electron into the lattice was spread equally over the six oxygen atoms that constitute the first coordination shell. As mentioned previously, rigid-ion models have been found to overestimate the reorganization energy as they lack the stabilization due to the lattice polarizability. This was found to also be the case for all four of the electron transfers considered in this work in rutile and anatase when compared with the DFT data set,19 as shown in Table 4. Therefore, the oxygen ions were modeled by a shell model and the oxygen polarizability was determined by minimizing the root mean-square difference between the reorganization energies calculated with the shell model and those obtained from the DFT calculations. Both sets of reorganization energies were obtained using energy minimization calculations as described in the Charge Transfer Model subsection of the Computational Methods section. The DFT calculations made use of 3 × 3 × 3 rutile and 3 × 3 × 2 anatase supercells, whereas the rigidion and shell model simulations employed 5 × 5 × 6 rutile and 6 × 6 × 2 anatase supercells. The oxygen core charge was selected arbitrarily and the shell charge was obtained by subtracting the core charge to the total oxygen charge (see O(e) entry in Table 1). Finally, the spring constant was determined from the shell charge and the optimized value of the polarizability using eq 6 (Table 1). The agreement with the values obtained from the DFT calculations is greatly improved upon
7684 J. Phys. Chem. C, Vol. 112, No. 20, 2008
Kerisit et al.
Figure 5. Electron (left) and hole (right) polaron free energy profiles at the rutile (110) surface. Lines are shown as a guide to the eye. A section of the unrelaxed rutile (110) surface (middle) is shown for reference.
TABLE 13: Reorganization Energies (λ) Obtained for Electron Transfers at the Anatase (001) Surface from Umbrella Sampling Calculations at 300 Ka electron transfer surface layer
λ-E
λ - C1
λ - C2
1 2 3 4 5 6 Bulk
1.15 1.03 1.04 1.04 1.04 1.00
1.68 1.12 1.10 1.09 1.05 1.12 1.08
1.88 1.09 1.11 1.08 1.06 1.09 1.08
a All reorganization energies are in eV. Transfer directions (e, c1, and c2) are illustrated in Figure 3. Reorganization energies in bold indicate transfers in the inverted region. Reorganization energies of equivalent transfers in the bulk are given as a comparison.
TABLE 14: Reorganization Energies (λ) Obtained for Hole Transfers at the Anatase (001) Surface from Umbrella Sampling Calculations at 300 Ka Figure 6. Electrostatic potential at oxygen sites at the rutile (110) surface before and after energy minimization. Both electrostatic potential profiles are with respect to the electrostatic potential of an oxygen site in the middle of the rutile slab. Lines are shown as a guide to the eye.
introducing the shells. The reorganization energy of each of the electron transfer decreases by approximately 50%; however, there is no change in the relative order of the reorganization energies in going from the rigid-ion to the shell model. In both cases, the order closely reflects that obtained from the electronic structure calculations. Next, we derived potential parameters for modeling hole polarons using the same simulation cells employed above. Hole polarons are modeled as small polarons localized on an oxygen atom and they therefore diffuse on the oxygen sublattice. Five directions for nearest-neighbor hole transfer were considered in rutile and six in anatase. All the transfer directions are illustrated in Figure 2. As for the electron transfers, the reference reorganization energies were obtained from DFT calculations and will be described in detail in an upcoming publication.50 Following the approach employed for the electron polaron model, the polaron charge is distributed over an oxygen atom
hole transfer surface layer
λ-A
λ-B
λ-D
λ-F
1 2 3 4 5 Bulk
1.86 1.93 1.93 1.96 1.91
2.43/2.24 2.10 2.03 2.06 2.11 2.03
1.80 1.64 1.66 1.12 1.64 1.57
2.05 1.84 1.72 1.81 1.76 1.78
a All reorganization energies are in eV. Transfer directions (A, B, D, and F) are illustrated in Figure 3. Reorganization energies in bold indicate transfers in the inverted region. Reorganization energies of equivalent transfers in the bulk are given as a comparison.
and its three titanium nearest-neighbors. The shell spring constant of the oxygen on which the hole polaron is localized was set to the value determined above and the distribution of the charge between the oxygen atom and its titanium neighbors was fitted to minimize the root-mean square difference between the reorganization energies obtained from the shell model and the DFT calculations. In introducing the shells, the decrease in reorganization energy ranges from 15 to 35%. Good overall
Atomistic Simulation of Charge Transfer in Titania
J. Phys. Chem. C, Vol. 112, No. 20, 2008 7685
Figure 7. Electron (left) and hole (left) free energy profiles at the anatase (001) surface. Lines are shown as a guide to the eye. A section of the unrelaxed anatase (001) surface (middle) is shown for reference.
agreement was obtained with the DFT calculations with respect to the relative ordering of the reorganization energies. To validate the polaron models, we compared the lattice distortion due to the electron and hole polarons as obtained from the shell model and the electronic structure calculations. As shown in Table 5, even though this property was not included in the fit, the agreement with the DFT calculations is satisfactory. For the electron polarons, which are localized on Ti sites, the difference in the Ti-O bond length change is at most 0.02 Å for both polymorphs. There is one noticeable discrepancy in anatase whereby the shell model predicted a symmetric increase of two sets of Ti-O bonds (∆Ti-O ) 0.047 Å for both) whereas the electronic structure calculations gave an asymmetric lengthening of the same two sets of bonds (∆Ti-O ) 0.068 and 0.035 Å). However, the average increase of the two sets as obtained from the two methods is in excellent agreement. The agreement is also close for the hole polarons, which are localized on oxygen sites, with differences in O-Ti bond length distortions of at most 0.04 Å. The large discrepancy between the two methods for two of the distances to second-shell titanium ions in rutile is likely due to the small size of the rutile cell used in the DFT calculations, which causes the two distances to be to the same titanium ion. We carried out an additional validation test which consisted in comparing the energetics of polaron formation in the shell model with those obtained from the DFT calculations. Here, we define the polaron formation energy as the energy difference between the energy minimized titania supercell without a polaron and the same supercell after energy minimization in the presence of a polaron. In the latter simulation the net charge was neutralized by a uniform background charge.51 In addition, the energy was corrected for the polaron-polaron Coulombic interactions due to the periodic boundary conditions using the approach described by Watson et al.52 and implemented in the computer code PARAPOCS.53 The polaron formation energy was calculated for both rutile (ERf ) and anatase (EAf ), and the difference between the rutile and anatase polaron formation energies (∆Ef ) ERf - EAf ) was compared with that obtained from the DFT calculations for both the electron and hole polarons (Table 6). In both cases, the energy calculated by the shell model is approximately 35-40% that of the DFT value.
However, the sign of the energy difference is the same indicating that the shell model reproduces correctly the relative stability of the polarons in the two polymorphs. These calculations suggest that electron polarons have a greater intrinsic affinity for the anatase phase whereas hole polarons are energetically favored in the rutile phase. It should be noted that the affinities calculated in this work differ from those predicted from the relative position of the rutile and anatase conduction bands54 because small polarons are electrons and holes self-trapped at lattice sites whose energy lie within the band gap. Using the newly developed potential model, we also computed the polaron binding energy which is defined as the difference between the energies obtained from an energy minimization of the titania supercell with a polaron whereby only the shells are relaxed and one of the same supercell with all particles allowed to relax. As shown in Table 7, for both the electron and hole polaron models, the polaron binding energies in anatase and rutile are essentially identical. Finally, we verified that the introduction of the shells on the oxygen ions did not diminish the quality of the original model with respect to the bulk properties of several titania polymorphs. Table 8 shows that the lattice parameters remain in excellent agreement with the experimental values for all polymorphs. In addition, for rutile and anatase, for which there is experimental and ab initio data of the elastic moduli and elastic constants, Table 9 shows that the values obtained with the shell model are in good accord. Charge Transfer at Titania Surfaces. In the second part of this work, we investigated electron and hole transfers at vacuum-terminated rutile (110) and anatase (001) surfaces at room temperature. Both surfaces were studied in their (1 × 1) configuration. The rutile (110) surface is the most stable and most studied rutile surface.55 The anatase (001) surface is not the most thermodynamically stable anatase surface but has received some attention55,56 and has a relatively high symmetry, which allows for a comprehensive investigation with a small number of polaron transfers. Rutile (110) Surface. To create the rutile slab, a constant pressure geometry optimization was performed on the rutile unit cell which was then oriented with its z axis along the [110] direction. The slab was grown to 8 unit-cell thick and the surface
7686 J. Phys. Chem. C, Vol. 112, No. 20, 2008 vectors were scaled 6 × 3 to generate a slab with a surface area of 18.01 × 19.08 Å2 and a thickness of approximately 25 Å. A 25 Å vacuum gap was introduced between the two faces of the rutile slab. All possible nearest-neighbor electron transfers in one-half of the slab were considered. There are two symmetrically nonequivalent Ti sites at the rutile (110) surface: the 6-fold and 5-fold coordinated sites (see Figure 3). Due to the symmetry of the surface there are two possible edge transfers (Figure 3): between two 6-fold Ti atoms (e1) (and equivalent transfers in other layers) and between two 5-fold Ti atoms (e2) (and equivalent transfers in layers below). In addition, there are three possible corner transfers (Figure 3): across the surface between a 6-fold and a 5-fold Ti (c1), down the surface immediately below a 6-fold surface Ti (c2), and down the surface immediately below a 5-fold surface Ti (c3). For the hole transfers, in an effort to reduce the number of umbrella sampling calculations to be performed, not all possible nearest-neighbor transfers were considered. The bulk transfers D and E were not considered for the surface case as they showed the largest reorganization energies by a significant margin for the bulk case and therefore are predicted to be the slowest transfers. All possible A, B, and C hole transfers in one-half of the slab were considered. Transfers A1, B1, and C1 are those that take place between oxygen atoms that belong to TiO6 octahedra with axial ligands in the direction parallel to the surface and, conversely, transfers A2, B2, and C2 are the equivalent transfer involving oxygen atoms from TiO6 octahedra whose axial ligands lie in the direction normal to the surface. Because the polaron charge is spread over multiple sites, some fraction of the polaron charge needs to be redistributed for those cases where the polaron is located at under-coordinated surface sites. To determine the optimum charge redistribution at those sites, we carried out charge transfer calculations for electron transfer e2 (i.e., between two under-coordinated surface titanium sites) and for hole transfer A1 (i.e., between two undercoordinated surface oxygen sites) at the rutile (110) surface using the electronic structure approach described previously.19 The increase in Ti-O bond length with respect to the perfect surface without a polaron and the reorganization energy of the two transfers were used to evaluate the optimum charge redistribution (see Table 1 for the derived charges). For the electron polaron, all of the charge due to the missing sixth oxygen ligand was placed on the titanium ion. Although the Ti-O axial bond lengthens to a lesser extent in the shell model, the increase in the Ti-O equatorial bond length is in good agreement with the DFT calculations (Table 10). The reorganization energy increases, as predicted by the electronic structure calculations, although not as extensively. In the case of the surface hole polaron, the best agreement was obtained when placing 57% of the charge due to the missing third titanium ligand on the oxygen ion. This charge distribution yielded an increase in Ti-O bond length and reorganization energy of 13 and 12%, respectively, in excellent agreement with the DFT data, as shown in Table 10. The free energy profiles were obtained as follows: for each charge transfer, a 200 ps simulation was carried out after a 5 ps equilibration period for four values of θ, namely, 0.000, 0.167, 0.333, and 0.500. A configuration was collected every 0.2 ps and a value of the energy gap calculated for each configuration. As a result, each free energy profile was generated from 4000 energy gap values. All calculations were performed at 300 K. Figure 4 shows examples of free energy profiles obtained following this approach for the c2 electron transfer and equivalent transfers in other layers at the rutile (110) surface
Kerisit et al. (Figure 3). Tables 11 and 12 show the reorganization energies obtained for all electron and hole transfers, respectively. Umbrella sampling calculations were also carried out for the equivalent transfers in the bulk at 300 K as the results reported so far were for transfers at zero Kelvin. The reorganization energies of electron transfers in the top surface layer increase slightly by 7-18% but convergence to bulk values is predicted to be rapid, that is, from the second layer down. Hole transfers show on average a larger increase in reorganization energy in the top surface layer with increases from 13 to 25%. Again, convergence to bulk values is calculated to be rapid and reorganization energies not significantly different from that in the bulk are obtained from the third layer. In addition to reorganization energies, free energies of charge transfer were also obtained from the same set of calculations from the difference between the energy minima of the reactants and products free energy curves (see Figure 4 for examples). Combining all the charge transfers considered, a free energy profile versus depth was generated for both electron and hole polarons, as shown in Figure 5. Both profiles show that the free energy only varies significantly from that in the bulk in the topmost layer. The six-coordinated and five-coordinated Ti surface sites and the topmost (or surface bridging) oxygen site were found to have free energies 1.5, 2.0, and 1.5 eV higher than the free energy of the equivalent sites in the middle of the slab, respectively. The high free energies of these surface sites suggest that the topmost layer of a free rutile (110) surface will be highly depleted in electron polarons in the absence of defects and that hole polarons are unlikely to occupy surface bridging oxygen sites. Interestingly, a subsurface oxygen site was found to be energetically more favored than the bulk by approximately 1 eV. A profile of the electrostatic potential at oxygen sites versus depth was calculated using the two-dimensional electrostatic summation of Parry57,58 as implemented in the computer code METADISE59 and shows that, after cutting and before relaxing the surface, the electrostatic potential at this site is already more negative than in the bulk (Figure 6). The more negative the electrostatic potential at an oxygen site, the more energetically favored it is for this site to accept a hole polaron, that is, a positive charge. After relaxation, the electrostatic potential becomes even more negative with respect to the bulk. These results suggest that the attractive character of this subsurface site is electrostatic in nature. According to Figure 6, the surface bridging oxygen site should be attractive to hole polarons; however, as this site is under-coordinated, the hole polaron charge has to adopt a different distribution, which affects the electrostatic potential at that oxygen site. Anatase (001) Surface. We will now discuss the results obtained at the anatase (001) surface. The geometry of the anatase unit cell was optimized at constant pressure and the z axis was aligned along the [001] direction. The slab was grown to 3 unit-cell thick and the surface vectors were scaled 5 × 5 to generate a slab with a surface area of 18.94 × 18.94 Å2 and a thickness of approximately 27 Å. A 27 Å vacuum gap was introduced between the two faces of the anatase slab. All possible nearest-neighbor electron transfers were considered in one-half of the slab. Each surface layer has one titanium environment and all titanium ions exposed at the surface are in 5-fold coordination (see Figure 3). Due to the symmetry of the surface, there are two nonequivalent corner transfers which both occur within a surface layer whereas there is only one edge transfer that always takes place between two layers (see Figure 3). As was done for the rutile surface, the two hole transfers that showed the highest reorganization energies in the bulk,
Atomistic Simulation of Charge Transfer in Titania namely C and E, were not considered for the surface case, which reduced the number of transfers to be considered to four. Each layer has two nonequivalent oxygen atoms, with the topmost oxygen atom being in 2-fold coordination (Figure 3). The same procedure as described above for the rutile surface was employed to generate the free energy surfaces and thus determine the reorganization energies and free energies of charge transfer. Tables 13 and 14 show the reorganization energies obtained for all electron and hole transfers, respectively. Considering first the electron transfers, as was observed for rutile, the reorganization energy only deviates from the bulk value in the topmost surface layer with calculated increases from 7 to 22%. Similarly to what was predicted for the rutile surface, reorganization energies for hole transfers are calculated to increase, on average, to a larger extent with increases from 15 to 20% in the topmost layer. Again, the reorganization energy is seen to converge to its bulk values rapidly, from the second layer down. Figure 7 shows the free energy profiles versus depth obtained for electron and hole polarons at the anatase surface. The electron and hole polaron profiles are similar in that only the topmost sites differ significantly from the bulk free energy. In both cases, the topmost site is found to be approximately 3 eV higher in energy than a bulk site. As a result, our calculations predict a depletion of electron polarons in the topmost layer with respect to the bulk, as was found at the rutile (110) surface. Although the under-coordinated surface site is highly repulsive to hole polarons, the second oxygen atom exposed at the surface is fully coordinated and has a free energy not significantly larger than that of a bulk site. This would suggest that, as predicted for the rutile (110) surface, hole polarons will be less depleted at the surface than electron polarons. Therefore, the two surfaces show similar characteristics; however, we did not observe any stable subsurface site as was found at the rutile (110) surface. Finally, it should be noted that, at both the rutile and anatase surfaces, the large free energies of polaron transfer calculated at the topmost layer are larger than the reorganization energies. This implies that those transfers are in the so-called Marcus inverted region. In the normal region, the diabatic activation energy decreases with increasing free energy of charge transfer; the activation energy becomes zero when the free energy equals the reorganization energy. Larger free energies while keeping the reorganization energy constant are in the inverted region, wherein the activation energy increases with increasing charge transfer free energy, that is, the charge transfer slows with stronger driving force. All transfers in the inverted region are shown in bold in Tables 11–14. These include electron transfers c2 and c3 and hole transfers B1 in the first layer of the rutile (110) surface and electron transfer e and hole transfers B and F in the first layer of the anatase (001) surface. Except for electron transfer c2 from the first layer of the rutile (110) surface, all of these transfers involve charge transfer from an under-coordinated surface site. Conclusions Potential parameters for modeling electron and hole polarons in titania polymorphs were fitted to a set of reorganization energies obtained from DFT calculations. The model is based on a polarizable version of the Matsui and Akaogi rigid-ion model32 for TiO2 polymorphs whereby the oxygen polarizability is described via a shell model. The model thus derived was validated against the extent of lattice distortion of two titania crystal lattices upon introduction of electron and hole polarons and the difference in polaron formation energies between rutile and anatase. Good agreement was obtained throughout. Both
J. Phys. Chem. C, Vol. 112, No. 20, 2008 7687 the DFT calculations and shell model simulations predict that electron polarons have a greater intrinsic affinity for the anatase phase whereas hole polarons favor the rutile phase. The newly derived potential model was employed in a series of molecular dynamics simulations to investigate the intrinsic surface effects on the energetics of polaron transfer using the rutile (110) and anatase (001) surfaces as examples. Both the reorganization energy and free energy as a function of depth indicate that the surface effects are limited to the first couple of surface layers at most. However, the topmost surface layer shows large deviations from bulk behavior with increases in reorganization energy varying from approximately 5 to 25% and free energies of surface sites up to 3 eV higher than equivalent bulk sites. The high free energy of these surface sites causes many surface charge transfers to be in the Marcus inverted region. Although, by and large, the polaron free energy increases as it approaches the surface, as one would expect for a charge approaching a medium of lower dielectric constant, the electrostatic potential of some subsurface sites can be significantly lower than that of equivalent sites in the bulk and thus act as traps for diffusing electron or hole polarons, as was found for hole polarons at the rutile (110) surface. The presence of these sites may have important implications for the availability of electrons and holes in photocatalytic reactions such as the photoinduced oxidation of TMA molecules at the rutile (110),4–9 which, as described in the Introduction section, involves the capture of holes by TMA molecules at the surface. Acknowledgment. This research was supported by the U.S. Department of Energy’s (DOE) Office of Basic Energy Sciences. The computer simulations were performed in part using the Molecular Science Computing Facility in the William R. Wiley Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory (PNNL). PNNL is operated for the DOE by Battelle Memorial Institute under Contract DE-AC05-76RL01830. References and Notes (1) Linsebigler, A. L.; Lu, G.; Yates, J. T., Jr Chem. ReV. 1995, 95, 735. (2) Fujishima, A.; Rao, T. N.; Tryk, D. A. J. Photochem. Photobiol. C 2000, 1, 1. (3) Kudo, A. Catal. SurVeys Asia 2003, 7, 31. (4) Henderson, M. A.; White, J. M.; Uetsuka, H.; Onishi, H. J. Am. Chem. Soc. 2003, 125, 14974. (5) White, J. M.; Szanyi, J.; Henderson, M. A. J. Phys. Chem. B 2004, 108, 3592. (6) White, J. M.; Henderson, M. A. J. Phys. Chem. B 2005, 109, 12417. (7) Uetsuka, H.; Onishi, H.; Henderson, M. A.; White, J. M. J. Phys. Chem. B 2004, 108, 10621. (8) Henderson, M. A.; White, J. M.; Uetsuka, H.; Onishi, H. J. Catal. 2006, 238, 153. (9) Lyubinetsky, I.; Yu, Z. Q.; Henderson, M. A. J. Phys. Chem. C 2007, 111, 4342. (10) Bogomolov, V. N.; Kudinov, E. K.; Firsov, Y. A. SoV. Phys. Solid State 1968, 9, 2502. (11) Bogomolov, V. N.; Firsov, Y. A.; Kudinov, E. K.; Mirlin, D. N. Phys. Status Solidi 1969, 35, 555. (12) Yagi, E.; Hasiguti, R. R.; Aono, M. Phys. ReV. B 1996, 54, 7945. (13) Hendry, E.; Wang, F.; Shan, J.; Heinz, T. F.; Bonn, M. Phys. ReV. B 2004, 69, 081101. (14) Heluani, S. P.; Comedi, D.; Villafuerte, M.; Juarez, G. Phys. B 2007, 398, 305. (15) Holstein, T. Ann. Phys. 1959, 8, 325. (16) Holstein, T. Ann. Phys. 1959, 8, 343. (17) Yamashita, J.; Kurosawa, T. J. Phys. Chem. Solids 1958, 5, 34. (18) Jacob, I.; Moreh, R.; Shahal, O.; Wolf, A. Phys. ReV. B 1987, 35, 8. (19) Deskins, N. A.; Dupuis, M. Phys. ReV. B 2007, 75, 195212.
7688 J. Phys. Chem. C, Vol. 112, No. 20, 2008 (20) Friedman, L.; Holstein, T. Ann. Phys. 1963, 21, 494. (21) Emin, D.; Holstein, T. Ann. Phys. 1969, 53, 439. (22) Austin, I. G.; Mott, N. F. AdV. Phys. 1969, 18, 41. (23) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265. (24) Marcus, R. A. J. Chem. Phys. 1956, 24, 966. (25) Iordanova, N.; Dupuis, M.; Rosso, K. M. J. Chem. Phys. 2005, 123, 074710. (26) Iordanova, N.; Dupuis, M.; Rosso, K. M. J. Chem. Phys. 2005, 122, 144305. (27) Rosso, K. M.; Smith, D. M. A.; Dupuis, M. J. Chem. Phys. 2003, 118, 6455. (28) Kerisit, S.; Rosso, K. M. J. Chem. Phys. 2005, 123, 224712. (29) Kerisit, S.; Rosso, K. M. J. Chem. Phys. 2007, 127, 124706. (30) Liechtenstein, A. I.; Anisimov, V. I.; Zaanen, J. Phys. ReV. B 1995, 52, R5467. (31) Kerisit, S.; Rosso, K. M. Geochim. Cosmochim. Acta 2006, 70, 1888. (32) Matsui, M.; Akaogi, M. Mol. Simul. 1991, 6, 239. (33) Deskins, N. A.; Kerisit, S.; Rosso, K. M.; Dupuis, M. J. Phys. Chem. C 2007, 111, 9290. (34) Sayle, T. X. T.; Catlow, C. R. A.; Sayle, D. C.; Parker, S. C.; Harding, J. H. Philos. Mag. A 1993, 68, 565. (35) Warshel, A. J. Phys. Chem. 1982, 86, 2218. (36) Hwang, J. K.; Warshel, A. J. Am. Chem. Soc. 1987, 109, 715. (37) King, G.; Warshel, A. J. Chem. Phys. 1990, 93, 8682. (38) Rosso, K. M.; Dupuis, M. J. Chem. Phys. 2004, 120, 7050. (39) Kuharski, R. A.; Bader, J. S.; Chandler, D.; Sprik, M.; Klein, M. L.; Impey, R. W. J. Chem. Phys. 1988, 89, 3248. (40) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: Oxford, UK, 1954. (41) Swamy, V.; Gale, J. D.; Dubrovinsky, L. S. J. Phys. Chem. Solids 2001, 62, 887. (42) Dubrovinskaia, N. A.; Dubrovinsky, L. S.; Ahuja, R.; Prokopenko, V. B.; Dmitriev, V.; Weber, H.-P.; Osorio-Guillen, J. M.; Johansson, B. Phys. ReV. Lett. 2001, 87, 275501. (43) Bandura, A. V.; Kubicki, J. D. J. Phys. Chem. B 2003, 107, 11072. (44) Dick, B. G.; Overhauser, A. W. Phys. ReV. 1958, 112, 90. (45) Smith, W.; Forester, T. R. DL_POLY, a package of molecular simulation routines; The Council for the Central Laboratory of the Research Councils: Daresbury Laboratory at Daresbury, Nr. Warrington., 1996.
Kerisit et al. (46) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (47) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (48) Ewald, P. P. Ann. Phys. 1921, 64, 253. (49) Mitchell, P. J.; Fincham, D. J. Phys.: Condens. Matter 1993, 5, 1031. (50) Deskins, N. A.; Dupuis, M. Submitted. (51) Leslie, M.; Gillan, M. J. J. Phys. C: Solid State Phys. 1985, 18, 973. (52) Watson, G. W.; Wall, A.; Parker, S. C. J. Phys.: Condens. Matter 2000, 12, 8427. (53) Watson, G. W.; Tschaufeser, P.; Wall, A.; Jackson, R. A.; Parker, S. C. Lattice energy and free energy minimization techniques. In Computer Modeling in Inorganic Crystallography; Catlow, C. R. A., Ed.; Academic Press: New York, 1997; pp 55. (54) Li, G.; Gray, K. A. Chem. Phys. 2007, 339, 173. (55) Diebold, U. Surf. Sci. Rep. 2003, 48, 53. (56) Diebold, U.; Ruzycki, N.; Herman, G. S. S., A. Catal. Today 2003, 85, 93. (57) Parry, D. E. Surf. Sci. 1975, 49, 433. (58) Parry, D. E. Surf. Sci. 1976, 54, 195. (59) Watson, G. W.; Kelsey, E. T.; de Leeuw, N. H.; Harris, D. J.; Parker, S. C. J. Chem. Soc., Faraday Trans. 1996, 92, 433. (60) Rice, C. E.; Robinson, W. R. Acta Crystallogr. 1977, B33, 1342. (61) Chi, T. C.; Sladek, R. J. Phys. ReV. B 1973, 7, 5080. (62) Rimai, D. S.; Sladek, R. J.; Nichols, D. N. Phys. ReV. B 1978, 18, 6807. (63) Abrahams, S. C.; Bernstein, J. L. J. Chem. Phys. 1971, 55, 3206. (64) Rao, K. V. K.; Naidu, S. V. N.; Iyengar, I. J. Am. Ceram. Soc. 1970, 53, 124. (65) Meagher, E. P.; Lager, G. A. Can. Mineral. 1979, 17, 77. (66) Simons, P. Y.; Dachille, F. Acta Crystallogr. 1967, 23, 334. (67) Isaak, D. G.; Carnes, J. D.; Anderson, O. L.; Cynn, H.; Hake, E. Phys. Chem. Minerals 1998, 26, 31. (68) Swamy, V.; Dubrovinsky, L. S. J. Phys. Chem. Solids 2001, 62, 673. (69) Iuga, M.; Steinle-Neumann, G.; Meinhardt, J. Eur. Phys. J. B 2007, 58, 127.
JP8007865
