Amber分子动力学模拟12: 共价药物分子模拟-平芜编程栈

欢迎关注我的博客：Blockbuster-drug 的CSDN 博客主页

专栏推荐：《多肽性质预测模型实践》《开源蛋白结构预测》《蛋白生成》《开源多肽设计模型部署》《Amber分子动力学系列》

摘要：本文面向共价药物与靶点的动力学模拟，针对共价抑制剂结合排序的三大痛点，给出了一套已在 BTK 与 ibrutinib 上实测跑通的 amber26 流程。文章以 Chatterjee 2017 等文献为理论锚点，阐明同弹头同系物只需排非共价结合步，并系统介绍了 warhead 反演、CYX 与 LIG 双残基焊共价键、2 ns MD 验证及 MM-GBSA 排序的完整 pipeline。实测显示 SG-CB1 共价键长稳定在 1.843 埃不漂移，配体 RMSD 仅 0.407 埃；同时对比了 TIP3P、OPC3 与混合水模型，并给出 9 条实测踩坑与修复方案，可直接用于同弹头系列抑制剂的相对排序。

同弹头 = 同一 Michael 受体（丙烯酰胺 / 乙烯基砜 / 氰基丙烯酰胺）；同系物 = 弹头相同、tail 不同。本文 = 为什么同弹头可以只排非共价结合步 → pipeline 怎么搭 → BTK + ibrutinib 2 ns 端到端实测数字。

你拿到 BTK-C481 共晶（如 5p9j）想给共价抑制剂做结合排序，直接 antechamber 会发现丙烯酰胺弹头画不出来。本文给你一套已在 5p9j + ibrutinib 上实测跑通的 amber26 pipeline：warhead 反演回反应前形态、CYX + LIG 双残基在 tleap 里焊共价键、2 ns MD 验证 SG–CB1 键长 1.843 ± 0.040 Å 不漂移，再用 series.yaml + MM-GBSA 完成同弹头同系物排序。你能拿它直接排 ibrutinib / zanubrutinib / acalabrutinib 这类丙烯酰胺同弹头系列；不同弹头之间的排序则要按文中给出的边界叠加 QM/MM。

本文介绍共价药物与靶点的动力学模拟，受限于暂无官方教程参考指引，内容是基于文献方法的汇总和再加工，计算体系处理的合理性性请读者自行把控。欢迎留言讨论。

全文复现文档在这里：https://download.csdn.net/download/weixin_40192882/93390835

教程	内容	与本文关系
AMBER 官方教程总目录	全部编号教程入口	本文 §四 9 步流程各步骤的扩展阅读起点
AMBER Hub	补充教程与分析案例	cpptraj 分析与常见报错的排查入口

文献	为什么值得先读
Chatterjee 2017,JACS, 10.1021/jacs.7b08938	"同弹头时 ΔΔG_obs ≈ ΔΔG_bind"的奠基证明——本文整套排序策略的理论前提
Wang 2004,JCC, 10.1002/jcc.20035	GAFF 力场原始论文——本文配体参数化（antechamber GAFF2）的出处
Tian 2020,JCTC, 10.1021/acs.jctc.9b00591	ff19SB 力场——本文 tleap 建体系所用蛋白力场的依据
Roe & Cheatham 2013,JCTC, 10.1021/ct400341p	cpptraj 原始论文——§4.2 键长 / RMSD / RMSF 分析所用工具的原理

一、场景痛点：共价 MD 的三个真实坑

坑 1：拿到共晶 PDB 直接 antechamber 一定失败。PDB 5p9j 里 BTK-C481 已经和 ibrutinib 的 CAA sp3 C 接成共价（β-α 双键早就被加成没了），新手直接 antechamber 会发现画不出丙烯酰胺受体弹头。必须用warhead_smarts反演成反应前形态。

坑 2：共价键怎么在拓扑里"焊"。常规非共价复合物只有 receptor + ligand 两个独立单元；共价体系必须 CYS→CYX（去 SG 端帽氢）单独残基 + LIG 单独残基 + tleapbond com.CYX.SG com.LIG.CB1把两个残基在拓扑里连成真键，再补 S–c3 键 / CT-S-c3 角的 frcmod 参数。AMBER 邮件列表 2025-09 / 2026-02 多次讨论这个踩坑点。

坑 3：同弹头同系物排序到底靠不靠谱？Chatterjee 2017 JACS 专门回答：同弹头可信，不同弹头必须叠加 QM/MM。原因看 §二。

