Oxidation-State Constrained Density Functional Theory for the Study

Subscriber access provided by BUFFALO STATE

Quantum Electronic Structure

Oxidation-State Constrained Density Functional Theory for the Study of Electron-Transfer Reactions Calvin Ku, and Patrick H.-L. Sit J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00281 • Publication Date (Web): 24 Jul 2019 Downloaded from pubs.acs.org on July 27, 2019

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

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

Page 1 of 27 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

Journal of Chemical Theory and Computation

Oxidation-State Constrained Density Functional Theory for the Study of Electron-Transfer Reactions Calvin Ku and Patrick H.-L. Sit* School of Energy and Environment, City University of Hong Kong, Hong Kong Special Administrative Region

We propose a new constrained density functional theory (CDFT) approach which directly controls the oxidation state of the target atoms. In this new approach called oxidation-state constrained density functional theory (OS-CDFT), the eigenvalues of the occupation matrix obtained from projecting the Kohn-Sham wavefunctions onto the valence orbitals are constrained to obtain the desired oxidation states. This approach is particularly useful to study electron transfer problems in transition metal-containing systems due to the multivalent nature of the transition metal ions. The calculation of the forces on the ions and of the coupling constant was implemented under the OS-CDFT scheme to allow efficient and accurate study of electron transfer reactions. We demonstrated the application of this method in the study of different electron transfer reactions including the aqueous ferrous-ferric self-exchange reaction, polaron hopping in the TiO2 anatase and bismuth vanadate, and photoexcited electron transfer in the sapphire.

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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 2 of 27

1. Introduction Density Functional Theory (DFT) has been widely used for the electronic ground-state study of molecular and extended systems. The Hohenberg-Kohn theorems1, the basis of DFT guarantees a ground-state electronic density, which can be obtained by the variation principle, for any external potential. The Kohn-Sham approach converts the many-electron Schrödinger equation to an iterative single-electron Kohn-Sham equation, with the introduction of the approximate Exchange-Correlation (XC) energy functionals. Such functionals usually adopt an analytic functional form of only the electronic density (local density approximation) or of both the electronic density and its gradient (generalized gradient approximation). Due to the mean-field nature of the typical XC functionals, DFT calculations often suffer from the self-interaction error (SIE). For example, given a H2+ system at infinite separation, H+–H, H–H+, and H0.5+–H0.5+ states are exactly degenerate. However, DFT underestimates the energy of the H0.5+–H0.5+ state compared to the H+–H, H–H+ states2-3. This results in DFT predicting the H0.5+–H0.5+ state to be the non-degenerate ground state due to the electron delocalization energy. In this regard, constrained DFT (CDFT) was introduced4-6 which provides a workaround to the SIE problem. It also allows the access to the excited-state properties which are not typically available in the DFT-based calculation6. In typical CDFT, a constraint term is added to the Kohn-Sham energy functional so that the partial charge (and optionally magnetization) around the target atom(s) is constrained to a predefined value. Weighting functions on the charge density are used to control where in the system the constraint is applied. For example, a weighting function can be set to 1 inside the Voronoi cell of the target atom(s) and 0 everywhere else. Several standard weighting functions are commonly implemented including the density-based functions, such as Becke7 and Hirshfeld8, and the atomic population-based functions, such as Löwdin9 and Mulliken10. When performing CDFT calculations, one must be careful with the choice of the weighting function 2 ACS Paragon Plus Environment

Page 3 of 27 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


and the system it applies on. For example, it was suggested that atomic population-based weighting function does not work well with the diffuse basis set6. Moreover, CDFT may fail when the atoms are too close to each other6. Therefore, more advanced CDFT-based approaches11-12 have been developed targeting the different issues. One particularly important property in the study of electron transfer reactions is the oxidation state (OS). However, the charges obtained from the density-based and atomic population-based charge analysis methods are often of non-integral values, which give ambiguous description of the OS9, 13-15. Therefore, a CDFT approach which targets explicitly the oxidation states would be desirable16-19. In this work, we offer a new CDFT approach specifically target the oxidation states of the atoms16. Similar control of the OS was previously proposed by one of the authors using a penalty functional15. Here, we adopt a constraint formalism, which provides simpler and more robust control of the OS. The focus here is the transition metal ions due to their ability to accommodate variable oxidation states. To ease with the notation, we name the new approach oxidation-state constrained density functional theory (OS-CDFT), which has been implemented in the

QUANTUM ESPRESSO

package20-21. Unlike the conventional CDFT4, 22-23

which uses the charge density as the basis for the constraint, OS-CDFT provides an alternative approach which allows direct control on the oxidation state of the target transition metal ion. The OS-CDFT scheme can also effectively generates the diabatic states of the electron transfer reactions involving transition metal ions through constraining the ions to their desired oxidation states. For OS-CDFT in the QUANTUM ESPRESSO package, the speed of OS-CDFT is about the same as that of the CDFT implementation in the QUANTUM ESPRESSO package. Depending on the systems studied, it is about 5-10 times of the time needed for standard DFT calculations. In Section 2, we first detail the formalism of the OS-CDFT implementation. We also discuss the calculation of the forces on the ions based on OS-CDFT which are useful for structural



Page 4 of 27

optimization and molecular dynamics simulations. Section 3 showcases the application of this approach in different electron-transfer problems. Finally, we conclude our work in Section 4.

2. Methodology In our previous work16, a simple, unambiguous approach was introduced for oxidation state determination of transition metal ions, which can accommodate a variable number of delectrons. Determining the number of d-electrons belonging to the ion allows us to determine the oxidation state of the ion. First, we project the occupied Kohn-Sham wavefunctions onto the atomic d-orbitals of the ion to obtain the 5  5 occupation matrix. The (h,k)-element of the occupation matrix (n) of a transition metal ion I is defined as 𝐼 𝐼 𝑒𝑙 [𝑛𝐼 ]ℎ𝑘 = ∑𝑁 𝑖 ⟨𝜙ℎ |𝜓𝑖 ⟩⟨𝜓𝑖 |𝜙𝑘 ⟩

(1)

