INVESTIGATION OF OXYGEN EVOLUTION REACTION CATALYSTS: ELECTRONIC PROPERTIES OF CR SUBSTITUTED COBALT (II, III) OXIDE By Guangyao Gao A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Chemistryโ€”Master of Science 2023 ABSTRACT Sustainable energy is an important topic with increasing energy demands, worldwide. Hydrogen energy is one of the promising directions in developing new energy resources. One of the promising methods for hydrogen generation is oxygen evolution reaction (OER) process. However, because of the high energy barrier associated with oxygen evolution reaction, novel catalysts are desired to promote catalytic activity. Among others, transition metal oxides have been one of the optimal materials for oxygen evolution reaction. Recent experimental studies have evaluated the catalytic performance of cobalt oxide structures. However, the fundamental details into the mechanism of OER is largely unknown, this study is focused on providing such insights of OER using first-principle calculation with density functional theory. To overcome the underestimation of the band gap, Hubbard correction (Hubbard values U = 3.3 eV for Co, 3.7 eV for Cr) is applied to the PBE functional throughout the calculation. The Hubbard correction values are validated by comparing to the experimental band gap values and the lattice parameters. The electronic properties of pristine and Cr-substituted cobalt oxides structures were investigated by computing their band structures, density of states, and projected density of states. This study clearly reports a strong correlation between the band gap and the concentration of Cr substituted in the substrate framework. The computed band gap stays relatively constant and was consistent with the experimental values of Co3O4 for Cr concentration less than 1.25. TABLE OF CONTENTS CHAPTER 1: INTRODUCTION .......................................................................................... 1 Motivation and Importance ............................................................................................... 1 Cobalt Oxide ..................................................................................................................... 5 CHAPTER 2: METHODOLOGY ....................................................................................... 10 Theoretical Method ......................................................................................................... 10 Plane-Wave Basis Set and Energy Cutoff ........................................................................ 14 Pseudopotential ............................................................................................................... 16 K-Points .......................................................................................................................... 18 Band Theory and Density of States (DOS) ...................................................................... 18 CHAPTER 3: COMPUTATIONAL DETAILS .................................................................... 23 CHAPTER 4: RESULTS AND DISCUSSIONS .................................................................. 25 Energy Cutoff ................................................................................................................. 25 K-Points Optimization .................................................................................................... 26 Convergence Threshold ................................................................................................... 26 Substitution Occupation Analysis .................................................................................... 27 Optimalization of Chromium Substituted Structures ....................................................... 28 Electronic Properties of Conventional Structures ............................................................ 30 Band Gap Calculation with Hubbard Correction & Band Structure ................................. 30 Band Structure & Density of States ................................................................................. 32 CHAPTER 5: ONGOING CALCULATIONS ..................................................................... 41 Super Cell with OER Intermediates ................................................................................ 41 CHAPTER 6: FUTURE WORK ......................................................................................... 43 REFERENCES ................................................................................................................... 44 iii CHAPTER 1: INTRODUCTION Motivation and Importance Electricity has become the major power source in the late nineteenth century. The consumption of electricity worldwide has increased to 23,398 million kWh per year.1 With this tremendous energy demand combined with depletion of fossil fuels and the environmental impact of these unrenewable sources,2 the importance of clean secondary energy is paramount. Consequently, hydrogen energy, which uses hydrogen as a fuel to generate energy either by combustion or fuel cells is a highly promising renewable energy resource and has the potential to be widely used in the future. Unlike other alternate sources, hydrogen energy is considered to be a clean secondary energy resource because of its high burning velocity and environmentally friendly combustion product (water). Moreover, hydrogen can also be used as a reactant in industrial process such as the Fischer-Tropsch process and ammonia synthesis.3 There are many methods employed in industry to generate hydrogen, including water electrolysis and natural gas reforming.4 Electrolytic cells are also used for hydrogen production using different kinds of oxides such as methanol, ethanol, urea in industry.5 Electrolysis of water is a promising future direction to generate hydrogen. Water electrolysis is a process that will decompose water into oxygen and hydrogen gas using external potential. Water electrolysis can be divided into two half-reactions an oxygen evolution reaction (OER) and a hydrogen evolution reaction (HER).6 A schematic representation of an electrolytic cell for water electrolysis is represented in Figure 1.6 The electrolysis process requires an external electrical energy to drive the nonspontaneous reactions. The OER reactions occur at the anode where the water molecules are oxidized to form oxygen molecules and protons. The protons pass through the membrane and 1 combine with electrons to form hydrogen gas at the cathode and hence termed as HER. One of the difficulties associated with the water electrolysis process is the high energy barrier that must be overcome to initiate OER and generate molecular oxygen Figure 1: Schematic representation of electrocatalytic water splitting. OER and HER occur under different pH conditions. The overall reaction involved in water splitting can be written as: ๐‘‰๐‘œ๐‘™๐‘ก๐‘Ž๐‘”๐‘’ 2๐ป2 ๐‘‚(๐‘™) โ†’ ๐‘‚2 (๐‘”) + 4๐ป + + 4๐‘’ โˆ’ (1) This is an endothermic reaction. Under ideal thermodynamic conditions, there are four electrochemical steps (discussed in the methodology section) involved in water electrolysis and each step will require a difference in Gibbs free energy of 1.23 eV under equilibrium condition (Temperature (T) = 298.15 K, Pressure (P) = 1 bar, pH = 0). However, this value is always higher under practical conditions and can vary significantly depending on the electrode surface where the reactions occur. The additional potential greater than 1.23 eV required to initiate the reaction is termed as overpotential. (A detailed discussion of the overpotential is presented in the 2 methodology section). In short, the overpotential is the energy difference between the rate determining step reaction and the thermodynamic equilibrium potential for this step in an ideal condition, which is 1.23 eV. It can be summarized as the equation: ๐œ‚ = ๐ธ๐‘š๐‘Ž๐‘ฅ โˆ’ 1.23๐‘’๐‘‰ (2) Thus, the lower the overpotential, the higher the efficiency of a reaction due to its lower energy cost. Subsequently, catalysts are used to reduce the overpotential value as much as possible by overcoming energy barriers as evidenced in various chemical reactions. There are two kinds of inorganic catalysts namely, homogeneous catalysts and heterogeneous catalysts examined for OER. Homogeneous catalysts are catalysts that are in the same physical phase as the reactants while heterogeneous catalysts are in a different physical phase than the reactants. One advantage for heterogeneous catalysis is that it is easy to separate the catalysts from the reactants and products, which typically results in an easier recycling process as compared to that for homogeneous catalysts and reduces the operational cost substantially.7 Hence, heterogeneous catalysts are more desirable as compared to homogeneous ones in OER. A number of heterogeneous catalysts examined in OER studies are based on novel metal elements, transition metals and their oxides. These systems are considered as ideal catalysts based on Sabatierโ€™s principle, which states that the binding energy of an ideal catalyst with a reactant should be neither too strong nor too weak.8 Importantly, novel metals such as Iridium, Ruthenium and Platinum have shown promising catalytic activity during OER studies.9 A recent density functional theory (DFT) computational study has shown that the rutile IrO2 has a theoretical overpotential of 0.38 V in an OER cycle, which is promising.10 However, the limited reserves of 3 those noble metal and their oxides inhibit the use of large scale catalytic application. At the same time, the use of transition metals and their metal oxides have shown highly promising behavior during OER by yielding lower overpotentials. Morales et al.11 has performed DFT calculations based on monoxides of various transition metals showed different catalytic activity with different metal oxides. Figure 2 illustrates a volcano diagram showing the OER catalytic activity for those monoxides. Figure 2: Volcano plot for transition metal monoxides for (001) surface. Following Sabatierโ€™s principle, the catalysts in the middle of the volcano plot (Figure 2) show the binding energy between the reactants and catalysts indicates the optimal catalytic activity. This also agrees with the theory of the theoretical overpotential, the smaller the overpotential the better OER catalysts efficiency. In Figure 2, oxides of Ni, Mn, Co and Fe all represent a strong potential to be the OER catalysts. Besides pure metal oxides, studies have indicated that foreign metal doping in transition metal oxides will further enhance the catalytic activity in OER. Sun et al.12 demonstrated the 4 influence of foreign elements in altering the catalytic activity of transition metal oxide using copper (Cu) doped IrO2 surface. The increased catalytic activity could be attributed to 2 important factors: (i) an oxygen vacant site is formed due to the substitution of Cu2+ for Ir4+; (ii) distortion of IrO2 lattice due to the presence of CuO6 octahedra. Consequently, the electronic configuration in t2g and eg orbitals changes significantly. Similarly, experimental studies from Lin and McCrory examined the influence of Cr substitution on the catalytic behavior of Co3O416. Despite the experimental catalytic properties, the mechanism for this catalyst in OER is still unknown. To have a better understanding of the Cr substituted Co3O4 and OER mechanism, this study will deeply investigate the importance of Cr substituted system and their electronic properties. Cobalt Oxide Among the transition metals mentioned above, cobalt is a common, earth abundant material. Due to the transition metal properties, cobalt has been widely used in many different industrial applications including but not limited to batteries and paints. The electronic properties of the cobalt oxide play a significant role in rechargeable batteries. For instance, lithium doped cobalt oxides are considered as a good positive electrode material because of the high structural reversibility and thermal stability.13 In recent years, increasing number of studies are focusing on the efficacy of cobalt oxides during electrocatalysis14,43. For instance, the catalytic efficiency of Co3O4 is greatly enhanced with incorporation of dopants in their structure as demonstrated by Banerjee et al.14 with three different foreign metals: nickel (Ni), chromium (Cr), and iron (Fe). as shown in Figure 3(a). In particular, the Cr dopant has the strongest positive effect on catalytic behavior of cobalt oxide because it requires lowest potential at the same current density. The higher efficiency of Co3O4 is evident on Tafel slope which is another important parameter to evaluate the catalytic behavior and is given by the following equation: 5 ๐‘– ฮท = A ร— ๐‘™๐‘œ๐‘”10(๐‘– ) (3) 0 where ฮท is the overpotential, A is the Tafel slope and ๐‘– is the current density. In general, a catalyst with a smaller Tafel slope results in an efficient catalytic property. The Cr substituted cobalt oxide shows a slope of 57.7 mV/dec, which is lower than any other catalysts in their study as shown in Figure 3(b). Figure 3: (a) Current densities for different cobalt oxide dopants. (b) Tafel slope for different cobalt oxide dopants. Moreover, recent DFT studies has shown that iron and nickel doped cobalt oxides have different overpotentials on different active metal sites.15 In their study, the Co3O4 surface were cleaved with two different terminations namely A and B. In a solid structure, the surface exposed can be different due to the atom arrangements. In cobalt oxide, the A and B layers are arranged alternately. Different atom arrangements will show different properties Both terminated surfaces were examined because different terminations have different active metal sites. Their calculation 6 clearly highlights that the Ni doped Co3O4 showed a lower overpotential comparing to pure Co3O4 for the octahedral cobalt in termination B (Table 1). Table 1: Overpotential for different active metal sites on different terminations. Catalyst (B-layer) Overpotential (V) Catalyst (A-layer) Overpotential (V) Co3O4 (CoOct) 0.46 Co3O4 (CoOct) 0.63 FexCo3-xO4 (CoOct) 0.45 FexCo3-xO4 (CoOct) 0.37 NixCo3-xO4 (CoOct) 0.34 NixCo3-xO4 (CoOct) 0.41 FexCo3-xO4 (Fe) 0.76 FexCo3-xO4 (Ni) 0.77 Importantly, recent experimental studies by Lin et al.16 demonstrated the influence of dopant concentration on the electrocatalytic properties of Co3O4. Chromium was used as the dopant metal in their study. They investigated a series of chromium doped cobalt oxide Co3-xCrxO4 and their associated catalytic activity of OER, where โ€˜xโ€™ represents the amount of chromium was doped into the spinel system. Their study clearly demonstrated that the catalytic performance greatly increases with the increasing amount of Cr doped in the cobalt oxides (when x โ‰ค 0.75).16 The lattice parameters (shown in Table 2) tend to increase with increase in the amount of chromium dopants. A high catalytic performance was reported when x = 0.75 (see Table 2). The high OER activity characteristics is primarily due to the lowest value for the overpotential. However, the fundamental insights are not completely known. For example, what factors dictate the catalytic activity increases at low Cr concentration while it decreases at high concentrations when x is 7 greater than 0.75. Besides, what drives the change in the electroactivity behavior if chromium occupied to a Co2+ site instead of Co3+ site, which is the primary basis of this study. Table 2: Experimental results for different chromium doped catalysts lattice parameters and overpotentials. Catalyst Lattice Parameters (โ„ซ) Overpotential (V) Co3O4 8.081 0.42 ยฑ 0.01 Co2.75Cr0.25O4 8.106 0.38 ยฑ 0.01 Co2.5Cr0.5O4 8.124 0.36 ยฑ 0.01 Co2.25 Cr0.75O4 8.135 0.35 ยฑ 0.01 Co2CrO4 8.203 0.37 ยฑ 0.01 Co1.75Cr1.25O4 8.209 0.38 ยฑ 0.01 Co1.5Cr1.5O4 8.233 0.39 ยฑ 0.01 Co1.25Cr1.75O4 8.267 0.38 ยฑ 0.01 CoCr2O4 8.285 0.40 ยฑ 0.01 In this project, the goal is to set up a detailed understanding about the conventional structure for chromium substituted structures, especially focusing on their electronic properties. There have been studies about nickel substituted cobalt oxide15, However, there havenโ€™t been theoretical studies about using chromium as the foreign metal to create the novel substituted structures. The cobalt in pristine Co3O4 structures have 2 different coordination sites, namely tetrahedral Co2+ and octahedral Co3+ sites. Therefore, it is important to determine the optimal crystallographic sites for chromium substitution. The Hubbard value in PBE+U method will also be tested with different U values. 8 Importantly, the current project will be focused on the electronic properties, which includes both band structures and the corresponding density of states. The electronic properties play an important role in solid state calculation which shows the conductivity of the material. Such insights into the electronic properties of the material will be valuable to understand the bonding environment and also lay a strong basis for predicting mechanism in a catalytic rection. The project will address the relationship between the calculated band gaps and the concentration of Cr in Co3O4. We also highlight the correlation of density of states with band gap and how it varies with increasing Cr concentration. Furthermore, understanding the electronic properties for the conventional structure will be constructive for the future direction. To investigate the OER catalysis on the metal surface, local density of states will be performed. These calculations will help experimentalists to understand the bonding properties with the small water molecule intermediates and potentially explain the mechanism of the OER process which is one of the major subjects for this collaboration. Therefore, understanding the energy states distribution from each atomic orbitals in the conventional structure will be a strong support and reference to compare with the local density of states. 9 CHAPTER 2: METHODOLOGY Theoretical Method This study will address the fundamental characteristics dictating the electrocatalytic behavior of Co3O4 using density functional theory. The principal idea of quantum mechanics is to solve the time independent Schrรถdinger equation and is given below: ฬ‚๐œ“ = ๐ธ๐œ“ ๐ป (4) where H represents the Hamiltonian, E represents eigenvalues of this equation, respectively.17 The wave functions is a set of functions to describe the quantum state for each particle in the space. According to Schrรถdinger equation, each particle has its own motion so they should be considered as different individual units. This is the many-body Schrรถdinger equation. One of the biggest concerns for many-body system is that every particle will have its own coordinates. Taking cobalt oxide as an example, the total number of electrons in the system is 113, resulting in a total of 339 coordinate functions and their system property is hard to obtain by solving the Schrรถdinger equation. On the other hand, density functional theory (DFT) considers electron density at any given position in space to calculate the ground state energy of the system rather than the many-body wave function problem. Thus, the ground state energy indicates the atoms all occupy the lowest energy level. Following with this idea, the DFT is based on 2 Hohenberg and Kohn theories:18 (1) The ground state energy of the system described by the Schrรถdinger equation is a unique functional of the electron density. 10 (2) The electron density that minimizes the total energy of the system is the exact ground state density of the system described by the Schrรถdinger equation. Since all the electrons are in one 3-dimensional space, we can reduce the many body coordinates into a 3-dimensional coordinate system. The Kohn and Sham equation in terms of electrons density can be written as follows: โ„2 2 (5) [โˆ’ ๐›ป + ๐‘‰(๐’“) + ๐‘‰๐ป (๐’“) + ๐‘‰๐‘‹๐ถ (๐’“) ] ๐›น๐‘– (๐’“) = โ„‡๐‘– ๐›น๐‘– (๐’“) 2๐‘š The terms in order from the left side are kinetic energy of the electrons, the coulombic nuclei-electron attraction ๐‘‰(๐’“) , electron-electron repulsion๐‘‰๐ป (๐’“) , and an unknown exchange- correlation potential ๐‘‰๐‘‹๐ถ (๐’“). The first 3 terms in eq.5 are well-defined. In addition, it is important to note that in the Hartree potential, considering only 1 electron, the electron interacts with their electric field/density generated by all electrons. However, in this approximation, the electric field contributed from all the electrons including the electron that is being observed. This is called self- interaction error and is corrected using a term namely exchange-correlation term. However, the exchange-correlation potential is still a complex functional derivative of the exchange-correlation energy. This term defines all the unknown effects in this equation that is not included from the single electron approximation. However, the true functional form of the exchange-correlation is not known. Consequently, there are different approximations used to define this functional. The 2 main approximations commonly employed are Local (Spin) Density Approximation (LDA or LSD)18 and Generalized Gradient Approximation (GGA).19 The simplest functional is LDA which considers the system filled with uniform electron gas. This results in a constant electron density at 11 any given point within system space. Under this condition, the density within the entire volume is constant. If the density is a constant, then the local density around the particle of interest is considered. The equation can be written as: ๐ฟ๐ท๐ด [๐‘›(๐’“)] ๐ธ๐‘‹๐ถ = โˆซ ๐‘›(๐’“) ๐œ€๐‘ฅ๐‘ [๐‘›(๐’“)] ๐‘‘ 3 ๐’“ (6) The exchange-correlation energy can be divided into 2 terms: ๐ธ๐‘‹๐ถ = ๐ธ๐‘‹ + ๐ธ๐ถ (7) Where ๐ธ๐‘‹ and ๐ธ๐ถ represents exchange and correlation term. The exchange term can be solved by using the idea of uniform electron gas. The energy density at a certain position will be the same as the energy density of the uniform electron gas. It is obvious to see that the LDA works fine for slowly varying densities.20 However, for the system where the electron densities change fast, LDA tend to fail to predict the values comparing to experimental values. Then it is appropriate to introduce a functional that includes an adaptive density distribution: ๐บ๐บ๐ด [๐‘›(๐’“)] ๐ธ๐‘‹๐ถ = โˆซ ๐‘“(๐‘›(๐’“), โˆ‡๐‘›(๐’“)) ๐‘‘3 ๐’“ (8) This is the generalized gradient approximation (GGA). The GGA functionals are based on the local electron density and the gradient of the electron density. Two main functional namely Perdew- Wang (PW91) and PBE are widely used in GGA. Comparing to PW91, PBE is a more simple GGA 12 functional developed by John P. Perdew, Kieron Burke and Matthias Ernzerhof.21 Although GGA functional represents a more sophisticated description of electron density compared to LDA, this does not guarantee that GGA is always better than LDA for solid-state calculations. These functionals are all widely used in different field in chemistry. In addition, there are many other functionals that have been developed to improve the accuracy of approximation of the exchange- correlation functional such as meta-GGA22 and hybrid-functionals. However, DFT have its limitations on underestimating the electron-electron interactions. When examining transition metals, because of the strong correlation effect of the 3d electrons, DFT cannot provide an accurate result. The main reason of this underestimation is because of the electron-electron interactions. To overcome this, the Hubbard model has been used in this situation. The Hubbard Hamiltonian can be written as:23 โ€  ๐ป = โˆ’๐‘ก โˆ‘ โˆ‘ ๐ถ๐‘–๐œŽ ๐ถ๐‘—๐œŽ + ๐‘ˆ โˆ‘ ๐‘›๐‘–๐œŽ ๐‘›๐‘–๐œŽโˆ— (9) โŸจ๐‘–,๐‘—โŸฉ ๐œŽ ๐‘– Where t stands for the energy needed to transfer the electron to its neighbor orbitals. โŸจ๐‘–, ๐‘—โŸฉ stands โ€  for the 2 nearest neighbor orbitals. And the ๐ถ๐‘–๐œŽ , ๐ถ๐‘—๐œŽ are the operators that create and eliminate a certain electron in state โ€˜iโ€™ without changing the spin. ๏ณ and ๏ณ* stands for the opposite electron โ€  spins. The โ€˜nโ€™ operator is the counting operators which satisfies ๐‘›๐‘–๐œŽ = ๐ถ๐‘–๐œŽ ๐ถ๐‘–๐œŽ . Because of the strong repulsion between the โ€˜dโ€™ orbital electrons, the DFT result will give a shorter bandwidth. This is the reason why sometimes DFT fails to predict the conductivity of a transition metal material. The solution for the failure of DFT is the Hubbard model. An additional term that describes the โ€˜dโ€™ electrons repulsion can be added to overcome this insufficient description. This 13 is the DFT + U. In a more specific way, this method can be written as LDA + U or GGA + U depends on the functional has been used to approximate the exchange-correlation term. Plane-Wave Basis Set and Energy Cutoff The KS orbital can be represented using basis set and can be written with a linear combination form: ๐พ ๐œ“๐‘— (๐‘ฅ) = โˆ‘ ๐›ผ๐‘—,๐‘– ๐œ‘๐‘– (๐‘ฅ) (10) ๐‘–=1 where ๐›ผ is the expansion coefficients, ๐‘– stands for the number of electrons, ๐‘— is number of the number of spin orbitals. This set of functions is called basis set. It is easy to understand that choosing a larger basis set will have a more accurate calculation result. However, this will also increase the computational cost. There are many ways to expand the basis set, for example, the local basis set: Slater type orbitals (STO),24 Gaussian type orbitals (GTO)25, etc. At the same time, plane wave basis set is widely used for periodic system and provides leverage over basis set superposition error. Based on Blochโ€™s theorem, the wave function of a periodic system can be written as:26 ๐œ“๐‘› (๐‘Ÿ) = ๐œ‡๐‘› (๐’“)๐‘’ ๐‘–๐’Œโ‹…๐’“ (11) where ๐‘˜ is the wave vector in first Brillouin zone or reciprocal space, ๐‘Ÿ is the wave vector in the real space. To obtain the reciprocal space, we can apply the Fourier transform to the real space. Reciprocal vectors ๐’ƒ๐’‹ are defined and satisfies ๐’‚๐’Š โˆ™ ๐’ƒ๐’‹ = 2๐œ‹ when ๐‘– = ๐‘—. The ๐’‚๐’Š are vectors in real 14 space. The space that consists of all the reciprocal vectors is called the reciprocal space. In the reciprocal space, the volume formed by the planes that are perpendicular to reciprocal vectors are called Brillouin zone. The space which has the smallest distance to a center point of a lattice is called the first Brillouin zone. In solids state chemistry that focused on electrocatalytic activity, all calculations are performed in reciprocal space. The reason to use Blochโ€™s theorem is because this theorem reduces the area that needs to be calculated. The ๐œ‡๐‘› is the periodic function that have the same periodicity as the lattice. Because of this, we do not have to calculate for the entire lattice but only calculate the smallest repeating unit of the lattice, which is the Brillouin zone. Also, the bigger the lattice in the real space will have a smaller lattice in the reciprocal space. This is the main reason plane-wave basis set is best suited for calculation on solid systems. The periodicity of ๐œ‡๐‘› (๐’“) can be expanded as: ๐œ‡๐‘› (๐’“) = โˆ‘ ๐‘๐‘ฎ ๐‘’ ๐‘–๐‘ฎโˆ™๐’“ (12) ๐‘ฎ where G is a set of vectors in reciprocal space and satisfies ๐‘ฎ โˆ™ ๐’‚๐’Š = 2๐œ‹๐‘š๐‘– for any given real space vector ๐’‚๐’Š . Combine the 2 equations above can get: ๐œ“๐‘› (๐‘Ÿ) = โˆ‘ ๐‘๐‘ฎ+๐’Œ ๐‘’ ๐‘–(๐’Œ+๐‘ฎ)โ‹…๐’“ (13) ๐‘ฎ However, in this situation, the ๐‘ฎ has infinite terms to values. To make this wavefunction possible to solve, we can choose an energy cutoff for the kinetic energy and make sure: 15 โ„2 2 โ„2 2 (14) ๐ธ= |๐‘˜ + ๐บ| < ๐ธ๐‘๐‘ข๐‘ก = ๐บ 2๐‘š 2๐‘š ๐‘๐‘ข๐‘ก In this case, a finite number of plane-waves are quired. This is the reason why energy cutoff must be defined in plane-wave method. Pseudopotential A schematic representation of the pseudopotential can be viewed in Figure 4. It clearly emphasizes that when the electrons are closer to the nuclei, the wavefunction will be more complicated to calculate because they are highly oscillatory.27 A normal plane wave basis set will require a large number of plane waves to describe the system. This is primarily due to the contribution of the core electrons to the chemical or physical property is small. In this case, the true core electrons wave functions and their potential becomes less important. For this reason, a pseudopotential or pseudo wave function can be used to smoothen the electron density around the core region of an electron. The core electron and nuclei are termed as a pseudo core, hence referred as frozen core approximation. 16 Figure 4: Comparison between pseudo wavefunction (blue) and pseudopotential and full potential (red). Another advantage about this approximation is it reduces the total number of plane waves needed in a calculation as it focuses on the valence electrons.28 There are different kinds of pseudopotentials. The most basic one, norm-conserving pseudopotential (NCPP).29 The core idea for NCPP is that the total charge density stay the same comparing to all electronsโ€™ situations when the radius is greater than the cutoff radius. NCPP can give accurate results in many aspects, especially dealing with scattering patterns.29 However, NCPP need higher cut-off energy which can be computational demanding. Consequently, ultrasoft pseudopotentials (USPP)30 have been developed which represent a smoother wavefunction than the NCPP. Therefore, the basis set will 17 be smaller. Moreover, the cut-off energy will be smaller than the norm-conserving pseudopotentials. K-Points As discussed above, the Brillouin zone is used for periodic system. However, integrating the first Brillouin still needs lots of computational resources. In order to simplify the calculation, k-points are used. k-points stands for the high symmetry points selected in the reciprocal space. The k-points are commonly generated using Monkhorst-Pack method.31 For example, a 4 ๏‚ด 4 ๏‚ด 4 k-points grid indicates there are 4 k-points in each direction. Instead of calculating the entire Brillouin zone, only the selected k-points need to be calculated. It is also obvious to conclude that more k-points will better represent the system in the Brillouin zone however, this will increase the computational cost. Because of this, it is necessary to determine the appropriate k-point mesh for a system. Band Theory and Density of States (DOS) To understand the electronic property of a solid-state structure, one of the important methods is to analyze density of states. To better understand density of states, the band theory will be introduced. To interpretate the band theory as a chemist and starting from the molecular theory, 2 hydrogen atoms will form 2 bonds, bonding, and anti-bonding. Imagine that hydrogen atoms will from a hydrogen atom chain, when the โ€œmoleculeโ€ grows bigger, the number of molecular orbitals will also increase. When the molecule is big enough and still to form crystal, the orbitals are dense enough to form a band. 18 Figure 5: Illustrations of band theory from chemistry perspective. White and shaded circles represent imaginary hydrogen atoms forming a chain. In MO theory, molecular orbitals are formed by a linear combination of atomic orbitals. In crystal structure, the orbitals followed the same pattern. As discussed above, the periodic crystal orbitals can be described by Bloch function, we can rewrite equation 11 as: โˆž ๏น(x) = โˆ‘ ๐‘’ ๐‘–๐‘˜๐‘ฅ ๏ฃ(x) (15) ๐‘›=1 Where x is the position to the origin. Since the system is a crystal structure, x can only be within the boundary a, which a is the lattice constant or distance between two neighbor atoms. Because of the periodicity of the crystal orbitals, the equation can be written as: 19 โˆž ๏น(x + Na) = โˆ‘ ๐‘’ ๐‘–๐‘˜(๐‘ฅ+๐‘๐‘Ž) ๏ฃ(x + Na) (16) ๐‘›=1 N is the Nth atom on the infinite lattice chain. The exponential portion is expanded with the Eulerโ€™s equation: ๐‘›(2๐œ‹) ๐‘˜= (17) ๐‘๐‘Ž Since n can be either positive or negative value, the maximum of n will be half of the N. Therefore, ๐œ‹ ๐œ‹ k can only choose from the range between - ๐‘Ž to ๐‘Ž . This is the range of the first Brillouin zone. To visualize the energy state at different k value, band structure diagram was introduced. Even through the k value is quantized, however, because of the number of atomic orbitals are big enough, the band diagram can be considered as continuous49 as shown in the figure below: Figure 6: Illustration of band structure diagram and corresponding density of states. 20 when k = 0 is not always the lowest energy. This depends on the different type of atomic orbitals. The density of states (DOS) is closely related to the band structure. Density of states describe the number of energy states per unit volume over the energy range E and E + ๏คE. As the diagram shows, when the band structure is flat from the left, it indicates more states in a small energy change. This phenomenon is shown on the diagram on the right when the flat portion of the band structure will experience a higher peak on the density of states diagram. Density of states can be reviewed in a more detailed way. The projected density of states (PDOS) can project the energy states onto each of the atomic orbitals. One of the most important features for PDOS is to show the energy states contribution from different atomic orbitals, then potentially will show the bonding property information. For example, Mina Yoon et al. showed the bonding properties of MC60 (M = Be, Mg, Ca, Sr) complexes using PDOS50. Figure 7: PDOS of different foreign metal with the interaction of C60. As shown in Figure 7, the black peaks are the total PDOS from C60. It is obvious that the first 2 element shown in the figure, the orbitals from the metal atoms have no major interactions 21 with the carbon orbitals. This indicates there is no chemical bond formation between the carbon and the Be and Mg metal. On the other hand, the energy states from Ca and Sr s orbitals and p orbitals interact with the carbon orbitals from C60. This indicates a strong hybridization between the carbon orbitals and the metals orbitals. Moreover, the peak in Ca and Sr is shorter than the Be and Mg s orbitals. This depletion indicates the electrons from the metals transferred to the C60. 22 CHAPTER 3: COMPUTATIONAL DETAILS The structure of Co3O4 used in this current study was obtained from the material project database (ID: mp-18748). The Co3O4 belongs to the spinel structure like iron oxide35 where there is 2:1 distribution of octahedral (16) and tetrahedral (8) of cobalt sites with different oxidation state namely Co3+ and Co2+, respectively. The octahedral sites are filled by Co3+ while the tetrahedral sites are filled by Co2+. A schematic representation of the conventional unit cell of Co3O4 is shown in Figure 8. The 3 unpaired electrons in the d-orbital of Co2+ will half fill the t2g level with higher energy than eg while Co3+ will fully fill all the t2g orbitals and the eg at higher energy level are unfilled. The cobalt oxide behaves as a semiconductor under room temperature. Cobalt oxide is a p-type semiconductor with band gap value 1.52 to 2.13 eV.36 The Co3O4 structure represents face- centered cubic formed with oxygen having a lattice parameter of 8.08 ร… determined experimentally.37 The distance between the Co3+ and O2- is 1.89 ร… while between Co2+ and O2- is 1.99 ร….37 All the KS-DFT calculation were performed with Quantum Espresso package. The common GGA functionals namely the Perdew-Burke-Ernzerhof (PBE) were used to approximate the exchange-correlation functional. Ultrasoft pseudopotential were used through the calculation. The structure was geometry optimized to obtain the optimal values for energy cutoff, convergence threshold and k-points were calculated. Those parameters were evaluated to calculate with the efficient usage computational resources. All the doped conventional structures were relaxed to optimized structures including optimized atomic positions and cell parameters. 23 Oxygen Cobalt Figure 8: Conventional Co3O4 structure, red represents oxygen, blue represents Co2+ and Co3+. In order to determine the energetically stable doping sites in the Co3O4 lattice, a series of calculations were performed where the ratio of Cr doping varies from complete octahedral to tetrahedral substituted sites. After identifying the favorable dopants sites, a series of Cr doped Co3O4 structure was created with the structural formula of Co3-xCrxO4. To investigate the effect of Cr dopant concentration on the catalytic properties, the value of โ€˜xโ€™ varies from 0 to 2 with steps of 0.25 and their lattice parameters were obtained by allowing to change the cell dimension. To better understand electronic properties of the structure, Co3O4 band structure and corresponding density of state were determined. To overcome the underestimation of the band gap when calculating transition metal materials, there are variety studies have shown the importance of GGA + U method15, 43. The PBE+U (U = 3.3 eV for cobalt element, 3.7 eV for chromium element) calculation were performed to overcome the underestimation for DFT dealing with strong correlation system. The value of cobalt and chromium were obtained semi-empirically from literature46 and validated by comparing the cell parameters with the experimental value after the Hubbard correction. It should be noted that for all the calculations, the Fermi energies were as reference and kept at 0 eV. 24 CHAPTER 4: RESULTS AND DISCUSSIONS When calculating a novel structure in solid states, 3 parameters need to be tested: Energy cutoff, convergence threshold and k-points. It is important to test these parameters at the early stage to obtain an accurate result which is important to use the computational resources appropriately. Energy Cutoff Figure 9 shows the variation in the total energy of the system with increasing energy cutoff. Figure 9: (a) Variation of total energy as functions of energy cutoff for Co3O4 (b) Computational time for Co3O4 conventional structure when changing the energy cutoff from 50 Ry to 120 Ry (20 cores). As shown in Figure 9, the total energy of the Co3O4 decreases with increase in energy cutoff and reaches a shallow energy minimum at ~80 Ry. With further increase in the energy cutoff, the change in the total energy does not vary significantly. Notice that the total energies are relatively big in their absolute values. This indicates that the increase in energy cutoff does not impact the total energy of the system calculation result. Moreover, the total energy change per unit atom in this conventional structure between each calculation is smaller than 0.001eV/atom. The 25 inconsistency of the reported computational time could be attributed to the calculations being run on different nodes. K-Points Optimization Figure 10 clearly indicates that the increase in k-points beyond k = 4 have negligible impact on the total energy of the Co3O4 system. In addition, the computational time cost for 4x4x4 grid is 105 minutes which is substantially lower than larger number of k-points. Both results show that k = 4 should be appropriate size for studying the structural properties. It is also significant to point out that changing the cut off energy will not change the choice of the k-points. Figure 10: (a) different sets of k-points in different cut-off energy. (b) Computational time with different sets of k-points (20 cores). Convergence Threshold The convergence threshold parameter defines how accurate the self-consistent calculation result can be considered as converged. The calculations shows that the convergence threshold is 1.0๏‚ด10-8 is sufficient enough to address catalytic performance. The reason is the total energy reaches the lowest value for this variable. Moreover, as illustrated in Figure 11 b), the computational time will increase when decreasing the threshold value. 26 Figure 11: (a) Total energy with respect to convergence threshold. (b) Computational time with respect to convergence threshold. Substitution Occupation Analysis As mentioned earlier, there are two types of cobalt structures, Co3+ and Co2+ present in this spinel structure, the energetic favorable occupation sites need to be tested to verify the structures with foreign substitution have the lowest energy. Different types of doping will result in a different total energy of the system. A series of calculations were performed to determine the ideal location in the structure for dopant substitution. When x = 0.75, the optimizations were performed with chromium dopants at different sites: Firstly, chromium atoms were exclusively distributed in tetrahedral sites, later the distribution of tetrahedral chromium sites will be gradually decrease until all chromium atoms were doped in octahedral sites. It is evident from Figure 9 that the chromium dopants structures show high stability only when occupying the octahedral sites, based on the total energy profiles. Besides, evident from figure 12 b), the lattice parameters reach the lowest value with all octahedral occupation. 27 Figure 12: (a) The number of octahedral sites with respect to total energy. (b) Lattice parameters. Optimalization of Chromium Substituted Structures Based on the lattice parameters and occupation results, a set of geometry relaxation calculations with different amounts of chromium dopants were performed to obtain the optimized structures. The lattice parameters obtained of the optimized structures were compared with experimental values in Table 3. It is evident from Figure 13 that the magnitude of the lattice parameter increases with increase in Cr concentration both consistent with reported experimental data and computational results. Table 3: Lattice parameters with different fractions of chromium for calculation results (left) and experimental results (right). Catalyst Comp. Lattice Parameters (โ„ซ) Exp. Lattice Parameters (โ„ซ) Co3O4 8.164 8.081 Co2.75Cr0.25O4 8.201 8.106 Co2.5Cr0.5O4 8.242 8.124 Co2.25 Cr0.75O4 8.285 8.135 28 Table 3 (contโ€™d) Co2CrO4 8.321 8.203 Co1.75Cr1.25O4 8.369 8.209 Co1.5Cr1.5O4 8.409 8.233 Co1.25Cr1.75O4 8.496 8.267 CoCr2O4 8.518 8.285 The reported lattice parameters are larger than the experimental values, however, it is acceptable for the computational value is slightly different than the experimental values in a certain limit48. Moreover, the differences in lattice parameters between the two consecutive concentrations is consistent with different amount of chromium substitution (about 1% to 2%). This indicates the calculation results are acceptable comparing to experimental values, which is reasonable to validation of the calculations. Figure 13: (a) lattice parameters with magnetization & Hubbard correction (b) lattice parameters from experimental results. 29 At the same time Figure 14 clearly advocates the need to use Hubbard correction. The reference calculation performed without magnetization and Hubbard correction are completely not in correlation with experiments. this clearly emphasizes the importance of employing magnetization in our calculation. Figure 14: Lattice parameters without magnetization and Hubbard correction. Electronic Properties of Conventional Structures The electronic band properties also play a significant role in understanding the behavior of a catalyst. The two parameters that is widely used to determine the electronic properties of the materials is band gap and density of states. The band structure is strongly related to both density of states and the band gap. For instance, the band structure will show the band gap to determine whether the gap is a direct or indirect band gap. As for density of states, it showed how number of states change with the change of energy instead of the positions in reciprocal space. Band Gap Calculation with Hubbard Correction & Band Structure When calculating the band structure for a transition metal system, it is important to overcome the underestimation tendencies of the band gap with transition metals by using the Hubbard 30 correction23. However, to find the values of the Hubbard correction, band structure calculations need to be performed. The Hubbard correction will be considered in this step and the calculated band gap will be compared to results from calculations using the PBE functionals without the Hubbard correction. Theoretically, the Hubbard correction value was calculated by using the linear response approach.40 However, identifying Hubbard values in empirically is more common.41,45 To start the calculations, Co3+ and Co2+ were set values at 6.7 eV and 4.4 eV, respectively. To simplify the calculation, the Hubbard correction was averaged out by fraction to be 5.9 eV from the literature.41 However, the band structure showed that the Hubbard value from the literature was too big in the current conventional system, in a result of band gap about 2.5 eV (about 1.6 eV for experimental value47). To fit the band gap in the range of experimental values, different Hubbard values were tested. Our calculations clearly showed that the band gap will match the experimental value when U = 3.3 eV, as shown in figure 15. Figure 15: Band gap comparison between PBE method (left) and PBE + U (right). 31 For the calculation without Hubbard correction (left), the band gap around Fermi energy is relatively very small. The smallest gap at Gamma point is about 0.25 eV. This does not match with the experimental values nor the physical property since Co3O4 is a semi-conductor. On the right side of Figure 15, the band gap is opened up about 1.5 to 1.6 eV. Band Structure & Density of States The density of states for CrxCo3-xO4 structures are calculated. The band structure diagrams, and corresponding density of states are shown in Figure 16-21 (a). The Fermi energy is set to 0 eV and is shown as the dashed line. The x-axis shows the K-point path in the first Brillouin zone. Co3O4 are semiconductor material, the Fermi energy is right in between the valence band and the conduction band indicates the semiconductor properties with the Hubbard correction. The projected density of states is also shown in Figure 16-23 (right). The energy states are projected to different atomic orbitals to show the contribution from each orbital at the same energy level. In these structures, only 3d orbitals in Co3+, Co2+, Cr3+ and 2 p orbitals from O2- were shown. The rest of the atomic orbitals will have the contribution far from the Fermi energy, which is not significant in this study and hence not included in the plots. 32 (a) (b) Figure 16: (a) Band structure and density of states for Co3O4. (b) Projected density of states for Co3O4. In Figure 16 right, peaks above and below 0 in the y-axis are corresponding to the majority spin and minority spin (spin up and spin down) respectively. The majority and minority states are perfectly symmetric, this can be explained by the anti-ferromagnetic property of the Co3O4 at 0 k. 33 The peaks and pattern of the PDOS is consistent with the reported studies by Chen about the bulk Co3O4 property49. In the valence band close to the Fermi energy (0 eV), the 3d electrons from the octahedral sites contribute most of the energy states, the 3d electrons from the tetrahedral sites only have minor contribution at this energy region. Similarly, the octahedral 3d orbitals contribute the most around the conduction band close to the Fermi energy. In Co3O4, the total magnetism is zero because the electrons from octahedral cobalt sites are all paired. (a) Figure 17: (a) band structure and density of states. (b) Projected density of states for Co2.75Cr0.25O4. 34 Figure 17 (contโ€™d) (b) When Co3+ is substituted with Cr3+, there are unpaired electrons in the system which provides them with ferrimagnetic characteristics. The peak around 4 eV from Cr contribution has a larger majority spin than the minority spin. In Figure 18, the energy states around the Fermi level from Co3+ is relatively large comparing to the Cr3+ with band gap around 1.53 eV. 35 (a) (b) Figure 18: (a) Band structure and density of states (b). Projected density of states for Co2.25 Cr0.75O4. The PDOS in Figure 18 clearly illustrates that the energy states above the Fermi level gets smaller and the shape between spin up and spin down showed a different shape. This is primarily 36 due to the unpaired electrons from Cr3+ have a higher concentration in the structure which increase the magnetization properties of the substrate. (a) (b) Figure 19: (a) Band structure and density of states. (b) Projected density of states for Co1.75 Cr1.25O4. 37 (a) (b) Figure 20: (a) Band structure and density of states. (b) Projected density of states for Co1.5 Cr1.5O4. In Figure 20, the band structure and density of states stays relatively similar to the previous systems having less Cr concentration. However, for the projected density of states on the right side, though the energy states from the oxygen stays the same, the contribution from Co3+ have reduced 38 significantly to about 5 states/eV. For the Cr3+ states, the energy states have increased form 47 states/eV to about 57 states/eV. (a) (b) Figure 21: (a) Band structure and density of states. (b) Projected density of states for CoCr2O4. Figure 21 (a) and (b) showed the structure when all the cobalt atoms on the octahedral site (Co3+) substituted by chromium atoms. The energy states from Co2+ atoms stay relatively stable 39 comparing to the PDOS in figure 20. Because the contribution of the octahedral Co showed most of the energy levels around the Fermi region, it is reasonable to state that the band gap will be increasing while the scale of the octahedral cobalt decreasing. The band gap for this group of systems starts around 1.50 eV and the structure with all octahedral sites having Cr3+ (when x = 2) CoCr2O4, the band gap has increased to 1.87 eV. The trend of the band gap changing with different magnitude of chromium in a unit cell is shown in Figure 22. The band gap stays relatively the same when the quantity of Cr atoms is small. However, when x > 1.25, the band gap starts to increase with increasing amounts of Cr atoms in the conventional structure. This trend was discussed in the previous figure. Since the Cr quantity increases, the energy states closer from the Fermi level is decreasing. However, there is a threshold for the number of octahedral cobalt to show the change on band gap significantly. It is when x > than 1.25. Figure 22: Band gap trend with different amount of Cr substitution. 40 CHAPTER 5: ONGOING CALCULATIONS Super Cell with OER Intermediates The supercell shown in figure 23 with different intermediates will be used for computing the Gibbs free energy along each OER step. It is important to highlight that there are two different terminations of the unit cell.43 The surface can be defined as termination A and the following layer beneath can be defined as termination B. The termination A has Co3+, Cr3+ and Co2+ exposed to the vacuum however termination B only has Co3+ and Cr3+, which means the tetrahedral coordinated cobalt only exits on termination B. Since termination A has all three kinds of metal sites, it was chosen to be the exposed surface in order to investigate the effect of chromium doping properties. Notice that for the surface layers, there both cobalt and chromium are exposed. Theoretically, the water molecule can potentially bind on either site. To study the difference binding properties, three different sets of intermediates will be prepared: intermediates with cobalt active site and intermediates with chromium active site. Gibbs free energy will be calculated for both cycles for OER and compare the results. 41 (a) (b) (c) Figure 23: Different intermediates for the (110) cleaved supercell. (a) Cr*O; (b) Cr*OH; (c) Cr*OOH. In the supercell calculation, because the bottom layers are far from the active surface, the effect from those layers will be ignored. In this case, the bottom layers will be fixed during the intermediateโ€™s optimization. The Hubbard values for a doped Co3O4 system is valuable and will be tested and used in future calculations. Because the values are empirical, the Hubbard values for slab models The active layers will be tested to figure out the best terminations between termination A and B in order to confirm the active sites. 42 CHAPTER 6: FUTURE WORK Once the overpotentials are calculated, the goals are to study as many OER systems as possible. For example, instead of examining the chromium doped system, the iron doped, and nickel doped cobalt oxide will be considered. In figure 3 showed that nickel and iron doped Co3O4 will significantly enhance the catalytic activity. The reason for those dopants to enhance the catalytic activity might be the dopants will also prefer to substitute the octahedral sites, which has been considered to be the actual sites interacting with oxygen molecule.15 This will be confirmed by performing the dopants occupation analysis calculations as mentioned above, and the overpotential for different active sites will also be calculated to support the hypothesis. The methodology will be similar: Calculate the Gibbs free energy for each step and determine the overpotential. The dopants effect on the OER mechanism will be investigated by comparing the overpotentials of different active sites. Moreover, the bond length effect in OER especially the bond length between Co โ€“ Co will be investigated. The bond effect will be beneficial to understand the catalytic trend about the changes in catalytic activity when changing the concentration of Cr โ€˜xโ€™. As mentioned in the previous section, the local density of states on certain small particles will also be useful to understand the binding property on the exposed surface. The result from this project in terms of the projected density of states for the conventional structures will be a strong reference when analyzing the surface bindings for different intermediates. 43 REFERENCES 1. N. Sรถnnichsen, Statista, https://www.statista.com/statistics/280704/world-power- consumption (Accessed May 11, 2021). 2. Baz, K.; Cheng, J.; Xu, D.; Abbas, K.; Ali, I.; Ali, H.; Fang, C. Asymmetric Impact of Fossil Fuel and Renewable Energy Consumption on Economic Growth: A Nonlinear Technique. Energy (Oxf.) 2021, 226 (120357), 120357. 3. Handbook of Nanomaterials for Wastewater Treatment: Fundamentals and Scale up Issues; Bhanvase, B. A., Sonawane, S., Pawade, V. B., Pandit, A. B., Eds.; Elsevier Science Publishing: Philadelphia, PA, 2021. 4. Pashchenko, D.; Nikitin, M. Forging Furnace with Thermochemical Waste-Heat Recuperation by Natural Gas Reforming: Fuel Saving and Heat Balance. Int. J. Hydrogen Energy 2021, 46 (1), 100โ€“109. 5. Yuan, C.; Liu, Z.; Gu, W.; Teng, F.; Hao, W.; Hussain, S. A.; Jiang, W. Hydrogen Production Performance of Novel Glycerin-Based Electrolytic Cell. Renew. Energy 2021, 167, 862โ€“868. 6. Shang, X.; Tang, J.-H.; Dong, B.; Sun, Y. Recent Advances of Nonprecious and Bifunctional Electrocatalysts for Overall Water Splitting. Sustain. Energy Fuels 2020, 4 (7), 3211โ€“3228. 7. Fadhel, A. Z.; Pollet, P.; Liotta, C. L.; Eckert, C. A. Combining the Benefits of Homogeneous and Heterogeneous Catalysis with Tunable Solvents and Nearcritical Water. Molecules 2010, 15 (11), 8400โ€“8424. 8. Sabatier, P. La Catalyse En Chimie Organique, Encyclopรฉdie de Science Chimique Appliquรฉe; Ch. Bรฉranger: Dunham, QC, 1913. 9. Reier, T.; Oezaslan, M.; Strasser, P. Electrocatalytic Oxygen Evolution Reaction (OER) on Ru, Ir, and Pt Catalysts: A Comparative Study of Nanoparticles and Bulk Materials. ACS Catalysis 2012, 2 (8), 1765โ€“1772. 10. Seitz, L. C.; Dickens, C. F.; Nishio, K.; Hikita, Y.; Montoya, J.; Doyle, A.; Kirk, C.; Vojvodic, A.; Hwang, H. Y.; Norskov, J. K.; Jaramillo, T. F. A Highly Active and Stable IrOx/SrIrO3 Catalyst for the Oxygen Evolution Reaction. Science 2016, 353 (6303), 1011โ€“1014. 11. Diaz-Morales, O.; Ledezma-Yanez, I.; Koper, M. T. M.; Calle-Vallejo, F. Guidelines for the Rational Design of Ni-Based Double Hydroxide Electrocatalysts for the Oxygen Evolution Reaction. ACS Catalysis 2015, 5 (9), 5380โ€“5387. 44 12. Sun, W.; Song, Y.; Gong, X.-Q.; Cao, L.-M.; Yang, J. An Efficiently Tuned D-Orbital Occupation of IrO2 by Doping with Cu for Enhancing the Oxygen Evolution Reaction Activity. Chem. Sci. 2015, 6 (8), 4993โ€“4999. 13. Arai, H.; Hayashi, M. SECONDARY BATTERIES โ€“ LITHIUM RECHARGEABLE SYSTEMS โ€“ LITHIUM-ION | Positive Electrode: Lithium Cobalt Oxide. In Encyclopedia of Electrochemical Power Sources; Elsevier, 2009; pp 258โ€“263. 14. Banerjee, S.; Debata, S.; Madhuri, R.; Sharma, P. K. Electrocatalytic Behavior of Transition Metal (Ni, Fe, Cr) Doped Metal Oxide Nanocomposites for Oxygen Evolution Reaction. Appl. Surf. Sci. 2018, 449, 660โ€“668. 15. Peng, Y.; Hajiyani, H.; Pentcheva, R. Influence of Fe and Ni Doping on the OER Performance at the Co3O4(001) Surface: Insights from DFT+U Calculations. ACS Catal. 2021, 11 (9), 5601โ€“5613. 16. Lin, C.-C.; McCrory, C. C. L. Effect of Chromium Doping on Electrochemical Water Oxidation Activity by Co3โ€“XCrxO4 Spinel Catalysts. ACS Catal. 2017, 7 (1), 443โ€“451. 17. Sholl, D.; Steckel, J. A. Density Functional Theory: A Practical Introduction; Wiley- Blackwell: Hoboken, NJ, 2009. 18. Hohenberg, P., & Kohn, W. (1964). Inhomogeneous Electron Gas. Physical Review, 136(3B), B864โ€“B871. doi:10.1103/physrev.136.b864 19. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78 (7), 1396โ€“1396. 20. Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B Condens. Matter 1992, 46 (11), 6671โ€“6687. 21. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865โ€“3868. 22. Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91 (14), 146401. 23. Hubbard, J. Electron Correlations in Narrow Energy Bands. Proc. R. Soc. Lond. 1963, 276 (1365), 238โ€“257. 24. Slater, J. C. Atomic Shielding Constants. Physical Review 1930, 36 (1), 57โ€“64. 45 25. Electronic Wave Functions - I. A General Method of Calculation for the Stationary States of Any Molecular System. Proc. R. Soc. Lond. 1950, 200 (1063), 542โ€“554. 26. N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College Publishing, Philadelphia, PA, 1976). 27. Rรฅsander, Mikael. (2010). A Theoretical Perspective on the Chemical Bonding and Structure of Transition Metal Carbides and Multilayers. 28. Sachs, E. S. Frozen Core Approximation, a Pseudopotential Method Tested on Six States of NaH. J. Chem. Phys. 1975, 62 (9), 3393. 29. Hamann, D. R.; Schlรผter, M.; Chiang, C. Norm-Conserving Pseudopotentials. Phys. Rev. Lett. 1979, 43 (20), 1494โ€“1497. 30. Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B, 41, 7892-7895 (1990). 31. Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. rev. 1976, 13 (12), 5188โ€“5192. 32. Man, I. C.; Su, H. Y.; Calleโ€Vallejo, F.; Hansen, H. A.; Martรญnez, J. I.; Inoglu, N. G.; Kitchin, J.; Jaramillo, T. F.; Nรธrskov, J. K.; Rossmeisl, J. Universality in Oxygen Evolution Electrocatalysis on Oxide Surfaces. ChemCatChem 2011, 3 (7), 1159โ€“1165. 33. Liang, Q.; Brocks, G.; Bieberle-Hรผtter, A. Oxygen Evolution Reaction (OER) Mechanism under Alkaline and Acidic Conditions. J. Phys. Energy 2021, 3 (2), 026001. 34. 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 Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108 (46), 17886โ€“17892. 35. Roth, W. L. Journal of Physics and Chemistry of Solids 1964, 25, 1โ€“10. 36. Gulino, A.; Dapporto, P.; Rossi, P.; Fragalร , I. A Novel Self-Generating Liquid MOCVD Precursor for Co3O4Thin Films. Chem. Mater. 2003, 15 (20), 3748โ€“3752. 37. SMXTn, B. W. L.; HoBsoy, A. D.; Warren Spring Laboratory, P. O. The Structure of Cobalt Oxide, Co~O4. SG 1 2. 38. Zhou, J.; Zhang, L.; Huang, Y.-C.; Dong, C.-L.; Lin, H.-J.; Chen, C.-T.; Tjeng, L. H.; Hu, Z. Voltage- and Time-Dependent Valence State Transition in Cobalt Oxide Catalysts during the Oxygen Evolution Reaction. Nat. Commun. 2020, 11 (1), 1984. 46 39. Suen, N.-T.; Hung, S.-F.; Quan, Q.; Zhang, N.; Xu, Y.-J.; Chen, H. M. Electrocatalysis for the Oxygen Evolution Reaction: Recent Development and Future Perspectives. Chem. Soc. Rev. 2017, 46 (2), 337โ€“365. 40. Cococcioni, M.; de Gironcoli, S. A Linear Response Approach to the Calculation of the Effective Interaction Parameters in the LDA+U Method. arXiv [cond-mat.mtrl-sci], 2004. 41. Chen, J.; Wu, X.; Selloni, A. Electronic Structure and Bonding Properties of Cobalt Oxide in the Spinel Structure. Phys. Rev. B Condens. Matter Mater. Phys. 2011, 83 (24). 42. Maldonado, F.; Novillo, C.; Stashans, A. Ab Initio Calculation of Chromium Oxide Containing Ti Dopant. Chem. Phys. 2012, 393 (1), 148โ€“152. 43. Chen, J.; Selloni, A. Electronic States and Magnetic Structure at the Co3O4 (110) Surface: A First Principles Study. arXiv [cond-mat.mtrl-sci], 2012. 44. Lhermitte, C. R.; Sivula, K. Alternative Oxidation Reactions for Solar-Driven Fuel Production. ACS Catal. 2019, 9 (3), 2007โ€“2017. 45. Walsh, A.; Wei, S.-H.; Yan, Y.; Al-Jassim, M. M.; Turner, J. A.; Woodhouse, M.; Parkinson, B. A. Structural, Magnetic, and Electronic Properties of the Co-Fe-Al Oxide Spinel System: Density-Functional Theory Calculations. Phys. Rev. B Condens. Matter Mater. Phys. 2007, 76 (16). 46. Wang, L.; Maxisch, T.; Ceder, G. Oxidation Energies of Transition Metal Oxides within TheGGA+Uframework. Phys. Rev. B Condens. Matter Mater. Phys. 2006, 73 (19). 47. Shinde, V. R.; Mahadik, S. B.; Gujar, T. P.; Lokhande, C. D. Supercapacitive cobalt oxide (Co3O4) thin filmsby spray pyrolysis. Appl. Surf. Sci. 2006, 252, 7487โˆ’7492. 48. Oyewande, O. E.; Atsue, T.; Ogunniranye, I. B.; Usikalu, M. Prediction of Lattice Constants of Some Transition Metal Nitrides Using Different Functionals and Pseudo- Potentials. IOP Conf. Ser. Earth Environ. Sci. 2021, 655 (1), 012045. 49. Zhandun, V. S.; Nemtsev, A. Ab Initio Study of the Magnetic, Optical and Electronic Properties of Spinel Co3O4 within DFT and GW Approaches. J. Magn. Magn. Mater. 2020, 499 (166306), 166306. 50. Yoon, M.; Yang, S.; Hicke, C.; Wang, E.; Geohegan, D.; Zhang, Z. Calcium as the Superior Coating Metal in Functionalization of Carbon Fullerenes for High-Capacity Hydrogen Storage. Phys. Rev. Lett. 2008, 100 (20), 206806. 47