4. 原子特征工程：原子类型、杂化状态、电荷、手性等特征提取与编码

分子里每个原子都是一个“节点”，但节点本身不是只是符号。你得把它的真实属性——比如它是什么元素、杂化方式、带多少电荷、是不是手性中心——变成数字，GNN才能听懂。这个步骤我习惯叫它“原子本地特征工程”。

你想想看，一个碳原子在甲烷里是sp³，在乙烯里是sp²，在乙炔里是sp，这三种碳在分子里的行为完全不同。要是把它们的特征向量设成一样的，GNN学了等于白学。所以我们必须把原子的这些细微差异编码进去。

4.1 原子类型特征：从元素符号到One‑hot编码

最常见的做法是把原子序数或元素符号映射成整数，然后做 One‑hot。但这里有个坑——分子数据集里常见的元素只有十来种（C, N, O, S, F, Cl, Br, I, P, 再加上H），像过渡金属很少出现。

我自己的做法是预先定义一张元素表，只覆盖常见元素，把稀有元素统一归为“其他”。这样编码长度可控，模型不会因为极端罕见元素而浪费参数。

# 常见元素列表（按原子序数排序）
ATOM_TYPES = ['H', 'C', 'N', 'O', 'F', 'P', 'S', 'Cl', 'Br', 'I', 'OTHER']

def one_hot_element(atomic_number):
    symbol = element_symbol(atomic_number)  # 假设已有映射函数
    if symbol not in ATOM_TYPES:
        symbol = 'OTHER'
    idx = ATOM_TYPES.index(symbol)
    vec = [0.0] * len(ATOM_TYPES)
    vec[idx] = 1.0
    return vec

还有一种做法是直接用原子序数归一化：用原子序数除以最大支撑数（比如118），得到一个[0,1]的浮点数。这种连续编码在某些GNN任务里效果也不错，但我个人建议保留One‑hot，因为原子类型本质上是离散范畴，One‑hot给梯度提供的信息更清晰。

4.2 杂化状态：体现原子成键环境

杂化状态直接影响原子在3D空间里的取向和轨道交叠方式。RDKit为每个原子提供了 GetHybridization() 方法，返回值类型是 rdkit.Chem.rdchem.HybridizationType，常见的有：SP、SP2、SP3、SP3D、SP3D2，以及 UNSPECIFIED。

我在项目中遇到过这样的情况：同一个sp²杂化的碳原子，在芳香环中与在烯烃中其实可以通过“是否在环上”进一步区分。但杂化特征本身已经提供了一个很强的先验信号。

杂化枚举值	对应成键特点	编码方式建议
UNSPECIFIED	未定义（例如氢原子）	One‑hot 第1位
SP	直线型，如炔烃	One‑hot 第2位
SP2	平面三角，如烯烃、芳香环	One‑hot 第3位
SP3	四面体，如烷烃、醇	One‑hot 第4位
SP3D	三角双锥，如PCl₅	One‑hot 第5位
SP3D2	八面体，如SF₆	One‑hot 第6位

def encode_hybridization(atom):
    hyb = atom.GetHybridization()
    # 使用RDKit的枚举值 -> 索引映射
    hyb_map = {
        rdchem.HybridizationType.UNSPECIFIED: 0,
        rdchem.HybridizationType.SP: 1,
        rdchem.HybridizationType.SP2: 2,
        rdchem.HybridizationType.SP3: 3,
        rdchem.HybridizationType.SP3D: 4,
        rdchem.HybridizationType.SP3D2: 5,
    }
    idx = hyb_map.get(hyb, 0)  # 默认UNSPECIFIED
    vec = [0.0] * 6
    vec[idx] = 1.0
    return vec

4.3 电荷特征：形式电荷与部分电荷

电荷特征分两种：一种是显式的形式电荷（Formal Charge），大部分RDKit分子中原子带的整数电荷；另一种是计算化学方法算出的部分电荷（如Mulliken电荷、NPA电荷等）。

形式电荷很容易：直接用 atom.GetFormalCharge() 得到一个整数。我的经验是不要把它作为连续值直接喂给网络，因为+1和+2之间不是线性关系。我一般是用One‑hot编码覆盖−1, 0, +1, +2等几个常见值，超出范围的归到“其他”。

部分电荷则需要预先用半经验或DFT方法计算。如果是大型数据集，我倾向于省略部分电荷，只用形式电荷。但在精细化的GNN模型中（例如反应过渡态预测），部分电荷几乎是必选项。

# 形式电荷 One‑hot（支持范围 -2 到 +2）
def encode_formal_charge(atom):
    fc = atom.GetFormalCharge()
    offset = 2  # 索引0对应 -2
    idx = fc + offset
    if idx < 0 or idx > 4:
        # 超出范围，归为OTHER（索引5）
        idx = 5
    vec = [0.0] * 6
    vec[idx] = 1.0
    return vec

