mmol/g from 10 to 100 bar. Higher pressure fills more of the available pore volume.
Humidity-dependent CO₂ capture in ultraporous MOF-177
Peer-reviewed findings paired with an inspectable LAMMPS walkthrough. Published values and archive-derived teaching examples are identified separately throughout.
PublishedRunnableArchive-derivedTutorial
Adhikari, B., Bhusal, A., Sun, Q., & Adhikari, K. Computational and Theoretical Chemistry, 1253.
View DOI ↗On this page
The article and research question
MOF-177 has a very large internal surface area and can accommodate substantial CO₂. Real adsorption environments, however, are rarely perfectly dry. The study therefore asks a practical question: how do pressure, temperature, and pre-adsorbed water change the amount of CO₂ held by MOF-177, where the molecules prefer to sit, and how freely they move?
Answering that question requires more than a single uptake value. The work connects equilibrium loading with molecular structure, transport, and framework–adsorbate interactions so that the observed trends can be interpreted physically.
Study design: what was done
The study combines grand-canonical Monte Carlo with molecular dynamics. GCMC exchanges CO₂ molecules with a pressure-controlled reservoir until the loading reaches equilibrium; MD then follows the loaded population to measure spatial organisation, mobility, and interaction energy.
- Material
- A replicated, periodic MOF-177 framework with the framework atoms held fixed.
- Conditions
- Pressure: 10–100 bar at 308 K. Temperature: 298–318 K at 10 bar. Humidity: pre-adsorbed water at 308 K and 10 bar.
- Molecular models
- UFF-style framework interactions, MEPO-ML charges, TraPPE CO₂, and SPC/E water for the humidity study.
- Measurements
- Gravimetric uptake, radial distribution functions, mean-square displacement and diffusion, plus interaction energies.
Study map · published methodology
From experimental variable to physical evidence
Each condition changes one physical influence. GCMC establishes equilibrium loading; MD reads the organisation and motion of the loaded state. The map makes clear which measurement answers which question.
- Pressure
- At 308 K: how does a more compressed CO₂ reservoir change loading?
- Temperature
- At 10 bar: how does thermal motion change affinity and occupancy?
- Pre-adsorbed water
- At 308 K and 10 bar: how does a wet pore compete with CO₂?
- Uptake
- How much CO₂ is held at equilibrium, reported in mmol/g.
- RDF + interaction energy
- Where molecules sit and how strongly the framework–adsorbate system interacts.
- MSD + diffusion
- How mobility evolves after loading; the late-time MSD slope gives diffusion.
Published findings · DOI 10.1016/j.comptc.2025.115419
What the article reports
Relative uptake change from 298 to 318 K; increased thermal motion weakens adsorption.
At 15 wt% pre-adsorbed water relative to the dry case; water competes for sites and pore space.
All figures above are peer-reviewed article values. They are not derived from the archived inputs. Open the publisher DOI ↗
Accessible chart explorer · default: published pressure values
Compare the three condition studies
At 308 K, the published uptake rises from 17.19 mmol/g at 10 bar to 26.92 mmol/g at 100 bar.
Data status: published endpoints; the intervening line is a visual connection, not a fitted isotherm.
My contribution
As a co-author, I contributed across the computational workflow and its communication within the broader collaborative study.
- MethodologyHybrid GCMC/MD setup and condition sweeps.
- SoftwareLAMMPS inputs, analysis workflow, and reproducibility documentation.
- ValidationLoading, RDF, diffusion, and interaction-energy checks.
- WritingCo-authorship and this technical research guide.
File manifest
- Publication
doi:10.1016/j.comptc.2025.115419 ↗
Peer-reviewed record available through the publisher DOI. - Input
in.lmp
GCMC/MD loading protocol. - Molecule
CO2.mol·h2o.mol
TraPPE CO₂ and SPC/E water templates. - Data
MOF177.data
Framework geometry and force-field data. - Analysis
rdf.lmp
RDF, MSD, and interaction-energy calculations. - Collector
collect_md_data.py
Condition-level summary extraction.
1. From question to workflow
"How much CO₂ does the framework hold, and how freely does it move?" is really two different measurements. Adsorption is an equilibrium loading problem — counting molecules once the framework balances a gas reservoir — while diffusion is a time-dependent transport problem. No single simulation technique answers both, which is why every condition in this archive runs through the same four-stage pipeline:
The rest of this guide follows that pipeline in order: the model is built in section 3, loaded with GCMC in section 4, swept across pressure, temperature, and humidity in section 5, and turned into structure, transport, and energy measurements with MD in section 6.
2. Read the folder before you run anything
The archive is organised as three studies. Each one has an adsorption input (in.lmp), an MD/RDF input (rdf.lmp), a framework data file, a CO₂ molecule template, and a small Python collector for the MD summaries.
Pressure-308K
10–100 bar at 308 K. This isolates crowding: higher pressure inserts more CO₂ into the same pore network.
Temperature-10bar
298–318 K at 10 bar. This tests the balance between adsorption interactions and molecular motion.
Humidity-308K-10bar
Pre-adsorbed water at 308 K and 10 bar. Water is placed before CO₂ is sampled, so it competes for space and sites.
Core files to inspect
Run a study from inside its own folder so its relative paths resolve. Replace /path/to/your/project with the directory where you saved or cloned the files:
cd /path/to/your/project/Pressure-308K lmp -in in.lmp # adsorption/loading stage lmp -in rdf.lmp # MD, RDF, MSD, and energy stage
What you need installed
LAMMPS is one executable plus optional packages, and this workflow needs four of them: MOLECULE (bonded atoms and .mol templates), KSPACE (the PPPM long-range electrostatics solver), RIGID (the rigid/nvt/small integrator), and MC (the fix gcmc command). The easiest way to get a build with all four is conda:
conda create -n lammps -c conda-forge lammps conda activate lammps lmp -h | grep -A8 "Installed packages"
If a run stops immediately with Unknown fix style gcmc or Invalid kspace style, your executable was built without the package. That is a setup failure, not a parameter problem — nothing is wrong with the input file.
Do a five-minute smoke test first. The full input runs ten pressures × 500,000 steps — hours of computer time. Before committing, open in.lmp and temporarily change variable stepTimes equal 100 to equal 1. That shrinks each pressure to 5,000 steps: the numbers will be nowhere near converged, but in a few minutes you will see molecules being inserted, every output file being created, and the loop advancing to the next pressure — which proves your build and the files work before you spend real time.
3. Build the physical model
Open Pressure-308K/in.lmp next to this section — every block below appears in that file, in the same order, and nothing is skipped. The file has four jobs: define the sweep, set the global rules, load the framework, and define the CO₂ molecule that GCMC will insert.
The header: variables and the pressure loop
variable tfac equal 5.0/3.0
variable fug_coeff index 0.9512 0.9036 0.8571 ... 0.5347
variable pressBar index 10 20 30 40 50 60 70 80 90 100
variable temp equal 308
variable stepTimes equal 100
variable equilSteps equal 4000*${stepTimes} # 400,000 steps
variable prodSteps equal 1000*${stepTimes} # 100,000 steps
variable totalSteps equal ${equilSteps}+${prodSteps}
label next_pressure
...the whole simulation...
next fug_coeff
next pressBar
clear
jump SELF next_pressure
- An
equalvariable holds one value. Anindexvariable holds a list and starts at its first entry — so the first pass runs at 10 bar with fugacity coefficient 0.9512. The two lists are parallel: each pressure is paired with its own Peng-Robinson coefficient. label…jump SELFis how LAMMPS writes a loop. At the bottom,clearwipes the whole system from memory, the twonextcommands advance both lists in step, andjumprestarts the script at the label. When the lists run out, the jump is skipped and the script simply ends. One file — ten complete, independent simulations.tfac = 5/3fixes a subtle temperature bug you would otherwise never notice: a rigid linear molecule has five degrees of freedom (three translations, two rotations), but GCMC hands new molecules velocities assuming three. Without the correction, freshly inserted CO₂ starts slightly cold.
Global rules: units, style, boundaries
units real # Å, femtoseconds, kcal/mol, atmospheres atom_style full # atoms carry a charge plus bonds and angles boundary p p p # periodic in x, y, and z
units real decides what every number means for the rest of the run: distances in ångströms, time in femtoseconds (and a default timestep of 1 fs, which this input keeps), energies in kcal/mol, pressures in atmospheres. atom_style full is needed because both the framework and CO₂ carry partial charges and bonded topology. Periodic boundaries mean a molecule leaving one side of the box re-enters from the opposite side — the small cell behaves like an infinite crystal.
The framework: MOF-177
pair_style lj/cut/coul/long 12.0 pair_modify mix arithmetic bond_style harmonic angle_style harmonic read_data MOF177.data replicate 2 2 1 kspace_style pppm 1.0e-4
pair_style lj/cut/coul/long 12.0— every pair of atoms interacts through a Lennard-Jones term (shape and packing) plus a Coulomb term (partial charges). Both are computed directly out to 12 Å; the Coulomb tail beyond the cutoff is handed to the k-space solver below.pair_modify mix arithmetic— the data file defines Lennard-Jones parameters per atom type, but never for pairs of unlike types. Lorentz-Berthelot mixing builds those cross terms automatically: arithmetic mean for σ, geometric mean for ε.read_data MOF177.data— loads the 936-atom triclinic unit cell (64 Zn, 216 H, 456 C, 200 O) with its box vectors, masses, partial charges, and per-type Lennard-Jones coefficients. Skim the file's header once: it declares six atom types — four for the framework and two reserved for CO₂.replicate 2 2 1— duplicates the cell into a 2×2×1 supercell of 3,744 atoms. This is not cosmetic: with periodic boundaries the box must be at least twice the 12 Å cutoff in every direction, and a larger pore network gives better loading statistics.kspace_style pppm 1.0e-4— the particle-particle particle-mesh solver for long-range electrostatics, tuned so the estimated error is 10−4 of a reference force. Charged periodic systems cannot simply truncate Coulomb interactions; this is the standard fix.
The force-field choices come from the paper: UFF Lennard-Jones parameters for the framework atoms and MEPO-ML partial charges. You do not have to enter any of them by hand — they are already written into MOF177.data.
CO₂: a rigid TraPPE molecule
GCMC needs a template describing what to insert. CO2.mol is a three-site TraPPE model: carbon at the origin, an oxygen 1.16 Å away on each side, +0.70e on carbon and −0.35e on each oxygen. The molecule is neutral overall but carries the quadrupole — the charge pattern — that makes CO₂ stick to the framework.
# CO2.mol: local molecule types 1 and 2
Coords: C 0.00 0.00 0.00
O -1.16 0.00 0.00
O 1.16 0.00 0.00
Charges: C +0.70e; O −0.35e each
Angle: O–C–O = 180°
Three lines in in.lmp wire the template into the simulation:
molecule CO2mol CO2.mol toff 4 bond_coeff 1 1000 1.16 angle_coeff 1 100 180
toff 4shifts the template's local atom types 1 and 2 up to global types 5 and 6 — the two slots the data file reserved. Forget this offset and your CO₂ would be inserted as framework atoms, with the wrong mass, charge, and interactions.bond_coeff 1 1000 1.16andangle_coeff 1 100 180define stiff harmonic springs for the C=O bond (rest length 1.16 Å) and the 180° O–C–O angle. During this GCMC stage they never actually act — the rigid-body integrator holds each CO₂ perfectly rigid — but the MD analysis stage in section 6 integrates the molecules flexibly, and there these springs are what keep CO₂ linear.
Groups, bookkeeping, and one big speed trick
group MOF177 type 1 2 3 4 group CO2 type 5 6 neighbor 2.0 bin neigh_modify exclude group MOF177 MOF177
- The two
groupcommands give names to sets of atoms, so later commands can say "the framework" or "the CO₂" instead of listing types. neighbor 2.0 binadds a 2 Å skin to the neighbor lists — routine bookkeeping that avoids rebuilding the lists every step.- The
excludeline is the trick: this workflow treats the framework as rigid and frozen — no integrator ever moves a MOF atom. Since frozen atoms cannot exert forces on each other that matter, all framework–framework pair computations are skipped. In a 3,744-atom framework, that removes the vast majority of the pair interactions and makes the run dramatically faster.
Water: SPC/E, when humidity is studied
The humidity input creates water molecules before the GCMC stage. The SPC/E template has one oxygen (−0.8472e) and two hydrogen atoms (+0.4236e each), 1 Å O–H bonds, and a 109.47° H–O–H angle. This archive models pre-adsorbed water; it does not perform GCMC exchange moves for water.
Type offsets matter. A molecule file numbers its own types from 1. The toff option shifts those numbers into the types defined in the LAMMPS data file. Always make a four-column map—template type, shifted type, group membership, and pair coefficient—before running a multicomponent system.
4. Understand — then run — the hybrid GCMC/MD stage
Grand Canonical Monte Carlo samples a system at fixed temperature and chemical potential (specified here through pressure and fugacity). Its special power is that it can change the number of CO₂ molecules, which ordinary molecular dynamics cannot do — MD only moves the atoms it already has. This input interleaves the two: MD lets the adsorbed molecules wander and relax; GCMC periodically tries to insert new molecules from the gas reservoir or delete existing ones. The balance of accepted insertions and deletions is the adsorption equilibrium.
Counting molecules while they appear and disappear
compute temp_CO2 CO2 temp
compute_modify temp_CO2 dynamic/dof yes
variable CO2Car atom "type==5"
group CO2Car dynamic all var CO2Car
variable nCO2Mol equal count(CO2Car)
variable pressAtm equal ${pressBar}*0.986923
- Every CO₂ has exactly one carbon (type 5), so counting type-5 atoms counts molecules. The group must be declared
dynamicbecause GCMC changes its membership constantly — a static group would keep the members it had when it was defined. dynamic/dof yeshandles a bookkeeping subtlety: temperature is kinetic energy per degree of freedom, and the number of degrees of freedom changes with every insertion and deletion. Without this flag, the printed CO₂ temperature drifts away from reality as loading grows.- The last line converts the human-friendly bar value to atmospheres, because in
units realthe GCMC fix expects atmospheres.
The two fixes that do the work
fix myrigid CO2 rigid/nvt/small molecule mol CO2mol temp 308 308 100
fix mygcmc CO2 gcmc 100 30 0 0 65899 308 -1.0 0.5 mol CO2mol &
tfac_insert ${tfac} group CO2 rigid myrigid &
fugacity_coeff ${fug_coeff} pressure ${pressAtm}
The first fix is the integrator: it advances every CO₂ as one rigid body — translating and tumbling, never bending — with a Nosé-Hoover thermostat at 308 K (the trailing 100 is the thermostat's damping time, in femtoseconds). The second fix is the molecule exchanger, and its long argument list is worth decoding once:
| Argument | Value | Meaning |
|---|---|---|
N | 100 | Attempt a GCMC cycle every 100 MD timesteps. |
X | 30 | Each cycle makes 30 exchange attempts — insert a molecule at a random position and orientation, or delete a random existing one. |
M | 0 | Zero Monte Carlo displacement moves. The rigid-body MD between cycles handles all the motion instead. |
| type | 0 | Unused placeholder — the inserted species comes from the mol CO2mol template instead. |
| seed | 65899 | Random-number seed. Keep it if you want a bitwise-repeatable test; change it for an independent statistical sample. |
| T | 308 | Reservoir temperature in kelvin, matching the thermostat. |
| μ | −1.0 | A chemical potential that is ignored, because pressure is given below — the pressure route is easier to reason about. |
| displacement | 0.5 | Maximum MC move size — unused here since M = 0. |
rigid myrigid | — | Newly inserted molecules are handed to the rigid-body integrator, and group CO2 adds them to the CO₂ group. |
fugacity_coeff, pressure | per condition | Together they define the gas reservoir. The fugacity coefficient (from the Peng-Robinson equation of state) corrects for the gas being non-ideal — essential at 100 bar, where CO₂ is far from ideal. |
So the rhythm of the simulation is: 100 femtoseconds of rigid-body dynamics → 30 insertion/deletion attempts, accepted or rejected with grand-canonical probabilities → repeat, 5,000 times per pressure.
Watching it run
thermo_style custom step temp c_temp_CO2 press pe ke density v_nCO2Mol
thermo 1000
dump dmp all custom 1000 Outputs/dumps/dump_${pressBar}.lammpstrj ...
fix loading all ave/time 100 20 2000 v_nCO2Mol &
file Outputs/numbers/no_of_mol_${pressBar}.profile
run ${totalSteps}
write_data Outputs/dumps/afterGCMC_${pressBar}.data
Start it from inside the study folder:
cd Pressure-308K lmp -in in.lmp # serial mpirun -np 8 lmp -in in.lmp # or parallel, if your build supports MPI
A thermo line prints every 1,000 steps, and the column to watch is the last one, v_nCO2Mol: it starts at 0 and climbs quickly as GCMC floods the empty pores. Three output streams are written as it runs:
Outputs/numbers/no_of_mol_*.profile— the molecule count, sampled every 100 steps, averaged in blocks of 20, written every 2,000 steps (ave/time 100 20 2000). The averaging smooths the insert/delete noise into a readable curve.Outputs/dumps/dump_*.lammpstrj— a trajectory snapshot every 1,000 steps. Drag it into OVITO and you can literally watch the pores fill.Outputs/dumps/afterGCMC_*.data— the final loaded configuration, written once at the end. This is the input for the MD analysis stage in section 6.
Reading the result: from a count to an isotherm point
The values below are taken directly from Pressure-308K/no_of_mol_10.profile. The final 100,000-step window runs from 400,000 to 500,000 steps and contains 51 recorded averages.
# TimeStep v_nCO2Mol 2000 50.60 ← early loading into empty pores 16000 397.55 36000 612.30 ← still approaching the late-time regime ... 400000 891.98 ← start of the final 100,000-step window 492000 900.28 500000 903.23
The curve climbs steeply, then fluctuates around a late-time mean. Averaging only the final 100,000 steps gives ⟨N⟩ = 899.43 CO₂ molecules (51 profile records; standard deviation 16.36). Converting that profile-derived mean to standard units takes one formula and the mass of the replicated framework:
- ⟨NCO₂⟩
- Mean number of adsorbed CO₂ molecules in the production window.
- Mframework
- Molar mass of the replicated empty framework, in g mol−1.
- 1000
- Conversion from mol g−1 to mmol g−1.
One unit cell weighs 64 × 65.409 (Zn) + 216 × 1.008 (H) + 456 × 12.011 (C) + 200 × 15.999 (O) ≈ 13,081 g/mol, and the 2×2×1 supercell is four of those: 52,323 g/mol. Thus the profile average gives 1000 × 899.43 / 52,323 ≈ 17.19 mmol/g, matching the paper’s reported 10-bar value.
Cross-check complete. The final-100,000-step mean from no_of_mol_10.profile converts to the paper’s reported 17.19 mmol/g at 10 bar. The article remains the authoritative scientific record for the reported result.
Two habits worth building immediately: read the plateau from a late-time window, never from the first flat-looking point; and if the count is still trending upward at the end of equilibration, extend the run rather than averaging a slope.
Equilibration reader · archive-derived teaching example
Find the stable window before calculating uptake
A count that is still drifting upward is not an equilibrium result. The useful average begins only after the trace fluctuates around a stable mean.
Fig. 2 · loading count versus simulation step
Teaching boundary: this trace translates the archived 10-bar output into a reading guide. Its 400,000–500,000-step average is ⟨N⟩ = 899.43, which converts to the paper’s 17.19 mmol/g result.
5. Change one variable at a time
| Study | Controlled condition | Variable | What it answers |
|---|---|---|---|
| Pressure | 308 K | 10, 20, …, 100 bar | How loading, crowding, and CO₂ mobility change as the reservoir is compressed. |
| Temperature | 10 bar | 298, 303, 308, 313, 318 K | How thermal motion weakens adsorption but can increase transport. |
| Humidity | 308 K, 10 bar | Pre-adsorbed water content | How water occupies sites/free volume before CO₂ is sampled. |
Each pressure and temperature condition has its own Peng-Robinson fugacity coefficient in the input. Do not reuse a coefficient from a neighbouring condition: fugacity is the correction that connects the non-ideal gas reservoir to the GCMC acceptance probability.
For the humidity experiment described in the paper, water is randomly introduced before CO₂, at 3–15 wt% water. The important scientific distinction is that this is a preloading protocol: it measures CO₂ uptake in a wet pore network, rather than equilibrating CO₂ and water against two separate grand-canonical reservoirs.
6. Turn a loaded snapshot into physical measurements
The GCMC stage answered "how much?" — this stage answers "sitting where, moving how fast, bound how strongly?" The companion script rdf.lmp loops over the ten afterGCMC_*.data snapshots, relaxes each one, and runs plain molecular dynamics — no more insertions or deletions, just watching a fixed population move. Run it the same way, after in.lmp has finished:
lmp -in rdf.lmp
Restart, relax, reset
read_data Outputs/dumps/afterGCMC_${pressBar}.data
group MOF177 type 1 2 3 4
group co2 type 5 6
group c type 5 # one carbon per molecule
fix nvt_co2 co2 nvt temp 308 308 100.0
thermo 1000
run 50000
reset_timestep 0
- The snapshot already contains the loaded CO₂ with correct global types, so no
toffgymnastics are needed here — everything comes straight from the data file. - This stage integrates CO₂ with a plain
nvtthermostat, not the rigid-body fix. That means the harmonic bond and angle springs from section 3 are now doing real work: the 1,000 kcal/mol/Ų bond and 100 kcal/mol/rad² angle constants are stiff enough to keep each molecule effectively linear. The framework remains frozen, exactly as before. - 50,000 steps of relaxation let the population settle after the switch of integrator, and then
reset_timestep 0restarts the clock — important because the diffusion analysis below divides displacement by time, and that time must start when the measurement starts.
Adsorption amount
Uptake was already computed from the GCMC plateau in section 4; the snapshot gives the same answer, since it freezes that final population. One sanity check worth knowing: replication increases both the molecule count and the framework mass by the same factor, so a correctly normalised uptake does not depend on the supercell size. If yours does, one of the two numbers is from the wrong cell.
RDF: where molecules prefer to sit
compute myRDF all rdf 100 1 5 2 5 3 5 4 5 5 5 6 5 6 6
fix rdf_out all ave/time 500 1 500 c_myRDF[*] ave running &
file Outputs/md/rdf_${pressBar}.txt mode vector overwrite
A radial distribution function g(r) asks: relative to random chance, how likely am I to find an atom of type B at distance r from an atom of type A? The compute uses 100 distance bins, and the numbers after it are read in pairs — "type A, type B":
| Pair | Atoms | What a peak means |
|---|---|---|
1 5 | Zn ↔ C(CO₂) | CO₂ sitting near the zinc-oxide corner clusters. |
2 5, 3 5 | H, C(MOF) ↔ C(CO₂) | CO₂ against the organic linkers. |
4 5 | O(MOF) ↔ C(CO₂) | The pair the paper highlights: a peak near 4 Å shows the framework oxygens anchoring CO₂. |
5 5, 6 5, 6 6 | CO₂ ↔ CO₂ | How the adsorbed gas packs against itself — these grow with pressure as the pores crowd. |
ave running means the histogram keeps improving for the whole run, and the file is rewritten with the latest average. In rdf_*.txt, each block lists the bin's distance, then alternating g(r) and coordination-number columns for each pair, in the order above. Read it like this: a sharp peak marks a preferred separation; g(r) settling toward 1 at long range means the structure has faded into uniform, bulk-like behaviour.
MSD and diffusion
compute msd1 c msd com yes
variable dcoeff4 equal c_msd1[4]/(step*dt+0.0001)/6
fix msd_out all ave/time 100 1 100 c_msd1[1] c_msd1[2] c_msd1[3] c_msd1[4] &
file Outputs/md/msd_${pressBar}.txt
fix dcoeff_out all ave/time 25000 10 250000 v_dcoeff1 v_dcoeff2 v_dcoeff3 v_dcoeff4 &
ave running file Outputs/md/dcoeff_${pressBar}.txt overwrite
- The mean-square displacement is tracked for group
c— the carbons, one per molecule, sitting at each molecule's centre.com yessubtracts any drift of the group's centre of mass, so you measure diffusion rather than collective sliding. c_msd1is a four-component vector: ⟨Δx²⟩, ⟨Δy²⟩, ⟨Δz²⟩, and their sum. The Einstein relation turns the total into a diffusion coefficient by dividing by 6t (or 2t per axis, useful for spotting anisotropic transport along particular pore directions). The+0.0001only prevents a divide-by-zero at step 0.
Units matter here: in units real the MSD is in Ų and time in fs, so the variable's value is in Ų/fs — multiply by 10−5 to get m²/s. The archived 10-bar file ends with v_dcoeff4 = 0.0040, i.e. D ≈ 4.0×10−8 m²/s as the mechanical reading. The averaging window (25000 10 250000) deliberately samples only late times, because the early MSD is curved — molecules rattling inside a cage look fast before they have actually gone anywhere. If you fit D yourself, fit only the late, linear region.
Interaction energy and collected data
compute engkspace MOF177 group/group co2 kspace yes
fix energy_out all ave/time 500 1 500 c_engkspace[*] c_engkspace &
ave running file Outputs/md/energyKspace_${pressBar}.txt overwrite
compute group/group isolates the interaction energy between two groups — here, everything the framework exerts on the CO₂, with kspace yes including the long-range electrostatic part. In the output file the final column is that energy in kcal/mol; more negative means more strongly bound. The humidity study's script adds the water–framework term the same way. One caution before comparing conditions: a 100-bar system holds far more molecules than a 10-bar one, so always divide by ⟨NCO₂⟩ and compare energy per adsorbed molecule.
When the final pressure finishes, the input calls collect_md_data.py, which sweeps the per-condition files into combined_energy_data.txt and combined_dcoeff_data.txt — one row per pressure, ready to plot an isotherm-style trend.
7. Archive audit: where the files disagree with the paper
Use the paper as the authoritative methodology. The supplied archive is valuable for learning the file structure, but it should not be treated as a validated, rerunnable reproduction dataset without a review.
- The humidity input’s water/CO₂
toffmapping conflicts with the displayed pair-coefficient ordering and group labels. - The humidity input lists five water counts (28–144), while the paper describes 3–15 wt% water. Recalculate water weight percent from the replicated framework mass before using those labels.
- The stored humidity and temperature output snapshots are duplicates of the pressure outputs and contain no water atom types, so they are not suitable examples of the wet-system results.
- The archived RDF scripts use 50,000 relaxation steps before a 500,000-step run, whereas the paper describes 500,000 equilibration plus 500,000 production steps.
This is not a reason to discard the archive. It is exactly the kind of audit a reproducible workflow should make possible: inspect type maps, compare counts to the method, and validate a small case before a long production campaign.
Reproducibility matrix · last audit 2026-07-13
What is public, inspectable, and still limited
| Component | Status | What is available | Boundary |
|---|---|---|---|
| Article | Published | DOI and peer-reviewed scientific record. | Publisher record is authoritative for reported values. |
| Inputs + data | Public archive | Three condition folders, LAMMPS inputs, templates, and framework data. | Not end-to-end verified as a fresh reproduction. |
| Environment | Specified | Required LAMMPS packages are named in this guide. | Exact historical binary and package versions are not recorded. |
| Analysis | Inspectable | rdf.lmp and collection script are provided. | Stored outputs contain documented inconsistencies. |
| Verification | Audited | Archive differences are stated in section 7. | No unqualified “reproduce” claim is made. |
| Rights | Bounded | Public methods and inputs only. | Publisher PDF, private collaboration records, and unpublished data are not redistributed. |
What to verify before running a corrected or new humidity case
- Confirm each molecule template’s shifted atom types against the
Masses,Pair Coeffs, and group definitions. - Calculate water wt% after the
replicate 2 2 1command, not from the unreplicated cell. - Write a fresh post-GCMC data file and check that water atom types and the intended number of water molecules are present.
- Visualise one dry and one wet snapshot in OVITO before analysing RDF or diffusion.
- Keep an immutable copy of the original input and record every corrected parameter in a run log.
8. Before you trust a result
- Does the CO₂ count plateau over the chosen averaging window?
- Do your bar, atmosphere, temperature, and LAMMPS-unit conversions agree with the input?
- Does the snapshot contain the species and atom types you think it contains?
- Is the MOF intentionally fixed or intentionally flexible, with the integration commands matching that choice?
- Are RDF, MSD, and energy normalised in a way that permits comparison across loadings?
- Can a short test reproduce expected qualitative trends before you spend cluster time on the full sweep?
Context glossary
The minimum vocabulary for this record
GCMC
Grand-canonical Monte Carlo inserts and deletes molecules to sample a system controlled by temperature and chemical potential or its pressure/fugacity representation.
Adsorption
Equilibrium accumulation of molecules in a material; here it is reported as CO₂ uptake per framework mass.
RDF
The radial distribution function compares the probability of a pair separation with a random reference, revealing preferred molecular neighbourhoods.
MSD
Mean-square displacement tracks average molecular motion. Its late-time linear slope can yield self-diffusion.
Interaction energy
The energy associated with interactions between selected groups, useful for interpreting relative binding only with consistent normalization.
MOF-177
An ultraporous metal–organic framework. Here it is represented as a replicated periodic framework for adsorption and transport analysis.
Research horizon · not completed work
Questions this record makes possible
- QuestionHow would a verified water-loading protocol change the humidity trend when water and CO₂ are both checked against independently regenerated snapshots?
- HypothesisNormalising framework–adsorbate energy per adsorbed CO₂ may separate stronger local binding from simple high-loading effects.
- In progressImprove the evidence boundary between peer-reviewed values and archive-derived teaching calculations across the public research record.
The goal is not to memorise every LAMMPS command. It is to keep a clear chain from the physical question to the molecular model, to the sampled ensemble, to the quantity plotted in a paper.