GROMACS Protein Cofactor Simulation: Heme, NAD+ and FAD

To run a GROMACS protein cofactor simulation, first check whether your force field already defines the cofactor as a residue. The CHARMM36 port for GROMACS ships heme (HEME), NAD+ (NAD), NADH (NADH and NAI) and FAD (FAD and FADR). If it does, match the atom names to the .rtp and let pdb2gmx build it; if not, parameterize with CGenFF or GAFF.

This guide is written by the StemSkills Lab team, who have spent more than 10 years on structural bioinformatics, drug design and molecular modeling projects, many of them on enzymes that refuse to give up their cofactor. It is a spoke of our GROMACS molecular dynamics learning path and sits at step 6 (MD setup) of the computational biology skills roadmap.

Why does pdb2gmx stop when my protein has a cofactor?

gmx pdb2gmx builds a topology only for residues that exist in the force field’s residue topology database (the .rtp files inside the .ff folder). A heme, NAD or FAD in your PDB file is a HETATM residue, and if its name is not in any loaded .rtp, the run ends with Residue ‘XXX’ not found in residue topology database. The GROMACS error documentation explains that this happens when the database “does not contain the residue at all, or because the name is different.”

That second cause is the one most students miss. The cofactor may already be in the force field under a different name or with different atom names. Deleting it, which many first tutorials do for drug-like ligands, is the wrong default here: a dehydrogenase without NAD+ or a cytochrome without heme is not the enzyme you docked into.

Cofactor-bound structures are common. A search of the RCSB Protein Data Bank on 11 October 2026 returned 6,485 entries containing heme (chemical component HEM), 3,343 with FAD, 2,221 with NAD+ (NAD) and 476 with NADH (NAI). If you work on enzymes, you will meet this problem.

How do I find which cofactor is in my PDB file?

List every HETATM residue name. Columns 18 to 20 of a PDB record hold the residue name, so this works even when columns run together:

grep '^HETATM' 1A6M.pdb | cut -c18-20 | sort | uniq -c

For sperm whale oxymyoglobin (PDB 1A6M, 1.0 Å resolution) you will see HEM, the bound oxygen OXY, a sulfate ion SO4 and the waters HOH. Ions like that sulfate usually come from the crystallization buffer and are normally removed. Look each code up on its RCSB ligand page to confirm the chemical state. The codes matter: in the PDB, NAD is the oxidized form (NAD+), NAI is the reduced form (NADH), FAD is oxidized flavin adenine dinucleotide and FDA is the dihydro (reduced) form.

Then decide what stays. Keep the cofactor if the protein needs it to hold its fold or if your question involves the active site. Small ligands bound to the cofactor, like the O2 in 1A6M, need their own decision: OXY has no entry in standard protein force fields, so you either parameterize it or remove it and simulate the deoxy state, and you say which in your methods.

Tap to explore in 3D

As downloadedHeme and proximal histidineThe heme pocket

Does my force field already include heme, NAD+ or FAD?

Check before you parameterize anything. We grepped the residue databases of the force fields that ship with GROMACS and of the current CHARMM36 port from the MacKerell lab CHARMM36 page (the February 2026 release, CGenFF 5.0 build). Here is what is defined as a ready residue:

Run the same check on your own copy, because releases change:

