Standalone enzyme–PHA trajectory preparation

The preparation workflow now uses Bash and GROMACS, with two templates in src/md_simulation_scripts/trajectory_preparation/. Copy them into a simulation’s analysis/ folder, edit the plain KEY=value instructions, then run bash process_trajectory.sh. No Python package or preparation notebook is required.

• process_trajectory.sh

• instrcution.txt

The complete setup, settings, atom groups, error handling, VMD inspection and contact-analysis instructions are maintained in the README below.

Trajectory preparation needs only Bash and GROMACS. The research analysis and enzyme-contact notebooks use Python analysis dependencies.

Parameterisation provenance

The GK13/ANC45 × P3HO4/P3HB4 production simulations analysed here used the CHARMM/GROMACS route: CHARMM36m + CGenFF (via CHARMM-GUI Ligand Reader & Modeler, from iPHASimulator’s R-configured SDF/PDB) + CHARMM TIP3P + SOD/CLA in GROMACS. No GAFF2 charges enter them.

src/md_simulation_scripts/
├── README.md
├── trajectory_preparation/
│   ├── prepare_md_files.ipynb
│   ├── templates/                    # tested combine + protein-only SH/TXT files
│   ├── process_trajectory.sh         # older manual SOLU/index-based template
│   └── instrcution.txt
├── enzyme_contacts/
│   ├── enzyme_contacts.ipynb
│   ├── run_enzyme_contacts.py
│   └── GK13_P3HO_4.yaml
└── enzyme_pha_analysis/
    └── enzyme_pha_analysis.ipynb

Automatic preparation notebook

Open trajectory_preparation/prepare_md_files.ipynb. Set SIMULATION_DIR to your raw simulation folder and Run All:

1. Detect and combine production EDR parts in the raw folder.

2. Detect and combine production XTC parts in the raw folder.

3. Centre on Protein only with -pbc mol -ur compact -center, retain System, and write the full processed XTC and a -skip 100 preview in analysis/.

The notebook defaults to ANC55_P3HO4 and its step8_production_1us continuations. Edit the filename settings for other systems. It installs the tested SH/TXT templates from trajectory_preparation/templates/ only when local files are missing; existing local settings must agree with the notebook settings. Combined files remain in the raw folder. Non-empty existing outputs are reused if their inputs are not newer; this does not revalidate their integrity. New continuation parts require fresh output names. The notebook stops on failure and records incomplete stages and GROMACS logs in analysis/. Set DRY_RUN=True to inspect the plan without writes or GROMACS execution.

This workflow uses standard Python, Bash and GROMACS; it does not require MDAnalysis. It does not modify teaching notebooks or research analysis notebooks.

The manual instructions below describe the older SOLU/index-based template at trajectory_preparation/process_trajectory.sh. For the current protein-only workflow, use the notebook or its templates/ SH/TXT files.

1. Prepare a new simulation’s analysis folder

Keep the combined XTC in the simulation folder. Copy the two templates and matching production TPR, viewing GRO and index into a new analysis/ folder. For example, after setting your repository and simulation paths:

REPO=/path/to/iPHASimulator_v2
SIM=/path/to/GK13_P3HO_4_gromacs
mkdir "$SIM/analysis"
cp "$REPO/src/md_simulation_scripts/trajectory_preparation/process_trajectory.sh" "$SIM/analysis/"
cp "$REPO/src/md_simulation_scripts/trajectory_preparation/instrcution.txt" "$SIM/analysis/"
cp "$SIM/step7_production.tpr" "$SIM/analysis/"
cp "$SIM/step6.2_npt.gro" "$SIM/analysis/"
cp "$SIM/index.ndx" "$SIM/analysis/"
cd "$SIM/analysis"

Use filenames from your system. The production TPR must describe the same atoms in the same order as the combined XTC; the GRO is for viewing, and must also match that order. Matching filenames or atom counts alone do not prove this. If analysis/ already exists, inspect its contents and preserve previous files before copying templates into it. These examples are setup instructions, not an instruction to replace your existing settings.

simulation/
├── production_combined_1us.xtc       # raw combined trajectory stays here
└── analysis/
    ├── process_trajectory.sh
    ├── instrcution.txt
    ├── step7_production.tpr
    ├── step6.2_npt.gro
    └── index.ndx

2. Edit instrcution.txt

