Intro to adsorption energies

Intro to adsorption energies#

To introduce OCP we start with using it to calculate adsorption energies for a simple, atomic adsorbate where we specify the site we want to the adsorption energy for. Conceptually, you do this like you would do it with density functional theory. You create a slab model for the surface, place an adsorbate on it as an initial guess, run a relaxation to get the lowest energy geometry, and then compute the adsorption energy using reference states for the adsorbate.

You do have to be careful in the details though. Some OCP model/checkpoint combinations return a total energy like density functional theory would, but some return an “adsorption energy” directly. You have to know which one you are using. In this example, the model we use returns an “adsorption energy”.

Intro to Adsorption energies#

Adsorption energies are always a reaction energy (an adsorbed species relative to some implied combination of reactants). There are many common schemes in the catalysis literature.

For example, you may want the adsorption energy of oxygen, and you might compute that from this reaction:

1/2 O2 + slab -> slab-O

DFT has known errors with the energy of a gas-phase O2 molecule, so it’s more common to compute this energy relative to a linear combination of H2O and H2. The suggested reference scheme for consistency with OC20 is a reaction

x CO + (x + y/2 - z) H2 + (z-x) H2O + w/2 N2 + * -> CxHyOzNw*

Here, x=y=w=0, z=1, so the reaction ends up as

-H2 + H2O + * -> O*

or alternatively,

H2O + * -> O* + H2

It is possible through thermodynamic cycles to compute other reactions. If we can look up rH1 below and compute rH2

H2 + 1/2 O2 -> H2O  re1 = -3.03 eV, from exp
H2O + * -> O* + H2  re2  # Get from UMA

Then, the adsorption energy for

1/2O2 + * -> O*

is just re1 + re2.

Based on https://atct.anl.gov/Thermochemical%20Data/version%201.118/species/?species_number=986, the formation energy of water is about -3.03 eV at standard state experimentally. You could also compute this using DFT, but you would probably get the wrong answer for this.

The first step is getting a checkpoint for the model we want to use. UMA is currently the state-of-the-art model and will provide total energy estimates at the RPBE level of theory if you use the “OC20” task.

Need to install fairchem-core or get UMA access or getting permissions/401 errors?

  1. Install the necessary packages using pip, uv etc

  1. Get access to any necessary huggingface gated models

    • Get and login to your Huggingface account

    • Request access to https://huggingface.co/facebook/UMA

    • Create a Huggingface token at https://huggingface.co/settings/tokens/ with the permission “Permissions: Read access to contents of all public gated repos you can access”

    • Add the token as an environment variable using huggingface-cli login or by setting the HF_TOKEN environment variable.

If you find your kernel is crashing, it probably means you have exceeded the allowed amount of memory. This checkpoint works fine in this example, but it may crash your kernel if you use it in the NRR example.

This next cell will automatically download the checkpoint from huggingface and load it.

from __future__ import annotations

from fairchem.core import FAIRChemCalculator, pretrained_mlip

predictor = pretrained_mlip.get_predict_unit("uma-s-1p1")
calc = FAIRChemCalculator(predictor, task_name="oc20")

WARNING:root:device was not explicitly set, using device='cuda'.

