Protein-Ligand Complex PDB from an AutoDock Vina Pose

To make a protein-ligand complex PDB from an AutoDock Vina result, extract the top pose (MODEL 1) from the output PDBQT, convert it to SDF or PDB with correct bond orders and hydrogens (Meeko’s mk_export.py or Open Babel), then combine it with the original receptor PDB in PyMOL, ChimeraX or a text editor. Rename the ligand residue and give it its own chain before saving.

Docking ends with a file called something like out.pdbqt, and the next step (an MD run, a 2D interaction diagram, a figure for your thesis) needs a single PDB holding the protein and the ligand together. Nothing in Vina writes that file for you. This guide, from the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modeling), gives the exact commands, explains why a PDBQT pose is not ready for MD as it stands, and lists the errors that show up later in GROMACS when this step is done carelessly. It sits in our molecular docking learning path, right between interpreting your docking results and running a protein-ligand MD simulation.

What is inside a Vina output PDBQT file?

The output PDBQT holds only the ligand, not the protein. With default settings Vina writes up to 9 poses (num_modes defaults to 9) within 3 kcal/mol of the best one (energy_range defaults to 3), each wrapped in a MODEL and ENDMDL block and ranked from best to worst. The first line inside each model is a remark such as:

MODEL 1
REMARK VINA RESULT:    -8.4      0.000      0.000

The three numbers are the predicted affinity in kcal/mol and two RMSD values (lower and upper bound) from the best pose, which is why MODEL 1 always shows 0.000. Your values will differ; the format will not.

Three properties of this file matter for what comes next:

• No bond orders. PDBQT stores coordinates, partial charges and AutoDock atom types, but not double or aromatic bonds. Any program reading it has to guess them from distances.
• Polar hydrogens only. Vina uses a united-atom scoring model, so the ligand preparation step merged the non-polar hydrogens into their carbons. A force field for MD needs every hydrogen.
• The receptor did not move. In a standard rigid-receptor run, Vina only moves the ligand. The pose coordinates are in the same frame as the receptor file you prepared, so a pose drops straight into the pocket of the original PDB. (Flexible side chains are the exception, covered below.)

The official AutoDock Vina documentation describes the input and output formats, and the program itself is described in Trott and Olson, J. Comput. Chem. 2010, which reported “approximately two orders of magnitude speed-up” over AutoDock 4. Version 1.2 and its Python bindings are described in Eberhardt et al., J. Chem. Inf. Model. 2021.

How do I extract the top pose from out.pdbqt?

Pick the pose first. The best score is MODEL 1, but it is not automatically the pose you should simulate. If you validated your protocol by redocking, or if a known binding mode exists for a related ligand, check that the pose you carry forward makes the same key contacts. Write down which model you chose and why; reviewers ask.

There are three reliable ways to pull one pose out.

Option 1: Meeko’s mk_export.py (best if you prepared the ligand with Meeko)

If you made the ligand PDBQT with mk_prepare_ligand.py, Meeko stored the original molecule’s chemistry in the file, so it can write the poses back out as SDF with correct bond orders:

mk_export.py out.pdbqt -s docked_poses.sdf

This is the form given in the Meeko README. The SDF contains all poses in rank order; the first record is MODEL 1. Because bond orders survive, this is the cleanest route to a ligand you can parameterise for MD.

Option 2: Open Babel (works on any PDBQT)

Open Babel reads each MODEL as a separate molecule, and its -f and -l flags set the first and last molecule to convert. To take only pose 1:

obabel out.pdbqt -O pose1.sdf -f 1 -l 1

To split every pose into its own numbered file instead, use -m:

obabel out.pdbqt -O pose.sdf -m

Open Babel has to perceive bond orders from geometry here, so open the result and check every ring and every carbonyl, carboxylate and nitro group. Our Open Babel format conversion guide covers the flags in more depth, and the tool itself is described in O’Boyle et al., J. Cheminform. 2011.

Option 3: vina_split

Vina ships a small utility that writes each model to its own PDBQT file:

vina_split --input out.pdbqt

Run ls afterwards to see the file names it produced. The pieces are still PDBQT, so you will convert one with Open Babel or Meeko anyway. If you docked with flexible side chains, the output also contains those residues, and vina_split separates the ligand from the flexible residues for you.

How do I add back the hydrogens and fix bond orders?

For a figure, polar hydrogens are often enough. For MD, the ligand needs every hydrogen at the protonation state you intend to simulate. With Open Babel:

obabel pose1.sdf -O pose1_H.sdf -p 7.4

-p 7.4 adds hydrogens appropriate for pH 7.4 using Open Babel’s built-in rules; -h adds them without any pH model. Both rely on correct bond orders, which is why Option 1 above is worth the effort: adding hydrogens to a ring that was read as non-aromatic produces a wrong molecule with a plausible look.

If you only have the PDBQT and Open Babel guesses badly, take the bond orders from the original ligand file instead. RDKit can assign bond orders from a template molecule (its AssignBondOrdersFromTemplate function takes the original SMILES or SDF as the template and the docked pose as the target). Whatever route you use, compare the final ligand with the molecule you started docking with: same formula, same charge, same stereocentres.

How do I merge the pose and the receptor into one PDB?