where Nel is the number of electrons in the system (or the number of bands if smearing is used, with the terms multiplied by the smearing weights), ψi are the Kohn-Sham wavefunctions, and 𝜙ℎ𝐼 are the d-orbitals of the target ion. For the calculations with multiple k-points, the occupation matrix of each k-point is summed to obtain the resulting occupation matrix. Furthermore, for planewave calculations using ultrasoft pseudopotentials, the overlap operator S is included to account for the augmented part of the electronic density. Next, the 5  5 occupation matrix is diagonalized to obtain the eigenvalues. We refer to such eigenvalues as occupation numbers. Occupation numbers which are close to one are labelled as “fully-occupied”. The electrons in the “fully-occupied” d-orbitals are the d-electrons belonging to the transition metal ion as the full occupation originates from orbital mixing of the occupied d atomic orbitals and the orbitals of the surrounding ligands or atoms. Occupation numbers which are significantly smaller than one are labelled as “partially-occupied”. The partial occupation originates from electron donation from the surrounding ligands or atoms to the originally empty d-orbital. Therefore, these electrons do not belong to the transition metal


Page 5 of 27 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


ion. Knowing the number of d-electrons associated with the transition metal ion, the oxidation state can be determined. It should be noted that this approach can be considered as a projection method. However, instead of computing the trace of the occupation matrix as in Löwdin or Mulliken population analysis, this approach counts the number of “fully-occupied” orbitals, and thus determining the number of d-electrons belonging to the ion and the oxidation state. Based on the same principle, the oxidation state can be controlled by constraining the values of the occupation numbers, using the constraint terms added to the KS energy functional as follows min 𝐸[𝜌(𝑟⃗)] → min max 𝑊[𝜌(𝑟⃗), 𝑉] = min max[𝐸[𝜌(𝑟⃗)] + ∑𝐼 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )] 𝜌

𝜌

𝜌

𝑉

𝑉

(2)

where VI are the Lagrange multipliers, f I are the calculated occupation numbers, and 𝑓0𝐼 are the target occupation numbers. It should be note that while OS-CDFT target specifically the oxidation states, choosing the target occupation numbers still introduces certain ambiguity in the application of the approach. As will be discussed later, we introduced procedures to determine the target occupation numbers which greatly reduces the ambiguity.

2.1.

Constraint

Hamiltonian

and

Self-consistent

Field

Calculations With the modification to the KS energy functional, the Hamiltonian becomes 𝐻|𝜓⟩ → 𝐻|𝜓⟩ +

𝛿 [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )] |𝜓⟩ 𝛿𝜌

(3)

𝐼

𝛿𝜌

Using 𝜌 = ∑𝑖⟨𝜓𝑖 |𝜓𝑖 ⟩ ⇒ δ⟨𝜓| = |𝜓⟩ for all 𝜓 ∈ {𝜓𝑖 }, the second term on the right-hand side of Equation 3 becomes 𝛿 𝛿 𝛿𝜌 [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 ) ] |𝜓⟩ = [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )] 𝛿𝜌 𝛿𝜌 𝛿⟨𝜓| 𝐼

𝐼



=

Page 6 of 27

𝛿 [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 ) ] 𝛿⟨𝜓| 𝐼

𝛿𝑓 𝐼 = ∑ 𝑉𝐼 𝛿⟨𝜓|

(4)

𝐼

The derivative of the eigenvalue f I can be expressed as 𝐼 𝛿𝑓 𝐼 𝛿𝑛ℎ𝑘 = ∑ 𝑢ℎ𝐼 𝑢𝑘𝐼 𝛿⟨𝜓| 𝛿⟨𝜓|

(5)

ℎ𝑘

Therefore, Equation 4 becomes: 𝐼 𝛿 𝛿𝑛ℎ𝑘 [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )] |𝜓⟩ = ∑ 𝑉𝐼 ∑ 𝑢ℎ𝐼 𝑢𝑘𝐼 𝛿𝜌 𝛿⟨𝜓| 𝐼

𝐼

ℎ𝑘

= ∑ 𝑉𝐼 ∑ 𝑢ℎ𝐼 𝑢𝑘𝐼 ⟨𝜙ℎ𝐼 |𝜓⟩|𝜙𝑘𝐼 ⟩ 𝐼

(6)

ℎ𝑘

The Hamiltonian (Equation 3) then becomes 𝐻|𝜓⟩ +

𝛿 [∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )] |𝜓⟩ = 𝐻|𝜓⟩ + ∑ 𝑉𝐼 ∑ 𝑢ℎ𝐼 𝑢𝑘𝐼 ⟨𝜙ℎ𝐼 |𝜓⟩|𝜙𝑘𝐼 ⟩ 𝛿𝜌 𝐼

𝐼

(7)

ℎ𝑘

where 𝑢ℎ𝐼 are the unit eigenvectors corresponding to the eigenvalues f I. As the occupation matrix is symmetric, the derivative of the unit eigenvectors cancels out as shown in Section S1 of the Supporting Information (SI). Following the proof in Wu and Van Voorhis5, the existence of a maximum point in the KS energy with respect to the Lagrange multiplier can be proved as shown in Section S2 of the SI. The use of constraints on the eigenvalues of the d-orbital occupation matrix has the advantage that no restriction is placed on which d-orbital (or combination of the d-orbitals) is to have full or partial occupation. In the current implementation, the eigenvalues are sorted in the ascending order. This however introduces complication in some cases. Suppose we are to constrain the largest eigenvalue to a target value lower than its unconstrained value. The constraint will force the electrons in the orbital corresponding to that eigenvalue to migrate to other orbitals or to the neighboring atoms. At the next optimization step of the Lagrange multiplier, this eigenvalue 6 ACS Paragon Plus Environment

Page 7 of 27 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


goes down, but other eigenvalues are either unaffected or may potentially increase. This could lead to the abrupt change of the order of the eigenvalues. As the constraint is still set to apply to the largest eigenvalue, the orbital it applies to changes abruptly. This swapping could repeat indefinitely leading to non-convergence in the calculation. To handle this problem, the implementation compares the orbital (or eigenvector) to which the constraint applies in the previous optimization step with the new set of eigenvectors. The constraint is then applied to the most similar eigenvector by considering the dot products between the saved eigenvector and the new set of eigenvectors. The new eigenvector is then saved for the next step. This prevents the abrupt change of the constrained orbital, while still allowing it to change gradually. This remedy is not needed when constraining the highest eigenvalue up or when constraining the lowest eigenvalue down as no orbital swapping occurs in these cases.

2.2.

Forces on the Ions

Under the OS-CDFT scheme, we can obtain the forces on the ions for geometry optimization and molecular dynamics simulations. Taking advantage of the Hellmann-Feynman theorem, the force on the ion J is 𝜕 (𝐻 + ∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )) 𝜕𝜏⃗𝐽 𝐼 𝜕𝐻 OS-CDFT =− − 𝐹⃗𝐽 𝜕𝜏⃗𝐽

