| Property | Value |
|---|---|
| Difficulty | Beginner |
| Time | 15-30 minutes |
| Prerequisites | Basic Python, familiarity with ASE |
| Goal | Calculate adsorption energies using UMA models |
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.
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-ODFT 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* + H2It 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 UMAThen, the adsorption energy for
1/2O2 + * -> O*is just re1 + re2.
Based on https://atct.anl.gov/Thermochemical Data/version 1.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
! pip install fairchem-core fairchem-data-oc fairchem-applications-cattsunamiGet 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 loginor by setting the HF_TOKEN environment variable.
# Login using the huggingface-cli utility
! huggingface-cli login
# alternatively,
import os
os.environ['HF_TOKEN'] = 'MY_TOKEN'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-1p2")
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_9829/3752951811.py:17: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
WARNING:root:Model is being compiled this might take a while for the first time
W0807 05:54:57.583000 9829 site-packages/torch/_logging/_internal.py:1345] [0/0] Profiler record function <class 'torch.autograd.profiler.record_function'> will be ignored
Step Time Energy fmax
BFGS: 0 05:55:49 -104.694017 0.695051
BFGS: 1 05:55:50 -104.750470 0.597035
BFGS: 2 05:55:50 -104.896307 0.382718
BFGS: 3 05:55:50 -104.926053 0.441386
BFGS: 4 05:55:50 -105.022118 0.447653
BFGS: 5 05:55:50 -105.082313 0.322309
BFGS: 6 05:55:50 -105.111678 0.162520
/tmp/ipykernel_9829/3752951811.py:22: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
WARNING:root:The UMA fast path (merge_mole + compile) is only available for fixed composition, task, charge, and spin. This is optimized for MD applications. Falling back to a less optimized version for subsequent evaluations. Reason: 'Compositions differ from merged model'.
Use inference_settings='batch' for heterogeneous batched evaluations.
BFGS: 7 05:55:50 -105.122290 0.038242
Step Time Energy fmax
BFGS: 0 05:55:58 -110.077201 1.746972
BFGS: 1 05:55:58 -110.258222 0.993461
BFGS: 2 05:55:58 -110.405519 0.740253
BFGS: 3 05:55:59 -110.453435 0.792029
BFGS: 4 05:55:59 -110.570003 0.602221
BFGS: 5 05:55:59 -110.638271 0.491858
BFGS: 6 05:55:59 -110.695131 0.598342
BFGS: 7 05:56:00 -110.741332 0.612749
BFGS: 8 05:56:00 -110.772910 0.428314
BFGS: 9 05:56:00 -110.788240 0.192009
BFGS: 10 05:56:01 -110.791604 0.095433
BFGS: 11 05:56:01 -110.792269 0.094849
BFGS: 12 05:56:01 -110.793021 0.085641
BFGS: 13 05:56:01 -110.793671 0.071266
BFGS: 14 05:56:02 -110.794200 0.053573
BFGS: 15 05:56:02 -110.794439 0.042699
-1.4981483382228116It 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. Xu & Kitchin (2014)
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 eVThen, 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.263842000000002Next, 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)/tmp/ipykernel_9829/647904475.py:10: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
Step Time Energy fmax
BFGS: 0 05:56:05 -82.881492 1.012517
BFGS: 1 05:56:06 -82.940117 0.758967
BFGS: 2 05:56:06 -83.035745 0.334363
BFGS: 3 05:56:06 -83.039932 0.304527
BFGS: 4 05:56:06 -83.049373 0.206965
BFGS: 5 05:56:07 -83.054436 0.140383
BFGS: 6 05:56:07 -83.057107 0.076567
BFGS: 7 05:56:07 -83.057955 0.064835
BFGS: 8 05:56:07 -83.058538 0.066526
BFGS: 9 05:56:08 -83.058830 0.045645
/tmp/ipykernel_9829/647904475.py:14: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
Step Time Energy fmax
BFGS: 0 05:56:08 -88.773356 0.334878
BFGS: 1 05:56:08 -88.777414 0.290913
BFGS: 2 05:56:08 -88.789784 0.119410
BFGS: 3 05:56:09 -88.791837 0.124366
BFGS: 4 05:56:09 -88.795399 0.130636
BFGS: 5 05:56:09 -88.797986 0.118737
BFGS: 6 05:56:10 -88.800127 0.085202
BFGS: 7 05:56:10 -88.801291 0.091586
BFGS: 8 05:56:10 -88.802145 0.065007
BFGS: 9 05:56:10 -88.802709 0.042341
-4.149879082903078
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]/tmp/ipykernel_9829/1356342052.py:33: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
Step Time Energy fmax
BFGS: 0 05:56:11 -48.890191 0.646801
BFGS: 1 05:56:11 -48.913215 0.542334
BFGS: 2 05:56:11 -48.978855 0.272943
BFGS: 3 05:56:12 -48.980963 0.248402
BFGS: 4 05:56:12 -48.989843 0.142195
BFGS: 5 05:56:12 -48.994133 0.109508
BFGS: 6 05:56:12 -48.995943 0.057373
BFGS: 7 05:56:13 -48.996412 0.052332
BFGS: 8 05:56:13 -48.996839 0.050481
BFGS: 9 05:56:13 -48.997190 0.037852
Step Time Energy fmax
BFGS: 0 05:56:13 -55.183791 0.317002
/tmp/ipykernel_9829/1356342052.py:37: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
BFGS: 1 05:56:13 -55.186079 0.260545
BFGS: 2 05:56:13 -55.194759 0.163272
BFGS: 3 05:56:13 -55.196862 0.156293
BFGS: 4 05:56:14 -55.200264 0.089518
BFGS: 5 05:56:14 -55.202051 0.085341
BFGS: 6 05:56:14 -55.203480 0.085305
BFGS: 7 05:56:14 -55.204715 0.106556
BFGS: 8 05:56:14 -55.206073 0.098622
BFGS: 9 05:56:14 -55.206943 0.055586
100%|██████████| 1/1 [00:04<00:00, 4.11s/it]100%|██████████| 1/1 [00:04<00:00, 4.11s/it]
BFGS: 10 05:56:15 -55.207300 0.041434
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:15 -48.915497 0.555938
BFGS: 1 05:56:15 -48.933106 0.473539
BFGS: 2 05:56:15 -48.987176 0.208265
BFGS: 3 05:56:15 -48.988350 0.196053
BFGS: 4 05:56:15 -48.996556 0.040859
Step Time Energy fmax
BFGS: 0 05:56:15 -55.087818 0.314616
BFGS: 1 05:56:16 -55.089884 0.253627
BFGS: 2 05:56:16 -55.096718 0.155815
BFGS: 3 05:56:16 -55.098611 0.158030
BFGS: 4 05:56:16 -55.101906 0.102147
BFGS: 5 05:56:16 -55.103301 0.061643
BFGS: 6 05:56:16 -55.104082 0.058925
BFGS: 7 05:56:17 -55.104754 0.080596
BFGS: 8 05:56:17 -55.105736 0.089870
BFGS: 9 05:56:17 -55.106587 0.063183
100%|██████████| 1/1 [00:02<00:00, 2.90s/it]100%|██████████| 1/1 [00:02<00:00, 2.91s/it]
BFGS: 10 05:56:18 -55.106974 0.033877
Ag
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:18 -33.015774 0.626056
BFGS: 1 05:56:18 -33.034755 0.545928
BFGS: 2 05:56:18 -33.103217 0.188722
BFGS: 3 05:56:19 -33.104919 0.179778
BFGS: 4 05:56:19 -33.106624 0.166504
BFGS: 5 05:56:19 -33.109909 0.127307
BFGS: 6 05:56:19 -33.113531 0.109334
BFGS: 7 05:56:20 -33.115653 0.053176
BFGS: 8 05:56:20 -33.116105 0.039454
Step Time Energy fmax
BFGS: 0 05:56:20 -38.158732 0.127418
BFGS: 1 05:56:20 -38.159781 0.119851
BFGS: 2 05:56:20 -38.170047 0.074136
BFGS: 3 05:56:21 -38.171022 0.083949
BFGS: 4 05:56:21 -38.174079 0.098329
BFGS: 5 05:56:21 -38.176468 0.091288
BFGS: 6 05:56:21 -38.178595 0.073994
BFGS: 7 05:56:22 -38.179925 0.083156
BFGS: 8 05:56:22 -38.180628 0.065178
100%|██████████| 1/1 [00:04<00:00, 4.66s/it]100%|██████████| 1/1 [00:04<00:00, 4.66s/it]
BFGS: 9 05:56:22 -38.180974 0.037563
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:23 -33.037809 0.552148
BFGS: 1 05:56:23 -33.052341 0.486542
BFGS: 2 05:56:23 -33.108823 0.155455
BFGS: 3 05:56:23 -33.109719 0.145457
BFGS: 4 05:56:23 -33.111363 0.118068
BFGS: 5 05:56:24 -33.113395 0.079860
BFGS: 6 05:56:24 -33.115417 0.053419
BFGS: 7 05:56:24 -33.116023 0.030041
Step Time Energy fmax
BFGS: 0 05:56:24 -38.073489 0.119078
BFGS: 1 05:56:24 -38.074444 0.115486
BFGS: 2 05:56:25 -38.085485 0.074974
BFGS: 3 05:56:25 -38.086458 0.080957
BFGS: 4 05:56:25 -38.088451 0.081874
BFGS: 5 05:56:26 -38.089967 0.070834
BFGS: 6 05:56:26 -38.091909 0.042921
100%|██████████| 1/1 [00:03<00:00, 3.63s/it]100%|██████████| 1/1 [00:03<00:00, 3.63s/it]
Pd
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:26 -70.174812 0.646926
BFGS: 1 05:56:27 -70.200771 0.520374
BFGS: 2 05:56:27 -70.253580 0.195411
BFGS: 3 05:56:27 -70.255196 0.188046
BFGS: 4 05:56:27 -70.262548 0.137571
BFGS: 5 05:56:28 -70.265013 0.106889
BFGS: 6 05:56:28 -70.266655 0.074815
BFGS: 7 05:56:28 -70.267604 0.061562
BFGS: 8 05:56:28 -70.268390 0.035482
Step Time Energy fmax
BFGS: 0 05:56:29 -76.139648 0.221609
BFGS: 1 05:56:29 -76.142815 0.197612
BFGS: 2 05:56:29 -76.157140 0.181487
BFGS: 3 05:56:29 -76.159369 0.159807
BFGS: 4 05:56:30 -76.163719 0.132238
BFGS: 5 05:56:30 -76.166338 0.105136
BFGS: 6 05:56:30 -76.168581 0.099640
BFGS: 7 05:56:30 -76.169880 0.098468
BFGS: 8 05:56:30 -76.170694 0.077171
BFGS: 9 05:56:31 -76.171128 0.049108
100%|██████████| 1/1 [00:04<00:00, 4.83s/it]100%|██████████| 1/1 [00:04<00:00, 4.84s/it]
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:31 -70.208058 0.465459
BFGS: 1 05:56:31 -70.222887 0.381518
BFGS: 2 05:56:32 -70.257808 0.181345
BFGS: 3 05:56:32 -70.259016 0.170180
BFGS: 4 05:56:32 -70.266016 0.073034
BFGS: 5 05:56:32 -70.266661 0.070113
BFGS: 6 05:56:32 -70.267869 0.048762
Step Time Energy fmax
BFGS: 0 05:56:33 -75.957684 0.183791
BFGS: 1 05:56:33 -75.960869 0.164287
BFGS: 2 05:56:33 -75.970372 0.169078
BFGS: 3 05:56:33 -75.972327 0.164641
BFGS: 4 05:56:33 -75.977419 0.120429
BFGS: 5 05:56:34 -75.979833 0.110607
BFGS: 6 05:56:34 -75.981659 0.082725
BFGS: 7 05:56:34 -75.982830 0.078026
BFGS: 8 05:56:34 -75.983552 0.046149
100%|██████████| 1/1 [00:03<00:00, 3.59s/it]100%|██████████| 1/1 [00:03<00:00, 3.60s/it]
Pt
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:35 -82.881491 1.012517
BFGS: 1 05:56:35 -82.940117 0.758968
BFGS: 2 05:56:35 -83.035745 0.334363
BFGS: 3 05:56:35 -83.039932 0.304531
BFGS: 4 05:56:36 -83.049371 0.206955
BFGS: 5 05:56:36 -83.054436 0.140386
BFGS: 6 05:56:36 -83.057105 0.076576
BFGS: 7 05:56:37 -83.057954 0.064860
BFGS: 8 05:56:37 -83.058538 0.066522
BFGS: 9 05:56:37 -83.058830 0.045636
Step Time Energy fmax
BFGS: 0 05:56:37 -88.773355 0.334878
BFGS: 1 05:56:38 -88.777414 0.290913
BFGS: 2 05:56:38 -88.789784 0.119410
BFGS: 3 05:56:38 -88.791837 0.124366
BFGS: 4 05:56:38 -88.795399 0.130638
BFGS: 5 05:56:39 -88.797988 0.118733
BFGS: 6 05:56:39 -88.800128 0.085211
BFGS: 7 05:56:39 -88.801291 0.091583
BFGS: 8 05:56:39 -88.802145 0.065008
BFGS: 9 05:56:39 -88.802709 0.042340
100%|██████████| 1/1 [00:04<00:00, 4.92s/it]100%|██████████| 1/1 [00:04<00:00, 4.92s/it]
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:40 -82.968454 0.688065
BFGS: 1 05:56:40 -82.995521 0.558978
BFGS: 2 05:56:40 -83.049826 0.200180
BFGS: 3 05:56:40 -83.051294 0.185780
BFGS: 4 05:56:40 -83.057463 0.066862
BFGS: 5 05:56:41 -83.057908 0.055435
BFGS: 6 05:56:41 -83.058686 0.031712
Step Time Energy fmax
BFGS: 0 05:56:41 -88.396682 0.203983
BFGS: 1 05:56:41 -88.400062 0.174390
BFGS: 2 05:56:41 -88.408718 0.136409
BFGS: 3 05:56:41 -88.410332 0.134069
BFGS: 4 05:56:42 -88.414336 0.095108
BFGS: 5 05:56:42 -88.415784 0.090170
BFGS: 6 05:56:42 -88.417002 0.104754
BFGS: 7 05:56:42 -88.417817 0.094188
BFGS: 8 05:56:42 -88.418373 0.051483
BFGS: 9 05:56:43 -88.418641 0.031991
100%|██████████| 1/1 [00:03<00:00, 3.61s/it]100%|██████████| 1/1 [00:03<00:00, 3.62s/it]
Rh
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:43 -100.191090 0.703650
BFGS: 1 05:56:44 -100.219546 0.608665
BFGS: 2 05:56:44 -100.286016 0.178447
BFGS: 3 05:56:44 -100.289172 0.138491
BFGS: 4 05:56:44 -100.299063 0.074319
BFGS: 5 05:56:45 -100.300012 0.067005
BFGS: 6 05:56:45 -100.301061 0.062656
BFGS: 7 05:56:45 -100.301787 0.049559
Step Time Energy fmax
BFGS: 0 05:56:45 -106.949006 0.238786
BFGS: 1 05:56:46 -106.954055 0.199552
BFGS: 2 05:56:46 -106.963647 0.066008
BFGS: 3 05:56:46 -106.963942 0.058097
BFGS: 4 05:56:46 -106.964805 0.029320
100%|██████████| 1/1 [00:03<00:00, 3.49s/it]100%|██████████| 1/1 [00:03<00:00, 3.49s/it]
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:47 -100.169507 0.774354
BFGS: 1 05:56:47 -100.204661 0.634293
BFGS: 2 05:56:47 -100.287440 0.228397
BFGS: 3 05:56:48 -100.290399 0.177250
BFGS: 4 05:56:48 -100.298467 0.080649
BFGS: 5 05:56:48 -100.299792 0.067136
BFGS: 6 05:56:48 -100.300911 0.052046
BFGS: 7 05:56:49 -100.301597 0.051987
BFGS: 8 05:56:49 -100.302239 0.043631
Step Time Energy fmax
BFGS: 0 05:56:49 -106.904528 0.271607
BFGS: 1 05:56:50 -106.909937 0.214886
BFGS: 2 05:56:50 -106.920133 0.083912
BFGS: 3 05:56:50 -106.920388 0.076302
BFGS: 4 05:56:51 -106.921195 0.030130
100%|██████████| 1/1 [00:04<00:00, 4.37s/it]100%|██████████| 1/1 [00:04<00:00, 4.38s/it]
Ir
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:51 -124.226220 1.208047
BFGS: 1 05:56:51 -124.303068 0.944243
BFGS: 2 05:56:52 -124.413829 0.177421
BFGS: 3 05:56:52 -124.417069 0.150391
BFGS: 4 05:56:52 -124.423699 0.053157
BFGS: 5 05:56:53 -124.424242 0.050764
BFGS: 6 05:56:53 -124.424757 0.044990
Step Time Energy fmax
BFGS: 0 05:56:53 -130.642842 0.410127
BFGS: 1 05:56:53 -130.656118 0.294537
BFGS: 2 05:56:54 -130.672540 0.084505
BFGS: 3 05:56:54 -130.673412 0.070182
BFGS: 4 05:56:54 -130.674099 0.062829
BFGS: 5 05:56:54 -130.675015 0.055398
BFGS: 6 05:56:55 -130.675473 0.050562
BFGS: 7 05:56:55 -130.675726 0.043154
100%|██████████| 1/1 [00:04<00:00, 4.29s/it]100%|██████████| 1/1 [00:04<00:00, 4.29s/it]
0%| | 0/1 [00:00<?, ?it/s] Step Time Energy fmax
BFGS: 0 05:56:55 -124.219231 1.178119
BFGS: 1 05:56:56 -124.302855 0.920683
BFGS: 2 05:56:56 -124.415676 0.214575
BFGS: 3 05:56:56 -124.418174 0.200550
BFGS: 4 05:56:56 -124.424014 0.075790
BFGS: 5 05:56:56 -124.424558 0.045441
Step Time Energy fmax
BFGS: 0 05:56:56 -130.530228 0.471262
BFGS: 1 05:56:57 -130.546421 0.333597
BFGS: 2 05:56:57 -130.565046 0.073569
BFGS: 3 05:56:57 -130.566054 0.078632
BFGS: 4 05:56:57 -130.566520 0.069134
BFGS: 5 05:56:58 -130.567331 0.073884
BFGS: 6 05:56:58 -130.567673 0.065453
100%|██████████| 1/1 [00:03<00:00, 3.29s/it]BFGS: 7 05:56:58 -130.567876 0.043100
100%|██████████| 1/1 [00:03<00:00, 3.29s/it]
'Elapsed time = 47.78350043296814 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 the adsorption energies section 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_9829/338101817.py:5: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9829/338101817.py:14: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
nlayers = 3: -1.64 eV
nlayers = 4: -1.47 eV
nlayers = 5: -1.50 eV
nlayers = 6: -1.48 eV
nlayers = 7: -1.49 eV
nlayers = 8: -1.49 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_9829/1426773950.py:8: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9829/1426773950.py:18: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
nlayers = 3: -1.54 eV
nlayers = 4: -1.35 eV
nlayers = 5: -1.38 eV
nlayers = 6: -1.37 eV
nlayers = 7: -1.38 eV
nlayers = 8: -1.38 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_9829/3371624330.py:7: FutureWarning: Please use atoms.calc = calc
slab.set_calculator(calc)
/tmp/ipykernel_9829/3371624330.py:17: FutureWarning: Please use atoms.calc = calc
adslab.set_calculator(calc)
(1x1): -0.22 eV
(2x2): -1.38 eV
(3x3): -1.43 eV
(4x4): -1.45 eV
(5x5): -1.46 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.
- Xu, Z., & Kitchin, J. R. (2014). Probing the Coverage Dependence of Site and Adsorbate Configurational Correlations on (111) Surfaces of Late Transition Metals. The Journal of Physical Chemistry C, 118(44), 25597–25602. 10.1021/jp508805h