Introduction to using OCP in molecular simulations#
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”.
Calculating adsorption energies#
Calculating adsorption energy with OCP Adsorption energies are different than you might be used to in OCP. For example, you may want the adsorption energy of oxygen, and conventionally you would compute that from this reaction:
1/2 O2 + slab -> slab-O
This is not what is done in OCP. It is referenced to a different 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
H2O + * -> O* + H2 re2 # Get from OCP as a direct calculation
Then, 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. You could also compute this using DFT.
The first step is getting a checkpoint for the model we want to use. eSCN is currently the state of the art model arXiv. This next cell will download the checkpoint if you don’t have it already.
The different models have different compute requirements. 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.
%run ocp-tutorial.ipynb
checkpoint = get_checkpoint('eSCN-L6-M3-Lay20 All+MD')
Next we load the checkpoint. The output is somewhat verbose, but it can be informative for debugging purposes.
from ocpmodels.common.relaxation.ase_utils import OCPCalculator
calc = OCPCalculator(checkpoint_path=os.path.expanduser(checkpoint), cpu=False)
# calc = OCPCalculator(checkpoint_path=os.path.expanduser(checkpoint), cpu=True)
amp: true
cmd:
checkpoint_dir: /home/jovyan/shared-scratch/jkitchin/tutorial/ocp-tutorial/checkpoints/2023-08-19-16-44-48
commit: 3973c79
identifier: ''
logs_dir: /home/jovyan/shared-scratch/jkitchin/tutorial/ocp-tutorial/logs/tensorboard/2023-08-19-16-44-48
print_every: 100
results_dir: /home/jovyan/shared-scratch/jkitchin/tutorial/ocp-tutorial/results/2023-08-19-16-44-48
seed: null
timestamp_id: 2023-08-19-16-44-48
dataset: null
gpus: 1
logger: tensorboard
model: escn
model_attributes:
basis_width_scalar: 2.0
cutoff: 12.0
distance_function: gaussian
hidden_channels: 384
lmax_list:
- 6
max_neighbors: 20
mmax_list:
- 3
num_layers: 20
num_sphere_samples: 128
otf_graph: true
regress_forces: true
sphere_channels: 160
use_pbc: true
noddp: false
optim:
batch_size: 3
clip_grad_norm: 20
ema_decay: 0.999
energy_coefficient: 4
eval_batch_size: 3
eval_every: 5000
force_coefficient: 100
loss_energy: mae
loss_force: l2mae
lr_gamma: 0.3
lr_initial: 0.0008
lr_milestones:
- 433166
- 541460
- 649750
max_epochs: 24
num_workers: 8
optimizer: AdamW
optimizer_params:
amsgrad: true
warmup_factor: 0.2
warmup_steps: 100
slurm:
additional_parameters:
constraint: volta32gb
constraint: volta32gb
cpus_per_task: 9
folder: /checkpoint/zitnick/ocp_logs/4486283
gpus_per_node: 8
job_id: '4486283'
job_name: eSCN-L6-M3-Lay20-All-MD
mem: 480GB
nodes: 4
ntasks_per_node: 8
partition: ocp
time: 4320
task:
dataset: trajectory_lmdb
description: Regressing to energies and forces for DFT trajectories from OCP
eval_on_free_atoms: true
grad_input: atomic forces
labels:
- potential energy
metric: mae
primary_metric: forces_mae
relax_opt:
alpha: 70.0
damping: 1.0
maxstep: 0.04
memory: 50
name: lbfgs
traj_dir: traj_id
relaxation_steps: 200
train_on_free_atoms: true
type: regression
write_pos: true
trainer: forces
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 fcc111, add_adsorbate
from ase.optimize import BFGS
re1 = -3.03
slab = fcc111('Pt', size=(2, 2, 5), vacuum=10.0)
add_adsorbate(slab, '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()
slab_e + re1
Step Time Energy fmax
BFGS: 0 16:45:46 1.734571 1.6863
BFGS: 1 16:45:46 1.526617 0.9688
BFGS: 2 16:45:46 1.375589 0.7231
BFGS: 3 16:45:46 1.282416 0.7910
BFGS: 4 16:45:46 0.859390 0.6855
BFGS: 5 16:45:47 0.804168 0.4963
BFGS: 6 16:45:47 0.778587 0.6211
BFGS: 7 16:45:47 0.765790 0.7337
BFGS: 8 16:45:47 0.749063 0.4882
BFGS: 9 16:45:47 0.730627 0.1268
BFGS: 10 16:45:47 0.729829 0.0355
-2.300171225070953
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()
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.6 eV difference (we predicted -2.3 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 OCP
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
atoms = Atoms(sfcc['symbols'],
positions=sfcc['pos'],
cell=sfcc['cell'],
pbc=True)
atoms.set_calculator(calc)
opt = BFGS(atoms)
opt.run(fmax=0.05, steps=100)
re2 = atoms.get_potential_energy()
nO = 0
for atom in atoms:
if atom.symbol == 'O':
nO += 1
re2 += re1 + re3
print(re2 / nO)
Step Time Energy fmax
BFGS: 0 20:12:11 0.604887 0.3249
BFGS: 1 20:12:11 0.609885 0.2843
BFGS: 2 20:12:12 0.596275 0.0778
BFGS: 3 20:12:12 0.600325 0.0716
BFGS: 4 20:12:12 0.589127 0.0940
BFGS: 5 20:12:12 0.592868 0.0914
BFGS: 6 20:12:12 0.592103 0.0874
BFGS: 7 20:12:12 0.587785 0.0681
BFGS: 8 20:12:12 0.590783 0.0244
-5.019216940402984
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.
from tqdm import tqdm
import time
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]
atoms = Atoms(entry['symbols'],
positions=entry['pos'],
cell=entry['cell'],
pbc=True)
atoms.set_calculator(calc)
opt = BFGS(atoms, logfile=None) # no logfile to suppress output
opt.run(fmax=0.05, steps=100)
re2 = atoms.get_potential_energy()
nO = 0
for atom in atoms:
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
100%|██████████| 1/1 [00:01<00:00, 1.42s/it]
100%|██████████| 1/1 [00:00<00:00, 1.00it/s]
Ag
100%|██████████| 1/1 [00:00<00:00, 1.01it/s]
100%|██████████| 1/1 [00:00<00:00, 2.34it/s]
Pd
100%|██████████| 1/1 [00:01<00:00, 1.14s/it]
100%|██████████| 1/1 [00:00<00:00, 1.17it/s]
Pt
100%|██████████| 1/1 [00:01<00:00, 1.28s/it]
100%|██████████| 1/1 [00:00<00:00, 1.41it/s]
Rh
100%|██████████| 1/1 [00:00<00:00, 2.34it/s]
100%|██████████| 1/1 [00:00<00:00, 2.32it/s]
Ir
100%|██████████| 1/1 [00:00<00:00, 2.31it/s]
100%|██████████| 1/1 [00:00<00:00, 1.75it/s]
'Elapsed time = 9.732725858688354 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('OCP 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)', 'OCP']);
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)
add_adsorbate(slab, 'O', height=1.2, position='fcc')
slab.set_calculator(calc)
opt = BFGS(slab, logfile=None)
opt.run(fmax=0.05, steps=100)
slab_e = slab.get_potential_energy()
print(f'nlayers = {nlayers}: {slab_e + re1:1.2f} eV')
nlayers = 3: -2.40 eV
nlayers = 4: -2.29 eV
nlayers = 5: -2.30 eV
nlayers = 6: -2.28 eV
nlayers = 7: -2.28 eV
nlayers = 8: -2.30 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)
add_adsorbate(slab, 'O', height=1.2, position='fcc')
slab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in slab]))
slab.set_calculator(calc)
opt = BFGS(slab, logfile=None)
opt.run(fmax=0.05, steps=100)
slab_e = slab.get_potential_energy()
print(f'nlayers = {nlayers}: {slab_e + re1:1.2f} eV')
nlayers = 3: -2.25 eV
nlayers = 4: -2.17 eV
nlayers = 5: -2.16 eV
nlayers = 6: -2.18 eV
nlayers = 7: -2.18 eV
nlayers = 8: -2.18 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)
add_adsorbate(slab, 'O', height=1.2, position='fcc')
slab.set_constraint(FixAtoms(mask=[atom.tag > 1 for atom in slab]))
slab.set_calculator(calc)
opt = BFGS(slab, logfile=None)
opt.run(fmax=0.05, steps=100)
slab_e = slab.get_potential_energy()
print(f'({size}x{size}): {slab_e + re1:1.2f} eV')
(1x1): -0.74 eV
(2x2): -2.16 eV
(3x3): -1.38 eV
(4x4): -1.30 eV
(5x5): -1.30 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.