系综平均与时间平均：分子动力学、集平均与误差棒实操

1. 别急着按计算器：系综平均到底在平均什么

前阵子有个做分子模拟的朋友来找我，说他跑了三条轨迹，同一个体系、同一套力场，算出来的扩散系数差了百分之三十，问我是不是机器出了问题。我看了眼他的输出文件就明白了——他把三个不同初始条件、不同随机种子的结果直接做了算术平均，然后把那个数当成系综平均（Ensemble Average，也叫集平均）报了出去。这个做法在某些条件下是对的，在另一些条件下会错得离谱。系综平均本质上是对一个随机变量在某个概率分布下的期望值，而不是“把几组数加起来除以几”，这两件事看着像，实际差着十万八千里。

这篇东西我想把这件事掰开揉碎讲清楚：系综平均的物理直觉是什么，它在不同领域里长什么样，落到代码和参数上该怎么算，算完了怎么给出一个可信的误差棒。写这个的动机很直接，我自己在这个概念上翻过车，也见过太多人把“样本平均”和“系综平均”混着用。文章适合做分子动力学、蒙特卡洛采样、随机信号处理，以及在机器学习里靠多次独立实验取平均的读者；也适合只是被这个词卡住、想弄明白它和时间平均区别的人。不需要多深的统计底子，跟着走就行。

1.1 系综不是“很多条轨迹”，而是一个概率分布

先把最容易混淆的地方挑明。系综这个词是吉布斯提出来的，他的想法很巧妙：一个宏观体系在给定约束下（比如固定体积和能量，或者固定温度和压强），微观状态有无数种可能。我们没法也不需要跟踪其中某一个具体状态，而是想象有无穷多个思维中的副本系统，它们宏观条件完全一样，微观构型各不相同，散布在整个相空间里。这无穷多个副本构成的集合，就是系综。

关键在于，这些副本不是随便撒的，它们按照一个确定的概率密度 ρ(Γ) 分布，Γ 代表相空间里的一个点。系综平均的正式定义就是一个积分：

⟨A⟩ = ∫ A(Γ) ρ(Γ) dΓ

这里的 A(Γ) 是你关心的观测量，ρ(Γ) 是那个分布。所以系综平均是分布上的期望。你实际拿到手的是有限个样本，用样本平均去估计这个期望——这是估计，不是定义本身。这个区分很重要，因为它决定了后面所有的误差分析逻辑：样本平均和真值之间永远有偏差，偏差有多大、怎么量化，才是实操里的核心问题。

举个生活化的类比。你想知道全国成年人的平均身高。真值是那个“系综平均”，但你没可能量遍所有人。你只能抽样，抽若干个城市、若干个年龄段，然后算样本均值。样本均值能不能代表真值，取决于你的抽样方式是否对应那个真实的概率分布。如果你只在篮球队门口抽样，哪怕抽十万人，估计也是有偏的。分子模拟里“只跑一条轨迹、只取一段构型”就是典型的偏差抽样。

1.2 时间平均和系综平均：两把尺子量同一件事

如果只跑一条轨迹，我们能算的是时间平均：

Ā = lim(T→∞) (1/T) ∫₀ᵀ A(t) dt

它是对一条具体轨迹上的时间序列求平均。而系综平均是对无穷多个副本在同一时刻求平均。这两把尺子量的是同一个物理量，但它们相等是有条件的——这个条件叫遍历性（ergodicity）。直观地说，遍历性意味着一条足够长的轨迹最终会走遍所有能量允许的微观状态，并且在每个状态附近停留的时间比例恰好等于那个状态的系综概率。满足遍历性，时间平均就等于系综平均，你跑一条长轨迹就够了。

但现实往往是：遍历性假设在某些体系上直接失效。最经典的例子是低温下的双势阱体系或者玻璃态体系。温度低，轨迹可能长时间困在一个势阱里，看起来能量、温度都很“平稳”，像是平衡了，但它的时间平均只反映了那一个势阱的统计，完全看不到另一个势阱的贡献。这时候你跑再久也没用，十分钟和一小时算出来的结果几乎一样，因为轨迹根本没跳出去。解决办法有两个：一是加温、加偏置势、做副本交换来帮助跨越势垒；二就是老老实实跑多条独立轨迹，用跨轨迹的系综平均来替代时间平均。

