基于能量训练和自适应回火的生成式粗粒化方法：一种无数据的多模态玻尔兹曼分布逼近框架

《Journal of Chemical Theory and Computation》：Energy-Based Coarse-Graining in Molecular Dynamics: A Flow-Based Framework without Data

本文提出了一种创新的无数据生成式粗粒化框架，通过结合可逆变换（f?）和归一化流（qθ），直接利用势能函数（U(x)）逼近多模态玻尔兹曼分布（p(x)）。该方法采用反向KL散度（KL(qθ||p?)）作为训练目标，并引入自适应回火方案（adaptive tempering scheme）有效克服了模式锁定（mode-seeking）问题。其核心创新在于将全原子坐标（x）通过双射映射分解为慢变量（z，代表粗粒化自由度）和快变量（X，代表局部涨落），从而实现对复杂能量景观（energy landscape）中所有相关亚稳态（metastable states）的准确捕获。研究通过双阱势、高斯混合模型和丙氨酸二肽（alanine dipeptide）等基准测试表明，该方法能够在不依赖预采样分子动力学（MD）数据的情况下，高效生成热力学一致的全原子构象样本，为分子模拟中的采样和反向映射（back-mapping）问题提供了统一的解决方案。

　　

1. 引言

传统的分子动力学（MD）模拟在探索复杂分子系统的平衡态性质时，常常受限于高维构象空间中的能垒跨越和稀有事件采样。数据驱动的粗粒化（Coarse-Graining, CG）方法通过将全原子细节映射到更低维度的表示来加速采样，但其性能严重依赖于预先生成的高质量全原子数据，这导致了“先有鸡还是先有蛋”的困境。本文提出了一种全新的、完全无数据的生成式粗粒化框架。该框架的核心思想是直接利用系统的势能函数（U(x)）来训练一个生成模型，使其能够近似目标玻尔兹曼分布（p(x) ∝ exp(-βU(x))），而无需任何先验的采样数据。

该方法的关键创新在于引入了一个结构化的潜空间。具体而言，通过一个可逆的双射变换（f_?），将全原子坐标（x）分解为两部分：一部分是低维的粗粒化坐标（z），用于捕捉系统的慢变量和多重亚稳态（即分布的多模态特性）；另一部分是剩余的自由度（X），代表在给定z条件下的局部快速涨落。与此对应，生成模型（q_θ(X, z)）也由两部分构成：一个用于建模粗粒化坐标边缘分布（q_θ(z)）的多模态模型（如归一化流，Normalizing Flow），以及一个给定z条件下描述X分布的、相对简单的单模态条件分布（q_θ(X|z)，例如高斯分布）。模型的训练通过最小化近似分布q_θ与目标玻尔兹曼分布p之间的反向Kullback-Leibler（KL）散度来实现。为了克服反向KL散度固有的模式寻求（mode-seeking）行为，即容易忽略低概率密度区域的问题，本文设计了一种自适应的回火（tempering）方案。该方案从高温（β较小）开始训练，此时目标分布接近均匀分布，易于探索整个构象空间；然后逐步、自适应地降低温度（增加β），确保在每一步，模型都能平滑地跟踪分布的变化，最终准确捕获所有相关的亚稳态。

2. 方法论

2.1. 问题表述与概率框架

目标是在仅给定势能函数U(x)和逆温度β的情况下，学习一个生成模型q_θ(x)，使其尽可能接近真实的玻尔兹曼分布p(x) = exp(-βU(x)) / Z。直接学习高维的q_θ(x)是困难的。本文的解决方案是引入一个可逆的、参数化的双射映射f_?，将原空间x与一个新的坐标空间(X, z)联系起来：x = f_?(X, z)。这里，z是低维的粗粒化坐标（dim(z) << dim(x)），而X包含了剩余的自由度（dim(X) = dim(x) - dim(z)）。根据变量变换定理，在(X, z)空间中的概率密度为：

q_θ(X, z) = q_θ(X|z) q_θ(z) K_?(X, z)

其中K_?(X, z) = |det(?f_?/?(X, z))|是变换的雅可比行列式。相应地，目标分布在(X, z)空间中变为p_?(X, z) = p(f_?(X, z)) K_?(X, z)。模型训练的目标是最小化q_θ(X, z)与p_?(X, z)之间的反向KL散度：

L(θ, ?) = KL(q_θ(X, z) || p_?(X, z)) = β〈U_?(X, z)〉_qθ - 〈log q_θ(X|z)〉_qθ - 〈log q_θ(z)〉_qθ + log Z

由于配分函数Z与参数(θ, ?)无关，在优化中可以忽略。

对于映射f_?，本文考虑了一种特殊的线性形式：x = A_? [z; X]^T，其中A_?是一个右随机矩阵（每行元素之和为1）。这种选择保证了映射对刚体运动的等变性（equivariance）。一个特例是当A_?是置换矩阵时，相当于对原子坐标进行了一种划分。本文的方法通过学习得到这个映射，而不是预先指定。

2.2. 训练框架

对于条件分布q_θ(X|z)，假设其为高斯分布：q_θ(X|z) = N(X | μ_θ(z), diag(σ_θ²(z)))。这允许使用重参数化技巧（reparametrization trick）来采样和计算梯度。对于边缘分布q_θ(z)，使用归一化流模型，通过一系列可逆变换将一个简单的基分布（如高斯分布）转换为复杂的多模态分布。

