Last weeks I have been trying to learn how to use first-principles calculations to study electrochemical reactions. I am specifically using software for Density Functional Theory (DFT) calculations with the goal to apply it to interesting electrocatalytic reactions such as the hydrogen evolution reaction (HER), oxygen evolution reaction (OER) or carbon dioxide reductio reaction (CO2RR) among others. Even if I am not using DFT for my current research, which is more applied, I realize the growing importance of computational chemistry and how this could be very relevant in my future research.
During this learning path I am encountering some difficulties since I am coming from a quite applied background. For instance, there are lots of different computational software packages able to perform DFT calculations, but it is not clear the differences between them and which one could be more useful in your specific studies. It is also very difficult to find examples of real code employed in calculations reported in the literature, even in the current era of open source, open data and open publishing. So, there are a few roadblocks to overcome for a complete beginner trying to learn on his own. I plan to write a series of posts in the blog that will serve me to clarify my ideas during my learning path and, who knows, it may also be useful for someone out there. Please, bear with me since I am a beginner and feel free to correct me if I write something incorrect. FEEDBACK IS MUCH APPRECIATED. I will try to write in the most simple terms in order to make this tutorial helpful for non-experienced users (like me). However I will deal with some theory, so a few equations are mandatory. Since I will perform the calculations with my personal computer, I will use very simplistic molecular models, and sometimes I will skip the contribution of some parameters to make shorter posts even if it is not scientifically correct in that moment. I will deal with those specific parameters and calculations later in specific posts.
This will be a very practical tutorial, so I recommend to read some seminal literature in the field such as:
- Essentials of Computational Chemistry by Chris Cramer
- Introduction to Computational Chemistry by Frank Jensen
- Electrocatalysis – Theoretical Foundations and Model Experiments edited by Richard C. Alkire, Dieter M. Kolb, Ludwig A. Kibler and Jacek Lipkowski
- Handbook of Materials Modelling edited by Wanda Andreoni and Sidney Yip
Density Functional Theory
In an oversimplified way, the Schrödinger equation contains all the information about the energetic state of any chemical system, and therefore, it describes the properties of that system. Since solving the Schrödinger equation is not practically possible for almost any system (just for very basic ones such as the hydrogen atom), there are several approximations that enable the solving of this equation. Density Functional Theory (DFT) is one of these approximations. Its foundation lies basically in that the electron density of any system determines all the ground-state properties of the system (Hohenburg-Kohn theorem), and by using electron density functionals, the total energy of the system can be obtained. You could find more accurate definitions in any of the literature mentioned above. DFT computational methods have been gaining a lot of popularity in research in the past few years thanks to the balance between low computational cost and accuracy, enabling its use to study molecules, solids, heterogeneous reactions and so on.
Software packages for DFT
There are many computational packages able to perform DFT calculations. Some of the most known are Gaussian, Quantum Expresso, Orca, Nwchem, PSI, among many other. Some of them are available for free, some are also open source, and some are commercial. For electrochemical applications, maybe the most used is the Vienna Ab Initio Simultation Package (VASP), but this is a commercial software. Other alternatives that have been used to study electrochemical reactions are the Atomic Simulation Environment (ASE) integrated with packages such as GPAW, which uses the projector-aumented wave method. Dacapo is also a good option integrated with ASE for computational calculations.
I do not have very clear the advantages of one package over others, but I will be using GPAW/ASE since together they make a very powerful combo for designing molecular systems and perform several types of calculations. Many literature about simulated electrocatalytic reactions using GPAW/ASE have been reported, which is also a plus to compare results. Furthermore, they are open-source, which means the code is available, and they are written in Python, which I know how to code. ASE also supports other software to make the calculations, so it would be possible to reuse quite a bit of your code in case you would like to change the calculator at some point.
In this post, I will try to describe the necessary steps to generate a free energy diagram for an electrochemical reaction at different potentials in order to visualize if the reaction is thermodynamically favorable. This is also one of the most useful strategies to study the feasibility of reaction pathways and mechanisms using specific materials as electrode surfaces, and it is even possible to anticipate the performance of a tailored designed material.
Many important electro electrocatalytic reactions such as the HER, OER, CO2RR or oxygen reduction reaction among others have reaction steps the involve the transfer of electrons and protons, as in the following model reaction:
X* + H+(aq) + e- → XH
where * indicates that the species is adsorbed on the electrode surface. In order to generate the free energy diagram, we would need to calculate the free energy of the intermediates and products of every step of the reaction. For this specific case, the change in the free energy for this reaction would be:
ΔG = G(XH) - G(X*) - G(H+) - G(e-)
Calculating the free energy of electrons and protons is not an easy task, so there are several strategies to facilitate the computational calculations. One of the most employed is the use of the Computational Hydrogen Electrode (CHE) model, which is useful for reactions involving steps of electron and proton transfers such as the mentioned ones. The CHE model deals with the free energy issues of protons and electrons by considering the standard hydrogen electrode (SHE) as reference. The SHE involves the following reaction in equilibrium:
H+(aq) + e- ⇄ 1/2 H2 (g)
Since we are considering standard and equilibrium conditions, i.e. 1 M H+ concentration (pH = 0), 1 bar for the H2 gas, 25 ºC and 0 V, then the chemical potentials of both reaction terms are the same:
µ(H+) + µ(e-) = 1/2 µ(H2(g))
Therefore, it is possible to employ the free energy of H2(g) instead of the free energies of protons and electrons. The change of the free energy of the studied reaction step would then be:
ΔG = G(XH) - G(X*) - 0.5 G(H2(g))
This is for standard conditions (0 V vs SHE and pH=0). In order to consider the effect of the applied potential and pH, two factors should be added to the equation:
-eU (where e is the electron charge, and U is the applied potential vs SHE)
+ kB T log(10) pH (where kB is the Boltzmann constant, and T the temperature)
The CHE model assumes several situations: 1) the adsorption energy of the intermediates on the electrode surface is independent of the electrostatic field; 2) the bulk of the electrolyte is in equilibrium; 3) the surface where the reactions take place is in equilibrium with the electrodes and the bulk of the electrolyte. Therefore, it does not take into account the effect of the electrostatic potential at the double layer, which may have some influence on adsorption energies specially for large molecules with strong dipoles. It does also not consider activation barriers, i.e. kinetics, and the data is purely thermodynamic. And it does not consider molecule solvation effects. All these situations complicate a bit more the calculations, and therefore, I will address these specific situations in following posts.
Making a DFT calculation
Let’s apply this theory and make some DFT calculations to obtain the free energy diagram of the OER at a simple metallic surface. OER pathway is usually described by the following steps:
* + H2O → OH* + H+ + e- (step 1)
OH* → O* + H+ + e- (step 2)
O* + H2O → OOH* + H+ + e- (step 3)
OOH* → O2 + H+ + e- (step 4)
The Gibss free energy for these steps can be expressed as:
ΔG1 = ΔGOH
ΔG2 = ΔGO - ΔGOH
ΔG3 = ΔGOOH - ΔGO
ΔG4 = 4.92 eV - ΔGOOH
We avoid the use of free energies of protons and electrons, since we are using the CHE model and ΔGf(H2(g)) would be zero. The value 4.92 eV considers the thermodynamic value of the overall reaction, which is 4 electron transfers and 1.23 V as equilibrium potential. This way, calculation of O2 energy is avoided, which seems to be difficult to calculate under some conditions.
To calculate the change in the free energy of any of these terms, we use the famous thermodynamic equation:
ΔG = ΔH - TΔS
Now, let’s make some assumptions and utterly incorrect simplifications for the sake of make this post simple. Firstly, let’s assume that the enthalpy term comes defined only by adsorption energies of the adsorbate on the electrode surface (ΔE) and the change in zero point energy (ZPE). This is not totally correct, but sometimes it can be applied to minimize the complexity of calculations. Then, we would have:
ΔGi = ΔEi + ΔZPEi - TΔSi
Entropic terms can be obtained from thermodynamic tables since the contribution of adsorbed species are not usually considered (ΔS for *, *OH, *O, *OOH is 0 eV). Zero point energies can be obtained from DFT calculations using vibrational properties, but I will deal with those calculations in a future post. In order to avoid making this post utterly long and complicated, let’s assume values described in literature for other system. This is scientifically INCORRECT, and the values obtaned this way will probably be meaningless. To calculate the free energy of the reaction steps we have to consider the reactants and products of those reactions (see equations above). Just remember that we convert protons and electrons to hydrogen due to the CHE mode. I think that it is also possible to exclude the contribution of the ZPE(*) (I remember to have read about this but do not remember where). So, for the first step of the reaction:
ΔGOH = ΔEOH + ΔZPEOH - TΔSOH
ΔEOH = E(*OH) - E(*) - [E(H2O) - 1/2 E(H2)]
ΔZPEOH = ZPE(*OH)
- ZPE(*) - [ZPE(H2O) - 1/2 ZPE(H2)]
TΔS(*OH) - TΔS(*) - [TΔS(H2O) - 1/2 TΔS(H2)] = -TΔS(H2O) + 1/2 TΔS(H2)
To determine the energy change for OH adsorption on the surface, we would make the DFT calculations to calculate the energy of the whole system, E(*OH), the energy of the surface, E(*) and the total energy of a gas-phase reference system (H2O and H2). We then calculate the free energy change for all the intermediates of the OER reaction in a similar way, and we can generate the free energy diagram from the results. As I said our results are physically meaningless because I did too many simplifications in the calculations, but this would be the practical procedure to generate this useful representation. I also did not consider the surface coverage of the adsorbate on the metal, which can also have influence on the free energy calculations. I created an ASE/GPAW python script to calculate the energies for these systems (*, *OH, *O, *OOH, H2O and H2) using a Ni(111) surface. You can find it in a github repository. In order to avoid high computational costs, I have designed the surface with only 8 atoms, seven of them were fixed and one was allowed to relax for the optimization. Also, the vacuum chosen to account for the periodicity was only 2.0 Å since higher values make the calculations a lot slower. In a real calculation, the surface should be larger and the most similar possible to the hypothetical experimental situation. But, obviously, you would need a more powerful machine to perform the calculations. So, the free energy diagram obtained at 0 V and the optimized geometries of the intermediates were the following ones:
Calculating the free energy diagram at different potentials
We already have a free energy diagram at 0 V (vs SHE) and pH=0, since we are using the CHE model. To calculate the diagram at different potentials, we consider the eU term for each reaction step, which is straightforward:
ΔGi (U vs RHE) = ΔGi (0 V vs RHE) - eU
This way, we could calculate free energy diagrams at different applied potentials:
Here, we can see how the reaction become thermodynamically favorable when the applied potential is +2.87 V, having a theoretical overpotential of +1.64 V (η = E – Eeq, with Eeq = 1.23 V for the OER). This overpotential can be calculated beforehand using the following expression from the free energies at 0 V:
η = (max (ΔG1, ΔG2, ΔG3, ΔG4) / e) - Eeq
For this system, the formation of *OOH (step 3) is the limiting step of the reaction with ΔG3 = 2.87 eV, which will be responsible for this overpotential. As previously mentioned, this is only considering the final free energies (thermodynamics) and the activation barriers (kinetics) would also affect the overpotential needed to carry out the reaction experimentally. Activation barriers will be addressed in a following post.
In this post, I have described the typical approach considered to study an electrochemical reaction using DFT computational calculations and how to represent the free energy diagram of the reaction. I have used very simplistic models and avoided some calculations to make this post easier to understand. I will deal with all the possibilities and more correct approaches in future posts.