Molecular Dynamics Simulation with GROMACS

# GROMACS workflow
# 1. Prepare structure
gmx pdb2gmx -f protein.pdb -o structure.gro -water tip3p

# 2. Define simulation box
gmx editconf -f structure.gro -o box.gro -bt cubic -d 1.0

# 3. Solvate
gmx solvate -cp box.gro -cs spc216.gro -o solvated.gro

# 4. Add ions (neutralize)
gmx grompp -f ions.mdp -c solvated.gro -o ions.tpr
gmx genion -s ions.tpr -o system.gro -neutral

# 5. Energy minimization
gmx grompp -f em.mdp -c system.gro -o em.tpr
gmx mdrun -v -deffnm em

# 6. Production MD
gmx grompp -f md.mdp -c em.gro -o md.tpr
gmx mdrun -v -deffnm md -nt 8  # 8 CPU cores

FORCE_FIELDS = {
    'AMBER': 'Proteins, nucleic acids (biomolecules)',
    'CHARMM': 'Proteins, lipids, carbohydrates',
    'OPLS': 'Organic molecules, proteins',
    'GROMOS': 'Aqueous solutions, proteins',
}

# Force field parameters: bonds, angles, dihedrals, non-bonded
def calculate_forces(atom_positions, force_field):
    """
    F = -∇V where V = V_bonded + V_non-bonded

    V_bonded = V_bonds + V_angles + V_dihedrals
    V_non-bonded = V_electrostatic + V_van_der_Waals
    """
    energy = force_field.calculate_potential(atom_positions)
    forces = -gradient(energy)
    return forces

import MDAnalysis as mda

# Load trajectory
u = mda.Universe('system.gro', 'md.xtc')

# RMSD (structural stability)
from MDAnalysis.analysis import rms
R = rms.RMSD(u, select='backbone')
R.run()

# RMSF (flexibility per residue)
from MDAnalysis.analysis.rms import RMSF
rmsf = RMSF(u.select_atoms('backbone')).run()

# Hydrogen bonds
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
hbonds = HydrogenBondAnalysis(u).run()