MD Simulation of a Protein-Protein Complex in GROMACS

To simulate a protein-protein complex in GROMACS, run the normal pdb2gmx, editconf, solvate, genion, minimisation, NVT, NPT and production pipeline on one PDB that holds both chains, with clean chain IDs. Give the whole complex at least 1.0 to 1.2 nm of solvent to the box edge, split the chains into index groups with make_ndx, repair periodic images before analysis, and judge the interface by inter-chain RMSD, contacts and hydrogen bonds.

A docking server or AlphaFold-Multimer gives you a static guess at how two proteins fit together. Molecular dynamics is how you ask whether that guess survives in water at 300 K. The core GROMACS workflow is the same one you would use for a single protein (our GROMACS molecular dynamics guide walks through it end to end), but a two-chain system adds five specific places where students lose days: chain handling in pdb2gmx, box size, index groups, periodic boundary artefacts, and choosing analyses that actually describe the interface. This guide, written by the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modelling), covers each one with the exact commands.

What changes when you simulate two chains instead of one?

Very little in the physics, and quite a lot in the bookkeeping. GROMACS treats each chain as its own molecule by default, so the topology, the index groups and the trajectory processing all need to know there are two. The table below sets the single-protein habits you already have next to what a complex needs.

How do you prepare the complex PDB before pdb2gmx?

Start from the model you actually want to test, usually the top-ranked or most-populated cluster from ClusPro, HDOCK or HADDOCK, or the best-ranked AlphaFold-Multimer model judged by its interface PAE. Then check three things in a text editor or PyMOL:

1. Chain IDs. The two partners must carry different chain IDs (for example A and B), or be separated by a TER record. Some docking outputs reuse one chain ID for both partners. pdb2gmx’s -chainsep option decides where one chain ends: by default (id_or_ter) a new chain starts when the ID changes or a TER record appears.
2. Heteroatoms. Delete waters, crystallisation additives and ligands you do not intend to parameterise. Ions or cofactors that sit at the interface need their own parameters; see our guide on preparing a protein for MD.
3. Missing atoms and residues. Docked crystal structures often have gaps. A missing loop at the interface will distort exactly the region you are studying, so model missing residues before you start.

How do you run pdb2gmx on a two-chain complex?

Run it once on the whole file. The force field and water model choices are the same as for a single protein:

gmx pdb2gmx -f complex_clean.pdb -o complex.gro -p topol.top -ignh -ff amber99sb-ildn -water tip3p

With the default -merge no, GROMACS writes one include file per chain, typically topol_Protein_chain_A.itp and topol_Protein_chain_B.itp, each with its own position-restraint file, and lists both chains under [ molecules ] in topol.top. The official gmx pdb2gmx documentation describes -merge as controlling how multiple chains are merged into one moleculetype, with the choices no, all and interactive.

When should you use -merge all? Bonds can only exist inside a moleculetype. If your complex has an inter-chain disulfide (antibody heavy and light chains are the classic example), separate moleculetypes cannot hold that S-S bond. Use -merge all (or -merge interactive to merge only the chains involved) together with -ss so pdb2gmx asks you about each candidate disulfide. If there is no covalent link between the partners, keep the default; separate chains make the topology easier to read and the analysis groups easier to build. For a deeper tour of the flags, see our pdb2gmx tutorial.

How big should the box be for a protein complex?

Big enough that the complex never sees its own periodic image within the cut-off. A complex is elongated, and it rotates during the run, so size the box around the whole assembly:

gmx editconf -f complex.gro -o box.gro -c -d 1.2 -bt dodecahedron

-d 1.2 places at least 1.2 nm between the complex and the box edge, so two periodic copies are at least 2.4 nm apart at the start. A rhombic dodecahedron needs roughly 71% of the volume of a cubic box with the same image distance, which cuts the number of water molecules you pay for, and that matters more for a large complex than for a small protein. Use 1.0 nm as a floor; with CHARMM36, which uses a 1.2 nm cut-off, 1.2 nm or more is the safer choice. Our editconf tutorial explains the box types in detail.

Solvation and ions are unchanged from the single-protein recipe:

gmx solvate -cp box.gro -cs spc216.gro -o solv.gro -p topol.top
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15

Choose the SOL group when genion asks which group to replace.

After the run, verify that the box held. gmx mindist -pi reports the shortest distance between the protein and its periodic images; select your Protein group:

gmx mindist -s md.tpr -f md.xtc -pi -od mindist_pi.xvg -tu ns

If that minimum drops below your cut-off at any point, the complex interacted with itself through the boundary and the box was too small.

Do minimisation and equilibration change for a complex?

