GCMC Simulation

Brief introduction about GCMC simulation

Grand Canonical Monte Carlo (GCMC) is a stochastic simulation technique performed in the μVT ensemble, meaning the system’s chemical potential (μ) (usually through a fixed pressure, p), volume (V), and temperature (T) are held constant. Unlike conventional molecular simulations that maintain a fixed number of particles, GCMC introduces dynamic insertion and deletion moves to simulate the exchange of molecules with an infinite external reservoir. This open-system approach makes GCMC the gold standard for modeling gas adsorption in nanoporous materials, such as zeolites, porous carbons, and metal-organic frameworks (MOFs) and covalent organic frameworks (COFs). During these simulations, the porous material is typically modeled as a rigid host lattice, and the algorithm determines the equilibrium density of guest molecules (the adsorbate) that populate the pores at a given temperature and pressure.

---

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/3-Rigid_GCMC 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) and 1 bar.

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.gcmc import GCMC

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

TTemperature = 298.0  # in Kelvin
pressure = 100_000  # in Pa = 1 bar
MCSteps = 30_000

gcmc = GCMC(
    model=model,
    framework_atoms=framework,
    adsorbate_atoms=adsorbate,
    temperature=Temperature,
    pressure=pressure,
    device=device,
    vdw_radii=vdw_radii,
    vdw_factor=0.6,
    save_frequency=1,
    debug=False,
    output_to_file=True,
    criticalTemperature=304.1282,
    criticalPressure=7377300.0,
    acentricFactor=0.22394,
    random_seed=42,
    cutoff_radius=6.0,
    automatic_supercell=True,
)

gcmc.logger.print_header()

gcmc.run(MCSteps)
gcmc.logger.print_summary()
gcmc.save_results()

---

Breaking down the input script

The first block of the code is similar to the one used in the Widom insertion tutorial, where we import the necessary libraries, set up the device, load the model and read the framework and adsorbate structures.

Running the simulation

Temperature = 298.0
pressure = 100_000  # in Pa = 1 bar
MCSteps = 30_000

gcmc = GCMC(
    model=model,
    framework_atoms=framework,
    adsorbate_atoms=adsorbate,
    temperature=Temperature,
    pressure=pressure,
    device=device,
    vdw_radii=vdw_radii,
    vdw_factor=0.6,
    save_frequency=1,
    debug=False,
    output_to_file=True,
    criticalTemperature=304.1282,
    criticalPressure=7377300.0,
    acentricFactor=0.22394,
    random_seed=42,
    cutoff_radius=6.0,
    automatic_supercell=True,
)

gcmc.logger.print_header()

gcmc.run(MCSteps)
gcmc.logger.print_summary()
gcmc.save_results()

This block defines the gcmc simulator and runs the code. You can define the Temperature, pressure, and the number of Monte Carlo steps MCSteps.

Warning

Here we are using 100.000 steps as an example, but it is recommended that evaluate the best value for the simulations you are interested in, and check if in fact your results are equilibrated. Usually 100.000 is a good starting point, but it strongly depends on the system.

Note

FLAMES uses MC steps as basic step unit, which means that in each step the code will attempt either an insertion, deletion, or translation/rotation/reinsertion move. This differs from the implementation on other codes, such as RASPA, that uses cycles as basic step unit, where in each cycle the code will attempt N moves, being N=max(100, number of molecules in the system). So if you want to compare results with RASPA, make sure to use the same number of moves, which means that if you have 100 molecules in the system, you should use 100 MC steps in FLAMES to be equivalent to 1 cycle in RASPA.

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 twice the cutoff radius.

The criticalTemperature, criticalPressure, and acentricFactor are used to calculate the fugacity coefficient, which is used to convert the pressure to fugacity, that is the correct thermodynamic variable to be used in GCMC simulations. If no value is provided, the code will assume ideal gas behavior and use the pressure as fugacity, which is a good approximation at low pressures but can lead to significant errors at high pressures.

Analyzing the output

The system information

System Information:
Framework: C144H36Mg36O108
Framework: 324 atoms,
Framework mass: 4368.743999999999 g/mol, 7.254469967765757e-24 kg
Framework energy: -2351.243408203125 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
C     6.1334863    3.7363823    2.4891905
...
C    25.9933210   22.3518885    9.0813834
H    25.8864161   22.0824139    8.0335555

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

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