我在实践中通常的做法是混合使用：跑若干条独立轨迹（不同初始速度种子、不同初始构型），每条轨迹先各自平衡、各自采样，得到每条轨迹的时间平均，再把这若干条轨迹的时间平均当作同一分布下的独立样本，做跨轨迹平均。这样既利用了时间平均的采样效率，又用独立轨迹的数量来弥补单条轨迹遍历性不足的问题。至于跑几条合适，后面第 4 节会具体说，一般 5 到 20 条是个比较务实的区间。

2. 系综家族选型：先定哪几个量守恒，再谈怎么平均

系综平均不是凭空做的，你得先选一个系综。选系综这件事，本质上是在回答一个问题：这个体系在实际场景下，哪几个宏观量是固定的、哪几个是自由涨落的。选错了，不只是数值差一点，有时候连物理结论都会反过来。这一节把常见系综的适用边界讲清楚，顺便说说为什么“选系综”这一步决定了你后面能不能从涨落里提取信息。

2.1 四种常用系综及其适用场景

下面这张表是我自己整理时常看的版本，把约束条件、典型用途和常见坑标出来：

看这张表的时候，我建议先问自己一句：我关心的这个量，在实验上是在什么条件下测的？实验在常压下测密度，你就得用 NPT；实验在真空里测团簇振动，那 NVE 或 NVT 更自然。模拟条件和实验条件的对应关系没理清，后面的平均做得再精细也是白搭。

还有一个细节值得强调：系综平均不只是均值。方差、涨落同样是宝贵信息。在 NVT 系综下，定容热容可以直接从能量涨落里拿到：

C_v = (⟨E²⟩ − ⟨E⟩²) / (k_B T²)

这就是涨落-耗散定理的一个具体体现。注意它要求你用的是能产生正确正则分布的温控器，并且采样足够长、涨落被正确保留。如果你用的是 Berendsen 温控器，它会把能量涨落人为压小，这个公式算出来的热容会系统性偏低。这是个非常隐蔽的坑，很多人算完热容发现和文献对不上，查了半天才发现问题出在温控器上。

2.2 选系综时的三个判断顺序

我一般按这个顺序来判断，能避免大部分返工。

第一，看实验条件。常压常温的实验过程，模拟里优先 NPT。超高真空、孤立体系，优先 NVE。有物质交换的过程（吸附、渗透），才考虑巨正则。

第二，看你要算什么量。如果目标量是均值类的（密度、结构因子、扩散系数），NVT 和 NPT 都能给，差异主要体现在体积是否涨落上；如果目标量是涨落类的（热容、压缩系数、介电常数），那对系综和温控器的要求就严格得多，必须保证涨落统计是物理正确的。

第三，看计算成本。NPT 比 NVT 贵，巨正则蒙特卡洛又比分子动力学在某些场景下贵得多，因为粒子插入删除的接受率经常低得让人抓狂。在能回答问题的前提下选最省的那个，这是工程思维，不是偷懒。

提示：换系综之后，一定要重新平衡。从 NVT 切到 NPT，盒子会有一个体积弛豫过程，这段时间的数据不能拿去平均。

3. 分子动力学实操：把系综平均真正算出来

这一节进入具体操作。我拿分子动力学举例，因为它是系综平均这个概念最“落地”的场景之一，参数多、坑也多，讲透了其他领域可以类推。整条链路的顺序是：建模、能量最小化、平衡、采样、后处理。真正决定系综平均质量的，是平衡和采样这两步。

3.1 先判断平衡，再谈采样

新手最常犯的错误是：跑完模拟直接从头平均到尾，把前期的弛豫过程也算进去了。体系从初始构型出发，能量、温度、密度都在漂移，这段数据根本不属于平衡分布的样本。把它们混进去，均值直接被拉偏。

判断平衡的常用手段有这么几个，我一般同时看：

