3. 分子动力学核心原理
分子动力学（MD）模拟的本质，是把原子或分子看成在势能面上运动的经典粒子，用牛顿第二定律逐步推进它们的轨迹。对扩散与迁移研究而言，我们真正关心的不是某一瞬间的原子位置，而是原子在长时间尺度上如何从一个格点或自由体积跳到另一个位置。因此，理解“方程—势函数—系综—积分—控温控压”这条链条，并据此在精度与效率之间做出取舍，是开展可靠模拟的前提。
牛顿运动方程与势函数：模拟的物理内核
对第 i 个原子，牛顿方程写作 m_i·a_i = F_i = -∇_i U(r_1, r_2, …, r_N)。其中 U 是体系总势能，F_i 是作用在该原子上的合力。势函数决定了原子间如何相互作用，也决定了扩散势垒、迁移路径和缺陷稳定性。工程上最常用的是经验势（力场），例如金属体系中的嵌入原子法（EAM）、共价体系中的 Stillinger-Weber 或 Tersoff、有机与生物体系中的 AMBER/OPLS/CHARMM，以及离子晶体中的 Buckingham 势。
以 EAM 为例，总能量可写为 U = Σ_i F_i(ρ_i) + ½Σ_{i≠j} φ_ij(r_ij)，其中 φ_ij 是对势，F_i 是嵌入能，ρ_i 是电子密度。对势描述两体排斥与吸引，嵌入能则体现多体效应，这正是金属中空位迁移、表面扩散能被较好再现的原因。势函数的截断半径 r_cut 直接影响计算量：r_cut 越大，邻居数越多，单步耗时近似按 r_cut³ 增长。因此，在保证扩散势垒和迁移机制不失真的前提下，应优先选择经过验证、截断合理的势函数，而不是盲目追求复杂形式。
系综与边界条件：决定“统计什么”和“边界如何”
系综规定了模拟中哪些热力学量被固定。扩散研究常用 NVT（粒子数、体积、温度恒定）和 NPT（粒子数、压力、温度恒定）。NVT 适合研究给定体积下的自扩散系数和均方位移；NPT 则更适合模拟真实材料在常压下的热膨胀、相变和压力驱动的迁移过程。若只关心本征扩散，NVE 也可用于检验能量守恒和积分稳定性。
边界条件同样关键。三维块体扩散通常采用周期性边界条件（PBC），使原子从一侧离开后从对侧进入，从而消除表面效应。对于表面扩散、晶界扩散或纳米线中的迁移，则需在相应方向使用自由边界或固定边界，以保留真实表面和界面结构。PBC 下必须注意最小镜像约定：原子间距离取最近镜像，且截断半径应小于模拟盒长的一半，否则会出现自相互作用伪影。
时间步长与积分算法：稳定推进的关键
MD 最常用的积分器是 Verlet 算法及其速度形式。速度 Verlet 的更新为：
r(t+Δt) = r(t) + v(t)·Δt + ½·a(t)·Δt²
v(t+Δt) = v(t) + ½·[a(t) + a(t+Δt)]·Δt它时间可逆、能量长期漂移小，且只需存储一个时刻的位置、速度和加速度。时间步长 Δt 的选取必须小于体系最快振动周期的约 1/10。对含轻原子（如 H）的体系，O-H 或 C-H 伸缩振动周期约 10 fs，因此 Δt 常取 0.5–1 fs；对金属或无机氧化物，可取 1–2 fs。若 Δt 过大，会出现原子重叠、能量急剧漂移，甚至模拟崩溃；若过小，则计算量成倍增加，难以覆盖扩散所需的纳秒级时间。
为兼顾精度与效率，工程上常采用约束算法（如 SHAKE、LINCS）固定键长，从而允许使用 2 fs 甚至更大的步长；也可使用多时间步算法，对慢变力用大步长、快变力用小步长。但需注意，约束会改变体系动力学，对扩散系数和迁移势垒可能产生系统偏差，必须通过对比无约束模拟进行验证。
温度与压力控制：让系综真正成立
控温方法中，Berendsen 热浴通过速度标定使温度快速趋近目标值，适合平衡阶段，但它不产生严格的正则系综，温度涨落偏小。Nosé-Hoover 热浴引入热浴自由度，能产生正确的 NVT 系综，适合采样阶段。实际流程常先用 Berendsen 平衡 100–200 ps，再切换 Nosé-Hoover 采样 1–5 ns。控压方面，Berendsen 控压适合快速调整盒长，Parrinello-Rahman 或 MTK 控压则能给出更合理的体积涨落，适合 NPT 下的扩散与相变研究。
需要强调的是，控温控压的耦合强度不能过大。过强的热浴会使体系像被“强行搅拌”，破坏真实动力学；过弱则温度、压力漂移明显。一般建议温度耦合时间常数取 0.1–1 ps，压力耦合时间常数取 1–5 ps，并根据体系大小和扩散快慢微调。
面向扩散迁移的精度与效率权衡
扩散模拟的核心矛盾是：扩散事件发生在纳秒到微秒尺度，而原子振动在飞秒尺度。若直接用 Δt = 1 fs 模拟 1 μs，需要 10⁹ 步，计算量巨大。因此，实际工作中常采用以下策略：第一，用短时间、多温度的 MD 提取均方位移（MSD）和扩散系数，再通过 Arrhenius 关系外推；第二，用 NEB、CI-NEB 或过渡态搜索计算迁移势垒，与 MD 结果交叉验证；第三，对长时间扩散，可结合加速分子动力学（hyperdynamics）、温度加速动力学或动力学蒙特卡洛（kMC），把 MD 得到的势垒和频率作为输入。
判断模拟是否可靠，不能只看轨迹“看起来合理”。应检查：能量守恒是否稳定、温度是否在目标值附近合理涨落、MSD 是否进入线性区、扩散系数随温度是否满足 Arrhenius 关系、不同初始构型和不同势函数给出的结果是否一致。只有这些判据同时成立，扩散系数和迁移势垒才具有工程参考价值。