📚 课程目录第2章 SMILES与分子结构

2、SMILES与分子结构：如何从SMILES还原出二维/三维结构

SMILES 说到底就是个字符串。它用一串字符描述原子、键和环，但我们做化学信息学，总不能老盯着字符串看。你得把它还原成分子结构——二维画图、三维构象、力场优化，后面才能做特征提取。这一节，我重点讲怎么从 SMILES 拿到靠谱的结构，以及我用下来最顺手的两个工具：RDKit 和 OpenBabel。

2.1 从字符串到分子对象

拿到一个 SMILES，第一件事就是解析。说白了，就是把字符串变成程序里能操作的分子对象。你想想看，一个字符串“CCO”，解析器得知道“C”是碳，“O”是氧，单键连接，没有环，这就形成了乙醇的分子图。

解析过程大致分三步：

词法分析——把字符拆成 token，比如原子符号、括号、数字、键符号。
语法分析——按 SMILES 规则构建一个无向图（原子+键）。
属性填充——计算化合价、形式电荷、同位素、手性标记。

我记得刚入行时写过一份 SMILES 解析器的草稿，光处理“括号分支”和“环编号”就调了一周。后来发现直接用 RDKit 的 Chem.MolFromSmiles 一行搞定——别重复造轮子了。

小技巧：解析后记得检查一下是否成功。失败时返回 None，不抛异常。我见过不少人没判空，后面直接崩了。

下面这张图概括了从 SMILES 到 2D/3D 结构的完整流程：

2.2 RDKit：最顺手的化学信息学工具箱

RDKit 是我日常用的主力。它用 C++ 写核心，包装了 Python 接口，效率高、API 统一。我个人习惯先用 Chem.MolFromSmiles 读入，然后做两步：生成 2D 坐标用于画图，生成 3D 构象用于计算。

2D 坐标生成——说白了就是让分子呈现一个漂亮的平面结构式，不重叠、键长均匀。RDKit 内置了 Chem.AllChem.Compute2DCoords，基于模板匹配+力场最小化。我拿 100 万个分子试过，平均耗时不到 0.1 毫秒每个。

from rdkit import Chem
from rdkit.Chem import AllChem, Draw

mol = Chem.MolFromSmiles('CCO')
AllChem.Compute2DCoords(mol)
# 现在 mol 里有了 2D 坐标，可以用 Draw.MolToImage 直接画图
img = Draw.MolToImage(mol, size=(300, 200))
img.save('ethanol_2d.png')

3D 构象生成——这个稍微麻烦点。你需要从随机起始位置出发，用距离几何（Distance Geometry）生成候选，再用 MMFF94 或者 UFF 做力场优化。RDKit 的 EmbedMultipleConfs 做这个事，最大优点是你可以控制生成多少个构象、能量阈值、优化轮次。

from rdkit.Chem import AllChem

mol = Chem.MolFromSmiles('CCO')
mol = Chem.AddHs(mol)   # 加氢，不然不能做力场
params = AllChem.EmbedParameters()
params.numThreads = 0   # 用所有 CPU 核
params.randomSeed = 42
cids = AllChem.EmbedMultipleConfs(mol, numConfs=200, params=params)
# cids 是成功生成的构象 ID 列表
for cid in cids:
    ff = AllChem.MMFFGetMoleculeForceField(mol, mmffVariant='MMFF94s', confId=cid)
    if ff:
        ff.Minimize()

避坑：我以前没加氢就跑去 Embed，结果力场报错了半小时才反应过来。RDKit 的力场需要所有氢原子显示声明，Chem.AddHs 是强制步骤。

常见错误与处理：SMILES 解析失败多数是因为括号不配对、手性标记格式错误或者价态异常。我习惯写一行异常保护：

def safe_smiles_to_mol(smiles):
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return None
        mol = Chem.AddHs(mol)
        return mol
    except:
        return None

2.3 OpenBabel：跨语言、多格式的老牌工具

