FE-ToolKit方法学深度解析：推公式和物理意义

第一部分：使用ndfes构建和分析自由能形貌

伞形采样与集合变量（CVs）的原理

谐波偏置势 (Harmonic Biasing Potential) 的定义

• 在 $ndfes$ 中，用于伞形采样的谐波偏置势由以下函数形式定义：

  \[W(\xi) = \sum_{d=1}^{N_{dim}} k_{d} (\xi_{d} - \xi_{0,d})^2\]

  其中，$\xi_{d}$ 是第 $d$ 个集合变量（CV）的值，$\xi_{0,d}$ 是伞形窗（umbrella window）的中心位置，而 $k_{d}$ 是谐波系数。

• 一个必须注意的实践细节是，ndfes 所采用的公式（与 Amber 等主流模拟软件一致）省略了传统物理学中弹簧势能公式前导的 $1/2$ 因子。传统公式通常写作：

  \[W(\xi) = \sum_{d=1}^{N_{dim}} \frac{k_{\text{spr},d}}{2} (\xi_{d} - \xi_{0,d})^2\]

  这意味着在配置输入文件时，用户提供给 $ndfes$ 的谐波系数 $k_{d}$ 应该是传统物理意义上弹簧常数 $k_{\text{spr},d}$ 的一半。这种对实际应用细节的明确，体现了该工具包在设计上的严谨性，旨在帮助用户避免因定义不一致而导致的常见配置错误。

多态贝内特接受率（MBAR）方法

公式的通俗解释

我们的最终目标是得到无偏的自由能 $F_h(\xi)$，它与无偏概率分布 $\rho_h(\xi)$ 的关系由统计力学的基本公式定义：

其中，$k_B$ 是玻尔兹曼常数，$T$ 是温度。对于一个离散的箱（bin）$m$，其概率可以看作是所有样本点落入该箱的加权总和：

这里的 $\delta(\cdot)$ 函数判断样本是否在箱内，关键在于权重 $\omega_h$。这个权重告诉我们，一个在偏置模拟中采样到的点，在真实的、无偏的世界里应该有多“重要”。其公式为：

简单来说，这个权重是一个校正因子。它的物理意义是：将一个在人工偏置（biased）环境下得到的观测样本，其重要性（或贡献度）修正回它在真实物理（unbiased）环境下应有的水平。

为了理解这一点，我们需要从一个更简单的概念“重要性采样”说起，然后将其推广到MBAR的复杂情况。

1. 从一个简单的例子说起：重要性采样 (Importance Sampling)

想象一个思想实验：

• 我们的目标: 测量一个山脉（代表真实的、无偏的能量形貌）中，海拔低于1000米区域的平均温度。
• 我们的工具: 一个有故障的探测器，它更喜欢在海拔高的山峰上着陆采样（代表一个有偏置的模拟），而在山谷里采样很少。

如果我们直接平均所有采集到的温度数据，由于大部分数据来自寒冷的山峰，我们得到的平均温度一定会远低于真实的谷底平均温度。这就是偏置（bias）。

如何校正？对于每一个采集到的数据点，我们需要乘以一个“权重”：

• 如果在山峰（探测器喜欢去的地方）采集到一个数据点，我们需要降低它的权重，因为它被过采样了。
• 如果幸运地在山谷（探测器不喜欢去的地方）采集到一个数据点，我们需要大大增加它的权重，因为它被欠采样了。

这个权重具体应该是多少呢？直观上，它应该是：权重∝我们的探测器实际在该地点采样的概率/一个地点在真实世界中应该被采样的概率。这个比值就是“重要性权重”的核心思想。

变量解释

分子的详细解析

分子项是 $\exp[\beta F_{h} - \beta U_{h}(r_{hkn})]$。要理解它，我们首先需要明确 $F_h$ 的含义。

$F_h$ 是什么？

