RDKit 分子三维叠合完整指南：rdMolAlign 模块的 RMSD 对齐与 Open3DALIGN 实践-平芜编程栈

科学计算
科研
机器学习

【免费下载链接】rdkit

The official sources for the RDKit library

项目地址：https://gitcode.com/gh_mirrors/rd/rdkit

导读

rdkit.Chem.rdMolAlign是 RDKit 中负责分子三维叠合（alignment / superimposition）的官方 Python 模块，用于把探针分子（probe）的某个构象对齐到参考分子（reference）的某个构象上，并返回 RMSD 值、变换矩阵与原子映射。本文以 rdkit.Chem.rdMolAlign.rst 对应的自动文档页为骨架，结合仓库内 C++ 实现与 Python 测试，完整讲解AlignMol、GetBestRMS、CalcRMS、BestAlignmentParams、AlignMolConformers以及 Open3DALIGN（O3A）系列接口。读完本文，你将能够：对任意两个分子/构象做最小 RMSD 叠合并取回变换矩阵；在考虑分子对称性的前提下计算"最优"RMSD；批量对齐多构象分子并生成 RMSD 距离矩阵；以及使用 MMFF94 或 Crippen 原子性质驱动的 O3A 形状对齐。

一、模块概览：它解决什么问题

分子三维叠合是计算化学中高频出现的基础操作，典型场景包括：

比较两个配体的结合构象与对接姿势（docking pose）之间的偏差；
将多构象分子统一对齐到同一参考系，便于视觉比较与后续分析；
计算同一分子不同构象之间的 RMSD 距离矩阵，用于构象聚类；
基于分子形状与原子性质的配体形状相似性打分（Open3DALIGN）。

rdMolAlign模块的定位在源码注释中写得很明确："Module containing functions to align a molecule to a second molecule"（见 Wrap/rdMolAlign.cpp）。其 Python 绑定由 Code/GraphMol/MolAlign/Wrap/rdMolAlign.cpp 用 Boost.Python 生成，底层核心算法位于 Code/GraphMol/MolAlign/AlignMolecules.cpp 与 Code/GraphMol/MolAlign/O3AAlignMolecules.h，测试用例集中在 Wrap/testMolAlign.py。

模块内的函数大致分三组：基础 RMSD 对齐（AlignMol、GetAlignmentTransform、AlignMolConformers、RandomTransform）、对称感知最优对齐（GetBestRMS、GetBestAlignmentTransform、CalcRMS、BestAlignmentParams）、O3A 形状对齐（GetO3A、GetCrippenO3A、GetO3AForProbeConfs、GetCrippenO3AForProbeConfs）。

二、基础 RMSD 对齐 API

2.1 AlignMol：对齐并直接修改探针坐标

AlignMol是使用频率最高的函数：它计算将探针分子指定构象对齐到参考分子指定构象的 3D 变换，使指定原子集合的均方根距离最小，并把该变换直接应用到探针分子上，然后返回 RMSD 值。

from rdkit import Chem from rdkit.Chem import rdMolAlign rmsd = rdMolAlign.AlignMol(prbMol, refMol, prbCid=-1, refCid=-1, atomMap=(), weights=(), reflect=False, maxIters=50)

参数语义（与 Wrap/rdMolAlign.cpp 的文档字符串一致）：

参数	默认值	说明
`prbMol`	必填	待对齐的探针分子
`refMol`	必填	对齐所参照的参考分子
`prbCid`	`-1`	探针分子中参与对齐的构象 ID，`-1`表示第一个构象
`refCid`	`-1`	参考分子中作为对齐目标的构象 ID
`atomMap`	`()`	探针原子索引与参考原子索引的配对列表，格式为`((prbIdx, refIdx), ...)`；若不指定，则通过子结构匹配自动生成
`weights`	`()`	每个原子对的可选权重，长度必须与`atomMap`长度一致
`reflect`	`False`	若为`True`，先对探针构象做镜像反射再对齐
`maxIters`	`50`	最小化 RMSD 过程中的最大迭代次数

注意：AlignMol会就地修改prbMol的坐标（C++ 层的alignMol在计算变换后调用applyTransform写入构象，见 AlignMolecules.cpp）。若不想改动原分子，请先Chem.Mol(prbMol)复制一份。

2.2 GetAlignmentTransform：只取变换，不碰坐标

如果你只需要 RMSD 与变换矩阵（例如用于后续自建流水线），用GetAlignmentTransform。它不会修改任何分子的坐标，返回一个二元组：

rmsd, trans = rdMolAlign.GetAlignmentTransform(prbMol, refMol, prbCid=-1, refCid=-1, atomMap=(), weights=(), reflect=False, maxIters=50)