/home/runner/work/_tool/Python/3.12.12/x64/lib/python3.12/site-packages/pydantic/_internal/_generate_schema.py:2249: UnsupportedFieldAttributeWarning: The 'repr' attribute with value False was provided to the `Field()` function, which has no effect in the context it was used. 'repr' is field-specific metadata, and can only be attached to a model field using `Annotated` metadata or by assignment. This may have happened because an `Annotated` type alias using the `type` statement was used, or if the `Field()` function was attached to a single member of a union type.
  warnings.warn(
/home/runner/work/_tool/Python/3.12.12/x64/lib/python3.12/site-packages/pydantic/_internal/_generate_schema.py:2249: UnsupportedFieldAttributeWarning: The 'frozen' attribute with value True was provided to the `Field()` function, which has no effect in the context it was used. 'frozen' is field-specific metadata, and can only be attached to a model field using `Annotated` metadata or by assignment. This may have happened because an `Annotated` type alias using the `type` statement was used, or if the `Field()` function was attached to a single member of a union type.
  warnings.warn(

Next we can build a slab with an adsorbate on it. Here we use the ASE module to build a Pt slab. We use the experimental lattice constant that is the default. This can introduce some small errors with DFT since the lattice constant can differ by a few percent, and it is common to use DFT lattice constants. In this example, we do not constrain any layers.

from ase.build import add_adsorbate, fcc111
from ase.optimize import BFGS

# reference energies from a linear combination of H2O/N2/CO/H2!
atomic_reference_energies = {
    "H": -3.477,
    "N": -8.083,
    "O": -7.204,
    "C": -7.282,
}

re1 = -3.03

slab = fcc111("Pt", size=(2, 2, 5), vacuum=20.0)
slab.pbc = True

adslab = slab.copy()
add_adsorbate(adslab, "O", height=1.2, position="fcc")

slab.set_calculator(calc)
opt = BFGS(slab)
opt.run(fmax=0.05, steps=100)
slab_e = slab.get_potential_energy()

adslab.set_calculator(calc)
opt = BFGS(adslab)
opt.run(fmax=0.05, steps=100)
adslab_e = adslab.get_potential_energy()

# Energy for ((H2O-H2) + * -> *O) + (H2 + 1/2O2 -> H2) leads to 1/2O2 + * -> *O!
adslab_e - slab_e - atomic_reference_energies["O"] + re1

/tmp/ipykernel_3398/3752951811.py:17: FutureWarning: Please use atoms.calc = calc
  slab.set_calculator(calc)

      Step     Time          Energy          fmax
BFGS:    0 07:02:49     -104.695193        0.711414

BFGS:    1 07:02:49     -104.753188        0.608871

BFGS:    2 07:02:49     -104.906054        0.369514

BFGS:    3 07:02:49     -104.938753        0.439860

BFGS:    4 07:02:49     -105.016137        0.464029

BFGS:    5 07:02:49     -105.076559        0.356131

BFGS:    6 07:02:49     -105.112618        0.189549

BFGS:    7 07:02:50     -105.126757        0.045521

      Step     Time          Energy          fmax
BFGS:    0 07:02:50     -110.048787        1.757515

BFGS:    1 07:02:50     -110.230722        0.986430

/tmp/ipykernel_3398/3752951811.py:22: FutureWarning: Please use atoms.calc = calc
  adslab.set_calculator(calc)

BFGS:    2 07:02:50     -110.379365        0.731520

BFGS:    3 07:02:50     -110.430299        0.807331

BFGS:    4 07:02:50     -110.543682        0.692002

BFGS:    5 07:02:50     -110.617478        0.492741

BFGS:    6 07:02:50     -110.673464        0.666975

BFGS:    7 07:02:50     -110.721680        0.710481

BFGS:    8 07:02:50     -110.756904        0.443809

BFGS:    9 07:02:50     -110.769294        0.214548

BFGS:   10 07:02:51     -110.772451        0.091726

BFGS:   11 07:02:51     -110.772963        0.062470

BFGS:   12 07:02:51     -110.773256        0.046208

-1.4724996287652696

It is good practice to look at your geometries to make sure they are what you expect.

import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms

fig, axs = plt.subplots(1, 2)
plot_atoms(slab, axs[0])
plot_atoms(slab, axs[1], rotation=("-90x"))
axs[0].set_axis_off()
axs[1].set_axis_off()
../../_images/1e9975099b3d379bec3620d01733397931960485864b94438768f13e1551417d.png 
import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms

fig, axs = plt.subplots(1, 2)
plot_atoms(adslab, axs[0])
plot_atoms(adslab, axs[1], rotation=("-90x"))
axs[0].set_axis_off()
axs[1].set_axis_off()
../../_images/011ae04d7ed7e0cefdb42d07e652918fdcb49292eba19c886e2e255f3d81bddb.png

How did we do? We need a reference point. In the paper below, there is an atomic adsorption energy for O on Pt(111) of about -4.264 eV. This is for the reaction O + * -> O*. To convert this to the dissociative adsorption energy, we have to add the reaction:

1/2 O2 -> O   D = 2.58 eV (expt)

to get a comparable energy of about -1.68 eV. There is about ~0.2 eV difference (we predicted -1.47 eV above, and the reference comparison is -1.68 eV) to account for. The biggest difference is likely due to the differences in exchange-correlation functional. The reference data used the PBE functional, and eSCN was trained on RPBE data. To additional places where there are differences include:

  1. Difference in lattice constant

  2. The reference energy used for the experiment references. These can differ by up to 0.5 eV from comparable DFT calculations.

  3. How many layers are relaxed in the calculation

Some of these differences tend to be systematic, and you can calibrate and correct these, especially if you can augment these with your own DFT calculations.

See convergence study for some additional studies of factors that influence this number.

Exercises#

  1. Explore the effect of the lattice constant on the adsorption energy.

  2. Try different sites, including the bridge and top sites. Compare the energies, and inspect the resulting geometries.

Next steps#

In the next step, we consider some more complex adsorbates in nitrogen reduction, and how we can leverage OCP to automate the search for the most stable adsorbate geometry. See the next step.

Convergence study#

In Calculating adsorption energies we discussed some possible reasons we might see a discrepancy. Here we investigate some factors that impact the computed energies.

In this section, the energies refer to the reaction 1/2 O2 -> O*.

Effects of number of layers#

Slab thickness could be a factor. Here we relax the whole slab, and see by about 4 layers the energy is converged to ~0.02 eV.

for nlayers in [3, 4, 5, 6, 7, 8]:
    slab = fcc111("Pt", size=(2, 2, nlayers), vacuum=10.0)

    slab.pbc = True
    slab.set_calculator(calc)
    opt_slab = BFGS(slab, logfile=None)
    opt_slab.run(fmax=0.05, steps=100)
    slab_e = slab.get_potential_energy()

    adslab = slab.copy()
    add_adsorbate(adslab, "O", height=1.2, position="fcc")

    adslab.pbc = True
    adslab.set_calculator(calc)
    opt_adslab = BFGS(adslab, logfile=None)
    opt_adslab.run(fmax=0.05, steps=100)
    adslab_e = adslab.get_potential_energy()

    print(
        f"nlayers = {nlayers}: {adslab_e - slab_e - atomic_reference_energies['O'] + re1:1.2f} eV"
    )

/tmp/ipykernel_3398/338101817.py:5: FutureWarning: Please use atoms.calc = calc
  slab.set_calculator(calc)

/tmp/ipykernel_3398/338101817.py:14: FutureWarning: Please use atoms.calc = calc
  adslab.set_calculator(calc)

nlayers = 3: -1.58 eV

nlayers = 4: -1.47 eV

nlayers = 5: -1.47 eV

nlayers = 6: -1.46 eV

nlayers = 7: -1.47 eV

nlayers = 8: -1.47 eV

Effects of relaxation#

It is common to only relax a few layers, and constrain lower layers to bulk coordinates. We do that here. We only relax the adsorbate and the top layer.

This has a small effect (0.1 eV).

from ase.constraints import FixAtoms

for nlayers in [3, 4, 5, 6, 7, 8]:
    slab = fcc111("Pt", size=(2, 2, nlayers), vacuum=10.0)

    slab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in slab]))
    slab.pbc = True
    slab.set_calculator(calc)
    opt_slab = BFGS(slab, logfile=None)
    opt_slab.run(fmax=0.05, steps=100)
    slab_e = slab.get_potential_energy()

    adslab = slab.copy()
    add_adsorbate(adslab, "O", height=1.2, position="fcc")

    adslab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in adslab]))
    adslab.pbc = True
    adslab.set_calculator(calc)
    opt_adslab = BFGS(adslab, logfile=None)
    opt_adslab.run(fmax=0.05, steps=100)
    adslab_e = adslab.get_potential_energy()

    print(
        f"nlayers = {nlayers}: {adslab_e - slab_e - atomic_reference_energies['O'] + re1:1.2f} eV"
    )