$F_h$ 是整个无偏物理系统的亥姆霍兹自由能，它是一个描述系统整体热力学性质的常数，与任何一个具体的微观构象无关。它由系统的配分函数 $Z_h$ 决定：

配分函数 $Z_h$ 是对系统所有可能构象的玻尔兹曼因子求和：

什么是玻尔兹曼因子？对于一个特定的微观构象 $r_{hkn}$，其玻尔兹曼因子是 $\exp[-\beta U_{h}(r_{hkn})]$。这个值本身是一个未归一化的、相对的概率。能量越低的构象，其玻尔兹曼因子越大，出现的可能性也越大。

分子究竟是什么？

根据 $F_h$ 的定义，我们可以推导出：

因此，分子可以重写为：

这个完整的项 $\frac{\exp[-\beta U_h(r_{hkn})]}{Z_h}$ 才是该构象 $r_{hkn}$ 在无偏物理系统中出现的真实的、归一化的概率。所以，整个分子项代表的是真实概率，其核心是玻尔兹曼因子。

分母是什么？

通过将分母设置为对所有偏置实验的总和，MBAR确保了信息的最大化利用。任何一个构象，无论它是在哪个窗口被采样的，它的权重都是综合考虑了所有其他窗口的信息后计算出来的。这种“全局视野”使得MBAR在理论上比只考虑相邻窗口重叠的WHAM等方法更为精确和高效。

这是一种将所有信息源进行最佳组合的方式，以推断真实世界的分布。

这个过程的难点在于，权重 $\omega_h$ 的计算依赖于所有偏置态的自由能 $F_{hk’}$，而这些自由能本身又是待求量。因此，这是一个必须通过自洽迭代求解的方程组，MBAR算法的核心就是高效地解决这个问题。

变分自由能剖面（vFEP）方法

与MBAR不同，vFEP是一种参数化方法。它不直接计算离散点的概率，而是假设整个FES可以用一个全局的、连续光滑的解析函数 $F_h(\xi;p)$ 来建模。

vFEP的核心思想：用光滑函数拟合能量地貌

可以这样理解：如果说MBAR是在地图上测量并标注出一系列离散点的海拔高度，那么vFEP则是尝试找到一条单一的、平滑的数学曲线（或曲面）来完美地拟合这整个山脉的地形。vFEP通过最大化观测到所有偏置模拟样本的对数似然函数来找到最优的函数参数 $p$。

这里的核心是，通过调整参数 $p$ 来让我们的模型 $F_h(\xi;p)$ 变得最好。

基函数详解：基数B样条 (Basis Functions Explained: Cardinal B-splines)

为了构建这个全局函数 $F_h(\xi;p)$，vFEP使用了一种强大而灵活的数学工具——B样条。

什么是B样条？

通俗理解：可以把B样条想象成一种标准化的、柔性的“曲线积木”。每一块“积木”（一个B样条基函数）都是一个平滑的、钟形的局部曲线，它只在空间的一小块区域内有值（不为零），在其他地方都为零。

基数（Cardinal）的含义

特指这些“积木”的节点（knots）是等间距分布的。这使得它们特别适合用来描述在规则网格上定义的函数，比如FES计算的自由能表面。

B样条与CVs（距离、角度）的关系

这正是理解vFEP的关键。B样条基函数 $B_i(\xi)$ 的自变量 $\xi$ 就是我们定义的集合变量（CV）。

一维情况（如原子间距离）:

假设我们的CV $\xi$ 是原子A和原子B之间的距离。那么，一个一维的B样条基函数 $B_i(\xi)$ 就是一个关于这个距离的函数。例如，一个基函数 $B_5(\text{distance})$ 可能被中心设置在3.0 Å处，它的形状像一个平滑的“小山包”，只在距离为2.5 Å到3.5 Å的区间内有值，而在其他距离值上都为零。整个一维的FES就是由许多这样的、沿着距离坐标轴排布的“小山包”加权叠加而成。

