AutoDock Vina 分子对接教程：从 PDB 到可写进论文的结果

版本差异

当前工具链与旧教程的差异

多数中文教程仍基于 AutoDockTools（MGLTools 1.5.7）图形界面和 Vina 1.1.2。下表列出照抄这些教程时会在 Vina 1.2.x 上出问题的地方。

安装

安装 Vina 1.2.7、Meeko 与配套工具

Vina 的 pip 包只含 Python 绑定，不含 vina 命令。PyPI 上 vina 1.2.7 只有 Linux x86_64（Python 3.8–3.12）的预编译包，在 macOS 上 pip 会转为源码编译，缺少 Boost 时报 Boost library location was not found!。conda-forge 的 vina 1.2.7 在 Linux 和 macOS 上同时提供 vina、vina_split 命令和 Python 绑定，下面的环境用它。Windows 用户在 WSL 中执行，或下载 GitHub Releases 中的 win.exe。本页实测：下面第 1 段命令在 macOS arm64（Apple M2）上用 micromamba 从零安装成功。

openbabel 和 RDKit 在新流程中只做辅助：RDKit 负责读写 SDF、给 PDB 配体补键级和算 RMSD；Open Babel 用于格式互转和 PLIP 的依赖，不再用来生成对接用的 PDBQT。

# 1) 环境：vina 用 conda-forge 包，它同时提供 vina、vina_split 命令和 Python 绑定
#    （linux-64、linux-aarch64、osx-64、osx-arm64；Python 3.10-3.14）
micromamba create -n dock -c conda-forge python=3.11 rdkit gemmi prody pdbfixer openmm vina=1.2.7 -y
micromamba activate dock
# Meeko 0.8.0 与 molscrub 0.3.0 只在 PyPI 上有；molscrub 0.3.0 用到 joblib 但没有声明依赖，需要手动装
pip install meeko==0.8.0 molscrub==0.3.0 propka joblib
vina --version        # AutoDock Vina v1.2.7

# 2) 不用 conda 时：从 GitHub Releases 下载对应平台的二进制
#    文件名后缀：linux_x86_64、linux_aarch64、mac_aarch64、mac_x86_64、win.exe
curl -fL -o vina https://github.com/ccsb-scripps/AutoDock-Vina/releases/download/v1.2.7/vina_1.2.7_linux_x86_64
chmod +x vina && ./vina --version
#    pip install vina 只装 Python 绑定；PyPI 上 1.2.7 只有 Linux x86_64 的预编译包，
#    macOS 上会转为源码编译并报 Boost library location was not found!

# 3) 可选：口袋预测与互作分析
#    P2Rank 2.5.1：从 https://github.com/rdk/p2rank/releases 下载解压，需 Java，命令为 prank
#    fpocket：micromamba install -c conda-forge -c bioconda fpocket
#    PLIP：conda-forge 装 openbabel 后 pip install plip

受体准备

受体准备：去水、补残基、质子化态、金属与辅因子

Meeko 按模板逐个匹配残基，重原子必须与模板完全一致，氢可有可无。因此受体准备的重点是补全重原子并定好可滴定残基的状态。

配体准备

配体准备：3D 构象、质子化与互变异构

mk_prepare_ligand.py 要求输入含全部氢原子的 3D 结构，推荐 SDF。PDB 格式没有键级信息，官方明确不建议用来准备小分子。

直接用 SMILES 或 2D 结构对接会出错，原因有两个：Vina 在对接中只旋转可旋转键，环保持输入时的构象，椅式/船式错误无法纠正；2D 平面输入是官方 FAQ 列出的对接失败原因之一。molscrub 会用 ETKDGv3 生成 3D 构象、做力场优化，并把六元环的船式改成椅式。

Vina 用联合原子打分，输出中的氢位置没有意义；但输入中的氢决定哪些原子是氢键供体或受体，所以质子化态会直接改变打分。以伊马替尼为例，pH 5–9 范围内 molscrub 会给出吡啶中性和吡啶质子化两个状态，哌嗪胺都质子化。