No. Run energy minimisation, then NVT and NPT equilibration with position restraints on both chains (the -DPOSRES define switches on every posre_*.itp that pdb2gmx included), then unrestrained production. tc-grps = Protein Non-Protein still works because the default Protein group already contains both chains. The GROMACS tutorials by Justin Lemkul give tested .mdp files for each stage, and our NVT and NPT equilibration guide explains how to confirm temperature, pressure and density have settled before you start production.

The one complex-specific check: open npt.gro in VMD or PyMOL and confirm the two chains are still in contact. A docked pose with a bad clash can push the partners apart during minimisation, which tells you to revisit the docking model rather than run 100 ns on a broken starting point.

How do you make index groups for each chain?

Most interface analyses need chain A and chain B as separate groups. gmx make_ndx has a command for exactly this:

gmx make_ndx -f em.gro -o index.ndx
> splitch 1
> name 19 chainA
> name 20 chainB
> 19 & 4
> name 21 chainA_BB
> 20 & 4
> name 22 chainB_BB
> q

splitch 1 splits group 1 (Protein) into one group per chain, appended to the end of the list. The numbers 19 to 22 above are examples: read the actual numbers make_ndx prints on your system before naming them. Group 4 is Backbone in the default list; the & lines create backbone-only groups for each chain, which you will use for fitting. Our make_ndx tutorial covers the selection syntax if your system needs interface-residue groups too.

Why do the two chains look separated in the trajectory?

Almost always because of periodic boundary conditions, not because the complex fell apart. GROMACS keeps each molecule’s atoms inside the box, so chain A can be drawn in one periodic image and chain B in the neighbouring one. Your RMSD, distance and contact numbers will be nonsense until you fix this.

Try the simple route first, centring on the whole Protein group:

gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center -ur compact

Choose Protein for centring and System (or Protein) for output. If the chains still split in some frames, use the clustering route. The gmx trjconv documentation describes -pbc cluster as moving all atoms of a group so that they sit closest to the cluster’s centre of mass, with a caveat worth reading twice: “this will only give meaningful results if you in fact have a cluster.”

gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -pbc whole
gmx trjconv -s md.tpr -f md_whole.xtc -o first_cluster.gro -pbc cluster -dump 0
gmx grompp -f md.mdp -c first_cluster.gro -p topol.top -o md_cluster.tpr
gmx trjconv -s md_cluster.tpr -f md_whole.xtc -o md_nojump.xtc -pbc nojump

The first line makes broken molecules whole. The second clusters the complex in the first frame only (select Protein for clustering). The third builds a run input file whose reference coordinates are that clustered frame, and the fourth uses -pbc nojump so no atom jumps across the boundary relative to that reference. The complex then stays together for the rest of the trajectory, though it may drift out of the box, which is harmless for analysis. Finish by centring and fitting if you want a clean movie. Our trjconv tutorial explains each -pbc mode.

How do you tell whether the docked interface holds?

Use metrics that describe the relative position of the partners, not just the shape of each one. A classic two-chain complex shows what these metrics watch: bovine trypsin bound to its inhibitor BPTI (PDB 2PTC), where the chains are labelled E and I rather than A and B, and the inhibitor’s lysine 15 reaches into the protease’s binding pocket. Four analyses carry most of the evidence:

Inter-chain RMSD

gmx rms -s md.tpr -f md_center.xtc -n index.ndx -o rmsd_chainB_on_A.xvg -tu ns

Select chainA_BB for the least-squares fit and chainB_BB for the RMSD calculation. Because only chain A is superimposed, this number grows when chain B slides or rotates on the surface of chain A, which is the same idea as the ligand RMSD used to grade docking poses. In the CAPRI assessment framework (Méndez et al., Proteins, 2003), even an “acceptable” model must recover at least 10% of the native contacts and have a ligand RMSD under 10 Å or an interface RMSD under 4 Å, so drifts of several ångströms are meaningful. GROMACS reports nm; 1 nm is 10 Å. Also run the ordinary backbone RMSD of each chain fitted on itself, so you can tell “chain B moved” apart from “chain B unfolded”. Our RMSD and RMSF guide covers plotting and reading these curves.

Contacts and minimum distance

gmx mindist -s md.tpr -f md_center.xtc -n index.ndx -od mindist_AB.xvg -on numcont_AB.xvg -d 0.35 -tu ns

Select chainA, then chainB. -od writes the minimum distance between the chains over time and -on the number of atom pairs closer than -d (default 0.6 nm; 0.35 nm counts tighter, near-hydrogen-bond contacts). A minimum distance that jumps from around 0.2 to 0.3 nm to well over 1 nm and stays there means the partners have separated.

Inter-chain hydrogen bonds