𝐹⃗𝐽 = −

(8)

where, 𝐹⃗𝐽OS-CDFT = −

𝜕 (∑ 𝑉𝐼 (𝑓 𝐼 − 𝑓0𝐼 )) 𝜕𝜏⃗𝐽 𝐼

𝜕𝑓 𝐼 = − ∑ 𝑉𝐼 𝜕𝜏⃗𝐽 𝐼

𝑁𝑒𝑙

= − ∑ 𝑉𝐼 ∑ 𝑢ℎ𝐼 𝑢𝑘𝐼 ∑ [ 𝐼

ℎ𝑘

𝑖

𝜕⟨𝜙ℎ𝐼 |𝜓𝑖 ⟩ ⟨𝜓𝑖 |𝜙𝑘𝐼 ⟩ + 𝑐. 𝑐. ] 𝜕𝜏⃗𝐽

(9)



Page 8 of 27

⃗𝐹 ⃗⃗𝐽 is the total force with constraint, and 𝐹⃗𝐽OS-CDFT is the force due to the constraint term. 𝜏⃗𝐽 is the coordinate vector of the ion J. Further details on

𝐼 𝜕⟨𝜙ℎ |𝜓𝑖 ⟩ ⃗⃗𝐽 𝜕𝜏

are included in Section S3 of the

SI.

3. Results and Discussions In this section, we discuss the application of the OS-CDFT scheme to the study of different electron transfer problems. These include the aqueous ferrous-ferric self-exchange reaction, polaron hopping in the TiO2 anatase, photoexcited electron transfer in the sapphire and polaron hopping in the bismuth vanadate.

3.1.

Aqueous Ferrous-Ferric Self-Exchange Reaction

The aqueous ferrous-ferric ions self-exchange reaction is a prototypical example of the nonadiabatic long-range electron transfer problem6, 15, 24-25. Following Sit et al.15, we adopted in this study the simulation model which includes one hexa-aqua ferrous ion and one hexa-aqua ferric ion with the two iron centers fixed at 5.5 Å apart shown in Figure 1. Ultrasoft pseudopotentials with the Perdew-Burke-Ernzerhof (PBE)26 XC functional were used. The cutoffs for the wavefunctions and the augmented charge density are 30 Ry and 240 Ry respectively. The simulation cell dimensions are 20 Å × 20 Å × 20 Å with the Γ-point sampling.


Page 9 of 27 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 1: The Ferrous-Ferric hexa-aqua ions. Fe are pink; O are red; H are white Table 1 shows the occupation numbers of the 3d-orbitals of the Fe in the optimized Fe2+(H2 O)6, Fe3+(H2O)6, and Fe2+(H2O)6–Fe3+(H2O)6 without and with constraints. Note that the occupation numbers of the “unoccupied” d-orbitals are not zero due to the finite electron donation from the lone pairs of the water ligands to the empty d-orbitals. On the other hand, the occupation numbers for the “fully-occupied” orbitals are close to but not exactly one. The deviation from unity is due to the pi-back donation to the surrounding water, and the distortion of the d-orbitals from the atomic d-orbitals in the chemical environment16. For the Fe2+(H2O)6– Fe3+(H2O)6 system without constraints, it can be seen that the self-interaction error affects the largest eigenvalue of the spin down 3d-orbital of both Fe, leading to neither of them being “fully-occupied”. To restore the correct oxidation states, two constraints were applied as follows •

f 1 is the 5th eigenvalue of the spin down 3d-orbital of Fe2+, the target 𝑓01 = 0.94

•

f 2 is the 5th eigenvalue of the spin down 3d-orbital of Fe3+, the target 𝑓02 = 0.26



Page 10 of 27

The target values above were determined according to the procedure in Section S7.1 of the SI. With the use of the constraints, the correct oxidation states of the two Fe were obtained (Table 1). Further evidence can be seen from the average distance between the iron and the surrounding oxygen atoms. Without the constraints, the average Fe–O distances for the two Fe are both 2.11 Å. With the constraints, the average Fe–O distances for Fe2+ and Fe3+ are 2.17 Å and 2.06 Å respectively. We also compared the OS-CDFT implementation with the plane-wave CDFT scheme22-23. In the CDFT calculation, we used the PBE norm-conserving pseudopotentials with the wavefunction cutoff of 80 Ry. The Hirshfeld charges can be constrained in different ways. Here we adopted two ways to study the differences in their effects. For the first one, we constrained the first Fe to +0.41 and the second Fe to +0.75, where the two values are the Hirshfeld charges of Fe in an isolated Fe(H2O)6 cluster with the total charge of +2 and +3 respectively. For the second method, we constrained the difference between the Hirshfeld charges of the two Fe(H2O)6 cluster to +1. The resulting occupation numbers are also shown in Table 1. From the table, we can see that the latter way is able to produce the desired oxidation states of the two Fe, while the former cannot.

Table 1: Occupation numbers of the Fe2+ and Fe3+ 3d-orbitals for Fe2+(H2O)6, Fe3+(H2O)6, and Fe2+(H2O)6 – Fe3+(H2O)6 without constraints, with the OS-CDFT constraints and with the CDFT constraints22-23. The occupation number concerned is the 5th spin down eigenvalue shown in bold. System

Ion Spin

Fe2+(H2O)6

Fe2+ ↑

0.99

1.00

1.00

1.00

1.00

Fe2+ ↓

0.05

0.06

0.14

0.14

0.94

Fe3+ ↑

1.00

1.00

1.00

1.00

1.00

Fe3+ ↓

0.14

0.14

0.14

0.26

0.26

Fe3+(H2O)6

and Occupation Numbers


Page 11 of 27 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


Fe2+ ↑

1.00

1.00

1.00

1.00

1.00

2+ without Fe ↓

0.08

0.09

0.19

0.20

0.52

Fe3+ ↑

1.00

1.00

1.00

1.00

1.00

Fe3+ ↓

0.08

0.09

0.20

0.21

0.52

Fe2+ ↑

0.99

1.00

1.00

1.00

1.00

Fe2+ ↓

0.04

0.05

0.13

0.15

0.94

Fe3+ ↑

1.00

1.00

1.00

1.00

1.00

Fe3+ ↓

0.14

0.16

0.17

0.25

0.26

Fe2+ ↑

0.98

0.98

0.98

0.98

0.99