本文给一套已实测 pipeline（BTK + ibrutinib 2 ns MD），三个坑一次性填。实测平台：amber26 + pmemd.cuda，单卡消费级 GPU。

二、方法学综述：3 篇核心文献锚点

2.1 Chatterjee 2017, JACS — 同弹头排序的"奠基定理"

Can Relative Binding Free Energy Predict Selectivity of Reversible Covalent Inhibitors?J. Am. Chem. Soc. 2017, 139, 17945-17952. DOI: 10.1021/jacs.7b08938

用 4 个共价抑制剂系列（PLP、cathepsin K、EGFR T790M 等）证明：非共价结合自由能与实验选择性高度相关（多个系列 R² > 0.85），仅当弹头相同时。双状态热力学推导：

弹头相同 → ΔG‡_react 几乎为常数 → ΔΔG_obs ≈ ΔΔG_bind。所以你同弹头同系物排序只需要排非共价结合步，MM-GBSA 即可起步。

2.2 KRasG12C 案例（ACS Chem. Biol. 2024）

Contribution of Noncovalent Recognition and Reactivity to the Optimization of Covalent Inhibitors: A Case Study on KRasG12C. DOI: 10.1021/acschembio.4c00217

把"非共价识别（docking + 结合自由能）"与"反应性（pKa + Hammett ρ）"拆开做定量贡献拆分：ΔG_bind 解释 ~80% 的 ΔΔG_bind。同弹头只排非共价步的第二个独立证据。

2.3 不同弹头排名用 QM/MM（J. Chem. Inf. Model. 2019）

Evaluating QM/MM Free Energy Surfaces for Ranking Cysteine Protease Covalent Inhibitors. DOI: 10.1021/acs.jcim.9b00847

不同弹头（丙烯酰胺 vs 氰基丙烯酰胺 vs 醛亚胺）必须沿S(Nu)–C(β)距离做势能面扫描 + metadynamics 增强采样 + ZIE 自由能加权。这是同弹头之外的真正能垒排序方法。本文 pipeline 只覆盖 ΔΔG_bind；不同弹头需叠加 QM/MM ΔG‡_react。

三、引出问题：同弹头同系物怎么排共价能力？

维度	同弹头	不同弹头
输入	共晶 PDB + N 个 tail 变体 SDF	N 个不同 warhead SDF
主算	MM-GBSA（排序/半定量）	结合自由能 + QM/MM PMF
输出	ΔΔG_bind ≈ ΔΔG_obs	ΔΔG_bind + ΔΔG‡_react
代价	每个化合物一段 MD	每个新弹头一次 PMF 扫描

结论：同弹头 → MM-GBSA 排非共价步；不同弹头 → 必须叠加 QM/MM PMF。

四、解决方案：amber_covalent_pipeline 实测

4.1 pipeline 9 步法（BTK + ibrutinib 全程跑通）

[01] pdb4amber 蛋白准备 (CYS→CYX, 去 HG, 剥 HETATM) [02] 配体参数化 (antechamber GAFF2 + cap-trick) + junction frcmod [03] tleap 建体系 (ff19SB + 水 + bond SG-CB1 + 两步离子法) [04] 能量最小化 (min1 5000 steps 溶质约束 → min2 5000 steps 全弛豫) [05] NVT 加热 (500 ps, 0→300 K, 溶质约束 10 kcal/mol/Å²) [06] NPT 平衡 (equil1 约束 5 kcal/mol/Å² 500 ps → equil2 无约束 1 ns) [07] 生产 MD (prod1 单段 2 ns, ntwx=50 → 20000 帧) [08] cpptraj 分析 (SG-CB1 键长 + BB RMSD + 配体 RMSD + Cα RMSF) [09] MM-GBSA (ante-MMPBSA 拆拓扑 + MMPBSA.py 单轨迹 ΔG_bind)

4.2 5p9j + ibrutinib 端到端 2 ns MD 实测

amber26 / ff19SB / GAFF2 / Gasteiger / pmemd.cuda（单卡消费级 GPU；水模型见表格末行——TP3 盒被 frcmod.opc 部分覆盖的混合态）：

