说明:本文华算科技主要介绍 MD 模拟时长与目标物理量收敛的关系,以及初始弛豫、时间相关、有效样本数和稀有事件对统计结果的影响。


MD模拟时长与统计收敛的基本含义
分子动力学按照粒子受力积分运动方程,生成位置、速度等量随时间变化的轨迹。模拟的物理时长为 tsim = NstepΔt,其中 Nstep 是积分步数,Δt 是时间步长。计算机运行一天与材料演化一纳秒属于不同时间;壁钟耗时受体系大小、受力计算和硬件影响,轨迹中的物理时间由积分步长累计。
用于平衡性质统计的轨迹通常分为初始弛豫段和生产段。初始结构的密度、局部配位或链构象可能偏离所研究系综的典型分布,这段演化的数据会使均值产生偏置。剔除初始时段后,生产时间 tprod 才用于估计密度、配位数等平衡平均;初始弛豫所占比例随模型和目标量改变。

对于目标量 A,生产段的估计值可写为 Ā = N−1Σi=1NAi。统计收敛针对某个量及其允许误差:均值在轨迹延长或独立重复后仍处于可接受的误差范围,才适合用于比较材料。原子坐标继续涨落并不妨碍平均密度收敛;平均密度已稳定时,扩散系数仍可能具有很大的统计误差。
经典力场 MD、AIMD 和机器学习势 MD 使用不同受力模型,势垒高低与势阱形状会改变模型中的弛豫时间。AIMD 每一步求电子结构的代价高,少量皮秒轨迹仍只覆盖这段物理时间。增大 Δt 会改变积分误差,过大的步长可能扭曲振动或造成能量漂移;延长采样与控制积分误差分别约束统计结果。


温度能量平稳后目标物理量为何仍变
恒温器调节动能,使温度在设定值附近涨落;链构象交换、晶界重排和离子跨位点迁移具有各自的弛豫时间。这些过程尚未充分发生时,温度曲线仍可相当平稳。晶体局部键长可以围绕一个极小值持续振动,而相变所需的集体位移在整段轨迹中都未出现。
液氩从偏高密度的初始构型出发时,密度经历显著松弛。把这段初始松弛保留在累积平均中,会使均值长期偏高,即使后续瞬时密度已经在平衡区间涨落。原论文以多条重复轨迹的平均展示这种偏置;横轴的 τ 是其模型定义的约化时间单位,图中的数值没有通用于其他材料的皮秒含义。

改变舍弃的起始时间后,若目标均值仍发生超出统计误差的移动,生产段可能残留初始偏置。舍弃更多数据也会缩短可用生产段,固定删除前半条轨迹未必能得到足够样本。密度、配位和构象比例的平均值可以分别检查,因为它们对同一初始结构的松弛速度可能不同。
正丁烷的 C–C–C 键角持续快速涨落,C–C–C–C 主链二面角却可在某个构象附近停留较久,随后才发生交换。同一条轨迹内,振动与构象占据的时间尺度已经相差明显。聚合物链段、柔性框架和分子吸附层也可能存在局部振动与整体重排的分离,目标量若依赖慢构象占据,短时间内的快速涨落只统计了其中一部分构型。

局部平均键长、位点占据和自扩散系数分别描述振动分布、状态比例与长程位移。结构快照保持原有形貌,只记录观察窗口内的结构演化;计算平衡占据需要构型交换,计算扩散则需要足够的长程位移。三类结果对应的生产时间可以不同。


存帧很多,有效样本数为何仍少?
相邻构型从同一动力学过程连续产生,Ai 与 Ai+k 通常存在相关。归一化自相关函数 ρA(k) 描述间隔 k 次采样后涨落的相关程度。正丁烷构象占据的自相关衰减记录状态持续时间,振荡或零附近的尾部还会受到有限轨迹噪声影响。图中的 trans、g+ 和 g− 是构象标签。
对已经平稳、相关函数可积且生产段足够长的序列,统计低效因子可近似写为 gA = 1 + 2Σk≥1ρA(k),相应有效样本数 Neff ≈ N/gA。它表示这 N 帧对 A 的均值估计,相当于多少个独立样本。实际估计须控制相关尾部的噪声,直接把所有滞后点累加可能得到失真的相关时间。

均值的标准误差近似为 SE(Ā) ≈ σA/√Neff,其中 σA 是 A 本身的涨落宽度。轨迹标准差与均值标准误差描述不同对象:前者可在热平衡下保持非零,后者随独立信息积累而减小。把存帧间隔缩短,会增加高度相似的构型,N 与 gA 往往一起增大,有效样本数未必随文件帧数同比增加。
分块平均把生产轨迹分成 B 段,分别计算块均值,再估计 SE ≈ sblock/√B。短块之间仍相关时,误差估计通常偏小;分块长度超过主要相关时间后,均值误差随块长出现平台。块长继续增加又会使剩余块数变少,误差估计开始波动,因此平台与足够的块数须在同一段数据中出现。

在上述平稳采样条件下,若均值标准误差为 SEnow,目标为 ε,可用 treq ≈ tprod(SEnow/ε)2 估计所需生产时间。把均值标准误差减半,通常需要约四倍生产时间,前提是目标量的分布与相关时间未变。这个外推适用于已充分取样后的均值误差;轨迹从未访问的慢状态不会自动出现在该估计中。


独立重复和稀有事件改变所需时长
同一构型区域内可以获得许多去相关的振动样本,但不同区域的占据比例仍可能未采样。有效样本数随所选物理量而变,快速键长变量的相关时间可能很短,跨配位状态的指示变量却长期不变。构型空间中的状态划分有助于识别轨迹究竟访问过哪些区域,以及各区域之间发生了多少次转换。
从不同随机速度开始的轨迹可能给出不同的位点占据比例,共用同一初始结构的重复模拟仍可能保留相同的构型偏置。若材料存在不同链构象、缺陷排列或界面水配位状态,从这些候选状态分别开始演化,状态占据和目标量随生产时间的变化才有机会显示尚未消除的初始依赖。

MSD 中不同滞后时间共享同一批原子和时间起点,相邻点通常相关;线性拟合的点数不能替代独立扩散事件数。较长滞后时间剩余的可用起点更少,少数大位移即可显著改变斜率。扩散系数的误差估计应保留这些相关,并检查所选区间是否已进入普通扩散过程。
对罕见的缺陷跃迁或重构,若等待时间近似独立指数分布且事件速率 k 恒定,观察时间 tobs 内零次事件的概率为 P(0) = exp(−ktobs)。零次事件只约束该模型下的事件速率:单侧 95% 上限约为 k obs。这个结果依赖泊松过程假设;存在诱导期、相关跃迁或逐渐变化的结构时,相应速率模型须改变。

高温轨迹可以增加跨位点跃迁,但降温后的等待时间可能显著增长。用高温数据外推低温扩散系数,需要相同的物相和迁移机制;相变、缺陷重新分布或传输机制改变会使单一 Arrhenius 拟合失效。增强采样也会改变构型访问频率,偏置轨迹的经过时间通常须经过专门的动力学处理才能用于真实速率。
长轨迹降低采样误差时,势函数误差和有限盒尺寸效应仍可保留。对于机器学习势或经典力场,目标温度下的配位与跃迁构型必须处于模型经过检验的化学环境中。若某类离子迁移势垒持续偏低,增加轨迹长度会积累更多该类跃迁,其过快的跳跃速率仍会进入扩散系数。
