Widom insertion

Breafly introduction about Widom insertion method

The Widom insertion method is a powerful and computationally efficient technique in statistical mechanics used to calculate the excess chemical potential of a species at infinite dilution. The method operates by inserting a “ghost” or “test” particle at random positions and orientations within a static configuration of a host system, such as a porous material. For each insertion, the interaction energy between the ghost particle and the host is calculated, but the particle is not actually added, so the host’s configuration remains unchanged. By averaging the Boltzmann factor of these interaction energies over thousands or millions of trials, one can directly compute the Henry’s constant, which is fundamentally related to the material’s affinity for the adsorbate at low pressures, and the Isosteric Heat of Adsorption at zero coverage (Qst). Furthermore, the distribution of insertion energies reveals the potential energy landscape, identifying the coordinates of the most stable adsorption sites within the framework.

This makes the Widom method arguably the easiest and most direct way to test a new Machine Learning Potential (MLP) for adsorption applications. Unlike full Grand Canonical Monte Carlo simulations, it does not require lengthy equilibration or the simulation of multiple interacting guest molecules. The method yields two crucial, physically meaningful metrics: the most stable binding configurations and the enthalpy of adsorption. These two outputs can be directly compared with high-fidelity experimental data from microcalorimetry and diffraction techniques, or with results from expensive quantum mechanical calculations, providing a robust and computationally inexpensive first-pass validation of the MLP’s accuracy in describing host-guest interactions.

---

Setting the stage

In this example we will evaluate the Mg-MOF-74, a material well-known for having adsorption of CO ₂ on Open-Metal sites, which is very hard to model with classical force fields and usually require electronic structure methods, such DFT, to describe the interaction properly.

On the flames/examples/Basic/1-Widom folder you will find the cif file for the Mg-MOF-74 [1]. The simulation will be executed with the machine learning potential MACE [2], using the medium-0b2 pre-trained foundation model with D3(0) dispersion correction.

Warning

This model is used here just as an example, you should look carefully into the mace documentation to undestand how the potential work and which one is the best for your specific case.

The script below will run the simulation at 298 K (25°C).

import os

# Hide UserWarning and RuntimeWarning messages
import warnings

import ase
import torch
from ase.data import vdw_radii
from ase.io import read
from mace.calculators import mace_mp

from flames.widom import Widom

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)

os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"

device = "cuda" if torch.cuda.is_available() else "cpu"

model = mace_mp(
    model="medium-0b2",
    dispersion=True,
    damping="zero",  # choices: ["zero", "bj", "zerom", "bjm"]
    dispersion_xc="pbe",
    default_dtype="float32",
    device=device,
)

FrameworkPath = "mg-mof-74.cif"
AdsorbatePath = "co2.xyz"

# Load the framework structure
framework: ase.Atoms = read(FrameworkPath)  # type: ignore

# Load the adsorbate structure
adsorbate: ase.Atoms = read(AdsorbatePath)  # type: ignore

Temperature = 298.0

NSteps = 3000

widom = Widom(
    model=model,
    framework_atoms=framework,
    adsorbate_atoms=adsorbate,
    temperature=Temperature,
    device=device,
    vdw_radii=vdw_radii,
    debug=False,
    output_to_file=True,
    random_seed=42,
    cutoff_radius=6.0,
    automatic_supercell=True,
)

widom.logger.print_header()

widom.run(NSteps)
widom.logger.print_summary()
widom.save_results()

---

Breaking down the input script

Importing libraries

import os

# Hide UserWarning and RuntimeWarning messages
import warnings

import ase
import torch
from ase.data import vdw_radii
from ase.io import read
from mace.calculators import mace_mp

from flames.widom import Widom

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)

os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"

This first part only import the libraries and make the code to ignore a few non-important warnings from mace. The os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" command is ony important on Windows machines. If you are running on Linux, this command will just be ignored.

Here the vdw_radii from ase.data is imported to calculate the smaller distance allowed between two atoms.

---

Defining the calculator

device = "cuda" if torch.cuda.is_available() else "cpu"

