企业官网建设流程全解析

做分子性质预测、虚拟筛选或者化合物结构分析的时候，我最常用也最离不开的一个工具就是 RDKit。不管是手里拿到的还只是一行 SMILES，还是要处理已经带有三维坐标的文件，绝大部分工作最后都能落到这个工具包上。简单说，RDKit 就是化学信息学领域处理 2D 结构式和 3D 构象的一把瑞士军刀，它能帮你完成从分子读入、清洗、画图、子结构匹配、指纹相似度计算，到三维构象生成、力场优化、分子叠加对齐这一整条链路。

这篇文章主要是想聊聊我平时用 RDKit 处理 2D 和 3D 分子的完整工作流，包括每一步为什么要这么做、参数到底怎么选、以及那些文档里不会明说但实际一定会踩的坑。不管你是刚接触化学信息学的学生，还是已经在做药物设计、材料筛选的开发人员，这套流程应该都能直接拿去用。

1. 环境准备与分子读入：从一行 SMILES 开始的完整流程

1.1 装环境为什么优先选 conda-forge

RDKit 的安装方式其实有不少，pip 也能装，但我个人强烈建议用 conda 从 conda-forge 通道装。原因很简单：RDKit 的底层依赖了 Boost、NumPy 这些比较敏感的编译环境，用 pip 装预编译包偶尔会遇到 OpenMP 版本不匹配或者某些绘图功能缺失的问题。conda-forge 上的版本更新及时，而且会把 cairo、boost、numpy 这些配套依赖一起解决掉，尤其是后面要用到 2D 描绘功能时，没有 cairo 支持会特别难受。

conda create -n rdkit-env -c conda-forge rdkit python=3.11 conda activate rdkit-env

装好之后我通常建议先跑一个自检脚本，确认关键模块都能正常导入，同时看一眼版本号，方便后续排查问题：

from rdkit import Chem from rdkit.Chem import AllChem, Draw, rdMolDraw2D print(Chem.RDKFingerprint)

如果你的环境已经有其他深度学习依赖，也可以直接把 rdkit 加进去，不用刻意隔离。只是要注意一点：不要把 conda-forge 和 defaults 通道混得太乱，否则解算依赖时容易卡半天。

1.2 读入分子时的常见坑位

RDKit 读入 2D 分子最基础的方式就是Chem.MolFromSmiles，一行 SMILES 就能生成一个Mol对象。比如乙醇是CCO，苯是c1ccccc1，布洛芬大概是CC(C)Cc1ccc(C(C)C(=O)O)cc1。

但你如果直接拿这个函数去处理一批数据，很快就会发现它有个让人头疼的地方——解析失败时不会抛异常，而是默默返回None。这就意味着你写循环批量读 SMILES 的时候，不能想当然地认为每个分子都成功解析了，一定要做空值判断，并且把失败记录打出来，不然程序跑到一半才发现某个分子丢了，排查起来很麻烦。

from rdkit import Chem smiles_list = ['CCO', 'c1ccccc1', 'invalid_smiles'] for smi in smiles_list: mol = Chem.MolFromSmiles(smi) if mol is None: print(f'解析失败: {smi}') continue print(f'成功: {smi}, 原子数 {mol.GetNumAtoms()}')

还有一个新手容易忽略的点：RDKit 内部的Mol对象本质上是一张图，原子是节点，键是边。它默认不显式保存氢原子，所以你去数原子数时，乙醇CCO只有 3 个重原子，而不是 9 个原子。这一点在后续做 3D 构象生成时尤其重要，因为三维坐标计算时通常需要先通过Chem.AddHs把氢原子显式加回来。

2. 2D 层面到底能做什么：子结构匹配、指纹与相似度

2.1 子结构搜索不是简单的字符串匹配

处理 2D 分子最频繁的需求就是子结构搜索。比如你想在一批化合物里筛出所有含有羧基的分子，天然想法可能是直接在 SMILES 里找C(=O)O这种字符串，但这种方法很快会出问题，因为同一个官能团在不同的分子里写出来可能长得完全不一样。RDKit 的方式是让你定义 SMARTS 模式，它是一个描述“原子类型 + 连接关系 + 环境”的模式语言。