• 从 PubChem 下载的结构常带盐和反离子，scrub.py 默认只保留最大片段。
• 质子化态不确定时用 --ph_low 和 --ph_high 枚举多个状态，各自对接，报告最好的一个并注明状态。
• 手性中心未定义的 SMILES 要先确定立体构型。molscrub 0.3.0 没有立体异构体枚举选项，未定义的手性中心由 ETKDG 随机确定；需要枚举时先用 RDKit 的 EnumerateStereoisomers 展开成多条 SMILES，再交给 scrub.py。
• 大环分子：Meeko 默认让大环柔性，Vina 按 BRANCH 数计算构象熵罚，分数可能偏差 2 kcal/mol；比较时对所有分子用同一设置，必要时加 --rigid_macrocycles。
• 含硼、硅以外特殊元素或金属的配体，Vina 会报 Atom type ... is not a valid AutoDock type，需要自定义原子类型参数。

对接盒子

对接盒子怎么设：中心与尺寸

盒子越小，搜索越容易收敛；盒子之外的位置不会被搜索。Vina FAQ 的建议是“尽可能小，但不能更小”，超过 30×30×30 Å 时应同时提高 exhaustiveness。

盒子过小：配体部分原子被挤出盒子，Vina 只能给出被截断的构象或明显偏高的分数。盒子过大：同样的 exhaustiveness 下搜索不充分，结果在不同种子之间变化大。Feinstein 与 Brylinski 在 3659 个复合物上测得，立方盒子边长取配体回转半径的 2.857 倍时 Vina 位姿最准；“每维至少 22.5 Å”的通用盒子平均 RMSD 为 4.9 Å，优化后为 4.0 Å。

盲对接（盒子包住整个蛋白）在网络药理学论文中很常见。Che 与 Zhang 统计 2024 年 35 篇相关论文，17 篇使用或疑似使用盲对接；他们引用的 CASF-2016 数据中，Vina 盲对接的 RMSD < 2 Å 比例为 34–47%，对接到指定位点为 90.2%。

本页实测（P2Rank 2.5.1、fpocket 4.2.3，输入为上一节 PDBFixer 处理后的 1IEP A 链）：P2Rank 第 1 名口袋（probability 0.937）的中心距晶体配体中心 3.0 Å，fpocket 第 1 名口袋的中心距 3.4 Å，口袋残基包含 Thr315 和 Met318。--box_enveloping xtal_lig.sdf --padding 5 生成 18.7×26.7×23.5 Å 的盒子；按 2.857 × 回转半径计算，晶体构象给出 19.1 Å，scrub.py 生成的构象给出 18.7 Å。P2Rank 在 8 核机器上单个结构约 1 分钟，输出 CSV 的列名带前导空格，用 pandas 读取时加 skipinitialspace=True。

# 晶体结构
prank predict -f protA_H.pdb -o p2rank_out
# AlphaFold、NMR、冷冻电镜结构：默认模型把 B 因子当特征，AlphaFold 的 B 因子列存的是 pLDDT
prank predict -c alphafold -f AF-model.pdb -o p2rank_out
# 第一行是排名最高的口袋，center_x/y/z 列即盒子中心
head -3 p2rank_out/protA_H.pdb_predictions.csv

# 交叉验证：fpocket 输出在 protA_H_out/ 目录
fpocket -f protA_H.pdb

# 用口袋中心和配体大小定盒子（Feinstein & Brylinski 2015：边长 = 2.857 x 回转半径）
python - <<'EOF'
from rdkit import Chem
from rdkit.Chem import Descriptors3D
m = Chem.MolFromMolFile("lig_scrub.sdf")  # 默认去氢，只算重原子
print("cube edge (A):", round(2.857 * Descriptors3D.RadiusOfGyration(m), 1))
EOF
mk_prepare_receptor.py --read_pdb protA_H.pdb -o rec -p -v \
    --box_center 12.3 45.6 7.8 --box_size 22 22 22

参数

exhaustiveness、num_modes、energy_range 与随机种子

同一 seed 也不保证与别人的结果完全一致：随机数生成器来自 Boost，编译时的 Boost 版本不同，同一 seed 会得到不同轨迹。维护者的建议是各跑约 10 次，比较分数分布是否不同。Che 与 Zhang 引用的数据中，exhaustiveness 从 8 提到 64，中位 RMSD 从 3.37 Å 降到 2.21 Å，耗时约增加 8 倍。

