3. 分子动力学核心原理

分子动力学(MD)模拟的本质,是把材料中原子看作遵循经典力学的质点,通过数值积分求解牛顿运动方程,得到每个原子在相空间中的运动轨迹,再对轨迹做统计平均,得到宏观可观测的物理量。辐照损伤模拟中,级联碰撞、离位峰演化、缺陷迁移与复合,都是这一套力学与统计框架下的自然结果。本节围绕六个相互耦合的核心要素展开:牛顿运动方程、势函数、周期性边界条件、系综与温度控制、时间步长与积分算法、模拟尺度与限制。

牛顿运动方程:模拟的动力学引擎

对第 i 个原子,位置为 ri、质量为 mi,其运动由

mi d2ri/dt2 = Fi = −∇iU(r1, …, rN)

描述。这里 U 是体系总势能,Fi 是原子所受合力。方程右端只依赖原子位置,因此只要给定初始位置与初始速度,体系后续演化完全确定。辐照级联中,PKA(初级撞出原子)获得 keV 量级动能后,正是通过这一方程把能量逐级传递给周围原子,形成离位级联。

需要强调参照约定:本文所有位置、速度、力均在同一实验室坐标系下定义,位移方向以初始晶格位置为原点,向外为正。切换参照系时必须显式说明,否则缺陷迁移方向、级联取向等结论会失去可比性。

势函数:决定“原子如何相互作用”

势函数 U 是 MD 模拟中唯一承载材料物理信息的输入。它把原子构型映射为能量,进而给出力。辐照损伤常用势函数分三类:

  • 对势(如 Lennard-Jones):只依赖原子对距离,计算快,但无法正确描述金属键与缺陷形成能。
  • 多体势(如 EAM、Finnis-Sinclair):嵌入能加对势,能较好描述金属的弹性、空位与间隙形成能,是金属辐照模拟的主力。
  • 键级势(如 Tersoff、Stillinger-Weber):适合共价材料(Si、SiC),能描述键角与配位变化。

选择势函数时,必须核对其是否经过辐照相关量验证:离位阈能 Ed、缺陷形成能、迁移能垒、级联存活缺陷数。若势函数只拟合了平衡态性质,级联模拟结果可能严重失真。例如,某 EAM 势若未正确再现 Fe 的 Ed≈40 eV,则级联中产生的 Frenkel 对数量会系统性偏离实验与 DFT 参考值。

周期性边界条件:用小盒子模拟大块体

真实材料宏观上近似无限大,而 MD 只能模拟有限原子数。周期性边界条件(PBC)通过让模拟盒在三个方向重复,使边界原子与对侧镜像原子相互作用,从而消除表面效应。实现要点:

  • 原子坐标超出盒子时,按盒长取模折回;
  • 计算原子间距离时采用最小镜像约定,即取 rij 与 rij−L 中较短者;
  • 截断半径 rc 必须小于 L/2,否则同一原子会与自身镜像作用。

对辐照级联,PBC 还有一个隐含约束:级联不能跨越盒子边界“自我碰撞”。若 PKA 能量过高、级联尺寸接近盒长,热峰会与自身镜像重叠,导致缺陷数虚高。工程上通常要求级联最大尺寸不超过盒长的 1/3~1/2。

系综与温度控制:让模拟对应真实实验条件

MD 可在不同系综下运行。辐照损伤模拟最常用的是 NVE(微正则,能量守恒)与 NVT(正则,温度恒定)。级联初期通常用 NVE 让能量自然分配,避免恒温器把级联热峰“抽走”;级联冷却后或做长时间退火时,切换到 NVT 控制温度。

常见控温方法:

  • 速度标定:每步把速度缩放到目标温度,简单但破坏动力学,不推荐用于级联。
  • Berendsen 恒温器:以弛豫时间 τ 缓慢拉向目标温度,适合平衡态,但系综不严格。
  • Nosé-Hoover 恒温器:引入热浴自由度,能产生正确的正则系综,适合长时间退火与缺陷迁移研究。
  • Langevin 恒温器:加阻尼与随机力,适合含电子阻止或需要耗散的场景。

温度与动能的关系为 T = 2K/(3NkB),其中 K 为总动能。级联模拟中,初始温度通常设为 300 K,PKA 动能单独注入,不参与初始温度统计。

时间步长与积分算法:精度与稳定性的平衡

MD 用有限差分积分牛顿方程。最常用的是 Velocity-Verlet 算法:

v(t+Δt/2) = v(t) + F(t)Δt/(2m)
r(t+Δt)   = r(t) + v(t+Δt/2)Δt
F(t+Δt)   = -∇U(r(t+Δt))
v(t+Δt)   = v(t+Δt/2) + F(t+Δt)Δt/(2m)

该算法时间可逆、能量长期漂移小,是辐照 MD 的标准选择。时间步长 Δt 需满足:Δt < 0.1 Tmin,其中 Tmin 是体系最高振动频率对应的周期。对 Fe、W 等金属,典型 Δt = 1~2 fs;对含轻元素(H、He)或高能级联初期,需降到 0.1~0.5 fs。级联热峰阶段原子速度极高,若 Δt 过大,能量会不守恒,缺陷数出现虚假增长。

模拟尺度与限制:知道边界在哪里

MD 的时空尺度受计算资源硬约束。典型辐照级联模拟:原子数 105~107,时间 10~100 ps,对应级联从产生到冷却。这带来几个根本限制:

  • 时间尺度:无法直接模拟长时间扩散(秒级),需借助加速方法(如 ART、TAD)。
  • 尺寸效应:小盒子会抑制级联充分发展,缺陷数偏低或偏高,需做尺寸收敛测试。
  • 电子阻止:经典 MD 不含电子激发,高能 PKA 需额外加入电子阻止力。
  • 势函数精度:所有结论受限于势函数对辐照缺陷的再现能力。

下面用一张关系图把六个核心要素的耦合关系与关键参数约束集中呈现,便于在搭建模拟时快速核对。

分子动力学核心要素耦合关系与参数约束 牛顿运动方程 m·d²r/dt² = −∇U 轨迹演化引擎 势函数 U(r) EAM / Tersoff / LJ 决定力与缺陷能 需验证 Ed、迁移能垒 周期性边界条件 最小镜像约定 rc < L/2 级联尺寸 ≤ L/3 系综与温度控制 NVE / NVT Nosé-Hoover / Langevin T = 2K/(3NkB) 时间步长与积分 Velocity-Verlet Δt = 1~2 fs(金属) 级联初期 0.1~0.5 fs 模拟尺度与限制 10⁵~10⁷ 原子 · 10~100 ps 尺寸效应 · 电子阻止 · 势精度 提供 F = −∇U 边界镜像作用 轨迹统计 数值离散 控温影响热峰 Δt 限制总时长 关键约束 rc < L/2 · Δt < 0.1 Tmin 级联尺寸 ≤ L/3

把这张图当作搭建模拟时的检查清单:先确认势函数对辐照缺陷的再现能力,再根据级联最大尺寸确定盒长与截断半径,随后选择系综与恒温器,最后按元素种类和级联阶段设定时间步长。任何一项参数越界,都会让后续的缺陷统计失去物理意义。

在级联模拟中,建议先用 NVE 运行 5~10 ps 让热峰自然冷却,再切换到 Nosé-Hoover NVT 做退火。这样既保留级联初期的真实动力学,又能在长时间尺度上获得稳定的温度统计。
不要用速度标定法在级联热峰阶段控温。热峰中原子速度分布本身偏离平衡态,强行标定会人为抽走能量,导致缺陷数被严重低估。