/tmp/ipykernel_3398/1426773950.py:8: FutureWarning: Please use atoms.calc = calc
  slab.set_calculator(calc)

/tmp/ipykernel_3398/1426773950.py:18: FutureWarning: Please use atoms.calc = calc
  adslab.set_calculator(calc)

nlayers = 3: -1.47 eV

nlayers = 4: -1.35 eV

nlayers = 5: -1.35 eV

nlayers = 6: -1.34 eV

nlayers = 7: -1.35 eV

nlayers = 8: -1.35 eV

Unit cell size#

Coverage effects are quite noticeable with oxygen. Here we consider larger unit cells. This effect is large, and the results don’t look right, usually adsorption energies get more favorable at lower coverage, not less. This suggests fine-tuning could be important even at low coverages.

for size in [1, 2, 3, 4, 5]:

    slab = fcc111("Pt", size=(size, size, 5), vacuum=10.0)

    slab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in slab]))
    slab.pbc = True
    slab.set_calculator(calc)
    opt_slab = BFGS(slab, logfile=None)
    opt_slab.run(fmax=0.05, steps=100)
    slab_e = slab.get_potential_energy()

    adslab = slab.copy()
    add_adsorbate(adslab, "O", height=1.2, position="fcc")

    adslab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in adslab]))
    adslab.pbc = True
    adslab.set_calculator(calc)
    opt_adslab = BFGS(adslab, logfile=None)
    opt_adslab.run(fmax=0.05, steps=100)
    adslab_e = adslab.get_potential_energy()

    print(
        f"({size}x{size}): {adslab_e - slab_e - atomic_reference_energies['O'] + re1:1.2f} eV"
    )

/tmp/ipykernel_3398/3371624330.py:7: FutureWarning: Please use atoms.calc = calc
  slab.set_calculator(calc)

/tmp/ipykernel_3398/3371624330.py:17: FutureWarning: Please use atoms.calc = calc
  adslab.set_calculator(calc)

(1x1): -0.12 eV

(2x2): -1.35 eV

(3x3): -1.38 eV

(4x4): -1.41 eV

(5x5): -1.42 eV

Summary#

As with DFT, you should take care to see how these kinds of decisions affect your results, and determine if they would change any interpretations or not.

Contents