二维情况（如距离和角度）:

假设我们的CVs是距离 $\xi_1$ 和角度 $\xi_2$。此时，基函数就是一个二维的B样条 $B_{i,j}(\xi_1, \xi_2)$。现在，这个“积木”不再是线上的“小山包”，而是一个在（距离，角度）平面上的“平滑土堆”。例如，一个基函数 $B_{5,10}(\text{distance}, \text{angle})$ 可能中心位于 (3.0 Å, 90°)，并且只在一个小的矩形区域（比如距离在2.5-3.5 Å之间，同时角度在80°-100°之间）有值。

vFEP的应用：构建最终的FES

vFEP的最终目标，就是将整个自由能表面表示为这些B样条“积木”的线性组合：

• $B_i(\xi)$：是我们预先定义好的、固定形状和位置的B样条基函数（“积木”）。
• $c_i$：是待优化的系数，可以理解为每一块“积木”所需要的高度（权重）。

vFEP通过复杂的优化算法，找到一组最优的系数 $c_i$，使得由这些“积木”搭建起来的最终FES曲面，能够最好地拟合我们从分子动力学模拟中得到的所有采样数据。

vFEP vs. MBAR：互补的方法论

工具包中同时包含这两种方法并非冗余，而是反映了数据建模中的一个根本性权衡。

• MBAR：非参数化，更忠实于原始数据，但在数据稀疏的区域，其结果可能不连续或充满噪声。
• vFEP：参数化，假设了FES的平滑性，因此总能给出一个平滑、连续的表面，并且具有解析形式，便于后续分析（如求导寻找MFEP），但可能因函数形式选择不当而引入模型偏见。

提供这两种方法使用户能够根据具体问题选择最合适的工具，或同时使用两者进行交叉验证，这是严谨科学研究的关键实践。

重权熵的解释

重权熵（reweighting entropy）是一个关键的质量评估指标。对于一个给定的箱 $m$，其重权熵 $S_t(\xi_m)$ 定义为： \(S_{t}(\xi_{m}) = -\frac{\sum_{h,k,n} \delta(\xi_{m} - \xi(r_{hkn})) \frac{\omega_{t}(r_{hkn})}{s_{tm}} \ln \frac{\omega_{t}(r_{hkn})}{s_{tm}}}{\ln \left(\sum_{h,k,n} \delta(\xi_{m} - \xi(r_{hkn}))\right)}\)

其中，$s_{tm}$ 是对箱 $m$ 的总权重贡献。这个公式本质上是计算了对箱 $m$ 有贡献的所有样本权重的香农熵，并进行了归一化。

• 直观解释：可以把它理解为对该箱有贡献的偏置模拟（伞形窗）的有效数量。
• 值接近1：表示来自多个不同伞形窗的样本都对这个箱做出了均匀的贡献。这意味着不同伞形窗之间在该区域的相空间重叠良好，重权结果非常可靠。
• 值接近0：表示这个箱的自由能值几乎完全由来自少数几个甚至一个伞形窗的样本决定。这意味着该区域的相空间重叠很差，重权结果可能存在较大偏差，统计上不可靠。

第二部分：使用edgembar进行网络范围的炼金术计算

1. 边目标函数 $F_{(ab)}$ 到底是什么？

边目标函数 $F_{(ab)}$ 是为网络中一条特定的边（即一个从配体 a 到配体 b 的炼金术转换）构建的总目标函数，用于估计该边所有中间 $\lambda$ 状态的自由能。其构建方式如下：

1.1 最小单元：单次重复试验的目标函数 $F_{(ab)est}$

对于一次炼金术模拟（例如，在特定环境 e、特定阶段 s 的一次重复试验 t），我们会在多个 $\lambda$ 插入值下进行采样。MBAR 方法为这一次重复试验构建了一个目标函数 $F_{(ab)est}$。这个函数输入是该次试验中所有 $\lambda$ 状态的待求自由能（向量 $G_{(ab)est}$），输出是一个标量值。当这个函数达到最小值时，对应的输入 $G_{(ab)est}$ 就是 MBAR 方法给出的对这些状态自由能的最佳估计。