grep -nE '^\[ (HEME|NAD|NADH|NAI|FAD|FADR) \]' charmm36-feb2026_cgenff-5.0.ff/*.rtp

The CHARMM36 NAD+ and NADH parameters come from Pavelites, Gao, Bash and MacKerell, “A Molecular Mechanics Force Field for NAD+, NADH, and the Pyrophosphate Groups of Nucleotides” (J. Comput. Chem. 1997, 18, 221 to 239). The heme residue comes from the CHARMM stream file toppar_all36_prot_heme.str, and the flavins from CGenFF. If you are still choosing a force field, our guide on how to choose a force field in GROMACS covers the trade-offs; for cofactor work, built-in coverage is a strong argument for CHARMM36.

Which redox and protonation state should I simulate?

Pick the state the enzyme is in for the step you study, not the state the crystallographers happened to model. A dehydrogenase caught before hydride transfer holds NAD+; after transfer it holds NADH. Use the matching residue: NAD for NAD+ and NADH or NAI for NADH in the CHARMM36 port, FAD or FADR for oxidized or fully reduced flavin.

Know the net charges, because they decide how many counter-ions gmx genion adds. Summing the partial charges in the February 2026 CHARMM36 port gives:

• NAD (NAD+): −1 (two phosphates at −2, the nicotinamide ring at +1)
• NADH and NAI: −2
• FAD and FADR: −2
• HEME: −2

Also check the protein side: the histidine that binds the heme iron and any residue that hydrogen-bonds the nicotinamide or isoalloxazine ring. Set those histidines explicitly (HSD, HSE or HSP in CHARMM36) instead of trusting automatic assignment.

How do I build the topology when the cofactor is a built-in residue?

This is path one, and it uses no external parameters at all.

1. Match the residue name. The PDB residue name must match the .rtp entry. NAD, NAI and FAD already match. For heme, the CHARMM36 port’s aminoacids.r2b maps HEM to HEME, so pdb2gmx handles that rename for you.
2. Compare atom names. Open the [ NAD ] (or [ FAD ], [ HEME ]) block in the .rtp and compare it atom by atom with your PDB. The heavy-atom names largely follow PDB conventions, but check every one; fix a mismatch with a text edit of the PDB, never by changing the .rtp.
3. Let pdb2gmx add hydrogens. The port has hydrogen database entries for HEME, NAD, NADH, NAI, FAD and FADR, so -ignh can drop crystal hydrogens and rebuild them consistently.
4. Run pdb2gmx with the cofactor in the file:

gmx pdb2gmx -f complex.pdb -o complex.gro -p topol.top \
  -ff charmm36-feb2026_cgenff-5.0 -water tip3p -ignh -his

-his makes pdb2gmx ask you about each histidine, which is where you set the iron-binding one. Every other flag is explained in our pdb2gmx tutorial. From here, solvation, ions and minimization proceed exactly as for any protein.

What do I do about the bond between heme iron and histidine?

Decide whether the iron is bonded to the protein or held only by non-bonded forces. pdb2gmx creates the special bonds listed in specbond.dat. The default file shipped with GROMACS includes lines such as:

HIS	NE2	1	HEM 	FE	1	0.2	HIS1	HEME
CYS	SG	1	HEM 	FE	2	0.25	CYS2	HEME
MET	SD	1	HEM 	FE	1	0.24	MET	HEME

Each line names the two atoms, the distance in nm used to detect the link, and the residue names applied afterwards. Two things trip students up. First, both residues must be in the same molecule type, so you need -merge all (or -merge interactive) when the heme sits in its own chain. Second, the rename column says HIS1, a name from GROMOS-style force fields; if your .rtp has no HIS1, copy specbond.dat into your working directory, where GROMACS looks first, and edit the rename to a histidine your force field defines.

A bond is only as good as its angle and dihedral parameters. If the force field lacks them, grompp stops on missing bonded types. That is the point to stop guessing: our post on metal ions in docking and MD covers bonded, non-bonded and restrained models of metal sites, and the safest choice is the one a published study of your protein family used.

What if my force field has no residue for the cofactor?

This is path two. You parameterize the cofactor like a ligand, with extra care for its size and charge. Two standard routes exist, and you pick the one that matches your protein force field.

CHARMM36 + CGenFF. Submit a mol2 with all hydrogens and the correct protonation to the CGenFF program, then convert the .str file with the MacKerell lab’s cgenff_charmm2gmx_py3_nx2.py:

python cgenff_charmm2gmx_py3_nx2.py LIG lig.mol2 lig.str charmm36-feb2026_cgenff-5.0.ff

The CGenFF parameter version in your .ff folder must match the CGenFF program version that wrote the .str; the MacKerell page ships separate 4.6 and 5.0 ports for exactly this reason. Read the penalty scores in the .str. The CGenFF documentation treats penalties below 10 as acceptable, 10 to 50 as needing basic validation, and above 50 as poor. Large cofactors with phosphates and fused rings often land in the higher range, so check them. The force field itself is described in Vanommeslaeghe et al., J. Comput. Chem. 2010, 31, 671 to 690.

AMBER + GAFF2 via ACPYPE. With an AMBER protein force field, either take the cofactor from the AMBER parameter database (its cofactor list includes heme, NAD+, NADH, NADP+, NADPH, FMN and FADH−) or generate GAFF2 parameters:

acpype -i lig.mol2 -c bcc -n -2 -a gaff2

Set -n to the cofactor’s real net charge. Never mix GAFF cofactor parameters into a CHARMM36 protein, or CGenFF into AMBER. Our full walkthrough of both routes is in how to parameterize a ligand for GROMACS.

Which route should I choose?

How do I merge the cofactor topology and set index groups?

For path one, pdb2gmx writes the cofactor into topol.top for you. For path two, follow the same merge as for any ligand, covered step by step in our protein-ligand MD tutorial: include the cofactor’s parameter file after the force field, include its .itp after the protein, and add one line to [ molecules ].

Then fix the coupling groups. A cofactor treated as its own group of a few dozen atoms gives noisy temperature coupling, so couple it with the protein. Build the group with make_ndx:

gmx make_ndx -f em.gro -o index.ndx
> 1 | r HEME
> q

Group 1 is Protein by default; check the numbers printed on your screen. Then in the .mdp set tc-grps = Protein_HEME Water_and_ions, using the exact name make_ndx printed. Our make_ndx index groups tutorial covers the selection syntax.

Finally, restrain the cofactor during NVT and NPT. If it is part of the protein molecule, pdb2gmx’s posre.itp covers it. If it is a separate molecule, generate restraints:

gmx genrestr -f lig.gro -o posre_lig.itp -fc 1000 1000 1000

and include them in the cofactor’s .itp inside an #ifdef POSRES_LIG block, with define = -DPOSRES -DPOSRES_LIG in the equilibration .mdp files.

What are the common errors and how do I fix them?

• Residue ‘HEM’ (or ‘NAD’, ‘FAD’) not found in residue topology database. The name does not match any .rtp entry, or the force field you picked lacks it. Grep the .rtp files as shown above; rename, switch to the CHARMM36 port, or parameterize.
• Atom X in residue NAD not found in rtp entry. An atom name in your PDB differs from the .rtp. Rename it in the PDB. If it is a hydrogen, use -ignh.
• Missing atoms in the cofactor. Crystal structures sometimes omit part of the adenosine or ribityl tail. Rebuild them from the ideal coordinates on the cofactor’s RCSB ligand page before running pdb2gmx.
• Net charge not an integer after merging a CGenFF or ACPYPE topology. The cofactor charges were rounded or the wrong protonation state was submitted. Fix the input structure and regenerate; do not hand-edit a charge to force an integer.
• No default Bond, Angle or Proper Dih. types around the iron. A special bond was created without matching bonded parameters. Remove the bond and use a non-bonded model, or take parameters from a published heme model.
• The cofactor drifts out of the pocket in NVT. Restraints were not applied (missing -DPOSRES_LIG) or the protonation state is wrong. Compare cofactor RMSD with and without restraints, and check the starting contacts.

More general fixes are collected in our list of common GROMACS errors and how to fix them.

How do I check that the cofactor behaves after the run?

Three checks catch most setup mistakes. Compute the cofactor’s RMSD after fitting on the protein backbone (gmx rms with a backbone fit group and a cofactor output group). Track the distance between the iron and the histidine NE2, or between the nicotinamide C4 and the substrate, with gmx distance. And watch the hydrogen bonds that hold the adenine or flavin ring. A cofactor that keeps its crystal contacts over the production run is a good sign; one that loses them within the first nanoseconds usually points back to a charge or state error, not to interesting biology.

FAQ

Can I just delete the cofactor and simulate the apo protein?

Only if your question is about the apo state. Removing heme or FAD from a protein that needs it creates an empty cavity, and the pocket you care about may collapse during the run.

Does the CHARMM36 port in GROMACS include heme?

Yes. The February 2026 port defines HEME in aminoacids.rtp, taken from the CHARMM heme stream file, and maps the PDB name HEM to it. The bundled CHARMM27, GROMOS 54A7 and GROMOS 43A1 also define HEME.

Which PDB code is NADH?

NAI. NAD is the oxidized form, NAD+. The CHARMM36 port has entries for both NAD and NAI, so the PDB names can be used directly.

What CGenFF penalty is too high for a cofactor?

The CGenFF documentation calls penalties above 50 poor and suggests validation for 10 to 50. For cofactors that already exist in CHARMM36, use the built-in residue instead of a CGenFF estimate.

Should the cofactor be in its own temperature coupling group?

No. Couple it together with the protein, for example Protein_HEME, and keep solvent and ions as the second group.

Written by the StemSkills Lab team. Force-field coverage was checked against the CHARMM36 February 2026 port and the force fields bundled with GROMACS 2024; re-check if you use a different release.