分子对接自由能计算

分子对接做了这么多年，最大的感受是：结合自由能不是算出来的，是选出来的。同样的体系，用不同的对接软件、不同的评分函数、不同的采样参数，跑出来的结果可能完全不一样。

所以对接的核心不是选一个”看起来最对的”结果，而是建立一套稳定的打分体系，然后批量筛选。MM/GBSA是这套体系中性价比最高的方法之一——比量子力学精确，比纯经验打分可靠。

PDB下载的晶体结构往往带有多余的配体和水分子。处理流程：去掉水分子(距离活性位点8Å以外的水可以全删);加氢(PDB2PQR或H++自动处理);检查残基质子化状态(特别是His、Asp、Glu);保存为PDBQT格式。

配体分子的处理更繁琐：构象优化(PM7或PM6半经验方法);电荷计算(Gasteiger或AM1-BCC);可旋转键定义(Vina自动处理，但自定义时注意关键药效基团不要设成可旋转)。

Vina是目前最常用的对接工具，速度快，精度也不错。参数设置：exhaustiveness = 32(搜索穷尽度，越高越全面);num_modes = 20(输出构象数);energy_range = 3(与最优构象的能量差上限)。

搜索空间(box)要刚好包住活性口袋，四周边界留5Å左右。太大了搜索效率低，太小了漏掉真实结合模式。

Vina打分(kcal/mol)只是初筛依据，不能直接当结合自由能用。一般流程是：先Vina批量对接，筛选top 20%构象，再用MM/GBSA或MM/PBSA精打。判断结合模式是否合理，看两点：关键氢键是否形成(尤其是配体与蛋白主链NH/CO的氢键);疏水基团是否嵌入疏水口袋。

MM/GBSA(分子力学/广义玻恩表面积)是精打分的标准方法。计算公式：ΔG_bind = ΔE_MM + ΔG_sol + ΔG_SA。ΔE_MM是分子力学能(键合+静电+范德华);ΔG_sol是溶剂化自由能(极性部分用GB模型，非极性部分用表面积);ΔG_SA是脱溶剂化自由能。

用AmberTools的MM_PBSA.py：

python $AMBERHOME/bin/MM_PBSA.py -i mmgbsa.in -o results.dat -f vina_results.csv

输入文件里需要指定：igb=2(GB模型，ICM方案);saltcon=0.15(离子浓度)。

看ΔG_bind的绝对值意义不大，关键是相对排序。同一系列化合物，ΔG_bind越负，结合能力越强。自由能分解(per-residue)能告诉你每个氨基酸残基对结合的贡献：哪几个残基是主要驱动力，哪几个是障碍。

对接分数好的构象，但氢键位置不对，检查受体活性位点的质子化状态。His、Asp、Glu的质子化状态对静电势影响极大，质子化状态错了，氢键预测必然不准。MM/GBSA算出来的值偏正或偏负太多，检查介电常数设置——igb=2的介电常数默认是1.新体系可能需要调到2来更好地描述静电屏蔽。

回过头看，分子对接自由能计算的本质是在海量候选构象中，用分级筛选的方式找到最有价值的那些。Vina负责广度，MM/GBSA负责精度——两者配合，才是性价比最优的策略。