gmx hbond -s md.tpr -f md_center.xtc -n index.ndx -r 'group "chainA"' -t 'group "chainB"' -num hbnum_AB.xvg -tu ns

This is the selection-based gmx hbond in current GROMACS releases. The gmx hbond documentation notes that the reference and target selections must be either identical or non-overlapping, which two chains always are. On releases older than 2023, or with gmx hbond-legacy, you pick the two groups at an interactive prompt instead, and there -r means a distance cut-off, not a selection.

Buried surface area

Run gmx sasa three times with -surface set to chainA, chainB and Protein, then compute buried surface area as SASA(A) + SASA(B) minus SASA(complex). In their survey of 75 protein-protein complexes, Lo Conte, Chothia and Janin (J. Mol. Biol., 1999) found that most “standard-size” interfaces bury 1600 ± 400 Ų in total. GROMACS reports nm², so that range is roughly 12 to 20 nm². A docked interface far below that range, or one that shrinks steadily during the run, deserves scepticism. For a residue-level view of the same interface on a single frame, run it through PDBePISA.

Tap to explore in 3D

Two chainsInterface residuesLysine in the pocket

How long should you run, and how many replicas?

There is no universal length. Many student projects run 100 ns per replica; what matters is whether the interface metrics above reach a plateau and whether independent replicas agree. A single run in which the complex stays bound is weak evidence, because one trajectory samples one path. Run at least three replicas from the same equilibrated structure with different starting velocities (set gen_vel = yes with a different gen_seed in each NVT run), and report mean and spread. Our replicate simulation guide shows how to set these up and present them.

If you want an energy estimate on top of the structural evidence, MM/PBSA on the stable part of the trajectory is the usual next step for students; treat its absolute values with caution and use it to compare complexes rather than to quote a binding affinity.

What errors come up most often, and how do you fix them?

• pdb2gmx joins both partners into one chain and reports long bonds. The two partners share a chain ID and there is no TER record between them, so GROMACS bonds the last residue of A to the first residue of B. Relabel one chain, or add a TER record, then rerun.
• “Atom … in residue … was not found in rtp entry”. Usually hydrogens named in a different convention by the docking server. Use -ignh so pdb2gmx rebuilds hydrogens. If a heavy atom is named in the message, the residue is incomplete; rebuild it before pdb2gmx.
• The inter-chain disulfide is missing from the topology. The chains were kept as separate moleculetypes. Rerun with -merge all (or interactive) and -ss.
• grompp: “The cut-off length is longer than half the shortest box vector”. The box is too small for the cut-off. Rebuild it with a larger -d in editconf, then re-solvate.
• grompp: number of coordinates does not match topology. The [ molecules ] section in topol.top no longer matches the .gro file, often after a manual edit or after deleting a chain. Make sure both protein chains, SOL, NA and CL appear in the same order and counts as in the coordinate file.
• Chains appear to fly apart in VMD. Check the PBC-corrected trajectory before concluding anything. If they stay apart after the clustering route and the inter-chain minimum distance stays above 1 nm, the complex really dissociated; check whether other replicas agree.

Our list of common GROMACS errors covers the general ones that are not specific to complexes.

Frequently asked questions

Can I simulate an AlphaFold-Multimer model directly in GROMACS?

Yes, once it is cleaned like any other model: confirm distinct chain IDs, remove low-confidence tails that would flop around the box if they are not part of your question, and check that the interface residues have reasonable confidence. The GROMACS steps are identical to those for a docked complex.

Should I keep the two chains as separate molecules or merge them?

Keep them separate (the default -merge no) unless they share a covalent bond such as an inter-chain disulfide. Separate moleculetypes make per-chain topology and analysis groups simpler, and they do not change the physics.

Does it matter which chain I fit on for the inter-chain RMSD?

Fit on the larger or more rigid partner, usually the receptor, and calculate on the smaller one. That mirrors how docking poses are graded and makes the curve easier to read. Report which chain you fitted on in your methods section.

How do I know the complex really dissociated and it is not a PBC artefact?

Process the trajectory with the clustering and nojump route, then plot the inter-chain minimum distance from gmx mindist. A true dissociation shows a sustained distance well above 1 nm in the corrected trajectory; a PBC artefact disappears after correction.

Which analyses belong in a thesis results section for a complex?

Per-chain backbone RMSD, inter-chain RMSD fitted on the receptor, RMSF per chain, inter-chain hydrogen bonds and contacts over time, buried surface area, and, if you ran it, MM/PBSA. Show replicas on the same plot so readers can judge agreement.

Protein-protein complex simulation sits in the molecular dynamics stage of our computational biology skills roadmap, right after docking and structure prediction.