公式如下：

1.2 构建完整的边目标函数 $F_{(ab)}$

一条边 (ab) 的计算通常包含多个组成部分：两种环境（如靶标环境和参考环境）、多个阶段（如去电荷、范德华变换、加电荷）和多次重复试验。边目标函数 $F_{(ab)}$ 是将这条边所包含的所有最小单元的目标函数 $F_{(ab)est}$ 加权求和：

这里 $G_{(ab)}$ 是一个包含该边所有环境中、所有阶段、所有重复试验里、所有 $\lambda$ 状态自由能的巨大向量。

2. 核心问题一：目标函数 $F_{(ab)est}$ 的意义是什么？目标是什么？

这个目标函数的意义源于统计学中的最大似然估计（Maximum Likelihood Estimation）原理。

2.1 目标：找到最“可信”的自由能

简单来说，edgembar（以及其底层的 MBAR 方法）的目标是：寻找一组自由能值（$G_{(ab)est}$），使得我们观测到的所有模拟数据出现的总概率最大。

换句话说，我们在问这样一个问题：“假设真实的自由能是这个样子的（由 $G_{(ab)est}$ 描述），那么我能做实验采集到我手上这堆数据的可能性有多大？” 我们要做的就是不断调整对“真实自由能”的猜测（即调整 $G_{(ab)est}$ 的值），直到这个可能性达到最大。

2.2 $F_{(ab)est}$：负对数似然函数

在数学上，直接最大化一个由许多概率连乘的函数（似然函数）很困难，通常我们会转而最小化它的负对数形式。目标函数 $F_{(ab)est}$（如您截图中的 eq. 21）正是这个负对数似然函数。

• 似然函数 (Likelihood): $L(G) \propto P(\text{观测到所有数据} \mid \text{给定自由能 } G)$

• 目标函数 (Objective Function): $F(G) = -\ln(L(G))$

因此，最小化目标函数 $F_{(ab)est}$ 就等价于最大化我们观测到这组模拟数据的概率。这个最小化问题的解，就是 MBAR 方法给出的对各个 $\lambda$ 状态自由能的最优统计估计。

3. 核心问题二：函数的输入是什么？

要构建并求解目标函数 $F_{(ab)est}$，edgembar 需要以下两类直接来自模拟的原始数据：

3.1 势能矩阵 (The Potential Energy Matrix)

这是最关键的输入。对于在某个炼金术状态 $\lambda_j$ 下进行模拟得到的一个构象（样本）$r_{jk}$，我们不仅需要知道它在自身势能函数 $U(r_{jk}; \lambda_j)$ 下的能量，还需要计算它在所有其他炼金术状态 $\lambda_l$ 的势能函数下的能量 $U(r_{jk}; \lambda_l)$。这正是 eq. 21 中 $\ln$ 函数内部求和项 $\exp[-\beta U_{(ab)es}(r_{(ab)estjk}; \lambda_l)]$ 所需要的输入。这个交叉计算的势能矩阵是实现重权和信息组合的基础。

3.2 每个状态的样本数 (Number of Samples per State)

即 eq. 22 中的 $N_{\mathrm{s},(ab)esti}$。这告诉算法每个状态采样了多少数据，用于正确的统计加权。

公式如下：

4. 核心问题三：为什么最小化 $F_{(ab)}$ 得到的是“无约束”相对自由能 $g_{(ab)}$？

“无约束”在这里的含义是，计算这条边 (ab) 的自由能时，完全不考虑网络中任何其他边的信息或限制。例如，对于一个 a→b→c→a 的三角形环路，在计算 a→b 这条边时，完全不考虑 b→c 和 c→a 的存在，也不强制要求三条边的 $\Delta\Delta G$ 之和必须为零（即循环闭合条件）。

