Constrained-Orbital Density Functional Theory. Computational Method

Subscriber access provided by NEW YORK UNIV

Article

Constrained-Orbital Density Functional Theory – Computational Method and Applications to Surface Chemical Processes Craig Patrick Plaisance, Rutger A. van Santen, and Karsten Reuter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00362 • Publication Date (Web): 28 Jun 2017 Downloaded from http://pubs.acs.org on June 30, 2017

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

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

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

Constrained-Orbital Density Functional Theory – Computational Method and Applications to Surface Chemical Processes Craig P. Plaisance,∗,† Rutger A. van Santen,‡,¶ and Karsten Reuter† †Chair for Theoretical Chemistry and Catalysis Research Center, Technische Universit¨ at M¨ unchen, Lichtenbergstr. 4, D-85747 Garching, Germany ‡Institute for Complex Molecular Systems, Technische Universiteit Eindhoven, Ceres building, P.O. Box 513, 5600 MB Eindhoven, Netherlands ¶Laboratory of Inorganic Materials Chemistry, Schuit Institute of Catalysis, Technische Universiteit Eindhoven, Helix building, P.O. Box 513, 5600 MB Eindhoven, Netherlands E-mail: [email protected]

Abstract

1

We present a method for performing densityfunctional theory (DFT) calculations in which one or more Kohn-Sham orbitals are constrained to be localized on individual atoms. This constrained-orbital DFT (CO-DFT) approach can be used to tackle two prevalent shortcomings of DFT – the lack of transparency with regards to the governing electronic structure in large (planewave based) DFT calculations and the limitations of semi-local DFT in describing systems with localized electrons or a large degree of static correlation. CODFT helps to address the first of these issues by decomposing complex orbital transformations occurring during elementary chemical processes into simpler and more intuitive transformations. The second issue is addressed by using the CO-DFT method to generate configuration states for multi-configuration Kohn-Sham calculations. We demonstrate both of these applications for elementary reactions involved in the oxygen evolution reaction.

Semi-local density-functional theory (DFT) is presently an indispensable tool for performing ground-state electronic structure calculations on a wide variety of systems of practical interest. It yields relatively accurate energetics for many systems at a computational cost much less than wave-function-based methods or non-local versions of DFT. Despite this success, semi-local DFT clearly also has severe shortcomings. A first shortcoming, which we will refer to as the accuracy problem, concerns the poor accuracy of semi-local DFT observed for certain systems. It is well know that large errors can be present in systems with localized electronic states (due to the self-interaction error) or significant static correlation (the inability of a single Kohn-Sham determinant to approximate the electronic structure). 1–3 There are many applications for which this can be a significant problem, for example when calculating reaction barriers 2 or fundamental band gaps, 2,4 or when performing calculations on systems such as transition metal oxides with localized electrons in the ground state. 1,5 A second shortcoming of semi-local DFT, or rather all flavors of DFT, we refer to as the transparency problem. This denotes the

Introduction

ACS Paragon Plus Environment

1


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

difficulty in extracting useful, i.e. conceptual, information about the influence of the electronic structure on elementary chemical processes. Here, the advantage of quantitatively describing and explicitly resolving the full electronic structure turns into a disadvantage. In contrast to qualitative approaches treating only a few states, DFT calculations often include hundreds of delocalized electronic states seemingly involved in an elementary chemical process, each making only a small contribution. In DFT calculations based on planewave basis sets, the situation is even worse since the wave functions of the electronic states are expressed in terms of millions of planewaves. As a result, conceptual analyses of bonding contributions or extraction of qualitative mechanistic aspects, as so successfully done in the early days of empirical theories like extended H¨ uckel, 6,7 are rarely done in present-day practical DFT work. Instead, DFT is primarily employed to produce total energies, which then determine thermochemical or kinetic quantities. This is unfortunate as the electronic structure contains a wealth of information that could potentially help us to better understand how structure and chemistry are related. Despite the seemingly very different nature of the transparency and accuracy problems, we have developed an approach that could potentially be utilized to address both of them, particularly with respect to studying chemical reactions on surfaces. This approach involves performing DFT calculations wherein certain Kohn-Sham orbitals are constrained to lie on single atoms. These orbitals, which we refer to as atomically-localized orbitals, are constructed from local basis functions centered on a particular atom and therefore resemble atomic orbitals that are not allowed to hybridize with their surroundings. We now discuss separately the potential we see for using this approach to address the transparency and accuracy problem. DFT calculations with atomically-localized orbitals can be used to address the transparency problem by allowing one to decompose complex bonding and orbital transformations occurring during elementary chemical processes into several simpler and more intuitive trans-

Page 2 of 25

formations. Elementary chemical processes in extended systems typically involve the simultaneous formation and breaking of multiple chemical bonds. The individual bonds are the result of hybridization between atomic orbitals on the two bonded atoms. By atomically-localizing the orbitals involved in a bond, we can artificially break the bond, allowing us to deconvolute the different bond transformations occurring during the elementary chemical process. Atomically-localized orbitals can be utilized to address the accuracy problem associated with semi-local DFT because of a particular property of these functionals. Despite its shortcomings, semi-local DFT actually works quite well for systems in which the exchangecorrelation hole is localized around the reference electron. 8 This happens to be fulfilled in isolated atoms, and by extension, in atoms with atomically-localized orbitals embedded in a larger system. This can be exploited to perform multi-configurational Kohn-Sham calculations 9,10 in which the configurations consist of Kohn-Sham determinants formed from an active space of atomically-localized orbitals. This is especially useful for treating only a small part of the system at a higher level since the active space can be restricted to the atomicallylocalized orbitals on the atoms in this subsystem, while the orbitals in the remainder of the system are not constrained. We have developed a constrained-orbital DFT (CO-DFT) method for performing calculations containing atomically-localized orbitals that is suitable for implementation in a planewave DFT code. It works by solving the Kohn-Sham equations under the constraint that one or more orbitals are forced to lie within restricted subspaces, such as a set of local functions centered on a single atom. The details of this method are presented in Sec. 2, along with how it was implemented in the Vienna Ab-initio Simulation Package (VASP), 11 a widely used planewave DFT code. In Sec. 3, we further discuss how this method can be used to address the transparency and accuracy problems afflicting semilocal DFT and illustrate this for reactions occurring during the oxygen evolution reaction (OER) on 3d transition metal oxides.


2

Page 3 of 25

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

2


Method for performing constrained-orbital DFT calculations

code, and to our knowledge, this has not been done before. Another deficiency of BL-DFT is that the partitioning is fixed during the calculation. In some situations, it may be desirable to optimize the local subspace itself. For instance, if we want to constrain the two hybrid atomic orbitals forming a covalent bond between two atoms, one would want to converge to the ‘best’ hybrid orbitals for this purpose (we discuss such an example in Sec. 3.1 of this manuscript). The CO-DFT method that we present herein combines aspects of both CDFT and BL-DFT. It works by solving the Kohn-Sham equations under the constraint that one or more orbitals are forced to lie within restricted subspaces. Like BL-DFT, these orbitals do not hybridize with the rest of the system; however, the remainder of the restricted subspace not spanned by the constrained orbitals can freely hybridize with the rest of the system, similar to CDFT where there is no explicit restriction on hybridization between fragment subspaces. The CO-DFT method has certain technical aspects that make it more difficult to implement in a planewave DFT code (and likely also in an atomic orbital basis code) than CDFT and BL-DFT. In CDFT, the constraint is enforced by adding a constraining potential (it can be local or non-local) to the Kohn-Sham Hamiltonian. 13,14 This potential applies equally to all Kohn-Sham orbitals in the system so that they are still eigenstates of the modified Hamiltonian. This way, the Kohn-Sham equations can still be cast in the form of an eigenvalue equation as is done in standard DFT. 11 A similar concept applies to BL-DFT methods since the Hamiltonian can be diagonalized separately in each fragment subspace. 19,20 In the proposed method, however, the constraint only applies to certain orbitals and thus requires an orbitaldependent potential. As a result, the KohnSham equations can no longer be cast in the form of an eigenvalue equation. Since the standard methods to solve the Kohn-Sham equations are based on this form, 11 we must develop an alternative method of solving them in the presence of orbital constraints. Finally, we would like to note that the current implementation of the CO-DFT method does

