CO2 hydrate nucleation (2026) - GROMACS REPPTIS =============================================== Replica Exchange Partial Path Transition Interface Sampling (REPPTIS) of CO2 gas hydrate nucleation at 270 K and 100 bar with GROMACS. The order parameter is a smoothed form of the Mutually Coordinated Guest (MCG) count of Barnes et al. [1]. Paper: in preparation. Files ----- retis.toml the PyRETIS 4 input of the study (task "repptis"); there is no PyRETIS 3 (retis.rst) counterpart. The current PyRETIS 4 refuses it as it stands: see "Running it with the current PyRETIS 4" below. orderp_soft.py order parameter MCG3Soft (needs numpy, scipy, mdtraj) gromacs_input/ GROMACS system: conf.gro configuration; also the topology of the order parameter ([orderparameter] topology) basinA_270K.gro a second configuration of the same system (same atoms, other coordinates, a 4.88 nm box); not read by retis.toml topol.top topology grompp.mdp the mdp template the engine reads forcefield/ TraPPE CO2 (reoptimised for aqueous solubility) and TIP4P/Ice water pyretis.mdp, files written by an earlier run, kept topol.tpr for reference; PyRETIS writes its own in the run directory System ------ 512 CO2 and 2944 TIP4P/Ice water molecules, 14336 atoms (conf.gro: a cubic box of 4.78 nm). md-vv with dt = 2 fs, Nose-Hoover at 270 K (tau-t = 2.0 ps), C-rescale at 100 bar (tau-p = 10.0 ps), PME for electrostatics and for van der Waals at r = 1.4 nm, no dispersion correction. One path step is subcycles = 1000 MD steps (2 ps); the 25 interfaces run from 17.2 to 250.0 in MCG3Soft, with zero_left = 9.0. Order parameter --------------- MCG3Soft is based on the MCG order parameter of Barnes, Beckham, Wu and Sum [1], introduced for quantifying clathrate hydrate nucleation and growth. MCG counts guest molecules whose neighbours are separated by shared water cages; MCG3 is the variant requiring at least three qualifying bonds per guest. MCG3Soft keeps that definition and replaces its yes/no criteria with a one-sided cubic smoothstep, so that the count is continuous in the atomic coordinates. The switch is exactly 1 wherever the hard criterion accepts and softens only the rejection: widening the gates can add near-misses but never discounts a true member. The widths shipped are the class defaults. The counting criteria are those of [1]; the smoothing is the authors' addition. Requirements ------------ * PyRETIS 4 (the repptis task, TOML input) * GROMACS with "gmx" on the PATH, built with PME and md-vv * the Python packages numpy, scipy and mdtraj Status / caveats ---------------- The input starts from initial paths ([initial-path] method = "load", load_folder = "load"). The load/ folder is NOT included: provide your own initial paths in load/, or change [initial-path] (see the PyRETIS user guide on initial paths) before running. These are the inputs of the study, written for a development version of PyRETIS 4. This is a production study of a large system: reproducing its statistics needs HPC resources. Edit [runner] for your machine first. It ships with three workers and one mdrun command per worker: workers = 3 wmdrun = [ "gmx mdrun -ntomp 8 -ntmpi 1", "gmx mdrun -ntomp 8 -ntmpi 1", "gmx mdrun -ntomp 8 -ntmpi 1", ] Running it with the current PyRETIS 4 ------------------------------------- As they stand, the current PyRETIS 4 refuses these inputs when it sets up the run, and its message names what to set: * gromacs_input/grompp.mdp sets tcoupl = Nose-Hoover. A path-sampling run refuses this thermostat: its friction variables carry a memory from one time step to the next, which starts again with every new piece of trajectory, and a deterministic thermostat run on reversed velocities does not retrace the path it ran forward. Set tcoupl = v-rescale nsttcouple = 1 (stochastic velocity rescaling, CSVR). This changes the thermostat the study used. * grompp.mdp sets no gen-temp. GROMACS draws the velocities of a run at gen-temp, and at 300 K when it is absent. Add gen-temp = 270 With both changes the set-up passes, and PyRETIS warns that pcoupl = C-rescale carries a memory between coupling steps that starts again with every piece of trajectory; for exact path sampling it advises pcoupl = no, at the box of an equilibrated constant-pressure run. The run then needs its initial paths in load/ (see "Status / caveats"). How to run ---------- PyRETIS 4: pyretis run -i retis.toml -p pyretis analyse -i retis.toml References ---------- [1] B. C. Barnes, G. T. Beckham, D. T. Wu and A. K. Sum, "Two-component order parameter for quantifying clathrate hydrate nucleation and growth", J. Chem. Phys. 140, 164506 (2014). https://doi.org/10.1063/1.4871898