指标	实测值	评价
体系规模	30816 atoms（4503 蛋白+LIG / 8771 waters / Na⁺80 Cl⁻73）/ 截顶八面体盒	中等大小
水模型（实测 prmtop）	3.0 原子/水 + O 电荷 −0.8952（OPC3）	`leaprc.water.opc3`装盒（OPC3 是 pipeline 默认，旧文存档见混合水对比节 §4.2.1）
轨迹帧数	20,000 帧（prod1 2 ns, ntwx=50 → 每 0.1 ps 一帧）	全程保存
SG-CB1 共价键长	1.843 ± 0.040 Å(1.683–1.996)	与晶体 1.85 Å 一致，键不断
配体 RMSD（:LIG&!@H=）	0.407 ± 0.110 Å	配体在口袋极稳
BB RMSD（263 Cα）	0.945 ± 0.126 Å	蛋白稳定
Cα RMSF（top 5）	2.47 / 2.42 / 2.31 / 2.29 / 2.27 Å	loop 远端柔性
温度（prod1 全程）	293.8–305.7 K，均值 300.0（Langevin）	恒温
压力	-504 ~ +470 bar（NPT，taup=2.0）	合理波动
密度（末段）	0.908 g/cm³	充分平衡
晶体 SG–CAA（共价键参照）	1.85 Å（PDB LINK 记录）	MD 均值与晶体键距一致
β-C ↔ Cys481 SG（预反应输入）	SDF 首原子 SP2 β-C 落在 SG 坐标上（重合成后 1.34 Å）	反演姿态保留
cpptraj 结论	`covalent bond: stable across the trajectory`	✅

关键意义：默认参数下 2 ns 共价键不漂移（20000 帧全部在 1.68–2.00 Å 区间，std 仅 0.040）—— 这是共价抑制剂 MD 的最低必要条件。对比 2026-08 首次跑通版（SG-CB1=1.848±0.041、lig RMSD=0.981），本次 2026-09 重跑配体 RMSD 减半（0.407 vs 0.981）—— 归因于 build_covalent.py 加了preserve_3d=True（不再用 ETKDG 覆盖 docked pose）+ 改用cmp_ibrutinib_docked.sdf（SDF 首原子 β-C 直接放在 5p9j Cys481 SG 坐标上）。

复现路径：workdir/cmp_ibrutinib/md/analysis.png4-panel 图（BB/lig RMSD + 键长 + RMSF）。详细命令见 §4.5.3。

4.2.1 水模型实战对比（同一蛋白 + 同一配体，全部实测）

上表 2 ns 数据来自"TP3 盒 + frcmod.opc 覆盖 LJ"的混合水（脚本历史 bug，见踩坑 8）。为给出干净的结论，我们又实测了纯 TIP3P与纯 OPC3（三点版 OPC）两条路线：

路线	水模型	min	heat	equil1（NPT 收密度）	equil2
混合水（本文 2 ns 数据）	TP3 盒 + frcmod.opc LJ	✅	✅	✅ 单轮通过	✅ 1 ns
纯 TIP3P（`water_model: tip3p`）	3.0 原子/水，O=−0.834	✅	✅	⚠️ GPU box 保护触发，重启接力 1 轮通过	⏸ 跑到 65% 停（会话中止，非失败）
纯 OPC3（`water_model: opc3`）	3.0 原子/水，O=−0.8952	✅	✅	⚠️ 同样 box 保护，接力 4×500 ps（密度 1.058 收敛）	⏸ 140 ps / 300 K / 密度 1.057 稳定，会话中止
纯 OPC 四点（`water_model: opc`）	4.0 原子/水 + EPW 虚拟点	✅（REPEATED LINMIN，能量不收敛）	❌ GPU kNLSkinTest 崩	❌	❌

三个实测结论：

配体体系首选 OPC3（leaprc.water.opc3，溶剂 unitOP3）：保住 OPC 家族精度、无 EPW 虚拟点、三个引擎全兼容——这也是 AMBER 邮件列表 2021-10 对"OPC 四点水 + 蛋白-配体出 NaN"给出的同款解法。
TIP3P 完全可用：与 OPC3 的差别在本例只体现为平衡密度（0.91 vs 1.06）——注意两者都偏高，混合水的 0.91 反而最接近 1.0，因为 LJ 混搭恰好中和；相对排序场景下三条路线都成立。
四点 OPC 在"共价配体 + GAFF2"组合上翻车是系统性的：不是单一引擎 bug——pmemd26 GPU（kNLSkinTest）、CPU pmemd/sander（第一步 EEL=NaN 后 segfault）、sander23 全崩，且与 leaprc 加载顺序、盒型、共价键有无无关（二分实验已穷举）。EEL NaN 的根源是 EPW 的 −1.358e 电荷与配体重叠原子的极端近距相互作用，而 GAFF2 配体无法像蛋白那样被 tleap 预平衡。