• 累计平均值曲线。把某个量（比如势能）的累计平均画出来，横轴时间、纵轴累计均值。如果曲线在前期剧烈变化后逐渐走平，说明进入了平衡区。走平的那一段起点，就是采样的起点。
• 时间序列的滑动平均。看温度和密度有没有系统性漂移。只有随机涨落、没有趋势，才算稳。
• 结构指标。径向分布函数、均方根偏差这类量，如果它们在两条不同初始构型的轨迹里收敛到同一条曲线，说明采样比较充分了。

具体时长没有通用答案。小分子液体在室温下，通常 100 ps 到 500 ps 的平衡就够了；蛋白质这类大体系可能要几纳秒甚至更长。判断依据永远是数据本身，不是文献里抄来的数字。我见过有人拿别人论文里的“平衡 1 ns”当圣旨，结果自己的体系大了十倍，1 ns 连局部结构都没松弛开。

3.2 温控器和压控器的选择与参数设置

系综平均的物理正确性，很大程度上取决于温控器。这里给一个实用的选型判断：

Berendsen 温控器：弛豫快、鲁棒，适合做初步平衡。但它不产生正确的正则分布，会系统性压制涨落。所以平衡阶段可以用，采样阶段绝对不能用。如果你要算热容、压缩系数这类依赖涨落的量，用它的结果会偏小一大截。

Nose-Hoover 温控器：能给出正确的正则系综分布，是做采样时的标准选择。麻烦在于它的时间常数 τ_T 需要调。τ_T 太小，温度会剧烈振荡，甚至出现能量不守恒的伪影；τ_T 太大，温度和体系耦合不足，采样效率低。经验值是 τ_T 取 100 到 1000 倍的积分步长。如果步长 dt = 1 fs，那 τ_T 大概在 0.1 到 1 ps。对水这类体系，我一般取 0.5 到 2 ps，实测下来比较稳。

Langevin 温控器：加了摩擦项和随机力，也能给出正确分布，摩擦系数 γ 一般取 1 到 5 ps⁻¹。它的好处是对大体系收敛稳，坏处是摩擦会略微影响动力学性质，比如扩散系数。如果算动力学量，γ 要取小一点，或者干脆做多组 γ 然后外推到零。

压控器这边，Berendsen 压控器同样只适合平衡；采样阶段推荐 Parrinello-Rahman，时间常数 τ_P 一般取 1 到 5 ps，同时要注意它对盒子形状的处理——各向同性体系用 iso，膜体系或者晶体用 semiisotropic 或 anisotropic，选错会让盒子被压成奇怪形状。

给一段实际用的配置示意，下面是 LAMMPS 的写法：

# 积分步长 1 fs timestep 0.001 # 平衡阶段：Berendsen，快但只用于弛豫 fix eq_t all temp/berendsen 300.0 300.0 0.1 fix eq_p all press/berendsen iso 1.0 1.0 1.0 run 500000 # 500 ps 平衡 # 采样阶段：Nose-Hoover + Parrinello-Rahman unfix eq_t unfix eq_p fix prod_t all temp/nose-hoover 300.0 300.0 0.5 fix prod_p all press/parrinello-rahman 1.0 1.0 5.0 run 2000000 # 2 ns 采样

注意：从 Berendsen 切到 Nose-Hoover 之后，最好再空跑一段不采样的过渡期，让温控器的额外自由度自己弛豫到稳态。这个细节很多教程都不提，但不做的话，前一百皮秒的数据会带上切换瞬间的伪影。

3.3 从轨迹到平均值：后处理脚本怎么写

采样跑完，你会得到一堆能量文件、轨迹文件。以 GROMACS 为例，提取某个量的时间序列：

# 提取势能时间序列 gmx energy -f md.edr -o potential.xvg # 提取均方根偏差 gmx rms -s topol.tpr -f traj.xtc -o rmsd.xvg -tu ns

拿到 xvg 之后，用 Python 做平均和误差分析。下面这段代码我改过很多次，现在算是个趁手的工具：

