Reaction analysis on heterogeneous catalysts

In the previous section, we learned how to use the NEB method to calculate activation energies using organic molecules as the example. In this section, we will try to analyze the reaction of NO dissociation on the heterogeneous catalyst Rh as a more practical application.

The workflow is as follows. Create the Rh slab structure that will be the catalyst. and adsorb the NO which is the reactant molecule. Then, the states before and after the reaction are created and the activation energy is calculated using the NEB method.

Note that the content of this tutorial is based on the example published in matlantis-contrib.

Setup

[1]:

import numpy as np
import pandas as pd
from IPython.display import display_png
from IPython.display import Image as ImageWidget
import ipywidgets as widgets
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
import glob
from pathlib import Path
from PIL import Image, ImageDraw

import ase
from ase import Atoms, units
from ase.units import Bohr,Rydberg,kJ,kB,fs,Hartree,mol,kcal
from ase.io import read, write
from ase.build import surface, molecule, add_adsorbate
from ase.constraints import FixAtoms, FixedPlane, FixBondLength, ExpCellFilter
from ase.neb import NEB
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
from ase.visualize import view
from ase.build.rotate import minimize_rotation_and_translation
from ase.optimize import BFGS, LBFGS, FIRE
from ase.md import MDLogger
from ase.io import read, write, Trajectory
from ase.build import sort

from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode


estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL_PLUS_D3, model_version="v3.0.0")
calculator = ASECalculator(estimator)

1. Create slab from bulk structures

Rh is used as a catalyst in this example. The slab structure can be created from the bulk structure.

1-1 Bulk structure loading

The cif file of the bulk structure is downloaded from Materials Project and saved into the input folder.

Input cif file is from
A. Jain, S.P. Ong, G. Hautier, W. Chen, W.D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K.A. Persson (*=equal contributions)
The Materials Project: A materials genome approach to accelerating materials innovation APL Materials, 2013, 1(1), 011002.
doi:10.1063/1.4812323
[bibtex]
Licensed under CC BY 4.0
[2]:

bulk = read("../input/Rh_mp-74_conventional_standard.cif")
print("Number of atoms =", len(bulk))
print("Initial lattice constant =", bulk.cell.cellpar())


bulk.calc = calculator
opt = LBFGS(ExpCellFilter(bulk))
opt.run()
print ("Optimized lattice constant =", bulk.cell.cellpar())

Number of atoms = 4
Initial lattice constant = [ 3.843898  3.843898  3.843898 90.       90.       90.      ]
       Step     Time          Energy         fmax
*Force-consistent energies used in optimization.
LBFGS:    0 07:14:41      -25.513042*       3.6306
LBFGS:    1 07:14:42      -24.835918*      13.6274
LBFGS:    2 07:14:42      -25.584866*       0.6765
LBFGS:    3 07:14:42      -25.587189*       0.1175
LBFGS:    4 07:14:43      -25.587259*       0.0013
Optimized lattice constant = [ 3.793282    3.79328225  3.79328165 90.00000028 90.00000252 90.00001485]

The sort method in ASE sorts the order of atoms according to the atomic number.

[3]:

from pfcc_extras.visualize.view import view_ngl

bulk = bulk.repeat([2, 2, 2])
bulk = sort(bulk)
# Shift positions a bit to prevent surface to be cut at wrong place when `makesurface` is called
bulk.positions += [0.01, 0, 0]

view_ngl(bulk, representations=["ball+stick"])

1-2 Create slab structures

Create slab structures from bulk structures with any Miller index. We can specify the miller_indices=(x,y,z). The surface method is used in the makesurface to create a surface structure.

[5]:

from ase import Atoms
from ase.build import surface


def makesurface(
    atoms: Atoms, miller_indices=(1, 1, 1), layers=4, rep=[4, 4, 1]
) -> Atoms:
    s1 = surface(atoms, miller_indices, layers)
    s1.center(vacuum=10.0, axis=2)
    s1 = s1.repeat(rep)
    s1.set_positions(s1.get_positions() - [0, 0, min(s1.get_positions()[:, 2])])
    s1.pbc = True
    return s1

[6]:

slab = makesurface(bulk, miller_indices=(1, 1, 1), layers=2, rep=[1, 1, 1])
slab = sort(slab)
# adjust `positions` before `wrap`
slab.positions += [1, 1, 0]
slab.wrap()

view_ngl(slab, representations=["ball+stick"])

1-3 Check the z-coordinate of the created slab

Check the highest coordinates of the slab along the z axis (it is necessary when creating the adsorption structure). Check the coordinates of each layer of the slab (it is necessary to determine how many layers to fix).

[7]:

import plotly.express as px

# Check z_position of atoms
atoms = slab
df = pd.DataFrame({
    "x": atoms.positions[:, 0],
    "y": atoms.positions[:, 1],
    "z": atoms.positions[:, 2],
    "symbol": atoms.symbols,
})
display(df)

coord = "z"
df_sorted = df.sort_values(coord).reset_index().rename({"index": "atom_index"}, axis=1)
fig = px.scatter(df_sorted, x=df_sorted.index, y=coord, color="symbol", hover_data=["x", "y", "z", "atom_index"])
fig.show()
x y z symbol
0 1.007071 0.995918 4.458544e-08 Rh
1 10.394965 1.770218 2.190052e+00 Rh
2 6.371580 8.738923 2.190052e+00 Rh
3 2.348199 1.770218 2.190052e+00 Rh
4 14.418346 8.738923 2.190052e+00 Rh
... ... ... ... ...
59 5.030453 3.318820 6.570156e+00 Rh
60 1.007071 0.995918 6.570156e+00 Rh
61 5.030453 4.867420 4.380104e+00 Rh
62 6.371581 2.544519 4.380104e+00 Rh
63 7.712709 4.867420 4.380104e+00 Rh

64 rows × 4 columns