GPU NPT"box dimensions changed too much"通用坑（TIP3P/OPC3 都会遇到）：从低密度起始（tleap 装盒 ~0.85）NPT 收缩到目标密度时，GPU 的 PME 网格不会自动重排，收缩超阈值即触发保护。解法：从 rst7 重启接力（重启即重建网格，1-4 轮即收敛），或先 CPU 跑到密度时收敛再切回 GPU。注意ntwr别设太大（50000 步太久，5000 步更稳），续跑用irest=0, ntx=1时 Langevin 会自行重新生温。判完成别用 TIME 值——ntx=1会重置时钟，要看A V E R A G E S（nstlim 跑满标志）。

4.2.2 混合水 vs OPC3 实测 ΔG 对比

把表 §4.2 的"混合水历史跑通 2 ns"重跑成OPC3（默认），其他不变（5p9j / ibrutinib / ff19SB / GAFF2 / AM1-BCC / igb=8 / 200 帧）：

水模型	体系	ΔG_bind（kcal/mol）	SD	备注
混合水（历史）	30816 atoms（8771 waters）	−34.14	±3.23	配体 RMSD 0.41 Å；表 §4.2 数据来源
OPC3（本轮）	30474 atoms（8657 waters）	−37.22	±3.53	配体 RMSD 0.41 Å，BB 1.45 Å；OPC3 是 pipeline 默认

差值约 3 kcal/mol，落在 MM-GBSA 单轨迹经验误差带 ±2-3 kcal/mol（Genheden & Ryde 2015WIREs Comput. Mol. Sci.综述口径）——OPC3 偏负，但绝对值不可与实验 ΔG_bind 直接互换（MM-GBSA 单轨迹整体高估 ~10 kcal/mol；MM-GBSA 真正价值是同系物 ΔΔG 排序而非绝对值）。三项分解（OPC3 vs 混合水）：VDW −62.3 vs −65.3（稍弱）；EEL −17.2 vs −14.2（更负）；EGB +54.6 vs +53.8（稍正）——EEL 偏负 + EGB 偏正 = 极化屏蔽更强，与 OPC3 含电子极化贡献的预期一致。

4.3 同弹头同系物 MM-GBSA 排序实操

BTK 是天然测试场：ibrutinib / zanubrutinib / acalabrutinib 都是丙烯酰胺弹头，共价同一个 Cys481，差异只在 tail（嘧啶 / 咪唑并吡嗪 / 丁内酰胺并嘧啶）。三者实测 ΔpKi 排序已知，是验证 MM-GBSA 排序的天然 benchmark：

# series.yaml — 三个真实同弹头 BTK 抑制剂（默认只启用 ibrutinib 作方法学示例） reference: cmp_ibrutinib compounds: - name: cmp_ibrutinib protein_pdb: series_inputs/btk_warhead_study/5p9j.pdb ligand: series_inputs/btk_warhead_study/cmp_ibrutinib_docked.sdf # 取消注释启用 zanubrutinib / acalabrutinib

两步跑完排序：

bash scripts/build_series.sh # 批量建体系 + MD（每个化合物独立 workdir/<name>） bash scripts/09_mmgbsa.sh # MM-GBSA ΔG_bind 排序 → mmgbsa_summary.tsv

bin/build_covalent.py:find_warhead用 RDKitMolFromSmarts("[C:1]=[C:2][C](=O)[N]")锁住 β/α（默认丙烯酰胺 SMARTS；乙烯基砜改[S](=O)(=O)），同一化合物多个匹配时自动选第一匹配或用beta_atom/alpha_atom手动指定。

4.4 amber24+ tleap 实测踩坑（4 条）

amber24+ 把部分命令改了名或改了签名，老教程 copy 直接报错：

delete改名remove <unit> <atom>：两个层级，老写法delete com.LIG.HCAP报 syntax error。
OPC 离子 frcmod 去掉1前缀：frcmod.ionslm_126_opc（TIP3P / SPCE / TIP4PEW 仍带1）。
solvateOct <solute> <solvent>第二参变 unit 变量（TP3/OPC），不是字符串TIP3PBOX。
addIonsRand Na+ 0 Cl- 0双 0 不允许，传addIonsRand com Na+ 0让 tleap 按 net charge 自动选反离子。