import numpy as np def load_xvg(path): """读取 GROMACS xvg 文件，跳过注释和标题行""" data = [] with open(path) as f: for line in f: if line.startswith(('#', '@')): continue parts = line.split() if len(parts) >= 2: data.append((float(parts[0]), float(parts[1]))) return np.array(data) def block_error(x, n_blocks=10): """块平均法估计标准误""" x = np.asarray(x, dtype=float) N = len(x) L = N // n_blocks blocks = x[:L * n_blocks].reshape(n_blocks, L) bm = blocks.mean(axis=1) # 块均值近似独立，标准差除以 sqrt(块数) return bm.std(ddof=1) / np.sqrt(n_blocks) t, e = load_xvg('potential.xvg').T # 丢掉前 20% 作为平衡段 n_cut = int(0.2 * len(e)) e_prod = e[n_cut:] mean_e = e_prod.mean() err_e = block_error(e_prod, n_blocks=10) print(f"势能系综平均估计 = {mean_e:.3f} ± {err_e:.3f} kJ/mol")

这里的关键点在于：丢掉平衡段再算平均，以及误差用块平均而不是朴素标准差除以根号 N。为什么不能用朴素公式？因为时间序列的相邻帧高度相关，它们不算独立样本。直接套 σ/√N 会严重低估误差，我见过低估三到五倍的例子。下一节细讲。

4. 误差棒才是灵魂：自相关、块平均与有效样本数

一个没有误差棒的系综平均值，基本没有参考价值。原因很简单：你报出的数字是估计值，它和真值之间差多少，你必须有办法说清楚。这一节讲三个工具：自相关时间、块平均法、有效样本数。它们解决的是同一个问题——时间序列里的样本不独立。

4.1 自相关时间怎么估

先建立直觉。你每隔 1 ps 记录一次能量，但能量的实际“记忆长度”可能是 10 ps。也就是说，第 1 个点和第 11 个点的相关性已经衰减得差不多了，而第 1 个点和第 2 个点几乎完全相关。那这 1000 个数据点里，真正独立的样本可能只有 100 个左右。样本量虚高，误差自然被低估。

量化这个记忆长度的是自相关函数：

C(t) = ⟨δA(0) δA(t)⟩ / ⟨δA²⟩

其中 δA = A − ⟨A⟩。C(0) = 1，随着 t 增大逐渐衰减。定义积分自相关时间为：

τ_int = ∫₀^∞ C(t) dt

实际算的时候，用离散求和，并且积分到 C(t) 第一次过零为止，避免尾部噪声的贡献：

def integrated_autocorr_time(x, dt=1.0): x = np.asarray(x, float) x = x - x.mean() n = len(x) # 用 FFT 算自相关，比双重循环快得多 f = np.fft.rfft(x, n=2 * n) acf = np.fft.irfft(f * np.conj(f))[:n] acf /= acf[0] # 找第一次过零的位置 zero_cross = np.where(acf < 0)[0] cut = zero_cross[0] if len(zero_cross) else n tau_int = dt * (0.5 + acf[1:cut].sum()) return tau_int, acf

拿到 τ_int 之后，统计效率因子g = 1 + 2τ_int/dt，有效样本数就是：

N_eff = N / g

如果 τ_int = 10 ps、dt = 1 ps，那么 g = 21，一千个数据点的有效样本只有不到 50 个。这个数字能立刻解释为什么有些人明明跑了几百万步，误差还那么大——采样量看着多，独立信息量很少。

4.2 块平均法：不用算自相关也能给误差

自相关时间算起来要挑积分截止点，有点主观。块平均法是另一条路，更省事，也更鲁棒。

思路很直接：把长度 N 的时间序列切成 M 个长度为 L 的连续块，每块算一个均值。只要 L 远大于相关时间 τ_int，这些块均值之间就近似独立了。于是总均值的标准误可以写成：

误差 = std(块均值) / √M

实操上最关键的判断是块长要足够大。我的做法是：让块长从很小开始逐步增加，画一条“块长 vs 误差估计”的曲线。误差估计会先上升，然后趋于一个平台。平台出现的位置就说明块长已经够大，取平台区的数值作为最终误差。用 4.3 的代码可以一次把这条曲线扫出来：

def block_error_curve(x, max_blocks=64): x = np.asarray(x, float) N = len(x) out = [] m = 2 while m <= max_blocks: L = N // m if L < 2: break bm = x[:L * m].reshape(m, L).mean(axis=1) out.append((L, bm.std(ddof=1) / np.sqrt(m))) m *= 2 return out for L, err in block_error_curve(e_prod): print(f"块长 {L:6d} 误差估计 {err:.4f}")

