3POZ / 配体 03P:一次真实短分子动力学与轨迹分析

module-07 报告 · TeachOpenCADD T020 案例


0. 摘要

用输入的 3POZ / 配体 03P(上游由 T019 准备)体系,在本机 CPU 上实跑了一段 真实分子动力学,然后把新轨迹与输入包自带的上游轨迹用完全相同的口径 重新分析、对照,出 A–F 六组图。

⚠️ 一句话边界:"跑完了、看着稳" ≠ "结合"、"有效"或"收敛"。 本报告只交付可核验的坐标与计数,以及它们能和不能支持的说法。


1. 做了什么(以及为什么起点不是晶体结构)

输入包提供的是已经准备好的溶剂化体系坐标(topology.pdb:蛋白 + 03P + TIP3P 水 + Na⁺/Cl⁻),不是原始晶体文件。这带来两个必须交代的事实:

  1. 上游坐标被复用为起点 —— 新模拟从上游客观给定的同一套准备坐标出发, 而不是从 3POZ 重新做一遍 PDBFixer 加水加离子。
  2. 上游的 System XML 没有随包提供 —— PDB 只存坐标、不存力场。 所以必须重新建立 OpenMM System,也就是重新生成配体力场参数。

因此:力场"名字"一致,参数集并不逐比特一致。 这条差异是本报告最关键的诚实声明, 它决定了后面所有跨轨迹比较的措辞(见 §5、§6)。


2. 关键选择(写给初学者)

2.1 配体的化学态必须先"认下来",不能靠坐标猜

PDB 坐标不包含键级和质子化状态。一个只有坐标的小分子无法直接参数化。 本目录的做法是分两步、并留下证据:

核验项 结果
组分名 / 别名 N-{2-[4-({3-chloro-4-[3-(trifluoromethyl)phenoxy]phenyl}amino)-5H-pyrrolo[3,2-d]pyrimidin-5-yl]ethyl}-3-hydroxy-3-methylbutanamide / TAK-285
分子式 CCD C26 H25 Cl F3 N5 O3 ↔ 感知 C26H25ClF3N5O3 ✅
InChIKey ZYQXEVJIFYIBHZ-UHFFFAOYSA-N 完全一致 ✅
原子映射 63/63(按连接性图同构,非按原子名)✅
元素 / 形式电荷 63/63 / 63/63 ✅
键级分歧 0 处真实分歧(66 键中 54 精确一致、12 为同一芳环的 Kekulé 等价写法)✅
立体中心 无(与 UHFFFAOYSA 无立体层一致)

注意一个容易踩的坑:准备坐标里原子名被重新编号了(N1…N5 / C1…C26 / F1…F3…), 与 CCD 和晶体副本的原始命名(N/N1/N2/…、C/C1/…、F/F1/F2、CL)不同。 所以映射必须用图同构,用原子名匹配会得到 47/63 的假失败。

结论:配体身份与键级已被独立证据确证,不是"我说它是"。

2.2 力场与积分:贴着原教程,但把"改动"写出来

项目 原章节(T019) 本次运行 关系
力场 amber14-all + amber14/tip3pfb 同 一致
配体力场 GAFF 2.2.20 + AM1-BCC 同 一致(参数重新生成)
非键 PME,1.0 nm 同 一致
积分器 Langevin,300 K,摩擦 1/ps 同 一致
步长 2 fs 同 一致(标准氢质量,未用 HMR)
恒压 无 → NVT 无 → NVT 一致
溶质 X–H 约束 未显式设置 HBonds 改动,已披露
平台 (原文为 GPU/Colab) CPU 12 线程 本机无 GPU

两点需要解释清楚:

2.3 为什么可以用 CPU

本机无 GPU。技能要求:CPU 做小规模烟雾/短模拟是合法的, 但绝不能把 CPU 回退伪装成 GPU 生产运行。因此本报告在 A1、B 组图和执行记录里 显式标注 CPU 12 线程与实测 ns/day,并给出实测时间预算:

线程数 ns/day
4 1.25
8 1.47
12 1.56(实测最好)
16 1.48

这是按证据选的线程数,不是靠假设"线程越多越快"。

2.4 一个真实踩坑:AM1-BCC 首次尝试被中断与单线程重试(完整记录)

AM1-BCC 由 AmberTools 的 sqm 做半经验 QM,是本流程最贵的一步。