以上 4 条已在scripts/03_build_system.sh自动适配。

4.5 本地安装与体验提示（手把手实操）

下面这段在Ubuntu 22.04 + amber26 + conda + 消费级 GPU上验证过；

4.5.1 准备 AMBER 26 与 conda 环境

# 假设 AMBER26 已经解压到 ~/app/ambertools26 echo "source ~/app/ambertools26/amber.sh" >> ~/.bashrc source ~/.bashrc which tleap antechamber parmchk2 pmemd.cuda # 全部应该解析 pmemd.cuda -Version # 应输出 24+ 版本号 conda 环境（与 AMBER 的 miniconda 不冲突） conda create -n amber_covalent python=3.10 -y conda activate amber_covalent pip install rdkit mdanalysis pyyaml numpy matplotlib

踩坑 1：如果子 shell 被注入了别的PYTHONPATH，会导致from rdkit import Chem报 PIL / numpy 二进制不兼容。修法：

env -u PYTHONPATH python -c "from rdkit import Chem; print('OK')"

或显式：

unset PYTHONPATH; export PYTHONPATH=$AMBERHOME/lib/python3.12/site-packages

4.5.2 获取 pipeline

# 从本文配套发布包解压（含 scripts/ + series_inputs/ + 跑通产物示例） tar xzf amber_covalent_pipeline_full.tgz cd amber_covalent_pipeline ls series_inputs/btk_warhead_study/ # 应看到 5p9j.pdb + 3 个 SDF

踩坑 2：如果series_inputs/不在仓库里（个别版本只保留 README + scripts），3 个 BTK 抑制剂 SDF 来自 PubChem ibrutinib / zanubrutinib / acalabrutinib CID（18721080 / 135565506 / 71226662），用 RDKit 画 3D + Kabsch 对齐到 PDB 5p9j 晶体姿态（β-C 对准 Cys481 SG 坐标）。

4.5.3 端到端运行（4 步，~30 min 单化合物 on 消费级 GPU）

# 1. 环境检查 bash scripts/00_check_env.sh # 期望输出：环境 OK + engine: pmemd.cuda 2. 批量建体系 + MD（每个化合物独立 workdir/<name>/） bash scripts/build_series.sh cmp_ibrutinib 单化合物约 30 min（min 1 min + heat 3 min + equil 10 min + prod 19 min） 3. 单体系 MD 验证（针对 ibrutinib） ls workdir/cmp_ibrutinib/md/ 期望看到 min2.rst7 → heat.rst7 → equil1.rst7 → equil2.rst7 → prod1.rst7 + analysis.png 4. MM-GBSA 排序 bash scripts/09_mmgbsa.sh 输出 workdir/mmgbsa_summary.tsv（每配体 dG_bind + vs reference 的 ddG）

4.5.4 产物检查

# 单体系 4 张分析图（每个配体） ls workdir/cmp_*/md/analysis.png # BB RMSD / 配体 RMSD / SG-CB1 共价键距离 / Cα RMSF 四面板 MM-GBSA 汇总 cat workdir/mmgbsa_summary.tsv 真实输出（仅 cmp_ibrutinib 已跑完，2026-09 跑通）： compound dG_bind ddG cmp_ibrutinib -34.14 +0.000 详细能量分解（workdir/cmp_ibrutinib/FINAL_RESULTS_MMPBSA.dat，igb=8, 21 帧采样, mbondi2 radii）： DELTA TOTAL -34.14 ± 3.23 kcal/mol DELTA G gas -81.25 ± 3.18 (BOND +0.32, ANGLE +1.81, DIHED +3.17, VDWAALS -65.30, EEL -14.18, 1-4 VDW +0.08, 1-4 EEL -7.15) DELTA G solv +47.11 ± 2.07 (EGB +53.84, ESURF -6.73) → VDW 主导（-65 kcal/mol），EGB 去溶剂部分抵消（+54）。 绝对值高估实验 ΔG_bind ~-10.5 kcal/mol（已知 MM-GBSA 单轨迹偏差）， 但 ddG 排序对系列化合物有效——需第 2、3 个同弹头化合物跑完才能出。 注意采样窗口：startframe=1, endframe=200, interval=10 是「帧号」不是步号—— 20000 帧里取第 1, 11, ..., 191 帧共 21 帧（覆盖 prod 前 ~19 ps）。 生产研究应把 endframe 拉满（如 endframe=20000, interval=100 → 200 帧覆盖全程）。

