gromacs

name

gromacs

description

Create, run, and analyze molecular dynamics simulations using GROMACS. Use this skill when the user asks for: (1) System setup and topology preparation, (2) Ligand parameterization and small molecule setup, (3) .mdp parameter file design, (4) Simulation workflows (EM → NVT → NPT → production), (5) HPC job scripts and performance tuning, (6) Trajectory analysis (energy, density, MSD, SASA, RDF), (7) Free energy calculations (umbrella sampling, FEP, PMF), (8) Troubleshooting simulation failures (LINCS, exploding systems, barostat instability, PLUMED errors), (9) Polymer-solvent MD simulations for COSMO-RS analysis.

GROMACS Molecular Dynamics Skill

Guide for reproducible, production-grade GROMACS MD workflows: system building → equilibration → production → analysis.

Design Thinking

Before writing commands, identify:

• Goal: Stability, sampling, kinetics, thermodynamics, free energy?
• Model: Force field, water model, constraints strategy
• Risk: LINCS failures, bad contacts, barostat instability, PBC artifacts
• Deliverable: Command sequence, job script, .mdp blocks, analysis recipe

Standard Workflow

1. System Preparation

Protein Systems

# Generate topology (example for protein)
gmx pdb2gmx -f input.pdb -o processed.gro -water tip3p -ff charmm27

# Define box and solvate
gmx editconf -f processed.gro -o boxed.gro -c -d 1.2 -bt dodecahedron

# Add WATER using gmx solvate (for SPC, TIP3P, TIP4P water models)
gmx solvate -cp boxed.gro -cs spc216.gro -o solvated.gro -p topol.top

# Add ions to neutralize
gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o system.gro -p topol.top -pname NA -nname CL -neutral

Note: Use gmx solvate for water, gmx insert-molecules for all other solvents (see below).

Ligand/Small Molecule Parameterization

# Option 1: From compound name (PubChem lookup)
python scripts/ligand_setup.py --name "hexane" --resname HEX

# Option 2: From SMILES string
python scripts/solvent_to_gmx.py --smiles "CCCCCC" --resname HEX --base hexane

# Option 3: Complex molecules from SMILES
python scripts/solvent_to_gmx.py \
  --smiles "O=C([O-])C1=C(O)C(/N=N/C2=CC=C(C)C=C2S(=O)([O-])=O)=C3C=CC=CC3=C1" \
  --resname DYE --base dye_molecule

# Outputs:
# - HEX.acpype/HEX_GMX.itp (topology)
# - HEX.acpype/HEX_GMX.gro (structure)
# - HEX.acpype/atomtypes.itp (force field parameters)

# Include in topology:
# #include "HEX.acpype/atomtypes.itp"
# #include "HEX.acpype/HEX_GMX.itp"

See references/ligand-parameterization.md for detailed ligand setup.

Solvation Guide

Use gmx solvate for WATER:

# Standard water solvation (SPC, TIP3P, TIP4P models)
gmx solvate -cp protein.gro -cs spc216.gro -o solvated.gro -p topol.top

# With custom box
gmx solvate -cp protein.gro -cs spc216.gro -box 8 8 8 -o solvated.gro -p topol.top

# Add water shell around protein
gmx solvate -cp protein.gro -cs spc216.gro -shell 1.0 -o solvated.gro -p topol.top

Use gmx insert-molecules for NON-WATER solvents:

# Add organic solvent (e.g., ethanol, DMSO, ionic liquids)
gmx insert-molecules -f protein.gro -ci ethanol.gro -nmol 500 -o solvated.gro

# Create pure solvent box
gmx insert-molecules -ci hexane.gro -nmol 100 -box 5 5 5 -o hexane_box.gro

# Mixed solvents (water + organic)
gmx solvate -cp protein.gro -cs spc216.gro -o water_solvated.gro -p topol.top
gmx insert-molecules -f water_solvated.gro -ci dmso.gro -nmol 50 -o mixed_solvent.gro

2. Energy Minimization

gmx grompp -f em.mdp -c system.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em

Validation: Check em.log for Fmax < emtol and potential energy < 0.

3. Equilibration (NVT → NPT)

# NVT (temperature coupling)
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt

# NPT (pressure coupling)
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt

Validation: Temperature stable at target ±5K, density ~1000 kg/m³ for aqueous.

4. Production MD

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md

Core Tools Reference

gmx grompp (Preprocessor)

Creates .tpr run input from topology + coordinates + parameters.

Key flags:

• -f .mdp parameters, -p topology, -c coordinates
• -t checkpoint for continuation (velocities, coupling variables)
• -r restraint reference coordinates (required for position restraints)
• -n custom index groups
• -po mdout.mdp shows actual interpreted parameters
• -pp processed.top shows expanded topology (debug macros/includes)

gmx mdrun (Engine)

Executes simulation from .tpr.

Key flags:

• -deffnm sets base name for all outputs
• -cpi state.cpt continues from checkpoint
• -maxh N wall-time limit (writes checkpoint before exit)
• -nb gpu -pme gpu -bonded gpu -update gpu GPU acceleration
• -ntmpi N -ntomp M MPI ranks × OpenMP threads

Analysis Tools

Tool	Purpose	Key Options
gmx energy	Energy/pressure/temperature extraction	Interactive selection
gmx density	Density profiles along axis	-center, -symm for bilayers
gmx msd	Mean-square displacement, diffusion	-mol for per-molecule, -beginfit/-endfit
gmx sasa	Solvent-accessible surface area	-surface, -output selections
gmx rdf	Radial distribution functions	-ref, -sel groups
gmx trjconv	Trajectory processing/PBC fixing	-pbc mol, -center, -fit