Fe2+(H2O)6 – Fe3+(H2O)6 using CDFT22-23, Fe2+ ↓ constraining the first Fe to +0.41 and the Fe3+ ↑ second Fe to +0.75

0.06

0.08

0.19

0.21

0.71

0.98

0.98

0.99

0.99

0.99

Fe3+ ↓

0.05

0.07

0.18

0.20

0.31

Fe2+ ↑

0.99

1.00

1.00

1.00

1.00

Fe2+(H2O)6 – Fe3+(H2O)6 using CDFT22-23, Fe2+ ↓ constraining the difference in the charges Fe3+ ↑ of the two Fe(H2O)6 cluster to be +1

0.05

0.06

0.16

0.18

0.98

0.99

1.00

1.00

1.00

1.00

Fe3+ ↓

0.09

0.10

0.12

0.27

0.29

Fe2+(H2O)6–Fe3+(H2O)6 constraints

Fe (H2O)6 – Fe (H2O)6 using OS-CDFT 2+

3+

According to the Marcus theory, the rate for non-adiabatic electron transfer can be obtained from the following formula 𝑘𝐸𝑇 =

2𝜋 |𝐻𝑎𝑏 |2 ℏ √4𝜋𝜆𝑘𝐵 𝑇

exp [−

(𝜆+𝛥𝐺 𝑜 )2 4𝜆𝑘𝐵 𝑇

]

(10)

The key parameters are the driving force (ΔG0), the reorganization energy (λ), and the electronic coupling between the initial and final diabatic states (Hab). For the ferrous-ferric selfexchange reaction, the driving force is zero. To obtain the reorganization energy, we first performed geometry optimization with the two hexa-aqua Fe ions constrained to the electronic structure of the initial diabatic state (e.g. Fe2+ on the left and Fe3+ on the right). Then we used the resulting geometry to run a self-consistent field calculation with the two hexa-aqua Fe ions constrained to the final diabatic state (i.e. Fe3+ on the left and Fe2+ on the right). The energy 11 ACS Paragon Plus Environment


Page 12 of 27

difference of the two diabatic states involving the two hexa-aqua Fe ions is the inner-sphere reorganization energy. The outer-sphere reorganization energy was calculated using the continuum theory and the details are provided in Section S4 of the SI. The inner-sphere and outer-sphere reorganization energies were calculated to be 0.83 eV and 1.45 eV respectively. The total reorganization energy is therefore 2.28 eV. This is in good agreement with the experimental values of 2.10 eV27 The calculation of the electronic coupling (Hab) follows the procedure in Wu and Van Voorhis28, which is summarized in Section S5 of the SI. The geometry of the electrontransferring state was taken to be the relaxed geometry without any constraints. Then constraints were applied to obtain the initial and final diabatic states. The Hab term was calculated to be 16 meV. Using Equation 10, kET was found to be 667 s-1, which is also in reasonable agreement with the experimental value27 of 790 s-1. For comparison, we also calculated the electron transfer parameters using CDFT22-23. The Hab is 436 meV when constraining the charge of the first Fe to +0.41 and the charge of the second Fe to +0.75. This is significantly larger the OS-CDFT result due to the inadequacy in restoring the desired oxidation states. When constraining the difference of charges between the two Fe(H2O)6 clusters to +1, the Hab is 1.75 meV, which is in closer agreement with the OS-CDFT result.

3.2.

Polaron Hopping in TiO2 Anatase

Titanium dioxide (TiO2) has attracted great interests due to its many potential photochemical and catalytic applications14, 29-31. The motion of the photo-excited electrons and holes in the TiO2 is key to these applications. Unlike the conventional picture in which the electrons are delocalized in the conduction band, the presence of the electrons often induces lattice distortion to form polarons which are trapped at the Ti lattice sites13, 32-34. The migration of the polarons from one site to another becomes an electron transfer problem under the Marcus theory. Such 12 ACS Paragon Plus Environment

Page 13 of 27 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


polaron migration has also been studied using CDFT in MgO35-36. Here, we demonstrate the use of OS-CDFT for efficient determination of the rate of polaron hopping in the bulk TiO2 anatase. In this study, the TiO2 anatase has a tetragonal unit cell with the lattice parameters37 3.78 Å × 3.78 Å × 9.51 Å. A 4×4×2 supercell with the tetragonal cell dimensions of 15.12 Å × 15.12 Å × 19.02 Å was adopted (Figure 2). We performed the calculations using the ultrasoft pseudopotentials with the PBE functional with Γ-point sampling. The wavefunction and the augmented charge density cutoffs are 30 Ry and 240 Ry respectively.

Figure 2: TiO2 anatase structure. The directions of electron transfer are indicated. Ti are pink; O are red.



Page 14 of 27

With one additional electron, the oxidation state of one of the Ti should become +3. However, without constraints, none of the 3d-orbitals of Ti has full occupation as shown in Table 2. This suggests that the additional electron delocalizes over the whole crystal. To localize the electron, a constraint was applied to the spin up d-orbitals of the left bottom Ti shown in Figure 2. The target occupation number is 0.92 as determined in Section S7.2 of the SI. Applying the constraint prescribed above, the OS-CDFT constraint effectively controls the oxidation state of the selected Ti to Ti(III) as shown in Table 2. Furthermore, the lattice distortion due to the additional electron can be observed in elongation of the average Ti–O bond length from 1.95 Å for the neutral TiO2 crystal to 2.01 Å around the Ti(III).

Table 2: Occupation numbers of the selected Ti in the TiO2 anatase with an additional electron per supercell. A constraint is applied to the largest eigenvalue of the Ti spin up occupation matrix with the target value of 0.92 (in bold). Without the constraint, the additional electron spreads over the whole crystal. The constraint localizes the electron on a single Ti atom. System

Ion and spin

Occupation numbers

No constraint

Ti ↑

0.23

0.23

0.24

0.41

0.41

Ti ↓

0.23

0.23

0.24

0.41

0.41

Ti ↑

0.17

0.17

0.31

0.35

0.92

Ti ↓

0.13

0.16

0.16

0.30

0.33

Constrained

To calculate the electron transfer rate, we first obtained the transmission probability to determine adiabaticity of the electron transfer process. The transmission probability κ is 𝜅=

𝑜 2𝑃12 𝑜 1 + 𝑃12

(11)

where 14 ACS Paragon Plus Environment

Page 15 of 27 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


𝑜 𝑃12 = 1 − exp [−

|𝐻𝑎𝑏 |2 𝜋3 √ ] ℎ𝑣𝑛 𝜆𝑘𝑏 𝑇

