UMA Intro Tutorial#
This tutorial will walk you through a few examples of how you can use UMA. Each step is covered in more detail elsewhere in the documentation, but this is well suited to a ~1-2 hour tutorial session for researchers new to UMA but with some background in ASE and molecular simulations.
Before you start / installation#
You need to get a HuggingFace account and request access to the UMA models.
You need a Huggingface account, request access to https://huggingface.co/facebook/UMA, and to create a Huggingface token at https://huggingface.co/settings/tokens/ with these permission:
Permissions: Read access to contents of all public gated repos you can access
Then, add the token as an environment variable (using huggingface-cli login:
# Enter token via huggingface-cli
! huggingface-cli login
or you can set the token via HF_TOKEN variable:
# Set token via env variable
import os
os.environ['HF_TOKEN'] = 'MYTOKEN'
Installation process#
It may be enough to use pip install fairchem-core. This gets you the latest version on PyPi (https://pypi.org/project/fairchem-core/)
Here we install some sub-packages. This can take 2-5 minutes to run.
! pip install fairchem-core fairchem-data-oc fairchem-applications-cattsunami x3dase
# Check that packages are installed
!pip list | grep fairchem
fairchem-applications-cattsunami 1.1.2.dev66+g85bd8353 /home/runner/work/fairchem/fairchem/packages/fairchem-applications-cattsunami
fairchem-core 2.9.1.dev4+g85bd8353 /home/runner/work/fairchem/fairchem/packages/fairchem-core
fairchem-data-oc 1.0.3.dev66+g85bd8353 /home/runner/work/fairchem/fairchem/packages/fairchem-data-oc
import fairchem.core
fairchem.core.__version__
'2.9.1.dev4+g85bd8353'
Illustrative examples#
These should just run, and are here to show some basic uses.
Critical points:
Create a calculator
Specify the task_name
Use calculator like other ASE calculators
Spin gap energy - OMOL#
This is the difference in energy between a triplet and single ground state for a CH2 radical. This downloads a ~1GB checkpoint the first time you run it.
We don’t set a device here, so we get a warning about using a CPU device. You can ignore that. If a CUDA environment is available, a GPU may be used to speed up the calculations.
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
WARNING:root:device was not explicitly set, using device='cuda'.
/home/runner/work/_tool/Python/3.12.12/x64/lib/python3.12/site-packages/pydantic/_internal/_generate_schema.py:2249: UnsupportedFieldAttributeWarning: The 'repr' attribute with value False was provided to the `Field()` function, which has no effect in the context it was used. 'repr' is field-specific metadata, and can only be attached to a model field using `Annotated` metadata or by assignment. This may have happened because an `Annotated` type alias using the `type` statement was used, or if the `Field()` function was attached to a single member of a union type.
warnings.warn(
/home/runner/work/_tool/Python/3.12.12/x64/lib/python3.12/site-packages/pydantic/_internal/_generate_schema.py:2249: UnsupportedFieldAttributeWarning: The 'frozen' attribute with value True was provided to the `Field()` function, which has no effect in the context it was used. 'frozen' is field-specific metadata, and can only be attached to a model field using `Annotated` metadata or by assignment. This may have happened because an `Annotated` type alias using the `type` statement was used, or if the `Field()` function was attached to a single member of a union type.
warnings.warn(
from ase.build import molecule
# singlet CH2
singlet = molecule("CH2_s1A1d")
singlet.info.update({"spin": 1, "charge": 0})
singlet.calc = FAIRChemCalculator(predictor, task_name="omol")
# triplet CH2
triplet = molecule("CH2_s3B1d")
triplet.info.update({"spin": 3, "charge": 0})
triplet.calc = FAIRChemCalculator(predictor, task_name="omol")
print(triplet.get_potential_energy() - singlet.get_potential_energy())
-0.5508372783660889
Example of adsorbate relaxation - OC20#
Here we just setup a Cu(100) slab with a CO on it and relax it.
This is an OC20 task because it is a slab with an adsorbate.
We specify an explicit device in the predictor here, and avoid the warning.
from ase.build import add_adsorbate, fcc100, molecule
from ase.optimize import LBFGS
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="oc20")
# Set up your system as an ASE atoms object
slab = fcc100("Cu", (3, 3, 3), vacuum=8, periodic=True)
adsorbate = molecule("CO")
add_adsorbate(slab, adsorbate, 2.0, "bridge")
slab.calc = calc
# Set up LBFGS dynamics object
opt = LBFGS(slab)
opt.run(0.05, 100)
print(slab.get_potential_energy())
WARNING:root:device was not explicitly set, using device='cuda'.
Step Time Energy fmax
LBFGS: 0 17:49:28 -89.596203 11.451673
LBFGS: 1 17:49:28 -92.497567 6.543859
LBFGS: 2 17:49:29 -92.624418 7.536267
LBFGS: 3 17:49:29 -93.000882 3.716082
LBFGS: 4 17:49:29 -93.158665 3.480264
LBFGS: 5 17:49:29 -93.264137 2.256967
LBFGS: 6 17:49:29 -93.505253 1.133176
LBFGS: 7 17:49:29 -93.595944 0.991884
LBFGS: 8 17:49:29 -93.705359 0.683701
LBFGS: 9 17:49:29 -93.791540 0.506523
LBFGS: 10 17:49:29 -93.837935 0.364215
LBFGS: 11 17:49:29 -93.856811 0.349524
LBFGS: 12 17:49:29 -93.881778 0.498544
LBFGS: 13 17:49:29 -93.900200 0.432719
LBFGS: 14 17:49:30 -93.910017 0.156140
LBFGS: 15 17:49:30 -93.915885 0.170016
LBFGS: 16 17:49:30 -93.922143 0.211887
LBFGS: 17 17:49:30 -93.929011 0.260653
LBFGS: 18 17:49:30 -93.935159 0.183819
LBFGS: 19 17:49:30 -93.938071 0.057453
LBFGS: 20 17:49:30 -93.938499 0.039144
-93.93849899606094
Example bulk relaxation - OMAT#
from ase.build import bulk
from ase.filters import FrechetCellFilter
from ase.optimize import FIRE
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="omat")
atoms = bulk("Fe")
atoms.calc = calc
opt = FIRE(FrechetCellFilter(atoms))
opt.run(0.05, 100)
print(atoms.get_stress()) # !!!! We get stress now!
WARNING:root:device was not explicitly set, using device='cuda'.
Step Time Energy fmax
FIRE: 0 17:49:32 -8.261158 0.651784
FIRE: 1 17:49:32 -8.271310 0.358118
FIRE: 2 17:49:32 -8.264588 1.650199
FIRE: 3 17:49:32 -8.273672 0.177967
FIRE: 4 17:49:32 -8.272633 0.269093
FIRE: 5 17:49:32 -8.272766 0.257552
FIRE: 6 17:49:32 -8.273009 0.234348
FIRE: 7 17:49:33 -8.273319 0.199201
FIRE: 8 17:49:33 -8.273635 0.151747
FIRE: 9 17:49:33 -8.273890 0.091454
FIRE: 10 17:49:33 -8.274015 0.017801
[ 1.5562733e-03 1.5562323e-03 1.5562797e-03 3.2104857e-08
-7.8731563e-09 4.2523460e-08]
Molecular dynamics - OMOL#
import matplotlib.pyplot as plt
from ase import units
from ase.build import molecule
from ase.io import Trajectory
from ase.md.langevin import Langevin
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="omol")
atoms = molecule("H2O")
atoms.info.update(charge=0, spin=1) # For omol
atoms.calc = calc
dyn = Langevin(
atoms,
timestep=0.1 * units.fs,
temperature_K=400,
friction=0.001 / units.fs,
)
trajectory = Trajectory("my_md.traj", "w", atoms)
dyn.attach(trajectory.write, interval=1)
dyn.run(steps=50)
# See some results - not paper ready!
traj = Trajectory("my_md.traj")
plt.plot(
[i * 0.1 * units.fs for i in range(len(traj))],
[a.get_potential_energy() for a in traj],
)
plt.xlabel("Time (fs)")
plt.ylabel("Energy (eV)");
WARNING:root:device was not explicitly set, using device='cuda'.
Catalyst Adsorption energies#
The basic approach in computing an adsorption energy is to compute this energy difference:
dH = E_adslab - E_slab - E_ads
We use UMA for two of these energies E_adslab and E_slab. For E_ads We have to do something a little different. The OC20 task is not trained for molecules or molecular fragments. We use atomic energy reference energies instead. These are tabulated below.
The OC20 reference scheme is this reaction:
x CO + (x + y/2 - z) H2 + (z-x) H2O + w/2 N2 + * -> CxHyOzNw*
For this example we have
-H2 + H2O + * -> O*. "O": -7.204 eV
Where "O": -7.204 is a constant.
To get the desired reaction energy we want we add the formation energy of water. We use either DFT or experimental values for this reaction energy.
1/2O2 + H2 -> H2O
Alternatives to this approach are using DFT to estimate the energy of 1/2 O2, just make sure to use consistent settings with your task. You should not use OMOL for this.
from ase.build import add_adsorbate, fcc111
from ase.optimize import BFGS
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="oc20")
WARNING:root:device was not explicitly set, using device='cuda'.
# 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 # Water formation energy from experiment
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.calc = calc
opt = BFGS(slab)
print("Relaxing slab")
opt.run(fmax=0.05, steps=100)
slab_e = slab.get_potential_energy()
adslab.calc = calc
opt = BFGS(adslab)
print("\nRelaxing adslab")
opt.run(fmax=0.05, steps=100)
adslab_e = adslab.get_potential_energy()
Relaxing slab
Step Time Energy fmax
BFGS: 0 17:49:41 -104.710392 0.708597
BFGS: 1 17:49:41 -104.767890 0.606681
BFGS: 2 17:49:41 -104.919426 0.369421
BFGS: 3 17:49:41 -104.952876 0.441471
BFGS: 4 17:49:41 -105.030008 0.467628
BFGS: 5 17:49:42 -105.091456 0.365310
BFGS: 6 17:49:42 -105.128722 0.195170
BFGS: 7 17:49:42 -105.143316 0.049006
Relaxing adslab
Step Time Energy fmax
BFGS: 0 17:49:42 -110.055657 1.762239
BFGS: 1 17:49:42 -110.239040 0.996664
BFGS: 2 17:49:42 -110.389559 0.747527
BFGS: 3 17:49:42 -110.441190 0.818345
BFGS: 4 17:49:42 -110.557369 0.688399
BFGS: 5 17:49:42 -110.631216 0.497355
BFGS: 6 17:49:42 -110.687289 0.690846
BFGS: 7 17:49:42 -110.737890 0.729422
BFGS: 8 17:49:43 -110.774875 0.435693
BFGS: 9 17:49:43 -110.786663 0.199915
BFGS: 10 17:49:43 -110.789560 0.080679
BFGS: 11 17:49:43 -110.790038 0.057993
BFGS: 12 17:49:43 -110.790288 0.044003
Now we compute the adsorption energy.
# Energy for ((H2O-H2) + * -> *O) + (H2 + 1/2O2 -> H2O) leads to 1/2O2 + * -> *O!
adslab_e - slab_e - atomic_reference_energies["O"] + re1
-1.4729716975518907
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.
It is always a good idea to visualize the geometries to make sure they look reasonable.
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()
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()
Molecular vibrations#
from ase import Atoms
from ase.optimize import BFGS
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="omol")
from ase.vibrations import Vibrations
n2 = Atoms("N2", [(0, 0, 0), (0, 0, 1.1)])
n2.info.update({"spin": 1, "charge": 0})
n2.calc = calc
BFGS(n2).run(fmax=0.01)
WARNING:root:device was not explicitly set, using device='cuda'.
Step Time Energy fmax
BFGS: 0 17:49:47 -2981.068009 1.645285
BFGS: 1 17:49:47 -2980.961841 6.601501
BFGS: 2 17:49:47 -2981.076753 0.203644
BFGS: 3 17:49:47 -2981.076882 0.024169
BFGS: 4 17:49:47 -2981.076883 0.000104
np.True_
vib = Vibrations(n2)
vib.run()
vib.summary()
---------------------
# meV cm^-1
---------------------
0 0.0i 0.0i
1 0.0i 0.0i
2 0.0 0.0
3 2.0 16.1
4 2.0 16.1
5 309.2 2494.2
---------------------
Zero-point energy: 0.157 eV
Bulk alloy phase behavior#
Adapted from https://kitchingroup.cheme.cmu.edu/dft-book/dft.html#orgheadline29
We manually compute the formation energy of pure compounds and some alloy compositions to assess stability.
from ase.atoms import Atom, Atoms
from ase.filters import FrechetCellFilter
from ase.optimize import FIRE
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
cu = Atoms(
[Atom("Cu", [0.000, 0.000, 0.000])],
cell=[[1.818, 0.000, 1.818], [1.818, 1.818, 0.000], [0.000, 1.818, 1.818]],
pbc=True,
)
cu.calc = FAIRChemCalculator(predictor, task_name="omat")
opt = FIRE(FrechetCellFilter(cu))
opt.run(0.05, 100)
cu.get_potential_energy()
WARNING:root:device was not explicitly set, using device='cuda'.
Step Time Energy fmax
FIRE: 0 17:49:50 -3.756933 0.161999
FIRE: 1 17:49:50 -3.757594 0.110083
FIRE: 2 17:49:50 -3.758130 0.020766
-3.758130096599536
pd = Atoms(
[Atom("Pd", [0.000, 0.000, 0.000])],
cell=[[1.978, 0.000, 1.978], [1.978, 1.978, 0.000], [0.000, 1.978, 1.978]],
pbc=True,
)
pd.calc = FAIRChemCalculator(predictor, task_name="omat")
opt = FIRE(FrechetCellFilter(pd))
opt.run(0.05, 100)
pd.get_potential_energy()
Step Time Energy fmax
FIRE: 0 17:49:50 -5.211726 0.240057
FIRE: 1 17:49:50 -5.213070 0.131580
FIRE: 2 17:49:50 -5.213503 0.060259
FIRE: 3 17:49:51 -5.213528 0.051644
FIRE: 4 17:49:51 -5.213565 0.035872
-5.213564696643232
Alloy formation energies#
cupd1 = Atoms(
[Atom("Cu", [0.000, 0.000, 0.000]), Atom("Pd", [-1.652, 0.000, 2.039])],
cell=[[0.000, -2.039, 2.039], [0.000, 2.039, 2.039], [-3.303, 0.000, 0.000]],
pbc=True,
) # Note pbc=True is important, it is not the default and OMAT
cupd1.calc = FAIRChemCalculator(predictor, task_name="omat")
opt = FIRE(FrechetCellFilter(cupd1))
opt.run(0.05, 100)
cupd1.get_potential_energy()
Step Time Energy fmax
FIRE: 0 17:49:51 -9.202820 0.142029
FIRE: 1 17:49:51 -9.203042 2.754047
FIRE: 2 17:49:51 -9.201167 0.106312
FIRE: 3 17:49:51 -9.201227 0.104707
FIRE: 4 17:49:51 -9.201344 0.101568
FIRE: 5 17:49:51 -9.201514 0.097055
FIRE: 6 17:49:51 -9.201730 0.091412
FIRE: 7 17:49:51 -9.201985 0.084958
FIRE: 8 17:49:52 -9.202272 0.078093
FIRE: 9 17:49:52 -9.202583 0.071408
FIRE: 10 17:49:52 -9.202951 0.069736
FIRE: 11 17:49:52 -9.203377 0.067345
FIRE: 12 17:49:52 -9.203862 0.063554
FIRE: 13 17:49:52 -9.204403 0.057514
FIRE: 14 17:49:52 -9.204991 0.051090
FIRE: 15 17:49:52 -9.205607 0.054900
FIRE: 16 17:49:52 -9.206231 0.060124
FIRE: 17 17:49:52 -9.206858 0.063154
FIRE: 18 17:49:53 -9.207512 0.060215
FIRE: 19 17:49:53 -9.208219 0.050452
FIRE: 20 17:49:53 -9.209001 0.038185
-9.209001477498369
cupd2 = Atoms(
[
Atom("Cu", [-0.049, 0.049, 0.049]),
Atom("Cu", [-11.170, 11.170, 11.170]),
Atom("Pd", [-7.415, 7.415, 7.415]),
Atom("Pd", [-3.804, 3.804, 3.804]),
],
cell=[[-5.629, 3.701, 5.629], [-3.701, 5.629, 5.629], [-5.629, 5.629, 3.701]],
pbc=True,
)
cupd2.calc = FAIRChemCalculator(predictor, task_name="omat")
opt = FIRE(FrechetCellFilter(cupd2))
opt.run(0.05, 100)
cupd2.get_potential_energy()
Step Time Energy fmax
FIRE: 0 17:49:53 -18.126594 0.181633
FIRE: 1 17:49:53 -18.127546 0.162952
FIRE: 2 17:49:53 -18.129066 0.127293
FIRE: 3 17:49:53 -18.130534 0.078068
FIRE: 4 17:49:53 -18.131347 0.021486
-18.13134669450632
# Delta Hf cupd-1 = -0.11 eV/atom
hf1 = (
cupd1.get_potential_energy() - cu.get_potential_energy() - pd.get_potential_energy()
)
hf1
-0.23730668425560086
# DFT: Delta Hf cupd-2 = -0.04 eV/atom
hf2 = (
cupd2.get_potential_energy()
- 2 * cu.get_potential_energy()
- 2 * pd.get_potential_energy()
)
hf2
-0.18795710802078425
hf1 - hf2, (-0.11 - -0.04)
(-0.04934957623481662, -0.07)
These indicate that cupd-1 and cupd-2 are both more stable than phase separated Cu and Pd, and that cupd-1 is more stable than cupd-2. The absolute formation energies differ from the DFT references, but the relative differences are quite close. The absolute differences could be due to DFT parameter choices (XC, psp, etc.).
Phonon calculation#
This takes 4-10 minutes. Adapted from https://wiki.fysik.dtu.dk/ase/ase/phonons.html#example.
Phonons have applications in computing the stability and free energy of solids. See:
https://www.sciencedirect.com/science/article/pii/S1359646215003127
https://iopscience.iop.org/book/mono/978-0-7503-2572-1/chapter/bk978-0-7503-2572-1ch1
from ase.build import bulk
from ase.phonons import Phonons
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="omat")
# Setup crystal
atoms = bulk("Al", "fcc", a=4.05)
# Phonon calculator
N = 7
ph = Phonons(atoms, calc, supercell=(N, N, N), delta=0.05)
ph.run()
# Read forces and assemble the dynamical matrix
ph.read(acoustic=True)
ph.clean()
path = atoms.cell.bandpath("GXULGK", npoints=100)
bs = ph.get_band_structure(path)
dos = ph.get_dos(kpts=(20, 20, 20)).sample_grid(npts=100, width=1e-3)
WARNING:root:device was not explicitly set, using device='cuda'.
WARNING, 3 imaginary frequencies at q = ( 0.00, 0.00, 0.00) ; (omega_q = 9.460e-09*i)
WARNING, 3 imaginary frequencies at q = ( 0.00, 0.00, 0.00) ; (omega_q = 9.460e-09*i)
# Plot the band structure and DOS:
import matplotlib.pyplot as plt # noqa
fig = plt.figure(figsize=(7, 4))
ax = fig.add_axes([0.12, 0.07, 0.67, 0.85])
emax = 0.04
bs.plot(ax=ax, emin=0.0, emax=emax)
dosax = fig.add_axes([0.8, 0.07, 0.17, 0.85])
dosax.fill_between(
dos.get_weights(),
dos.get_energies(),
y2=0,
color="grey",
edgecolor="k",
lw=1,
)
dosax.set_ylim(0, emax)
dosax.set_yticks([])
dosax.set_xticks([])
dosax.set_xlabel("DOS", fontsize=18);
Transition States (NEBs)#
Nudged elastic band calculations are among the most costly calculations we do. UMA makes these quicker!
Get initial state
Get final state
Construct band and interpolate the images
Relax the band
Analyze and plot the band.
We explore diffusion of an O adatom from an hcp to an fcc site on Pt(111).
Initial state#
from ase.build import add_adsorbate, fcc111, molecule
from ase.optimize import LBFGS
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1")
calc = FAIRChemCalculator(predictor, task_name="oc20")
# Set up your system as an ASE atoms object
initial = fcc111("Pt", (3, 3, 3), vacuum=8, periodic=True)
adsorbate = molecule("O")
add_adsorbate(initial, adsorbate, 2.0, "fcc")
initial.calc = calc
# Set up LBFGS dynamics object
opt = LBFGS(initial)
opt.run(0.05, 100)
print(initial.get_potential_energy())
WARNING:root:device was not explicitly set, using device='cuda'.
Step Time Energy fmax
LBFGS: 0 17:50:01 -141.329801 3.509945
LBFGS: 1 17:50:01 -141.719920 3.515035
LBFGS: 2 17:50:01 -142.980940 2.978375
LBFGS: 3 17:50:01 -143.684060 0.968228
LBFGS: 4 17:50:02 -143.787057 1.271701
LBFGS: 5 17:50:02 -143.858783 0.874640
LBFGS: 6 17:50:02 -143.933930 0.170905
LBFGS: 7 17:50:02 -143.937122 0.152452
LBFGS: 8 17:50:02 -143.944600 0.122233
LBFGS: 9 17:50:02 -143.948829 0.109291
LBFGS: 10 17:50:02 -143.952238 0.069969
LBFGS: 11 17:50:02 -143.953716 0.080101
LBFGS: 12 17:50:02 -143.955175 0.083476
LBFGS: 13 17:50:02 -143.956799 0.066261
LBFGS: 14 17:50:02 -143.958311 0.031498
-143.95831134909758
Final state#
# Set up your system as an ASE atoms object
final = fcc111("Pt", (3, 3, 3), vacuum=8, periodic=True)
adsorbate = molecule("O")
add_adsorbate(final, adsorbate, 2.0, "hcp")
final.calc = FAIRChemCalculator(predictor, task_name="oc20")
# Set up LBFGS dynamics object
opt = LBFGS(final)
opt.run(0.05, 100)
print(final.get_potential_energy())
Step Time Energy fmax
LBFGS: 0 17:50:02 -141.282602 3.340432
LBFGS: 1 17:50:03 -141.659778 3.323330
LBFGS: 2 17:50:03 -142.891396 2.596631
LBFGS: 3 17:50:03 -143.418904 1.225889
LBFGS: 4 17:50:03 -143.484098 0.977163
LBFGS: 5 17:50:03 -143.606345 0.136741
LBFGS: 6 17:50:03 -143.610846 0.118666
LBFGS: 7 17:50:03 -143.613317 0.100143
LBFGS: 8 17:50:03 -143.615008 0.078485
LBFGS: 9 17:50:03 -143.616457 0.051599
LBFGS: 10 17:50:03 -143.617140 0.033150
-143.61714008444915
Setup and relax the band#
from ase.mep import NEB
images = [initial]
for i in range(3):
image = initial.copy()
image.calc = FAIRChemCalculator(predictor, task_name="oc20")
images.append(image)
images.append(final)
neb = NEB(images)
neb.interpolate()
opt = LBFGS(neb, trajectory="neb.traj")
opt.run(0.05, 100)
Step Time Energy fmax
LBFGS: 0 17:50:04 -143.193978 3.039294
LBFGS: 1 17:50:04 -143.360197 1.460895
LBFGS: 2 17:50:04 -143.411350 0.450304
LBFGS: 3 17:50:05 -143.423402 0.447563
LBFGS: 4 17:50:05 -143.443052 0.476380
LBFGS: 5 17:50:05 -143.459315 0.378513
LBFGS: 6 17:50:05 -143.469906 0.211691
LBFGS: 7 17:50:06 -143.474781 0.177275
LBFGS: 8 17:50:06 -143.475806 0.183713
LBFGS: 9 17:50:06 -143.477287 0.178490
LBFGS: 10 17:50:07 -143.478786 0.167189
LBFGS: 11 17:50:07 -143.479388 0.094911
LBFGS: 12 17:50:07 -143.479539 0.096953
LBFGS: 13 17:50:07 -143.479789 0.100272
LBFGS: 14 17:50:08 -143.480349 0.124775
LBFGS: 15 17:50:08 -143.481008 0.092786
LBFGS: 16 17:50:08 -143.481431 0.050532
LBFGS: 17 17:50:08 -143.481734 0.040580
np.True_
from ase.mep import NEBTools
NEBTools(neb.images).plot_band();
This could be a good initial guess to initialize an NEB in DFT.
Ideas for things you can do with UMA#
FineTuna - use it for initial geometry optimizations then do DFT
a. https://iopscience.iop.org/article/10.1088/2632-2153/ac8fe0
b. https://iopscience.iop.org/article/10.1088/2632-2153/ad37f0
AdsorbML - prescreen adsorption sites to find relevant ones
a. https://www.nature.com/articles/s41524-023-01121-5
CatTsunami - screen NEBs more thoroughly
a. https://pubs.acs.org/doi/10.1021/acscatal.4c04272
Free energy estimations - compute vibrational modes and use them to estimate vibrational entropy
a. https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07477
Massive screening of catalyst surface properties (685M relaxations)
a. https://arxiv.org/abs/2411.11783
Advanced applications#
These take a while to run.
AdsorbML#
It is so cheap to run these calculations that we can screen a broad range of adsorbate sites and rank them in stability. The AdsorbML approach automates this. This takes quite a while to run here, and we don’t do it in the workshop.
Expert adsorption energies#
This tutorial reproduces Fig 6b from the following paper: Zhou, Jing, et al. “Enhanced Catalytic Activity of Bimetallic Ordered Catalysts for Nitrogen Reduction Reaction by Perturbation of Scaling Relations.” ACS Catalysis 134 (2023): 2190-2201 (https://doi.org/10.1021/acscatal.2c05877).
This takes up to an hour with a GPU, and much longer with a CPU.
CatTsunami#
The CatTsunami tutorial is an example of enumerating initial and final states, and computing reaction paths between them with UMA.
Acknowledgements#
This tutorial was originally compiled by John Kitchin (CMU) for the NAM29 catalysis tutorial session, using a variety of resources from the FAIR chemistry repository.