model = mace_mp(
    model="medium-0b2",
    dispersion=True,
    damping="zero",  # choices: ["zero", "bj", "zerom", "bjm"]
    dispersion_xc="pbe",
    default_dtype="float32",
    device=device,
)

This is an extremelly important part, as the calculator will be used to calculate the energy of interaction between the adsorbate (CO2) and the adsorbent (Mg-MOF-74). If you are running on GPU (which is extremelly recommended) you need to define the device correctly, so the calculator use the GPU. Currently, basically all MLPs only use a single GPU.

Here the pre-trained model medium-0b2 is used, but currently there are several other variations available. For a complete list, check the mace-foundation repository: https://github.com/ACEsuit/mace-foundations.

For a list of available models, check https://github.com/ACEsuit/mace/blob/100a29149d90a5945eddec1f0940bd88a6e3b363/mace/calculators/foundations_models.py#L19. You can also run a local model, just passing the path of the model instead of the name to the variable model.

Reading the structure

FrameworkPath = "mg-mof-74.cif"
AdsorbatePath = "co2.xyz"

# Load the framework structure
framework: ase.Atoms = read(FrameworkPath)  # type: ignore

# Load the adsorbate structure
adsorbate: ase.Atoms = read(AdsorbatePath)  # type: ignore

This part of the code will read the framework and the adsorbate files. It uses the read function from ase, which is very powerful and can read most file formats. Here, just make sure that the framework file is in a file format that is compatible with periodic boundary conditions and has a unt cell correctly defined.

In case you want to make sure that your framework was read properly, you can create a separate python script and visualize the structure with the ase visualizator:

from ase.io import read
from ase.visualize import view

FrameworkPath = "mg-mof-74.cif"

# Load the framework structure
framework: ase.Atoms = read(FrameworkPath)  # type: ignore

# Visualize the framework structure
view(framework)

Running the simulation

Temperature = 298.0
NSteps = 3000

widom = Widom(
    model=model,
    framework_atoms=framework,
    adsorbate_atoms=adsorbate,
    temperature=Temperature,
    device=device,
    vdw_radii=vdw_radii,
    debug=False,
    output_to_file=True,
    random_seed=42,
    cutoff_radius=6.0,
    automatic_supercell=True,
)

widom.logger.print_header()

widom.run(NSteps)
widom.logger.print_summary()
widom.save_results()

Finally, this block defines the widom simulator and runs the code. You can define the Temperature and the number of Monte Carlo steps NSteps.

Here we are using 3.000 for example, but it is recommended that you use way more steps for a serious simulations, and check if in fact your results are equilibrated. Usually 100.000 is a good starting point, but it strongly depends on the system.

You should change the cutoff_radius according to your system and potential used, usually 6.0 Å is a good starting point but check carefully the cutoff of your potential. The automatic_supercell option will create a supercell big enough to fit the cutoff radius, so you don’t need to worry about that.

Analyzing the output

The header

===========================================================================
         _______  __          ___      .___  ___.  _______     _______.
        |   ____||  |        /   \     |   \/   | |   ____|   /       |
        |  |__   |  |       /  ^  \    |  \  /  | |  |__     |   (----`
        |   __|  |  |      /  /_\  \   |  |\/|  | |   __|     \   \
        |  |     |  `----./  _____  \  |  |  |  | |  |____.----)   |
        |__|     |_______/__/     \__\ |__|  |__| |_______|_______/

        Flexible Lattice Adsorption by Monte Carlo Engine Simulation
                        powered by Python + ASE
                    Author: Felipe Lopes de Oliveira
===========================================================================

Code version: 0.4.4
Simulation started at 2025-11-17 23:06:46
Hostname: felipe-Dell-G16-7630
OS type: Linux
OS release: 6.7.10-060710-generic
OS version: #202403151538 SMP PREEMPT_DYNAMIC Fri Mar 15 20:14:18 UTC 2024

Python version: 3.12.11
Numpy version: 1.26.4
ASE version: 3.25.0

