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_5220/3752951811.py:17: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
Step Time Energy fmax
BFGS: 0 00:25:52 -104.695193 0.709590
BFGS: 1 00:25:53 -104.753187 0.607748
BFGS: 2 00:25:53 -104.906057 0.369388
BFGS: 3 00:25:53 -104.938755 0.439737
BFGS: 4 00:25:53 -105.016144 0.463999
BFGS: 5 00:25:53 -105.076556 0.356078
BFGS: 6 00:25:53 -105.112622 0.189419
BFGS: 7 00:25:53 -105.126757 0.045374
Step Time Energy fmax
BFGS: 0 00:25:53 -110.048788 1.757515
BFGS: 1 00:25:53 -110.230720 0.986298
/tmp/ipykernel_5220/3752951811.py:22: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 00:25:53 -110.379353 0.731484
BFGS: 3 00:25:54 -110.430292 0.807280
BFGS: 4 00:25:54 -110.543670 0.692022
BFGS: 5 00:25:54 -110.617482 0.492731
BFGS: 6 00:25:54 -110.673463 0.666971
BFGS: 7 00:25:54 -110.721677 0.710499
BFGS: 8 00:25:54 -110.756906 0.443841
BFGS: 9 00:25:54 -110.769296 0.214536
BFGS: 10 00:25:54 -110.772453 0.091731
BFGS: 11 00:25:54 -110.772965 0.062450
BFGS: 12 00:25:54 -110.773259 0.046179
-1.4725024897882188
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 00:25:57 -82.862062 1.000262
BFGS: 1 00:25:57 -82.918113 0.752921
/tmp/ipykernel_5220/647904475.py:10: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
BFGS: 2 00:25:57 -83.011695 0.335877
BFGS: 3 00:25:57 -83.016496 0.303125
BFGS: 4 00:25:57 -83.024498 0.224966
BFGS: 5 00:25:57 -83.030036 0.152582
BFGS: 6 00:25:57 -83.033360 0.088642
BFGS: 7 00:25:57 -83.034469 0.065780
BFGS: 8 00:25:58 -83.035207 0.068697
BFGS: 9 00:25:58 -83.035564 0.046775
Step Time Energy fmax
BFGS: 0 00:25:58 -88.770343 0.312349
BFGS: 1 00:25:58 -88.773629 0.271819
/tmp/ipykernel_5220/647904475.py:14: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 00:25:58 -88.784474 0.104054
BFGS: 3 00:25:58 -88.786016 0.101495
BFGS: 4 00:25:58 -88.788548 0.106339
BFGS: 5 00:25:58 -88.790503 0.089104
BFGS: 6 00:25:58 -88.791969 0.096346
BFGS: 7 00:25:58 -88.792657 0.099425
BFGS: 8 00:25:58 -88.793177 0.077477
BFGS: 9 00:25:59 -88.793677 0.036098
-4.164113379223033
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 00:25:59 -48.910949 0.646351
BFGS: 1 00:25:59 -48.934031 0.539302
/tmp/ipykernel_5220/1356342052.py:33: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
BFGS: 2 00:25:59 -49.000045 0.280117
BFGS: 3 00:25:59 -49.002375 0.251316
BFGS: 4 00:25:59 -49.008803 0.158983
BFGS: 5 00:25:59 -49.013981 0.134494
BFGS: 6 00:25:59 -49.016868 0.071556
BFGS: 7 00:25:59 -49.017633 0.051237
BFGS: 8 00:25:59 -49.018012 0.047572
Step Time Energy fmax
BFGS: 0 00:26:00 -55.204221 0.339047
BFGS: 1 00:26:00 -55.206662 0.267936
/tmp/ipykernel_5220/1356342052.py:37: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 2 00:26:00 -55.213715 0.160209
BFGS: 3 00:26:00 -55.215870 0.166138
BFGS: 4 00:26:00 -55.219684 0.116282
BFGS: 5 00:26:00 -55.221818 0.079620
BFGS: 6 00:26:00 -55.223364 0.075241
BFGS: 7 00:26:00 -55.224369 0.090439
BFGS: 8 00:26:00 -55.225523 0.098844
BFGS: 9 00:26:00 -55.226532 0.068284
BFGS: 10 00:26:00 -55.226978 0.043064
100%|██████████| 1/1 [00:01<00:00, 1.93s/it]
100%|██████████| 1/1 [00:01<00:00, 1.93s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:01 -48.936503 0.552695
BFGS: 1 00:26:01 -48.954084 0.470533
BFGS: 2 00:26:01 -49.008577 0.211009
BFGS: 3 00:26:01 -49.009768 0.197681
BFGS: 4 00:26:01 -49.017746 0.041106
Step Time Energy fmax
BFGS: 0 00:26:01 -55.119070 0.305598
BFGS: 1 00:26:01 -55.120918 0.239016
BFGS: 2 00:26:01 -55.125707 0.129256
BFGS: 3 00:26:01 -55.127340 0.147662
BFGS: 4 00:26:01 -55.130401 0.128546
BFGS: 5 00:26:02 -55.132308 0.103050
BFGS: 6 00:26:02 -55.133380 0.051931
BFGS: 7 00:26:02 -55.133850 0.056103
BFGS: 8 00:26:02 -55.134357 0.071276
BFGS: 9 00:26:02 -55.135042 0.069433
BFGS: 10 00:26:02 -55.135574 0.039954
100%|██████████| 1/1 [00:01<00:00, 1.55s/it]
100%|██████████| 1/1 [00:01<00:00, 1.55s/it]
Ag
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:02 -33.034314 0.634820
BFGS: 1 00:26:02 -33.053470 0.551048
BFGS: 2 00:26:02 -33.120905 0.172022
BFGS: 3 00:26:03 -33.122777 0.157332
BFGS: 4 00:26:03 -33.124180 0.143388
BFGS: 5 00:26:03 -33.126944 0.112465
BFGS: 6 00:26:03 -33.130771 0.110663
BFGS: 7 00:26:03 -33.133994 0.067524
BFGS: 8 00:26:03 -33.134834 0.034877
Step Time Energy fmax
BFGS: 0 00:26:03 -38.173509 0.098015
BFGS: 1 00:26:03 -38.174435 0.094366
BFGS: 2 00:26:03 -38.184910 0.084180
BFGS: 3 00:26:03 -38.185730 0.085524
BFGS: 4 00:26:03 -38.189101 0.085542
BFGS: 5 00:26:04 -38.190996 0.067885
BFGS: 6 00:26:04 -38.192281 0.087741
BFGS: 7 00:26:04 -38.193202 0.076349
BFGS: 8 00:26:04 -38.193673 0.041782
100%|██████████| 1/1 [00:01<00:00, 1.74s/it]
100%|██████████| 1/1 [00:01<00:00, 1.74s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:04 -33.056750 0.557366
BFGS: 1 00:26:04 -33.071228 0.489413
BFGS: 2 00:26:04 -33.126399 0.123838
BFGS: 3 00:26:04 -33.127308 0.112962
BFGS: 4 00:26:04 -33.128516 0.097381
BFGS: 5 00:26:04 -33.130259 0.079550
BFGS: 6 00:26:05 -33.133159 0.076352
BFGS: 7 00:26:05 -33.134743 0.040729
Step Time Energy fmax
BFGS: 0 00:26:05 -38.102025 0.092717
BFGS: 1 00:26:05 -38.102856 0.088832
BFGS: 2 00:26:05 -38.111092 0.134761
BFGS: 3 00:26:05 -38.112013 0.120351
BFGS: 4 00:26:05 -38.115299 0.064070
BFGS: 5 00:26:05 -38.116157 0.056498
BFGS: 6 00:26:05 -38.116716 0.058454
BFGS: 7 00:26:05 -38.117585 0.073179
BFGS: 8 00:26:05 -38.118456 0.057569
BFGS: 9 00:26:06 -38.118835 0.026185
100%|██████████| 1/1 [00:01<00:00, 1.74s/it]
100%|██████████| 1/1 [00:01<00:00, 1.74s/it]
Pd
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:06 -70.165062 0.631446
BFGS: 1 00:26:06 -70.189910 0.508128
BFGS: 2 00:26:06 -70.241648 0.191936
BFGS: 3 00:26:06 -70.243346 0.182565
BFGS: 4 00:26:06 -70.248577 0.144562
BFGS: 5 00:26:06 -70.251931 0.115929
BFGS: 6 00:26:06 -70.254308 0.074838
BFGS: 7 00:26:06 -70.255205 0.059459
BFGS: 8 00:26:06 -70.255945 0.041102
Step Time Energy fmax
BFGS: 0 00:26:07 -76.130241 0.213656
BFGS: 1 00:26:07 -76.133242 0.190340
BFGS: 2 00:26:07 -76.148720 0.169651
BFGS: 3 00:26:07 -76.150811 0.152250
BFGS: 4 00:26:07 -76.154316 0.130500
BFGS: 5 00:26:07 -76.157060 0.106088
BFGS: 6 00:26:07 -76.159304 0.088151
BFGS: 7 00:26:07 -76.160431 0.080476
BFGS: 8 00:26:07 -76.161129 0.064477
BFGS: 9 00:26:07 -76.161614 0.042955
100%|██████████| 1/1 [00:01<00:00, 1.83s/it]
100%|██████████| 1/1 [00:01<00:00, 1.84s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:08 -70.197712 0.454636
BFGS: 1 00:26:08 -70.211776 0.373891
BFGS: 2 00:26:08 -70.245976 0.178188
BFGS: 3 00:26:08 -70.247162 0.164584
BFGS: 4 00:26:08 -70.253355 0.078574
BFGS: 5 00:26:08 -70.254514 0.075442
BFGS: 6 00:26:08 -70.255273 0.058063
BFGS: 7 00:26:08 -70.255845 0.036861
Step Time Energy fmax
BFGS: 0 00:26:08 -75.892388 0.169440
BFGS: 1 00:26:08 -75.895234 0.150801
BFGS: 2 00:26:08 -75.906635 0.144091
BFGS: 3 00:26:09 -75.908426 0.135659
BFGS: 4 00:26:09 -75.912426 0.128988
BFGS: 5 00:26:09 -75.915494 0.108566
BFGS: 6 00:26:09 -75.917927 0.075207
BFGS: 7 00:26:09 -75.918955 0.067267
BFGS: 8 00:26:09 -75.919479 0.044597
100%|██████████| 1/1 [00:01<00:00, 1.85s/it]
100%|██████████| 1/1 [00:01<00:00, 1.85s/it]
Pt
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:09 -82.862062 1.000262
BFGS: 1 00:26:09 -82.918114 0.752921
BFGS: 2 00:26:10 -83.011695 0.335878
BFGS: 3 00:26:10 -83.016496 0.303125
BFGS: 4 00:26:10 -83.024498 0.224962
BFGS: 5 00:26:10 -83.030036 0.152581
BFGS: 6 00:26:10 -83.033361 0.088645
BFGS: 7 00:26:10 -83.034470 0.065780
BFGS: 8 00:26:10 -83.035207 0.068698
BFGS: 9 00:26:10 -83.035563 0.046775
Step Time Energy fmax
BFGS: 0 00:26:10 -88.770343 0.312348
BFGS: 1 00:26:10 -88.773630 0.271819
BFGS: 2 00:26:11 -88.784475 0.104064
BFGS: 3 00:26:11 -88.786016 0.101495
BFGS: 4 00:26:11 -88.788548 0.106339
BFGS: 5 00:26:11 -88.790503 0.089101
BFGS: 6 00:26:11 -88.791968 0.096347
BFGS: 7 00:26:11 -88.792658 0.099429
BFGS: 8 00:26:11 -88.793177 0.077473
BFGS: 9 00:26:11 -88.793677 0.036099
100%|██████████| 1/1 [00:01<00:00, 1.93s/it]
100%|██████████| 1/1 [00:01<00:00, 1.93s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:11 -82.946997 0.675644
BFGS: 1 00:26:11 -82.973148 0.546707
BFGS: 2 00:26:12 -83.026505 0.203008
BFGS: 3 00:26:12 -83.028072 0.183984
BFGS: 4 00:26:12 -83.033698 0.083431
BFGS: 5 00:26:12 -83.034396 0.061089
BFGS: 6 00:26:12 -83.035121 0.048421
Step Time Energy fmax
BFGS: 0 00:26:12 -88.350562 0.234982
BFGS: 1 00:26:12 -88.353606 0.203899
BFGS: 2 00:26:12 -88.362973 0.117048
BFGS: 3 00:26:12 -88.364342 0.121272
BFGS: 4 00:26:12 -88.368490 0.090166
BFGS: 5 00:26:12 -88.370474 0.073296
BFGS: 6 00:26:13 -88.371781 0.075100
BFGS: 7 00:26:13 -88.372314 0.079555
BFGS: 8 00:26:13 -88.372729 0.061518
BFGS: 9 00:26:13 -88.373020 0.027294
100%|██████████| 1/1 [00:01<00:00, 1.65s/it]
100%|██████████| 1/1 [00:01<00:00, 1.65s/it]
Rh
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:13 -100.192566 0.688933
BFGS: 1 00:26:13 -100.218239 0.600213
BFGS: 2 00:26:13 -100.286835 0.174749
BFGS: 3 00:26:13 -100.289767 0.129383
BFGS: 4 00:26:13 -100.298130 0.066657
BFGS: 5 00:26:13 -100.299265 0.053489
BFGS: 6 00:26:14 -100.299778 0.045085
Step Time Energy fmax
BFGS: 0 00:26:14 -106.921793 0.203149
BFGS: 1 00:26:14 -106.924541 0.168768
BFGS: 2 00:26:14 -106.930976 0.052117
BFGS: 3 00:26:14 -106.931583 0.045576
100%|██████████| 1/1 [00:01<00:00, 1.10s/it]
100%|██████████| 1/1 [00:01<00:00, 1.10s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:14 -100.171741 0.714774
BFGS: 1 00:26:14 -100.203867 0.626030
BFGS: 2 00:26:14 -100.287846 0.263179
BFGS: 3 00:26:14 -100.290544 0.196018
BFGS: 4 00:26:14 -100.297167 0.080393
BFGS: 5 00:26:15 -100.299031 0.056461
BFGS: 6 00:26:15 -100.299772 0.043867
Step Time Energy fmax
BFGS: 0 00:26:15 -106.876840 0.168839
BFGS: 1 00:26:15 -106.879420 0.139489
BFGS: 2 00:26:15 -106.885548 0.071258
BFGS: 3 00:26:15 -106.886238 0.064287
BFGS: 4 00:26:15 -106.886579 0.054985
BFGS: 5 00:26:15 -106.886842 0.043047
100%|██████████| 1/1 [00:01<00:00, 1.28s/it]
100%|██████████| 1/1 [00:01<00:00, 1.28s/it]
Ir
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:15 -124.235610 1.189926
BFGS: 1 00:26:15 -124.310840 0.937185
BFGS: 2 00:26:16 -124.422189 0.140814
BFGS: 3 00:26:16 -124.425270 0.133727
BFGS: 4 00:26:16 -124.429597 0.060625
BFGS: 5 00:26:16 -124.430259 0.038680
Step Time Energy fmax
BFGS: 0 00:26:16 -130.605623 0.371646
BFGS: 1 00:26:16 -130.616146 0.265891
BFGS: 2 00:26:16 -130.630315 0.110766
BFGS: 3 00:26:16 -130.633614 0.094357
BFGS: 4 00:26:16 -130.634276 0.075691
BFGS: 5 00:26:16 -130.634989 0.065703
BFGS: 6 00:26:16 -130.636198 0.083030
BFGS: 7 00:26:17 -130.636657 0.075134
BFGS: 8 00:26:17 -130.637048 0.047437
100%|██████████| 1/1 [00:01<00:00, 1.47s/it]
100%|██████████| 1/1 [00:01<00:00, 1.47s/it]
0%| | 0/1 [00:00<?, ?it/s]
Step Time Energy fmax
BFGS: 0 00:26:17 -124.227592 1.165685
BFGS: 1 00:26:17 -124.311475 0.915619
BFGS: 2 00:26:17 -124.422556 0.165254
BFGS: 3 00:26:17 -124.425479 0.200141
BFGS: 4 00:26:17 -124.428759 0.168918
BFGS: 5 00:26:17 -124.430254 0.088490
BFGS: 6 00:26:17 -124.430560 0.041839
Step Time Energy fmax
BFGS: 0 00:26:17 -130.496159 0.405482
BFGS: 1 00:26:18 -130.507201 0.266989
BFGS: 2 00:26:18 -130.518615 0.126413
BFGS: 3 00:26:18 -130.522315 0.080271
BFGS: 4 00:26:18 -130.523181 0.066540
BFGS: 5 00:26:18 -130.523452 0.056600
BFGS: 6 00:26:18 -130.524510 0.055660
BFGS: 7 00:26:18 -130.524677 0.042718
100%|██████████| 1/1 [00:01<00:00, 1.47s/it]
100%|██████████| 1/1 [00:01<00:00, 1.47s/it]
'Elapsed time = 19.582616090774536 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_5220/338101817.py:5: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_5220/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_5220/1426773950.py:8: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_5220/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_5220/3371624330.py:7: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_5220/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.