细节见 ../outputs/new_md/charge_assignment_history.md。电荷结果已缓存到 ../inputs/ligand_03p_charges.json(键 = SDF sha256 + 方法 + SMILES), 所以后续任何无关的小错误都不必重跑 QM。


3. 执行记录(全部来自引擎自身)

项目 数值
体系 58,047 原子 · 317 蛋白残基(701–1017)· 17,587 水 · 55 Na⁺ · 48 Cl⁻ · 1 × 03P
盒子 8.489 nm 立方,NVT 全程不变
约束数 55,360(HBonds)
自由度 受约束 118,778(=3N−3−约束数);朴素 3N = 174,141
最小化 投影 +188.9 / 最小化 −26.5 kJ/mol;原子最大位移 0.0033 / 0.0042 nm
平衡 限制性 NVT 5 ps(2,584 个重原子受限,T = 297.8 K)+ 自由 NVT 3 ps(302.7 K)
正式采样 9,000 步 × 2 fs = 18.00 ps,36 帧,每 0.5 ps 一帧
温度(正式采样) 均值 299.93 K,SD 1.29 K(n = 36 帧,ddof = 1,与 ../tables/summary.json 的 temperature_sd_K = 1.2935 一致)
势能(正式采样) 均值 −932,160 kJ/mol,SD 988(0.11%);全程未出现能量发散(有限值)。此 SD 仅为描述量,不作为平衡的充分证据
墙钟 总墙钟 ≈33.6 分钟(started_utc→finished_utc = 2,014 s);所列计算阶段合计 ≈33.3 分钟(2,000 s)= 最小化 19 s + 限制平衡 370 s + 自由平衡 233 s + 生产 1,378 s;生产 ≈23 分钟
吞吐 1.13 ns/day(12 CPU 线程)
平台 OpenMM 8.6.0.dev · CPU

时间的三条独立来源互相吻合,不是照抄某个文件名:

  1. 执行设置:integrator.getStepSize() 实测 2.0 fs × 9,000 步 = 18.00 ps;
  2. 引擎状态日志 prod.log:末行 time = 18.000 ps;
  3. DCD 头部的步长信息:ISTART=250、NSAVC=250、delta=0.040910 AKMA = 0.002 ps → 首帧 0.50 ps、帧间隔 0.50 ps,与 run_report.json 记录的上报间隔一致 (脚本会比对,不一致就报错退出)。
  4. 上游 XTC 头自带真实时间:10–1000 ps。

关于 DCD 的正确说法:DCD 没有每帧时间数组,但带有步长/步数信息, 所以物理时间是可以恢复的,不能说成"完全没有物理时间"。 另外头部 NSET 帧数会随写入更新,对已写完的文件它是正确的, 只是不能在读一个仍在写入的文件时把它当作最终帧数。

3.1 限制力不是摆设:跑了引擎级测试

限制性平衡用 periodicdistance(...) 表达式的谐振子。这个 OpenMM 构建版的 CustomExternalForce 没有 setUsesPeriodicBoundaryConditions 设值器 (只有取值器),所以最小镜像只能靠表达式函数。这是关于引擎行为的断言, 因此写了 5 项测试(../outputs/new_md/restraint_test.json,全过):

测试 验证内容 结果
T1 参考坐标处能量与受力为 0 ✅ 0 / 0
T2 已知位移 0.1 nm → 能量 5.000、回复力 −100.0(解析值 5 / −100) ✅
T3 整体平移一个盒矢量后不变 ✅ 差 1e−14
T4 跨盒面时按最小镜像(0.1 nm)而非直线穿越(2.9 nm) ✅ 周期式 5.0 vs 非周期式 4205.0
T5 真实体系 2,584 个受限原子在自身参考坐标处为 0 ✅

T4 还量化了写错的代价:非周期式的能量放大 841×、力放大 29× —— 所以这不是形式主义。

3.2 周期边界重建:用可证伪的判据,而不是"看起来还行"

两条轨迹的键长都做过检查:5,176 根蛋白键,最大 0.1932 nm, 直线距离与最小镜像距离之差 ≈ 1.6–1.8e−8 nm,超过 0.01 nm 的键 = 0。 也就是说蛋白在两套文件里本来就没断,重建是恒等操作(坐标最大改动 1.5e−8 nm), 下游用周期距离计算是安全的。

