构象生成入门：2D 分子如何变成 3D 构象

对接、三维相似性、MD、三维模型都需要分子的三维构象。而 SMILES 只有二维信息——构象生成的质量，直接决定了这些下游计算的可靠性。

核心挑战

# 【为什么这个问题不简单】：
#
# 1) 【构象空间随可旋转键指数增长】
#    n 个可旋转键，每个 3 个优势构象
#    → 3^n 种组合
#    → 10 个可旋转键 = 约 6 万种
#
# 2) 【最低能量构象未必是结合构象】
#    分子结合到蛋白时，
#    常常采取【比全局最低能量高几 kcal/mol】的构象
#    → 只生成低能构象可能漏掉正确的
#    → 【这是构象生成中最重要的认识】
#
# 3) 【溶液中是构象系综】
#    分子在溶液中不断变换构象
#    → 「一个构象」是简化
#
# 4) 【环的构象】
#    柔性环（如环己烷的椅式/船式）
#    大环的构象空间尤其复杂（见 373）

ETKDG：RDKit 的标准方法

# 【ETKDG = 距离几何 + 实验扭转角知识】
#
# 步骤：
#   1) 【距离几何】：
#      根据键长、键角、环约束，
#      构建原子间距离的上下界矩阵
#   2) 【嵌入】：
#      从距离矩阵生成三维坐标
#   3) 【扭转角修正】：
#      用从晶体结构统计的扭转角偏好，
#      调整生成的构象
#      → 【这是 ETKDG 相对纯距离几何的关键改进】
#   4) 【力场优化】：
#      用 MMFF 或 UFF 优化几何

from rdkit import Chem
from rdkit.Chem import AllChem

def generate_conformers(smiles, n_confs=50, seed=0xf00d,
                        prune_rms=0.5, optimize=True):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    # 【必须先加氢】
    mol = Chem.AddHs(mol)

    params = AllChem.ETKDGv3()
    params.randomSeed = seed              # 【固定种子保证可复现】
    params.pruneRmsThresh = prune_rms     # 【去掉过于相似的构象】
    params.numThreads = 0                 # 用所有核心
    params.useSmallRingTorsions = True    # 小环的扭转角知识
    params.useMacrocycleTorsions = True   # 【大环专用】
    params.enforceChirality = True        # 保持手性

    cids = AllChem.EmbedMultipleConfs(mol, numConfs=n_confs, params=params)
    if not cids:
        print("【构象生成失败】")
        return None

    if optimize:
        # MMFF 优化（比 UFF 好，但不是所有分子都支持）
        results = AllChem.MMFFOptimizeMoleculeConfs(
            mol, maxIters=1000, numThreads=0)
        # results: [(收敛标志, 能量), ...]
        energies = [e for conv, e in results]
    else:
        energies = [None] * len(cids)

    return mol, list(cids), energies

# 【关键参数】：
#
# numConfs（构象数）：
#   经验规则：
#     可旋转键 <= 7  → 50 个
#     8 ~ 12         → 200 个
#     > 12           → 300+ 个
#   【或用公式】：min(300, 显式计算)
#
# pruneRmsThresh：
#   去掉 RMSD 小于阈值的重复构象
#   0.5 Å 是常用值
#   → 【减少冗余，提高有效构象的多样性】
#
# randomSeed：
#   【必须固定】，否则结果不可复现

构象数量的选择

from rdkit.Chem import Descriptors

def suggest_n_confs(mol):
    """根据柔性建议构象数"""
    n_rot = Descriptors.NumRotatableBonds(mol)
    if n_rot <= 3:
        return 20
    elif n_rot <= 7:
        return 50
    elif n_rot <= 12:
        return 200
    else:
        return 300

# 【但更好的做法是检查收敛】：
#   逐步增加构象数，看新增的构象
#   是否还在贡献新的构象空间

import numpy as np
from rdkit.Chem import rdMolAlign