rmsd：浮点数，对齐后的均方根距离；
trans：numpy.ndarray，形状(4, 4)的齐次变换矩阵（底层由RDGeom::Transform3D导出，见 Wrap/rdMolAlign.cpp）。

测试 testMolAlign.py 中明确断言了trans是numpy.ndarray，且AlignMol与GetAlignmentTransform在同一对分子上得到的 RMSD 基本一致（0.6578与0.6579，差异来自浮点舍入）。

2.3 AlignMolConformers：分子内部多构象对齐

当你想把同一个分子的多个构象统一叠到第一个构象上时，使用AlignMolConformers：

rdMolAlign.AlignMolConformers(mol, atomIds=(), confIds=(), weights=(), reflect=False, maxIters=50, RMSlist=None)

以mol的第一个构象为参考，逐个对齐其余构象；
atomIds：参与对齐的原子索引列表，默认全部原子；
confIds：要对齐的构象 ID 列表，默认全部构象；
RMSlist：若传入一个 list，函数会把"参考构象到各对齐构象"的 RMSD 追加进去（见 Wrap/rdMolAlign.cpp）。

测试 testMolAlign.py 展示了完整用法：先用 UFF 力场最小化big_aromat.sdf中的每个构象，再用AlignMolConformers(mol, aids)对齐指定原子，并断言这些原子在各构象中的位置趋于一致（容差 0.5 Å）；随后通过RMSlist收集 RMSD，并断言列表长度为构象数 - 1。

2.4 RandomTransform：随机变换探针

RandomTransform对分子施加一次随机刚体变换，常用于构造对齐算法的"挑战输入"（测试中用它打乱探针再验证能否找回正确叠合）：

rdMolAlign.RandomTransform(mol, cid=-1, seed=-1)

seed默认为-1（不设种子）；传入非负整数可复现随机结果（见 Wrap/rdMolAlign.cpp）。

三、考虑对称性的最优对齐：GetBestRMS 与 GetBestAlignmentTransform

基础对齐依赖固定的原子映射。但分子常常具有对称性（如苯环的 6 重旋转、羧基两个氧的等价性），此时"最小 RMSD"未必出现在你手工给定的映射上。为此，RDKit 提供了对称感知的最优对齐 API：它们会遍历所有等价的子结构匹配，找出 RMSD 最小的那个映射。

3.1 GetBestRMS

rmsd = rdMolAlign.GetBestRMS(prbMol, refMol, prbId=-1, refId=-1, map=(), maxMatches=1000000, symmetrizeConjugatedTerminalGroups=True, weights=(), numThreads=1)

参数	默认值	说明
`map`	`()`	可选的多组原子映射（list of lists of`(prbIdx, refIdx)`）；不提供时用子结构搜索自动生成
`maxMatches`	`1000000`	`map`为空时，`SubstructMatch()`最多返回的匹配数量上限
`symmetrizeConjugatedTerminalGroups`	`True`	若为`True`，共轭末端官能团（如硝基、羧酸根）按对称方式处理
`weights`	`()`	映射对应的权重
`numThreads`	`1`	计算线程数

副作用提醒：GetBestRMS的文档字符串明确写着"As a side-effect, the probe molecule is left in the aligned state"——即它同样会把探针留在对齐后的状态，并且官方建议："for some molecules it will lead to 'combinatorial explosion' especially if hydrogens are present"（见 Wrap/rdMolAlign.cpp）。因此：

需要保持原子顺序不变时，请用AlignMol（或GetAlignmentTransform）；
含氢较多、对称性高的分子，maxMatches的穷举可能急剧膨胀，应权衡或预置map。

3.2 GetBestAlignmentTransform

与GetBestRMS相同的最优搜索，但不修改坐标，且返回更完整的信息：

rmsd, bestTrans, bestMatch = rdMolAlign.GetBestAlignmentTransform( prbMol, refMol, prbCid=-1, refCid=-1, map=(), maxMatches=1000000, symmetrizeConjugatedTerminalGroups=True, weights=(), reflect=False, maxIters=50, numThreads=1)

返回三元组(最优 RMSD, 4×4 变换矩阵, 最佳原子映射)。测试 testMolAlign.py 验证了：bestTrans是numpy.ndarray，bestMatch长度等于参考分子原子数，且该 RMSD 与GetBestRMS的结果一致。

3.3 BestAlignmentParams：参数对象

除散参数写法外，两组函数都提供等价的参数对象重载。BestAlignmentParams由 Wrap/rdMolAlign.cpp 绑定，属性如下：

