Expand description
确定性的分子 3D 初始构型生成 —— 给后续力场优化提供一个起点。
§目标,以及不是目标
目标:给一个分子生成一个三维构型,键长键角基本合理、立体化学正确、 没有原子重叠,交给后续力场优化。
不是目标:构象系综、构象搜索、接近实验结构、MMFF 优化后的几何分布。
§与 RDKit(ETKDG)的分岔点:去掉随机采样,不是去掉嵌入
RDKit 的流程是:界矩阵 → 三角光滑化 → 在区间里随机取一组距离 →
度量矩阵嵌入 → 误差函数精修 → 不合格就换一组随机距离重来(最多 10×N 次)。
病在“随机取一组距离“这一步:每一对原子各自独立取值,取出来的表
往往任何空间里都摆不出来 —— 好比“A 离 B 三米、B 离 C 三米、A 离 C 十米“,
写得下来,画不出来。RDKit 的应对是作废整次尝试重掷,而当病因是结构性的时候,
10×N 次重试会以同样的方式全部失败(实测失败分子上的触发率是 100%)。
它的失败分子与最慢分子高度重合,原因就在这里。
本算法保留嵌入,换掉那一步随机采样。
关键在 smooth:Floyd–Warshall 之后的上限矩阵 U 本身就是一个度量
(U_ij ≤ U_ik + U_kj 按构造恒成立),它是一张“画得出来“的距离表。
直接拿它当参考距离表,在同一批界矩阵上与 RDKit 那样的随机取法对照
(400 个分子,由 examples/bounds_oracle.rs 的判据三每次跑出来):
| 谱:负份额 | 谱:前三占绝对谱 | 几何:越界 RMS | 几何:1-2 键上越界 | 几何:长程下越界 | |
|---|---|---|---|---|---|
| 区间内随机取 | 0.284 | 0.463 | 0.812 Å | 63.2% | 3.2% |
直接用 U | 0.042 | 0.889 | 0.323 Å | 26.2% | 0.4% |
U 每一栏都赢,而且全程没有随机数 —— 同一个分子永远同一个答案,不需要 seed。
§但这张表读的时候有两个坑,都是踩过的
一、别用“越界的对数“当几何指标。 按对数算,U 是 29.3%、随机取 36.7%,
只好 1.25 倍;按 RMS 算是 0.323 / 0.812,好 2.5 倍。差别在于
长程对占了全部原子对的 87%,而那一档两种取法都基本没问题 ——
计数把差距摊没了,而且 0.1 Å 的阈值本身是个悬崖。报量,别只报数。
二、错不在长程,在短程。 直觉会以为“U 让每对取上限,所以长程会崩“,
实测正相反:长程的下越界只有 0.4%,而 1-2 键有 25.1% 低于下限、
26.2% 高于上限,1-3 角更差(26.6% / 41.2%)。
原因在经典 MDS 的目标本身:它拟合的是平方距离,于是大距离天然占主导 —— 一对 10 Å 的长程对差 0.1 Å,在目标里的分量远大于一根 1.5 Å 的键差 0.1 Å。 嵌入是拿键长去换长程的。
这一条不必在嵌入里修:精修的误差函数罚的是相对越界
(d²/u² − 1),同样 0.1 Å 的偏差,键上比 10 Å 的长程对贵约 7 倍,
天然会先把键拉回来。但判据必须分档看,否则“总体越界率下降“可以是
拿键长换来的 —— 长程占 87%,足够把这笔账盖住。
§曾经走错的两条路(别再走)
- 按分子类别切分支(无环走构造法、有环另说、超配位拒绝……):覆盖率停在 14.3%,而且每来一类分子就得加一个分支。
- 换个便宜初值来砍掉嵌入(生成树 NeRF 走一遍):实测慢约 2 倍,不是快 —— 嵌入不是成本,它给优化器一个好起点,省下的迭代远超那次特征分解。
§通用性的判定标准
新来一类分子(笼状、金属配合物、大环)时,要改的必须是约束表,
不是代码分支。嵌入器里不允许出现 is_macrocycle / is_metal 这类分子类别谓词 ——
这条可以直接 grep 进闸门。
§RDKit 那条“作废整次尝试“的规则有多狠
computeInitialCoords 只要发现一个原子的 sqD0i 低于阈值就把整次尝试
判死(DistGeomUtils.cpp:110-115)。判据三顺手记了这笔账:照搬这条规则,
400 个分子里 U 表会被打掉 265 个。
这不说明 U 差 —— 它的谱与几何两栏都赢 —— 而说明那条规则本身太粗:
sqD0i 是“到质心的平方距离“,U 把长程距离整体抬高,居中的原子于是算出负值。
那是长程被抬高的症状,不是摆不出来的证据。本算法因此不作废,只记账。
§精修阶段:力场里放全部 N² 对 —— 而这一条 RDKit 自己也这么干
RDKit 建力场时只放 u − l ≤ basinThresh(默认 5.0)的原子对
(DistGeomUtils.cpp 的 constructForceField)。实测这条过滤在语料上丢掉
16.8% 的原子对,逐分子中位 3.5%、p90 27.3%、最高 51.1% ——
柔性大分子上近一半的原子对在力场里没有任何约束,而那正是自穿会发生的地方。
但这不是我们对 RDKit 的批评,更不是分岔点。 Embedder.cpp:866 里,
只要走随机坐标那一支,RDKit 自己就把 basinThresh 设成 1e8(等于全部对),
旁边的注释写的理由与上面一字不差:
The basin threshold just gets us into trouble when we're using
random coordinates since it ends up ignoring 1-4 (and higher)
interactions. This causes us to get folded-up (and self-penetrating)
conformations for large flexible molecules上面那几个百分比是这段注释的定量版。真正的差别只有一处:RDKit 把
basinThresh = 1e8 与四维绑在一起放(同一个 useRandomCoords 分支也把
fourD 设 true),“全部对 + 只有三维“是它从不跑的组合。
Modules§
- bounds
- 界矩阵 —— 把化学事实翻译成“每一对原子之间必须离多远“。
- chiral
- 手性:有符号四点体积,以及嵌入之后那一次全局定向。
- embed
- 度量矩阵嵌入 —— 把一张距离表变成一组三维坐标。
- field
- 距离几何的误差函数 —— 精修阶段要最小化的那个目标。
- linalg
- 对称矩阵的特征分解 —— 手写循环 Jacobi 旋转。
- optimize
- L-BFGS —— 有限内存的拟牛顿下降,手写,零外部依赖。
- params
- 键长与键角的查表,以及查不到时的兜底。
- pipeline
- 端到端:分子进,一组三维坐标出。
- smooth
- 三角光滑化 —— 把一张自相矛盾的距离区间表收紧成自洽的。
- spread
- 确定性地把重合的原子分开 —— 必须在优化器之前跑。
- threading
- 自穿检测 —— 链有没有从环里穿过去、两根键有没有交叉。