Need Advice on GCMC Simulation Settings for MOFs

[AubrielleCruz]

Hi everyone, war the knights

I'm currently running a GCMC simulation for gas adsorption in a metal–organic framework (MOF), but I'm unsure whether my simulation parameters are appropriate. The adsorption results seem lower than expected, and I'd like to verify that my setup is correct.

Has anyone experienced something similar? Are there any key parameters or best practices I should check to improve the accuracy of my simulation?

Any suggestions would be greatly appreciated.

Thanks!

[dubbelda]

This is very common, and "too low" is usually a comparison problem (units, absolute vs excess, unfinished sampling) or a setup problem (no working swap moves, missing charges, force field / cutoff / box size). A single parameter almost never explains it on its own.

RASPA2 uses simulation.input (not JSON). A minimal GCMC setup looks like:

SimulationType MonteCarlo
NumberOfCycles 50000
NumberOfInitializationCycles 10000
PrintEvery 1000

Forcefield YourForceField
CutOffVDW 12.0

Framework 0
FrameworkName YourMOF
UnitCells 2 2 2
HeliumVoidFraction 0.7
ExternalTemperature 298.0
ExternalPressure 1e5
ChargeMethod Ewald
UseChargesFromCIFFile yes

Component 0 MoleculeName CO2
 MoleculeDefinition TraPPE
 TranslationProbability 0.5
 RotationProbability 0.5
 ReinsertionProbability 0.5
 SwapProbability 1.0
 CreateNumberOfMolecules 0


1. Confirm you are actually running GCMC, and that insertions work

You need insertion/deletion: SwapProbability (CBMC). Translation / rotation / reinsertion alone is NVT: the number of molecules never changes.

In the output, check the swap statistics. If insertions are almost never constructed or accepted, the loading will sit too low. That is typical at high loading or in tight pores with ordinary CBMC. Switch to CFCMC (CFCMCProbability and/or CBCFCMCProbability) and give NumberOfEquilibrationCycles so the lambda biasing can be learned.

Also check that the five block averages of the loading have actually flattened. If they are still climbing, the run is not equilibrated.

CreateNumberOfMolecules 0 is correct for GCMC starting from an empty framework. Starting with a guessed loading is optional; it does not replace swap moves.


2. Check units and which loading you are comparing

- ExternalPressure is in pascal. 1 bar is 1e5, not 1.0.
- RASPA reports both absolute and excess loading, in several units, for example:

  Average loading absolute [molecules/unit cell]
  Average loading absolute [mol/kg framework]
  Average loading absolute [milligram/gram framework]
  Average loading absolute [cm^3 (STP)/gr framework]
  Average loading absolute [cm^3 (STP)/cm^3 framework]

  and the same set for excess.

- Experiments almost always report excess adsorption. Absolute is the actual number of molecules in the pores; excess subtracts the amount that would already be in the pore volume at the bulk-gas density.
- Excess adsorption needs HeliumVoidFraction from a separate helium Widom run (the Talu–Myers probe, typically helium at room temperature). Without it, excess is not meaningful. At MOF-relevant pressures (often several bar), absolute and excess already differ; at high pressure excess can even fall while absolute stays high.
- molecules/unit cell is per crystallographic unit cell, not per simulation supercell. With UnitCells 2 2 2 the supercell is 8 unit cells, but RASPA already reports per unit cell in that line. Still, do not mix mol/kg, mg/g, cm3 STP/g, and molecules/uc without converting.


3. Force field, charges, cutoff, and box size

The usual suspects for mismatch with literature are simulation length, system size, cutoff, shifted vs truncated+tail, electrostatics, the crystal structure, and unit conversion.

For polar adsorbates (CO2, N2, H2O, ...):

- ChargeMethod Ewald (not None)
- framework charges either from the CIF (UseChargesFromCIFFile yes, via _atom_site_charge) or from pseudo_atoms.def (UseChargesFromCIFFile no). Use the same choice as the published model.
- guest charges from pseudo_atoms.def
- check Movies/System_0/Framework_0_initial_P1.cif: the _atom_site_charge column is what was actually used. The output should also show a near-zero net framework charge.

A typical Van der Waals setup is CutOffVDW 12.0, with shifted yes / tailcorrections no (or truncated with tail corrections) in force_field_mixing_rules.def. The supercell must satisfy shortest box length > 2 × cutoff. A 1 1 1 MOF cell is often too small.

CIF atom labels are a frequent silent failure: O1, O2, C1 often do not match O, C in pseudo_atoms.def, so those atoms get no Lennard-Jones parameters and/or charge 0, and uptake comes out far too low. Use RemoveAtomNumberCodeFromLabel yes, or map the labels explicitly, and confirm every framework type appears in the force-field output.

