How to Extract a Ligand from a PDB File (PyMOL, ChimeraX)

To extract a ligand from a PDB file, select its HETATM residue by its three-letter code and save only those atoms: save lig.pdb, resn MK1 in PyMOL, save lig.pdb selectedOnly true in ChimeraX, or an awk filter on columns 18 to 20. A PDB file stores no bond orders, so restore them from the RCSB SDF before docking or MD.

This guide comes from the StemSkills Lab team, who have more than 10 years of experience in structural bioinformatics, drug design and molecular modeling. Every command below was checked against one real entry, PDB 1HSG (HIV-1 protease bound to the drug indinavir, ligand code MK1, 2.0 Å resolution), so you can reproduce each number yourself. It sits inside our molecular docking learning path and covers the ligand half of step 4 of the computational biology skills roadmap.

Why would I extract the ligand from a PDB file?

The co-crystallized ligand is the only experimentally observed binding pose you have for that pocket, so it does three jobs in a docking or MD project:

• Redocking validation. You dock the extracted ligand back into its own receptor and measure the RMSD between the docked pose and the crystal pose. A top pose within 2 Å is the usual success threshold, and it tells a reviewer your docking settings can reproduce reality before you trust them on new compounds.
• Grid box centring. The ligand’s coordinates tell you exactly where the pocket is. Its centre becomes the centre of your docking box, which is far more defensible than eyeballing it.
• Ligand parameterization for MD. To simulate the complex in GROMACS you need topology files for the ligand, and tools such as CGenFF or ACPYPE need a clean, correctly bonded ligand file as input.

It is also a very common need. A query of the RCSB PDB search API on 13 October 2026 returned 197,631 of 260,626 experimental entries (about 76%) with at least one non-polymer entity, meaning a ligand, ion, cofactor or crystallization additive sits next to the macromolecule in roughly three structures out of four.

Where is the ligand stored inside a PDB file?

The wwPDB format specification (v3.3, coordinate section) is explicit: “Non-polymer chemical coordinates use the HETATM record type.” Protein atoms are ATOM lines; ligands, waters, ions and additives are all HETATM lines. The header also lists every hetero group in HET and HETNAM records, which is the fastest way to learn the ligand’s code.

grep "^HET " 1HSG.pdb
HET    MK1  B 902      45

That one line tells you the residue name (MK1), the chain (B), the residue number (902) and the atom count (45). In 1HSG the 172 HETATM lines split into 45 MK1 atoms plus 127 water oxygens (44 in chain A, 83 in chain B). Each HETATM line has fixed columns: the residue name sits in columns 18 to 20, the chain in column 22 and the x, y, z coordinates in columns 31 to 54.

HETATM 1517  N1  MK1 B 902       9.280  23.763   3.004  1.00 28.25           N

If you have not downloaded the file yet, our guide on how to download a protein structure from RCSB PDB covers the formats and the command-line route. You can also look up any ligand code on its RCSB page, for example rcsb.org/ligand/MK1, to confirm it is the molecule you expect.

Tap to explore in 3D

As downloadedOnly the HETATM ligandPocket around the ligand

How do I extract a ligand in PyMOL?

Type these into the PyMOL command line. They use standard PyMOL commands documented on the PyMOL Wiki save page.

fetch 1hsg
select lig, resn MK1
save mk1_crystal.pdb, lig
save mk1_crystal.sdf, lig

resn selects by residue name. If a ligand code appears more than once (two copies in two chains is common), narrow the selection so you keep one copy: select lig, resn MK1 and chain B. PyMOL picks the output format from the file extension, so the same selection can go to .pdb, .sdf or .mol2.

The receptor side is a separate job. In PyMOL, remove resn HOH deletes waters and save receptor.pdb, polymer writes the protein only. Our posts on removing water molecules from a PDB file and extracting a single chain cover the edge cases, and the full command set is in our PyMOL commands for docking and MD cheat sheet.

How do I extract a ligand in UCSF ChimeraX?

ChimeraX uses : for residue names and a selectedOnly option on save, documented in the ChimeraX save command reference.

open 1hsg
select :MK1
save mk1_crystal.pdb selectedOnly true
save mk1_crystal.mol2 selectedOnly true