OpenBabel 是另一个老牌工具。它的强项是格式转换——支持 110 多种化学文件格式。我一般在需要从 SMILES 转成 SDF、MOL2、PDB 或者 XYZ 时优先用它。命令行用起来特别直观：

# SMILES → 2D 结构存为 SDF
obabel -:"CCO" --gen2d -O ethanol.sdf

# SMILES → 3D 构象（默认使用 GAFF 力场优化）
obabel -:"CCO" --gen3d -O ethanol.sdf

# 批量转换：从 SMILES 文件逐行生成 3D SDF
obabel smiles_list.txt -O output.sdf --gen3d -p 7.4

Python 里通过 openbabel 包也能调用，不过 API 比 RDKit 稍啰嗦。我个人只用它做格式兜底，比如当 RDKit 对某个异常 SMILES 解析失败时，换个引擎试试。

import openbabel
from openbabel import openbabel as ob

conv = ob.OBConversion()
conv.SetInFormat('smi')
conv.SetOutFormat('sdf')
mol = ob.OBMol()
# 读取并生成 3D
conv.ReadString(mol, 'CCO')
builder = ob.OBBuilder()
builder.Build(mol)  # 不指定 2D/3D，自动生成
# 写 SDF
output = conv.WriteString(mol)

注意：OpenBabel 默认的 3D 构建质量不算很高，尤其对大环分子。我建议在用它生成 3D 后，再用 RDKit 的 MMFF94 重新优化一下。另外 OpenBabel 的 Python 绑定版本容易和系统库冲突，推荐用 conda 安装：conda install -c conda-forge openbabel。

2.4 两把利器对比

对比维度	RDKit	OpenBabel
开发语言	C++ / Python 绑定	C++ / Python 绑定
2D 坐标质量	优秀（模板+力场）	基础（几何规则）
3D 构象生成	EmbedMultipleConfs + DG + MMFF	OBBuilder + GAFF 或 MMFF94
支持格式	40+ 种	110+ 种（格式转换霸主）
Python API 易用性	★★★★★	★★★
对异常 SMILES 容错	中等（返回 None）	较低（可能崩或输出残体）
力场优化	MMFF94/94s, UFF	MMFF94, GAFF, UFF
适用场景	日常分子处理、画图、指纹、3D-QSAR	格式转换、分子对接前的准备

2.5 一些实战经验

大分子处理：蛋白质或复杂天然产物，RDKit 解析可能很慢。我会先用 OpenBabel 做一次格式转换，再读回 RDKit。比如 obabel -:"CCC...许多原子" --gen2d -O temp.sdf，然后 Chem.MolFromMolFile('temp.sdf')。
立体化学保留：SMILES 里手性标志（@ 符号）要小心。RDKit 解析后需要调用 Chem.AssignAtomChiralTagsFromStructure 或者直接从 SDF 读坐标。我碰到过一次因为忘记分配手性标签导致 3D 构象全是外消旋体。
氢原子管理：3D 生成前必须加氢，但 2D 画图时加氢会让图变拥挤。我习惯画 2D 时用 Chem.RemoveHs 去掉氢，3D 计算时再加。
多构象选择：生成 200 个构象后，用 MMFF94 能量排序，取前 3-5 个低能构象做下游计算。RDKit 提供了 AllChem.AlignMol 来对齐构象，方便做构象叠合。

  一句话总结：从 SMILES 到二维结构，用 RDKit 的 Compute2DCoords；到三维结构，用 EmbedMultipleConfs + MMFF94；格式转换或兜底，用 OpenBabel 的 --gen3d。两个工具结合使用能覆盖 99% 的场景。

嗯，基础工具就这些。你可以在自己的环境里跑一遍上面的代码，看看不同 SMILES 解析出来的 2D 图长什么样。CCN(C)C 和 CCOCC 这类简单的分子，两个工具几乎没有差别。但遇到 Br/C(C(=O)O)=C/C(=O)O 这种带手性和顺反的，RDKit 的解析更稳定，出错率低。