There are already several methods in the literature for performing DFT calculations on systems with constraints. The most wellknown of these methods is constrained DFT (CDFT), introduced by Dederichs et al. 12 and later refined and extensively applied by Wu, Van Voorhis, and coworkers. 13–15 In CDFT, the total charge and/or spin in certain regions of space is constrained to a particular value. The most commonly employed approach defines the constrained regions using a real-space weighting function, although projector-based versions have also been developed based on the Löwdin population. 13,14,16 In principle, CDFT prevents hybridization between orbitals on different fragments of the system by constraining the number of electrons (and spin) on each fragment to an integer number. While this usually works for preventing hybridization between well-separated fragments, it would not work when the fragments are strongly interacting because it is possible to constrain the number of electrons (and spin) on the fragment to an integer number and still have the fragment hybridize with the rest of the system. Nonetheless, CDFT has been used to construct configuration states for multi-configurational KohnSham calculations with some success, 17,18 although it is unlikely that the states calculated in this way are completely localized when the fragments are strongly interacting. A more appropriate class of methods for obtaining atomically-localized states are the block-localized DFT (BL-DFT) methods. 19,20 The premise of these methods is to partition the Hilbert space between the different fragments and then force the one-particle density matrix to be block diagonal in terms of this partitioning, thereby completely preventing hybridization between the fragments. BL-DFT is rather simple to implement in DFT codes using atomcentered basis functions since the partitioning is done in this basis. It would be considerably more difficult to implement in a planewave DFT


3


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

2.1

not give atomic forces in the presence of orbital constraints so that geometry optimization cannot be performed. There are two technical hurdles to including atomic forces in CO-DFT calculations. The first is due to the presence of Pulay forces that arise since the quasiatomic orbitals (vide infra) that define the restricted subspace move with the atoms. The second is that the shapes of these orbitals also change with the geometry since they are adapted to the local chemical environment. In this section, we give a detailed overview of the CO-DFT method that we have developed and implemented in the VASP planewave DFT code. To start out, we discuss in Sec. 2.1 how the locality constraints on the Kohn-Sham orbitals are defined. We then state the Lagrange functional used to optimize the electronic structure in the presence of the orbital constraints in Sec. 2.2 and derive the governing equations and gradients used in the optimization, also discussing how the optimization of the constrained and unconstrained bands can be separated. In Sec. 2.3, we present the complete iterative algorithm for minimization of the electronic structure in the presence of orbital constraints. Before beginning, we would like to note that the implementation that follows is independent of the pseudopotential scheme used in the calculations – it works with both ultra-soft 21 and norm-conserving 22 pseudopotentials as well as with the projector-augmented wave method. 23 One simply has to use the appropriate overlap and Hamiltonian operators for the pseudopotential that is used. Furthermore, neither the method itself nor the implementation discussed below are in principle specific for calculations performed in a planewave representation. The implementation could easily be adapted to other representations such as a real-space or atomic orbital representation, although the optimization procedure for the unconstrained bands may have to be altered to achieve efficient convergence.

Page 4 of 25

Wave functions of the constrained bands

The goal of the CO-DFT method is to minimize the Kohn-Sham energy with the wave functions of one or more bands constrained to lie within restricted subspaces. The restricted subspaces are defined by a set of orbitals that in principle can take an arbitrary form. In the present implementation, we define the restricted subspaces in terms of orbitals centered on a given atom. The calculation can include multiple restricted subspaces, each localized on a different atom; as such, we can perform calculations with different orbitals constrained to different atoms. The atomic orbitals used to define the restricted subspaces are generated from an initial unconstrained DFT calculation using the quasiatomic orbital (QO) method of Qian et al. 24 This method generates a set of QOs that exactly reproduce the ground state of the original planewave calculation while maintaining maximum similarity to a fitting set of original atomic orbitals that are given as input to the procedure. In practice, the atomic orbitals of the fitting set are taken as the orbitals of the neutral atom in vacuum. The QOs that are generated by the method end up being very similar to these input atomic orbitals but are slightly modified so that they are able to exactly reproduce the occupied Bloch orbitals of the planewave calculation – essentially, they are adapted to the chemical environment. Furthermore, the generation of the QOs is only as computationally demanding as a few steps of the SCF cycle. Because of these favorable characteristics, we choose this approach to define the model subspace in our method. We have implemented the QO method in the VASP code, the details of which are given in the Supporting Information. We note that other methods for defining atomic orbital projections in planewave calculations such as the maximally localized Wannier functions of Marzari and Vanderbilt 25 or the LOBSTER method of Maintz et al. 26,27 could in principle also be used to define the restricted subspaces. In general, there is a set of atoms {I} in each unit cell R that each accommodates one or


4

Page 5 of 25

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


tion. The local functions |aµ,R i on each atom I are required to be orthonormal to simplify the subsequent operations, X ˆ ν,R0 i = δRR0 haµ,R |S|a A∗iµ Aiν = δµν δRR0 .

more constrained orbitals. The wave functions of the constrained orbitals |ψa,R i on atom I of this set must lie completely within the subspace spanned by a set of local functions |aµ,R i that are centered on this atom (the notation µ ∈ I beneath the sum implies that the sum is over all local functions µ on atom I), X |ψa,R i = Cµa |aµ,R i. (1)

i∈I

(5) Local functions on different atoms are automatically orthogonal since they have no QOs in common. To maintain periodicity, the constrained orbitals should have the same shape in each unit cell – therefore, the coefficients Cµa and Aiµ are the same in all unit cells. In planewave DFT codes, it is more convenient to transform the constrained orbitals |ψa,R i to constrained bands |ψa,k i that have the same periodicity as the crystal lattice, X X |ψa,k i = eik·R |ψa,R i = Cµa |aµ,k i. (6)

µ∈I

This constraint can be expressed mathematically by defining a projection operator PÎ for this subspace, XX PÎ = |aµ,R ihaµ,R | Sˆ R

=

µ∈I

XX k

ˆ |aµ,k ihaµ,k | S,

(2)

µ∈I

R

( |aµ,k i are the local functions in reciprocal space, defined later in Eq. (7)) and requiring that the wave functions of the constrained orbitals on atom I obey PÎ |ψa i = |ψa i .

µ∈I

These are defined in terms of a set of local functions |aµ,k i, obtained by converting the local functions |aµ,R i to reciprocal space, X X ˜ i,k i, (7) |aµ,k i = eik·R |aµ,R i = Aiµ |Q

(3)

R

The overlap operator Sˆ is included explicitly in the definition of the projection operator (and all inner products) because the planewave basis is not orthonormal when either ultra-soft pseudopotentials 21 or the projector-augmented wave method 23 is used. The local functions |aµ,R i are defined as linear combinations of the orthogonalized QOs ˜ i,R i centered on atom I (the notation i ∈ I |Q beneath the sum implies that the sum is over all QOs i on atom I), X ˜ i,R i. |aµ,R i = Aiµ |Q (4)

i∈I

˜ i,k i being the set of orthogonalized rewith |Q ciprocal space QOs on atom I.

2.2

Electronic structure minimization in the presence of orbital constraints

The electronic structure of the constrained system is obtained by minimizing the Kohn-Sham energy for a set of orthonormal bands |ψr i (throughout the remainder of this section, we group the band index, k-point index, and spin index for each band into a single index) subject to the orbital constraints in Eq. (3). This is done by minimizing the following Lagrange functional: i X h ˆ s i − δrs L=F− γsr hψr |S|ψ

i∈I

The procedure for constructing the orthogonalized QOs is detailed in the Supporting Information. Defining a basis of local functions in terms of the QOs rather than using the QOs directly allows one to constrain orbitals to a subspace of the space spanned by all of the QOs on atom I. The coefficients Aiµ defining the restricted subspace are specified in the input to the calcula-

r,s

−

Xh a∈C


5

i hγa |(1 − PÎ )|ψa i + c.c. .

(8)


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

gradient in the absence of orbital constraints while the second term is the contribution due to the orbital constraints. This is apparent when one considers that PÎ = 1 in the absence of orbital constraints. The expressions for the gradients in Eqs. (9) and (10) can be rearranged into a more convenient form,

The first term is the electronic free energy of the system, given e.g. in Ref. 11, while the double sum in the second term of the definition of L contains the orthonormality constraints between all of the bands, with Lagrange multipliers γsr . These first two terms comprise the Lagrange functional in the absence of orbital constraints. The orbital constraints specified by Eq. (3) are contained in the last term of the Lagrange functional, with C defining the set of all constrained bands. Since the constraint must be satisfied at every point in the Hilbert space, a Lagrange multiplier function |γa i is required to specify this constraint. Differentiating L with respect to the wave functions and solving for the Lagrange multipliers that satisfy the constraints gives the gradients and governing equations (details are given in the Supporting Information),

