Messages

The warning is not that UFF/DREIDING/OPLS are incomplete. RASPA2 never matched your framework atoms to the types in the force field, so CO2 has no van der Waals interaction with the MOF. The results are not usable until this is fixed. Do not silence those pairs with "none".

RASPA2 takes the force-field type from _atom_site_label, not from _atom_site_type_symbol. Your CIF is:

Cu1  Cu     ...
O1   O_MOF  ...
C1   C_MOF  ...

so the types RASPA actually uses are Cu1, O1, C1, C2, C3, H1. That is exactly what the warning prints (O_CO2-Cu1, O_CO2-O1, ...). The mixing-rule types Cu, C_MOF, O_MOF, H_MOF are never applied. Unknown CIF labels are added as new pseudo-atoms with no Lennard-Jones parameters (and, if UseChargesFromCIFFile is no, charge 0). Check the pseudo-atom table in the output: you will see Cu1, O1, C1 listed separately from Cu, O_MOF, C_MOF.

You do not need a unique charge on every atom of every MOF. For a typed UFF/DREIDING-style model, give each chemical type one epsilon, sigma, and charge in pseudo_atoms.def / force_field_mixing_rules.def, then make the CIF labels match those type names. UseChargesFromCIFFile no is the right choice if conversion to IL-modified structures is wiping _atom_site_charge.

Two practical ways to make the names match, depending on how you generate CIFs.

1. Keep C_MOF / O_MOF / H_MOF / Cu in the force field (good if you want MOF atoms distinct from CO2, water, and the IL). Rewrite the CIF labels to those names. Labels do not have to be unique:

loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
Cu     Cu  0.2853 0.2853 0.0000
O_MOF  O   0.3166 0.2431 0.9478
C_MOF  C   0.2968 0.2032 0.9313
H_MOF  H   0.3802 0.2280 0.8802

Keep O_CO2 and C_CO2 only for the guest.

2. Better for many CoRE-style CIFs that already use Cu1, O1, C1, C2. Rename the framework types in the force field to Cu, O, C, H and put this in simulation.input:

RemoveAtomNumberCodeFromLabel yes

That maps Cu1 → Cu, O1 → O, C1/C2/C3 → C, H1 → H. Leave the guests as O_CO2 / C_CO2 so they do not collide with framework O and C. Do not use RemoveAtomNumberCodeFromLabel with types named O_MOF: O1 becomes O, which still will not match O_MOF.

After the names match, the warning should disappear for CO2–framework pairs. Framework–framework missing VDW (Cu–O, C–H, ...) is normal for a rigid framework and can be left undefined, or listed as none if you want a clean log.

Charges for high-throughput: one charge per type in pseudo_atoms.def is enough; you do not assign a charge per MOF. UseChargesFromCIFFile no. Make sure every framework type you actually use is in pseudo_atoms.def, otherwise RASPA adds it with charge 0. The CIF charges you showed (Cu 1.248, O −0.624, ...) are not the same as the pseudo_atoms.def values (Cu 0.836, O_MOF −0.834, C_MOF 0.417, H_MOF 0.417). Pick one consistent set. H_MOF = 0.417 looking identical to C_MOF is worth checking; that is not a typical DREIDING hydrogen charge.

Also fix the mixing-rules header. RASPA2 expects:

# general rule for shifted vs truncated
shifted no
# general rule tailcorrections
no

not the word "truncated" on its own.

Once types match, confirm in the output that (1) the warning is gone for O_CO2/C_CO2 vs Cu/C/O/H, (2) the framework net charge is what you expect, and (3) Movies/System_0/Framework_0_initial_P1.cif has the charges you intended. Mixing UFF/DREIDING (framework), OPLS (IL), and TraPPE-style CO2 is then a force-field choice: Lorentz–Berthelot will run, but Cu-BTC open-metal-site CO2 uptake is often still low with generic Cu parameters. That is a separate issue from this warning.

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.

This is a windows problem plus Qt problem. The windows port has now been rewritten in WinUI3 and DirectX12.

These are two different "grids." Energy interpolation grids only speed up host–guest energy evaluation. Adsorbate density plots come from a separate histogram, ComputeDensityGrid, which is sampled over many production configurations. A restart file is one snapshot; it is a good starting point for a new production run, not a density profile by itself.

Yes — restart from the JSON restart and run production with density grids on. You do not need energy/interpolation grids for that.

A density cube is an ensemble average: every SampleDensityGridEvery production cycles, adsorbate positions are binned onto a 3D grid and written as Gaussian cube files under density_grids/. One final configuration (or one restart) only shows where the molecules happened to sit at that moment.

Energy interpolation grids (UseInterpolationGrids in force_field.json) are optional and unrelated. They precompute framework–adsorbate energies. Turning them on later does not create density plots, and you can compute density grids without them.

