Part 6: Oxygen Evolution Reaction

Definition of the Problem Space

Finally, we are ready to study some reaction steps of the oxygen evolution reaction (OER).

The "canonical" hypothesis for what might be the key reaction steps is summarized in Ref. https://dx.doi.org/10.1021/ja301567f. We adopt this hypothesis here.
We furthermore pursue the hypothesis that a particular site at our hematite-derived nanocluster model created earlier could be a possible reaction site.
We additionally do not consider transition states in our analysis below (due to time limitations), hypthesizing that the sequence of equilibrium structures will give us an initial idea of how the reaction might proceed.
We finally employ the hypothesis of a cluster geometry in vacuum, without considering solvation effects implicitly or explicitly. There is abundant literature on the limits of this hypothesis and how to overcome it, but as an important advantage, this hypothesis allows us to consider a limited set of local minimum-energy geometries, rather than requiring us to pursue statistical averages over a potentially large solvent conformational ensemble even for a single hypothetical reaction site.

Whether or not all of the above hypotheses are meaningful or sufficient is an excellent question, but for the present purpose of demonstrating a set of specific simulations, the above hypotheses provide us with a well defined starting point.

Reaction Steps

Water splitting is an endothermic chemical reaction. By providing the right amount of energy, water can be decomposed into its constituent atoms, oxygen and hydrogen, through the following chemical reaction:

\(2H_{2}O \rightarrow O_{2} + 2H_{2}\)

The measured minimum energy required for this reaction is \(1.23\) eV.

We assume the following reaction scheme for the oxygen-evolution half-reaction of water splitting at hematite clusters:

(A) \(H_{2}O + ^{*} \rightarrow ^{*}{OH}_{2}\) ,

(B) \(^{*}OH_{2} \rightarrow ^{*}OH + H^{+} + e^{-}\) ,

(D) \(H_{2}O + ^{*}O \rightarrow ^{*}OOH+ H^{+} + e^{-}\) ,

(E) \(^{*}OOH \rightarrow ^{*} + O_{2} + H^{+} + e^{-}\) .

Here, \(^{*}\) represents the bare cluster. The symbols \(^{*}OH_{2}\), \(^{*}OH\), \(^{*}O\), and \(^{*}OOH\) represent the cluster with corresponding chemisorbed species. During each oxidation step (B-E), a proton is removed into the electrolyte and an electron is created. To close the water splitting reaction, a proton will eventually combine with an electron in the hydrogen-evolution half-reaction (not shown) of the initial net reaction equation.

Simulation Tasks

As a first step, the potential effectiveness of hematite clusters for water splitting can be evaluated by examining the energetics of above reactions steps. To proceed, we use the cluster structure that was relaxed with help of hydrogen saturation (cf. Tutorial-5/solutions/relaxation-H-saturated). We pick one (of potentially many possible) reaction site and calculate total energies of the adsorbed reaction intermediates on the right and left hand sides of the reaction equations (A)-(E).

As before, we need to identify a local minimum-energy structure of each system. This task is accomplished by first preparing a reasonable initial structure, charge, and spin initialization, and then relaxing that structure until it reaches a minimum energy configuration (for which the magnitudes of the any residual forces on each atom is below a given threshold, specified in control.in).

The above systems can be divided into to categories:

Isolated molecules, i.e. \(H_{2}O\), \(O_{2}\), and \(H_{2}\)
Cluster plus chemisorbed species

Now, we deal with each category one after another.

Isolated molecules

See Tutorial Basics of Running FHI-aims for a more specific discussion of simple molecules, such as:

H\(_{2}\)
O\(_{2}\)
H\(_{2}\)O

and how to deal with them.

The solutions for molecular structures and total energies ("light" settings) for Part 6 can be found in solutions/molecules.

Cluster with chemisorbed species

Let's begin with the bare cluster. Our first task is to identify a reactive Fe atom which might have a high tendency to participate in chemical reaction with corresponding groups in the above steps. In principle, although not trivially, this can be done systematically by using further information (e.g., local electrostatic potentials, local bonding geometries, local densities of states) extracted from first-principles simulations.

However, here, by visually analyzing the structure using GIMS, we see that atom # 32 has fewer number of neighboring atoms compared to other Fe atoms. Therefore, one possible hypothesis is that this atom might have a tendency to participate in chemical reactions. We adopt it as our reaction site, keeping in mind that we might have to address other potential reaction sites in an analogous manner later.

Next, we test two separate approaches to find meaningful structures of the cluster with chemisorbed reaction intermediates.

In our first approach, we derive the structure sequentially from step A to E, i.e. we initialize the structure for the next step, using the relaxed structure obtained from the previous step. This amounts to the hypothesis that individual reaction steps happen in sequence and the system retains a "memory" of its history from one step to the next.
In our second approach, we initialize each structure separately, using the initial,bare cluster. This amounts to the hypothesis that the partial reactions occur rather independently of one another, i.e. the "history" of successive reaction steps is not remembered and the system is thermodynamically "reset" between individual reaction steps.

Approach 1: Reaction Steps with "Memory"