(12)

hvn is the typical frequency of nuclear motion and was set to 0.11 eV for the anatase38. A transmission probability close to one indicates an adiabatic electron transfer which can be treated classically. A value close to zero indicates a non-adiabatic electron transfer which can be studied using the Marcus electron transfer rate formula (Equation 10)39. We studied the electron transfer at the room temperature for the directions [021] and [010]. Table 3 shows the calculated reorganization energy , coupling constant Hab, activation energy, and transmission probability. The results from the previous study38 which uses DFT+U to generate the diabatic states are also shown in Table 3. The transmission probability, 0.03 for [021] and 0.03 for [010], confirms that the electron transfer in both directions are non-adiabatic. While the DFT+U work38 concludes that the electron transfer is non-adiabatic, their transmission probability value is significantly higher. Finally, we calculated the rate of electron transfer at the room temperature and the results are shown in Table 3. Our values of 𝑘𝐸𝑇 are in good agreement the results in the previous work38. However, our values of both the Hab and reorganization energy are lower compared to the DFT+U work for both directions. This difference could be due to their use of TiO2 clusters instead of a periodic crystal for the coupling constant calculations. Furthermore, we also compare our value of the activation energies, defined in Section S6 of the SI, with the previous work using the RPA method 40. The RPA calculation with 35% hybrid exchange functional found the activation energy of the anatase to be 0.17 eV for the [201] direction, and 0.14 eV for the [100] direction. Due to the symmetry of anatase, [021] and [201] are the same direction, and [010] and [100] are the same direction. Therefore, their results are in good agreement with the activation energies for the [021] and [010] directions in our work.



Page 16 of 27

Table 3: The Hab, reorganization energy, transmission probability and kET obtained from our OS-CDFT of the TiO2 anatase using OS-CDFT in this work and using DFT+U38 for the directions [021] and [010]. Direction Hab (meV)

Reorganization Transmission Activation kET (s-1) energy (eV) Probability κ Energy (eV)

[021]

6.38

0.70

0.03 (non- 0.168 adiabatic)

9.09×108

This work

[010]

6.48

0.73


6.73×108

This work

[021]

30

1.19


1.73×108

Deskins and Dupuis38

[010]

20

1.21


3.68×107

Deskins and Dupuis38

3.3.

Ref.

Photo-excitation Energy in Blue Sapphire

Sapphire is a gem stone of the corundum Al2O3 structure. While pure sapphire is colorless, the presence of impurities introduces color to it. For example, titanium and iron impurities turn the crystal blue. It is generally believed that that the blue color is caused by intervalence electron transfer of the red-light absorbing Fe−Ti pairs41-43. However, alternative mechanisms were suggested44-45. Moreover, the oxidation states of the ground-state Fe−Ti pair are still debated41, 43

. Here, we apply the OS-CDFT method to provide insights into this photo-excited

intervalence electron transfer process. In this study, the sapphire has a hexagonal unit cell with the lattice parameters of 4.76 Å × 4.76 Å × 12.99 Å. We adopted a 2×2×1 supercell with the hexagonal cell dimensions of 9.52 Å × 9.52 Å × 12.99 Å (Figure 3). We performed the calculations using the ultrasoft pseudopotentials with the PBE XC functional and 1×1×2 k-point sampling. The wavefunction 16 ACS Paragon Plus Environment

Page 17 of 27 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


and augmented charge density cutoffs are 30 Ry and 240 Ry respectively. One Fe−Ti pair per supercell were introduced by replacing two neighboring Al atoms.

Figure 3: Structure of the Al2O3 corundum with the Fe−Ti pair impurity. Fe are purple; Ti are pink; Al are green; O are red We first performed geometry optimization without the OS-CDFT constraint on the Al2O3 with one Fe–Ti pair impurity. Table 4 shows the occupation numbers of the 3d-orbitals of the Fe and Ti. As shown in the table, there are 6 “fully-occupied” 3d-orbitals in the Fe and no “fully-occupied” 3d-orbitals in the Ti. Therefore, the oxidation states of the Fe and Ti ions in the ground state are +2 and +4, respectively, suggesting that photo-excited transition is from Fe2+–Ti4+ to Fe3+–Ti3+. To generate the photo-excited Fe3+–Ti3+ pair, we constrained the largest spin down occupation number of the Ti to 0.96. We also applied the constraints on the 4th and 5th largest spin down occupation numbers of the Fe to 0.98 and 1.00 respectively. The details in 17 ACS Paragon Plus Environment


Page 18 of 27

determining the target occupation number is outlined in Section S7.3 of the SI. The constrained state was obtained without structural optimization for the study of this photo-excitation process. The excitation energy is the energy difference between the Fe3+–Ti3+ and the Fe2+–Ti4+ states, which was found to be 1.56 eV. This gives the absorption wavelength of around 795 nm in the red color region of the visible spectrum. This explains the distinct blue color of the corundum crystal with Fe−Ti pair impurities. This is also in reasonable agreement with the 580-710 nm absorption range found in the photoluminescence experiments of the blue sapphire45-46.

Table 4: Occupation numbers of the 3d-orbitals of the Fe-Ti pair in the Al2O3 corundum without and with the constraints. The constraints were applied to the largest eigenvalue of the spin down 3d-orbitals occupation matrix of the Ti with the target value of 0.96 (in bold), and to the largest and second largest eigenvalues of the spin down occupation matrix of the Fe with the target values of 0.98 and 1.00 respectively (in bold). System

No constraint Fe2+-Ti4+

Constrained Fe3+-Ti3+

Ion and spin

Occupation numbers

Fe ↑

0.39

0.39

0.84

0.99

0.99

Fe ↓

0.39

0.39

0.84

0.99

0.99

Ti ↑

0.19

0.19

0.36

0.40

0.40

Ti ↓

0.19

0.19

0.36

0.40

0.40

Fe ↑

0.48

0.48

0.98

1.00

1.00

Fe ↓

0.30

0.42

0.42

0.98

1.00

Ti ↑

0.15

0.15

0.17

0.35

0.35

Ti ↓

0.16

0.16

0.36

0.37

0.96


Page 19 of 27 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


3.4.

Polaron Hopping in Bismuth Vanadate (BiVO4)

