How to Handle Metal Ions in Molecular Docking and MD (Zinc, Magnesium and Iron, Step by Step)

Keep the metal ion whenever it sits in or beside the binding site, because deleting it removes the strongest electrostatic feature of the pocket. For docking, keep the ion in the rigid receptor and switch to the AutoDock4 scoring function, or to AutoDock4Zn for zinc. For MD, decide between a nonbonded point-charge ion and a bonded metal-centre model before you build the topology.

Metalloproteins break the standard tutorial workflow. Every docking and molecular dynamics walkthrough assumes a protein made of the twenty standard residues plus water, so the moment your structure has a zinc, a magnesium or a haem iron sitting in the active site, the instructions stop matching what you see on screen. The two failure modes are predictable. Students either delete the HETATM line because a preparation script complained about it, or they keep the ion and let a force field treat it as a bare point charge that wanders out of the pocket by nanosecond three.

This guide from the StemSkills Lab team (10+ years in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling) walks through the decision and both workflows with the real commands. It sits inside our pillar guide on learning molecular docking, and the upstream structure-preparation steps are covered in our staged computational biology skills roadmap.

How common are metal ions in the structures you will actually dock?

Common enough that you will meet one early. A search of the RCSB Protein Data Bank on 31 August 2026 returns 27,809 released entries containing a zinc ion and 28,673 containing magnesium, out of 259,129 released structures in total. Roughly one structure in ten in the PDB carries a zinc. Kinases (magnesium plus ATP), matrix metalloproteinases and carbonic anhydrases (zinc), and cytochromes (haem iron) are among the most heavily studied drug targets in the database, so the odds that your thesis target is metal-free are lower than most students expect.

That matters because the ion is not decoration. A divalent cation in a pocket is a formal +2 charge with a short, geometrically strict coordination shell, usually 2.0 to 2.2 Angstrom to the coordinating nitrogen, oxygen or sulfur atoms. If your ligand is designed to chelate that ion, as hydroxamates do in metalloproteinase inhibitors and sulfonamides do in carbonic anhydrase inhibitors, then the metal contact is the binding mode. Remove it and you are docking into a different protein.

Should you keep or delete the metal ion?

Ask one question: is the ion structural or catalytic, or is it a crystallisation artefact? Three checks settle it.

1. Distance to the binding site. Open the structure and measure. An ion within about 6 Angstrom of the co-crystallised ligand or the pocket you defined is part of the site. Keep it.
2. Coordination environment. A catalytic or structural metal is held by protein side chains, typically three or four of His, Cys, Asp or Glu, in a recognisable tetrahedral or octahedral arrangement. An ion coordinated only by water and a single surface carboxylate, sitting far from any pocket, is almost always from the crystallisation buffer. Delete that one.
3. The literature and the PDB entry. Read the paper the structure came from and the entry’s own annotation. If the enzyme is described as a zinc metalloenzyme, the zinc stays, even if your docking software would rather it did not.

Do this check before anything else in preparation, because the ion also changes how you protonate the site. Our guide on setting protonation states and adding hydrogens covers the general rules, but the metal-specific one is worth stating plainly here: a histidine that coordinates a metal must not carry a hydrogen on the coordinating nitrogen. If it does, you have put a proton and a +2 cation on the same lone pair, and both docking scores and MD stability will suffer for it.

Why does AutoDock Vina handle metal ions badly by default?

Because the default Vina scoring function was never built for them. It is an empirical function fitted to typical organic protein-ligand interactions, with hydrophobic and hydrogen-bond terms and no explicit electrostatic term. Metal ions are typed as metal donors and scored with the same short-range polar term used for hydrogen bonds. That approximation is reasonable for a distant surface ion and poor for a ligand whose whole binding mode is a coordination bond.

You have three better options, in increasing order of effort.

How do you dock into a zinc metalloprotein with AutoDock4Zn?

AutoDock4Zn is the specialised extension of the AutoDock4 force field for zinc sites, published by Santos-Martins, Forli, Ramos and Olson in Journal of Chemical Information and Modeling in 2014 (doi:10.1021/ci500209e). It adds tetrahedral zinc pseudo atoms, typed TZ, at the positions a ligand donor atom would occupy if it completed the zinc’s tetrahedral coordination shell, so both the energy and the geometry of the interaction are described. The official AutoDock Vina documentation records that the force field was, in its words, “calibrated on a data set of 292 crystal complexes containing zinc”.

The workflow below follows the official Vina tutorial on docking with zinc metalloproteins, which uses PDB entry 1S63. It needs the ADFR software suite, the Meeko package and the AD4Zn.dat parameter file shipped in the AutoDock-Vina example directory.

