The Surface Phase Diagram tab is an app on the Lykoi dashboard that allows users to load a .csv file containing information regarding the surface energies of various facets and terminations purely as a function of chemical potential and plots a phase diagram based on the user's desired reaction conditions. In this tutorial, I will go over how to use the Surface Phase Diagram tab. I will also go through some fundamental thermodynamics to give context to how this app works, however, for a deeper dive into the science, I suggest giving the following a read:
Otherwise, a future tutorial regarding surface thermodynamics will be released in the future.
To follow along with this tutorial, you can download the neccesary .csv file in the Tutorials page.
Opening up our Lykoi dashboard and going to the "Surface Phase Diagram" tab (1) will bring us to an empty Plotly graph (2). We will begin by loading our .csv file either by dragging and dropping or clicking the tab at the bottom (3). As an example, we will be creating a .csv file containing surface energy information regarding the surface stability of $SrTiO_3$ reconstructions on the (001) facet. Data for this example will be taken from:
Table 2 taken from ref 2.
To help you follow along the tutorial, I've provided a .csv file you can download here that lists the surface energy of each termation in Table 2 as a function of chemical potential. More on this later.
The surface energy ($\phi$) is defined as the Landau (or grand) potential per unit of surface area:
$\Omega = U - TS - \Sigma_iN_i\mu_i$ (Landau potential)
$\phi = \Omega / A = $ (surface energy)
Here, $U$ is the internal energy, $T$ is temperature, $S$ is entropy, $i$ is the component making up the system (typically an atom), $\mu_i$ is the chemical potential of component $i$, $N_i$ is the number of components $i$, and $A$ is the surface area. We'll simplify this in the context of ab-initio density functional theory calculations and assume we are at 0 K with no vibrational components, making the internal energy, $U$, the total DFT energy of the system, as such:
This of course assumes our system models the surface of a material as a slab where the denominator of 2 accounts for the fact that a slab is terminated by 2 surfaces. For a detailed explanation of how to model surfaces, see the following link:
https://matgenb.materialsvirtuallab.org/2017/04/03/Slab-generation-and-Wulff-shape.html
We'll reference $\mu_i$ to 0 eV where $\Delta \mu_i = \mu_i-E^{ref}_i < 0$ eV. Here $E^{ref}_i$ is the reference to the chemical potential. Typically, to account for components $i$ in the bulk, we use the DFT energy per formula unit of the bulk as our reference ($\mu_{SrTiO_3} = \mu_{Sr}+\mu_{Ti}+3\mu_{O} = E^{DFT}_{SrTiO_3}$).
For systems where the stoichiometry of system $U$ is inconsistent with the bulk, we will need to account for exccess or defficient species by incorporating the chemical potential of individual atomic components while the reference for the adsorbate is some molecular species as a gas (e.g. $\mu_{Sr} = E^{ref}_{Sr} +\Delta \mu_{Sr}$ and $\mu_{O} = E^{ref}_{O} +\Delta \mu_{O}$ ). And so Eqn. 3 becomes:
Now you may have notice in Table 2 ref2 above gives us surface energy when the chemical potential terms are saturated ($\Delta \mu_i = 0$ eV). We can retrieve the general surface energy as a function of chemical potential by using the following equations from ref2.
*Note here that the surface energy is not normalized by surface area, this is becasue all surfaces considered are modelled with the same surface area and as such, normalization will not change the results.
where $E_{Sr}^{ref}$ is the energy per atom of fcc Sr, $E_{O}^{ref}$ is the energy per atom of gaseous $O_{2(g)}$ where $\Delta \mu_{O} = \Delta \mu_{O_2}/2$, and $\Gamma_{A}^{Ti}$ is the excess amount of atoms A given by:
With $\phi(0)=\phi(\Delta \mu_O=0, \Delta \mu_{Sr}=0)$ and $\Gamma_{A}^{Ti}$, we can rewrite Eqn5 as:
$\phi(\Delta \mu_O, \Delta \mu_{Sr})=\phi(0)-\Gamma_{Sr}^{Ti}\Delta \mu_{Sr}-\Gamma_O^{Ti}\Delta\mu_{O}$
Also keep in mind that the chemical potentials of a bulk system are capped by a limit that is determined by the chemical potential phase diagram of the bulk. The surface is only stable in a chemical potential window where bulk SrTiO3 is stable and as such, the limits should be set based on this. Going outside the chemical potential window will lead to the stabilization of other bulk phase. For more information about how to derive the limits of $\Delta \mu_i$, see ref2. For the Surface Phase Diagram app, we will use the Materials Project API to query the phase diagram of a three component (Sr-Ti-O) system to identify the chemical potentials where bulk $SrTiO_3$ is stable.
We now have a way to rewrite our surface energy as a function of chemical potentials. In your .csv file, each row will correspond to one surface energy and each column will correspond to the coefficient of each chemical potential, the reduced formula of the bulk structure from which the slab is created from, and an additional column corresponding to the constant. The bulk column will tell the MPAPI what phase diagram to query for in order to calculate the appropriate chemical potential limits. We also have one additional column '#' for indexing. This is enough to create your phase diagram, however additional columns will appear on each facet of the phase diagram as hovertext.
import pandas as pd # Import pandas for creating csv files
heifets_SrTiO3_se = [] # Create a list of surface energies to be written to .csv
# simple function to handle parsing tabulated data
def make_se_dict(se0, Sr_cov, O_cov, name):
d = {}
d['structure'] = name
d['const'] = se0 # constant
d['delu_Sr'] = Sr_cov # coefficient for Sr chemical potential
d['delu_O'] = O_cov # coefficient for O chemical potential
d['bulk'] = 'SrTiO3' # coefficient for O chemical potential
return d
# add each surface energy
heifets_SrTiO3_se.append(make_se_dict(4.84, 1/2, 1/2, 'Regular TiO2'))
heifets_SrTiO3_se.append(make_se_dict(13.07, 3/2, 3/2, 'DL TiO2'))
heifets_SrTiO3_se.append(make_se_dict(25.14, 3, 3, '(2x1) DL TiO2'))
heifets_SrTiO3_se.append(make_se_dict(-2.63, -1/2, -1/2, 'Regular SrO'))
heifets_SrTiO3_se.append(make_se_dict(-8.49, -3/2, -3/2, 'DL SrO'))
# lastly add a # to each row
for i, d in enumerate(heifets_SrTiO3_se):
heifets_SrTiO3_se[i]['#'] = i
df = pd.DataFrame(heifets_SrTiO3_se)
df.to_csv('SrTiO3001_surface_energy.csv') # save the dataframe as a .csv file
df
structure | const | delu_Sr | delu_O | bulk | # | |
---|---|---|---|---|---|---|
0 | Regular TiO2 | 4.84 | 0.5 | 0.5 | SrTiO3 | 0 |
1 | DL TiO2 | 13.07 | 1.5 | 1.5 | SrTiO3 | 1 |
2 | (2x1) DL TiO2 | 25.14 | 3.0 | 3.0 | SrTiO3 | 2 |
3 | Regular SrO | -2.63 | -0.5 | -0.5 | SrTiO3 | 3 |
4 | DL SrO | -8.49 | -1.5 | -1.5 | SrTiO3 | 4 |
Finally we can pass our csv file into the Surface Phase Diagram app. Upon doing so, you will notice that the chemical potential table (4) is now occupied by the chemical potential of O and Sr
You will want to fill out each parameter to create the phase diagram. "min" and "max" is the minimum and maximum value of your parameter respectively. x-axis and y-axis are booleans used to designate what axis you wish to assign to what axis. Increment is how fine you would like to sample your phase diagram.
Note that while you are doing this, you will recieve warning messages at the bottom:
After inputting all the appropriate parameters as such:
This will output messages instructing what needs to be inputted into the table. Once you see the following message:
You can proceed to clicking the "Run Phase Diagram" button
Here, we have our phase diagram generated by Lykoi on the left and the phase diagram from ref2 on the right. In the left plot, shaded in regions indicate the region of stability for $SrTiO_3$ while everything outside of this region is grayed out (unstable). This is in contrast to Heifets et al's figure which shows all terminations will stabilize before spilling into the chemical potential region of bulk instability. This is due to the difference in formation energies for the different Sr-Ti-O phases where MP shows a more pessimistic window of stability for $SrTiO_3$. For more information about bulk stability, the user is encourage to explore the DFT methods used in ref2 and compare them to the methods used in MP.
You will see a button here:
Pressing this button will rewrite the surface grand potential as a function of available reaction conditions, providing two new parameters that can be plotted, the partial pressure and temperature. The conversion between chemical potential and reaction conditions is given by the following:
where $k_B$ is Boltzmann's constant, $p_i$ is the partial pressure of species $i$, $p_o$ is the standard pressure, and $\Delta G$ can is obtained from available thermodynamic tables from NIST/JANAF: https://janaf.nist.gov/
Keep in mind this conversion only works for chemical potential species listed on the table.
Subsequently, this makes our phase diagram a function of three parameters (T, $P_{O}$, and $\Delta \mu_{Sr}$). I don't know about you, but visualizing things in three dimensions on a two dimensional plane makes my eyes hurt, so at the moment, we will only ever plot as a function of two parameters. However, we also have the option to set a slider (5) for the third parameter. Doing so will provide us with a way to change the phase diagram as we move the slider.
*Note: Future implementations will allow for more than one slider, however at the moment, we only allow for one. Any remaining parameter will remain constant with the minimum value.
Recall that the surface energy is the grand potential normalized by surface area. At the end of the day, this is just a normalized formation energy. So one might ask "can I just use this to plot the phase diagram for bulk formation energies instead of surface energies so long as I provide the formation energy as a function of chemical potentials?". And the answer to this is yes. You technically can build phase diagrams for bulk structures or anything for that matter, so long as the units of your Landau potential are consistent for all phases you are interested in.