属性	默认值	说明
`maxMatches`	`1000000`	子结构匹配数量上限
`symmetrizeConjugatedTerminalGroups`	`True`	共轭末端官能团对称化
`ignoreHs`	`True`	对齐时忽略氢原子（注意：显式出现在`map`中的原子始终会被使用）
`numThreads`	`1`	线程数
`map`	`()`	原子映射集合
`weights`	`None`	权重向量

注意一个细节：散参数版本的GetBestRMS/GetBestAlignmentTransform内部会把ignoreHs固定为False（见 AlignMolecules.h 内联实现），而BestAlignmentParams默认ignoreHs=True——两者默认行为并不完全相同，测试 testMolAlign.py 也展示了ignoreHs对结果的影响（同为默认参数对象时 RMSD 为1.8100，而ignoreHs=False时回到2.43449）。

参数对象用法：

params = rdMolAlign.BestAlignmentParams() params.map = matchesPruned # 预置原子映射 params.weights = weights # 预置权重 rmsd = rdMolAlign.GetBestRMS(prbMol, refMol, params, prbId=0, refId=0)

3.4 CalcRMS：原地计算 RMSD

CalcRMS与GetBestRMS的最大区别是不做对齐："probe molecules are not aligned to the reference ahead of the RMS calculation"，直接以两分子当前坐标计算对称感知的最优 RMSD。官方文档特别指出它适合计算"对接姿势与共结晶配体之间的 RMSD"（见 Wrap/rdMolAlign.cpp）。

rmsd = rdMolAlign.CalcRMS(prbMol, refMol, prbId=-1, refId=-1, map=(), maxMatches=1000000, symmetrizeConjugatedTerminalGroups=True, weights=())

测试 testMolAlign.py 给出同一对分子的三个典型值，直观展示了三者的差异：CalcRMS（不移动探针）为2.6026，注释说明AlignMol会得到2.50561，而对称感知的GetBestRMS得到最优值2.43449。

四、构象集合间的 RMSD 矩阵

4.1 GetAllConformerBestRMS：单分子内部距离矩阵

返回同一分子任意两个构象之间的最优 RMSD（内部调用getBestRMS），结果是一个对称距离矩阵的上三角展平：

vals = rdMolAlign.GetAllConformerBestRMS(mol, numThreads=1, map=(), maxMatches=1000000, symmetrizeConjugatedTerminalGroups=True, weights=())

返回顺序为[(1,0),(2,0),(2,1),(3,0),(3,1),(3,2),...]，共n×(n-1)/2个值（n 为构象数）。测试 testMolAlign.py 验证了返回值长度，并断言origVals[0]与GetBestRMS(mcp, mcp, 0, 1)一致，同时确认numThreads=4的多线程结果与单线程完全一致。该结果可直接喂给聚类算法做构象聚类。

4.2 GetAllConformerBestRMSToRef：两组分子间的交叉矩阵

当探针与参考分别是两个多构象分子时，GetAllConformerBestRMSToRef返回参考所有构象对探针所有构象的 RMSD：

vals = rdMolAlign.GetAllConformerBestRMSToRef(prbMol, refMol, params=None)

返回顺序为[(0,0),(0,1),(0,2),(1,0),(1,1),(2,0),...]，其中第一个索引是 refMol 的构象 ID，第二个是 prbMol 的构象 ID（见 Wrap/rdMolAlign.cpp）。测试 testMolAlign.py 使用butane_ref.sdf/butane_prb.sdf校验了数值（如首个值为0.19474）。

五、Open3DALIGN（O3A）形状对齐

O3A（Open3DALIGN）是一种不依赖显式子结构匹配的叠合方法：它基于原子性质（MMFF94 原子类型与电荷，或 Crippen logP 原子贡献）构建"成本矩阵 + 加权 + 打分"体系，通过线性分配问题（LAP）与 Sinkhorn 距离矩阵（SDM）迭代寻找原子对应关系，从而对骨架不同但形状相似的分子进行形状对齐。核心类O3A与相关数据结构（O3AConstraintVect、LAP、SDM、MolHistogram）都在 O3AAlignMolecules.h 中。

5.1 GetO3A（MMFF94 版）与 GetCrippenO3A

o3a = rdMolAlign.GetO3A(prbMol, refMol, prbPyMMFFMolProperties=None, refPyMMFFMolProperties=None, prbCid=-1, refCid=-1, reflect=False, maxIters=50, options=0, constraintMap=(), constraintWeights=()) o3a = rdMolAlign.GetCrippenO3A(prbMol, refMol, prbCrippenContribs=(), refCrippenContribs=(), prbCid=-1, refCid=-1, reflect=False, maxIters=50, options=0, constraintMap=(), constraintWeights=())