4.4 手性特征：原子立体化学编码

手性中心是原子级别的重要属性。一个碳原子连四个不同的基团，它就是手性中心，有R或S构型。在分子图里，手性信息通常存储在 atom.GetChiralTag() 中。

RDKit定义了几种标签：CHI_UNSPECIFIED, CHI_TETRAHEDRAL_CW, CHI_TETRAHEDRAL_CCW。对于其他手性类型（如螺旋性、顺反异构），通常用 atom.GetPropsAsDict() 获取额外属性。

def encode_chiral(atom):
    chiral_tag = atom.GetChiralTag()
    if chiral_tag == rdchem.ChiralType.CHI_TETRAHEDRAL_CW:
        return [1.0, 0.0, 0.0]  # R
    elif chiral_tag == rdchem.ChiralType.CHI_TETRAHEDRAL_CCW:
        return [0.0, 1.0, 0.0]  # S
    else:
        return [0.0, 0.0, 1.0]  # 非手性

注意：不是所有sp³碳都是手性中心。我曾犯过一个错——看到分子中有四个不同基团的碳就直接标手性，结果没考虑该碳在对称环境中（比如内消旋化合物）实际上是非手性的。RDKit的ChiralTag已经帮你处理了这些情况，直接用它就好。

4.5 综合特征向量构建

上面几种特征可以拼接成一个稠密向量。我常用的做法是先确定总维度，然后按顺序填充。比如：元素类型(11维) + 杂化(6维) + 形式电荷(6维) + 手性(3维) = 26维。你也可以加入原子质量、范德华半径、在环中/环大小等特征。

def build_atom_feature(atom):
    feat = []
    feat.extend(one_hot_element(atom.GetAtomicNum()))      # 11
    feat.extend(encode_hybridization(atom))                # 6
    feat.extend(encode_formal_charge(atom))                # 6
    feat.extend(encode_chiral(atom))                       # 3
    # 可选：原子次数、在环中指示等
    feat.append(float(1.0 if atom.IsInRing() else 0.0))   # 1
    return feat  # 总共 27 维

我习惯把所有这些特征的值域限定在[0,1]之间，避免大幅度的数值差异。像One‑hot编码天然就是0/1，而连续值（如原子质量）需要额外归一化。

4.6 一个完整流程：从SMILES到原子特征矩阵

这里是一个完整的处理函数，你可以直接放进你的数据预处理管线里。

from rdkit import Chem
from rdkit.Chem import rdchem, AllChem

def smiles_to_atom_features(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    mol = Chem.AddHs(mol)          # 显式加氢
    AllChem.ComputeGasteigerCharges(mol)
    Chem.SanitizeMol(mol)          # 确保杂化、手性等被计算
    
    atom_feat_list = []
    for atom in mol.GetAtoms():
        feat = build_atom_feature(atom)
        atom_feat_list.append(feat)
    return atom_feat_list  # shape: (num_atoms, feat_dim)

在我的实际项目中，这类特征提取函数通常被封装成 DataLoader 的一部分，在读取SMILES时即时计算，不预先存储到磁盘——因为特征维度固定，而且计算量不大，反而比加载大文件快。

知识结构总览

下面这张图总结了原子特征工程的整体逻辑，从原始SMILES到多维特征向量的处理流程。

图中的四个分支分别对应原子类型、杂化状态、电荷和手性，它们被独立编码后再拼接成固定维度的特征向量。我个人建议最后加上一个“在环中”的二值位，这个信息对芳香性识别和键角预测特别有用，而且成本几乎为零。

4.7 经验总结与避坑

我从几个实际项目里总结了一些教训，这里列出来供你参考：

• 特征维度不宜过大。有的框架把原子特征做到50+维，但数据量若不够大，GNN很容易过拟合。我常用的上限是32维。
• 千万不要忘记H原子。很多SMILES（如“CCO”）不写H，但分子图里每个C上是有H的。你要么用 Chem.AddHs(mol) 显式添加，要么在处理特征时把H的信息隐式包含在重原子的度里。我问过一些同行，他们一开始没加H，模型收敛很慢。
• 手性标签依赖于预处理。如果SMILES来自不同来源（比如ChEMBL、PubChem），有些分子本身的立体化学可能缺失或错误。我曾经遇到一个分子，RDKit标记为CHI_UNSPECIFIED，但实际是有手性的。这时候需要手动检查立体键是否完整。
• 部分电荷计算要配置文件缓存。如果你用Gasteiger电荷，RDKit算得比较快；如果用Mulliken电荷（需要额外调用Open Babel或Pyscf），就会慢很多。我建议只在验证集上用小分子子集试一下效果再决定是否加入。