举个例子，比较稳健的羧基 SMARTS 可以写成[CX3](=O)[OX2H1]。这里CX3表示碳原子连接了 3 个其他原子，=O是羰基氧，OX2H1是连了 2 个原子的羟基氧。比单纯找字符串靠谱得多。

carboxyl = Chem.MolFromSmarts('[CX3](=O)[OX2H1]') mol = Chem.MolFromSmiles('CC(C)Cc1ccc(C(C)C(=O)O)cc1') match = mol.GetSubstructMatch(carboxyl) print(match) # 返回匹配到的重原子索引元组

实际操作中还有一个非常实用的细节：如果同一个模式在一个分子里出现多次，你希望全部找出来而不是只取第一组，就要用GetSubstructMatches，它返回的是所有匹配结果。做构象叠加或骨架对齐时，这个功能特别有用，因为一个分子里可能有两个甚至更多的候选匹配位点，需要枚举后取最优的那一个。

2.2 指纹计算：怎么选 ECFP 还是 MACCS

相似度计算是 2D 层面另一个核心操作。做虚拟筛选、聚类、去重都离不开分子指纹。RDKit 里最常用的两类指纹，一类是 Morgan 指纹，对应大家常说的 ECFP；另一类是 MACCS keys。

Morgan 指纹的本质是把每个原子周围的环境做“圆形展开”，半径越大，囊括的原子环境和键合信息越远。常用的 ECFP4 其实就是 Morgan 指纹半径等于 2、再按哈希方式映射到固定位长的位向量。我的经验是：做药化分子相似性搜索和活性预测时，Morgan 指纹通常表现最好，因为它对分子的局部环境描述细腻；而 MACCS keys 只有 166 位，更多用于快速粗筛、骨架多样性分析。

from rdkit import DataStructs from rdkit.Chem import AllChem mol1 = Chem.MolFromSmiles('CCO') mol2 = Chem.MolFromSmiles('CCOC') fp1 = AllChem.GetMorganFingerprintAsBitVect(mol1, radius=2, nBits=2048) fp2 = AllChem.GetMorganFingerprintAsBitVect(mol2, radius=2, nBits=2048) sim = DataStructs.TanimotoSimilarity(fp1, fp2) print(sim)

批量计算相似度时，千万别写两层 for 循环然后逐个调TanimotoSimilarity，数据量一大就慢得离谱。RDKit 提供了DataStructs.BulkTanimotoSimilarity，一次算出一个指纹和一堆指纹的相似度列表，速度提升非常明显。

2.3 2D 描绘：让别人一眼看懂你的分子

很多时候你不需要把分子讲得多复杂，直接画一张 2D 结构图就够了。RDKit 的rdMolDraw2D模块支持生成 PNG 或 SVG。我最常用的是MolDraw2DCairo，因为它能直接输出位图，放在论文和 PPT 里很方便。

from rdkit.Chem import rdMolDraw2D mol = Chem.MolFromSmiles('CC(C)Cc1ccc(C(C)C(=O)O)cc1') drawer = rdMolDraw2D.MolDraw2DCairo(400, 300) rdMolDraw2D.PrepareAndDrawMolecule(drawer, mol) drawer.FinishDrawing() drawer.WriteDrawingText('ibuprofen.png')

如果目的不是示意图而是讲解机制，我会把子结构匹配到的原子高亮出来，或者把原子索引号画出来，这样在一些远程协作场景里描述“第 7 号原子发生了什么反应”就会清楚得多。

3. 3D 构象生成的完整链路：从加氢到坐标输出

3.1 为什么不能直接把 2D 结构拉成 3D

很多刚接触的人以为 3D 结构就是给 2D 坐标加一个 Z 值就行，这是个很大的误解。分子的三维构象不是平面结构简单抬高，而是要在满足键长、键角、二面角优先倾向以及非键原子不重叠的前提下，找到一套合理的空间坐标。