def conformer_coverage(mol, cids, threshold=1.0):
    """评估构象集合的多样性"""
    n = len(cids)
    # 计算两两 RMSD
    rmsds = np.zeros((n, n))
    for i in range(n):
        for j in range(i + 1, n):
            r = rdMolAlign.GetBestRMS(mol, mol, prbId=cids[i], refId=cids[j])
            rmsds[i, j] = rmsds[j, i] = r

    # 【贪心聚类：有多少个「独立」的构象】
    unassigned = set(range(n))
    clusters = []
    while unassigned:
        seed = min(unassigned)
        cluster = {k for k in unassigned if rmsds[seed, k] < threshold}
        clusters.append(cluster)
        unassigned -= cluster

    print(f"{n} 个构象 → {len(clusters)} 个独立构象簇")
    return clusters

# 【判读】：
#   簇数远小于构象数 → 生成的构象冗余，可以减少
#   簇数接近构象数   → 【构象空间可能未充分覆盖，应增加】

能量窗口与构象筛选

# 【常见做法】：只保留能量窗口内的构象
#   如：最低能量 + 10 kcal/mol 以内
#
# 【但要谨慎】：
#   1) 力场能量不准确
#      MMFF 的相对能量误差可达数 kcal/mol
#      → 【能量窗口太窄可能丢掉正确的构象】
#
#   2) 【结合构象常常不是最低能量构象】
#      → 用于对接时，窗口应该宽一些（10~20 kcal/mol）
#
#   3) 真空 vs 溶液
#      MMFF 优化通常在真空中
#      → 分子内氢键被过度稳定
#      → 【溶液中的构象分布可能不同】

def filter_by_energy(mol, cids, energies, window=10.0):
    """按能量窗口筛选构象"""
    if not energies or energies[0] is None:
        return cids
    e_min = min(energies)
    kept = [cid for cid, e in zip(cids, energies)
            if e - e_min <= window]
    print(f"能量窗口 {window} kcal/mol: {len(cids)} → {len(kept)}")
    return kept

# 【实践建议】：
#   用于对接：宽窗口（或不筛），让对接程序自己搜索
#   用于三维相似性：中等窗口
#   用于三维描述符计算：可以只用最低能量构象

特殊情况

• 大环化合物：构象空间复杂，标准 ETKDG 可能效果不佳。必须开启 useMacrocycleTorsions=True，且需要更多构象（见 373《大环化合物 Macrocycle》）；
• 手性中心未指定：会随机分配——应该先枚举所有立体异构体，或明确指定；
• 生成失败：某些高度约束的结构（如稠合的多环）可能嵌入失败。可以尝试 useRandomCoords=True 或增加 maxAttempts；
• MMFF 不支持的原子：某些元素 MMFF 没有参数，会退回 UFF（质量较低）；
• 与实验构象比较：如果有该分子或类似物的晶体结构，应该检查生成的构象集合中是否包含接近的构象——这是最直接的质量验证。

质量验证

# 【最有说服力的验证：能否复现晶体构象】

from rdkit.Chem import rdMolAlign

def validate_against_crystal(smiles, crystal_sdf, n_confs=200):
    """检查生成的构象中是否包含接近晶体的构象"""
    crystal = Chem.RemoveHs(Chem.SDMolSupplier(crystal_sdf)[0])
    mol, cids, energies = generate_conformers(smiles, n_confs=n_confs)
    mol_noh = Chem.RemoveHs(mol)

    best_rmsd, best_cid = float("inf"), None
    for cid in cids:
        try:
            r = rdMolAlign.GetBestRMS(mol_noh, crystal, prbId=cid)
            if r < best_rmsd:
                best_rmsd, best_cid = r, cid
        except Exception:
            continue

    print(f"最接近晶体构象的 RMSD: {best_rmsd:.2f} Å")
    # 【判据】：
    #   < 1.0 Å  很好
    #   1.0~2.0  可接受
    #   > 2.0    【构象生成未覆盖到结合构象】
    #            → 需要增加构象数或换方法

    # 【同时看这个构象的能量排名】
    if best_cid is not None and energies[0] is not None:
        rank = sorted(energies).index(energies[cids.index(best_cid)])
        print(f"该构象的能量排名: 第 {rank+1} / {len(cids)}")
        # 【常常发现结合构象的能量排名并不靠前】
        # → 这印证了「不要只用最低能量构象」

延伸资源

• Omega：279《OpenEye Omega》；分子图：034《分子图表示》；大环：373《大环化合物 Macrocycle》；
• 对接准备：106《Meeko》；Uni-Mol：142《Uni-Mol 论文精读》；等变网络：160《等变神经网络》。