本页实测（1IEP 重对接，--box_enveloping 加 5 Å 外扩得到 18.7×26.7×23.5 Å 盒子，其余条件见下一节）：exhaustiveness 8、16、32、64 各跑 3 个种子（8 和 32 各补到 10 个种子），最佳构象 RMSD 都在 0.77–0.85 Å，分数在 −12.74 至 −12.82 kcal/mol 之间。单次运行的 CPU 时间（user）中位数分别为 109、281、362、821 s，墙钟时间中位数为 70、175、94、280 s（负载干扰大，只作量级参考）。这个盒子紧贴晶体配体，默认 8 已经 10/10 进入 2 Å；官方文档提到默认 8 偶尔找不到正确构象，大盒子或柔性更大的配体仍建议从 32 起步。同一台机器、同一个 vina 二进制、同一 seed 重复运行时，分数和构象完全相同。

完整命令

完整流程：1IEP 伊马替尼重对接

以 Vina 官方示例体系 1IEP（c-Abl 激酶 + 伊马替尼）为例，从 PDB 文件到分数与 RMSD 表。配体从 SMILES 出发，检验的是整套准备流程能否复现晶体构象。官方示例在 Vina 打分下最佳构象约 −13 kcal/mol。

本页实测（AutoDock Vina 1.2.7（conda-forge）、Meeko 0.8.0、molscrub 0.3.0、PDBFixer 1.12、RDKit 2025.09 与 2026.03；macOS arm64，Apple M2 8 核 16 GB；测试期间机器上同时运行其他任务，系统负载约 25–70）：上面 4 个代码块完整跑通。PDBFixer 的输出交给 Meeko 时，A 链所有残基模板匹配成功；scrub.py 在 pH 7.4 只给出一个状态（两个哌嗪 N 质子化，净电荷 +2）。mk_export.py 为每个构象写一条 SDF 记录，不是一个分子带多个构象。

重对接结果：从 SMILES 出发、exhaustiveness 32、种子 1–10，最佳构象分数 −12.74 至 −12.82 kcal/mol，与晶体构象的原位 RMSD 为 0.80–0.85 Å，10 次全部小于 2 Å；第 2 名构象是首尾翻转的结合方式（RMSD 约 13 Å，分数高约 1.4 kcal/mol）。PLIP 显示第 1 名构象再现了与 Met318 主链和 Thr315 侧链的氢键。晶体构象直接打分（--score_only）为 −12.51 kcal/mol，原位局部优化（--local_only）后为 −13.24 kcal/mol，RMSD 0.21 Å。

起始构象会改变分数：在 RDKit 2026.03 的环境中重跑同一脚本，scrub.py 生成的起始构象不同，3 个种子的最佳分数都是 −12.41 至 −12.42 kcal/mol，RMSD 0.84–0.85 Å。交叉对接确认差异来自配体起始构象，与受体加氢无关。Vina 在对接中不改变键长、键角和环构象，所以比较分数时要固定配体准备环境，或把起始构象文件一起保存。

Vina 终端表格最多列出 num_modes 个构象，输出 PDBQT 只写入与最佳分数相差不超过 energy_range 的构象：本例终端列出 9 个，文件中每次只有 2–7 个。rmsd_table.py 按 PDBQT 中的 REMARK VINA RESULT 读分数，与 SDF 一一对应。

set -euo pipefail
SMI='Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1'   # 伊马替尼

# 1. 下载 1IEP，取 A 链蛋白（只保留无 altloc 或 altloc A）和 A 链共晶配体 STI
curl -fsSLO https://files.rcsb.org/download/1IEP.pdb
awk 'substr($0,1,4)=="ATOM" && substr($0,22,1)=="A" && (substr($0,17,1)==" " || substr($0,17,1)=="A")' 1IEP.pdb > protA.pdb
echo END >> protA.pdb
awk 'substr($0,1,6)=="HETATM" && substr($0,18,3)=="STI" && substr($0,22,1)=="A"' 1IEP.pdb > xtal_lig.pdb

# 2. 补缺失重原子，按 pH 7.4 加氢，去掉水和所有异质分子
pdbfixer protA.pdb --output=protA_H.pdb --add-atoms=all --keep-heterogens=none --ph=7.4