RDKit 生成 3D 构象的核心逻辑是距离几何法。它先根据 2D 分子图的连接关系估算出原子两两之间的距离上下界，然后在一个高维距离矩阵里随机采样，再把采样结果投影到三维空间。这一步做完后，分子的原子坐标已经有了，但键长键角通常不够理想，所以还要用力场优化一下。

一个容易踩的坑：做 3D 之前一定先调用Chem.AddHs(mol)把氢原子补上。因为 2D 表示里的氢通常是隐式的，如果不加氢就生成构象，你会得到一个没有氢、看起来好像也合理但实际物理上很怪的结构。而且后面很多力场参数是依赖氢原子的，少了氢就没法正常算。

3.2 EmbedMolecule 的关键参数与随机种子

RDKit 里实际生成单构象的函数是AllChem.EmbedMolecule。推荐的做法是使用带经验知识约束的 ETKDG 方法，也就是设置useExpTorsionAnglePrefs=True和useBasicKnowledge=True，它会参考已知的扭转角偏好和常见环构象库，生成的构象合理度明显更高。

mol = Chem.MolFromSmiles('CCO') mol = Chem.AddHs(mol) AllChem.EmbedMolecule(mol, randomSeed=42, useExpTorsionAnglePrefs=True, useBasicKnowledge=True) AllChem.MMFFOptimizeMolecule(mol) conformer = mol.GetConformer() pos = conformer.GetAtomPosition(0) print(pos.x, pos.y, pos.z) # 单位是埃

这里有个很重要的经验：一定要设置randomSeed。如果不设置，程序内部会取系统随机数，结果就是每次运行生成的构象可能都不一样。做实验、写论文、跑流程的时候，如果你不固定随机种子，别人根本没法复现你的结果。我一般习惯固定为 42 或者某个跟数据集编号有关的值。

MMFFOptimizeMolecule这一行是力场优化，默认选用 MMFF94 力场。MMFF94 对有机分子尤其是类药分子的几何优化效果通常不错。如果 MMFF 因为元素类型或参数缺失报错，可以退一步用UFFOptimizeMolecule，通用性更强但几何精度略差。

3.3 多构象生成并不是越多越好

有时候你需要的不只是一个构象，而是尽量覆盖分子可能采取的多种空间形状，比如做构象搜索或者 3D 相似性分析。这时可以用EmbedMultipleConfs一次性生成多个构象，然后逐一优化。

confs = AllChem.EmbedMultipleConfs(mol, numConfs=50, randomSeed=1, useExpTorsionAnglePrefs=True, useBasicKnowledge=True) res = AllChem.MMFFOptimizeMoleculeConfs(mol)

但我要提醒一句：生成多少个构象，取决于你的下游任务。虚拟筛选的 3D 对接，每个配体生成 50 到 200 个构象通常就够用了；如果你只是想要一个合理的展示结构，生成一个就完了。我见过有人一上来就生成 500 个构象，最后发现大部分都聚成一类，白白浪费时间。

生成之后还有一个重要问题：你怎么判断这些构象是不是真的“不一样”？这时候不能只看坐标，而要算两两之间的 RMSD。RMSD 低于某个阈值（比如 0.5 埃）的两个构象就可以认为是重复的，后续用聚类去冗余才有意义。

3.4 3D 结构如何写出到文件

得到 3D 构象后，最常用的输出格式是 SDF。SDF 的好处是不仅能保存分子连接关系，还能保存每个原子的三维坐标以及分子的属性字段。

from rdkit.Chem import SDWriter mol = Chem.MolFromSmiles('CCO') mol = Chem.AddHs(mol) AllChem.EmbedMolecule(mol, randomSeed=42) AllChem.MMFFOptimizeMolecule(mol) writer = SDWriter('ethanol_3d.sdf') writer.write(mol) writer.close()

写文件之前一定要确保分子确实带有坐标。判断方法很简单：mol.GetNumConformers() > 0就说明有构象。如果分子没有任何构象，写出来也会报错或丢失 3D 信息。