===========================================================================
Equation of State Parameters:

    Critical temparure [K]: 304.128200
    Critical pressure [Pa]: 7377300.000000
    Acentric factor [-]:    0.223940

    Vapour=stable, Liquid=metastable

    MolFraction:           1.0000000000 [-]
    Compressibility:       0.994513 [-]
    Fugacity coeff.:       0.9945357225 [-]
    Bulk phase pressure:   99453.572250 [Pa]

    Density of the bulk fluid phase:      1.785996 [kg/m^3]

    Amount of excess molecules:        0.0000000000 [-]


===========================================================================
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:
            99453.572250000000000 Pascal
                0.994535722500000 bar
                0.981530444000000 atm
                0.001291487000000 Torr
===========================================================================

Movement statistics:
    Insertion  : 0.200
    Deletion   : 0.200
    Translation: 0.200
    Rotation   : 0.200
    Reinsertion: 0.200

The header and simulation summary are the same as in the Widom insertion output. The system information section, however, now shows the Equation of State parameter data and the movement weights used on the simulation. The values are normalizes so the sum is 1.

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 GCMC simulation
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 Iteration |  Number of  |  Uptake  |    Tot En.   |Av. Ads. En.|  Pacc  |  Pdel  |  Ptra  |  Prot  |  Prin  |  Time
     -     |  Molecules  | [mmol/g] |     [eV]     |  [kJ/mol]  |    %   |    %   |   %    |   %    |   %    |   [s]
---------- | ----------- | -------- | ------------ | ---------- | ------ | ------ | ------ | ------ | ------ | -------
     1     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.31
     2     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.29
     3     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.00
     4     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.41
     5     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.01
     6     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.00
     7     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.01
     8     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.01
     9     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.00
    10     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.09
    11     |      0      |     0.00 |   -2351.2434 |     0.0000 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.01
    12     |      1      |     0.23 |   -2374.0972 |    -6.5381 |   8.33 |   0.00 |   0.00 |   0.00 |   0.00 |     0.09
    13     |      1      |     0.23 |   -2374.1006 |    -6.8679 |   8.33 |   0.00 |   0.00 | 100.00 |   0.00 |     0.09
    14     |      0      |     0.00 |   -2351.2432 |     0.0000 |   8.33 | 100.00 |   0.00 | 100.00 |   0.00 |     0.08
    15     |      0      |     0.00 |   -2351.2432 |     0.0000 |   7.69 | 100.00 |   0.00 | 100.00 |   0.00 |     0.00
    16     |      0      |     0.00 |   -2351.2432 |     0.0000 |   7.14 | 100.00 |   0.00 | 100.00 |   0.00 |     0.00
    17     |      0      |     0.00 |   -2351.2432 |     0.0000 |   6.67 | 100.00 |   0.00 | 100.00 |   0.00 |     0.11
...

Here the main simulation output is presented. For each GCMC step, the uptake (in mmol/g) total energy and average adsorption energy (in kJ/mol) are reported. THe acceptance probabilities for each type of move are also printed, along with the time taken for each step.

The results

...
179989   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.00
179990   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.97
179991   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     1.00
179992   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.99
179993   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     1.00
179994   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     1.02
179995   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.99
179996   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.00
179997   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.97
179998   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     1.00
179999   |     33      |     7.55 |   -3116.1934 |   -38.0439 |   0.00 |   0.00 |   0.00 |   0.00 |   0.00 |     0.97
180000   |     33      |     7.55 |   -3116.1948 |   -38.0481 |   0.00 |   0.00 |   0.00 |  14.29 |   0.00 |     0.99

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

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Finishing GCMC simulation
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    pyMSER Equilibration Results:
    ------------------------------------------------------------------------------
    Start of equilibrated data:          53850 of 180000
    Total equilibrated steps:            125653  (70.00%)
    Equilibrated:                        True
    Average over equilibrated data:      33.3205 ± 0.9215 molecules/unit cell
    Number of uncorrelated samples:      39.8
    Autocorrelation time:                3155.0
    ------------------------------------------------------------------------------

    Average properties of the system:
    ------------------------------------------------------------------------------
    Average loading absolute [molecules/unit cell]           33.32055 +/-      0.92151 [-]
    Average loading absolute [mol/kg framework]               7.62703 +/-      0.21093 [-]
    Average loading absolute [mg/g framework]               335.65804 +/-      9.28296 [-]
    Average loading absolute [cm^3 (STP)/gr framework]      170.95200 +/-      4.72785 [-]
    Average loading absolute [cm^3 (STP)/cm^3 framework]    150.46957 +/-      4.16139 [-]
    Average loading absolute [%wt framework]                 33.56580 +/-      0.92830 [-]

    Average excess absolute [molecules/unit cell]            33.32055 +/-      0.92151 [-]
    Average loading absolute [mol/kg framework]               7.62703 +/-      0.21093 [-]
    Average loading absolute [mg/g framework]               335.65804 +/-      9.28296 [-]
    Average loading absolute [cm^3 (STP)/gr framework]      170.95200 +/-      4.72785 [-]
    Average loading absolute [cm^3 (STP)/cm^3 framework]    150.46957 +/-      4.16139 [-]
    Average loading absolute [%wt framework]                 33.56580 +/-      0.92830 [-]


    Enthalpy of adsorption: [kJ/mol]                        -21.25549 +/-      3.05431 [kJ/mol]

