Protein Molecular Dynamics Simulation

1. Download & Inspect the PDB File

Get your structure from the RCSB Protein Data Bank. Before doing anything else, open the raw .pdb file in a text editor and check for:

• MISSING entries — residues or atoms absent from the crystal structure. Terminal regions are common culprits. Any missing internal residues must be modeled in with external software (Modeller, Swiss-Model) before proceeding — pdb2gmx cannot handle gaps.
• Alternate conformations — pick one and remove the others.
• Non-standard residues or ligands — pdb2gmx only knows residues defined in the force field .rtp files. Anything else needs separate parameterization.

2. Remove Crystal Waters

Crystal waters are artifacts of the crystallization process. You will be adding a proper explicit solvent box yourself. Exception: if a water molecule is tightly bound in the active site and functionally important, keep it.

3. Run pdb2gmx

This step adds missing hydrogen atoms, assigns protonation states, and generates three output files:

4. Optional Flags

• -ter — interactively set N- and C-terminal charge states
• -inter — manually set protonation for Glu, Asp, Lys, Arg, His and disulfide bonds
• -ignh — ignore H atoms in PDB (useful for NMR structures)

Define the Simulation Box

• -c — centers the protein in the box
• -d 1.2 — places the protein at least 1.2 nm from any box edge
• -bt cubic — cubic box type. Rhombic dodecahedron is ~71% the volume and is more efficient for globular proteins

Solvate the Box

Uses spc216.gro, a generic 3-point solvent configuration compatible with SPC, SPC/E, and TIP3P water models. The topol.top file is updated automatically with the correct number of SOL molecules added.

Generate the Run Input File

grompp (GROMACS pre-processor) assembles the coordinate file, topology, and simulation parameters into a binary .tpr file that contains the complete atomic description of the system.

Add Counter-ions with genion

When prompted, choose group 13 "SOL" — you want to replace water molecules, not protein atoms.

• -pname NA / -nname CL — positive and negative ion names (always elemental symbol in all-caps)
• -neutral — adds the minimum number of ions to reach net zero charge
• -conc 0.15 — optionally add physiological ionic concentration (0.15 M) on top of neutralization

Assemble the Input

Run Energy Minimization

The -v flag makes mdrun verbose (prints progress). Steepest descent iteratively moves atoms downhill on the potential energy surface until the maximum force drops below the threshold.

Success Criteria

Check Potential Energy

Why NVT First?

After energy minimization, the solvent is geometrically relaxed but not thermally equilibrated. Jumping directly to a full production run would cause the system to collapse or explode. NVT uses position restraints on protein heavy atoms so the solvent equilibrates around the protein without structural changes.

Key .mdp Parameters

• gen_vel = yes — random velocity generation from Maxwell-Boltzmann distribution at target T
• tcoupl = V-rescale — stochastic velocity rescaling thermostat (Bussi et al.)
• pcoupl = no — pressure coupling is OFF at this stage
• define = -DPOSRES — activates position restraints from posre.itp

Monitor Temperature Convergence

Why NPT After NVT?

The NVT step fixed the temperature but not the density. The NPT ensemble (constant Number, Pressure, Temperature) most closely resembles real experimental conditions. The barostat adjusts the box volume until the pressure equilibrates.

Key .mdp Changes from NVT

• continuation = yes — continuing from NVT phase
• gen_vel = no — velocities read from NVT checkpoint
• pcoupl = C-rescale — C-rescale barostat (Bernetti & Bussi) for pressure control
• ref_p = 1.0 — target pressure 1 bar

Monitor Pressure & Density

Prepare and Launch

The checkpoint file from NPT carries the fully equilibrated coordinates, velocities, and barostat state. Position restraints are absent from the production .mdp file — the protein moves freely.

GPU Acceleration

GROMACS offloads PME (particle mesh Ewald) calculations and nonbonded interactions to GPU. As of GROMACS 2025, a Titan Xp GPU can run this lysozyme system at ~196 ns/day. Bonded forces remain on CPU cores.

Minimum GPU requirements: CUDA SDK + compute capability ≥ 2.0.

Key Production Parameters

Step 1 — Correct for Periodicity (PBC)

Proteins diffuse across periodic boundaries and can appear "broken" or jumping across the box. trjconv re-images the trajectory so the protein is intact and centered throughout.

RMSD — Structural Stability

Measures average displacement of backbone atoms from the reference. For lysozyme: 0.09 ± 0.01 nm — indicating a stable, compact fold. Note: RMSD is not a convergence criterion on its own.

Radius of Gyration (Rg) — Compactness

A stably folded protein maintains a near-constant Rg. For lysozyme: 1.409 ± 0.008 nm — confirming no unfolding or elongation over 10 ns.

Secondary Structure — DSSP

Invokes the DSSP algorithm to assign α-helices, β-strands, 3₁₀-helices, turns, and bends to each residue at each frame. The per-residue time series reveals which secondary structural elements are stable vs. transiently fluctuating.

Hydrogen Bonds

Criteria: donor-acceptor distance ≤ 0.35 nm AND donor-acceptor-H angle ≤ 30°. Lysozyme maintains ~50 backbone H-bonds, ~20 sidechain H-bonds, and ~280 protein-water H-bonds throughout the simulation.