This spelling is retained for compatibility. It is plain KEY=value data, not YAML. Do not use quotes, spaces around =, or inline comments; paths containing spaces are otherwise supported. Settings are read as data, never sourced as code. Paths resolve from the script/instruction folder, including when the script is invoked from another working directory.

INPUT=../production_combined_1us.xtc
TPR=step7_production.tpr
INDEX=index.ndx
CENTER_GROUP=SOLU
OUTPUT_GROUP=SYSTEM
STRIDE=100
PROCESSED=processed.xtc
PREVIEW=processed_every100.xtc
GMX=gmx

Change the input filenames and groups for a new system. GMX is an executable name or path (for example gmx_mpi), not a command with extra arguments. Output settings must be distinct .xtc filenames inside analysis/.

In GK13–P3HO₄, the inspected index contains:

Group	Atoms	Use
SOLU	3,800	Enzyme (3,701 atoms) + LIG (99 atoms); default centre
SOLV	133,560	Solvent and ions
SYSTEM	137,360	Entire simulation; default output

Use actual names from your index, not guessed group numbers. Retain the full system in OUTPUT_GROUP, so the same TPR/index also match the processed XTC in step 2. A subset needs its own matching topology/index and is outside this simple two-pass template. There is no automatic choice of centring group or fitting. For enzyme-only centring, explicitly create an enzyme group in your copied index first; Protein and protein Cα groups are not present in this example’s index. The previously inspected protein has 245 Cα atoms, but this script does not create selections from Python expressions.

3. Check and run

From the analysis folder:

bash process_trajectory.sh --dry-run
bash process_trajectory.sh

The dry run checks settings, executable availability, input file readability and output conflicts, then prints the two commands and selections. It executes no GROMACS command and writes no files. It cannot verify trajectory integrity, topology compatibility or index-group contents; check those for your system.

The normal run performs, in this order:

1. Read the combined XTC from the parent folder; use the production TPR and index with -pbc mol -ur compact -center. Select CENTER_GROUP, then OUTPUT_GROUP. Write every processed frame to processed.xtc, with no frame reduction or fitting in this step.

2. Read processed.xtc and write every 100th frame to processed_every100.xtc using -skip 100 (or the configured STRIDE). This keeps zero-based frames 0, 100, 200, …; it is a frame stride, not a time interval.

The script stops on errors and refuses existing output files, including symlinks. Rename earlier outputs or choose new PROCESSED/PREVIEW names to rerun. A failed run can leave a partial file: inspect it before deciding what to keep. Run only one preparation job per analysis folder at a time. Raw inputs are not modified. The full processed XTC may be roughly as large as the raw file; allow space for both full and reduced outputs.

4. Analyse the full trajectory; view the reduced one

Use processed.xtc for analysis. In MDAnalysis you can sample without creating another stored trajectory:

import MDAnalysis as mda
from MDAnalysis.lib.distances import distance_array

u = mda.Universe("step7_production.tpr", "processed.xtc")
enzyme = u.select_atoms("protein and not name H*")
pha = u.select_atoms("resname LIG and not name H*")
for ts in u.trajectory[::100]:
    distances = distance_array(enzyme.positions, pha.positions, box=ts.dimensions)
    # Consume the sampled distances here (MDAnalysis length units are Angstrom).

This is a simple sampling/PBC illustration; the contact workflow below provides more complete heavy-atom selection, summaries and plots. Enzyme–PHA distances must account for periodic boundaries using each frame’s box, even after centring. Centring does not bind a separated ligand to the enzyme or guarantee that the displayed pair shares the nearest periodic image. Do not apply rotational fitting and then calculate periodic distances using an unrotated box.

In VMD, load the matching step6.2_npt.gro, then add processed_every100.xtc to the same molecule. The GRO’s initial coordinates need not be the first processed frame; start viewing the added XTC frames (or remove the initial GRO frame). Display protein or resname LIG for the pair, or resname LIG for the PHA. Scrub through the beginning, middle and end, checking molecular continuity, periodic images and the chosen centre. Visual inspection is a manual step; running this script does not perform it.

Do not infer time coverage from 1us in a filename. Inspect the actual times. Do not concatenate the combined XTC with continuation files it already contains.

For Python analysis, continue with the enzyme-contact guide. The complete script-folder README also covers the moved entry points.