=============== 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 :sub:`2` 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. .. code:: python 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 ---------------------- .. code:: python 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 ---------------------- .. code:: none 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 -------------- .. code:: none +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 ----------- .. code:: none ... 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. .. code:: python { "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. .. raw:: html
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