1 ˆ ρˆS] ˆ |ψa i , |ga i = PÎ† [H, 2 1 ˆ ρˆS] ˆ |ψn i |gn i = (1 − PˆC )[H, 2 X − Sˆ |ψa i hga |ψn i,

(13)

(14)

a∈C

where PˆC is the projection operator onto the constrained bands, X ˆ PˆC = |ψa ihψa |S. (15)

δL =0 hδψ a | 1 ˆ ˆ 1 ˆ ρˆS] ˆ |ψa i , = [H, ρˆS] |ψa i − (1 − PÎ† )[H, 2 2

|ga i =

a∈C

The gradient in Eq. (14) for the unconstrained bands consists of two terms. The first term is associated with the energy of the unconstrained band in the absence of any update of the constrained bands. The second term accounts for the changes in the unconstrained bands resulting from updating the constrained bands – it essentially keeps the unconstrained bands orthogonal to the constrained bands when the latter are updated. While it is possible in principle to optimize the constrained and unconstrained bands simultaneously, we have found that this leads to severe oscillations in the charge density when the charge density is determined using a selfconsistent field (SCF) procedure as in Ref. 11. This problem can be solved by optimizing these two sets of bands separately. The optimization algorithm we have implemented thus alternates between the two sets of bands. First, one or more optimization steps are performed for the unconstrained bands using an SCF procedure to update the charge density while keeping the constrained bands fixed. Then, an optimization step is performed for the constrained bands using a so-called “direct” method (cf. Ref. 28) that avoids the problems with the SCF proce-

(9) δL =0 hδψ n | 1 ˆ ˆ = [H, ρˆS] |ψn i 2 1Xˆ ˆ ρˆS] ˆ (1 − PÎ )|ψn i, − S |ψa i hψa |[H, 2 a∈C

|gn i =

(10) ˆ is the Kohn-Sham Hamiltonian given where H in Ref. 11 and ρˆ is the density operator corresponding to the Kohn-Sham orbitals |ψr i with occupancies fr , X ρˆ = fr |ψr ihψr |. (11) r

ˆ ρˆS] ˆ indicates the commutator The notation [H, ˆ ˆ between H and ρˆS, ˆ ρˆS] ˆ =H ˆ ρˆSˆ − SˆρˆH. ˆ [H,

Page 6 of 25

(12)

The first term in both Eqs. (9) and (10) is the


6

Page 7 of 25

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


dure. These two steps alternate until convergence is achieved. When the optimizations of the constrained and unconstrained bands are separated, the expression for the gradients of the latter is simplified. During the optimization of the unconstrained bands, the constrained bands are frozen so that the second term appearing in Eq. (14) that contains the overlap between the gradients of the constrained bands and the wave function for band n vanishes, 1 ˆ ρˆS] ˆ |ψn i . |gn i = (1 − PˆC )[H, 2

nient to rewrite this gradient in the form 1 ˆ ˆ p,I |ψa i , [HI , ρˆCI SÎ ] + H |ga i = (21) 2 ˆ I and SÎ are the projections of the where H Hamiltonian and overlap operators onto the space of local functions on atom I to which band a is constrained,

(16)

ˆ − PˆC ). SÛ = (1 − PˆC† )S(1

(19)

a∈CI

while the density operator for the unconstrained bands is X ρÛ = fn |ψn ihψn |, (26) n∈U

Finally, the gradient for the unconstrained bands during the optimization of the constrained bands are given by the second term in Eq. (14), X ˆ n i. |gn i = − |ψa i hga |S|ψ (27)

With this substitution, we can see that the wave functions satisfy the generalized eigenvalue equation ˆ U |ψn i = εn SÛ |ψn i . H

(23)

The density operator for the constrained bands on atom I (contained in set CI ) is X fa |ψa ihψa |, (25) ρˆCI =

where εn is the energy eigenvalue of band n. By noting that |ψn i is already orthogonal to PˆC , we can then insert a factor of (1 − PˆC ) to the right of the Hamiltonian and overlap operators and define effective Hamiltonian and overlap operaˆ U and SÛ for the unconstrained bands, tors H (18)

SÎ = PÎ† SˆPÎ ,

ˆ p,I = Pˆ † (1 − SˆρÛ )H(1 ˆ − ρÛ S) ˆ PÎ − H ˆ I . (24) H I

(17)

ˆ U = (1 − Pˆ † )H(1 ˆ − PˆC ), H C

(22)

ˆ p,I is a term that accounts for the changes and H in the unconstrained bands required to maintain orthogonality with the constrained bands,

Using Eq. (11) to expand ρˆ in terms of the wave functions, along with the property that the wave functions are mutually orthogonal, it can be seen that this gradient vanishes when the wave functions of the unconstrained bands obey the equation ˆ − εn S) ˆ |ψn i , (1 − PˆC )(H

ˆ I = Pˆ † H ˆ PÎ , H I

a∈C

(20)

As mentioned before, these correspond to the changes in the unconstrained bands required to keep them orthogonal to the constrained bands as the latter are updated.

Since this equation has the same form as the typical Kohn-Sham equations, we can use the usual methods employed in planewave DFT codes (cf. Ref. 11) to optimize the unconˆ U and SÛ instead of H ˆ strained bands, using H ˆ and S. The gradient for the constrained bands given by Eq. (13) is unchanged by optimizing the two sets of bands separately. It is however conve-

2.3

Complete electronic structure minimization algorithm

The complete electronic structure minimization algorithm consists of the following steps:


7


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

(1) SCF procedure for unconstrained bands (perform up to NSCF times)

Page 8 of 25

subspace. These two methods and the modifications required to accommodate constrained bands are discussed in the Supporting Information. At the end of each SCF cycle for the unconstrained bands, the occupancies of these bands are recomputed (1b) and the charge density ρout,1 is constructed from the new wave functions (1c). This charge density is then mixed with charge densities from previous iterations to determine a new input charge density ρ using the Pulay method as discussed in Ref. 11 (1d). This method builds up an approximation of the inverse charge dielectric matrix G that is used to update the charge density according to

(a) Iterative update of unconstrained bands (b) Calculate occupancies of unconstrained bands (c) Calculate charge density ρout,1 for current wave functions (d) Update charge density ρ using the Pulay method (Eq. (28)) (e) Check convergence, go to (2) if converged (2) Direct procedure for constrained bands (a) Check convergence of total energy, exit electronic structure minimization if converged (b) Iterative update of constrained bands

ρ → ρ + G · (ρout,1 − ρ).

(28)

In order for the above alternating update scheme to work efficiently, the matrix G is retained during the update of the constrained bands for use as the initial G in the next set of SCF cycles for the unconstrained bands. Finally, convergence of the unconstrained bands is checked (1e) using the criteria for the original electronic structure minimization scheme in VASP (Refs. 11 and 28). The so-called “direct” procedure for optimizing the constrained bands involves a simultaneous update of all constrained bands using a conjugate gradient procedure followed by the update of the charge density. After convergence of the total energy is checked (2a), an iterative update (2b) is performed for the constrained bands, with the wave functions of the unconstrained bands also being updated to maintain orthogonality to the former. This update is performed in the basis of local functions |aµ,k i (Eq. (7)) and is discussed in further detail in the Supporting Information. The charge density ρout,2 of the new wave functions is then calculated (2c) and the input charge density for the SCF procedure (1) is updated (2d) according to

(c) Calculate charge density ρout,2 for current wave functions (d) Update charge density using Eq. (29) The procedure for the unconstrained bands (1) is repeated until convergence is reached or NSCF iterations have been performed. The algorithm then moves to the procedure for the constrained bands (2), which is executed only once before returning to (1). We have found that NSCF = 4 gives good convergence at minimal computation cost for a variety of different systems. The entire algorithm is terminated when the change in the total energy of the system between successive executions of (2) is below a specified threshold. The SCF procedure for the unconstrained bands is essentially the same procedure (without orbital constraints) used already in the VASP code, 11 the only difference being that the iterative update of the wave functions (1a) is modified for the presence of the constrained bands. Two of the procedures used in VASP for performing the iterative update of the unconstrained bands have been modified to include orbital constraints – the Davidson-block iteration scheme and the residual minimization method – direct inversion in the iterative

ρ → ρ + ρout,2 − ρout,1 .


8

(29)

Page 9 of 25

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


Using atomically-localized orbitals to improve the transparency and accuracy of DFT via the CODFT method

and Hamiltonian and only include the meanfield electron-electron interactions. One class of commonly used approaches that fall into this category is based on the projected density of states, three of these being the crystal orbital overlap population (COOP) analysis, 6,7,33 the related crystal orbital Hamiltonian population (COHP) analysis, 34 and the d-band center model. 35 While these methods have been quite successful in rationalizing a variety of chemical phenomena, 36–40 they do not capture changes in the electron-electron interactions that occur during the course of an elementary chemical process. As a result, they can fail for processes that involve large changes in electron density such as reactions involving oxidation/reduction of metal centers in transition metal oxides. Analysis of elementary chemical processes involving large changes in electron-electron interactions requires an interacting analysis technique. The energy decomposition analysis (EDA) techniques are a widely employed class of such methods that are typically used to analyze intermolecular interactions. These approaches decompose the interaction energy by performing calculations in which the wave functions are restricted to certain variational subspaces (e.g. using block-localized DFT); as such, these calculations have some similarities to our CO-DFT method. While the EDA techniques have been successful for describing intermolecular interactions, they were not intended for analyzing bonds with strong covalent character, particularly those on surfaces like we discuss below. We now discuss how calculations with atomically-localized orbitals performed using the CO-DFT method can be used to decompose the energies of elementary chemical processes on surfaces. Unlike those discussed so far, this technique is applicable to processes involving strong chemical bonds and large changes in electron populations. To motivate this approach, we first consider the case of chemisorption of an adsorbate on a surface. This phenomenon lies at the heart of all heterogeneous catalytic processes, thus a thorough theoretical understanding of it is vital for gaining fundamental insights into these processes.

We now discuss how the constrained orbital DFT method can be used to address the transparency and accuracy problems afflicting semilocal DFT. The first of these involves the deconstruction of orbital transformations occurring during chemical reactions. The second is for performing multi-configurational KohnSham calculations with a local active space. All of the constrained and unconstrained DFT calculations in this section are performed in the VASP code 11 using the RPBE generalized gradient approximation (GGA) exchangecorrelation functional. 29 Unless otherwise specified, a +U correction 30 is applied to the 3d electrons on transition metal atoms, with the values of the U parameter (listed in the Supporting Information) determined using a linear response technique that we have employed in previous work. 31 In all calculations, the projector augmented wave method 23,32 was used to describe the wave functions in the core region. Further computational details are given in Ref. 31.

3.1

Improving the transparency of DFT – deconstructing orbital transformations occurring during elementary chemical processes

Many techniques have evolved over the years for peering inside the “electronic structure black box” of complex DFT calculations. The general goal of these techniques is to break down the energy of an elementary chemical process into components that have intuitive connections with the electronic structure. These can be divided into non-interacting and interacting techniques. The former are based on an analysis of the one-electron Kohn-Sham orbitals


9

Journal of Chemical Theory and Computation A

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

M

M M

M

M M

A

M

M M

M M

M M

Δୟୢୱ

M

M M

M

M M

M

M

Δ୪୭ୡ

M M

M

M

M

orbitals) that directly bond to the adsorbate. This can be done by constraining this orbital to be atomically-localized in the “localized” state. We refer to this as an in-place localization of the orbital since the metal atom is not physically removed from the surface. We have carried out this decomposition of chemisorption energies to transition metal surfaces using the CO-DFT method and will present this study in a future publication. For now, we focus on using the CO-DFT method to decompose chemisorption energies to transition metal oxide surfaces. For chemisorption on transition metal oxide surfaces, the decomposition in Fig. 1 is not as useful, because for these systems it is more intuitive to cleave the surface metal-oxygen bonds heterolytically rather than homolytically. This would lead to charged surfaces in the intermediate states in Fig. 1, which are not welldefined in periodic calculations. Using the CODFT method, these bonds can be heterolytically cleaved in place so that no long-range charge separation is introduced.

