分子动力学模拟全解析：微观状态、时间推进、轨迹统计、时间相关函数及计算条件限制

说明：本文华算科技主要介绍分子动力学模拟中的微观状态、离散时间推进、轨迹统计、时间相关函数，以及模型、时长和体系尺寸对计算结论的限制。

基本概念：微观状态怎样随时间推进？

分子动力学模拟，简称 MD，描述一组粒子在给定相互作用和外界条件下随时间演化的过程。这里的“粒子”可以是原子、离子，也可以是粗粒化模型中的珠子。某一时刻的计算状态由全部粒子的位置 r_i 和动量 p_i 组成，常写成 Γ(t) = {r_i(t), p_i(t)}。一帧坐标只给出构型；加入速度或动量，才足以继续推进经典动力学。

经典 MD 先用势能函数 U(R) 计算每个粒子的受力 F_i = −∂U/∂r_i，再数值求解 m_id²r_i/dt² = F_i。这一步产生下一个时刻的位置和速度，随后重新计算受力。循环执行后得到按时间排序的微观状态序列，也就是轨迹。势能可以来自经验势、机器学习势或即时电子结构计算；它决定相同构型下采用什么力。

图1的红线表示连续运动，黑色虚线表示数值积分给出的离散近似。积分时间步 Δt、轨迹保存间隔和总模拟时长是三个不同量：时间步控制方程推进的分辨率，保存间隔决定写入多少帧，总时长决定轨迹可能访问多慢的过程。每一步都保存坐标并不会增加已经计算的物理时间，稀疏保存也不会把较大的时间步变小。

时间步过大时，最快振动不能得到足够分辨，积分误差会积累，严重时轨迹发散。恒能量条件下的总能量漂移可用于检查数值推进。图2采用钨 BCC 超胞，在 3000 K 下比较多种机器学习势与 SNAP 势；小时间步区间内，能量偏差近似随 Δt² 增长，这与所用 Verlet 类积分的二阶误差相符。稳定推进只说明离散积分在当前模型中表现正常，它不检验势能函数能否描述熔化、缺陷迁移或成键变化。

单条轨迹为何不能等同于真实录像？

初始坐标、初始速度和温压控制方式共同选定一次模拟的起点。即使受力模型完全相同，稍微改变初速度，经过一段时间后也会得到另一条微观路径；含随机热浴时，随机数序列还会改变每次推进。原子在某一帧中的具体去向通常不是可重复的宏观结论。可重复比较的是指定模型与统计系综下的分布、平均量、涨落和相关时间。

图3给出甲烷在无序多孔碳中的运动。彩色折线是一条分子的代表性位移轨迹，它包含停留和跳跃；背景色表示全部分子与多个时间起点统计得到的自相关分布，黑色虚线才是总体均方位移对应的尺度。单条路径适合识别“困在孔腔后跳到相邻区域”的机制，扩散系数则来自足够多粒子、时间起点与长时间区间的平均。挑一条看起来平滑或跳跃频繁的路径去代表全部粒子，会把轨迹选择引入结论。

平衡态统计要先排除初始构型带来的暂态偏差。生产段中的温度、密度、势能或目标结构变量应围绕稳定分布波动；若它们持续单向漂移，直接计算平均值会混入松弛过程。去掉开头若干皮秒只是数据处理动作，判断依据应是观测量是否进入稳定分布，不能由一个固定时长替所有体系作决定。

相邻帧由前一帧推进而来，数据天然相关。保存一万帧时，有效样本数仍由该观测量的相关时间决定；若某个结构变量的自相关时间很长，有效样本数可能远小于帧数。多条采用不同初始速度或不同初始构型的独立轨迹，可以暴露状态占比是否依赖起点。轨迹越长只会扩大已访问区域内的统计量；一个从未跨越的高能垒状态不会凭平均运算自动出现。

相关轨迹怎样得到热力学和输运量？