Step 1. Prepare the receptor and add the pseudo atoms. Start from a receptor file that already has all hydrogens, then add TZ atoms around the zinc.

mk_prepare_receptor.py -i proteinH.pdb -o protein -p
pythonsh <script_directory>/zinc_pseudo.py -r protein.pdbqt -o protein_tz.pdbqt

The second command reports how many pseudo atoms it wrote, for example Wrote 1 TZ atoms on protein_tz.pdbqt. If it reports zero, the zinc never made it into the PDBQT and everything downstream is meaningless.

Step 2. Prepare the ligand. Add hydrogens, then convert to PDBQT with Meeko.

scrub.py 1s63_ligand.sdf -o 1s63_ligandH.sdf
mk_prepare_ligand.py -i 1s63_ligandH.sdf -o 1s63_ligand.pdbqt

Step 3. Generate the zinc-aware grid maps. A dedicated script writes a grid parameter file that points at AD4Zn.dat and carries the pairwise potentials for the TZ and Zn atom types.

pythonsh <script_directory>/prepare_gpf4zn.py -l 1s63_ligand.pdbqt -r protein_tz.pdbqt \
  -o protein_tz.gpf -p npts=40,30,50 -p gridcenter=18,134,-1 \
  -p parameter_file=AD4Zn.dat
autogrid4 -p protein_tz.gpf -l protein_tz.glg

Step 4. Dock with the AutoDock4 scoring function. Once the force field is baked into the maps, you pass the maps rather than a receptor file.

vina --ligand 1s63_ligand.pdbqt --maps protein_tz --scoring ad4 \
  --exhaustiveness 32 --out 1s63_ligand_ad4_out.pdbqt

For this tutorial system the documentation states the best pose should come out at roughly -13 kcal/mol and reproduce the crystallographic pose. Treat that as your sanity check: if you can redock the native ligand of a metalloprotein and recover the crystal pose, your setup is sound. The same redocking logic applies to every docking protocol, metal or not, and is the basis of the validation step every examiner asks about.

What about magnesium, iron and other metals in docking?

There is no shipped equivalent of AutoDock4Zn for magnesium or iron, so the practical answer is the AutoDock4 scoring function without pseudo atoms. Keep the ion in the receptor, build AutoGrid maps so the electrostatic and desolvation terms are active, and run with --scoring ad4. Two habits keep this honest.

• Check the atom type in the PDBQT. Open the receptor file and find the ion line. The final column is the AutoDock atom type and must read Zn, MG, FE or MN, not C. A metal silently typed as carbon is the single most common cause of nonsense metal-site docking results.
• Never dock the ion as part of the ligand. Magnesium in a kinase belongs to the receptor, together with any coordinating waters you have decided to keep, and the ATP-competitive ligand is docked against that assembly.

Haem iron deserves its own warning. A haem group is a large organic cofactor with an iron at its centre, and preparation scripts frequently mangle it or drop it. Confirm that the whole porphyrin survived preparation, and the FE line with it, before you trust any pose. If your preparation step is dropping residues in general, our guide on fixing missing residues and loops in a PDB structure covers the upstream repair work.

How do you handle a metal ion in a GROMACS MD simulation?

Here the choice is explicit and you must make it deliberately: a nonbonded model or a bonded model.

In the nonbonded model, the metal is a simple ion, a +2 point charge with Lennard-Jones parameters and no bonds to the protein. This is what you get for free when the ion is defined in your force field’s ion parameters. It is fast, it needs no extra work, and it allows ligand exchange in principle. Its weakness is that nothing enforces coordination geometry, so the ion can and often does drift out of the site during an unrestrained run.

In the bonded model, explicit bonds, angles and charges describe the metal centre. You build these with a parameter builder rather than by hand. MCPB.py, published by Li and Merz in Journal of Chemical Information and Modeling in 2016 (doi:10.1021/acs.jcim.5b00674), generates AMBER parameters for a metal centre, which you then convert to GROMACS format with ACPYPE. For zinc specifically, the Zinc AMBER Force Field of Peters and colleagues (doi:10.1021/ct1002626) provides pre-derived parameters for common coordination motifs.

For a first MSc project on a metalloprotein, the restrained nonbonded model is usually the defensible middle path. State the choice and its limitation in your methods section rather than hiding it.

Building the topology

Before running gmx pdb2gmx, check what your force field actually knows about your ion instead of assuming. Naming differs between ports: the CHARMM36 port names the zinc ion residue ZN2, while AMBER ports generally use ZN.

grep -riE 'ZN|ZN2' /path/to/charmm36.ff/ions.itp /path/to/charmm36.ff/merged.rtp

