4. 扩散迁移模拟的势函数选择与构建

在扩散与迁移的分子动力学模拟中，势函数决定了原子间如何相互作用，也直接决定了扩散能垒、迁移路径、缺陷形成能以及温度依赖关系。换句话说，势函数不是“背景参数”，而是整个模拟的物理模型本身。选择或构建一个不合适的势函数，即使积分步长、系综和模拟时间都设置正确，得到的扩散系数和迁移机制也可能完全偏离真实材料行为。因此，本节围绕五类常用势函数——Lennard-Jones、EAM、MEAM、ReaxFF 和机器学习势——讨论它们的适用边界、对扩散能垒与迁移路径的影响，并给出工程材料势函数验证与参数化的可执行流程。

一、五类势函数的能力边界

Lennard-Jones（LJ）势只描述范德华式的短程排斥与长程吸引，形式简单、计算极快，适合惰性气体、简单液体和表面吸附等体系。但它没有方向性键合，也无法区分不同金属的电子密度效应，因此对金属中的空位迁移、晶界扩散等问题的描述能力有限。EAM（嵌入原子法）势通过嵌入能加电子密度项，能够较好描述金属键和密堆积金属的弹性、缺陷与扩散行为，是金属扩散模拟中最常用的经验势之一。MEAM（修正嵌入原子法）在 EAM 基础上引入角度依赖的电子密度，能描述共价性较强的金属间化合物、半导体和部分合金中的方向性键合，对迁移路径的各向异性更敏感。

ReaxFF 属于反应力场，允许键断裂与形成，适合描述氧化物、碳氢体系、界面反应和化学扩散过程，能够给出反应路径与能垒，但计算成本显著高于 EAM/MEAM。机器学习势（MLP）以第一性原理数据为训练集，通过高维势能面拟合获得接近 DFT 的精度，同时保持远低于 DFT 的计算成本，适合复杂合金、界面和缺陷扩散研究，但其可靠性完全取决于训练集覆盖度与外推能力。

二、势函数如何影响扩散能垒与迁移路径

扩散能垒 Eₐ 是原子从一个稳定位置迁移到相邻稳定位置所需跨越的能量势垒。势函数的排斥项斜率、吸引项深度和截断半径都会改变势能面形状，从而改变 Eₐ。以面心立方金属中的空位迁移为例，EAM 势通常给出 0.6～0.8 eV 的迁移能垒，而 LJ 势由于缺乏电子密度嵌入项，往往低估能垒并给出过于平滑的势能面，导致迁移路径趋于直线化。MEAM 引入角度项后，在 B2 或 L1₂ 结构合金中会打开不同的迁移通道，使最近邻与次近邻路径的能垒差增大，迁移路径呈现明显的各向异性。

ReaxFF 对氧化物中氧空位迁移的描述更接近 DFT，因为它允许金属—氧键在迁移过程中部分断裂与重组，能垒曲线会出现多个鞍点，对应协同迁移机制。机器学习势则可以在训练集覆盖的构型空间内复现 DFT 级别的能垒和路径，但一旦迁移路径涉及训练集未包含的局部环境，预测能垒可能出现系统性偏差。因此，势函数选择的核心判断是：目标扩散过程是否落在该势函数的“可信描述域”内。

三、工程材料势函数验证与参数化流程

工程中不能直接拿一个文献势函数就投入生产模拟，必须经过验证与必要的参数化。下面给出一套可执行的流程，并以“Cu 中空位扩散”为教学示例，所有数值均为示例性合理假设，不代表实际计算结果。

第一步：明确目标量与参考数据。确定需要复现的物理量，例如空位形成能 E_f、迁移能垒 Eₐ、扩散前因子 D₀ 和 300～900 K 的扩散系数。参考数据优先来自实验或高精度 DFT/NEB 计算。示例目标：E_f ≈ 1.27 eV，Eₐ ≈ 0.70 eV。

第二步：候选势函数筛选。对金属体系优先选 EAM/MEAM；若涉及氧化物或界面反应则选 ReaxFF；若要求 DFT 精度且训练数据充足则选 MLP。示例中选取一个已发表的 Cu EAM 势和一个 MEAM 势作为候选。

第三步：基础物性验证。在 0 K 下计算晶格常数 a₀、内聚能 E_coh、弹性常数 C₁₁、C₁₂、C₄₄ 和空位形成能。示例判据：a₀ 与实验偏差小于 1%，E_f 偏差小于 0.1 eV。若偏差过大，说明势函数不适合该体系，应更换或重新拟合。

第四步：扩散能垒与路径计算。用 NEB 或 CI-NEB 计算空位迁移能垒，同时提取迁移路径上的原子位移。示例命令片段如下：

# LAMMPS 中 NEB 计算空位迁移能垒（示例）
min_style       fire
neb             0.0 0.1 100 500 5.0
neighbor        2.0 bin
fix             neb all neb 1.0
# 初始态与末态构型分别读取 initial.data 和 final.data
# 输出 neb.log 中提取能垒曲线

计算完成后，检查能垒曲线是否光滑、鞍点是否唯一。若出现多个鞍点，说明迁移路径可能涉及协同机制，需要结合 ReaxFF 或 MLP 进一步确认。

第五步：有限温度扩散系数验证。在 300～900 K 范围内运行 NVT 或 NPT 系综，计算均方位移 MSD，通过 Einstein 关系 D = lim(t→∞) MSD/(6t) 提取扩散系数，再用 Arrhenius 关系 D = D₀ exp(−Eₐ/k_BT) 拟合。示例判据：拟合得到的 Eₐ 与 NEB 结果偏差小于 0.05 eV，且 lnD 与 1/T 线性度 R² > 0.98。

第六步：与实验或 DFT 对比迭代。若 Eₐ 或 D₀ 偏差超过阈值，则回到势函数参数化阶段，调整截断半径、电子密度函数或嵌入能曲线，重新执行第三至第五步。对于 MLP，则需要补充训练集并重新训练，直到验证指标满足要求。

四、选择策略与常见误区

实际工程中，选择势函数应遵循“先匹配物理机制，再匹配精度与成本”的原则。金属体扩散优先 EAM/MEAM；氧化物和反应扩散优先 ReaxFF；复杂合金或界面扩散且已有高质量 DFT 数据时优先 MLP。常见误区包括：用 LJ 势研究金属空位扩散；用 EAM 势研究氧化物中的氧迁移；用 MLP 外推到训练集未覆盖的缺陷构型；以及只验证晶格常数而不验证扩散能垒。后者尤其危险，因为一个势函数可以给出正确的晶格常数，却给出完全错误的迁移能垒。

此外，势函数的截断半径和光滑处理也会影响扩散结果。截断半径过小会导致势能面不连续，使 NEB 能垒曲线出现伪振荡；截断半径过大则增加计算量。建议在验证阶段同时检查能量守恒和力收敛，确保势函数在目标温度下数值稳定。只有把势函数选择、验证和参数化作为一个整体工程流程来执行，扩散与迁移模拟的结果才具备可解释性和可复现性。