Solvation Tools

Tool	Use Case	Example
gmx solvate	Water only (SPC, TIP3P, TIP4P)	gmx solvate -cp protein.gro -cs spc216.gro -o solvated.gro -p topol.top
gmx insert-molecules	All non-water solvents (organic, ionic liquids, etc.)	gmx insert-molecules -f protein.gro -ci ethanol.gro -nmol 500 -o solvated.gro

Key difference: gmx solvate intelligently fills space and updates topology; gmx insert-molecules randomly inserts molecules.

Default Parameters

Unless specified otherwise, use:

• Thermostat: v-rescale (stochastic velocity rescaling, tau_t = 0.1-1.0 ps)
• Barostat: C-rescale (stochastic cell rescaling, tau_p = 2.0-5.0 ps)
• Cutoffs: 1.2 nm for vdW, PME for electrostatics
• Constraints: h-bonds with LINCS (or all-bonds if needed)
• Time step: 2 fs (with h-bond constraints), 1 fs (unconstrained)

Utility Scripts

• scripts/ligand_setup.py: Automated ligand parameterization from compound name (PubChem → RDKit → ACPYPE)
• scripts/solvent_to_gmx.py: Convert SMILES to GROMACS topology (SMILES → RDKit → ACPYPE)
• scripts/generate_swarm_job.sh: Generate SLURM job scripts for swarm cluster
• scripts/check_equilibration.py: Validate NVT/NPT equilibration quality
• scripts/analyze_density.py: Track density evolution during MD simulations (requires MDAnalysis)
• scripts/gmx_wrapper.sh: Auto-detect gmx_mpi vs gmx

Specialized Workflows

Polymer-Solvent MD for COSMO-RS

Automated end-to-end workflow for polymer oligomer simulations in organic solvents, designed for COSMO-RS conformational analysis.

Quick Start:

python workflows/polymer-cosmo/polymer_md_workflow.py \
  --polymer PVDF \
  --solvent NMP \
  --work-dir ./pvdf_nmp

Features:

• Auto-generates polymer oligomers with backbone atom validation (ideal: 10-12 atoms, max: 20)
• Built-in polymer database: PVDF, PEO, PMMA, PS
• Built-in solvent database: NMP, DMF, DMSO, THF, water, etc.
• Intelligent solvent insertion (timeout-based retry with automatic target reduction)
• Merged atomtypes with [ defaults ] section (prevents duplicate atomtype errors)
• Generates swarm cluster and workstation submission scripts
• Creates solvated boxes with 10,000-15,000 solvent molecules in 11×11×11 nm box

What it does:

1. ✅ Validates polymer + solvent SMILES (with user confirmation)
2. ✅ Parameterizes both using ACPYPE (via solvent_to_gmx.py)
3. ✅ Creates solvated simulation box
4. ✅ Tests energy minimization
5. ✅ Generates submission scripts for 5ns NPT + 10ns NVT

Documentation:

• workflows/polymer-cosmo/README.md: Complete user guide
• workflows/polymer-cosmo/WORKFLOW_SUMMARY.md: Technical details

Example usage:

# PVDF 4mer in NMP
python workflows/polymer-cosmo/polymer_md_workflow.py \
  --polymer PVDF \
  --solvent NMP \
  --oligomer-size 4 \
  --work-dir ./test_pvdf_nmp

# Custom polymer and solvent SMILES
python workflows/polymer-cosmo/polymer_md_workflow.py \
  --polymer "C(C(F)(F))([H])[H]" \
  --solvent "CN1CCCC1=O" \
  --work-dir ./custom_sim

Future modules (in development):

• Conformational sampling (Rg-SASA grid extraction from trajectory)
• DFT input generation (ORCA COSMO-RS files)

Reference Files

For detailed information, see:

• references/swarm-workflow.md: Workstation preparation + swarm cluster submission workflow
• references/ligand-parameterization.md: Ligand topology generation, ACPYPE, CGenFF, ATB
• references/mdp-templates.md: Complete .mdp examples for EM, NVT, NPT, production, and specialized runs
• references/free-energy.md: Umbrella sampling, steered MD, PMF calculations, FEP
• references/troubleshooting.md: Common errors (LINCS, exploding systems, NaN) and fixes
• references/mpi-configuration.md: gmx vs gmx_mpi setup, PLUMED issues

Note: Create cluster-specific workflow files (like swarm-workflow.md) for each HPC resource you use, as requirements vary significantly between clusters.

Quick Troubleshooting

Symptom	Likely Cause	Fix
LINCS warning	Bad contacts, large forces	Longer/stricter EM, reduce dt
System explodes	Overlapping atoms, bad topology	Check topology, visual inspection
NaN in energy	Unstable simulation	Reduce dt, check parameters
Density too low/high	Wrong pressure coupling	Check barostat settings, ref_p
Temperature drift	Thermostat misconfigured	Check tc-grps match index
Atom count mismatch	Topology molecule count ≠ .gro	Check if diff = atoms/molecule, fix [ molecules ]
PLUMED symbol error	PLUMED library not loaded	Set LD_LIBRARY_PATH or use gmx

See references/troubleshooting.md and references/mpi-configuration.md for detailed diagnostics.

Note: Examples use gmx for simplicity. On workstation, use gmx_mpi. On swarm cluster, use gmx (see swarm-workflow.md and mpi-configuration.md).