Bismuth vanadate (BiVO4) has attracted great attention in recent years as the photoelectrode material for solar fuel production due to its high stability and optical absorption in the visible spectrum47-49. However, the main bottleneck for its development lies in its low electron mobility50. Understanding its electron mobility is crucial for the development of this material for photochemical applications. BiVO4 has a tetragonal unit cell with lattice parameters 5.15 Å × 5.15 Å × 11.72 Å. A 2×2×2 supercell and a 3×3×2 supercell with tetragonal cell dimensions of 10.30 Å × 10.30 Å × 23.44 Å and 15.45 Å × 15.45 Å × 23.44 Å, respectively, were adopted in this study (Figure 4). We performed the calculations using the GBRV ultrasoft pseudopotentials with the PBE XC functional with Γ-point sampling. The wavefunction and the augmented charge density cutoffs are 40 Ry and 240 Ry respectively as in Garrity et al.51 Table 5 shows the occupation numbers of the 3d-orbitals of a selected vanadium with one added electron to the BiVO4. Without constraints, all the 3d-orbitals are not “fully-occupied” suggesting that the added electron delocalizes over the whole crystal like the case in the anatase. To localize the additional electron, we constrained the largest spin up occupation number of the V 3d-orbital occupation matrix to 0.98. The details of determining the target occupation number is found in Section S7.4 of the SI. The occupation numbers of the 3d-orbitals of the V are shown in Table 5.



Page 20 of 27

Figure 4: Structure of the 2×2×2 BiVO4 supercell. Bi are pink, V are green, O are red.

Table 5: The occupation numbers of the 3d-orbitals of the selected V in the 2×2×2 BiVO4 supercell with an additional electron per supercell without and with constraints. Without the constraint, the additional electron spreads over the whole crystal. The constraint is applied to the V ↑ 5th eigenvalue at 0.98 (in bold). System

Ion and spin

Occupation numbers

No constraint

V↑

0.42

0.45

0.48

0.48

0.51

V↓

0.42

0.45

0.48

0.48

0.51

V↑

0.40

0.41

0.41

0.47

0.98

Constrained


Page 21 of 27 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


V↓

0.24

0.35

0.38

0.38

0.45

We study here the electron transfer between two adjacent vanadium separated at 3.9 Å apart. Table 6 shows the reorganization energy () and the coupling constant (Hab) obtained. For the 2×2×2 supercell, we obtained a Hab of 54 meV and  of 1.38 eV, in good agreement with the results from the CDFT study39 of 91 meV and 1.20 eV for the Hab and , respectively, despite the difference in the constraint approach. For the 3×3×2 supercell, the same reorganization energy of 1.37 eV was obtained. However, the Hab value is 19 meV, suggesting a dependency on the size of the system. However, the difference between the Hab values for the two supercell sizes is not significant given the delicate dependence of the value on the underlying electronic structure. Figures S5 and S6 in Section S8 of the SI show the spin densities of the initial and final states of the electron transfer reaction in BiVO4 for the two supercell sizes. As shown in the figures, the smaller size of the 2×2×2 supercell leads to a different final state electronic structure. This contributes to the difference in the Hab for the two calculations. Next, we calculated the electron transmission probability 𝜅 using Equation 11 and 12, with the typical frequency of nuclear motion hvn = 276 meV39. We obtained a transmission probability of 0.42 for the 2×2×2 supercell and 0.07 for the 3×3×2 supercell. While this indicates non-adiabatic electron transfer for the 3×3×2 supercell, the value of the 2×2×2 supercell is neither close to one nor zero. Hence, neither the classical theory for adiabatic electron transfer nor the Marcus rate equation for non-adiabatic electron transfer can be used. Instead, we used the Landau-Zener theory52-54 with nuclear quantum effects55-56 to calculate the electron transfer rate. The details of the Landau-Zener theory are detailed in Section S6 of the SI. As the Landau-Zener theory becomes the Marcus rate equation in the non-adiabatic limit (κ ≈ 0)56, we also used the Landau-Zener theory for the 3×3×2 supercell.



Page 22 of 27

The activation energy and the electron transfer rate are shown in Table 6. We obtained an activation energy of 0.29 eV and kET of 3 × 108 s-1 for the 2×2×2 supercell and an activation energy of 0.32 eV and kET of 2× 107 s-1 for the 3×3×2 supercell. For comparison, we also show in Table 6 the activation energy (0.25 eV) and kET (4 × 109) from the previous CDFT work39. Due to the similar values of the Hab and λ, the results for both 2×2×2 supercell calculations are in good agreement. However, the kET value for the 3×3×2 supercell is lower.

Table 6: The Hab, reorganization energy, transmission probability, activation energy, and kET obtained from our OS-CDFT work and work of Wu and Ping39 Reorganization Transmission energy (eV) Probability

Activation Energy (eV)

kET (s-1)

2×2×2 supercell 54 constrained with OS-CDFT

1.38

0.42

0.29

3 × 108

2×2×2 supercell 91 constrained with CDFT from Wu and Ping39

1.20

0.76

0.25

4 × 109

3×3×2 supercell 19 constrained with OS-CDFT

1.37

0.07

0.32

2 × 107

System

Hab (meV)

4. Conclusion In this work, we introduced and implemented a new approach of constrained DFT called oxidation-state constrained density functional theory (OS-CDFT). This method directly controls the oxidation states of the target atoms through constraining the orbital occupations. This is particularly useful for the study of redox chemistry in transition metal systems. We demonstrated the use of the OS-CDFT scheme in the study of different electron transfer


Page 23 of 27 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


problems, including the aqueous ferrous-ferric self-exchange reaction, polaron hopping in the TiO2 anatase, photoexcited electron transfer in the sapphire and polaron hopping in the bismuth vanadate. This approach allows efficient and robust generation of the diabatic states important in electron transfer studies. The electronic coupling constant from the diabatic states can also be readily obtained. Based on OS-CDFT, the ionic forces can be calculated for structural optimization and molecular dynamics simulations. Through the calculation of the key parameters for the electron transfer reactions using OS-CDFT, important properties of electron transfer reactions including the excitation energy, reaction rates and adiabaticity can be determined. This new constrained DFT approach is a powerful tool for the study of electron transfer problems, redox reactions and excited-state properties in transition metal-containing systems. ASSOCIATED CONTENT Proof on the derivative of the eigenvalue of occupation matrix, proof on the existence of maximum of the constraint term, derivation on the forces of the constraint term, details on the outer-sphere reorganization energy of the aqueous ferrous-ferric system, obtaining the electronic coupling from OS-CDFT, equations regarding Landau-Zener theory, procedures in obtaining the target occupation numbers, and the spin densities of the initial and final states of the electron transfer reaction in BiVO4 are included in the Supporting Information. This information is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author Patrick H.-L. Sit – School of Energy and Environment, City University of Hong Kong, Hong Kong Special Administrative Region ([email protected]) ACKNOWLEDGEMENT 23 ACS Paragon Plus Environment