结论： 最小化 $F_{(ab)}$ 得到的是无约束自由能 $g_{(ab)}$，因为 $F_{(ab)}$ 的数学构造决定了它是一个只包含单条边信息的“局部”优化问题。

graph LR
    %% 定义节点样式
    classDef inputNode fill:#FEF9E7,stroke:#F39C12,stroke-width:2px,font-family:SimHei
    classDef processNode fill:#EAF2F8,stroke:#5499C7,stroke-width:2px,font-family:SimHei
    classDef outputNode fill:#E8F8F5,stroke:#16A085,stroke-width:2px,font-family:SimHei

    %% 节点定义
    A("<b>第一步：输入</b><br/>提供所有模拟轨迹<br/>计算交叉势能矩阵与各状态样本数")
    B["<b>第二步：构建目标函数</b><br/>edgembar利用输入数据<br/>构建F<sub>(ab)est</sub>函数"]
    C["<b>第三步：最小化</b><br/>程序求解 G<sup>*</sup><sub>(ab)est</sub> = argmin F<sub>(ab)est</sub><br/>得到一组最优的中间状态自由能"]
    D("<b>第四步：计算物理量</b><br/>利用中间态自由能 G<sup>*</sup><br/>计算最终结果，如 ΔG 和 ΔΔG")

    %% 流程连接
    A --> B
    B --> C
    C --> D

    %% 为节点应用样式
    class A inputNode
    class B,C processNode
    class D outputNode

第三部分：辅助工具与实用工作流程 (fetkutils)

3.1 优化模拟效率

核心问题：模型函数 $O(\lambda_i, \lambda_j)$ 的形式是什么？

为了能够预测任意 $\lambda$ 对之间的交换接受率$O(\lambda, \lambda’)$，而不仅仅是“预烧”模拟中实际计算过的那些点，FE-ToolKit 构建了一个连续的、可解析的数学模型。这个模型的构建过程精巧，分为几个步骤。

核心思想与目标

模型的目标是创建一个连续函数 $O(\lambda, \lambda’)$，它必须满足两个基本条件：

1. 再现性：对于所有在“预烧”模拟中已计算过的 $\lambda$ 对 $(\lambda_i, \lambda_j)$，模型预测值必须等于观测到的平均接受率，即 $O(\lambda_i, \lambda_j) = O_{ij}$。
2. 同一性：任何状态与自身的交换接受率必须为 1，即 $O(\lambda, \lambda) = 1$。

简化问题的坐标变换

为了更容易地满足上述的“同一性”条件，程序首先进行了一个坐标变换（平面直角坐标系旋转45度），从 $(\lambda, \lambda’)$ 变换到新的坐标系 $(u, v)$： \(u(\lambda, \lambda') = \frac{\lambda - \lambda'}{\sqrt{2}}\) \(v(\lambda, \lambda') = \frac{\lambda + \lambda'}{\sqrt{2}}\) 在这个新坐标系中，当 $\lambda = \lambda’$ 时，总有 $u = 0$，这使得处理 $O(\lambda, \lambda) = 1$ 这个条件变得非常方便。

模型函数的核心结构：指数衰减

模型函数 $O(\lambda, \lambda’)$ 的核心结构是一个关于状态间“距离” $u$ 的指数衰减函数： \(O(\lambda, \lambda') = \exp\left(\dfrac{-z(u(\lambda, \lambda'), v(\lambda, \lambda'))}{|u(\lambda, \lambda')|}\right)\)

• 这个形式巧妙地保证了当 $\lambda = \lambda’$ 时，$u = 0$，整个指数项为 $e^0 = 1$，自动满足了“同一性”条件。
• 物理本质：交换成功率随自由能差ΔG增大而指数下降
• 函数衰减的快慢，即交换接受率如何随着 $\lambda$ 的差异而降低，则完全由指数函数 $z(u, v)$ 决定。