两者的区别仅在于驱动对齐的原子性质：

MMFF 版：基于 MMFF94 原子类型与电荷。prbPyMMFFMolProperties/refPyMMFFMolProperties可由ChemicalForceFields.MMFFGetMoleculeProperties(mol)提前算好并复用；若传None，函数会在内部调用MMFFMolProperties现场计算，此时若分子缺少 MMFF94 参数会抛出ValueError: missing MMFF94 parameters for probe/reference molecule（见 Wrap/rdMolAlign.cpp，测试 test7 也专门覆盖了该错误路径）。
Crippen 版：基于 Crippen logP 原子贡献。prbCrippenContribs/refCrippenContribs可通过rdMolDescriptors._CalcCrippenContribs(mol, True)预先计算（见 Wrap/rdMolAlign.cpp）。

5.2 O3A 对象方法

GetO3A/GetCrippenO3A返回O3A对象（Python 类rdMolAlign.O3A），提供五个方法（见 Wrap/rdMolAlign.cpp）：

方法	说明
`Align()`	执行叠合，返回对齐后的 RMSD；探针坐标被就地更新
`Trans()`	返回`(RMSD, 4×4 变换矩阵)`，不修改坐标
`Score()`	返回 O3A 形状打分（越高越好）
`Matches()`	返回 O3A 找到的原子映射`(prbIdx, refIdx)`列表
`Weights()`	返回 O3A 计算出的权重向量

典型用法（与测试 test5/test6 一致）：

pyO3A = rdMolAlign.GetO3A(prbMol, refMol) score = pyO3A.Score() rmsd = pyO3A.Align() # 探针被对齐

测试 testMolAlign.py 对ref_e2.sdf中全部分子逐个做 MMFF O3A 叠合，累计打分约6942、平均 RMSD 约0.345；Crippen 版（test10/test11）累计打分约4918、平均 RMSD 约0.304。这些数值可作为功能自检的参照。

5.3 约束对齐：constraintMap 与 constraintWeights

O3A 支持把若干原子对作为"硬约束"注入对齐：约束原子对会优先匹配，权重越大越被强制靠近。约束的底层类型是 O3AAlignMolecules.h 中的O3AConstraint/O3AConstraintVect，默认约束权重常量为O3_DEFAULT_CONSTRAINT_WEIGHT = 100.0。

pyO3A = rdMolAlign.GetO3A(m3, m1, constraintMap=[[cIdx, nIdx]]) pyO3A.Align()

约束使用有三条硬性规则（Python 层会直接抛ValueError，见 Wrap/rdMolAlign.cpp）：

约束索引必须在各自分子的原子数范围内（Constrained atom idx out of range）；
约束原子必须是重原子（氢原子会报Constrained atoms must be heavy atoms）；
constraintWeights的长度必须与constraintMap一致（The number of weights should match the number of constraints）。

测试 testMolAlign.py 给出一个精巧的验证：两份完全相同的共平面 4-苯基吡啶，正常叠合时对应原子距离为 0；而强制把吡啶氮与苯环对位碳配对（constraintMap=[[cIdx, nIdx]]）后，对应原子距离变为 7 Å——说明约束确实"强制"了错误的配对。带权约束与局部优化的组合（test9/test13）则验证了权重越大约束越紧（如权重100.0时距离0.3-0.4 Å，权重10.0/0.1时距离2.7-3.2 Å）。

5.4 options 位标志与批量构象 API

options参数是一个位标志整数，语义（见 Wrap/rdMolAlign.cpp）：

最低 2 位编码计算精度：0= 最大精度（默认），3= 最低精度；
第 3 位（值4）触发仅局部优化：跳过成本矩阵计算，直接在当前位置附近做局部寻优（测试 test9/test13 中options=4的用法）。

批量构象 API 用于一次性把探针的所有构象分别对齐到参考：

o3as = rdMolAlign.GetO3AForProbeConfs(prbMol, refMol, numThreads=1, prbPyMMFFMolProperties=None, refPyMMFFMolProperties=None, refCid=-1, reflect=False, maxIters=50, options=0, constraintMap=(), constraintWeights=()) o3as = rdMolAlign.GetCrippenO3AForProbeConfs(prbMol, refMol, numThreads=1, ...)

返回一个 O3A 对象列表（长度等于探针构象数）。numThreads只有 RDKit 以线程支持编译时生效；设为0时使用系统支持的最大线程数（见 Wrap/rdMolAlign.cpp）。测试 test15 断言批量结果与逐个GetO3A(..., prbCid=i)的 RMSD/Score 逐位一致（容差 1e-6），test16 则用numThreads=4验证多线程路径返回的构象数正确。

