← Selected Work

Humidity-dependent CO₂ capture in ultraporous MOF-177

Complete research record · 2025 · GCMC/MD · MOF-177 · LAMMPS

Peer-reviewed findings paired with an inspectable LAMMPS walkthrough. Published values and archive-derived teaching examples are identified separately throughout.

PublishedRunnableArchive-derivedTutorial

Published article · 2025

Adhikari, B., Bhusal, A., Sun, Q., & Adhikari, K. Computational and Theoretical Chemistry, 1253.

View DOI ↗
Illustration of the MOF-177 framework with carbon dioxide molecules
Fig. 1 — MOF-177 under adsorption conditions. The porous framework is investigated using pressure, temperature, and humidity studies coupled to molecular analysis.

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.
Method boundary: GCMC and MD answer complementary questions; neither substitutes for the other.

Published findings · DOI 10.1016/j.comptc.2025.115419

What the article reports

Pressure · 308 K17.19 → 26.92

mmol/g from 10 to 100 bar. Higher pressure fills more of the available pore volume.

Temperature · 10 bar−30.1%

Relative uptake change from 298 to 318 K; increased thermal motion weakens adsorption.

Humidity · 308 K, 10 bar−12.5%

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

  • Publicationdoi:10.1016/j.comptc.2025.115419 ↗
    Peer-reviewed record available through the publisher DOI.
  • Inputin.lmp
    GCMC/MD loading protocol.
  • MoleculeCO2.mol · h2o.mol
    TraPPE CO₂ and SPC/E water templates.
  • DataMOF177.data
    Framework geometry and force-field data.
  • Analysisrdf.lmp
    RDF, MSD, and interaction-energy calculations.
  • Collectorcollect_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:

01Build the systemMOF-177, molecule templates, charges, Lennard-Jones parameters, and periodic cell.
02Load with GCMCInsert and delete CO₂ until the loading fluctuates around equilibrium at a chosen pressure and temperature.
03Relax with MDStart from the loaded snapshot and sample positions, energies, and mean-square displacement.
04InterpretTurn molecule counts into uptake, and trajectories into RDF, diffusion, and interaction-energy trends.

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 study

Pressure-308K

10–100 bar at 308 K. This isolates crowding: higher pressure inserts more CO₂ into the same pore network.

Temperature study

Temperature-10bar

298–318 K at 10 bar. This tests the balance between adsorption interactions and molecular motion.

Humidity study

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.

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 equal variable holds one value. An index variable 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.
  • labeljump SELF is how LAMMPS writes a loop. At the bottom, clear wipes the whole system from memory, the two next commands advance both lists in step, and jump restarts 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/3 fixes 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 4 shifts 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.16 and angle_coeff 1 100 180 define 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 group commands give names to sets of atoms, so later commands can say "the framework" or "the CO₂" instead of listing types.
  • neighbor 2.0 bin adds a 2 Å skin to the neighbor lists — routine bookkeeping that avoids rebuilding the lists every step.
  • The exclude line 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 dynamic because GCMC changes its membership constantly — a static group would keep the members it had when it was defined.
  • dynamic/dof yes handles 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 real the 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:

ArgumentValueMeaning
N100Attempt a GCMC cycle every 100 MD timesteps.
X30Each cycle makes 30 exchange attempts — insert a molecule at a random position and orientation, or delete a random existing one.
M0Zero Monte Carlo displacement moves. The rigid-body MD between cycles handles all the motion instead.
type0Unused placeholder — the inserted species comes from the mol CO2mol template instead.
seed65899Random-number seed. Keep it if you want a bitwise-repeatable test; change it for an independent statistical sample.
T308Reservoir temperature in kelvin, matching the thermostat.
μ−1.0A chemical potential that is ignored, because pressure is given below — the pressure route is easier to reason about.
displacement0.5Maximum MC move size — unused here since M = 0.
rigid myrigidNewly inserted molecules are handed to the rigid-body integrator, and group CO2 adds them to the CO₂ group.
fugacity_coeff, pressureper conditionTogether 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:

Gravimetric uptake
qCO2=1000NCO2Mframeworkmmol g−1
⟨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.
Eq. 1 — Gravimetric uptake. Avogadro’s constant cancels when molecule count and framework molar mass refer to the same simulation cell. Source: study analysis convention.

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

N(CO₂)0production beginsstepslate mean = 899.43transient loading
transient: do not average production fluctuations averaging window

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

StudyControlled conditionVariableWhat it answers
Pressure308 K10, 20, …, 100 barHow loading, crowding, and CO₂ mobility change as the reservoir is compressed.
Temperature10 bar298, 303, 308, 313, 318 KHow thermal motion weakens adsorption but can increase transport.
Humidity308 K, 10 barPre-adsorbed water contentHow 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 toff gymnastics are needed here — everything comes straight from the data file.
  • This stage integrates CO₂ with a plain nvt thermostat, 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 0 restarts 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":

PairAtomsWhat a peak means
1 5Zn ↔ C(CO₂)CO₂ sitting near the zinc-oxide corner clusters.
2 5, 3 5H, C(MOF) ↔ C(CO₂)CO₂ against the organic linkers.
4 5O(MOF) ↔ C(CO₂)The pair the paper highlights: a peak near 4 Å shows the framework oxygens anchoring CO₂.
5 5, 6 5, 6 6CO₂ ↔ 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 yes subtracts any drift of the group's centre of mass, so you measure diffusion rather than collective sliding.
  • c_msd1 is 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.0001 only prevents a divide-by-zero at step 0.
D=limt|r(t)r(0)|26t
Eq. 2 — Einstein relation. Self-diffusion follows from the long-time slope of the three-dimensional mean-square displacement. Source: standard statistical-mechanics relation.

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₂ toff mapping 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

ComponentStatusWhat is availableBoundary
ArticlePublishedDOI and peer-reviewed scientific record.Publisher record is authoritative for reported values.
Inputs + dataPublic archiveThree condition folders, LAMMPS inputs, templates, and framework data.Not end-to-end verified as a fresh reproduction.
EnvironmentSpecifiedRequired LAMMPS packages are named in this guide.Exact historical binary and package versions are not recorded.
AnalysisInspectablerdf.lmp and collection script are provided.Stored outputs contain documented inconsistencies.
VerificationAuditedArchive differences are stated in section 7.No unqualified “reproduce” claim is made.
RightsBoundedPublic 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
  1. Confirm each molecule template’s shifted atom types against the Masses, Pair Coeffs, and group definitions.
  2. Calculate water wt% after the replicate 2 2 1 command, not from the unreplicated cell.
  3. Write a fresh post-GCMC data file and check that water atom types and the intended number of water molecules are present.
  4. Visualise one dry and one wet snapshot in OVITO before analysing RDF or diffusion.
  5. 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.