Current directory: /home/felipe/PRs/flames/examples/Basic/1-Widom
Random Seed: 42

Model: sumcalculator
Running on device: cuda

The header will print some information about the code version, the host machine, the python and ASE versions, the current directory, the random seed used, the model and the device used. It is always important to check this information to make sure everything is correct.

The simulation summary

===========================================================================

Constants used:
Boltzmann constant:     8.617330337217213e-05 eV/K
Beta (1/kT):            38.941 eV^-1
Fugacity coefficient:   0.000000000 (dimensionless)

===========================================================================

Simulation Parameters:
Temperature: 298.0 K
Pressure: 0.00000 bar
Fugacity: 0.000 Pa
Fugacity: 0.00000e+00 eV/m^3
(1/kB.T) * V * f = 0.0 [-]

===========================================================================

This block presents the constants used in the simulation, mostly for reference and debugging purposes. It also summarizes the simulation parameters, including temperature, pressure, and fugacity values. For Widom simulations, there is no pressure or fugacity involved, so these values are zero.

The system information

System Information:
Framework: C144H36Mg36O108
Framework: 324 atoms,
Framework mass: 4368.743999999999 g/mol, 7.254469967765757e-24 kg
Framework energy: -2365.8930653513903 eV
Framework volume: 8.24197326482685e-27 m^3
Framework density: 880.186059171615 kg/m^3, 0.8801860591716151 g/cm^3
Framework cell:
    26.1701000    0.0000000    0.0000000
    13.0850500   22.6639714    0.0000000
    0.0000000    0.0000000   13.8960000

Perpendicular cell:
    22.6639714    0.0000000    0.0000000
    0.0000000   22.6639714    0.0000000
    0.0000000    0.0000000   13.8960000

Atomic positions:
Mg   16.5873945   14.4378564    3.7075223
O     4.8947247    5.7482631    2.3059717
C     5.6694905    4.8079349    1.7135852
...

 H    25.8864161   22.0824139    8.0335555

===========================================================================
Adsorbate: CO2
Adsorbate: 3 atoms, 7.307866261136e-26 kg
Adsorbate energy: -22.788235797845456 eV

Atomic positions:
C     0.0000000    0.0000000    0.0000000
O     1.1606430    0.0000000    0.0000000
O    -1.1606430    0.0000000    0.0000000

===========================================================================
Shortest distances:
O  - H : 1.632 Å
O  - Mg: 1.950 Å
O  - C : 1.932 Å
H  - Mg: 1.758 Å
H  - C : 1.740 Å
Mg - C : 2.058 Å

===========================================================================
Conversion factors:
    Conversion factor molecules/unit cell -> mol/kg:         0.228898741
    Conversion factor molecules/unit cell -> mg/g:           10.073604679
    Conversion factor molecules/unit cell -> cm^3 STP/gr:    5.130527703
    Conversion factor molecules/unit cell -> cm^3 STP/cm^3:  4.515818960
    Conversion factor molecules/unit cell -> %wt:            1.007360468

Partial pressure:
                0.000000000000000 Pascal
                0.000000000000000 bar
                0.000000000000000 atm
                0.000000000000000 Torr
===========================================================================

This section provides detailed information about the framework and adsorbate, including their compositions, masses, energies, volumes, densities, and atomic positions.

It also lists the shortest distances between different atom types in the system, which can be useful for understanding potential interactions. Always check these distances to make sure there are no unphysical overlaps or if it does not prevent desired interactions.

Additionally, conversion factors are provided to facilitate the interpretation of simulation results in various units. For Widom simulations these factores are not used, but they are printed for reference.

