OpenFE 做 RBFE：先导优化中的自由能计算工作流

OpenFE（Open Free Energy）是开源的自由能计算工作流框架。它把此前需要大量手工配置的 RBFE 计算标准化成可复现的流程——这对让 FEP 从「专家专用」变成「团队可用」很关键。

OpenFE 的组成

完整流程

conda install -c conda-forge openfe
# 或 pip install openfe（依赖较多，建议用 conda）

# ---- 1) 准备输入 ----
import openfe
from openfe import SmallMoleculeComponent, ProteinComponent, SolventComponent
from rdkit import Chem

# 配体：必须是【已对齐到同一结合模式】的三维结构
#   通常来自共晶结构，或以共晶配体为模板的约束嵌入
supplier = Chem.SDMolSupplier("ligands_aligned.sdf", removeHs=False)
ligands = [SmallMoleculeComponent.from_rdkit(m) for m in supplier if m]
print(f"载入 {len(ligands)} 个配体")

protein = ProteinComponent.from_pdb_file("protein_prepared.pdb")
solvent = SolventComponent(ion_concentration=0.15 * openfe.units.molar)

# 【最关键的前提】：
#   所有配体必须以【相同的结合模式】叠合好
#   如果叠合错了，整个计算无意义
#   → 用共晶结构或约束对接生成（见 100 的约束对接）

# ---- 2) 生成原子映射 ----
from openfe.setup import KartografAtomMapper
# 或 LomapAtomMapper

mapper = KartografAtomMapper(atom_map_hydrogens=True)

# ---- 3) 规划微扰网络 ----
from openfe.setup.ligand_network_planning import generate_lomap_network
from openfe.setup.atom_mapping import lomap_scorers

network = generate_lomap_network(
    molecules=ligands,
    mappers=[mapper],
    scorer=lomap_scorers.default_lomap_score,
)

print(f"网络包含 {len(network.edges)} 条边")
for edge in network.edges:
    print(f"  {edge.componentA.name} → {edge.componentB.name}")

# 可视化网络（检查是否有闭合环路）
network.draw_to_file("network.png")

# ---- 4) 创建计算任务 ----
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

settings = RelativeHybridTopologyProtocol.default_settings()
settings.simulation_settings.equilibration_length = 1.0 * openfe.units.nanosecond
settings.simulation_settings.production_length = 5.0 * openfe.units.nanosecond
settings.lambda_settings.lambda_windows = 11
settings.engine_settings.compute_platform = "CUDA"
settings.protocol_repeats = 3        # 【重要】独立重复次数

protocol = RelativeHybridTopologyProtocol(settings)

# ---- 5) 导出任务并执行 ----
from openfe.setup import RBFEAlchemicalNetworkPlanner

# 通常用命令行工具执行：
#   openfe plan-rbfe-network -M ligands.sdf -p protein.pdb -o network/
#   openfe quickrun network/transformations/xxx.json -o results/
#
# 在集群上：每条边一个作业，可完全并行

命令行工作流（推荐）

# OpenFE 的 CLI 更适合生产使用

# 1) 规划网络
openfe plan-rbfe-network \
  -M ligands_aligned.sdf \
  -p protein_prepared.pdb \
  -o network_setup/

# 输出：
#   network_setup/transformations/*.json   每条边一个任务
#   network_setup/ligand_network.graphml   网络结构

# 2) 执行（每条边独立，可并行）
for tf in network_setup/transformations/*.json; do
    name=$(basename "$tf" .json)
    openfe quickrun "$tf" -o "results/${name}.json" -d "work/${name}"
done

# 集群上用作业数组：
#   sbatch --array=1-40 run_edge.sh

# 3) 汇总结果
openfe gather results/ --report dg -o final_results.tsv

# --report 选项：
#   dg      绝对自由能估计（需要至少一个实验锚点）
#   ddg     相对自由能（每条边）
#   raw     原始数据

网络设计的关键考量

# Lomap 评分考虑的因素：
#   - 最大公共子结构的大小（越大越好）
#   - 变化的原子数（越少越好）
#   - 是否有电荷变化（有则扣分）
#   - 是否有环开合（有则扣分 —— 这类变换很难算准）
#
# 手工检查网络的要点：
#
# 1) 【必须有闭合环路】
#    环路上 ΔΔG 之和应为 0
#    → 这是唯一的内部一致性检验
#    → 纯星形网络无法自检，不推荐
#
# 2) 每个分子至少 2 条边
#    单条边连接的分子，其结果无冗余验证
#
# 3) 避免困难的变换
#    - 净电荷变化（需特殊处理）
#    - 环的开合
#    - 大的骨架变化
#    → 如果 Lomap 给出这类边，考虑手工调整
#
# 4) 实验锚点
#    网络中应有 1~2 个已知实验值的分子
#    → 把相对值转换为绝对值
#
# 手工添加/删除边：
from openfe.setup import LigandNetwork

# 检查网络连通性
import networkx as nx
g = network.graph
print("连通:", nx.is_connected(g.to_undirected()))
print("环路数:", len(nx.cycle_basis(g.to_undirected())))
for node in g.nodes:
    print(f"{node.name}: 度数 {g.degree(node)}")

结果分析与质量检查

import json
import numpy as np

# 读取单条边的结果
with open("results/edge_A_B.json") as f:
    res = json.load(f)

# 关键指标：
#   estimate      ΔΔG 估计值（kcal/mol）
#   uncertainty   统计误差
#
# 【质量检查清单】：
#
# 1) 各重复之间的一致性
#    protocol_repeats=3 的三次结果差异
#    → 差异 > 0.5 kcal/mol 说明采样不足
#
# 2) λ 窗口之间的重叠
#    相邻窗口的能量分布应有足够重叠
#    → 重叠不足 = 需要增加窗口数
#    → OpenFE 会输出重叠矩阵
#
# 3) 【环路闭合误差】
#    这是最重要的整体检验
def cycle_closure_error(edges_ddg, cycles):
    """计算每个环路的闭合误差"""
    errors = []
    for cycle in cycles:
        total = sum(edges_ddg[e] for e in cycle)
        errors.append(abs(total))
    return errors

#    判据：
#      < 0.5 kcal/mol   良好
#      0.5 ~ 1.0        可接受
#      > 1.0            该区域不可信
#
# 4) 与实验值的比较（如果有）
#    计算 RMSE、Spearman 相关系数、Kendall tau
#    → RMSE < 1.5 kcal/mol 是可用的水平
#
# 5) 收敛性
#    把模拟分成前后两半，各算一次
#    → 差异大说明未收敛，需延长模拟

成本与实践建议

• 单条边：约 5～20 GPU 小时（3 次重复、11 个 λ 窗口、5 ns/窗口）；
• 20 个分子的网络：约 25～35 条边 → 数百 GPU 小时；
• 先做验证集：用该体系已知活性的 8～10 个分子先跑一遍，看能否复现实验排序。这一步不能省——它告诉你 FEP 在这个体系上能不能用；
• 逐步扩大：验证通过后再用于预测新分子；
• 结合模式是前提：没有可靠的共晶结构就不要做 FEP，先去拿结构；
• 记录所有设置：力场版本、水模型、λ 窗口数、模拟长度——这些都影响结果，必须可复现。

延伸资源

• FEP 原理：127《FEP、RBFE 与 ABFE》；MM/GBSA：126《MM/GBSA》；OpenMM：128《OpenMM 跑 10–100 ns MD》；
• 力场参数化：202《OpenFF Toolkit》；先导优化：405《Lead Optimization》。