4. 让 2D 和 3D 协同工作：子结构对齐与构象比较

4.1 子结构匹配如何转成 3D 操作

前面说 2D 子结构匹配是找原子索引，而一旦加了氢、生成了构象，原子的编号顺序就会变：原来 2D 里第 0 号原子是碳，加氢之后第 1 号原子可能是氢。很多人在这里栽跟头，原因就是没搞清楚索引映射。

正确做法是，在加氢之前先做子结构匹配，记录重原子索引，再对Mol的副本做 3D 操作。如果你要在生成 3D 结构后继续用 2D 的匹配结果，最好先把匹配到的原子索引存下来，因为那组索引在AddHs之前是稳定的，之后就不一定了。

mol = Chem.MolFromSmiles('CC(C)Cc1ccc(C(C)C(=O)O)cc1') substructure = Chem.MolFromSmarts('[OH1]') match = mol.GetSubstructMatch(substructure) mol3d = Chem.AddHs(mol) AllChem.EmbedMolecule(mol3d, randomSeed=42)

这里你会发现match里保存的是重原子的索引，在mol和mol3d上如果原子顺序没变，前几个重原子的索引还是可以对应的。但如果你之后要移除氢、再做凯库勒化之类的操作，索引可能就不可靠了，所以最稳妥的策略是优先处理完所有 2D 逻辑之后，再统一转向 3D。

4.2 分子对齐与 RMSD 计算

3D 层面最常用的操作之一是把两个分子的构象放到同一坐标系下比较。比如你想知道分子 A 的活性构象和分子 B 的某种构象在空间上差多少，就可以找一个公共骨架做对齐，然后算 RMSD。

RDKit 的AlignMol可以完成这件事：你提供参考分子和目标分子，以及两组原子索引的映射关系，它会通过最小二乘法旋转平移坐标，让对应原子的空间距离尽量小，然后返回 RMSD 值。

from rdkit.Chem import AllChem core = Chem.MolFromSmarts('c1ccccc1') mol1 = Chem.MolFromSmiles('c1ccc(C)cc1') mol2 = Chem.MolFromSmiles('c1ccc(C(=O)O)cc1') for m in [mol1, mol2]: m = Chem.AddHs(m) AllChem.EmbedMolecule(m, randomSeed=42) AllChem.MMFFOptimizeMolecule(m) match1 = mol1.GetSubstructMatch(core) match2 = mol2.GetSubstructMatch(core) rmsd = AllChem.AlignMol(mol1, mol2, atomMap=list(zip(match1, match2))) print(rmsd)

这里我要特别强调一个细节：如果公共骨架在分子里出现不止一次（比如分子里有两个苯环），GetSubstructMatch只返回第一组匹配，不一定是最优的。更稳妥的做法是用GetSubstructMatches枚举所有匹配，分别算对齐后的 RMSD，取最小值。

4.3 构象聚类的去冗余思路

做构象搜索后，经常需要对几十甚至几百个构象做去冗余。RMSD 是天然的距离度量，有了构象间的 RMSD 矩阵之后，可以用 RDKit 自带的 Butina 聚类或者直接用层次聚类。Butina 聚类的思路是选择周围邻居最多的构象作为类代表，然后不断移除已聚类的点。这个算法速度很快，而且不要求你预先指定类别数，适合构象去重。

实际操作上，你只需要从每个构象里取坐标，计算所有构象对之间的 RMSD。批量计算 RMSD 可以用AllChem.GetConformerRMSD或者GetBestRMSD。后者会先做对齐再算 RMSD，更适合比较柔性分子的不同构象，因为你去掉整体旋转和平移的影响后，才能真正看出构象本身的差别。

5. 高频报错与排查实录

5.1 SMILES 解析失败

这是最容易遇到也最好解决的问题。MolFromSmiles返回None时，第一反应不应该是怀疑 RDKit 坏了，而是检查 SMILES 本身的合法性。常见原因包括：元素符号写错、价态不满足、方括号里的原子写法有问题。