Page 24 of 27

We acknowledge the support by the CityU SRG Funds (Project No. 7004691 and 7004926) and the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU21305415 and CityU11304818). This work utilized the High-Performance Computer Cluster managed by the School of Energy and Environment (SEE) and the College of Science and Engineering (CSE) of the City University of Hong Kong. Use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC0206CH11357. TOC GRAPHIC

REFERENCES 1. Kohn, W.; Sham, L. J., Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133-A1138. 2. Vydrov, O. A.; Scuseria, G. E.; Perdew, J. P., Tests of functionals for systems with fractional electron number. J. Chem. Phys. 2007, 126, 154109. 3. Cohen, A. J.; Mori-Sánchez, P.; Yang, W., Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792-794. 4. Wu, Q.; Van Voorhis, T., Direct optimization method to study constrained systems within density-functional theory. Phys. Rev. A 2005, 72, 024502. 5. Wu, Q.; Van Voorhis, T., Constrained Density Functional Theory and Its Application in Long-Range Electron Transfer. J. Chem. Theory Comput. 2006, 2, 765-774.


Page 25 of 27 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


6. Kaduk, B.; Kowalczyk, T.; Van Voorhis, T., Constrained Density Functional Theory. Chem. Rev. 2012, 112, 321-370. 7. Becke, A. D., A multicenter numerical integration scheme for polyatomic molecules. J. Chem. Phys. 1988, 88, 2547-2553. 8. Hirshfeld, F. L., Bonded-atom fragments for describing molecular charge densities. Theor. Chim. Acta 1977, 44, 129-138. 9. Löwdin, P. O., On the Non‐Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals. J. Chem. Phys. 1950, 18, 365-375. 10. Mulliken, R. S., Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833-1840. 11. Plaisance, C. P.; van Santen, R. A.; Reuter, K., Constrained-Orbital Density Functional Theory. Computational Method and Applications to Surface Chemical Processes. J. Chem. Theory Comput. 2017, 13, 3561-3574. 12. Evangelista, F. A.; Shushkov, P.; Tully, J. C., Orthogonality Constrained Density Functional Theory for Electronic Excited States. J. Phys. Chem. A 2013, 117, 7378-7392. 13. Austin, I. G.; Mott, N. F., Polarons in crystalline and non-crystalline materials. Adv. Phys. 2001, 50, 757-812. 14. Bak, T.; Nowotny, M. K.; Sheppard, L. R.; Nowotny, J., Mobility of Electronic Charge Carriers in Titanium Dioxide. J. Phys. Chem. C 2008, 112, 12981-12987. 15. Sit, P. H. L.; Cococcioni, M.; Marzari, N., Realistic Quantitative Descriptions of Electron Transfer Reactions: Diabatic Free-Energy Surfaces from First-Principles Molecular Dynamics. Phys. Rev. Lett. 2006, 97, 028303. 16. Sit, P. H. L.; Car, R.; Cohen, M. H.; Selloni, A., Simple, Unambiguous Theoretical Approach to Oxidation State Determination via First-Principles Calculations. Inorg. Chem. 2011, 50, 10259-10267. 17. Jansen, M.; Wedig, U., A Piece of the Picture—Misunderstanding of Chemical Concepts. Angew. Chem., Int. Ed. 2008, 47, 10026-10029. 18. Raebiger, H.; Lany, S.; Zunger, A., Charge self-regulation upon changing the oxidation state of transition metals in insulators. Nature 2008, 453, 763-766. 19. Resta, R., Charge states in transition. Nature 2008, 453, 735. 20. Paolo, G.; Stefano, B.; Nicola, B.; Matteo, C.; Roberto, C.; Carlo, C.; Davide, C.; Guido, L. C.; Matteo, C.; Ismaila, D.; Andrea Dal, C.; Stefano de, G.; Stefano, F.; Guido, F.; Ralph, G.; Uwe, G.; Christos, G.; Anton, K.; Michele, L.; Layla, M.-S.; Nicola, M.; Francesco, M.; Riccardo, M.; Stefano, P.; Alfredo, P.; Lorenzo, P.; Carlo, S.; Sandro, S.; Gabriele, S.; Ari, P. S.; Alexander, S.; Paolo, U.; Renata, M. W., QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. 21. Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M. B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; Colonna, N.; Carnimeo, I.; Corso, A. D.; Gironcoli, S. d.; Delugas, P.; R. A. DiStasio, J.; Ferretti, A.; Floris, A.; Fratesi, G.; Fugallo, G.; Gebauer, R.; Gerstmann, U.; Giustino, F.; Gorni, T.; Jia, J.; Kawamura, M.; Ko, H. Y.; Kokalj, A.; Küçükbenli, E.; Lazzeri, M.; Marsili, M.; Marzari, N.; Mauri, F.; Nguyen, N. L.; Nguyen, H. V.; Otero-de-la-Roza, A.; Paulatto, L.; Poncé, S.; Rocca, D.; Sabatini, R.; Santra, B.; Schlipf, M.; Seitsonen, A. P.; Smogunov, A.; Timrov, I.; Thonhauser, T.; Umari, P.; Vast, N.; Wu, X.; Baroni, S., Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys.: Condens. Matter 2017, 29, 465901. 22. Goldey, M. B.; Brawand, N. P.; Vörös, M.; Galli, G., Charge Transport in Nanostructured Materials: Implementation and Verification of Constrained Density Functional Theory. J. Chem. Theory Comput. 2017, 13, 2581-2590.



Page 26 of 27

