GROMACS

What is GROMACS?

GROMACS (GROningen MAchine for Chemical Simulations) is one of the most widely used molecular dynamics engines in computational biology. Originally developed at the University of Groningen in 1991, it has grown into a community-maintained open-source project known for exceptional performance on both CPUs and GPUs.

Molecular dynamics (MD) simulates how atoms move over time by integrating Newton's equations of motion. For proteins, this reveals conformational flexibility, loop motions, domain rearrangements, and stability under different conditions. GROMACS handles the full pipeline: system setup, solvation, energy minimization, equilibration, production dynamics, and trajectory analysis.

How does GROMACS work?

A GROMACS simulation progresses through several stages, each building on the previous one.

Energy minimization

The input structure often contains steric clashes or suboptimal bond geometries from crystallography or homology modeling. Steepest-descent minimization relaxes these high-energy contacts before dynamics begins. The algorithm iteratively adjusts atomic positions to reduce the total potential energy until convergence or a step limit is reached.

NVT equilibration

With clashes resolved, the system is heated to the target temperature. During NVT (constant Number, Volume, Temperature) equilibration, position restraints hold the protein backbone in place while the solvent and ions settle around it. The V-rescale thermostat couples the system to a heat bath, generating a canonical ensemble with correct kinetic energy fluctuations.

NPT equilibration

Once the temperature is stable, the system transitions to NPT (constant Number, Pressure, Temperature) equilibration. A Berendsen barostat adjusts the simulation box volume until the density converges to the target pressure. Position restraints remain active during this phase to prevent premature structural drift while the box dimensions adjust.

Production MD

Restraints are released and the system evolves freely. The Parrinello-Rahman barostat replaces Berendsen for production runs because it produces a true isothermal-isobaric ensemble, which is important for accurate thermodynamic sampling. Coordinates are saved at regular intervals to form the trajectory.

Trajectory analysis

After production, MDTraj computes structural metrics across all saved frames: backbone RMSD, per-residue RMSF, radius of gyration, secondary structure assignment (DSSP), solvent-accessible surface area, hydrogen bond networks, Ramachandran angles, and inter-residue contact maps.

How to use GROMACS online

ProteinIQ runs GROMACS on cloud infrastructure with automated system preparation, solvation, and multi-stage equilibration. No command-line setup or force field file management is needed.

Input

The structure should contain standard amino acid residues. Non-standard residues, modified amino acids, or HETATMs not recognized by the selected force field will cause topology generation to fail. Pre-process structures with PDBFixer if needed.

Simulation settings

Environment settings

Advanced settings

Output

The simulation produces:

• Final structure: PDB file of the last frame, centered and processed
• Trajectory: XTC file containing all saved frames (protein-only if water removal is enabled)
• Analysis CSVs: Structural metrics computed across the trajectory (see interpreting results below)

Force fields

Each force field represents a different parameterization philosophy. The choice affects secondary structure propensities, sidechain dynamics, and solvation behavior.

For general protein dynamics, AMBER99SB-ILDN is the most common choice in recent literature. GROMOS96 54a7 offers faster simulations at the cost of losing explicit non-polar hydrogen detail.

Water models

All simulations use explicit solvation in a dodecahedral box with 1.0 nm padding around the solute.

The water model choice has little effect on folded protein structure but can influence unfolded-state conformational sampling and solvation dynamics.

Interpreting results

RMSD

Root mean square deviation of backbone atoms from the initial structure, measured in nanometers. A plateau in the RMSD trace indicates the protein has reached a stable conformation. Continuously rising RMSD suggests the simulation may need more time, or that the protein is undergoing a conformational transition.

RMSF

Root mean square fluctuation per residue captures local flexibility. High RMSF values correspond to mobile loops, termini, and disordered regions. Low RMSF indicates rigid core residues or those involved in stable secondary structure elements.

Radius of gyration

Measures overall compactness. A stable radius of gyration indicates the protein maintains its fold. A sharp increase suggests partial unfolding or domain separation.

Secondary structure (DSSP)

Per-residue secondary structure assignment across all frames. Tracks helix-to-coil transitions, beta-sheet stability, and transient structure formation. The summary reports percentage of helix, sheet, and coil content per frame.

Potential energy

Should decrease during minimization and stabilize during production. Large fluctuations during production may indicate simulation instability.

Additional analyses

• Hydrogen bonds: Persistent intramolecular H-bonds (frequency $\geq 0.1$) and per-frame H-bond counts
• Ramachandran angles: $\phi$/$\psi$ dihedral distributions per residue across all frames
• SASA: Per-residue solvent-accessible surface area over time (Shrake-Rupley algorithm)
• Contact map: $C_\alpha$-$C_\alpha$ distance-based contacts ($< 8$ angstrom cutoff) with occupancy frequencies
• B-factors: Crystallographic temperature factors estimated from RMSF ($B = \frac{8\pi^2}{3} \cdot \text{RMSF}^2$)

Limitations

• Only standard amino acid residues are supported. Non-standard residues, post-translational modifications, and ligands require manual parameterization that is not available through the web interface.
• Simulations run on CPUs. For the same wall-clock budget, GPU-accelerated engines like OpenMM can sample longer timescales.
• The maximum simulation duration is 200 ns, which may be insufficient for slow conformational transitions, protein folding, or large-scale domain motions that occur on microsecond timescales.
• Protein backbone is treated as fully flexible from the start of production MD. There is no option for partial restraints during production.