Generic MOF force fields (UFF/DREIDING plus TraPPE) also often underpredict uptake at open metal sites, especially for CO2 and water. That is a model limitation, not a RASPA setting.


4. Fugacity and Rosenbluth weight

- If FugacityCoefficient is omitted (or 0), RASPA converts pressure to fugacity with Peng–Robinson using Tc, Pc, and the acentric factor from the molecule file. FugacityCoefficient 1.0 treats the input as fugacity, i.e. the ideal-gas approximation. At high pressure that changes the imposed chemical potential.
- IdealGasRosenbluthWeight must be 1.0 for rigid molecules. For flexible chains it must be precomputed (Widom in an empty box at the same temperature). A wrong value shifts the chemical-potential reference and can systematically lower the loading.

Check in the output that the computed fugacity coefficient, bulk density, and "amount of excess molecules" look sensible, and that the fluid is still a gas (not past the vapor pressure).


5. Structure

Use a solvent-free, correctly bonded CIF. Block inaccessible pockets with BlockPockets / BlockPocketsFileName if the structure has them (sodalite cages, etc.). Leftover solvent, missing hydrogens, or a non-P1 cell read with the wrong symmetry will all change the pore volume and the loading.


6. Quick health checks in the output

- Energy drift should be ~1e-5 or better (running energy vs fully recomputed energy). Larger drift means the energies, and therefore the loading, are not trustworthy.
- Swap acceptance should be non-zero; at high loading prefer CFCMC.
- Measured chemical potential / fugacity from Widom should be close to the imposed value once the system is equilibrated.
- NumberOfCycles is the production run. A RASPA cycle is max(20, N) move attempts, so "50000 cycles" is not 50000 insertions. Trust the block averages and error bars, not the cycle count.


7. Validating against experiment

Match the quantity first, then the physical model.

What to plot against the isotherm
- Use excess loading, in the same units as the paper (usually mmol/g, mg/g, or cm3 STP/g). RASPA already prints these if HeliumVoidFraction is set.
- Compare at the same temperature. Pressure: at low P, fugacity ≈ pressure; at several bar you should convert experimental pressure to fugacity (Peng–Robinson, or the EOS the paper used) or let RASPA do that conversion and plot vs the experimental pressure consistently.
- The cleanest check is the Henry regime (low P): slope of loading vs pressure, or a separate Widom Henry coefficient. If Henry is already too low, the high-pressure isotherm will be too low as well.
- Do not compare above the adsorbate vapor pressure. Experimental gas-phase isotherms stop there; beyond it you are in a different regime and excess can even become negative.
- Convert units using the simulated framework mass, not a formula-unit mass that omits extra-framework species or includes solvent.

What must match the experimental sample
- Same crystal structure: solvent-free CIF, correct topology, no collapsed or desolvated phase unless that is what was measured. If the experimental BET area or pore volume is much larger/smaller than the perfect-crystal helium void volume, the loadings will not match even with a perfect force field.
- Activation, leftover solvent, defects, missing linkers, and extra-framework ions in the real sample are the most common reasons experiment sits above or below a perfect-crystal GCMC run.
- Open metal sites, humidity, and polar adsorbates: a generic force field without extra site–guest parameters will usually underpredict. That is expected, not a sign that SwapProbability is wrong.
- Framework flexibility: a rigid GCMC run can miss extra uptake from breathing or linker rotation. If the experimental isotherm is temperature- or pressure-stepped, you may need a flexible model.
- Mixture or impure feed gas in experiment (especially water) will not match a dry single-component GCMC run.

A practical validation sequence
1. Reproduce a published simulation for the same MOF + adsorbate + force field (not the experiment yet). If you cannot match the paper, the input is wrong.
2. Check Henry coefficient / low-pressure slope against experiment.
3. Check the full excess isotherm, with HeliumVoidFraction from a helium Widom run, in the experimental units.
4. Check a second observable if available (isosteric heat, preferential binding site). Loading alone can be right for the wrong reasons.
5. Only then decide whether remaining disagreement is force field, flexibility, or the real sample.

A practical order of checks: (1) pressure in Pa and excess vs absolute / matching units, (2) swap move actually accepting, (3) charges + Ewald for polar gases and CIF labels mapped to the force field, (4) box ≥ 2×cutoff, (5) force field and CIF matching the paper you compare to, (6) helium void fraction and experimental excess. If you post simulation.input, the force-field files, the molecule definition, and the loading table from the output, the failure mode is usually visible in a few lines.