这里也纠正一个我先前用错的判据:包围盒对角线不是"是否断裂"的判据。 一个完整的 317 残基激酶域本身张成 6.61 × 5.43 × 6.72 nm,对角线 10.59 nm > 盒边 8.489 nm, 但它是完整的。正确判据是逐轴坐标跨度 vs 盒边(6.61/5.43/6.72 < 8.489 全过), 加上逐键直线距离 vs 最小镜像距离。


4. 图组:每条曲线在回答什么

全部图为 PNG + SVG + PDF 同步导出;每张图旁的 CSV 只是该图所画面板的便捷副本, 完整数据一律在 ../tables/,逐面板对应关系见 ../tables/figure_index.md。

A 组 — 到底模拟了什么

A

B 组 — 引擎自身诊断

B

C 组 — 蛋白行为

C

D 组 — 配体在口袋里的状态

D

D2 与 D1 是两种不同的拟合(自拟合 vs 受体拟合), 本报告没有做过把位移"分解为内部形变 + 整体漂移"的计算, 因此只分别描述各指标,不对两者之差做归因。

E 组 — 口袋接触

E

轨迹 氢键 方向 占据率 平均 H…A 平均角
新跑 03P:O3 → Asp855:OD2 配体→蛋白 1.00 0.188 nm 158°
新跑 Met793:N → 03P:N4 蛋白→配体 0.75 0.239 nm 155°
新跑 Lys745:NZ → 03P:O3 蛋白→配体 0.28 0.283 nm 154°
上游 Met793:N → 03P:N4 蛋白→配体 0.95 0.218 nm 160°
上游 03P:O3 → Ser720:O 配体→蛋白 0.49 0.343 nm 114°

两条轨迹都保持 Met793 铰链氢键(0.75 / 0.95),这是 3POZ 系列体系里预期的 铰链结合特征。新跑在这段窗口内额外持续满足几何氢键判据的是 03P:O3 → Asp855:OD2(36/36 帧);上游则出现一个仅 0.49 占据率的 Ser720 接触。

这里不应称其为"盐桥":03P 的 O3 是中性羟基氧(形式电荷 0,已验证), 而盐桥需要两个带相反电荷的基团。上表按几何判据给出的是候选氢键, 不区分静电类型,也不构成对相互作用的化学分类。

F 组 — 短模拟不能说明什么

F


5. 两条轨迹各自说明了什么

上游轨迹(10–1000 ps,100 帧)说明: 在一个更长的窗口里,蛋白保持折叠、Met793 铰链氢键持续存在, 03P 在该窗口内未见脱离迹象;配体内部形变与净位移都处于亚埃到约 1.3 Å 量级。 它给出的是这段 1 ns 内的持续性,同样不构成亲和力或收敛结论。

新跑轨迹(0.5–18 ps,36 帧)说明: 在同一套上游准备坐标上、用重新生成的 GAFF/AM1-BCC 参数重跑, 得到的是在物理参数与采样都与上游不同的条件下,短窗口内再现了同样的定性图像: Met793 氢键保留,配体在该保存窗口内未见明显脱离(净质心位移峰值 0.077 nm = 0.77 Å)。

它同样说明了自己的边界。 本窗口实测(1 nm = 10 Å):

指标 新跑 上游
蛋白 Cα RMSD(各自首帧参考) 均值 0.064 / 最大 0.089 nm = 0.64 / 0.89 Å 均值 0.111 / 最大 0.135 nm = 1.11 / 1.35 Å
配体受体拟合 RMSD 均值 0.085 / 最大 0.108 nm = 0.85 / 1.08 Å 均值 0.133 / 最大 0.165 nm = 1.33 / 1.65 Å
配体净质心位移 均值 0.041 / 最大 0.077 nm = 0.41 / 0.77 Å 均值 0.052 / 最大 0.100 nm = 0.52 / 1.00 Å

这些量级(亚埃到约 1.6 Å)在这个窗口里既不能判定"比上游更稳", 也不能判定"结合更强"。

关键限制(必须一起读): 上游 System XML 未随包提供 → 配体参数重新生成 + 我加了 HBonds 约束。 所以两条轨迹的差异是"物理模型 + 采样"混合的结果, 不能只归因于随机采样。这是一个条件性复现,不是逐比特复现。