训练过程通过随机梯度下降（如ADAM优化器）最小化目标函数L(θ, ?)。梯度的计算涉及对势能U和模型参数(θ, ?)的求导，并通过蒙特卡洛采样进行估计。为了缓解反向KL散度的模式寻求问题，采用了自适应回火策略。该策略从初始的高温（β₀ ≈ 0）开始，此时目标分布平坦，易于训练。然后，根据KL散度的相对变化量（δKL_k）自适应地确定下一步的逆温度增量Δβ_k，确保分布变化不会太大，模型能够平稳过渡到更低的温度（更高的β），直至达到目标温度β_target。

2.3. 预测

一旦模型训练完成，可以从中学到的分布中直接生成样本。具体步骤是：1) 从基分布中采样噪声变量(?, ?_X)；2) 通过流模型计算z坐标：z = g_θ(?)；3) 从条件分布中采样X坐标：X ～ q_θ(X|z)；4) 通过逆映射还原到全原子坐标：x = f_?(X, z)。这样得到的样本近似服从目标玻尔兹曼分布。此外，还可以计算粗粒化空间中的自由能面A(z) = -β^-1 log q_θ(z)，用于分析过渡路径和能垒。对于任意物理观测量a(x)，可以使用重要性采样（Importance Sampling）基于学到的模型q_θ进行无偏估计。

2.4. 模型设定

在具体实现中，用于q_θ(z)的归一化流采用了基于单调有理二次样条（monotonic rational-quadratic splines）的耦合层（coupling layers），这些样条具有高度表达力且完全可微可逆。条件分布q_θ(X|z)的均值μ_θ(z)和方差σ_θ²(z)通过一个多层感知机（MLP）进行参数化。线性映射f_?对应的矩阵A_?通过行方向的softmax变换来保证其右随机性。

3. 数值示例

3.1. 双阱势

第一个示例是一个二维双阱势能函数，其形式为U(x) = (1/4)x₁⁴ - 3x₁² + x₁ + (1/2)x₂²。该势能在x₁方向上有两个明显的极小值（亚稳态），而在x₂方向上是简单的谐振子形式。目标是验证方法能否发现并准确表示这两个模式。模型设定dim(z)=1, dim(X)=1。学习到的变换矩阵A_?表明，z坐标主要与x₁相关联，而这正是系统中区分两个模式的慢变量。通过自适应回火训练，模型成功地在所有中间温度以及目标温度下捕获了两个模式的概率分布，生成的样本直方图与参考分布高度一致。

3.2. 高斯混合模型

第二个示例是一个由三个高斯分量混合而成的目标分布，用于测试方法处理更复杂多模态问题的能力。考虑了两种维度情况：低维（dim(x)=4）和高维（dim(x)=20）。在高维情况下，模式寻求问题更为突出。结果表明，采用自适应回火策略后，模型能够成功识别并表示所有三个模式，而未经回火的训练则会将样本锁定在其中一个模式上。学习到的线性变换A_?及其逆矩阵显示出对角占优的结构，说明模型自动地将粗粒化坐标z与原始空间中包含多模态信息的子空间（即x_z）对齐。

3.3. 丙氨酸二肽

第三个也是最重要的示例是丙氨酸二肽分子，这是一个经典的生物分子模型，其构象分布由两个二面角（Φ, Ψ）主导，并在Ramachandran图上呈现出多个明显的亚稳态区域。参考势能由经过训练的图神经网络DimeNet提供。系统包含10个重原子（已去除氢原子），通过固定部分原子的坐标来消除刚体运动，最终有效维度为24。模型设定dim(z)=15, dim(X)=9，即将其粗粒化为5个伪原子和3个伪原子。

训练过程采用了前述的自适应回火方案。结果分析表明：

1. 1.1.
   变换可视化： 学习到的逆变换矩阵A_?^-1显示，粗粒化坐标z主要与分子骨架原子和氧原子相关联，这些原子正是决定二面角（Φ, Ψ）的关键部分。
2. 2.2.
   构象分布： 在不同逆温度下，模型预测的Ramachandran图与参考模拟结果高度吻合，成功再现了所有主要的构象亚稳态（如α_R, β, P_II等）。
3. 3.3.
   结构细节： 生成的样本在键长、键角等局部结构细节的分布上与参考数据非常一致。
4. 4.4.
   物理观测量： 计算了回转半径（radius of gyration）、均方根偏差（RMSD）和能量等物理量的分布，模型预测与参考结果匹配良好。
5. 5.5.
   评分指标： 使用键合分数（bond score，衡量生成结构的化学合理性）和多样性分数（diversity score，衡量生成样本的覆盖度和真实性）进行评估，结果均显示出优异的性能。

4. 结论

本文提出并验证了一种全新的、无数据的生成式粗粒化方法。该方法通过将系统坐标分解为慢变量（z）和快变量（X），并利用可逆变换和归一化流模型，直接基于势能函数训练生成模型，从而避免了传统方法对预采样数据的依赖。引入的自适应回火策略有效克服了反向KL散度训练中的模式锁定问题。在包括双阱势、高斯混合模型和丙氨酸二肽在内的多个基准测试上，该方法均表现出色，能够准确捕获复杂能量景观中的多模态特性，并生成热力学一致的全原子构象样本。这项工作为分子模拟中的高效采样和粗粒化建模提供了一个强大而通用的框架，有望应用于更复杂的生物分子体系。未来的研究方向包括引入等变架构（如SE(3)-等变流）以嵌入物理对称性，以及开发自动确定最优粗粒化维度（dim(z)）的算法。