Create a new folder for the bare cluster. From folder Tutorial-5/solutions/relaxation-H-saturated/step-2-H-removed copy geometry.in.next_step into it. We will use this structure to prepare an initial structure for reaction A.
Before proceeding further, an appropriate initial_moment and initial_charge should be assigned to atoms in this file. To avoid doing this manually, we can use clims together with Tutorial-5/solutions/relaxation-H-saturated/step-2-H-removed/geometry.in. The geometry.in file has the initial_moment and initial_charge information but it is not the correct final atomic structure. Therefore, please note that we DO NOT read the atomic coordinates from this latter file.
So, let's copy the geometry.in file to geometry.in.initialized and then execute clims-reinitialize-geometry geometry.in.initialized geometry.in.next_step. This action will generate a new file, geometry-reinitialized.in, which has the correct initial atomic coordinates and initial_moment and initial_charge information of the bare cluster.
You can also copy Tutorial-5/solutions/relaxation-H-saturated/step-2-H-removed/aims.out into this folder. We will use this file later to extract total energy of the bare cluster.
Create a new folder for reaction A. Copy the geometry-reinitialized.in of the previous step to a file geometry.in. Then we can import this file to GIMS and add a H\(_{2}\)O molecule to atom # 32 (by adding to and editing the atom list). Finally, export the structure. We used the following coordinates for the H\(_{2}\)O molecule:
```
atom -1.643100 5.572560 2.710410 O
atom -1.366190 6.227120 2.061320 H
atom -2.476890 5.903780 2.454920 H
```
Prepare the control.in file for the relaxation (e.g., using GIMS). To reduce the computational workload, we use light basis sets and pbe exchange-correlation functional. Important parameters controlling this simulation are:
```
xc                                 pbe
relativistic                       atomic_zora scalar
relax_geometry                     bfgs         5e-3
spin                               collinear
```
Once the calculation is finished, set the initial_moment and initial_charge values in geometry.in.next_step using the clims-reinitialize-geometry geometry.in geometry.in.next_step command. Then create a folder for reaction B and take the relaxed structure from geometry-reinitialized.in and copy it into geometry.in in this new folder. Use GIMS to remove the outer H atom (atom # 82.
Export the new structure and prepare the control.in, same as before, and now you are ready to obtain relaxed structure of step B. The coordinates of our OH group are:
```
atom        -1.21797125        5.54596620        1.40418590  O
atom        -0.96502794        4.99126465        0.56822231  H
```
We proceed in the same way until final structures of all reaction steps are obtained.

GIMS has a powerful output analyzer feature, which can be used to see how cluster structure evolves throughout the relaxation process.

The solutions for relaxation of above systems are located in solutions/cluster/approach-1.

Approach 2: Reaction Steps without "Memory"

Repeat the steps listed in Approach 1, with the exception that this time, the initial structure for each reaction intermediate is prepared from the bare cluster structure, not from a preceding reaction step.

Relaxed structures using this approach can be found in solutions/cluster/approach-2.

Calculating the reaction energies

After obtaining the relaxed structures of all systems involved in reactions A-E, we are now ready to calculate reaction energies. The reaction (free) energy of \(A+B \rightarrow C\) is calculated by approximating \(\Delta G = E(C)-E(A)-E(B)\). Here, the thermodynamic related contributions, such as Zero point energy (ZPE) and entropic contributions, are not taken into account. Also, in practice, the proton and electron transfer process in reactions B-E takes place under an applied external bias \(\phi\), the effect of which could be added to \(\Delta G\) by including a \(-e.\phi\) term. In the following, we only focus on reaction energies calculated from total energies of our relaxed structures.

For starters, let's calculate free energy of \(2\textrm{H}_{2}\textrm{O} \rightarrow \textrm{O}_2 + 2\textrm{H}_2\):

Extract the total energy of corresponding relaxed structures. The final total energy can, for example, be done by grep uncorrected aims.out | tail -1 where aims.out is the output file of the corresponding relaxation simulation.
The free energy is calculated as \(\Delta G_{H_{2}O} = E(O_{2}) + 2E(H_{2}) - 2E(H_{2}O)\).

How does the DFT water splitting energy, \(\Delta G_{H_{2}O}/4\), compare to the measured \(1.23\) eV value?

Following similar steps, calculate the free energies of reaction A-E which can be obtained by:

\(\Delta G_{A} = E(^{*}OH_{2}) - E(^{*}) - E(H_{2}O)\)
\(\Delta G_{B} = E(^{*}OH) - E(^{*}OH_{2}) + E(H^{+} + e^{-})\)
\(\Delta G_{C} = E(^{*}O) - E(^{*}OH) + E(H^{+} + e^{-})\)
\(\Delta G_{D} = E(^{*}OOH) - E(^{*}O) - E(H_{2}O) + E(H^{+} + e^{-})\)
\(\Delta G_{E} = E(^{*}) + E(O_{2}) - E(^{*}OOH) + E(H^{+} + e^{-})\)

\(E(H^{+} + e^{-})\) is implicitly calculated using the standard hydrogen electrode (SHE) model, \(1/2H_{2} \rightarrow H^{+} + e^{-}\), so that \(E(H^{+} + e^{-}) = 1/2 E_{H_{2}}\).

Once you have all the reaction energies at hand, try to answer the following questions:

What is the estimated cumulative reaction energy for all reactions, \(\Delta G_{all} = \Delta G_{A} + \Delta G_{B} + \Delta G_{C} + \Delta G_{D} + \Delta G_{E}\) ?
How does \(\Delta G_{all}\) compare to \(\Delta G_{H_{2}O}\) calculated by DFT and the measured value, and what it the water splitting energy \(\Delta G_{all}/4\)?
Is there any evidence to indicate which of Approach 1 and Approach 2 provides a more realistic picture of water splitting on hematite clusters? (this is a pretty detailed question)

You can run solutions/reaction-energies.py to see various calculated reaction energies outlined above.