关键的指数函数 $z(u, v)$：多象限径向基函数

模型最精妙的部分在于 $z(u, v)$ 的构造。它并非一个简单的常数或函数，而是由多象限径向基函数 (multiquadric radial basis function, RBF) 叠加而成： \(z(u, v) = \sum_{ij} w_{ij} \varphi(\sqrt{(u - u_{ij})^2 + (v - v_{ij})^2})\)

• 可以把它看作是在每个原始数据点 $(u_{ij}, v_{ij})$ 上放置一个“影响力锥”，其形状由径向基函数 $\varphi(r) = \sqrt{1 + (\epsilon r)^2}$ 定义。
• 权重 $w_{ij}$：每个“影响力锥”的高度或强度由一个待定的权重 $w_{ij}$ 控制。
• 最终的 $z(u, v)$：在空间中任意一点 $(u, v)$ 的值，就是所有这些“影响力锥”在该点贡献的总和。

参数化：求解权重 $w_{ij}$

最后一步就是确定所有未知的权重 $w_{ij}$。这是通过建立并求解一个线性方程组来实现的： \(\sum_{kl} A_{(ij),(kl)} w_{kl} = z_{ij}\)

• 等式右边的 $z_{ij}$ 是根据“预烧”模拟中观测到的接受率 $O_{ij}$ 反算出来的目标指数值。
• $A_{(ij),(kl)}$ 是一个由径向基函数构成的系数矩阵。

成果：获得可预测任意λ对交换率的连续模型 → 生成最优λ调度表 🚀

🌀 3.2 约束优化：SSC多项式降维策略

痛点分析

在进行自由能计算时，通常需要优化一系列的λ值（λ-schedule），以确保模拟的效率和准确性。然而，标准的优化方法需要调整 $N-2$ 个λ点，这不仅计算量大，而且容易受到预烧模拟噪声的影响，导致过拟合问题。例如，如果预烧模拟的采样不够充分，可能会导致优化后的λ值分布不合理，影响最终的自由能计算结果。

SSC方案核心：以约束换鲁棒性

为了克服这些问题，提出了一种基于 SSC（Smoothstep Softcore）多项式 的优化策略。核心思想 是通过引入约束条件，将优化问题从高维空间（需要调整 $N-2$ 个λ值）降低到低维空间（只需调整1-2个参数），从而显著降低对噪声的敏感性，提高优化的鲁棒性。

数学工具：SSC多项式

SSC多项式是一种特殊的多项式函数，用于生成平滑的λ值分布。常用的SSC多项式包括三阶和五阶两种：

• 三阶多项式：保证在 $\lambda = 0$ 和 $\lambda = 1$ 时，函数值和一阶导数连续，确保端点平滑。
• 五阶多项式：除了上述特性外，还保证了更高阶导数的连续性，使得整个函数更加平滑。

优化实操：从高维到低维

通过引入SSC多项式，可以将复杂的λ值优化问题简化为调整少数几个参数。具体方法如下：

1. 对称调度（1参数优化）

   对于对称的λ值分布，可以使用一个参数 $\alpha$ 来生成整个λ调度表。公式如下： \(S(\lambda;\alpha) = (2-\alpha)S_1(\lambda) + (\alpha-1)S_2(\lambda)\)

   • 参数范围：$\alpha \in [1,2]$
   • 操作：只需调整 $\alpha$，即可自动生成对称的λ序列。这种方法特别适用于对称的自由能变化，能够显著减少优化的复杂度。

