Adsorption Energies#
Pre-trained ODAC models are versatile across various MOF-related tasks. To begin, we’ll start with a fundamental application: calculating the adsorption energy for a single CO2 molecule. This serves as an excellent and simple demonstration of what you can achieve with these datasets and models.
For predicting the adsorption energy of a single CO2 molecule within a MOF structure, the adsorption energy (\(E_{\mathrm{ads}}\)) is defined as:
Each term on the right-hand side represents the energy of the relaxed state of the indicated chemical system. For a comprehensive understanding of our methodology for computing these adsorption energies, please refer to our paper.
Loading Pre-trained Models#
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 loginor by setting the HF_TOKEN environment variable.
A pre-trained model can be loaded using FAIRChemCalculator. In this example, we’ll employ UMA to determine the CO2 adsorption energies.
from fairchem.core import FAIRChemCalculator, pretrained_mlip
predictor = pretrained_mlip.get_predict_unit("uma-s-1p1")
calc = FAIRChemCalculator(predictor, task_name="odac")
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(
Adsorption in rigid MOFs: CO2 Adsorption Energy in Mg-MOF-74#
Let’s apply our knowledge to Mg-MOF-74, a widely studied MOF known for its excellent CO2 adsorption properties. Its structure comprises magnesium atomic complexes connected by a carboxylated and oxidized benzene ring, serving as an organic linker. Previous studies consistently report the CO2 adsorption energy for Mg-MOF-74 to be around -0.40 eV [1] [2] [3].
Our goal is to verify if we can achieve a similar value by performing a simple single-point calculation using UMA. In the ODAC23 dataset, all MOF structures are identified by their CSD (Cambridge Structural Database) code. For Mg-MOF-74, this code is OPAGIX. We’ve extracted a specific OPAGIX+CO2 configuration from the dataset, which exhibits the lowest adsorption energy among its counterparts.
import matplotlib.pyplot as plt
from ase.io import read
from ase.visualize.plot import plot_atoms
mof_co2 = read("structures/OPAGIX_w_CO2.cif")
mof = read("structures/OPAGIX.cif")
co2 = read("structures/co2.xyz")
fig, ax = plt.subplots(figsize=(5, 4.5), dpi=250)
plot_atoms(mof_co2, ax)
ax.set_axis_off()
The final step in calculating the adsorption energy involves connecting the FAIRChemCalculator to each relaxed structure: OPAGIX+CO2, OPAGIX, and CO2. The structures used here are already relaxed from ODAC23. For simplicity, we assume here that further relaxations can be neglected. We will show how to go beyond this assumption in the next section.
mof_co2.calc = calc
mof.calc = calc
co2.calc = calc
E_ads = (
mof_co2.get_potential_energy()
- mof.get_potential_energy()
- co2.get_potential_energy()
)
print(f"Adsorption energy of CO2 in Mg-MOF-74: {E_ads:.3f} eV")
Adsorption energy of CO2 in Mg-MOF-74: -0.459 eV
Adsorption in flexible MOFs#
The adsorption energy calculation method outlined above is typically performed with rigid MOFs for simplicity. Both experimental and modeling literature have shown, however, that MOF flexibility can be important in accurately capturing the underlying chemistry of adsorption [1] [2] [3]. In particular, uptake can be improved by treating MOFs as flexible. Two types of MOF flexibility can be considered: intrinsic flexibility and deformation induced by guest molecules. In the Open DAC Project, we consider the latter MOF deformation by allowing the atomic positions of the MOF to relax during geometry optimization [4]. The addition of additional degrees of freedoms can complicate the computation of the adsorption energy and necessitates an extra step in the calculation procedure.
The figure below shows water adsorption in the MOF with CSD code WOBHEB with added defects (WOBHEB_0.11_0) from a DFT simulation. A typical adsorption energy calculation would only seek to capture the effects shaded in purple, which include both chemisorption and non-bonded interactions between the host and guest molecule. When allowing the MOF to relax, however, the adsorption energy also includes the energetic effect of the MOF deformation highlighted in green.
To account for this deformation, it is vital to use the most energetically favorable MOF geometry for the empty MOF term in Eqn. 1. Including MOF atomic coordinates as degrees of freedom can result in three possible outcomes:
The MOF does not deform, so the energies of the relaxed empty MOF and the MOF in the adsorbed state are the same
The MOF deforms to a less energetically favorable geometry than its ground state
The MOF locates a new energetically favorable geoemtry relative to the empty MOF relaxation
The first outcome requires no additional computation because the MOF rigidity assumption is valid. The second outcome represents physical and reversible deformation where the MOF returns to its empty ground state upon removal of the guest molecule. The third outcome is often the result of the guest molecule breaking local symmetry. We also found cases in ODAC in which both outcomes 2 and 3 occur within the same MOF.
To ensure the most energetically favorable empty MOF geometry is found, an addition empty MOF relaxation should be performed after MOF + adsorbate relaxation. The guest molecule should be removed, and the MOF should be relaxed starting from its geometry in the adsorbed state. If all deformation is reversible, the MOF will return to its original empty geometry. Otherwise, the lowest energy (most favorable) MOF geometry should be taken as the reference energy, \(E_{\mathrm{MOF}}\), in Eqn. 1.
H2O Adsorption Energy in Flexible WOBHEB with UMA#
The first part of this tutorial demonstrates how to perform a single point adsorption energy calculation using UMA. To treat MOFs as flexible, we perform all calculations on geometries determined by geometry optimization. The following example corresponds to the figure shown above (H2O adsorption in WOBHEB_0.11_0).
In this tutorial, \(E_{x}(r_{y})\) corresponds to the energy of \(x\) determined from geometry optimization of \(y\).
First, we obtain the energy of the empty MOF from relaxation of only the MOF: \(E_{\mathrm{MOF}}(r_{\mathrm{MOF}})\)
import ase.io
from ase.optimize import BFGS
mof = ase.io.read("structures/WOBHEB_0.11.cif")
mof.calc = calc
relax = BFGS(mof)
relax.run(fmax=0.05)
E_mof_empty = mof.get_potential_energy()
print(f"Energy of empty MOF: {E_mof_empty:.3f} eV")
Step Time Energy fmax
BFGS: 0 17:47:38 -1077.274062 0.206406
BFGS: 1 17:47:38 -1077.276781 0.152729
BFGS: 2 17:47:38 -1077.281942 0.169930
BFGS: 3 17:47:39 -1077.284772 0.155755
BFGS: 4 17:47:39 -1077.288836 0.108762
BFGS: 5 17:47:39 -1077.291022 0.086447
BFGS: 6 17:47:39 -1077.293362 0.093427
BFGS: 7 17:47:39 -1077.295417 0.100117
BFGS: 8 17:47:40 -1077.297827 0.102542
BFGS: 9 17:47:40 -1077.300017 0.091601
BFGS: 10 17:47:40 -1077.302006 0.079013
BFGS: 11 17:47:40 -1077.304136 0.105577
BFGS: 12 17:47:41 -1077.306726 0.087913
BFGS: 13 17:47:41 -1077.309515 0.086368
BFGS: 14 17:47:41 -1077.312262 0.086817
BFGS: 15 17:47:41 -1077.314703 0.106260
BFGS: 16 17:47:41 -1077.316987 0.106230
BFGS: 17 17:47:42 -1077.319483 0.085492
BFGS: 18 17:47:42 -1077.322264 0.109624
BFGS: 19 17:47:42 -1077.325133 0.148696
BFGS: 20 17:47:42 -1077.327765 0.125955
BFGS: 21 17:47:43 -1077.329922 0.069111
BFGS: 22 17:47:43 -1077.331956 0.087269
BFGS: 23 17:47:43 -1077.334273 0.125245
BFGS: 24 17:47:43 -1077.336833 0.166687
BFGS: 25 17:47:43 -1077.339541 0.145561
BFGS: 26 17:47:44 -1077.342152 0.087671
BFGS: 27 17:47:44 -1077.344539 0.076165
BFGS: 28 17:47:44 -1077.346893 0.148926
BFGS: 29 17:47:44 -1077.349782 0.170227
BFGS: 30 17:47:45 -1077.352529 0.109269
BFGS: 31 17:47:45 -1077.354748 0.070345
BFGS: 32 17:47:45 -1077.356778 0.089681
BFGS: 33 17:47:45 -1077.358657 0.124294
BFGS: 34 17:47:45 -1077.360601 0.108050
BFGS: 35 17:47:46 -1077.362467 0.068582
BFGS: 36 17:47:46 -1077.364168 0.070192
BFGS: 37 17:47:46 -1077.365711 0.105448
BFGS: 38 17:47:46 -1077.367260 0.104214
BFGS: 39 17:47:47 -1077.368770 0.062853
BFGS: 40 17:47:47 -1077.370140 0.057033
BFGS: 41 17:47:47 -1077.371380 0.061028
BFGS: 42 17:47:47 -1077.372447 0.064131
BFGS: 43 17:47:47 -1077.373413 0.057878
BFGS: 44 17:47:48 -1077.374358 0.046757
Energy of empty MOF: -1077.374 eV
Next, we add the H2O guest molecule and relax the MOF + adsorbate to obtain \(E_{\mathrm{MOF+H2O}}(r_{\mathrm{MOF+H2O}})\).
mof_h2o = ase.io.read("structures/WOBHEB_H2O.cif")
mof_h2o.calc = calc
relax = BFGS(mof_h2o)
relax.run(fmax=0.05)
E_combo = mof_h2o.get_potential_energy()
print(f"Energy of MOF + H2O: {E_combo:.3f} eV")
Step Time Energy fmax
BFGS: 0 17:47:48 -1091.565588 1.145036
BFGS: 1 17:47:48 -1091.585062 0.314149
BFGS: 2 17:47:49 -1091.590211 0.243429
BFGS: 3 17:47:49 -1091.608169 0.237244
BFGS: 4 17:47:49 -1091.614632 0.227934
BFGS: 5 17:47:49 -1091.625217 0.186791
BFGS: 6 17:47:49 -1091.632353 0.178910
BFGS: 7 17:47:50 -1091.640634 0.175138
BFGS: 8 17:47:50 -1091.648041 0.184518
BFGS: 9 17:47:50 -1091.656146 0.160905
BFGS: 10 17:47:50 -1091.663844 0.178452
BFGS: 11 17:47:51 -1091.672298 0.188715
BFGS: 12 17:47:51 -1091.682088 0.157512
BFGS: 13 17:47:51 -1091.692983 0.177286
BFGS: 14 17:47:51 -1091.704436 0.158161
BFGS: 15 17:47:52 -1091.715509 0.191699
BFGS: 16 17:47:52 -1091.725709 0.197822
BFGS: 17 17:47:52 -1091.735328 0.163740
BFGS: 18 17:47:52 -1091.745541 0.151478
BFGS: 19 17:47:52 -1091.754026 0.170816
BFGS: 20 17:47:53 -1091.761503 0.153435
BFGS: 21 17:47:53 -1091.767895 0.152626
BFGS: 22 17:47:53 -1091.774206 0.165996
BFGS: 23 17:47:53 -1091.780885 0.135401
BFGS: 24 17:47:54 -1091.788354 0.181010
BFGS: 25 17:47:54 -1091.794284 0.204536
BFGS: 26 17:47:54 -1091.800638 0.131373
BFGS: 27 17:47:54 -1091.806519 0.189985
BFGS: 28 17:47:54 -1091.812301 0.199070
BFGS: 29 17:47:55 -1091.817180 0.151565
BFGS: 30 17:47:55 -1091.822210 0.100094
BFGS: 31 17:47:55 -1091.826304 0.125548
BFGS: 32 17:47:55 -1091.832529 0.177355
BFGS: 33 17:47:56 -1091.837114 0.246707
BFGS: 34 17:47:56 -1091.842045 0.112575
BFGS: 35 17:47:56 -1091.845776 0.328954
BFGS: 36 17:47:56 -1091.850879 0.174288
BFGS: 37 17:47:57 -1091.858477 0.159779
BFGS: 38 17:47:57 -1091.865291 0.141226
BFGS: 39 17:47:57 -1091.872011 0.142019
BFGS: 40 17:47:57 -1091.878340 0.245916
BFGS: 41 17:47:57 -1091.881002 0.540811
BFGS: 42 17:47:58 -1091.885458 0.690840
BFGS: 43 17:47:58 -1091.893765 0.213875
BFGS: 44 17:47:58 -1091.899099 0.157215
BFGS: 45 17:47:58 -1091.917187 0.281985
BFGS: 46 17:47:59 -1091.924934 0.354937
BFGS: 47 17:47:59 -1091.941727 0.250463
BFGS: 48 17:47:59 -1091.956888 0.455844
BFGS: 49 17:47:59 -1091.972389 0.822167
BFGS: 50 17:48:00 -1091.963293 1.583929
BFGS: 51 17:48:00 -1092.007673 0.458918
BFGS: 52 17:48:00 -1092.022540 0.288069
BFGS: 53 17:48:00 -1092.070710 0.353344
BFGS: 54 17:48:00 -1092.087621 0.303882
BFGS: 55 17:48:01 -1092.119305 0.356177
BFGS: 56 17:48:01 -1092.131400 0.260545
BFGS: 57 17:48:01 -1092.142927 0.386514
BFGS: 58 17:48:01 -1092.156453 0.280821
BFGS: 59 17:48:02 -1092.168005 0.246276
BFGS: 60 17:48:02 -1092.184048 0.339130
BFGS: 61 17:48:02 -1092.191921 0.372564
BFGS: 62 17:48:02 -1092.206071 0.334997
BFGS: 63 17:48:03 -1092.216222 0.284884
BFGS: 64 17:48:03 -1092.226534 0.169215
BFGS: 65 17:48:03 -1092.233212 0.113595
BFGS: 66 17:48:03 -1092.239250 0.106260
BFGS: 67 17:48:03 -1092.244869 0.111908
BFGS: 68 17:48:04 -1092.250320 0.121593
BFGS: 69 17:48:04 -1092.255307 0.108935
BFGS: 70 17:48:04 -1092.260573 0.087000
BFGS: 71 17:48:04 -1092.265096 0.096234
BFGS: 72 17:48:05 -1092.269023 0.127351
BFGS: 73 17:48:05 -1092.272618 0.153027
BFGS: 74 17:48:05 -1092.276079 0.173345
BFGS: 75 17:48:05 -1092.279381 0.129082
BFGS: 76 17:48:05 -1092.282878 0.100085
BFGS: 77 17:48:06 -1092.285440 0.090062
BFGS: 78 17:48:06 -1092.288010 0.090081
BFGS: 79 17:48:06 -1092.290418 0.087247
BFGS: 80 17:48:06 -1092.292323 0.067653
BFGS: 81 17:48:07 -1092.294055 0.077681
BFGS: 82 17:48:07 -1092.295419 0.076206
BFGS: 83 17:48:07 -1092.296890 0.074594
BFGS: 84 17:48:07 -1092.298274 0.074564
BFGS: 85 17:48:08 -1092.299959 0.084227
BFGS: 86 17:48:08 -1092.301584 0.087670
BFGS: 87 17:48:08 -1092.303368 0.078295
BFGS: 88 17:48:08 -1092.305103 0.066904
BFGS: 89 17:48:08 -1092.306538 0.055481
BFGS: 90 17:48:09 -1092.307932 0.054582
BFGS: 91 17:48:09 -1092.309096 0.065436
BFGS: 92 17:48:09 -1092.310101 0.063511
BFGS: 93 17:48:09 -1092.310900 0.052848
BFGS: 94 17:48:10 -1092.311633 0.045233
Energy of MOF + H2O: -1092.312 eV
We can now isolate the MOF atoms from the relaxed MOF + H2O geometry and see that the MOF has adopted a geometry that is less energetically favorable than the empty MOF by ~0.2 eV. The energy of the MOF in the adsorbed state corresponds to \(E_{\mathrm{MOF}}(r_{\mathrm{MOF+H2O}})\).
mof_adsorbed_state = mof_h2o[:-3]
mof_adsorbed_state.calc = calc
E_mof_adsorbed_state = mof_adsorbed_state.get_potential_energy()
print(f"Energy of MOF in the adsorbed state: {E_mof_adsorbed_state:.3f} eV")
Energy of MOF in the adsorbed state: -1077.091 eV
H2O adsorption in this MOF appears to correspond to Case #2 as outlined above. We can now perform re-relaxation of the empty MOF starting from the \(r_{\mathrm{MOF+H2O}}\) geometry.
relax = BFGS(mof_adsorbed_state)
relax.run(fmax=0.05)
E_mof_rerelax = mof_adsorbed_state.get_potential_energy()
print(f"Energy of re-relaxed empty MOF: {E_mof_rerelax:.3f} eV")
Step Time Energy fmax
BFGS: 0 17:48:10 -1077.090753 0.985834
BFGS: 1 17:48:10 -1077.123073 0.873012
BFGS: 2 17:48:10 -1077.172403 0.825194
BFGS: 3 17:48:11 -1077.211440 0.533393
BFGS: 4 17:48:11 -1077.231029 0.437351
BFGS: 5 17:48:11 -1077.247319 0.284757
BFGS: 6 17:48:11 -1077.258201 0.259521
BFGS: 7 17:48:11 -1077.267536 0.245625
BFGS: 8 17:48:12 -1077.277196 0.220215
BFGS: 9 17:48:12 -1077.284077 0.156490
BFGS: 10 17:48:12 -1077.289014 0.141604
BFGS: 11 17:48:12 -1077.292999 0.139771
BFGS: 12 17:48:12 -1077.296471 0.155922
BFGS: 13 17:48:13 -1077.301476 0.167136
BFGS: 14 17:48:13 -1077.305704 0.148967
BFGS: 15 17:48:13 -1077.309996 0.134975
BFGS: 16 17:48:13 -1077.313950 0.159379
BFGS: 17 17:48:14 -1077.318260 0.164480
BFGS: 18 17:48:14 -1077.322905 0.147942
BFGS: 19 17:48:14 -1077.326557 0.123462
BFGS: 20 17:48:14 -1077.329109 0.117042
BFGS: 21 17:48:14 -1077.331457 0.112932
BFGS: 22 17:48:15 -1077.334006 0.123251
BFGS: 23 17:48:15 -1077.336923 0.117101
BFGS: 24 17:48:15 -1077.339723 0.097952
BFGS: 25 17:48:15 -1077.342011 0.086996
BFGS: 26 17:48:16 -1077.344059 0.082096
BFGS: 27 17:48:16 -1077.345704 0.072120
BFGS: 28 17:48:16 -1077.347485 0.082814
BFGS: 29 17:48:16 -1077.349038 0.100270
BFGS: 30 17:48:17 -1077.350775 0.078745
BFGS: 31 17:48:17 -1077.352068 0.052648
BFGS: 32 17:48:17 -1077.353380 0.063390
BFGS: 33 17:48:17 -1077.354704 0.080220
BFGS: 34 17:48:17 -1077.356125 0.087689
BFGS: 35 17:48:18 -1077.357774 0.074546
BFGS: 36 17:48:18 -1077.359291 0.071444
BFGS: 37 17:48:18 -1077.360566 0.073488
BFGS: 38 17:48:18 -1077.361785 0.070318
BFGS: 39 17:48:19 -1077.363226 0.086594
BFGS: 40 17:48:19 -1077.364365 0.069672
BFGS: 41 17:48:19 -1077.365477 0.051307
BFGS: 42 17:48:19 -1077.366552 0.063514
BFGS: 43 17:48:19 -1077.367685 0.069021
BFGS: 44 17:48:20 -1077.369029 0.080557
BFGS: 45 17:48:20 -1077.370319 0.049330
Energy of re-relaxed empty MOF: -1077.370 eV
The MOF returns to its original empty reference energy upon re-relaxation, confirming that this deformation is physically relevant and is induced by the adsorbate molecule. In Case #3, this re-relaxed energy will be more negative (more favorable) than the original empty MOF relaxation. Thus, we take the reference empty MOF energy (\(E_{\mathrm{MOF}}\) in Eqn. 1) to be the minimum of the original empty MOF energy and the re-relaxed MOf energy:
E_mof = min(E_mof_empty, E_mof_rerelax)
# get adsorbate reference energy
h2o = mof_h2o[-3:]
h2o.calc = calc
E_h2o = h2o.get_potential_energy()
# compute adsorption energy
E_ads = E_combo - E_mof - E_h2o
print(f"Adsorption energy of H2O in WOBHEB_0.11_0: {E_ads:.3f} eV")
Adsorption energy of H2O in WOBHEB_0.11_0: -0.689 eV
This adsorption energy closely matches that from DFT (–0.699 eV) [1]. The strong adsorption energy is a consequence of both H2O chemisorption and MOF deformation. We can decompose the adsorption energy into contributions from these two factors. Assuming rigid H2O molecules, we define \(E_{\mathrm{int}}\) and \(E_{\mathrm{MOF,deform}}\), respectively, as
\(E_{\mathrm{int}}\) describes host host–guest interactions for the MOF in the adsorbed state only. \(E_{\mathrm{MOF,deform}}\) quantifies the magnitude of deformation between the MOF in the adsorbed state and the most energetically favorable empty MOF geometry determined from the workflow presented here. It can be shown that
For H2O adsorption in WOBHEB_0.11, we have
E_int = E_combo - E_mof_adsorbed_state - E_h2o
print(f"E_int: {E_int}")
E_int: -0.9724213480947537
E_mof_deform = E_mof_adsorbed_state - E_mof_empty
print(f"E_mof_deform: {E_mof_deform}")
E_mof_deform: 0.28360462188720703
E_ads = E_int + E_mof_deform
print(f"E_ads: {E_ads}")
E_ads: -0.6888167262075466
\(E_{\mathrm{int}}\) is equivalent to \(E_{\mathrm{ads}}\) when the MOF is assumed to be rigid. In this case, failure to consider adsorbate-induced deformation would result in an overestimation of the adsorption energy magnitude.