Use the receptor PDB you made before converting it to PDBQT (the cleaned file with waters removed and hydrogens added, as in our protein and ligand preparation guide). Do not use the receptor PDBQT for MD: it has lost its non-polar hydrogens and carries AutoDock atom types instead of element-only records.

First convert the hydrogenated pose to PDB:

obabel pose1_H.sdf -O ligand.pdb

In PyMOL

load receptor_clean.pdb, receptor
load ligand.pdb, ligand
alter ligand, chain='L'
alter ligand, resn='LIG'
alter ligand, resi='1'
set retain_order, 1
create complex, receptor or ligand
save complex.pdb, complex

alter changes atom properties in place (see the PyMOL Wiki entry for alter), and retain_order stops PyMOL from re-sorting atoms when it writes the file, which matters later because the ligand atom order must match its topology. Choose a residue name of at most three characters that matches the one your ligand topology will use.

In UCSF ChimeraX

open receptor_clean.pdb
open ligand.pdb
changechains #2 L
combine #1,2 modelId 3 name complex
save complex.pdb models #3

The ChimeraX combine documentation explains how it handles chain ID clashes between models. Rename the ligand residue in ChimeraX or in a text editor afterwards if needed.

With plain text tools

A PDB file is plain text, so the quickest merge is three lines in a terminal:

grep -E '^(ATOM|TER)' receptor_clean.pdb > complex.pdb
grep HETATM ligand.pdb >> complex.pdb
echo END >> complex.pdb

This keeps the atom order exactly, but drops CONECT records, leaves duplicate atom serial numbers, and keeps whatever residue name Open Babel wrote (often UNL). Edit the residue name and chain columns by hand or in PyMOL before you use the file.

Which merge method should I use?

What if I docked with flexible side chains?

Then the receptor did change: the flexible residues have new side-chain positions that belong with that pose. If you prepared the receptor with Meeko and wrote the JSON file, export both parts in one step, as shown in the Meeko README:

mk_export.py out.pdbqt -j my_receptor.json -s lig_docked.sdf -p rec_docked.pdb

Otherwise, take the flexible residues from the split output and replace the same residues in the receptor PDB. Using the rigid receptor with a pose that was docked against moved side chains can leave atoms overlapping, which shows up as huge forces in energy minimisation.

How do I check the complex before using it?

Open complex.pdb and run four quick checks:

1. The ligand sits in the pocket. If it floats outside the protein, the receptor PDB is not the file you docked against (a different download, or a file you aligned or moved after preparing the PDBQT).
2. The chemistry is right. Rings are aromatic where they should be, the formal charge matches your intended protonation, and the hydrogen count matches the molecule you planned to simulate.
3. Chains are distinct. The ligand has its own chain ID and a residue name that is not also used by the protein.
4. No clashes. Measure the closest protein-ligand heavy-atom distances. Contacts much closer than about 2.5 Å between heavy atoms usually mean a receptor mismatch or a flexible-residue problem.

If the complex passes, you can draw it straight away: our guides to publication-quality PyMOL figures and 2D interaction diagrams with LigPlot+ and PLIP both start from this file.

How do I take the complex into GROMACS?

Do not give the whole complex to gmx pdb2gmx. The standard workflow, set out in Justin Lemkul’s GROMACS protein-ligand complex tutorial, is to split it again: run pdb2gmx on the protein only, build the ligand topology separately (CGenFF, ACPYPE or a similar tool, see our ligand parameterisation guide), then join the coordinates and add the ligand to topol.top. The complex PDB is what guarantees the two pieces sit in the same frame.

The single most useful habit: parameterise the docked pose itself, not the original ligand file. Then the atom names, atom order and hydrogen count in the ligand topology already match the coordinates you will place in the box.

Why does GROMACS fail on my docked complex? (troubleshooting)

For the wider list of errors you will meet after this point, see common GROMACS errors and how to fix them. To see where docking-to-MD fits in the whole skill set, follow our computational biology skills roadmap.

Frequently asked questions

Can I open a Vina output PDBQT directly in PyMOL?

Yes. PyMOL reads PDBQT and loads each pose as a separate state, which is fine for looking at poses. For an MD input or a clean figure file, convert the chosen pose and merge it with the receptor PDB as shown above.

Which pose should I use for MD, MODEL 1 or another one?

Start with MODEL 1 unless you have evidence against it, such as a known binding mode from a crystal structure of a related ligand. Many studies simulate the top two or three distinct poses and compare their stability, which is more convincing than relying on a single docking score.

Why does my ligand have so few hydrogens after docking?

Vina’s scoring uses only polar hydrogens, so ligand preparation merged the non-polar ones into their carbons. Add them back with Open Babel (-p 7.4 or -h) after bond orders are correct, or export with Meeko.

Do I need a complex PDB if I only want a 2D interaction diagram?

Yes. LigPlot+ and PLIP both expect one file with the protein and the ligand together, with the ligand as HETATM records under its own residue name.

Does merging change the docked coordinates?

No. Merging only puts the receptor and pose records into one file. As long as you use the receptor file you docked against, the ligand stays exactly where Vina placed it.

Written by the StemSkills Lab team, who teach docking and molecular dynamics to BSc and MSc students and early-career researchers.