6. 短模拟不能证明什么(逐条)

  1. 不能证明亲和力 / 结合自由能。RMSD 稳、接触多、氢键在,都不是 ΔG 或 Kd。 若需要,必须另跑 MM/GBSA 或自由能扰动,并单独交代假设。
  2. 不能证明收敛。18 ps 是极短窗口;F1 的块均值仍在缓升, 说明没有到达平稳的采样平台。上游 1 ns 也只覆盖单一随机实现的早期部分。
  3. 不能由本窗口推断解离行为。这里能说的只是:在保存的快照窗口内未见明显脱离。 本轨迹不能用来估计解离速率(或任何动力学速率):窗口只有 18 ps、 快照间隔 0.5 ps,且只有单条轨迹,既没有足够长的驻留时间采样, 也没有重复;"没有看到离开"与"不会离开"是两件不同的事。
  4. 不能证明筛选富集或成药性。本案例没有诱饵集、没有活性对照。
  5. 不能把两条轨迹的数值直接比大小。时间窗不同、物理参数不完全相同、 各自参考各自首帧。
  6. 不能把"跑通/文件哈希稳定/轨迹好看"当作成功标准。

7. 本次收尾做的修改(只列改动)

# 位置 问题 处理
1 figE · E1a/E1b E1b 直接继承 E1a 纵轴,裁掉上游的 22 和 30 共用纵轴改为覆盖两条数据实际范围(21.5–30.8);E2/E3 并列值加残基名次排序键,重绘次序确定
2 figA · A1 横轴类别文字重叠、底部说明压柱 改横向条形图 + 短标签;单位说明移到图注
3 figC · C4 标题过长被裁;"更大因为时间更长及首帧参考"无受控比较支持 换简短科学标题;RMSF 明确为对齐后围绕自身窗口均值的波动;总图注只对 C1/C2 声明"首帧为 0",并说明 Rg/RMSF 不适用;详细方法移入本报告图注
4 figF · F3 图例遮住蓝色直方图柱顶 图例移到绘图区外(下方)
5 文档 DCD 时间描述不准确 改为"DCD 无每帧时间数组,但带步长/步数信息,物理时间可恢复";NSET 随写入更新,不能断言关闭前一定是占位符
6 ../inputs/ligand_identity.json 该记录早于 CCD 核验,仍标 bond_orders_independently_verified: false 重跑仅感知/核验(不重参数化、不跑 MD),更新为 true 并附 CCD 明细;SDF 字节未变(sha256 一致),电荷缓存键仍有效

已确认无问题、未改动的部分:figD 标题换行、B3 温度标题、 C/D 各自的真实独立时间轴。

第二轮:纯文字精度修订(科学数据与图未动)

本轮不重算、不重绘(科学文件与图已在独立目录验证:七张表数值一致、 A–F 六张 PNG 逐像素一致、轨迹缺失时退出码 1)。只改下述可核实的表述:

# 位置 原文问题 处理
7 §4 · E1 接触判据写成 ≤ 0.45 nm 改回方法实际使用的 < 0.45 nm
8 §4 · E2 称两条轨迹接触残基数"中位数都是 26" 按 summary.json 改为新跑 26.5 / 上游 26
9 §4 · E3 称 03P:O3→Asp855:OD2 为"盐桥式氢键" 03P 的 O3 是中性羟基氧,不构成盐桥;只称满足几何判据的候选,并给出 36/36 帧
10 §5 "约 0.2 Å 级 RMSD"与实测 0.064–0.165 nm 不符 改为具体指标值并按 1 nm = 10 Å 正确换算
11 §4 · D2 称上游较大值"主要来自内部形变" D1/D2 是两种不同拟合,未做过位移分解;改为只分别描述各指标
12 §5 "口袋闭合""无解离" 无开闭指标 → 改为保存窗口内未见明显脱离
13 §3 / §4 · B1 "无漂移" 未做漂移检验 → 只保留未出现能量发散;并声明小 SD 不作为平衡的充分证据
14 §6 "0.5 ps 帧间隔无法解析快速的解离事件" 不宜泛称解离速度 → 改为本轨迹不能估计解离速率
15 §2.4 / charge_assignment_history.md 称首次尝试被"外层命令的超时杀掉" 与保存记录一致:编辑确认 PID 后主动 SIGTERM,非经证实的超时;并注明 logs/charge.log 只含第 2 次成功尝试
16 README.md §4 `MD=$(python3 scripts/md_env.py \ tail -1)实际得到 JSON 的}` 改用真实可用的函数调用(discover_md_python(verbose=False))或直接由用户提供环境 python 路径,并加 test -x 守卫
17 README.md §4 / §5 未区分验收范围 明确:已验收=轨迹重分析入口;未验收=新机器从零重跑的整链
18 report/REPORT.md 图/表链接指向会话绝对路径 改为相对 report/ 的 ../figures/、../tables/ 等,解压后可直接阅读