2. 非对称调度（2参数优化）

   如果自由能变化是非对称的，可以使用两个参数 $\alpha_0$ 和 $\alpha_1$ 来生成非对称的λ调度表。公式如下： \(S(\lambda;\alpha_0,\alpha_1) = (1-\lambda)S(\lambda;\alpha_0) + \lambda S(\lambda;\alpha_1)\)

   • 参数范围：$\alpha_0, \alpha_1 \in [1,2]$
   • 操作：通过调整 $(\alpha_0, \alpha_1)$，可以适应左右不对称的自由能变化。这种方法虽然比对称调度复杂一些，但仍然比直接优化 $N-2$ 个λ值要简单得多。

总结

通过引入SSC多项式，可以将复杂的λ值优化问题简化为调整少数几个参数，从而显著降低对噪声的敏感性，提高优化的鲁棒性。这种方法不仅适用于对称的自由能变化，还可以通过引入两个参数来适应非对称的变化，具有广泛的应用前景。

🌀 自动化平衡检测：确保数据质量的“守门员”

在分子模拟中，初始阶段的轨迹通常反映了系统从一个非平衡的初始状态逐渐弛豫到热力学平衡的过程。这部分“未平衡”的数据必须被丢弃，否则会严重影响自由能计算的准确性。FE-ToolKit 提供了一个自动化的算法来客观地确定需要丢弃多少数据。

核心思想

该算法的核心思想是，它不依赖于人眼观察，而是通过一个迭代式的、基于多重统计检验的投票系统来判断一个给定的数据段是否稳定。它从假设“0%的数据需要被丢弃”开始，对剩余的“生产区域”（production region）进行检验。如果检验不通过，它会增加需要丢弃的数据比例（例如，增加5%），然后对新的、更短的“生产区域”重复检验，直到该区域通过所有测试。

分析对象：什么数据被检验？

算法并不直接分析原子的三维坐标，而是分析一个能够反映系统能量状态的一维时间序列数据。

• 伞形采样：检验的是偏置势能随时间变化的数据。
• 炼金术模拟：检验的是相邻λ态之间的势能差。为了更加稳健，它会同时分析“前向”的能量差（即从 $\lambda_i$ 采样，在 $\lambda_{i+1}$ 下计算能量）和“后向”的能量差（即从 $\lambda_i$ 采样，在 $\lambda_{i-1}$ 下计算能量）。如果两种分析建议的平衡时间不同，算法会保守地选择更长的那一个。

“三局两胜”的迭代检验流程

对于每一个提议的“生产区域”，算法会执行以下三个统计检验。如果其中任意两个检验失败，则该数据段被判定为“未平衡”，需要增加舍弃的数据量，进入下一轮迭代。

• 检验一：Welch’s t-检验（检验统计显著性）
  • 方法：将提议的“生产区域”数据在时间上平均分成前后两半。Welch’s t-检验被用来判断这两半数据的均值是否存在统计学上的显著差异。
  • 判定：如果t检验给出的p值小于一个预设的阈值（默认为0.05），则认为前后两半的均值有显著不同，暗示数据尚未稳定，该检验失败。
• 检验二：均值差异容忍度（检验化学显著性）
  • 方法：同样将数据分成两半，但这次直接计算两个均值之差的绝对值。
  • 判定：如果这个差值小于一个用户定义的、在化学或物理上有意义的容忍度 $d_{tol}$（默认为0.1 kcal/mol），则认为这个差异是可以接受的，该检验通过。这个检验是t检验的一个重要补充，它防止了因数据方差极小而导致的、统计上显著但物理上无意义的微小差异被误判为未平衡。
• 检验三：线性回归漂移检测（检验系统性趋势）
  • 方法：对整个“生产区域”的数据进行线性回归，得到一个斜率。然后使用Wald卡方检验来判断这个斜率是否在统计上显著不为零。
  • 判定：如果有统计学证据表明斜率不为零（即存在系统性的能量漂移），则认为数据仍处于弛豫过程中，该检验失败。

这个三检验投票系统结合了对统计波动、实际差异幅度和系统性趋势的考量，提供了一种比单一检验更稳健、可重复且客观的方法来截断模拟数据，从而确保用于最终自由能分析的数据质量。