经验上，块长取 5 到 10 倍的 τ_int 保险。如果你实在懒得估 τ_int，就让块数控制在 10 到 20 之间，这是个折中，多数情况下够用。但如果你要发表的数字对精度敏感，还是老老实实扫一遍曲线。

4.3 独立轨迹：系综平均的“正牌”做法

前面讲的都是单条轨迹内部的时间平均分析。真正意义上的系综平均，是跑多条独立轨迹，每条轨迹给一个估计值，然后跨轨迹平均。这是我认为最稳妥的方案，尤其在遍历性存疑的体系里。

具体操作：用不同的初始速度种子（以及不同初始构型，如果能做到的话）跑 M 条轨迹，每条各自平衡、各自采样，得到 M 个估计值 a₁, a₂, ..., a_M。系综平均估计就是：

⟨A⟩ ≈ (1/M) Σ aᵢ

误差用轨迹之间的标准差：

误差 = std(aᵢ, ddof=1) / √M

这个 M 不需要很大，5 到 10 条往往就能把误差压到一个可接受的水平。因为误差随 √M 下降，从 1 条到 5 条能降一半多，从 5 条到 20 条只再降一半。性价比最高的区间就在 5 到 10。

提示：多条轨迹的初始构型最好真的不一样。如果都从同一个平衡构型出发、只改速度种子，轨迹之间的独立性会打折扣，尤其在短时间尺度上。

我在一个扩散系数项目里做过对比：单条 10 ns 轨迹，块平均给出的误差是 8%；换成 5 条 2 ns 的独立轨迹，跨轨迹误差降到了 4% 左右，而且结果和实验值的吻合度明显更好。总计算量一样，结论质量差一倍。

5. 集平均在信号处理里的另一种面孔

系综平均不只出现在物理模拟里。信号处理里的“集平均”，思路几乎一模一样，只是换了一套语言。理解这个对应关系，对跨领域工作的人特别有用，因为你会发现很多看起来不相关的技术，底层是同一个东西。

5.1 随机过程视角：固定时刻的期望

把信号 X(t) 看成一个随机过程，也就是说，对每一个固定的时刻 t，X(t) 本身是一个随机变量。那么集平均就是：

μ_X(t) = E[X(t)]

注意它一般依赖 t。如果这个过程是平稳的，μ_X(t) 变成常数，和 t 无关。平稳加上各态历经，时间平均才等于集平均。这套逻辑和分子模拟里的遍历性假设是完全对应的：分子模拟里的“一条长轨迹”对应时间平均，“多条轨迹”对应集平均。

工程上最经典的集平均应用是诱发响应提取，比如脑电的诱发电位、心电的锁时平均、雷达的脉冲积累。原理很简单：你对系统施加同一个刺激，重复 N 次，记录 N 段信号。真正的响应在不同次试验里是锁时、一致的；而噪声是零均值、与刺激不相关的随机的。把 N 段记录在时间上对齐后逐点平均，噪声被压低，信号被保留。

信噪比的改善是可以算出来的。假设单次记录里信号幅度为 S，噪声标准差为 σ。平均 N 次后，信号仍是 S（因为锁时叠加），而噪声的标准差变成 σ/√N（因为独立随机量求平均）。所以信噪比提升 √N 倍。想提升一倍信噪比，得采集四倍的试验次数。这个 √N 关系是所有集平均方法的共同特征，也是它最让人又爱又恨的地方——提升总是越来越贵。

5.2 实测数据里的集平均怎么做才对

理论上的 √N 很漂亮，实测里的坑主要在“对齐”和“样本筛选”上。

第一步是时间对齐。以诱发响应为例，你必须有一个精确的触发时间戳，所有试验的片段都从这个时间点开始截取。触发检测有抖动，抖动哪怕只有几毫秒，在高频成分上也会造成严重的相位错位，平均之后信号被自己抵消掉。我见过这样的案例：数据处理流程完全正确，就是触发对齐用了错误的通道，结果平均出来的波形平坦得像白噪声。