M M

A

M

M

M

Δୣ୫ୠ

A M

M

M M

M M

M M

Δୡ୭୴

M

M M

M

M M

M

M M

M M

Page 10 of 25

M M

Figure 1: Conceptual decomposition of a surface chemisorption process into three contributions: (1) localization, (2) metal-adsorbate covalent bond formation, and (3) re-embedding. Adapted from Ref. 41 Chemisorption typically involves the simultaneous formation of surface-adsorbate bonds and breaking (or weakening) of surface-surface bonds. 41 To isolate these different contributions, chemisorption on metal surfaces can be decomposed into the three separate steps shown in Fig. 1. 41 In the first step, a metal atom is removed from the surface, during which the bonds with neighboring metal atoms are broken and electrons are localized on this atom (referred to as localization). A covalent bond is then formed between the isolated metal atom and the adsorbate in the second step to form a covalently bonded surface molecule complex. In the third step, the surface molecule complex is re-embedded into the surface, whereby metalmetal bonds broken in the first step are partially re-formed. A drawback of this decomposition is that completely removing the metal atom from the surface introduces far more disruption to the electronic structure than occurs during chemisorption. For example, chemisorption on a transition metal atom may only involve a single d orbital whereas removing the atom from the surface (corresponding to localization in Fig. 1) ruptures the bonds associated with all five d orbitals. Therefore, the localization energy as defined in Fig. 1 is not a good measure of the strength of the surface bonds that must be broken to form the chemisorption bond. A better way to define the localization energy would be to localize only the metal orbital (or

X H

H

OH

O

H

H

X

O M4+

HO

OH

H

O M4+

O H

M3+

M3+ O H

Figure 2: Dissociation of a molecule H – X over two oxygen atoms of a dual-center active site on a reducible metal oxide surface To illustrate this, we examine the reaction depicted in Fig. 2, which involves the dissociation of a molecule H – X over two oxygen atoms on a reducible metal oxide surface, with the transfer of two electrons to metal cations. Catalytic processes occurring through such an elementary reaction step include the oxygen evolution reaction 31 (H2 O addition, X = OH) and various hydrocarbon dehydrogenation/partial oxidation reactions 36,42–44 (C – H bond activation, X = alkyl). This step often occurs by the abstraction of a hydrogen atom from the reactant molecule by one surface oxygen followed by the formation of a bond between the remaining radical fragment of the reactant molecule and a


10

Page 11 of 25

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


second surface oxygen H−X −−→ H• + X•

Although it would be possible to separate step {2c} into these two components, it does not provide much additional insight because the reembedding energy is very low in this reaction due to the large energy separation between the O – H bonding and antibonding orbitals and the Fermi level. The O – X bond formation process in step {3} can be decomposed in an analogous way. In the intermediate states formed by steps {2a} and {2b}, the pair of atomically-localized spin-orbitals is constrained to a subspace composed of the O 2s and three 2p quasiatomic orbitals. The pair of spin-orbitals is optimized within this subspace to the form that minimizes the total energy. This is an example of a case in which the CO-DFT method has an advantage over the simpler BL-DFT methods (discussed at the start of Sec. 2) due to the fact that optimization of the constrained pair of spin-orbitals would not be possible with the latter method. In the (O) state formed in step {2a}, both of these orbitals are occupied, while in the ( •O) state formed in step {2b}, both of these orbitals are half-filled (fa = 0.5). Alternatively, one of these orbitals could be filled and the other empty in this state, but this introduces a large magnetic moment on the oxo that is not present in any of the other states. Decomposing the O – H and O – X bond formation steps according to steps {2a} – {2c} leads to a set of three simple orbital transformations, illustrated in Fig. 3, each occurring within a subspace dominated by two atomic orbitals. In step {2a}, the coordinative bonds between an O sp orbital (sp refers to a hybrid orbital containing an arbitrary mix of s and p orbitals) and the neighboring metal cations are broken, with the lone pair localizing on the oxo. In step {2b}, an electron is transferred from the doubly-occupied atomically-localized O sp orbital to an empty state composed mainly of a 3d orbital on a metal cation. The singly-occupied O sp orbital then forms a covalent σ-bond with the singly-occupied orbital on the H• or X• . The energies of the overall O – H bond formation step (and similarly for the O – X bond formation step) can be decomposed into the three components corresponding to the decomposi-

{1}

H• + O−Mn −−→ H−O−Mn−1

{2}

X• + O−Mn −−→ X−O−Mn−1

{3}

During this reaction, two metal cations are formally reduced. It is also possible that the two oxygens are bound to the same metal cation and its oxidation state is decreased by two levels (cf. the single-Co site in Ref. 31). In steps {2} and {3}, a bond is formed between an H• or X• radical and a surface oxygen atom, accompanied by the transfer of an electron from the oxygen to a neighboring metal cation. Both of these steps involve complex bonding and orbital transformations that are easier to comprehend if these steps are decomposed into simpler sub-steps by employing atomically-localized orbitals with the CO-DFT method. By localizing a pair of spin-orbitals |ψa i on the oxygen atom, step {2} can be broken down into three sub-steps O−Mn

−−→

(O)−Mn

{2a}

−−→ ( •O)−Mn−1

{2b}

H• + ( •O)−Mn−1 −−→ H−O−Mn−1

{2c}

(O)−Mn

In the first sub-step, {2a}, the pair of spinorbitals |ψa i is atomically-localized onto the oxygen atom (indicated by placing () around the element symbol), with each occupied by one electron. This step corresponds to an in-place heterolytic cleavage of the bonds between |ψa i and the neighboring metal cations and localization of a lone pair into these spin-orbitals. In step {2b}, an electron is transferred from the localized lone pair to an empty d orbital on the Mn cation. The combination of steps {2a} and {2b} corresponds to the homolytic cleavage of the metal-oxygen bonds involving |ψa i (cf. ∆Eloc in Fig. 1). Finally, in step {2c}, the singly-occupied pair |ψa i forms a σ-bond with the H atom. This step combines both the surface molecule complex formation and reembedding steps in Fig. 1 (∆Ecov + ∆Eemb ).


11


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 12 of 25

ϮĐ

Δ୉୘

Δୡ୭୴

Figure 3: Decomposition (according to Rxns. {2a} – {2c}) of the orbital transformations occurring during O – H bond formation on a 3d transition metal oxide surface: (2a) localization of an electron pair on the oxo, (2b) transfer of an electron from the oxo to a metal cation, and (2c) formation of the O – H bond. Occupied (unoccupied) orbitals are dark (light) blue. Only the spin-β orbitals are shown, but similar transformations occur in steps 2a and 2c for the spin-α orbitals. perimental applied potentials on pure Co3 O4 . 31 Since the current implementation of the CODFT method does not allow for the calculation of atomic forces in the presence of orbital constraints, we use fixed geometries taken from unconstrained DFT calculations. The geometry of O – Mn is used for the CO-DFT calculation on state (O) – Mn produced by step {2a} since both states have metal cations in the same oxidation state. Likewise, the geometry of H – O – Mn – 1 , with the H atom removed, is used for the CODFT calculation on state ( •O) – Mn – 1 produced by step {2b}. The three energy contributions from Eq. (30) as well as the total energies for both the O – H and O – OH bond formation steps are plotted in Fig. 5 for the five metals examined. The total O – H bond formation energy follows a double-volcano trend with maxima at Mn and Co and a minimum at Fe that is characteristic of many properties of the 3d transition metal oxides. 36,45–47 The total O – OH bond formation energy has only a single maximum for Mn. We can immediately see that the covalent bond formation energies ∆Ecov (step {2c}) do not vary much between the different metals, meaning that the differences in the energies of the bond formation steps {2} and {3} (and the total reaction energy) between these sites arise almost entirely from differences in the energy to localize a hole on the O atom (∆Eloc +∆EET ). The qualitatively different behaviors of the O – H and O – OH bond formation energies are therefore unrelated to the nature of the adsorbate but are entirely due to intrinsic properties of the oxide surface – in particular,

tion of the orbital transformations, ∆EO−H = ∆Eloc + ∆EET + ∆Ecov .

(30)

Since each of these energetic components is associated with an intuitive orbital transformation, they can be used to associate differences in bond formation energies with differences in electronic structure. ;ĂͿ

;ďͿ

Figure 4: (a) Dual-center and (b) single-center active sites on the (311) and (110)-B surfaces of Co3 O4 We now demonstrate the approach just described for the water addition reaction occurring during the OER on active sites consisting of different 3d transition metal cations (M = Cr – Ni) embedded in the Co3 O4 surface. The active site we use consists of two redox-active 3d transition metal cations embedded in the (311) surface of Co3 O4 (dual-metal site), shown in Fig. 4a. In our previous work, this site was found to be the most active site at typical ex-


12

Page 13 of 25

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

Ϭ͘ϲ

Ϭ͘ϴ

Δ୪୭ୡ

Δ୓ିୌ

Ϭ͘ϰ

Ϳ s Ğ; Ǉ Őƌ Ğ Ŷ

Δ୪୭ୡ Δ୓ି୓ୌ

Ϭ͘ϲ

Ϭ͘ϰ

Ϳ Ϭ͘Ϯ s ;Ğ Ǉ Őƌ Ϭ͘Ϭ Ğ Ŷ ͲϬ͘Ϯ

Ϭ͘Ϯ

Ϭ͘Ϭ

Δୡ୭୴

ͲϬ͘Ϯ

Δ୉୘

H

H O

OH

H

Δୡ୭୴ Δ୉୘

O M4+

ͲϬ͘ϰ

H

M3+

ϯ ƌ;,^Ϳ

ϰ DŶ;,^Ϳ

ϱ &Ğ;,^Ϳ

ϲ Ž;>^Ϳ

O

H

O M3+

ͲϬ͘ϰ

O H ͲϬ͘ϲ

OH OH

M4+ O H

ͲϬ͘ϲ

ϳ Eŝ;>^Ϳ

ϯ ƌ;,^Ϳ

ϰ DŶ;,^Ϳ

ϱ &Ğ;,^Ϳ

ϲ Ž;>^Ϳ

ϳ Eŝ;>^Ϳ

Figure 5: Energy contributions (eV) in Eq. (30) for (left) O – H bond formation to an η-hydroxo and (right) O – OH bond formation to a μ3 -oxo on a dual-metal active site embedded in the Co3 O4 (311) surface. The spin state for each transition metal cation, high spin (HS) or low spin (LS), is indicated. All energies have been referenced to the Ni-doped site for better visualization. that the oxygen forming the O – H bond is part of an η-hydroxo while the oxygen forming the O – OH bond is a μ3 -oxo. The hole localization energy itself is composed of two parts, the energy to localize a lone pair on the oxo (∆Eloc ) and the energy to transfer one electron from this lone pair to an orbital centered on the neighboring metal cation (∆EET ). The lone pair localization energy is seen to increase monotonically from Ni to Cr, likely due to the increasing oxophilicity of the metal centers as the number of 3d electrons decreases. The extrema seen in the total O – H and O – OH bond formation energies are therefore due to the variations in the electron transfer energy (∆EET ), which displays extrema at the same positions. The electron transfer energy is primarily controlled by the reducibility of the metal centers, which is a complex function of the electronic structure. In particular, the reducibility is likely to be influenced by the spin state of the metal center, which is indicated in Fig. 5 where we can see that it changes from high spin to low spin between Fe and Co. Even for O – OH bond formation to the μ3 -oxo, we see a change in behavior of ∆EET between Cr – Fe and Co – Ni that is likely associated with this change in the spin state, with a volcano-shaped trend for the former and a flat trend for the

latter. Since the nature of the bond being formed with the oxide surface (O – H or O – OH) has almost no effect on the bond formation energies, we might expect these trends to be quite general for any O – X bond on the 3d transition metal oxide surfaces. This could explain the commonly observed double-volcano trend observed in many calculated and experimentally measured quantities across the 3d transition metals – that it is controlled by the reducibility of the metal center, supporting the hypothesis that this behavior results from the high spin state adopted by the 3d transition metal cations. 36,45–47 It would be interesting to investigate other reactions on these materials using this approach to confirm whether this is the case in general.

3.2

Improving the accuracy of DFT – a multi-configurational Kohn-Sham method

A second promising application of atomicallylocalized orbitals and the CO-DFT method is to be able to perform DFT calculations with a multi-determinantal Kohn-Sham wave function. As mentioned in the introduction, semilocal DFT gives very poor results for cer-


13


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

tain systems, due to the local model of the exchange-correlation hole on which the corresponding functionals are based. This model fails in systems where the actual exchangecorrelation hole is non-local, 3,4 two classic examples being stretched H2 and H2 + . 1,2,48 In both of these systems, the exchange hole is completely delocalized over both H atoms and is independent of the position of the reference electron. 49 In stretched H2 , the correlation hole is also completely delocalized over both H atoms and changes sign when the reference electron moves from one atom to the other. These phenomena, which we refer to as inter-atomic exchange and static correlation, are not described by semi-local exchange-correlation functionals, resulting in a qualitatively incorrect description of these species in the dissociation limit. In contrast, these phenomena are accurately treated by correlated wave function methods. In turn, the correlated wave function methods suffer from a high computational cost to describe intra-atomic (dynamic) correlation, which semi-local DFT describes quite well and efficiently. 10 It therefore seems logical to combine the two methods in a way so that the semi-local DFT captures the dynamic correlation and intra-atomic exchange, while the correlated wave function method captures the static correlation and inter-atomic exchange. One way of combining these methods is the multiconfigurational Kohn-Sham method (MC-KS). In this approach, the multi-electron wave function |Ψi is constructed as a linear combination of configuration states |Φi i, each being the interacting wave function corresponding to a single-determinant Kohn-Sham wave function |ΦKS i i, X |Ψi = ci |Φi i. (31)

Page 14 of 25

inter-atomic contributions, treating the first exclusively with semi-local DFT and the latter exclusively with the correlated wave function method. The separation is achieved by defining an active space consisting of atomicallylocalized orbitals a ∈ C which are frozen in all configurations. Each configuration i in Eq. (31) is defined by which of these atomically-localized orbitals are occupied and which are unoccupied. The Kohn-Sham configuration wave function can be expressed in terms of creation operators cˆ†a and cˆ†n,i for the occupied atomically-localized and unconstrained orbitals (in sets Ci and Ui ), respectively, ! ! Y Y † |ΦKS cˆ†a cˆn,i |0i , (32) i i = a∈Ci

n∈Ui

( |0i is the vacuum state). The method presented here was developed only for use with Γ-point DFT calculations, so the k-point index is omitted, although it could in principle be extended to calculations with multiple kpoints. The wave functions of the atomicallylocalized orbitals are the same in all configurations while those of the unconstrained orbitals differ between configurations. Additionally, the unconstrained orbitals are required to be orthogonal to both the occupied and unoccupied atomically-localized orbitals. We can see that this form of the configuration wave function corresponds to a CO-DFT calculation in which the occupied constrained orbitals (in set Ci ) have occupancies of fa = 1, while the unoccupied constrained orbitals have occupancies of fa = 0; also, in these CO-DFT calculations the constrained orbitals are frozen (there is no optimization of the constrained orbitals, they are given as input to the calculation). As demonstrated in the Supporting Information, the Kohn-Sham configuration wave functions defined by Eq. (32) are mutually orthogonal. The above approach works for separating intra- and inter-atomic exchange and correlation within the active space because the interatomic effects are due to hybridization between atomic orbitals on different atoms, which is suppressed in atomically-localized orbitals. As a result, the configuration states used for the

i

A major challenge in the MC-KS approach is devising a way to avoid double-counting any contributions to the exchange-correlation energy. 10 Here, we introduce an approach for carrying out these types of calculations that utilizes atomically-localized orbitals and the CODFT method. It involves separating the exchange and correlation energies into intra- and


14

Page 15 of 25

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


correlated wave function expansion in Eq. (31) contain most of the intra-atomic exchange and correlation, while containing little of the interatomic exchange and correlation. The latter effects arise from coherence between the different configuration states. In this way, there should be minimal double-counting of the exchangecorrelation energy when combining the two methods. This approach has similarities to the CDFTCI method of Wu et al., 17,18 the main difference being that their method employs configurations with constrained charge densities while ours employs configurations with constrained orbitals. As mentioned in Sec. 2, constraints on the charge density (CDFT) work reasonably well to localize orbitals on weakly-interacting fragments. On strongly-interacting atoms or fragments, however, we expect this method to perform less satisfactory, while the CO-DFT method will not be affected. The energy of the multi-configurational wave function in Eq. (31) is given by P ∗ ˆ hΨ|H|Ψi i,j ci cj Hij E= , (33) = P ∗ ˆ hΨ|S|Ψi i,j ci cj Sij

Kohn-Sham wave functions. The matrix elements of the overlap operator are approximated by

ˆ is the exact Hamiltonian and where H

We expect this to be a reasonable approximation since the Kohn-Sham wave functions are mutually orthogonal and atomically-localized; therefore, the terms arising from the twoelectron part of the exact Hamiltonian are expected to be small, as shown in the Supporting Information. The procedure used for evaluating the integrals in Eq. (39) is also given in the Supporting Information. Once the configuration states and the Hamiltonian matrix elements between them have been calculated, the expansion coefficients ci in Eq. (31) can be determined from the energy in Eq. (33) according to the variational principle, leading to the following eigenvalue equation: X Hij cj = E ci , (41)

ˆ ji , Hij = hΦi |H|Φ

(34)

ˆ ji . Sij = hΦi |S|Φ

(35)

ˆ KS Sij ≈ hΦKS i |S|Φj i = 0,

(36)

Sii = 1.

(37)

(38)

which is zero since the Kohn-Sham wave functions in Eq. (32) are mutually orthogonal. For the off-diagonal Hamiltonian matrix elements, we cannot follow their approach exactly because, unlike in CDFT, the interacting wave functions in CO-DFT are not eigenstates of ˆ + VˆC due to the fact that the constraining poH tential VˆC in CO-DFT is applied to the KohnSham orbitals and not the density. Instead, these matrix elements are evaluated using the Kohn-Sham Hamiltonians (rather than the constraining potentials used with CDFT) for the two configuration states i and j, ˆ ij KS Hij ≈ hΦKS i |HKS |Φj i ,

i 6= j,

(39)

ˆ ij is the average of the Kohn-Sham where H KS ˆ j for states i and j, ˆ i and H Hamiltonians H KS KS î + H ˆj ) ˆ ij = 1 (H H KS KS 2 KS

ˆ are equal to The diagonal matrix elements of H KS the Kohn-Sham energies Ei of the configurations while those of Sˆ are unity, Hii = EiKS ,

i 6= j,

It is not possible to evaluate the off-diagonal matrix elements directly since the interacting wave functions |Φi i are not known, only the noninteracting Kohn-Sham wave functions |ΦKS i i. We therefore follow an approach inspired by that of Wu et al. used within the framework of CDFT 17,50 and approximate the off-diagonal matrix elements in terms of the

(40)

j

We would like to point out that in this approach the configuration wave functions are determined by applying the variational principle


15


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

to the configuration energies EiKS , not the total energy E. This essentially corresponds to the perfect-screening limit where the unconstrained orbitals respond instantly to screen any changes in occupancy of the constrained orbitals. This limit becomes exact when the constrained orbitals are only weakly coupled amongst themselves and to the unconstrained orbitals. A more rigorous treatment would determine the configuration wave functions by applying the variational principle to the total energy in Eq. (33). This introduces coupling between the configuration states during the CO-DFT calculations. We have conceptually developed a method for performing these types of calculations and will report on this in a forthcoming publication. Finally, we apply the approach just described to improve on the semi-local DFT energies of the water addition step of the OER that was discussed in the preceding section. Semi-local DFT is known to perform very poorly for systems with localized 3d electrons such as the metal oxides we examined for this reaction. 5 This is believed to be due to the non-local form of the exchange hole. The GGA+U method that we utilized in the last section improves upon normal GGA calculations, but there are still significant differences when compared to more accurate CCSD(T) calculations. 51 Furthermore, there is considerable uncertainty as to what value to use for the effective on-site Coulomb integral U . These systems are therefore a good candidate for testing our multiconfigurational Kohn-Sham approach. We specifically apply this method to the two electrochemical oxidation steps preceding the water addition reaction on the dual-Co active site on the (311) surface discussed in the previous section (Fig. 4a). We also examine this reaction on a single-Co site on the Co3 O4 (110)-B surface (Fig. 4b). The two oxidation steps for the two surfaces are shown in Fig. 6 and involve the removal of a proton from either an H2 O or OH ligand and removal of an electron from a Co cation. On the single-Co site, these

Page 16 of 25

reactions are H2 O−Co3+ −OH −−→ HO−Co4+ −OH + H+ + e− {4} HO−Co4+ −OH −−→ HO−Co5+ −O + H+ + e− {5} and on the dual-Co site, they are H Co3+ −O−Co3+ −OH2 −−→ Co4+ −O−Co3+ −OH2 + H+ + e− {6} Co4+ −O−Co3+ −OH2 −−→ Co4+ −O−Co4+ −OH + H+ + e− {7} The second oxidation step on the single-Co site results in a terminal cobalt-oxo species, that is known to be particularly difficult to describe using both DFT and single-configuration-based correlated wave function techniques. 51 To calculate the electrochemical oxidation free energies shown in Fig. 6, we need to calculate the energies for three states of the active site, labeled S0 , S1 , and S2 in the figure. The free energies of the proton and electron are determined using the Computational Hydrogen Electrode. 53,54 The geometries of all states and the QOs used to define the restricted subspaces are determined using the unconstrained GGA+U method (vide supra). For both active sites, each electrochemical oxidation step involves the creation of a hole in the t2g subshell of a Co atom. In the initial state S0 , the Co center on the single-Co site (site B) and both Co centers on the dual-Co site (site A) have full t2g subshells. In S1 , a single hole is present in the t2g subshell of the single-Co site and on one of the Co centers of the dual-Co site. In the fully oxidized state S2 , two holes are localized in the t2g subshell on the single-Co site; on the dualCo site, one hole is localized in the t2g subshell of each of the Co centers. The constrained subspace on the dual-Co site consists of the 6 spin orbitals in the t2g sub-


16

Page 17 of 25

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 ;^ ϬͿ

H

OH

H O

Co3+

;^ ϭͿ HO Co3+

O H

H2O Co3+

OH

H

ϭ͘ϯϯ;DͲϲͲϭϬ ůŝŐĂŶĚƐƉĂĐĞ;ɴͿ

Figure 7: Atomically-localized orbitals comprising the constrained space, ligand space, and active space in state S2 of the (top) dual-Co and (bottom) single-Co sites on Co3 O4


18

Page 19 of 25

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


the Supporting Information. One may notice that different constrained spaces, active spaces, and ligand spaces have been used for the different configurations of each active site. In principle, equivalent subspaces should be used for all states in order to rigorously calculate energy differences. This is particularly required for constructing continuous potential energy surfaces, such as reaction trajectories between two states. There are two reasons that we instead use different subspaces for the different states in this application. The first is due to the fact that the orbitals excluded from the active spaces in states S0 and S1 of both sites are either non-bonding (occupied Co 3d or O 2p orbitals) or strongly bonding/antibonding (O – H bond). These types of orbitals are well-described by the GGA, so it is unlikely that including these orbitals in the correlated subspace would lead to more accurate energies. This can be seen for the CDFT-CI method of Wu et al., 17 for which it was observed that similar energies are obtained compared to the GGA at equilibrium bond distances of H2 and H2 + . The other reason for excluding from the active space the orbitals comprising the O – H bonds in states S0 and S1 is that we have found that our method, as currently implemented, significantly underestimates the strength of strong covalent bonds, such as in H2 at the equilibrium bond length. This is very likely due to our use of the orthogonalized QOs for constructing the configuration states. At short bond distances, the orthogonalized QOs are significantly distorted away from the interatomic region. This likely increases the self-interaction error in the configuration states describing the bound system relative to those describing the unbound system (where the QOs are not distorted by the orthogonalization), and as a result, the energy of the bound system is too high. The solution to this would be to use the non-orthogonalized QOs to construct the configuration states; however, this increases the algebraic complexity of the method. This will be a topic of a future publication. The configuration states are constructed by placing the holes in active orbitals or ligand or-

bitals. In states A(S0 ) and B(S0 ), there are no holes and no active space so that only one configuration is required, in which all orbitals in the constrained space are occupied. In state A(S1 ), the single hole can reside in either the t2g orbital in the active space or in one of the six ligand orbitals. This leads to seven configurations, indicated using the following notation in which the orbital containing the hole is specified for each configuration, (1)

dyz ,

(1)

Li

i = 1, . . . , 6. The same is also true for state B(S1 ), the configurations being dyz ,

Li

i = 1, . . . , 6. In state A(S2 ), we restrict the spin-α hole to Co(1) and its ligand orbitals (Li , i = 1, . . . , 6) and restrict the spin-β hole to Co(2) and its ligand orbitals (Lj , j = 7, . . . , 12). We did not consider configurations in which both holes reside in ligand orbitals since even the configurations with one ligand hole were found to have relatively low amplitudes in the wave function expansion. The resulting 13 configurations in state A(S2 ) are (1) (2)

dyz dyz ,

(1)

dyz Lj ,

i = 1, . . . , 6

(2)

Li dyz

j = 7, . . . , 12.

In state B(S2 ), the two spin-β holes can reside in the active space consisting of two t2g and two pπ orbitals (6 configurations) or one of them can reside in one of the 10 ligand orbitals (40 configurations; we also exclude configurations in which both holes reside in ligand orbitals), dxz dyz , px py , py dyz , py dxz , px dyz , px dxz Li dyz , Li dxz , Li py , Li px i = 1, . . . , 10. In the CO-DFT calculation performed for each configuration, the +U correction is applied to the 3d electrons of the Co atoms that are


19


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

not in the active site. A U value of 5.02 eV is used for the Co atoms in octahedral sites and a value of 4.30 for those in tetrahedral sites. No +U correction is applied to the Co atom(s) in the active site. The resulting electrochemical oxidation energies are listed in Fig. 6 along with the values from unconstrained GGA+U calculation reported in our previous work. 31 We can see that the GGA+U values are surprisingly close to the MC-KS values. We also list the energies calculated using a U parameter of 3.52 eV that has been used in another study of the OER on Co3 O4 . 52 These energies are slightly closer to the MC-KS values than those using U = 5.02 eV except for the reaction B(S1 ) → B(S2 ). For this latter reaction, the energy calculated using U = 3.52 eV is 0.34 eV lower than the MCKS energy while the energy calculated using U = 5.02 eV is only 0.14 eV lower.

4

Page 20 of 25

from a CO-DFT calculation can be used to address two particular shortcomings of DFT calculations. The first concerns the transparency of DFT calculations on large systems. Atomically-localized orbitals can be used to decompose complex orbital transformations occurring during elementary chemical processes into simpler and more intuitive steps. This facilitates the identification of the electronic structure properties that control these processes, which can significantly assist with the task of optimizing them at the molecular level. We demonstrate this approach for the water addition reaction that is a key step in the oxygen evolution reaction (OER) on transition metal oxide surfaces. By doing so, we show that this reaction is controlled by the process of transferring an electron from an oxygen 2p orbital to a metal cation. A second promising application of the CO-DFT method is to perform multiconfigurational Kohn-Sham calculations to improve the accuracy of semi-local DFT for certain systems. By constructing the configuration states from an active space of atomicallylocalized orbitals, it is possible to separate intra-atomic dynamic correlation from interatomic static correlation. The former (along with intra-atomic exchange) is then treated by a semi-local exchange-correlation functional, while the latter is explicitly contained in the correlated wave function expansion. By separating the correlation energy into intra- and inter-atomic components, we minimize any double-counting. We demonstrate the multiconfigurational Kohn-Sham approach for calculating the electrochemical oxidation free energies of active sites for the OER on a cobalt oxide surface, which are poorly described by semi-local DFT.

Conclusions

We have developed a method for performing density-functional theory (DFT) calculations in which certain Kohn-Sham orbitals are constrained to be localized on single atoms. This method, which we call constrained-orbital DFT (CO-DFT), was implemented in the Vienna Abinitio Simulation Package (VASP), a widelyused planewave DFT code. The fastest convergence is achieved when different electronic structure minimization techniques are used for the constrained and unconstrained bands. The unconstrained bands are determined using the iterative eigenvalue solvers (along with charge density mixing) already included in VASP, which were modified in order to maintain orthogonality with the constrained bands. The constrained bands are determined by a direct minimization procedure based on the conjugate gradient method, whereby the wave functions and charge density are updated simultaneously. Using this technique, the computational cost for performing the electronic structure minimization is not much greater than for a standard semi-local DFT calculation. The atomically-localized orbitals resulting

Supporting Information Available See Supporting Information for details of the quasiatomic orbital method, a full derivation of the CO-DFT governing equations, further details of the electronic structure minimiza-


20

Page 21 of 25

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


tion algorithm, values of the U parameter used in GGA+U calculations, the procedure for computing inter-configurational matrix elements in the MC-KS method, further details of the MC-KS calculations on cobalt oxide, and a list of symbols used. This information is available free of charge via the Internet at http://pubs.acs.org.

(7) Hoffmann, R. A chemical and theoretical way to look at bonding on surfaces. Rev. Mod. Phys. 1988, 60, 601–628. (8) Gritsenko, O. V.; Ensing, B.; Schipper, P. R. T.; Baerends, E. J. Comparison of the Accurate Kohn-Sham Solution with the Generalized Gradient Approximations (GGAs) for the SN2 Reaction F- + CH3F to FCH3 + F-: A Qualitative Rule To Predict Success or Failure of GGAs. J. Phys. Chem. A 2000, 104, 8558–8565.

Acknowledgement C. P. acknowledges partial support from the Solar Technologies Go Hybrid initiative of the State of Bavaria, NRSCCatalysis, and the Institute for Complex Molecular Systems at the Technical University of Eindhoven. We also acknowledge NWO for access to supercomputing facilities.

(9) Gräfenstein, J.; Cremer, D. The combination of density functional theory with multi-configuration methods CAS-DFT. Chem. Phys. Lett. 2000, 316, 569–577.

References

(10) Pollet, R.; Savin, A.; Leininger, T.; Stoll, H. Combining multideterminantal wave functions with density functionals to handle near-degeneracy in atoms and molecules. J. Chem. Phys. 2002, 116, 1250–1258.

(1) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science (80-. ). 2008, 321, 792–794. (2) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289–320.

(11) Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169–11186.

(3) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Staroverov, V. N.; Tao, J. Exchange and correlation in open systems of fluctuating electron number. Phys. Rev. A 2007, 76, 040501.

(12) Dederichs, P. H.; Bl¨ ugel, S.; Zeller, R.; Akai, H. Ground states of constrained systems: Application to cerium impurities. Phys. Rev. Lett. 1984, 53, 2512–2515. (13) Wu, Q.; Van Voorhis, T. Direct optimization method to study constrained systems within density-functional theory. Phys. Rev. A 2005, 72, 24502.

(4) Mori-Sánchez, P.; Cohen, A.; Yang, W. Localization and Delocalization Errors in Density Functional Theory and Implications for Band-Gap Prediction. Phys. Rev. Lett. 2008, 100, 146401.

(14) 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.

(5) Kulik, H. J. Perspective: Treating electron over-delocalization with the DFT+U method. J. Chem. Phys. 2015, 142, 240901.

(15) Kaduk, B.; Kowalczyk, T.; Van Voorhis, T. Constrained density functional theory. Chem. Rev. 2012, 112, 321–370.

(6) Hoffmann, R. How Chemistry and Physics Meet in the Solid State. Angew. Chemie Int. Ed. English 1987, 26, 846–878.

(16) Behler, J.; Delley, B.; Reuter, K.; Scheffler, M. Nonadiabatic potential-energy ACS Paragon Plus Environment

21


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

surfaces by constrained density-functional theory. Phys. Rev. B 2007, 75, 115409.

Page 22 of 25

Analytic projection from plane-wave and PAW wavefunctions and application to chemical-bonding analysis in solids. J. Comput. Chem. 2013, 34, 2557–2567.

(17) Wu, Q.; Cheng, C. L.; Van Voorhis, T. Configuration interaction based on constrained density functional theory: A multireference method. J. Chem. Phys. 2007, 127 .

(27) Maintz, S.; Deringer, V. L.; Tchougréeff, A. L.; Dronskowski, R. LOBSTER: A tool to extract chemical bonding from plane-wave based DFT. J. Comput. Chem. 2016, 37, 1030–1035.

(18) Wu, Q.; Kaduk, B.; Van Voorhis, T. Constrained density functional theory based configuration interaction improves the prediction of reaction barrier heights. J. Chem. Phys. 2009, 130 .

(28) Kresse, G.; Furthm¨ uller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15–50.

(19) Mo, Y.; Gao, J. An Ab Initio Molecular OrbitalValence Bond (MOVB) Method for Simulating Chemical Reactions in Solution. J. Phys. Chem. A 2000, 104, 3012– 3020.

(29) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved adsorption energetics within density-functional theory using revised Perdew-BurkeErnzerhof functionals. Phys. Rev. B 1999, 59, 7413–7421.

(20) Mo, Y.; Song, L.; Lin, Y. Block-localized wavefunction (BLW) method at the density functional theory (DFT) level. J. Phys. Chem. A 2007, 111, 8291–8301.

(30) Liechtenstein, A. I.; Anisimov, V. I.; Zaanen, J. Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators. Phys. Rev. B 1995, 52, R5467–R5470.

(21) Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B 1990, 41, 7892– 7895.

(31) Plaisance, C. P.; van Santen, R. A. Structure Sensitivity of the Oxygen Evolution Reaction Catalyzed by Cobalt(II,III) Oxide. J. Am. Chem. Soc. 2015, 137, 14660– 14672.

(22) Hamann, D. R.; Schl¨ uter, M.; Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. 1979, 43, 1494–1497. (23) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953– 17979.

(32) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758–1775.

(24) Qian, X.; Li, J.; Qi, L.; Wang, C.Z.; Chan, T.-L.; Yao, Y.-X.; Ho, K.-M.; Yip, S. Quasiatomic orbitals for ab initio tight-binding analysis. Phys. Rev. B 2008, 78, 245112.

(33) Hughbanks, T.; Hoffmann, R. Chains of trans-edge-sharing molybdenum octahedra: metal-metal bonding in extended systems. J. Am. Chem. Soc. 1983, 105, 3528– 3537.

(25) Marzari, N.; Vanderbilt, D. Maximally localized generalized Wannier functions for composite energy bands. Phys. Rev. B 1997, 56, 12847–12865.

(34) Dronskowski, R.; Bloechl, P. E. Crystal orbital Hamilton populations (COHP):

(26) Maintz, S.; Deringer, V. L.; Tchougréeff, A. L.; Dronskowski, R. ACS Paragon Plus Environment

22

Page 23 of 25

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


energy-resolved visualization of chemical bonding in solids based on densityfunctional calculations. J. Phys. Chem. 1993, 97, 8617–8624.

(43) Tyo, E. C.; Yin, C.; Di Vece, M.; Qian, Q.; Kwon, G.; Lee, S.; Lee, B.; DeBartolo, J. E.; Seifert, S.; Winans, R. E.; Si, R.; Ricks, B.; Goergen, S.; Rutter, M.; Zugic, B.; Flytzani-Stephanopoulos, M.; Wang, Z. W.; Palmer, R. E.; Neurock, M.; Vajda, S. Oxidative Dehydrogenation of Cyclohexane on Cobalt Oxide (Co3O4) Nanoparticles: The Effect of Particle Size on Activity and Selectivity. ACS Catal. 2012, 2, 2409–2423.

(35) Hammer, B. Electronic factors determining the reactivity of metal surfaces. Surf. Sci. 1995, 343, 211–220. (36) van Santen, R. A.; Tranca, I.; Hensen, E. J. Theory of surface chemistry and reactivity of reducible oxides. Catal. Today 2015, 244, 63–84.

(44) Chin, Y.-H. C.; Buda, C.; Neurock, M.; Iglesia, E. Consequences of Metal-Oxide Interconversion for C-H Bond Activation during CH4 Reactions on Pd Catalysts. J. Am. Chem. Soc. 2013, 135, 15425–42.

(37) van Santen, R. A.; Tranca, I. How Molecular is the Chemisorptive Bond? Phys. Chem. Chem. Phys. 2016, 18, 20868– 20894.

(45) Suntivich, J.; Gasteiger, H. a.; Yabuuchi, N.; Nakanishi, H.; Goodenough, J. B.; Shao-Horn, Y. Design principles for oxygen-reduction activity on perovskite oxide catalysts for fuel cells and metalair batteries. Nat. Chem. 2011, 3, 546–550.

(38) Greeley, J.; Nørskov, J. K. A general scheme for the estimation of oxygen binding energies on binary transition metal surface alloys. Surf. Sci. 2005, 592, 104– 111. (39) Stamenkovic, V.; Mun, B. S.; Mayrhofer, K. J. J.; Ross, P. N.; Markovic, N. M.; Rossmeisl, J.; Greeley, J.; Nørskov, J. K. Changing the Activity of Electrocatalysts for Oxygen Reduction by Tuning the Surface Electronic Structure. 2006, 2897–2901.

(46) Dowden, D. A. Crystal and Ligand Field Models of Solid Catalysts. Catal. Rev. 1972, 5, 1–32. (47) Kremenic, G.; Nieto, J. M. L.; Tascon, J. M. D.; Tejuca, L. G. Chemisorption and catalysis on LaMO3 oxides. J. Chem. Soc. Faraday Trans. 1 1985, 81, 939.

(40) Stamenkovic, V. R.; Mun, B. S.; Arenz, M.; Mayrhofer, K. J. J.; Lucas, C. A.; Wang, G.; Ross, P. N.; Markovic, N. M. Trends in electrocatalysis on extended and nanoscale Pt-bimetallic alloy surfaces. 2007, 6 .

(48) Fuchs, M.; Niquet, Y.-M.; Gonze, X.; Burke, K. Describing static correlation in bond dissociation by Kohn-Sham density functional theory. J. Chem. Phys. 2005, 122, 094116.

(41) van Santen, R. A.; Neurock, M. Molecular Heterogeneous Catalysis: A Conceptual and Computational Approach; John Wiley & Sons Inc (E), 2006; pp 116–117.

(49) Baerends, E. J.; Gritsenko, O. V. A Quantum Chemical View of Density Functional Theory. J. Phys. Chem. A 1997, 101, 5383–5403.

(42) Redfern, P. C.; Zapol, P.; Sternberg, M.; Adiga, S. P.; Zygmunt, S. A.; Curtiss, L. A. Quantum chemical study of mechanisms for oxidative dehydrogenation of propane on vanadium oxide. J. Phys. Chem. B 2006, 110, 8363–8371.

(50) Wu, Q.; Van Voorhis, T. Extracting electron transfer coupling elements from constrained density functional theory. J. Chem. Phys. 2006, 125 .


23


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

(51) Kwapien, K.; Piccinin, S.; Fabris, S. Energetics of Water Oxidation Catalyzed by Cobalt Oxide Nanoparticles: Assessing the Accuracy of DFT and DFT+U Approaches against Coupled Cluster Methods. J. Phys. Chem. Lett. 2013, 4, 4223– 4230. (52) Garc´ıa-Mota, M.; Bajdich, M.; Viswanathan, V.; Vojvodic, A.; Bell, A. T.; Nørskov, J. K. Importance of Correlation in Determining Electrocatalytic Oxygen Evolution Activity on Cobalt Oxides. J. Phys. Chem. C 2012, 116, 21077–21082. (53) Rossmeisl, J.; Qu, Z.-W.; Zhu, H.; Kroes, G.-J.; Nørskov, J. Electrolysis of water on oxide surfaces. J. Electroanal. Chem. 2007, 607, 83–89. (54) Nørskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; Jónsson, H. Origin of the Overpotential for Oxygen Reduction at a FuelCell Cathode. J. Phys. Chem. B 2004, 108, 17886–17892.


24

Page 24 of 25

Page 25 of 25

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


Graphical TOC Entry ĂƚŽŵŝĐĂůůǇͲůŽĐĂůŝǌĞĚŽƌďŝƚĂůƐ ;KͲ&dͿ

ůŽĐĂůŝǌĂƚŝŽŶ

ĞůĞĐƚƌŽŶ

ďŽŶĚ

ƚƌĂŶƐĨĞƌ

ĨŽƌŵĂƚŝŽŶ


25

Constrained-Orbital Density Functional Theory. Computational Method

Recommend Documents