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?
Install the necessary packages using pip, uv etc
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'.
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_9508/3752951811.py:17: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
Step Time Energy fmax
BFGS: 0 17:32:02 -104.695194 0.709590
BFGS: 1 17:32:02 -104.753187 0.607748
BFGS: 2 17:32:02 -104.906057 0.369388
BFGS: 3 17:32:02 -104.938755 0.439737
BFGS: 4 17:32:03 -105.016140 0.464003
BFGS: 5 17:32:03 -105.076559 0.356059
BFGS: 6 17:32:03 -105.112621 0.189429
BFGS: 7 17:32:03 -105.126757 0.045376
Step Time Energy fmax
BFGS: 0 17:32:03 -110.048789 1.757515
BFGS: 1 17:32:03 -110.230720 0.986298
/tmp/ipykernel_9508/3752951811.py:22: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 17:32:03 -110.379353 0.731484
BFGS: 3 17:32:03 -110.430292 0.807280
BFGS: 4 17:32:03 -110.543671 0.692034
BFGS: 5 17:32:03 -110.617482 0.492732
BFGS: 6 17:32:03 -110.673461 0.666936
BFGS: 7 17:32:04 -110.721682 0.710502
BFGS: 8 17:32:04 -110.756902 0.443860
BFGS: 9 17:32:04 -110.769296 0.214557
BFGS: 10 17:32:04 -110.772453 0.091736
BFGS: 11 17:32:04 -110.772964 0.062463
BFGS: 12 17:32:04 -110.773258 0.046184
-1.4725015361139024
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()

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()

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:
Difference in lattice constant
The reference energy used for the experiment references. These can differ by up to 0.5 eV from comparable DFT calculations.
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#
Explore the effect of the lattice constant on the adsorption energy.
Try different sites, including the bridge and top sites. Compare the energies, and inspect the resulting geometries.
Trends in adsorption energies across metals.#
Xu, Z., & Kitchin, J. R. (2014). Probing the coverage dependence of site and adsorbate configurational correlations on (111) surfaces of late transition metals. J. Phys. Chem. C, 118(44), 25597–25602. http://dx.doi.org/10.1021/jp508805h
These are atomic adsorption energies:
O + * -> O*
We have to do some work to get comparable numbers from OCP
H2 + 1/2 O2 -> H2O re1 = -3.03 eV
H2O + * -> O* + H2 re2 # Get from UMA
O -> 1/2 O2 re3 = -2.58 eV
Then, the adsorption energy for
O + * -> O*
is just re1 + re2 + re3.
Here we just look at the fcc site on Pt. First, we get the data stored in the paper.
Next we get the structures and compute their energies. Some subtle points are that we have to account for stoichiometry, and normalize the adsorption energy by the number of oxygens.
First we get a reference energy from the paper (PBE, 0.25 ML O on Pt(111)).
import json
with open("energies.json") as f:
edata = json.load(f)
with open("structures.json") as f:
sdata = json.load(f)
edata["Pt"]["O"]["fcc"]["0.25"]
-4.263842000000002
Next, we load data from the SI to get the geometry to start from.
with open("structures.json") as f:
s = json.load(f)
sfcc = s["Pt"]["O"]["fcc"]["0.25"]
Next, we construct the atomic geometry, run the geometry optimization, and compute the energy.
re3 = -2.58 # O -> 1/2 O2 re3 = -2.58 eV
from ase import Atoms
adslab = Atoms(sfcc["symbols"], positions=sfcc["pos"], cell=sfcc["cell"], pbc=True)
# Grab just the metal surface atoms
slab = adslab[adslab.arrays["numbers"] == adslab.arrays["numbers"][0]]
adsorbates = adslab[~(adslab.arrays["numbers"] == adslab.arrays["numbers"][0])]
slab.set_calculator(calc)
opt = BFGS(slab)
opt.run(fmax=0.05, steps=100)
adslab.set_calculator(calc)
opt = BFGS(adslab)
opt.run(fmax=0.05, steps=100)
re2 = (
adslab.get_potential_energy()
- slab.get_potential_energy()
- sum([atomic_reference_energies[x] for x in adsorbates.get_chemical_symbols()])
)
nO = 0
for atom in adslab:
if atom.symbol == "O":
nO += 1
re2 += re1 + re3
print(re2 / nO)
Step Time Energy fmax
BFGS: 0 17:32:06 -82.862062 1.000262
BFGS: 1 17:32:07 -82.918113 0.752921
/tmp/ipykernel_9508/647904475.py:10: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
BFGS: 2 17:32:07 -83.011695 0.335877
BFGS: 3 17:32:07 -83.016496 0.303125
BFGS: 4 17:32:07 -83.024498 0.224966
BFGS: 5 17:32:07 -83.030036 0.152583
BFGS: 6 17:32:07 -83.033361 0.088643
BFGS: 7 17:32:07 -83.034469 0.065779
BFGS: 8 17:32:07 -83.035207 0.068696
BFGS: 9 17:32:07 -83.035564 0.046773
Step Time Energy fmax
BFGS: 0 17:32:07 -88.770343 0.312349
BFGS: 1 17:32:07 -88.773631 0.271819
/tmp/ipykernel_9508/647904475.py:14: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 17:32:08 -88.784475 0.104053
BFGS: 3 17:32:08 -88.786016 0.101495
BFGS: 4 17:32:08 -88.788548 0.106341
BFGS: 5 17:32:08 -88.790503 0.089103
BFGS: 6 17:32:08 -88.791969 0.096363
BFGS: 7 17:32:08 -88.792658 0.099427
BFGS: 8 17:32:08 -88.793177 0.077464
BFGS: 9 17:32:08 -88.793675 0.036099
-4.164112902385875
Site correlations#
This cell reproduces a portion of a figure in the paper. We compare oxygen adsorption energies in the fcc and hcp sites across metals and coverages. These adsorption energies are highly correlated with each other because the adsorption sites are so similar.
At higher coverages, the agreement is not as good. This is likely because the model is extrapolating and needs to be fine-tuned.
import time
from tqdm import tqdm
t0 = time.time()
data = {"fcc": [], "hcp": []}
refdata = {"fcc": [], "hcp": []}
for metal in ["Cu", "Ag", "Pd", "Pt", "Rh", "Ir"]:
print(metal)
for site in ["fcc", "hcp"]:
for adsorbate in ["O"]:
for coverage in tqdm(["0.25"]):
entry = s[metal][adsorbate][site][coverage]
adslab = Atoms(
entry["symbols"],
positions=entry["pos"],
cell=entry["cell"],
pbc=True,
)
# Grab just the metal surface atoms
adsorbates = adslab[
~(adslab.arrays["numbers"] == adslab.arrays["numbers"][0])
]
slab = adslab[adslab.arrays["numbers"] == adslab.arrays["numbers"][0]]
slab.set_calculator(calc)
opt = BFGS(slab)
opt.run(fmax=0.05, steps=100)
adslab.set_calculator(calc)
opt = BFGS(adslab)
opt.run(fmax=0.05, steps=100)
re2 = (
adslab.get_potential_energy()
- slab.get_potential_energy()
- sum(
[
atomic_reference_energies[x]
for x in adsorbates.get_chemical_symbols()
]
)
)
nO = 0
for atom in adslab:
if atom.symbol == "O":
nO += 1
re2 += re1 + re3
data[site] += [re2 / nO]
refdata[site] += [edata[metal][adsorbate][site][coverage]]
f"Elapsed time = {time.time() - t0} seconds"
Cu
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:08 -48.910948 0.646351
BFGS: 1 17:32:08 -48.934031 0.539302
/tmp/ipykernel_9508/1356342052.py:33: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
BFGS: 2 17:32:09 -49.000045 0.280109
BFGS: 3 17:32:09 -49.002375 0.251316
BFGS: 4 17:32:09 -49.008803 0.158977
BFGS: 5 17:32:09 -49.013980 0.134498
BFGS: 6 17:32:09 -49.016868 0.071553
BFGS: 7 17:32:09 -49.017633 0.051239
BFGS: 8 17:32:09 -49.018012 0.047571
Step Time Energy fmax
BFGS: 0 17:32:09 -55.204221 0.339047
BFGS: 1 17:32:09 -55.206662 0.267936
/tmp/ipykernel_9508/1356342052.py:37: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 17:32:09 -55.213714 0.160209
BFGS: 3 17:32:09 -55.215870 0.166138
BFGS: 4 17:32:10 -55.219684 0.116282
BFGS: 5 17:32:10 -55.221818 0.079620
BFGS: 6 17:32:10 -55.223363 0.075240
BFGS: 7 17:32:10 -55.224369 0.090440
BFGS: 8 17:32:10 -55.225523 0.098844
BFGS: 9 17:32:10 -55.226532 0.068283
BFGS: 10 17:32:10 -55.226978 0.043064
100%|██████████| 1/1 [00:01<00:00, 1.91s/it]
100%|██████████| 1/1 [00:01<00:00, 1.91s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:10 -48.936503 0.552695
BFGS: 1 17:32:10 -48.954084 0.470533
BFGS: 2 17:32:10 -49.008577 0.211009
BFGS: 3 17:32:11 -49.009768 0.197680
BFGS: 4 17:32:11 -49.017746 0.041107
Step Time Energy fmax
BFGS: 0 17:32:11 -55.119069 0.305598
BFGS: 1 17:32:11 -55.120918 0.239016
BFGS: 2 17:32:11 -55.125707 0.129255
BFGS: 3 17:32:11 -55.127340 0.147662
BFGS: 4 17:32:11 -55.130401 0.128550
BFGS: 5 17:32:11 -55.132308 0.103047
BFGS: 6 17:32:11 -55.133380 0.051934
BFGS: 7 17:32:11 -55.133849 0.056110
BFGS: 8 17:32:11 -55.134356 0.071279
BFGS: 9 17:32:12 -55.135042 0.069427
BFGS: 10 17:32:12 -55.135573 0.039950
100%|██████████| 1/1 [00:01<00:00, 1.54s/it]
100%|██████████| 1/1 [00:01<00:00, 1.54s/it]
Ag
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:12 -33.034314 0.634820
BFGS: 1 17:32:12 -33.053470 0.551048
BFGS: 2 17:32:12 -33.120905 0.172022
BFGS: 3 17:32:12 -33.122777 0.157332
BFGS: 4 17:32:12 -33.124180 0.143388
BFGS: 5 17:32:12 -33.126944 0.112467
BFGS: 6 17:32:12 -33.130770 0.110662
BFGS: 7 17:32:12 -33.133994 0.067520
BFGS: 8 17:32:13 -33.134834 0.034877
Step Time Energy fmax
BFGS: 0 17:32:13 -38.173508 0.098015
BFGS: 1 17:32:13 -38.174435 0.094366
BFGS: 2 17:32:13 -38.184910 0.084181
BFGS: 3 17:32:13 -38.185730 0.085524
BFGS: 4 17:32:13 -38.189101 0.085540
BFGS: 5 17:32:13 -38.190997 0.067880
BFGS: 6 17:32:13 -38.192282 0.087748
BFGS: 7 17:32:13 -38.193201 0.076336
BFGS: 8 17:32:13 -38.193675 0.041787
100%|██████████| 1/1 [00:01<00:00, 1.72s/it]
100%|██████████| 1/1 [00:01<00:00, 1.72s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:14 -33.056750 0.557366
BFGS: 1 17:32:14 -33.071228 0.489413
BFGS: 2 17:32:14 -33.126399 0.123838
BFGS: 3 17:32:14 -33.127308 0.112963
BFGS: 4 17:32:14 -33.128516 0.097381
BFGS: 5 17:32:14 -33.130259 0.079550
BFGS: 6 17:32:14 -33.133159 0.076352
BFGS: 7 17:32:14 -33.134742 0.040729
Step Time Energy fmax
BFGS: 0 17:32:14 -38.102025 0.092717
BFGS: 1 17:32:14 -38.102856 0.088832
BFGS: 2 17:32:14 -38.111091 0.134762
BFGS: 3 17:32:15 -38.112013 0.120354
BFGS: 4 17:32:15 -38.115299 0.064067
BFGS: 5 17:32:15 -38.116158 0.056494
BFGS: 6 17:32:15 -38.116715 0.058457
BFGS: 7 17:32:15 -38.117585 0.073179
BFGS: 8 17:32:15 -38.118456 0.057555
BFGS: 9 17:32:15 -38.118835 0.026191
100%|██████████| 1/1 [00:01<00:00, 1.72s/it]
100%|██████████| 1/1 [00:01<00:00, 1.72s/it]
Pd
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:15 -70.165062 0.631446
BFGS: 1 17:32:15 -70.189910 0.508128
BFGS: 2 17:32:15 -70.241649 0.191936
BFGS: 3 17:32:16 -70.243346 0.182565
BFGS: 4 17:32:16 -70.248577 0.144559
BFGS: 5 17:32:16 -70.251930 0.115930
BFGS: 6 17:32:16 -70.254308 0.074845
BFGS: 7 17:32:16 -70.255205 0.059459
BFGS: 8 17:32:16 -70.255944 0.041104
Step Time Energy fmax
BFGS: 0 17:32:16 -76.130241 0.213656
BFGS: 1 17:32:16 -76.133242 0.190340
BFGS: 2 17:32:16 -76.148719 0.169654
BFGS: 3 17:32:16 -76.150811 0.152253
BFGS: 4 17:32:16 -76.154315 0.130500
BFGS: 5 17:32:17 -76.157061 0.106082
BFGS: 6 17:32:17 -76.159304 0.088152
BFGS: 7 17:32:17 -76.160431 0.080482
BFGS: 8 17:32:17 -76.161130 0.064463
BFGS: 9 17:32:17 -76.161613 0.042968
100%|██████████| 1/1 [00:01<00:00, 1.81s/it]
100%|██████████| 1/1 [00:01<00:00, 1.81s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:17 -70.197711 0.454636
BFGS: 1 17:32:17 -70.211775 0.373891
BFGS: 2 17:32:17 -70.245976 0.178188
BFGS: 3 17:32:17 -70.247162 0.164584
BFGS: 4 17:32:17 -70.253354 0.078572
BFGS: 5 17:32:18 -70.254513 0.075442
BFGS: 6 17:32:18 -70.255273 0.058065
BFGS: 7 17:32:18 -70.255846 0.036861
Step Time Energy fmax
BFGS: 0 17:32:18 -75.892387 0.169440
BFGS: 1 17:32:18 -75.895234 0.150801
BFGS: 2 17:32:18 -75.906635 0.144089
BFGS: 3 17:32:18 -75.908426 0.135658
BFGS: 4 17:32:18 -75.912424 0.128986
BFGS: 5 17:32:18 -75.915493 0.108571
BFGS: 6 17:32:18 -75.917927 0.075205
BFGS: 7 17:32:18 -75.918956 0.067266
BFGS: 8 17:32:19 -75.919478 0.044605
100%|██████████| 1/1 [00:01<00:00, 1.63s/it]
100%|██████████| 1/1 [00:01<00:00, 1.63s/it]
Pt
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:19 -82.862062 1.000262
BFGS: 1 17:32:19 -82.918113 0.752921
BFGS: 2 17:32:19 -83.011695 0.335878
BFGS: 3 17:32:19 -83.016496 0.303125
BFGS: 4 17:32:19 -83.024498 0.224966
BFGS: 5 17:32:19 -83.030036 0.152581
BFGS: 6 17:32:19 -83.033361 0.088643
BFGS: 7 17:32:19 -83.034469 0.065780
BFGS: 8 17:32:19 -83.035207 0.068697
BFGS: 9 17:32:20 -83.035563 0.046772
Step Time Energy fmax
BFGS: 0 17:32:20 -88.770343 0.312348
BFGS: 1 17:32:20 -88.773631 0.271819
BFGS: 2 17:32:20 -88.784475 0.104053
BFGS: 3 17:32:20 -88.786016 0.101495
BFGS: 4 17:32:20 -88.788549 0.106339
BFGS: 5 17:32:20 -88.790503 0.089105
BFGS: 6 17:32:20 -88.791969 0.096362
BFGS: 7 17:32:20 -88.792657 0.099420
BFGS: 8 17:32:20 -88.793177 0.077455
BFGS: 9 17:32:20 -88.793677 0.036096
100%|██████████| 1/1 [00:01<00:00, 1.90s/it]
100%|██████████| 1/1 [00:01<00:00, 1.90s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:21 -82.946995 0.675644
BFGS: 1 17:32:21 -82.973148 0.546707
BFGS: 2 17:32:21 -83.026506 0.203008
BFGS: 3 17:32:21 -83.028072 0.183984
BFGS: 4 17:32:21 -83.033698 0.083433
BFGS: 5 17:32:21 -83.034396 0.061089
BFGS: 6 17:32:21 -83.035121 0.048422
Step Time Energy fmax
BFGS: 0 17:32:21 -88.350562 0.234982
BFGS: 1 17:32:21 -88.353606 0.203898
BFGS: 2 17:32:21 -88.362973 0.117049
BFGS: 3 17:32:22 -88.364340 0.121262
BFGS: 4 17:32:22 -88.368488 0.090168
BFGS: 5 17:32:22 -88.370473 0.073296
BFGS: 6 17:32:22 -88.371781 0.075107
BFGS: 7 17:32:22 -88.372314 0.079561
BFGS: 8 17:32:22 -88.372727 0.061506
BFGS: 9 17:32:22 -88.373020 0.027284
100%|██████████| 1/1 [00:01<00:00, 1.83s/it]
100%|██████████| 1/1 [00:01<00:00, 1.83s/it]
Rh
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:22 -100.192566 0.688933
BFGS: 1 17:32:23 -100.218239 0.600213
BFGS: 2 17:32:23 -100.286835 0.174749
BFGS: 3 17:32:23 -100.289767 0.129383
BFGS: 4 17:32:23 -100.298132 0.066656
BFGS: 5 17:32:23 -100.299265 0.053487
BFGS: 6 17:32:23 -100.299778 0.045081
Step Time Energy fmax
BFGS: 0 17:32:23 -106.921794 0.203149
BFGS: 1 17:32:23 -106.924541 0.168768
BFGS: 2 17:32:23 -106.930976 0.052117
BFGS: 3 17:32:23 -106.931583 0.045576
100%|██████████| 1/1 [00:01<00:00, 1.09s/it]
100%|██████████| 1/1 [00:01<00:00, 1.09s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:24 -100.171741 0.714774
BFGS: 1 17:32:24 -100.203866 0.626030
BFGS: 2 17:32:24 -100.287845 0.263179
BFGS: 3 17:32:24 -100.290543 0.196018
BFGS: 4 17:32:24 -100.297166 0.080393
BFGS: 5 17:32:24 -100.299033 0.056462
BFGS: 6 17:32:24 -100.299773 0.043867
Step Time Energy fmax
BFGS: 0 17:32:24 -106.876840 0.168839
BFGS: 1 17:32:24 -106.879420 0.139489
BFGS: 2 17:32:24 -106.885548 0.071259
BFGS: 3 17:32:24 -106.886238 0.064287
BFGS: 4 17:32:25 -106.886579 0.054985
BFGS: 5 17:32:25 -106.886842 0.043047
100%|██████████| 1/1 [00:01<00:00, 1.27s/it]
100%|██████████| 1/1 [00:01<00:00, 1.27s/it]
Ir
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:25 -124.235610 1.189926
BFGS: 1 17:32:25 -124.310841 0.937185
BFGS: 2 17:32:25 -124.422189 0.140816
BFGS: 3 17:32:25 -124.425271 0.133724
BFGS: 4 17:32:25 -124.429598 0.060624
BFGS: 5 17:32:25 -124.430258 0.038669
Step Time Energy fmax
BFGS: 0 17:32:25 -130.605623 0.371646
BFGS: 1 17:32:25 -130.616146 0.265891
BFGS: 2 17:32:26 -130.630317 0.110764
BFGS: 3 17:32:26 -130.633614 0.094354
BFGS: 4 17:32:26 -130.634277 0.075691
BFGS: 5 17:32:26 -130.634990 0.065695
BFGS: 6 17:32:26 -130.636197 0.083044
BFGS: 7 17:32:26 -130.636656 0.075131
BFGS: 8 17:32:26 -130.637048 0.047413
100%|██████████| 1/1 [00:01<00:00, 1.45s/it]
100%|██████████| 1/1 [00:01<00:00, 1.45s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 17:32:26 -124.227592 1.165685
BFGS: 1 17:32:26 -124.311476 0.915619
BFGS: 2 17:32:26 -124.422556 0.165254
BFGS: 3 17:32:27 -124.425480 0.200142
BFGS: 4 17:32:27 -124.428759 0.168918
BFGS: 5 17:32:27 -124.430254 0.088490
BFGS: 6 17:32:27 -124.430560 0.041839
Step Time Energy fmax
BFGS: 0 17:32:27 -130.496158 0.405482
BFGS: 1 17:32:27 -130.507202 0.266989
BFGS: 2 17:32:27 -130.518616 0.126412
BFGS: 3 17:32:27 -130.522315 0.080272
BFGS: 4 17:32:27 -130.523180 0.066540
BFGS: 5 17:32:27 -130.523452 0.056600
BFGS: 6 17:32:27 -130.524509 0.055664
BFGS: 7 17:32:28 -130.524677 0.042715
100%|██████████| 1/1 [00:01<00:00, 1.45s/it]
100%|██████████| 1/1 [00:01<00:00, 1.45s/it]
'Elapsed time = 19.344956636428833 seconds'
First, we compare the computed data and reference data. There is a systematic difference of about 0.5 eV due to the difference between RPBE and PBE functionals, and other subtle differences like lattice constant differences and reference energy differences. This is pretty typical, and an expected deviation.
plt.plot(refdata["fcc"], data["fcc"], "r.", label="fcc")
plt.plot(refdata["hcp"], data["hcp"], "b.", label="hcp")
plt.plot([-5.5, -3.5], [-5.5, -3.5], "k-")
plt.xlabel("Ref. data (DFT)")
plt.ylabel("UMA-OC20 prediction");

Next we compare the correlation between the hcp and fcc sites. Here we see the same trends. The data falls below the parity line because the hcp sites tend to be a little weaker binding than the fcc sites.
plt.plot(refdata["hcp"], refdata["fcc"], "r.")
plt.plot(data["hcp"], data["fcc"], ".")
plt.plot([-6, -1], [-6, -1], "k-")
plt.xlabel("$H_{ads, hcp}$")
plt.ylabel("$H_{ads, fcc}$")
plt.legend(["DFT (PBE)", "UMA-OC20"]);

Exercises#
You can also explore a few other adsorbates: C, H, N.
Explore the higher coverages. The deviations from the reference data are expected to be higher, but relative differences tend to be better. You probably need fine tuning to improve this performance. This data set doesn’t have forces though, so it isn’t practical to do it here.
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_9508/338101817.py:5: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9508/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_9508/1426773950.py:8: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9508/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_9508/3371624330.py:7: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9508/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.