轨迹中的每个状态都能映射成一个观测量 A[Γ(t)]。对平衡生产段求时间平均，在遍历充分且统计平稳时，可估计对应系综平均。密度、内能、径向分布函数都属于这种“由许多状态计数或平均”的量。若选择反应坐标 q 并统计概率 P(q)，可写成 F(q) = −k_BT ln P(q) + C。常数 C 只移动自由能零点，状态间的自由能差保持不变。

图4以石墨炔孔道中的质子转移为例：二维概率分布记录质子转移坐标 δ 与 O—O 距离 d_OO 的联合出现频率，自由能曲线由概率分布换算而来。两个高概率区对应质子更常停留的构型，δ = 0 附近的低概率区对应跨越区。自由能谷深浅取决于状态占比，势垒取决于跨越区与稳定区的概率比。采样没有到达的区域只能得到空白或很大的统计误差，无法据此断言该区域的真实自由能。

输运量还保留时间先后关系。Green–Kubo 公式把平衡涨落的时间相关函数积分成黏度、热导率等宏观系数。以黏度为例，要比较应力张量分量 P_ij(t₀) 与延迟时间后的 P_ij(t₀ + t)，并对许多时间起点 t₀ 平均。同一条长轨迹可以提供多个时间起点，但相邻起点彼此重叠，误差评估仍需考虑相关长度。

图5b中的单次相关积分在长延迟处分散很宽，红线是时间起点平均；图5c的短时区间分散较小。这里展示了一个常见判断：积分平台必须早于噪声主导区出现，并对截断时间保持可接受的变化。均方位移的线性斜率、速度自相关积分和热流自相关积分也要检查相应的时间区间，不能由曲线末端一个数替代区间检验。

模型、时间和尺寸怎样限定计算结论？

MD 输出始终以相互作用模型为前提。固定成键拓扑的经典力场适合研究振动、构象和非反应扩散，却无法自行生成模型中没有定义的断键与成键；机器学习势在训练构型附近可以接近参考电子结构精度，轨迹进入训练集未覆盖的局部环境时，力误差可能迅速增大。数值稳定、温度平稳和能量守恒都不能替代势函数验证。验证对象应与研究问题一致，例如晶格常数、弹性响应、缺陷能、液体结构或迁移能垒。

轨迹能否覆盖目标构型，可借助集体变量分布和自相关函数检查。图6比较丙氨酸二肽的普通 MD 与不确定度偏置 MD：构型空间覆盖率随时间上升，位置自相关衰减速度反映相邻构型失去记忆所需的时间，Ramachandran 图记录不同扭转角区域的访问情况。较高覆盖率只支持指定坐标访问更广；若存在未纳入分析的慢变量，轨迹仍可能停留在有限状态集合中。偏置力还会改变原始动力学与状态权重，未经重加权的偏置轨迹不能直接给出平衡概率或真实速率。

有限周期盒会重复同一结构，盒长小于相关结构或协同运动的尺度时，涨落和输运过程会受到约束。图7中的铝 Σ3 晶界在较小横截面下出现阶梯状 MSD，不同速度种子的扩散率分散较大；横截面扩大后，多个局域扩散事件可在不同位置发生，MSD 更接近连续增长。体系尺寸改变的不只是原子数，还包括可容纳的波长、缺陷间距与协同区域。

读取一项 MD 结论时，应明确粒子表示、相互作用模型、周期盒设置、统计系综、积分时间步、保存间隔、平衡段、生产段、独立轨迹数，以及物性公式采用的拟合或积分区间。若研究的是扩散，要看 MSD 的线性区、方向分量和尺寸检验；若研究的是相变，要看升降温路径、成核尺度和多次重复；若研究的是界面吸附，要确认表面覆盖度、溶剂层和电荷状态。这些计算条件共同规定数值对应哪一类微观过程，也规定它能与哪种实验量比较。