# 3. 检查口袋内可滴定残基的 pKa（结果写入 protA_H.pka）
propka3 protA_H.pdb

# 4. 给共晶配体补键级，得到参考 SDF（定盒子和算 RMSD 都用它）
python make_ref.py "$SMI" xtal_lig.pdb xtal_lig.sdf

# 5. 受体 PDBQT + 盒子：以共晶配体为中心，各方向外扩 5 Å
#    需要改质子化态时加 -n，例如 -n A:<残基号>=HIP
mk_prepare_receptor.py --read_pdb protA_H.pdb -o rec -p -v \
    --box_enveloping xtal_lig.sdf --padding 5
#    输出 rec.pdbqt、rec.box.txt（Vina 配置）、rec.box.pdb（在 PyMOL 中查看盒子）

# 6. 配体：从 SMILES 生成 pH 7.4 质子化态和 3D 构象，再转 PDBQT
scrub.py "$SMI" -o lig_scrub.sdf --ph 7.4 --skip_tautomers
#    scrub.py 写出的分子没有名字；不加 --multimol_prefix 时单分子会写成隐藏文件 lig_pdbqt/.pdbqt
mk_prepare_ligand.py -i lig_scrub.sdf --multimol_outdir lig_pdbqt --multimol_prefix lig

# 7. 重对接：3 个随机种子 x exhaustiveness 32
mkdir -p out
for lig in lig_pdbqt/*.pdbqt; do
  name=$(basename "$lig" .pdbqt)
  for seed in 1 2 3; do
    vina --receptor rec.pdbqt --ligand "$lig" --config rec.box.txt \
      --exhaustiveness 32 --num_modes 9 --energy_range 3 --seed "$seed" \
      --out "out/${name}_s${seed}.pdbqt" | tee "out/${name}_s${seed}.log"
    mk_export.py "out/${name}_s${seed}.pdbqt" -s "out/${name}_s${seed}.sdf"
  done
done

# 8. 每个构象的分数与相对晶体构象的 RMSD
python rmsd_table.py xtal_lig.sdf out > redock_rmsd.tsv
sort -t$'\t' -k3,3g redock_rmsd.tsv | head

# make_ref.py：用 SMILES 模板给 PDB 中的配体补上键级
import sys
from rdkit import Chem
from rdkit.Chem import AllChem

smi, pdb_in, sdf_out = sys.argv[1:4]
template = Chem.MolFromSmiles(smi)
lig = Chem.MolFromPDBFile(pdb_in, removeHs=True)
ref = AllChem.AssignBondOrdersFromTemplate(template, lig)
Chem.MolToMolFile(ref, sdf_out)
print("heavy atoms:", ref.GetNumHeavyAtoms())

# rmsd_table.py：原位计算（不叠合）对称性感知的重原子 RMSD
import glob, os, sys
from rdkit import Chem
from rdkit.Chem import rdMolAlign

def heavy_neutral(m):
    # 去氢并去掉形式电荷，避免质子化态不同导致原子无法匹配
    m = Chem.RemoveHs(m)
    for a in m.GetAtoms():
        if a.GetFormalCharge() != 0:
            # 质子化的 N 去氢后保留了 1 个显式 H，只清电荷会报 valence 错误
            a.SetFormalCharge(0)
            a.SetNumExplicitHs(0)
            a.SetNoImplicit(False)
    Chem.SanitizeMol(m)
    return m

ref = heavy_neutral(Chem.MolFromMolFile(sys.argv[1]))
print("file\tpose\tvina_score\trmsd_to_xtal")
for sdf in sorted(glob.glob(os.path.join(sys.argv[2], "*.sdf"))):
    pdbqt = sdf[:-4] + ".pdbqt"
    scores = [float(l.split()[3]) for l in open(pdbqt) if l.startswith("REMARK VINA RESULT")]
    poses = []
    for m in Chem.SDMolSupplier(sdf, removeHs=False):
        if m is None:
            continue
        for conf in m.GetConformers():
            p = Chem.Mol(m, confId=conf.GetId())
            poses.append(p)
    for i, p in enumerate(poses):
        rmsd = rdMolAlign.CalcRMS(heavy_neutral(p), ref)
        print(f"{os.path.basename(sdf)}\t{i+1}\t{scores[i]:.2f}\t{rmsd:.2f}")

# 批量对接：地图只算一次，换配体时复用
from vina import Vina
import glob, os

v = Vina(sf_name="vina", cpu=0, seed=42)
v.set_receptor("rec.pdbqt")
os.makedirs("out_batch", exist_ok=True)
# 不先加载配体时，compute_vina_maps 为力场中的全部原子类型算地图，换配体时不用重算
v.compute_vina_maps(center=[15.19, 53.90, 16.92], box_size=[22, 22, 22])  # 换成 rec.box.txt 中的值

for lig in sorted(glob.glob("lib_pdbqt/*.pdbqt")):
    v.set_ligand_from_file(lig)
    v.dock(exhaustiveness=16, n_poses=20)
    v.write_poses(os.path.join("out_batch", os.path.basename(lig)), n_poses=9, energy_range=3.0, overwrite=True)
    print(lig, v.energies(n_poses=1)[0][0])

验证

结果可信吗：重对接、阳性对照与打分的局限

网络药理学论文常用的 −4.25、−5.0、−7.0 kcal/mol 阈值没有统一出处，不同对接软件的分数也不能套用同一阈值。比较不同大小的分子时，同时报告配体效率 LE = −分数 / 重原子数；0.3 kcal/mol/重原子常被用作参考线，Vina 分数会随分子变大而变得更负。

1. 01

   重对接 RMSD < 2 Å

   把共晶配体对接回它自己的结构，最佳构象与晶体构象的对称性感知重原子 RMSD 应小于 2 Å。这个阈值见于 Vina 原始论文（Trott & Olson 2010）和 PoseBusters 基准；PoseBusters 中 Vina 在 Astex Diverse 85 个体系上达到该标准的比例为 58%。RMSD 要在受体坐标系内原位计算：RDKit 的 CalcRMS 不叠合，GetBestRMS 会先把配体叠合到参考上，得到的数值偏小。

2. 02

   重对接不过时按顺序排查

   依次检查：盒子单位与位置；配体和受体的质子化态；增加 exhaustiveness 或换种子；配体环构象；受体结构质量。若晶体构象在打分函数中本身就不是最低点（用 vina --score_only 和 --local_only 对比；这两个模式要求配体坐标已在盒子内，要用带氢的晶体配体准备 PDBQT），继续加大搜索也没用，应换打分函数或改用其他方法。本页 1IEP 实测：晶体构象 --score_only 为 −12.51，--local_only 后 −13.24 kcal/mol，对接最佳构象 −12.8 kcal/mol，三者 RMSD 都小于 1 Å。

3. 03

   阳性对照

   用同一受体、同一盒子、同一参数对接一个已知活性分子（最好有实验 Kd 或 IC50），新分子的分数与它比较。Che 与 Zhang 统计的网络药理学论文中，没有一篇使用阳性对照。

4. 04

   多种子一致性

   3 个种子的第 1 名构象互相之间 RMSD < 2 Å，说明搜索已收敛。仅换种子分数就大幅变化时，先缩小盒子或提高 exhaustiveness。

结果分析

结果分析：选构象、PyMOL 作图与二维互作图

# 在 PyMOL 命令行中执行
load protA_H.pdb, rec
load out/lig-1_s1.sdf, poses
load xtal_lig.sdf, xtal
split_states poses, prefix=pose
delete poses
hide everything
show cartoon, rec
set cartoon_transparency, 0.5
select pocket, byres (rec within 4 of pose0001)
show sticks, pocket and not name N+C+O
show sticks, pose0001 or xtal
color grey70, xtal and elem C
color green, pose0001 and elem C
distance hb, pose0001, pocket, mode=2
label pocket and name CA, "%s%s" % (resn, resi)
orient pose0001
set ray_opaque_background, 0
png pose1.png, width=2400, height=1800, dpi=300, ray=1

# 导出复合物给 PLIP / LigPlot+：配体改成 HETATM 并统一残基名
alter pose0001, resn="LIG"
alter pose0001, chain="L"
alter pose0001, resi="900"
alter pose0001, type="HETATM"
save complex_pose1.pdb, rec or pose0001

# 互作表（txt + xml）和 PyMOL 会话；-p 额外输出图片
plip -f complex_pose1.pdb -t -x -y -o plip_pose1

1. 01

   读懂输出表

   affinity 是 Vina 打分（kcal/mol）。rmsd l.b. 和 rmsd u.b. 是该构象与本次第 1 名构象的距离，与晶体结构无关；u.b. 按原子一一对应、不考虑对称，l.b. 按最近的同元素原子匹配。“对接 RMSD 越小越好”的说法混淆了这两种 RMSD。

2. 02

   选构象

   不只取第 1 名。依次看：第 1 名是否在多个种子中都出现；前几名分数差远小于 Vina 打分误差（2.85 kcal/mol）时，优先选与已知关键互作一致的构象（1IEP 中伊马替尼与铰链区 Met318 主链、守门残基 Thr315 形成氢键）；排除配体大部分暴露在溶剂中或贴着盒子边缘的构象。

3. 03

   导出 SDF

   mk_export.py 依据 PDBQT 头部保存的 SMILES 还原键级、形式电荷和全部氢，导出的 SDF 可以直接给 PyMOL、PLIP、RDKit 和下游分子动力学使用。

4. 04

   三维图

   用下面的 PyMOL 命令画口袋、氢键和晶体构象叠加图。Vina 1.2.x 输出在 PyMOL 中显示不全时，多半是文件中混入了 NUL 字符，改看 SDF。

5. 05

   二维互作图

   PLIP 输出氢键、疏水、π 堆积、盐桥、卤键的残基和距离表；LigPlot+ 和 PoseView（proteins.plus）画二维示意图；Discovery Studio Visualizer 也能画二维图。这些工具都要求输入受体和配体在同一个 PDB 文件里，配体标记为 HETATM。

论文写作

在论文中报告对接结果

审稿人复现对接需要以下信息。缺少盒子中心和尺寸是网络药理学论文中最常见的问题。

• 软件与版本：AutoDock Vina 1.2.7、Meeko 0.8.0、molscrub 0.3.0、PDBFixer 等。
• 受体：PDB ID、链、分辨率；去除了哪些水、离子和配体，保留了哪些金属和辅因子；补了哪些残基。
• 质子化：目标 pH、所用工具（PDBFixer、PROPKA），手动指定的残基状态。
• 配体：来源（PubChem CID 等），质子化与 3D 生成方法，枚举了几个状态。
• 盒子：中心坐标和三维尺寸（Å），以及中心的依据（共晶配体、文献残基、P2Rank）。
• 参数：打分函数、exhaustiveness、num_modes、energy_range、随机种子与重复次数。
• 验证：共晶配体重对接 RMSD；阳性对照分子及其分数。
• 结果：每个分子在各种子中的最佳分数（均值 ± 标准差）、所选构象的依据、互作残基表与三维图。

排错

常见报错对照表

国内常用工具

国内常用工具的实际问题：CB-Dock2、金属蛋白、Discovery Studio 与对接后跑 MD

以下经验来自计算化学公社、Sobereva 博客与中文问答，已与工具文档或原始论文核对。

交给 Agent

把这套流程交给 Agent

指令示例：“用 AutoDock Vina 把这 5 个化合物（附 SMILES）对接到 PDB 1IEP 的 ATP 口袋，先用共晶伊马替尼做重对接验证，再以伊马替尼为阳性对照，每个分子 3 个随机种子，输出分数表、RMSD 表、PyMOL 图和 PLIP 互作表。”

智能体会在工作区安装 Vina 1.2.7、Meeko 和 molscrub（预装环境中已有 RDKit、OpenMM、AmberTools），按本文流程准备受体和配体、重对接并检查 RMSD、批量对接、导出 SDF、作图并整理互作表。交付文件包括 run_redock.sh 等脚本、rec.pdbqt 与 rec.box.txt、各种子的对接输出与日志、redock_rmsd.tsv、分数汇总表、PNG 图和 PLIP 报告。智能体会对结果做对抗审阅，例如检查重对接是否真的达到 2 Å 以内、多种子结果是否一致。

需要你自己核对的是：口袋内残基的质子化态是否符合实验条件，盒子位置是否与文献中的结合位点一致，所选构象的关键互作是否有实验依据，阳性对照是否合适。需要继续做分子动力学验证时，可以在同一工作区里用预装的 GROMACS 接着跑。