The simulation

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Starting Widom simulation
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Iteration  |  dE (eV)  |  dE (kJ/mol)  | kH [mol kg-1 Pa-1]  |  dH (kJ/mol) | Time (s)
---------------------------------------------------------------------------------------
    1      | -0.049759 |         -4.80 |           1.821e-06 |        -6.67 |     0.52
    2      | -0.118381 |        -11.42 |           1.657e-05 |       -13.37 |     0.61
    3      | -0.156906 |        -15.14 |           6.406e-05 |       -16.79 |     0.15
    4      | -0.031338 |         -3.02 |           5.156e-05 |       -16.72 |     0.15
    5      | -0.034721 |         -3.35 |           4.326e-05 |       -16.65 |     0.30
    6      | -0.117306 |        -11.32 |           4.339e-05 |       -16.24 |     0.15
    7      | -0.077221 |         -7.45 |           3.913e-05 |       -16.05 |     0.15
    8      | -0.076447 |         -7.38 |           3.578e-05 |       -15.88 |     0.30
    9      | -0.077040 |         -7.43 |           3.312e-05 |       -15.71 |     0.15
    10     | -0.142258 |        -13.73 |           4.073e-05 |       -15.84 |     0.15
    11     | -0.184081 |        -17.76 |           8.692e-05 |       -18.35 |     0.17
    12     | -0.074968 |         -7.23 |           8.089e-05 |       -18.28 |     0.15
    13     | -0.233352 |        -22.52 |           3.646e-04 |       -23.61 |     0.15
    14     | -0.084846 |         -8.19 |           3.412e-04 |       -23.58 |     0.16
    15     | -0.068756 |         -6.63 |           3.203e-04 |       -23.56 |     0.31
    16     | -0.313437 |        -30.24 |           5.694e-03 |       -32.23 |     0.15
...

Here the main simulation output is presented. For each Widom insertion step, the change in energy (dE) upon inserting the ghost particle is reported in both eV and kJ/mol. The Henry’s constant (kH) and isosteric heat of adsorption at zero coverage (dH) are calculated based on the average Boltzmann factor of the insertion energies up to that point. The time taken for each insertion step is logged, allowing users to monitor performance.

The results

...
2998    | -0.083372 |         -8.04 |           2.046e+00 |       -53.78 |     0.43
2999    | -0.146298 |        -14.12 |           2.045e+00 |       -53.78 |     0.42
3000    | -0.156466 |        -15.10 |           2.044e+00 |       -53.78 |     0.33

===========================================================================

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Finishing Widom simulation
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    Average properties of the system:
    ------------------------------------------------------------------------------
    Henry coefficient: [mol/kg/Pa]       2.04418e+00 +/-  3.33591e-01 [-]
    Enthalpy of adsorption: [kJ/mol]       -53.78391 +/-      0.34888 [-]

===========================================================================
Simulation finished successfully!
===========================================================================

Simulation finished at 2025-11-17 23:29:43
Simulation duration: 0:13:37.214066
===========================================================================

Finally the simulation result is printed, with the average values reported.

Also, a json file is saved with the main results, that can be easily read and processed for further analysis.

{
    "code_version": "0.4.4",
    "random_seed": 42,
    "enlapsed_time_hours": 0.22700391277777776,
    "total_insertions": 3001,
    "temperature_K": 298.0,
    "henry_coefficient_mol_kg-1_Pa-1": 2.044178342815808,
    "henry_coefficient_std_mol_kg-1_Pa-1": 0.3335905351383779,
    "enthalpy_of_adsorption_kJ_mol-1": -53.78390928504431,
    "enthalpy_of_adsorption_std_kJ_mol-1": 0.3488799686455756
}

References

[1]

Mason, J.A., Sumida, K., Herm, Z.R., Krishna, R. and Long, J.R., 2011. Evaluating metal-organic frameworks for post-combustion carbon dioxide capture via temperature swing adsorption. Energy & Environmental Science, 4(8), pp.3030-3040. https://pubs.rsc.org/en/content/articlelanding/2011/ee/c1ee01720a

[2]

Batatia, I., Benner, P., Chiang, Y., Elena, A.M., Kovács, D.P., Riebesell, J., Advincula, X.R., Asta, M., Avaylon, M., Baldwin, W.J. and Berger, F., 2025. A foundation model for atomistic materials chemistry. The Journal of Chemical Physics, 163(18). https://pubs.aip.org/aip/jcp/article/163/18/184110/3372267/A-foundation-model-for-atomistic-materials