如果你想看到具体错误信息，可以把sanitize=False传进去，然后手动调用Chem.SanitizeMol，它会告诉你到底是哪一步出了问题，是价态不对、还是芳香性定义有冲突、还是某个原子没有正确闭合。

mol = Chem.MolFromSmiles('C1CCCC1', sanitize=False) try: Chem.SanitizeMol(mol) except Chem.AtomValenceException as e: print(e)

还有一种情况是高通量筛选时遇到大量 SMILES 失败，这时最好的策略是直接跳过并记录，不要试图当场修复每个字符串。等流程跑完后再单独看失败列表，很多问题其实是同一个原因造成的，批量修起来更快。

5.2 3D 嵌入失败

EmbedMolecule偶尔会返回 -1，意思是在距离几何采样阶段没能找到合理构象。常见场景是分子含有苛刻的刚性环系、金属配位环境，或者某些原子的半径约束相互冲突。

解决思路不是死磕一个函数，而是多用EmbedMultipleConfs加randomSeed换种子多试几次。距离几何本质是随机采样，换个种子往往就成功了。如果多试几次还是失败，再考虑把问题分子用 UFF 力场重试优化，或者检查 SMILES 是否包含不常见元素导致距离约束计算失败。

5.3 加氢去氢后索引错乱

这是 2D/3D 协同处理里最让人头疼的问题。加了氢之后原子数量变多，索引变化是必然的。要避免混淆，最关键的原则是：想清楚你到底需要哪个阶段的信息。

如果你需要的是“某个官能团的位置”，就在 2D 阶段先匹配并保存索引；如果你需要的是“氢原子在 3D 空间里的位置”，就等加氢和构象生成之后再去遍历原子。不要在加氢之后还试图用原来的重原子索引去直接访问，也不要反过来。

5.4 力场优化报错

MMFFOptimizeMolecule报错通常跟元素类型有关系，特别是含有硼、硅、某些金属原子时，MMFF94 的参数不齐全。这时候要么改UFFOptimizeMolecule，要么检查分子是否需要先做电荷计算再优化，而不是硬扛。

6. 实用心得：流程设计比单点操作更重要

最后聊一点经验。处理 2D 和 3D 分子，单点 API 用熟只是第一步，真正考验人的是把整个流程串起来。比如一个典型任务：给定一批 SMILES，需要生成 3D 构象、计算指纹相似性、找出与参考分子最接近的候选。你需要明确每一步的数据格式转换：SMILES 转Mol，Mol加氢，加氢后算指纹还是生成构象，生成构象后怎么评估。只要一个环节输出格式不匹配，后面再强也会被卡住。

我自己在项目里常用的顺序是：

1. 先用 2D 操作完成清洗和筛选，比如剔除无效 SMILES、按子结构淘选感兴趣分子；
2. 再对筛选后的分子统一加氢、生成构象、力场优化；
3. 如果需要比较构象，先对齐再算 RMSD，避免直接比较绝对坐标；
4. 所有过程固定随机种子，并记录处理日志，方便追溯。

这个流程看起来很朴素，但真的能救大命。尤其是数据集很大的时候，提前在 2D 层面把不想要的分子筛掉，能省掉后面大量 3D 构象生成和优化的时间。

再提一个小技巧：很多人画 2D 结构图时喜欢直接用默认布局，但 RDKit 的 2D 坐标生成结果跟分子初始表示有关。如果你发现某个分子画出来特别奇怪、键长乱得没法看，可以试试先把分子做一次Chem.MolFromSmiles再转一次 canonical SMILES，很多时候重新读入一次就能拿到更规整的 2D 坐标。

RDKit 这套工具真正强大之处，不在于某个函数有多神奇，而在于它把 2D 的结构逻辑和 3D 的空间逻辑无缝地接在了一起。理解这一点之后，你再遇到“这个分子有没有羧基”“它的构象长什么样”“两个分子空间上像不像”这样的问题，基本都能在几分钟之内写出可复用的代码解决。