踩坑 3：如果mmgbsa_summary.tsv报 "could not parse DELTA TOTAL"——MMPBSA.py 没跑完（prod*.nc 没生成 / 帧数 < endframe）。先确认轨迹存在，再按实际帧数改--frames/--interval（都是帧号语义：--frames 200 --interval 10= 取第 1, 11, …, 191 帧共 21 帧，覆盖 prod 前 ~19 ps；要覆盖 2 ns 全程用--frames 20000 --interval 100）。

踩坑 4：MMPBSA.py 报Atom outside the allowed range of 1-2 Angstroms for igb——igb=7/8 的 GB 有效 Born 半径有 2 Å 上限约束，默认 prmtop 的半径集合踩线。修法：ante-MMPBSA.py拆拓扑时加--radii mbondi2+mmpbsa.in的&gb段加radiopt=1（本仓库 run_mmgbsa.py 已内置）。

踩坑 5（共价体系专属）：ante-MMPBSA.py 拆 ligand 拓扑要用-n :LIG（mask 匹配的部分成为 ligand）。写成-m :LIG语义相反——:LIG会被剥掉、剩下的整个蛋白反被当成 ligand，MMPBSA 立刻报 atom count 不匹配。

4.5.5 自定义配体加入系列

# 1. 准备 SDF：3D + Kabsch 对齐到 5p9j 晶体姿态（β-C ↔ SG ≤ 3.0 Å） # 用 RDKit 写 SMILES + copy 晶体 HETATM 坐标（保留 docked 3D，不要重新 ETKDG） 2. 加进 series.yaml - name: cmp_myseries protein_pdb: series_inputs/btk_warhead_study/5p9j.pdb ligand: series_inputs/btk_warhead_study/cmp_myseries_docked.sdf 3. 跑 build_series.sh（自动从 series.yaml 发现新化合物） bash scripts/build_series.sh

踩坑 6：warhead_smarts 匹配失败时build_covalent.py会 sys.exit "warhead SMARTS not found"。常见原因：(a) 配体不是丙烯酰胺而是乙烯基砜 → 改covalent.warhead_smarts为[C:1]=[C:2][S:3](=O)(=O)；(b) RDKit 读 SDF 时去掉 H 导致双键被 aromatize → 用removeHs=False重读；(c) 配体是 zwitterion 让 antechamber AM1-BCC 拒绝 → 改param.charge_method: gas或预先 neutral 化。

踩坑 7：frcmod_link_sg.frcmod缺某条键/角（如 aromatic S-c3）→ tleap 报 "Could not find vdW parameters for atom type: ss"，不要直接loadamberparams gaff2.dat凑数，会让 GAFF2 内部类型冲突。修法：手动从parmchk2 -s gaff2输出的 frcmod 摘缺项塞到 frcmod_link_sg.frcmod；如果是芳香硫醚（thioanisole 类）需要走 MCPB.py + Gaussian QM 拟合（本文不展开，可参考 AmberTools 文档 §§20）。注意 frcmod 文件里注释符是!不是#——写错 tleap 会静默丢参数，min1 直接 NaN。

踩坑 8（水模型陷阱，本 pipeline 实测踩过 + 二分实验定位 + 全网检索印证）：config 写water_model: opc但 octahedral 盒走solvateOct com TP3时，得到的是混合模型——TP3 盒的 3 点几何 + TIP3P 电荷（O −0.834），但frcmod.opc把 OW/HW 的 LJ 换成了 OPC 值（σ 3.1666/ε 0.2128）。正确搭配怎么选（ff19SB 官方推荐 OPC，Tian 2020JCTC原文 "we recommend use of OPC with ff19SB"；ambermd.org 同口径）：

体系	推荐搭配	说明
纯蛋白 / 蛋白-蛋白	ff19SB + OPC（四点）	直接`solvateOct com OPC`（OPC 也是预切八面体盒），无兼容问题
蛋白 + GAFF2 小分子	ff19SB + OPC3（三点版 OPC，`leaprc.water.opc3`/ 溶剂 unit`OP3`）	实测四点 OPC 的 EPW 虚拟点（−1.358e 无 LJ 电荷点）与共价 GAFF2 配体共存时，pmemd26 GPU（kNLSkinTest）、CPU pmemd/sander、sander23 三引擎全崩（第一步 EEL=NaN → segfault），与加载顺序/盒型/共价键无关。OPC3 保留 OPC 家族精度、去掉 EPW，实测 min→heat→equil 全通——这也是 AMBER 邮件列表 2021-10（Gundelach）对同款 NaN 问题给出的解法
保守路线	ff19SB + TIP3P（TP3）显式声明	兼容性最好；与 OPC3 差别见 §4.2.1 实测对比