Do not use RestartFromBinaryFile. That restores the original run's full state and ignores new simulation.json settings, so you cannot switch on density grids that way.

Use the JSON restart instead:

Copy the restart out of output/ so the next run does not overwrite it:
cp output/restart_<T>_<P>.s0.json .
Point the system at it with "RestartFileName", set "CreateNumberOfMolecules": 0, and enable density grids.
Because the configuration is already equilibrated, you can use few (or zero) initialization cycles, then enough production cycles for a converged histogram.
{
  "SimulationType": "MonteCarlo",
  "NumberOfInitializationCycles": 0,
  "NumberOfProductionCycles": 100000,
  "PrintEvery": 5000,
  "Systems": [
    {
      "Type": "Framework",
      "Name": "YourFramework",
      "NumberOfUnitCells": [2, 2, 2],
      "ExternalTemperature": 298.0,
      "ExternalPressure": 1.0e5,
      "RestartFileName": "restart_298_0.s0",
      "ComputeDensityGrid": true,
      "SampleDensityGridEvery": 10,
      "WriteDensityGridEvery": 5000,
      "DensityGridSize": [128, 128, 128]
    }
  ],
  "Components": [
    {
      "Name": "CO2",
      "TranslationProbability": 0.5,
      "RotationProbability": 0.5,
      "ReinsertionProbability": 0.5,
      "SwapProbability": 1.0,
      "CreateNumberOfMolecules": 0
    }
  ]
}
Output goes to density_grids/ as .cube files (view in VMD, iRASPA, etc.). Optional extras: "DensityGridPseudoAtomsList" for site-resolved maps (e.g. C_co2 vs O_co2), and "DensityGridBinning": "Equitable" for smoother histograms.

No — not from the final configuration alone. RASPA3 has no "make a density cube from this restart/PDB" step. If you already wrote a PDB movie (OutputPDBMovie) with many frames, you can histogram that trajectory yourself (for example VMD volmap). That is a post-processing substitute, not a built-in RASPA3 feature, and a single last frame is still not enough.

Practical recommendation: copy the JSON restart, enable ComputeDensityGrid, and run a production-only continuation. Add interpolation grids only if you want that continuation to run faster.

[you]

They are input, you set these. These probabilities are relative to each other, and they do not influence the accuracy of the computation, just the efficiency.

There is a 'help' in iRASPA ('help' menu). In principle, this is about understanding the three panes in the top-left: (I) projects, (II) scenes, (III) movie frames. If you select the second tab, you can add structures to the scene.
'File' -> 'New' -> 'Objects' -> 'Ellipsoid' and then edit the properties.

This works in raspa3 now.

I do not see the charges of the framework. How do you set these? In the input the charges for the framework appear to be zero.

First in RASPA do a Monte Carlo simulation with some particles and create the 3D density grid.
Page 172 manual:
ComputeDensityProfile3DVTKGrid yes
WriteDensityProfile3DVTKGridEvery 10000
DensityProfile3DVTKGridPoints 150 150 150

Than it creates files which can you import in iRASPA.

In RASPA3, it will write out a cube-file, which you can also import in iRASPA.

Use a restart-file:
Reads the positions, velocities, and force from the directory 'RestartInitial'. Any creation of molecules in the 'simulation.input' file will be in addition and after this first read from file. This is useful to loadinitial positions of cations for example, and after that create adsorbates. The restart file is written at 'PrintEvery' intervals.

Easiest is to first create a system with the desired number of molecules. Then run it to get a restart-file. Move the directory from 'Restart' to 'RestartInitial' and edit the positions. Then rerun with 'RestartFile yes' (and set 'CreateNumberOfMolecules 0' since they now come from the restart-file).

Different between what and what? What is the expected value and you're getting a wrong value? More info is needed.

Under "Appearances" you have "Volumetric data": click "Show grid data" and set rendering method to "isosurface".

No, you cannot get the energy per molecule, only the total energy of the component. Note however, that all molecules are equivalent, so you could just divide by the amount of molecules.

The pressure of the liquid-phase contains a large error (because of the strong molecular interactions). Usually people take the gas-phase value after the simulation is equilibrated.

I copied your input, and see no problems.
Have you tried to run it with zero cycles? Does it finish? is the output precisely (pseudo-atoms, force field etc) exactly what you put in?
If that works, run it with 10 cycles and see if it works.

One thing that could go wrong is that your input files are not ascci-files, but utf8 or windows-files that contain different newlines.
With utilities like 'dos2unix' and 'unix2dos' you convert between windows and linux/Mac.

You're swapping molecules in and out, and hence use the grand-canonical ensemble. In this ensemble, you fix the chemical potential (via the pressure/fugacity) and the amount of molecules fluctuates. So you're not able to control how many of where the molecules are. But that is the nature of adsorption.


Btw. you're missing a rotation move (although you could say this is also achieved via insertion of molecules, this is much less efficient).