六、底层实现原理：从 Python 到 C++

对齐链路分三层：

Python 绑定层（Wrap/rdMolAlign.cpp）：负责参数解析与校验（如translateAtomMap解析映射、translateDoubleSeq解析权重、原子索引范围检查），并在NOGIL块中释放 GIL 调用 C++ 核心，保证多线程场景下 Python 解释器不被阻塞。
对齐核心层（AlignMolecules.cpp）：
- getAlignmentTransform：若未提供atomMap，先执行SubstructMatch(refMol, prbMol, ...)自动生成原子映射（AlignMolecules.cpp），再调用 Kabsch 类最小二乘求解变换并迭代优化 RMSD；
- getBestRMSInternal：在SubstructMatch返回的所有匹配上逐一计算 RMSD，取最优值（AlignMolecules.cpp），这正是"对称感知"的来源；
- alignMol：getAlignmentTransform之后把变换写回探针构象。
O3A 形状层（O3AAlignMolecules.h）：O3A类通过LAP（线性分配，computeCostMatrix/computeMinCostPath）与SDM（Sinkhorn 距离矩阵迭代）在原子间建立软对应，再用o3aMMFFCostFunc/o3aCrippenCostFunc等成本函数驱动匹配。头文件开头的常量揭示了算法的关键超参，例如约束默认权重O3_DEFAULT_CONSTRAINT_WEIGHT = 100.0、电荷权重O3_CHARGE_WEIGHT = 10.0、Crippen 权重O3_CRIPPEN_WEIGHT = 10.0、RMSD 收敛阈值O3_RMSD_THRESHOLD = 1.0e-04等。

理解这一点有助于正确使用 API："最优"（Best）系列天然依赖子结构匹配，因此对含氢分子要注意maxMatches上限与ignoreHs的取舍；而 O3A 系列不要求子结构对应，更适合骨架位移较大的形状比对。

七、测试与数值验证

仓库自带的 Python 测试 Wrap/testMolAlign.py 与 C++ 测试 testMolAlign.cpp、testO3AAlign.cpp 覆盖了全部主要路径，可作为功能自检基准（测试数据位于 test_data）：

测试	场景	期望数值
test1	`AlignMol`（1oir.mol 对 1oir_conf.mol）	RMSD ≈`0.6578`
test2	指定 6 对`atomMap`后对齐	RMSD ≈`0.8525`
test3	指定映射 + 权重`(1,1,1,1,1,2)`	RMSD ≈`0.9513`
test4	`AlignMolConformers`+`RMSlist`	列表长度 = 构象数 - 1
test5/6	MMFF O3A 全库叠合	累计 Score ≈`6942`，平均 RMSD ≈`0.345`
test10/11	Crippen O3A 全库叠合	累计 Score ≈`4918`，平均 RMSD ≈`0.304`
test17	`CalcRMS`vs`GetBestRMS`	`2.6026`vs`2.43449`
test18	共轭末端对称化开关（丙酸根）	开：`0`；关：`0.747`
test20	`GetAllConformerBestRMSToRef`（丁烷）	首值`0.19474`

其中 test18 尤其直观：两个只有羧酸根书写顺序不同的丙酸根，默认symmetrizeConjugatedTerminalGroups=True时 RMSD 为0，关闭后变为0.747——这正是"共轭末端官能团对称化"参数价值的实证。

八、典型应用场景与选型建议

对接姿势比较：用CalcRMS原地计算对接姿势与共结晶配体的 RMSD，不引入额外刚体变换的干扰；
构象聚类 / 多样性分析：GetAllConformerBestRMS生成构象距离矩阵，配合GetBestRMS的对称感知特性保证距离度量不受原子编号顺序影响；
多构象统一叠合：AlignMolConformers一键把整套构象对齐到首构象，便于截图与目视比较；
骨架不同的形状比对：GetO3A（MMFF）或GetCrippenO3A输出Score()用于配体形状相似性排序，Matches()可反查形状等价的原子对应；
强制位点匹配：利用 O3A 的constraintMap/constraintWeights把药效团关键原子钉在参考位点上。

选型上：需要"不改变原子顺序、明确映射"时选AlignMol/GetAlignmentTransform；需要"对称无关的最小 RMSD"时选GetBestRMS/GetBestAlignmentTransform；需要"形状驱动、允许骨架差异"时选 O3A 系列。结合 AlignMolecules.h 与 O3AAlignMolecules.h 的源码注释，即可在 RDKit 提供的这套完整对齐工具箱中找到最适合场景的组合。