Virtual Screening with AutoDock Vina: Screen a Compound Library Against a Protein Target

Where to get compound libraries

Several databases provide free compound libraries in formats ready for docking. Choose based on your use case:

For this tutorial, download a small focused library from ZINC20. Go to zinc20.docking.org → Catalogs → In-Stock → filter for drug-like (MW 250–500, LogP −1 to 5) → download 500 compounds as a single SDF file. Save it as library_raw/library.sdf.

Batch-prepare ligands with Open Babel

You need a separate PDBQT file for every compound in the library. Open Babel can do this in a single command by splitting the SDF and converting each entry:

cd ~/docking/virtual_screen

# Convert entire library SDF to individual PDBQT files
obabel library_raw/library.sdf \
  -O library_pdbqt/compound_.pdbqt \
  --gen3d \
  -p 7.4 \
  --partialcharge gasteiger \
  -m \
  2>> library_pdbqt/conversion.log

# Check how many were created
ls library_pdbqt/*.pdbqt | wc -l

The -m flag tells Open Babel to split the input into individual output files, naming them compound_1.pdbqt, compound_2.pdbqt, etc. The 2>> conversion.log captures any conversion warnings without cluttering your terminal.

Inspect the log for failures — some compounds will fail 3D conformer generation (usually highly constrained ring systems or salts). This is normal. A 5–10% failure rate is expected; above 20% suggests a problem with your input file.

Copy your prepared receptor into place:

cp ~/docking/cox2_ibuprofen/receptor/5KIR_receptor.pdbqt receptor/

Create a config template. In virtual screening you’ll use the same receptor, grid box, and search settings for every compound — only the ligand path changes. Save this as config_template.txt:

# Receptor — same for every docking run
receptor = receptor/5KIR_receptor.pdbqt

# Grid box — binding site center from previous tutorial
center_x = 15.234
center_y = -8.441
center_z = 22.178
size_x   = 20
size_y   = 20
size_z   = 20

# Search settings — lower exhaustiveness for speed in VS
exhaustiveness = 8
num_modes      = 1
energy_range   = 3

Note num_modes = 1 — in virtual screening you only need the top-ranked pose per compound. Saving 9 poses per compound wastes disk space and slows down result parsing without adding value at the screening stage.

The core of virtual screening automation is a loop that: writes a per-compound config file, calls Vina as a subprocess, captures the output score, and logs everything. Save the following as run_vs.py:

import subprocess
import os
import glob
import csv
import re
from pathlib import Path
from tqdm import tqdm

# ── Configuration ────────────────────────────────────────────────
RECEPTOR    = "receptor/5KIR_receptor.pdbqt"
LIGAND_DIR  = "library_pdbqt"
OUTPUT_DIR  = "output"
RESULTS_CSV = "results/all_scores.csv"
VINA_BIN    = "vina"          # or full path: /usr/local/bin/vina
N_CPUS      = 4               # parallel Vina jobs

# Grid box — update with your binding site coordinates
GRID = {
    "center_x": 15.234, "center_y": -8.441, "center_z": 22.178,
    "size_x":   20,      "size_y":   20,      "size_z":   20,
    "exhaustiveness": 8,  "num_modes": 1,  "energy_range": 3,
}
# ─────────────────────────────────────────────────────────────────

def write_config(ligand_path, output_path, config_path):
    """Write a per-compound Vina config file."""
    lines = [
        f"receptor = {RECEPTOR}\n",
        f"ligand   = {ligand_path}\n",
        f"out      = {output_path}\n\n",
    ]
    for k, v in GRID.items():
        lines.append(f"{k} = {v}\n")
    Path(config_path).write_text("".join(lines))


def parse_score(output_pdbqt):
    """Extract the best binding affinity from a Vina output PDBQT."""
    try:
        text = Path(output_pdbqt).read_text()
        match = re.search(r"REMARK VINA RESULT:\s+([-\d.]+)", text)
        return float(match.group(1)) if match else None
    except:
        return None


def dock_one(ligand_pdbqt):
    """Run a single Vina docking job. Returns (compound_name, score)."""
    name        = Path(ligand_pdbqt).stem
    out_pdbqt   = f"{OUTPUT_DIR}/{name}_out.pdbqt"
    config_path = f"{OUTPUT_DIR}/{name}_config.txt"

    write_config(ligand_pdbqt, out_pdbqt, config_path)

    result = subprocess.run(
        [VINA_BIN, "--config", config_path],
        capture_output=True, text=True, timeout=300
    )

    if result.returncode != 0:
        return name, None  # docking failed

    score = parse_score(out_pdbqt)

    # Clean up config file to save space
    os.remove(config_path)
    return name, score


def main():
    os.makedirs(OUTPUT_DIR, exist_ok=True)
    os.makedirs("results", exist_ok=True)

    ligands = sorted(glob.glob(f"{LIGAND_DIR}/*.pdbqt"))
    print(f"Found {len(ligands)} ligands to dock.")

    scores = []
    failed = []

    for ligand in tqdm(ligands, desc="Docking"):
        name, score = dock_one(ligand)
        if score is not None:
            scores.append({"compound": name, "score_kcal_mol": score})
        else:
            failed.append(name)

    # Write results CSV
    scores.sort(key=lambda x: x["score_kcal_mol"])  # most negative first
    with open(RESULTS_CSV, "w", newline="") as f:
        writer = csv.DictWriter(f, fieldnames=["compound", "score_kcal_mol"])
        writer.writeheader()
        writer.writerows(scores)

    print(f"\nDone. {len(scores)} succeeded, {len(failed)} failed.")
    print(f"Results saved to {RESULTS_CSV}")
    print(f"\nTop 5 hits:")
    for row in scores[:5]:
        print(f"  {row['compound']:30s}  {row['score_kcal_mol']:.1f} kcal/mol")

    if failed:
        print(f"\nFailed compounds ({len(failed)}):")
        for name in failed[:10]:
            print(f"  {name}")


if __name__ == "__main__":
    main()

Run it:

python run_vs.py

You’ll see a progress bar counting through your library. On a modern laptop with exhaustiveness = 8, expect roughly 1–3 compounds per minute. A 500-compound library takes 3–8 hours on a single CPU — which is why this tutorial covers parallelization below.

Example virtual screening output. Your compounds and scores will differ.

The CSV gives you a ranked list of all docked compounds. Now you need to filter it down to a shortlist worth inspecting visually. The standard approach is a two-stage filter: a score cutoff, then manual pose inspection of the survivors.

Save this as filter_results.py:

import pandas as pd

# Load all results
df = pd.read_csv("results/all_scores.csv")
print(f"Total docked: {len(df)}")
print(f"Score distribution:\n{df['score_kcal_mol'].describe()}\n")

# Apply score cutoff — take top 10% or score below -9.0, whichever is smaller
cutoff  = min(df['score_kcal_mol'].quantile(0.10), -9.0)
hits    = df[df['score_kcal_mol'] <= cutoff].copy()
hits    = hits.sort_values('score_kcal_mol')

print(f"Score cutoff applied: {cutoff:.1f} kcal/mol")
print(f"Hits passing cutoff:  {len(hits)}\n")
print(hits.to_string(index=False))

# Save filtered hit list
hits.to_csv("results/hits.csv", index=False)
print("\nHit list saved to results/hits.csv")

python filter_results.py

Your ranked hit list will look something like this:

The score cutoff gets you from hundreds of compounds down to tens. The next cut — manual pose inspection — gets you from tens down to the handful you’d actually take to the lab. The funnel looks like this in practice:

Load your receptor and top-scoring compound outputs together in PyMOL. This script loads all your top hits at once for efficient visual comparison:

# Load receptor
load receptor/5KIR_receptor.pdbqt, receptor
hide everything, receptor
show cartoon, receptor
color gray80, receptor

# Load top hits — repeat for each compound in your hits.csv
load output/compound_142_out.pdbqt, hit_142
load output/compound_089_out.pdbqt, hit_089
load output/compound_311_out.pdbqt, hit_311

# Show all hits as colored sticks
show sticks, hit_142 or hit_089 or hit_311
color cyan,   hit_142
color magenta, hit_089
color yellow,  hit_311

# Show binding site residues within 5 Å of any hit
select bs, receptor and (byres (hit_142 or hit_089 or hit_311) around 5)
show sticks, bs
zoom bs

For each hit in your shortlist, work through the same checklist from the single-ligand tutorial: is the ligand inside the pocket, are there hydrogen bonds to known key residues, does the binding mode make chemical sense? Discard any compound whose top-ranked pose is outside the binding site or makes no meaningful contacts — these are false positives driven by the scoring function.