已确认无问题、未改动的部分:figD 标题换行、B3 温度标题、 C/D 各自的真实独立时间轴。

第三轮:审阅发现的问题(仅改文档)

独立审阅指出 §3 表与 §4·B2 的温度 SD 与保存产物不符。复核结果:审阅成立。

# 位置 问题 处理
19 §3 表 · §4·B2 温度"SD 0.94 K"在保存产物里找不到任何来源 改为 1.29 K(n = 36,ddof = 1),与 ../tables/summary.json 的 temperature_sd_K = 1.2935 及 ../tables/new_run_engine_log.csv 复算一致
20 §0 · §3 表 · ../README.md §2 总墙钟口径未写明,读起来像与阶段拆分矛盾 并非数据错误:33.6 分钟 = started_utc→finished_utc(2,014 s,进程启动至结束);33.3 分钟 = 各阶段计时之和(2,000 s)。现并列写出两个口径并注明差约 14 s 为导入/IO

关于 0.94 的核查:它在全部保存产物中只作为某帧温度值的一部分出现 (300.94…K 的子串),并非任何 SD;即该数字无来源,是转写错误。 图本身是对的:figB 的 B2 标题印的是 mean 299.93 K ± 1.29 (SD), 故本轮未重绘任何图。为避免再次发生,SD 现在连同 样本数 n 与自由度约定(ddof = 1)一起写出。

第 20 项不是数据更正,只是把原本没写明的口径补上: "总墙钟"是进程启动至结束(2,014 s),"阶段合计"是各阶段计时之和(2,000 s), 两者都对,并列写出即可。第一版修订曾误把 33.6 当作矛盾去"纠正",此处已撤回该判断。

同时按审阅意见补充:§8 的独立验证明确署名为"教程编辑独立验证"并给出记录路径, 不再表述为模块自行执行;并新增"读链接 vs 跑命令"的路径约定。


8. 复现与文件

独立验证的署名与出处

以下验证由教程编辑独立执行,不是本模块在本回合自行运行的。 记录保存在 outputs/teachopencadd-ten-20260911/validation/, 本报告只引用其结果:

验证记录 范围 结果
module07-preparation-verification.json 化学 / System / 拓扑(无 MD 重跑) passed: true;粒子数 58,047、约束 55,360、配体 63 原子;电荷缓存与粒子一致(最大差 0.0);总电荷 −8.5e−14;PME、无气压计;配体与 CCD 图/式一致
module07-run-verification.json 保存的 DCD / 日志 / PDB / System / checkpoint passed: true;步长与积分器一致(2.000000006 fs)、步号与日志一致、物理时间 0.5–18.0 ps、温度与受约束自由度一致、NVT 盒子恒定、checkpoint 可续跑一步
module07-analysis-verification.json 两轨迹独立复算(库成像 + NumPy Kabsch + 直接距离) passed: true;容差 0.5 pm 内全部指标吻合(最大差 ~1e−7 nm)、接触残基/占据率/计数一致、氢键几何逐一核对
module07-entry-verification.json 重分析入口在移动副本上重跑 passed: true;57.6 s;七张表行列一致、A–F PNG 逐像素一致、四格式齐全、轨迹缺失退出码 1

说明:module07-run-verification.json 只记录了温度的 min/max/mean, 没有记录温度 SD;上面 §3/§4 的 SD 1.29 K 来自 ../tables/summary.json 的 temperature_sd_K 与 ../tables/new_run_engine_log.csv 的 36 帧复算(ddof = 1),三者一致。

路径约定(读链接 vs 跑命令)

9. 数据来源与许可

TeachOpenCADD(Volkamer Lab)固定提交 4efdbd36723b5d07960142e15df393d7c6df2ba1,CC BY 4.0。 上游 topology.pdb / trajectory.xtc 是上游参考产物,本交付只作重分析; 新跑轨迹是本会话实际执行的真实动力学。配体身份取自 RCSB PDB 化学组分字典 (CCD 03P / TAK-285;3POZ 为其模型坐标来源)。