===========================================================================
GCMC simulation finished successfully!
===========================================================================

Simulation finished at 2026-03-11 02:48:41
Simulation duration: 6:32:28.669883
===========================================================================

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

Note

FLAMES uses the pyMSER library to analyze the equilibration of the system, and report the average values calculated for the equilibrated portion of the simulation. Check the pyMSER paper for more details about the method and how to interpret the results.

Also, a json file with the main results is saved with the name results_[TEMPERATURE]_[PRESSURE].json, that can be easily read and processed for further analysis.

{
    "simulation": {
        "code_version": "0.4.5A",
        "random_seed": 42,
        "temperature_K": 298.0,
        "pressure_Pa": 100000,
        "fugacity_coefficient": 0.9945357224994603,
        "fugacity_Pa": 99453.57224994604,
        "move_weights": {
            "insertion": 0.2,
            "deletion": 0.2,
            "translation": 0.2,
            "rotation": 0.2,
            "reinsertion": 0.2
        },
        "n_steps": 179503,
        "enlapsed_time_hours": 0.004653468055555556
    },
    "equilibration": {
        "LLM": true,
        "t0": 53850,
        "average": 33.3205494496749,
        "uncertainty": 0.9215130805969238,
        "equilibrated": true,
        "ac_time": 3155,
        "uncorr_samples": 39
    },
    "enthalpy": {
        "kJ_mol": {
            "mean": -21.25548795467983,
            "sd": 3.054310813299496
        }
    },
    "absolute_uptake": {
        "nmol": {
            "mean": 33.3205494496749,
            "sd": 0.9215130805969238
        },
        "mol/kg": {
            "mean": 7.627031810114402,
            "sd": 0.21093318373289396
        },
        "mg/g": {
            "mean": 335.6580428449786,
            "sd": 9.282958480512942
        },
        "cm^3 STP/gr": {
            "mean": 170.95200203068663,
            "sd": 4.727848388677955
        },
        "cm^3 STP/cm^3": {
            "mean": 150.469568974888,
            "sd": 4.1613862415913205
        },
        "% wt": {
            "mean": 33.565804284497865,
            "sd": 0.9282958480512944
        }
    },
    "excess_uptake": {
        "nmol": {
            "mean": 33.3205494496749,
            "sd": 0.9215130805969238
        },
        "mol/kg": {
            "mean": 7.627031810114402,
            "sd": 0.21093318373289396
        },
        "mg/g": {
            "mean": 335.6580428449786,
            "sd": 9.282958480512942
        },
        "cm^3 STP/gr": {
            "mean": 170.95200203068663,
            "sd": 4.727848388677955
        },
        "cm^3 STP/cm^3": {
            "mean": 150.469568974888,
            "sd": 4.1613862415913205
        },
        "% wt": {
            "mean": 33.565804284497865,
            "sd": 0.9282958480512944
        }
    }
}

Finally, it is always recommended to check the convergence of the results, by plotting the uptake as a function of the simulation steps, and check if it is equilibrated and if the results are stable. You can use the results_[TEMPERATURE]_[PRESSURE].json file to plot the results with your preferred plotting library, such as matplotlib, seaborn, or even excel.

Below there is the uptake as a function of the simulation steps. It is possible to see that the system is equilibrated after around 50.000 steps, as detected by pyMSER, allowing the precise separation between equilibration steps and production steps.

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