Then rename the residue and atom in your PDB to match exactly what you found, and build the topology as usual. The general pdb2gmx workflow is covered in our guide on generating a GROMACS topology with pdb2gmx, and the official gmx pdb2gmx documentation lists every flag, including -his for choosing histidine tautomers interactively. Use that flag on metalloproteins. It is how you stop pdb2gmx from protonating the nitrogen that is supposed to be coordinating your zinc.

If your ligand also needs parameters, that is a separate job, handled in our guide to parameterising a ligand for GROMACS.

Checking that the metal stayed put

Never assume the site survived. Measure a coordination distance across the trajectory with gmx distance, which computes distances between pairs of positions within a selection.

gmx distance -s md.tpr -f md_noPBC.xtc \
  -select 'resname ZN and name ZN plus resid 94 and name NE2' \
  -oall zn_his94.xvg

Plot the result. A well-behaved bonded or restrained model gives a flat line near 2.0 to 2.2 Angstrom. A nonbonded ion that has escaped shows an unmistakable staircase upward. To measure the minimum distance from the ion to a whole group of atoms rather than a fixed pair, the GROMACS documentation directs you to gmx pairdist instead.

What errors will you actually hit, and how do you fix them?

• Wrote 0 TZ atoms on protein_tz.pdbqt. The zinc is not in the receptor PDBQT. Re-run receptor preparation and confirm the ion line is present before adding pseudo atoms.
• AutoGrid cannot find the parameter file. The parameter_file keyword in the GPF is read relative to the working directory, so keep a local copy of AD4Zn.dat beside your input files or give the full path.
• Residue 'ZN' not found in residue topology database. The residue name in your PDB does not match the force field. Grep the force field’s ion definitions and rename to match, for example ZN2 for CHARMM36.
• Atom ZN in residue ZN2 was not found in rtp entry. The residue name is right but the atom name is not. Both columns must match the force field entry, and PDB column alignment matters.
• The system explodes near the metal during minimisation or NVT. Almost always protonation. Check that no coordinating histidine, cysteine or aspartate carries a hydrogen pointing at the ion, and rebuild hydrogens if it does.
• The ion drifts away during production MD. Expected behaviour for a nonbonded ion. Add distance restraints or move to a bonded model, and say which you did in your methods.
• Docking scores look implausibly good or bad. Check the atom type column in the receptor PDBQT first. A mistyped metal is more likely than a real discovery.

What does a defensible metalloprotein workflow look like end to end?

1. Download the structure and decide, with a reason you can write down, which metal ions and waters to keep. Our guide on finding a protein’s binding site helps you draw that boundary.
2. Repair missing residues, then set protonation states with the metal coordination in mind.
3. Prepare the receptor keeping the ion, and verify its atom type in the PDBQT.
4. Redock the native ligand with the AutoDock4 scoring function, or AutoDock4Zn for zinc, and confirm you recover the crystal pose.
5. Dock your compounds with the validated protocol, then inspect poses for sensible coordination geometry rather than for the best number alone.
6. For MD, choose nonbonded, restrained or bonded, build the topology against a force field you have checked, and monitor the coordination distance for the whole run.
7. Report the model and its limitation explicitly in your methods section.

Our step-by-step walkthrough of preparing a protein and ligand for docking covers steps 2 and 3 for the metal-free case, and this article is the layer you add on top when a metal is present.

FAQ

Can I just delete the zinc and dock anyway?

Only if the zinc is far from your binding site and clearly a crystallisation artefact. If it is in the pocket, deleting it removes the dominant electrostatic feature and your poses will not reflect the real binding mode, particularly for chelating ligands such as hydroxamates and sulfonamides.

Does AutoDock Vina support metals at all?

Yes, but weakly with the default scoring function, which types metals as polar donors and has no explicit electrostatic term. Vina 1.2 can run the AutoDock4 scoring function through --scoring ad4, which is the route to proper electrostatic and desolvation treatment, and to AutoDock4Zn for zinc.

Do I need to keep the water molecules coordinating the metal?

Often yes. If a water completes the coordination shell and mediates contact with the ligand in the crystal structure, keeping it usually improves the pose. Keep only waters you can justify from the structure, and state which ones you kept.

Is a bonded model always better than a nonbonded one?

No. A bonded model fixes the coordination pattern you built into it, so it cannot describe ligand exchange or a change in coordination number. It is the right choice when the metal centre is the subject of the study, and unnecessary overhead when the metal is simply holding a fold together.

How do I cite this properly in my thesis?

Cite the primary methods papers, not tutorials: the AutoDock4Zn paper for zinc docking, MCPB.py or ZAFF for bonded metal parameters, and the AutoDock Vina 1.2 paper (doi:10.1021/acs.jcim.1c00203) for the docking engine itself. Then state your model choice and its limitation in one sentence.