3. 分子动力学基本原理

分子动力学（MD）把材料看作由原子或离子组成的多体系统，用经典力学追踪每个粒子的运动轨迹，再通过统计力学把微观轨迹换算成宏观量。理解 MD 的关键，是同时抓住五件事：粒子按什么方程运动、粒子之间如何相互作用、系统处在什么统计环境、边界如何处理、时间如何离散推进。这五者共同决定一次模拟是否物理可信。

牛顿运动方程：轨迹的物理内核

对第 i 个原子，位置为 r_i、质量为 m_i，其运动遵循牛顿第二定律：

m_i · d²r_i/dt² = F_i = -∇_i U(r_1, r_2, …, r_N)

其中 U 是系统总势能，F_i 是其他所有原子施加给它的合力。求解这组耦合二阶微分方程，就得到每个原子的位置与速度随时间的演化。速度与位置的关系为 v_i = dr_i/dt，因此只要给定初始位置和初始速度，轨迹在经典近似下完全确定。MD 不追踪电子，原子间作用力全部由势函数提供，这是它区别于第一性原理方法的核心假设。

势函数：相互作用的来源

势函数决定模拟能描述什么材料、什么性质。工程中最常用的是嵌入原子法（EAM）势，其能量写成：

U = Σ_i F_i(ρ_i) + ½ Σ_{i≠j} φ_ij(r_ij)
ρ_i = Σ_{j≠i} f_j(r_ij)

式中 φ_ij 是原子对势，f_j 是电子密度贡献，F_i 是嵌入能，ρ_i 是原子 i 处的背景电子密度。EAM 能较好描述金属的弹性、缺陷与相变。共价材料常用 Stillinger-Weber 或 Tersoff 势，离子晶体用 Buckingham 势加库仑项，有机与生物体系用 CHARMM、AMBER 等力场。选择势函数时必须核对它的适用范围、截断半径和是否经过目标性质验证，不能只看名称。

系综：控制温度与压力

系综规定系统与外界交换什么。常见类型包括：

• NVE（微正则）：粒子数、体积、能量恒定，用于检验能量守恒与积分稳定性。
• NVT（正则）：温度恒定，通过 Nosé-Hoover 或 Berendsen 热浴实现，适合研究恒温下的结构演化。
• NPT（等温等压）：温度和压力恒定，通过 Parrinello-Rahman 或 Berendsen 控压实现，用于获得平衡晶格常数、热膨胀与相变。
• NPH：压力与焓恒定，较少单独使用。

温度由动能统计得到：T = 2K/(3N k_B)，压力由维里定理计算。系综选择错误会直接导致结果失真，例如用 NVE 研究热膨胀就得不到正确趋势。

边界条件：消除表面假象

块体材料模拟普遍采用周期性边界条件（PBC）：模拟盒在三个方向重复，原子从一侧穿出时从对侧穿入。这样少量原子即可代表无限大晶体，避免自由表面带来的额外能量与虚假表面效应。PBC 下计算原子间作用力采用最小镜像约定，即只考虑最近镜像，因此截断半径必须小于模拟盒最短边长的一半。对于纳米线、薄膜、裂纹或表面吸附等场景，则需在相应方向使用自由边界，或采用固定层加自由层的混合边界。边界条件与截断半径、盒子尺寸必须匹配，否则会出现原子与自身镜像相互作用的非物理结果。

时间积分算法：把连续方程离散化

MD 最常用的积分器是速度 Verlet 算法，其更新顺序为：

v(t+Δt/2) = v(t) + F(t)·Δt/(2m)
r(t+Δt) = r(t) + v(t+Δt/2)·Δt
计算 F(t+Δt)
v(t+Δt) = v(t+Δt/2) + F(t+Δt)·Δt/(2m)

该算法时间可逆、能量长期漂移小，是 NVE 模拟的首选。时间步长 Δt 通常取 1 fs 左右；若体系含氢原子或高温振动，应减小到 0.5 fs 或更小，否则高频振动会失稳。初始速度一般按目标温度从麦克斯韦-玻尔兹曼分布抽样，并扣除整体平动以保证总动量为零。模拟流程通常为：能量最小化消除原子重叠，然后在 NVT 或 NPT 下弛豫到平衡，再切换到目标系综采集数据。判断平衡的常用依据是温度、势能、压力随时间围绕稳定均值波动，且均方根偏差不再单调漂移。

一次可复现的示例设置

以面心立方铜块体为例，可按下表建立一次教学级模拟。以下参数为合理示例，用于说明配置逻辑，不代表唯一正确取值。

运行后应检查：势能与温度是否围绕稳定均值波动；压力是否在目标值附近；均方根偏差是否趋于平台。若能量持续单向漂移，优先减小时间步长或检查截断与邻居列表设置；若温度失控，检查热浴耦合强度与初始速度分布。只有确认平衡后，再统计径向分布函数、均方位移或弹性常数，结果才具有物理意义。