Without selectedOnly true ChimeraX writes the whole model, which is the most common reason students end up “docking” a ligand file that still contains the protein. To restrict to one chain, use select /B:MK1.

How do I extract a ligand from the command line?

On Linux, macOS or WSL2 you need no viewer at all. Because the PDB format is fixed-column, filter on the record type and the residue-name columns rather than grepping for the code anywhere in the line:

curl -O https://files.rcsb.org/download/1HSG.pdb
awk 'substr($0,1,6)=="HETATM" && substr($0,18,3)=="MK1"' 1HSG.pdb > mk1_crystal.pdb
wc -l mk1_crystal.pdb

The output file has 45 lines, one per heavy atom, matching the HET record. A plain grep MK1 also works for 1HSG, but on other entries the same three letters can occur in an atom name or a header line, so the column test is the safer habit.

If you prefer a dedicated tool, the pdb-tools suite from the Bonvin lab (pip install pdb-tools) does the same in one pipe:

pdb_fetch 1hsg | pdb_selresname -MK1 > mk1_crystal.pdb

Which extraction method should I use?

Why is my extracted ligand chemically wrong?

Because the PDB format records where atoms are, not how they are bonded. Most X-ray entries contain no hydrogens (at 2.0 Å they are not resolved), and there are no bond orders: no double bonds, no aromatic rings, no formal charges. Every program that reads mk1_crystal.pdb has to guess the chemistry from distances.

We tested this directly with RDKit 2026.09.1. Reading the 45-atom mk1_crystal.pdb with Chem.MolFromPDBFile gives a molecule with zero double or aromatic bonds. Indinavir has three aromatic rings and two amide carbonyls, so every one of those was silently lost. Feed that file to a docking program or a topology generator and you are working with a different molecule.

There are two reliable fixes, and both lean on the wwPDB Chemical Component Dictionary, which holds the curated chemistry for every ligand code.

Fix 1: download the crystal pose as SDF from RCSB

The RCSB ModelServer returns one ligand instance with its crystal coordinates and dictionary bond orders:

curl -o mk1_model.sdf "https://models.rcsb.org/v1/1hsg/ligand?auth_comp_id=MK1&encoding=sdf"

For 1HSG this file has 45 atoms and 49 bonds (11 of them double, in Kekulé form), and the coordinates match the PDB file exactly (the first atom, N1, sits at 9.280, 23.763, 3.004 in both). Add auth_asym_id=B to the query when the ligand appears in more than one chain. Hydrogens are still absent, so add them next.

Fix 2: transfer bond orders with RDKit

When you already have an extracted PDB ligand (from a docking output, an MD frame or an edited model), assign bond orders from the dictionary template using RDKit’s AssignBondOrdersFromTemplate:

curl -O https://files.rcsb.org/ligands/download/MK1_ideal.sdf

from rdkit import Chem
from rdkit.Chem import AllChem

pdb_lig  = Chem.MolFromPDBFile("mk1_crystal.pdb")
template = Chem.MolFromMolFile("MK1_ideal.sdf")
fixed    = AllChem.AssignBondOrdersFromTemplate(template, pdb_lig)
fixed_h  = Chem.AddHs(fixed, addCoords=True)
Chem.MolToMolFile(fixed_h, "mk1_fixed_h.sdf")

In our run the result had 92 atoms (45 heavy plus 47 hydrogens) and its canonical SMILES matched the RCSB ideal structure exactly, while the heavy atoms kept their crystal coordinates. Do not dock the ideal SDF itself for redocking: it has correct chemistry but computed coordinates, not the crystal pose.

Hydrogens added this way are neutral. If your ligand has an ionizable group, set its protonation state deliberately; Open Babel’s -p option adds hydrogens for a chosen pH (obabel mk1_model.sdf -O mk1_h.sdf -p 7.4), and our guide to protonation states and hydrogens for docking and MD explains when that matters.

How do I use the extracted ligand to centre a docking box?

Average the x, y and z coordinates of the ligand’s heavy atoms. One awk command does it on the extracted file:

awk '{x+=substr($0,31,8); y+=substr($0,39,8); z+=substr($0,47,8); n++}
     END {printf "%.3f %.3f %.3f\n", x/n, y/n, z/n}' mk1_crystal.pdb