第二步是伪迹剔除。不是所有试验段都干净。眨眼、运动、电极接触不良都会引入远超正常幅度的干扰。这些段如果混进去平均，会把结果往一边拽。通常按幅度阈值或方差阈值剔掉 10% 到 30% 的试验段，剩下的再平均。剔除标准要事先定好，不能看着结果不满意再回头调阈值，那就成了数据操纵。

第三步是加权平均。不同试验段的噪声水平可能不一样，此时用等权平均不是最优的，应该用方差倒数加权：

import numpy as np def weighted_ensemble_average(segments): """ segments: shape (n_trials, n_samples) 每段先用高频段估计噪声方差，再做方差倒数加权平均 """ segs = np.asarray(segments, float) # 用每段自身的高频能量作为噪声方差的粗估计 var = segs.var(axis=1) w = 1.0 / var w = w / w.sum() avg = np.tensordot(w, segs, axes=(0, 0)) return avg, w

这个加权平均在噪声水平不均匀时，等效样本数比等权平均高，误差能小一截。代价是要先估方差，估得不准反而会引入偏差，所以方差估计的稳定性要检查一下。

注意：集平均只能压制与刺激不相关的噪声。如果噪声本身也是锁时的（比如工频干扰恰好和刺激同步），平均非但压不掉，还会被强化。这类周期性干扰要在预处理阶段用陷波滤波去掉。

6. 常见问题与排查速查

前面讲的都是正向流程，这一节反着来，把我在实际操作中反复遇到的问题列成表，方便你对照症状找原因。

6.1 症状、原因、处理对照表

看这张表的时候，我想强调一点：先怀疑采样，再怀疑代码。我自己的经验是，八成以上的异常结果来自采样不足或系综选错，只有不到两成是脚本写错了。很多人一发现结果不对就去翻代码，翻半天没找出问题，其实问题在物理层面。

6.2 几条我在实操中踩出来的经验

第一条，永远先画累计平均曲线。这是成本最低的平衡诊断。花两分钟画一张图，可能省你两天返工。我现在已经养成习惯，任何模拟跑完第一件事就是画这条曲线，看走平点在哪。

第二条，报数必须带误差。不管是给同事看还是给自己存档，均值后面一定跟一个误差估计。哪怕只是块平均的粗略结果，也比裸数据强。裸数据会给人虚假的精确感，久了会养成坏习惯。

第三条，单位换算核对三遍。能量单位在 kJ/mol 和 kcal/mol 之间差 4.184 倍，这个坑我踩过。一次和实验组对数据，差了将近四倍，两边都以为对方算错了，最后发现是单位问题。现在我脚本里第一件事就是把单位统一，并在输出里显式标出来。

第四条，独立轨迹的初始条件要真独立。我早些年跑多轨迹，图省事只改随机种子，结果前几百皮秒的轨迹相关性很高，跨轨迹误差被低估。后来改成用不同的退火构型作为起点，误差估计立刻变得可靠了。

第五条，时间对齐的精度决定集平均的成败。信号处理里这条几乎是铁律。对齐误差要是超过信号最高频成分周期的十分之一，平均就是在自毁信号。有条件的话用互相关做亚采样对齐，效果比取整对齐好很多。

6.3 后续可以往哪扩展

如果你已经把这套流程跑通了，往下有几个方向值得试试。一个是多变量系综平均，同时平均多个关联量，然后看它们的协方差矩阵，能提取更多物理信息，比如耦合涨落。另一个是重加权方法，比如直方图重加权和各类自由能重加权技术，它们能在不重新模拟的前提下，把某个条件下的系综平均“搬”到另一个条件去，前提是样本重叠足够。还有自适应采样这类做法，让采样资源自动往需要的地方倾斜，对大体系和稀有事件特别有用。

不过这些扩展都有一个共同前提：基础的系综平均和误差估计你得做扎实。基础不牢，重加权出来的结果看着很美，实际站不住脚。我个人的建议是，先用一个简单体系（比如液态氩或者纯水）把整条链路跑顺，把平衡判断、块平均、独立轨迹这套流程都吃透，再去碰复杂体系。这个顺序别反过来。