23. Brawand, N. P.; Goldey, M. B.; Vörös, M.; Galli, G., Defect States and Charge Transport in Quantum Dot Solids. Chem. Mater. 2017, 29, 1255-1262. 24. Oberhofer, H.; Blumberger, J., Charge constrained density functional molecular dynamics for simulation of condensed phase electron transfer reactions. J. Chem. Phys. 2009, 131, 064101. 25. Migliore, A.; Sit, P. H. L.; Klein, M. L., Evaluation of Electronic Coupling in TransitionMetal Systems Using DFT: Application to the Hexa-Aquo Ferric−Ferrous Redox Couple. J. Chem. Theory Comput. 2009, 5, 307-323. 26. Perdew, J. P.; Burke, K.; Ernzerhof, M., Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865-3868. 27. Rosso, K. M.; Rustad, J. R., Ab Initio Calculation of Homogeneous Outer Sphere Electron Transfer Rates: Application to M(OH2)63+/2+ Redox Couples. J. Phys. Chem. A 2000, 104, 6718-6725. 28. Wu, Q.; Van Voorhis, T., Extracting electron transfer coupling elements from constrained density functional theory. J. Chem. Phys. 2006, 125, 164105. 29. Kalyanasundaram, K.; Grätzel, M., Applications of functionalized transition metal complexes in photonic and optoelectronic devices. Coord. Chem. Rev. 1998, 177, 347-414. 30. Ni, M.; Leung, M. K. H.; Leung, D. Y. C.; Sumathy, K., A review and recent developments in photocatalytic water-splitting using TiO2 for hydrogen production. Renewable Sustainable Energy Rev. 2007, 11, 401-425. 31. Perego, C.; Revel, R.; Durupthy, O.; Cassaignon, S.; Jolivet, J.-P., Thermal stability of TiO2-anatase: Impact of nanoparticles morphology on kinetic phase transformation. Solid State Sci. 2010, 12, 989-995. 32. Emin, D.; Holstein, T., Studies of small-polaron motion IV. Adiabatic theory of the Hall effect. Ann. Phys. 1969, 53, 439-520. 33. Holstein, T., Studies of Polaron Motion: Part II. The “Small” Polaron. Ann. Phys. 2000, 281, 725-773. 34. Nowotny, J.; Radecka, M.; Rekas, M., Semiconducting Properties of Undoped TiO2. J. Phys. Chem. Solids 1997, 58, 927-937. 35. Blumberger, J.; McKenna, K. P., Constrained density functional theory applied to electron tunnelling between defects in MgO. Phys. Chem. Chem. Phys. 2013, 15, 2184-2196. 36. McKenna, K. P.; Blumberger, J., Crossover from incoherent to coherent electron tunneling between defects in MgO. Phys. Rev. B 2012, 86, 245110. 37. Anthony, J. W.; Bideaux, R. A.; Bladh, K. W.; Nichols, M. C., Handbook of Mineralogy. Mineralogical Society of America. 2001 38. Deskins, N. A.; Dupuis, M., Electron transport via polaron hopping in bulk TiO2: A density functional theory characterization. Phys. Rev. B 2007, 75, 195212. 39. Wu, F.; Ping, Y., Combining Landau–Zener theory and kinetic Monte Carlo sampling for small polaron mobility of doped BiVO4 from first-principles. J. Mater. Chem. A 2018, 6, 20025-20036. 40. Spreafico, C.; VandeVondele, J., The nature of excess electrons in anatase and rutile from hybrid DFT and RPA. Phys. Chem. Chem. Phys. 2014, 16, 26144-26152. 41. Bristow, J. K.; Parker, S. C.; A. Catlow, C. R.; Woodley, S. M.; Walsh, A., Microscopic origin of the optical processes in blue sapphire. Chem. Commun. 2013, 49, 5259-5261. 42. Bristow, J. K.; Tiana, D.; Parker, S. C.; Walsh, A., Defect chemistry of Ti and Fe impurities and aggregates in Al2O3. J. Mater. Chem. A 2014, 2, 6198-6208. 43. Eigenmann, K.; Günthard, H. H., Valence states, redox reactions and biparticle formation of Fe and Ti doped sapphire. Chem. Phys. Lett. 1972, 13, 58-61. 44. Lehmann, G.; Harder, H., Optical spectra of Di- and trivalent iron in corundum. Am. Mineral. 1970, 55, 98-105. 26 ACS Paragon Plus Environment

Page 27 of 27 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


45. Wongrawang, P.; Monarumit, N.; Thammajak, N.; Wathanakul, P.; Wongkokua, W., Oxidation states of Fe and Ti in blue sapphire. Mater. Res. Express 2016, 3, 026201. 46. Bonizzoni, L.; Galli, A.; Spinolo, G.; Palanza, V., EDXRF quantitative analysis of chromophore chemical elements in corundum samples. Anal. Bioanal. Chem. 2009, 395, 20212027. 47. Abdi, F. F.; Firet, N.; van de Krol, R., Efficient BiVO4 Thin Film Photoanodes Modified with Cobalt Phosphate Catalyst and W-doping. ChemCatChem 2013, 5 (2), 490-496. 48. Park, Y.; McDonald, K. J.; Choi, K.-S., Progress in bismuth vanadate photoanodes for use in solar water oxidation. Chem. Soc. Rev. 2013, 42, 2321-2337. 49. Sivula, K.; van de Krol, R., Semiconducting materials for photoelectrochemical energy conversion. Nat. Rev. Mater. 2016, 1, 15010. 50. Rettie, A. J. E.; Lee, H. C.; Marshall, L. G.; Lin, J.-F.; Capan, C.; Lindemuth, J.; McCloy, J. S.; Zhou, J.; Bard, A. J.; Mullins, C. B., Combined Charge Carrier Transport and Photoelectrochemical Characterization of BiVO4 Single Crystals: Intrinsic Behavior of a Complex Metal Oxide. J. Am. Chem. Soc. 2013, 135, 11389-11396. 51. Garrity, K. F.; Bennett, J. W.; Rabe, K. M.; Vanderbilt, D., Pseudopotentials for highthroughput DFT calculations. Comput. Mater. Sci. 2014, 81, 446-452. 52. Landau, L. D., A Theory of Energy Transfer on Collisions. Phys. Z. Sowjetunion 1932, 1. 53. Landau, L. D., A Theory of Energy Transfer II. Phys. Z. Sowjetunion 1932, 2. 54. Zener, C.; Fowler Ralph, H., Non-adiabatic crossing of energy levels. Proc. R. Soc. London, Ser. A 1932, 137, 696-702. 55. Brunschwig, B. S.; Logan, J.; Newton, M. D.; Sutin, N., A semiclassical treatment of electron-exchange reactions. Application to the hexaaquoiron(II)-hexaaquoiron(III) system. J. Am. Chem. Soc. 1980, 102, 5798-5809. 56. Newton, M. D., Quantum chemical probes of electron-transfer kinetics: the nature of donoracceptor interactions. Chem. Rev. 1991, 91, 767-792.


Oxidation-State Constrained Density Functional Theory for the Study

Recommend Documents