For MK1 in 1HSG this prints 13.073 22.467 5.557. The ligand spans 14.9 × 10.2 × 8.6 Å, so a box that adds about 10 Å per axis still fits the pocket comfortably. In an AutoDock Vina config file that becomes:

center_x = 13.073
center_y = 22.467
center_z = 5.557
size_x = 25
size_y = 20
size_z = 19

Those sizes are a reasonable starting choice for this ligand, not a rule. Our full guide on setting the AutoDock Vina grid box centre and size explains how to size it for your own system, and the receptor-side steps are in how to prepare a protein and ligand for docking.

How do I use the extracted ligand for MD simulations?

For a protein-ligand simulation the complex leaves the PDB file in two parts. The protein (without the ligand) goes through gmx pdb2gmx, which cannot build an unknown residue like MK1 and stops with a “Residue ‘MK1’ not found in residue topology database” error if you leave it in. The ligand, with correct bond orders and hydrogens, goes to a parameterization tool. Start from mk1_fixed_h.sdf or a hydrogenated mk1_model.sdf, never from the raw 45-atom PDB lines. Our walkthrough on how to parameterize a ligand for GROMACS picks up from exactly this file, and the wider MD sequence lives in our GROMACS learning path.

What goes wrong, and how do I fix it?

• The output file is empty. The residue name is wrong or has different spacing. Run grep "^HET " file.pdb and copy the code exactly. Large entries are distributed only as mmCIF; download the .cif file and use PyMOL or ChimeraX, which read both.
• The ligand file has twice as many atoms as expected. The entry contains two copies (one per chain) or alternate conformations. Add a chain filter (and chain B in PyMOL, /B:MK1 in ChimeraX) and, for altlocs, keep one (remove not alt ''+A in PyMOL keeps blank and A).
• The ligand file still contains the protein. In ChimeraX you saved without selectedOnly true. Re-save with the option and check that the atom count matches the HET record.
• Aromatic rings look like cyclohexanes in the viewer, or the topology tool reports odd valences. Bond orders were lost from the PDB format. Use the ModelServer SDF or the RDKit template fix above.
• RDKit prints “molecule is tagged as 2D, but at least one Z coordinate is not zero”. This came from the ideal SDF header in our test and is harmless; RDKit marks the molecule as 3D and continues.
• AssignBondOrdersFromTemplate raises “No matching found”. The template and the extracted ligand do not have the same heavy atoms, usually because atoms are missing from the crystal structure (check the REMARK 470 missing-atom records) or the wrong copy was extracted. Rebuild or pick a complete copy.
• Redocking RMSD looks huge even though the poses overlap. The RMSD tool is matching atoms in different orders. Compare symmetry-aware, with hydrogens removed from both files.

Frequently asked questions

What is the difference between ATOM and HETATM in a PDB file?

ATOM lines hold standard amino acids and nucleotides. HETATM lines hold everything else: ligands, waters, metal ions, cofactors and crystallization additives such as glycerol or sulfate. Ligand extraction always starts from HETATM lines.

How do I find the three-letter code of the ligand?

Run grep "^HET " file.pdb or open the entry page on RCSB, where the small-molecules section lists every ligand code with its name. Use the code (MK1), not the drug name (indinavir), in your selection.

Should I extract the ligand as PDB, SDF or Mol2?

Extract to whatever is convenient, but hand downstream tools an SDF or Mol2 that carries correct bond orders and hydrogens. A raw PDB ligand is fine for centring a box and as an RMSD reference, not as docking input or for parameterization.

Can I just dock the RCSB ideal SDF instead?

For screening, yes, since docking re-samples the pose anyway. For redocking validation, no, because its coordinates are computed rather than the crystal pose, so an RMSD against it measures nothing useful.

Do I remove ions and cofactors along with the ligand?

Not automatically. A catalytic zinc or a bound NAD is part of the binding site and often belongs in the receptor. Decide residue by residue based on what the paper behind the entry says the site needs.

Once your extraction is clean, you have everything step 4 of the computational biology skills roadmap asks for: a validated pose to redock, a defensible box centre and a ligand ready for parameterization.