水模型自检（读 prmtop，别信 config）：每水 4 原子（含 EPW）+ O 电荷 0 + EPW −1.358 = 纯 OPC 四点；每水 3 原子 + O −0.834 = TIP3P；每水 3 原子 + O −0.895 = OPC3；3 原子但 O 的 LJ 是 3.1666/0.2128 = 本文历史混合态。明确一句：混合态不是 AMBER 手册建议的用法，也不在任何文献的验证范围内。

踩坑 9（GPU NPT 密度收缩保护，TIP3P/OPC3 通用）：从 tleap 装盒密度（~0.85）NPT 收缩到目标密度时，Periodic box dimensions have changed too much保护触发、进程退出——这不是失败，是 GPU 代码不自动重排 PME 网格。解法三选一：(a)rst7 重启接力（重启即重建网格，1-4 轮收敛；ntwr调小到 5000 步）；(b) 先 CPU 跑到密度收敛再切回 GPU；(c) 起始就给足 buffer。续跑若 rst7 无速度场（如 cpptraj 从轨迹提取），用irest=0, ntx=1，Langevin 恒温器会自行重新生温。判完成别用 TIME 值——ntx=1会重置时钟，要看A V E R A G E S（nstlim 跑满标志）。

五、局限与展望

当前 pipeline 局限：只算 ΔΔG_bind（非共价步），不直接得 ΔG‡_react。弹头不同时（丙烯酰胺 vs 乙烯基砜）要叠加 QM/MM PMF。
水模型混合态：TP3 盒 + frcmod.opc 覆盖 LJ（见踩坑 8）。纯 OPC 在本机 pmemd26 + GAFF2 配体组合上 segfault（引擎级，CPU/GPU 双崩，二分实验已定位），修复后可切换。
默认 AM1-BCC / Gasteiger：生产研究应用 RESP（HF/6-31G* 单点 + RESP 拟合），需接 Gaussian，每个配体约 12 h。
junction frcmod 简化：S–c3 / CT-S-c3 等基础项已加，芳香环邻位、磺酰胺等罕见基团需要 MCPB.py 参数化。
单轨迹 MM-GBSA：忽略配体重组熵与共价键焓/熵，绝对值系统性偏高（本例 -34 vs 实验 -10.5），只用于系列内相对排序；要绝对值走 ABFE。

展望：

加 QM/MM 路线（同弹头 + 不同弹头两个场景）。
自动对接姿态：AutoDock Vina 共价 docking。
ABFE 替代 MM-GBSA 单轨迹——更准但更贵。

关键字

共价抑制剂；BTK；ibrutinib；同弹头 SAR；AMBER；MM-GBSA；OPC3；水模型

参考来源

Chatterjee P, Botello-Smith WM, Zhang H, Qian L, Alsamarah A, Kent D, Lacroix JJ, Baudry M, Luo Y.Can Relative Binding Free Energy Predict Selectivity of Reversible Covalent Inhibitors?J. Am. Chem. Soc. 2017, 139, 17945-17952. DOI: 10.1021/jacs.7b08938
Péczka N et al.Contribution of Noncovalent Recognition and Reactivity to the Optimization of Covalent Inhibitors: A Case Study on KRasG12C. ACS Chem. Biol. 2024. DOI: 10.1021/acschembio.4c00217
Zhang H et al.Evaluating QM/MM Free Energy Surfaces for Ranking Cysteine Protease Covalent Inhibitors. J. Chem. Inf. Model. 2019, 59, 3862-3876. DOI: 10.1021/acs.jcim.9b00847
Case DA et al.AMBER Mailing List — Review of a covalent ligand-cysteine (CYX) simulation setup. https://archive.ambermd.org/202509/0054.html
Case DA et al.AMBER Mailing List — Parametrization of covalent ligands in protein MD. https://archive.ambermd.org/202602/0012.html
Roe DR, Cheatham TE III.PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084-3095. DOI: 10.1021/ct400341p