分子动力学模拟材料行为 - 第3节
第三章 原子间相互作用势：常见势函数的物理意义、适用材料类型及参数化方法
分子动力学（MD）模拟的精度与效率，从根本上取决于原子间相互作用势（Interatomic Potential）的选择。势函数是描述体系总能量与原子坐标之间关系的数学表达式，其梯度给出原子所受的力。本章将系统剖析五类最具代表性的势函数——Lennard-Jones（LJ）、EAM、MEAM、Tersoff 和 ReaxFF——从物理图像、数学形式、适用边界到参数化策略，帮助学习者在面对具体材料体系时，能够做出有理有据的“势函数选型决策”。
上图展示了势函数选型的核心决策逻辑。实际模拟中，没有任何一个势函数是“万能钥匙”。LJ 势适合描述闭壳层分子间的弱相互作用；EAM 及其扩展形式 MEAM 则针对金属键的“嵌入”本质；Tersoff 势抓住了共价键的方向性与饱和性；ReaxFF 则通过键级（Bond Order）概念实现了化学反应过程的动态描述。理解这些势函数的物理内涵，是正确使用 MD 模拟的前提。
3.1 Lennard-Jones 势：范德华相互作用的经典范式
3.1.1 核心原理：排斥与吸引的平衡
Lennard-Jones 势（又称 12-6 势）是形式最简单的双体势函数，其表达式为：
V(r) = 4ε [ (σ/r)^12 - (σ/r)^6 ]其中，ε（能量参数）代表势阱深度，反映两原子间结合强度的最大值；σ（距离参数）是势能为零时的原子间距，近似为原子碰撞直径。物理上，(σ/r)^12 项模拟电子云重叠时的强排斥力（Pauli 排斥），而 (σ/r)^6 项则源自瞬时偶极-诱导偶极相互作用（London 色散力）。该势函数仅依赖原子间距 r，不包含角度项，因此计算效率极高。
3.1.2 适用材料类型与局限性
LJ 势最成功的应用领域是惰性气体（Ar、Kr、Xe）、分子晶体（如 N₂、CO₂）以及层状材料层间相互作用的补充描述。对于这些体系，电子云分布近似球形对称，多体效应微弱，双体势假设成立。然而，LJ 势无法描述金属键、共价键或离子键的方向性。例如，将其用于模拟铁（Fe）的体心立方（BCC）结构，会错误地预测 FCC 结构更稳定，因为 LJ 势无法区分不同晶体结构的能量差异。
3.1.3 参数化方法：从实验数据到力场库
LJ 参数的获取通常有两种途径：一是拟合实验数据，如通过气体粘度、第二维里系数确定 ε 和 σ；二是采用混合规则（Lorentz-Berthelot 规则）从纯物质参数推算异种原子参数：
σ_ij = (σ_ii + σ_jj) / 2
ε_ij = sqrt(ε_ii * ε_jj)在分子力场（如 OPLS、CHARMM）中，LJ 参数通常与电荷参数协同优化，以复现液体的密度、汽化焓和径向分布函数。
3.2 EAM 势：金属体系的多体效应描述
3.2.1 核心原理：嵌入能 + 对势
嵌入原子方法（Embedded Atom Method, EAM）由 Daw 和 Baskes 于 1984 年提出，其核心思想是：将原子视为嵌入在由周围原子构成的电子密度背景中的“杂质”。总能量表达式为：
E_total = Σ_i F_i(ρ_i) + (1/2) Σ_i Σ_j≠i φ_ij(r_ij)其中，F_i 是嵌入能函数，依赖于原子 i 处的局域电子密度 ρ_i（由周围原子贡献叠加）；φ_ij 是传统的双体排斥项。嵌入能函数是非线性的，这自然引入了多体效应——当原子周围邻居较多时，每个键的强度会减弱，从而正确描述了金属的“饱和”特性。
3.2.2 适用材料类型：FCC 金属及其合金
EAM 势在面心立方（FCC）金属（如 Cu、Al、Ni、Au、Ag）及其合金体系中取得了巨大成功。它能够准确再现晶格常数、弹性常数、空位形成能、表面能以及层错能等关键物性。对于体心立方（BCC）金属（如 Fe、W、Mo），标准 EAM 形式往往不够精确，因为 BCC 金属的 d 电子成键具有明显的方向性。
3.2.3 参数化流程与实战案例
EAM 势的参数化是一个典型的“拟合-验证”闭环过程。以 Cu 的 EAM 势开发为例：
- 目标数据收集：从实验或第一性原理（DFT）计算获取晶格常数 a₀=3.615 Å、弹性常数 C₁₁/C₁₂/C₄₄、空位形成能 E_vf≈1.28 eV、层错能 γ_sf≈45 mJ/m²。
- 函数形式选择：采用 Johnson 分析型 EAM 或 Voter-Chen 型 EAM，确定嵌入函数 F(ρ) 和对势 φ(r) 的数学表达式。
- 数值优化：使用遗传算法或模拟退火算法，最小化目标函数
χ² = Σ w_i (P_i^calc - P_i^target)²，其中 w_i 为权重。 - 验证测试：计算未参与拟合的性质，如热膨胀系数、熔点、声子色散关系，确保势函数的可移植性。
3.3 MEAM 势：突破角度依赖的局限
3.3.1 核心原理：引入电子密度方向性
改进型嵌入原子方法（Modified EAM, MEAM）在 EAM 的基础上，将局域电子密度 ρ_i 从标量扩展为包含角度依赖的表达式。它通过引入 s、p、d、f 轨道的权重因子，使得电子密度分布具有方向性：
ρ_i = ρ_i^(0) * exp( -β * ( (ρ_i^(1)/ρ_i^(0))² + ... ) )这种形式使得 MEAM 能够描述共价键的取向性，同时保留金属键的嵌入特征。
3.3.2 适用材料类型：BCC 金属、共价金属及半导体
MEAM 势显著扩展了 EAM 的适用范围，尤其适用于BCC 金属（Fe、Cr、W）、共价金属（Si、Ge）以及金属间化合物（如 NiAl、TiAl）。例如，MEAM 势能够正确预测 Fe 的 BCC 结构稳定性，并再现其复杂的位错核结构。
3.3.3 参数化难点与注意事项
MEAM 势的参数数量远多于 EAM（通常有 15-20 个），包括电子云衰减系数 β、原子密度标度 ρ₀、嵌入函数参数 A 等。参数化过程高度依赖 DFT 参考数据，尤其是不稳定晶体结构的能量（如 BCC Fe 的 FCC 结构能量），这些数据无法从实验获得。初学者常犯的错误是直接使用文献中针对特定元素组合的参数，而忽略了这些参数对截断半径的敏感性。MEAM 的截断半径通常取 4.0-5.5 Å，过小的截断会导致虚假的力振荡。
3.4 Tersoff 势：共价键的方向性与键级
3.4.1 核心原理：键级依赖的键能
Tersoff 势（及其改进形式 Tersoff-Brenner、AIREBO）是描述共价体系的经典键级势。其核心思想是：原子 i 与 j 之间的键能不仅取决于 r_ij，还受到周围原子 k 的影响，通过键级参数 b_ij 实现：
V_ij = f_c(r_ij) [ A_ij exp(-λ_ij r_ij) - b_ij B_ij exp(-μ_ij r_ij) ]键级 b_ij 是键角 θ_ijk 的函数，当原子 i 与多个原子成键时，每个键的强度会因“竞争”而减弱，这自然描述了共价键的饱和性。例如，在金刚石中，每个 C 原子与 4 个邻居成键，每个键的键级约为 1；而在石墨中，每个 C 原子与 3 个邻居成键，键级约为 1.33。
3.4.2 适用材料类型：Si、Ge、C 及其化合物
Tersoff 势最经典的应用是半导体材料（Si、Ge、SiC）和碳材料（金刚石、石墨烯、碳纳米管）。它能够准确描述 sp³ 和 sp² 杂化态的转变，这是 LJ 势和 EAM 势完全无法实现的。例如，Tersoff 势可以模拟硅晶体在压痕下的非晶化相变，以及石墨烯在拉伸下的 Stone-Wales 缺陷形成。
3.4.3 参数化策略：从 DFT 到经验拟合
Tersoff 势的参数（A、B、λ、μ、n、c、d、h 等）通常通过拟合以下数据获得：
- 结合能与晶格常数：确保势函数正确预测晶体的内聚能。
- 弹性常数：特别是 C₁₁、C₁₂、C₄₄，保证力学行为准确。
- 小团簇能量：如 Si₂、Si₃ 的 DFT 能量，用于约束短程相互作用。
- 缺陷形成能：空位、间隙原子、表面重构能，提升缺陷模拟可靠性。
值得注意的是，Tersoff 势的原始参数（1988 年）对石墨烯的层间相互作用描述极差（无法给出层间距），后续的 AIREBO 势通过增加 Lennard-Jones 项解决了这一问题。因此，在模拟多层二维材料时，必须检查势函数是否包含层间范德华项。
3.5 ReaxFF 势：化学反应动力学的反应力场
3.5.1 核心原理：键级的动态演化
ReaxFF（Reactive Force Field）由 van Duin 等人于 2001 年提出，其革命性在于：键级（Bond Order, BO）不再是固定参数，而是根据原子间距离和局部环境连续变化的变量。总能量包含多项贡献：
E_system = E_bond + E_over + E_under + E_val + E_tors + E_vdWaals + E_Coulomb其中，E_bond 是键能（与 BO 相关），E_over/E_under 是过配位/欠配位校正项，E_val 是键角能，E_tors 是二面角能，E_vdWaals 和 E_Coulomb 分别描述非键相互作用。由于 BO 可以连续变化，ReaxFF 能够描述化学键的断裂与形成，而无需预先定义反应路径。
3.5.2 适用材料类型：含能材料、催化、燃烧
ReaxFF 特别适合模拟化学反应主导的过程：
- 含能材料分解：如 TATB、HMX 的热分解机理。
- 碳氢化合物燃烧：甲烷、正庚烷的氧化反应网络。
- 催化表面反应：金属（Pt、Ni）表面的 CO 氧化、水煤气变换反应。
- 材料腐蚀：金属在潮湿环境中的氧化膜生长。
3.5.3 参数化方法：大规模 DFT 数据拟合
ReaxFF 的参数化是极其繁重的任务，通常包含 50-100 个参数。标准流程为：
- 训练集构建：使用 DFT 计算大量构型的能量、力、电荷，包括反应物、产物、过渡态、中间体。一个典型的训练集包含 1000-5000 个构型。
- 电荷平衡：采用 QEq（电荷平衡）方法计算原子电荷，需拟合电负性 χ 和硬度 η 参数。
- 全局优化：使用并行遗传算法或共轭梯度法，最小化加权均方根误差。
- 验证：在分子动力学模拟中测试反应路径、活化能、反应热，确保动力学行为合理。
3.6 势函数选型对比与总结
| 势函数 | 物理基础 | 典型材料 | 计算成本 | 关键参数 | 主要局限 |
|---|---|---|---|---|---|
| LJ | 范德华双体作用 | Ar、N₂、分子晶体 | 极低 | ε, σ | 无方向性，不能描述化学键 |
| EAM | 嵌入能+对势 | FCC 金属（Cu、Al、Ni） | 低 | F(ρ), φ(r) | 无法处理强共价键 |
| MEAM | 角度依赖电子密度 | BCC 金属、Si、金属间化合物 | 中 | β, A, t 参数 | 参数多，拟合复杂 |
| Tersoff | 键级-键角关联 | Si、Ge、C、SiC | 中 | A, B, λ, μ, n, c, d, h | 层间范德华力缺失 |
| ReaxFF | 动态键级+多能项 | 含能材料、催化、燃烧 | 极高 | 50-100 个参数 | 参数移植性差，成本高 |
3.6.1 常见误区与复盘
- 误区一：盲目追求“高精度”势函数。对于纯金属的力学性能模拟，EAM 往往比 MEAM 更稳定，因为 MEAM 的复杂角度项可能导致能量曲面出现虚假局部极小值。
- 误区二：忽略截断半径的影响。所有势函数都需要设置截断半径。过小的截断会导致能量不连续，过大的截断则引入不必要的计算开销。建议对每个势函数进行截断半径收敛性测试。
- 误区三：直接使用文献参数而不验证。文献中的参数可能是在特定温度、压力或缺陷状态下拟合的。使用前应至少计算晶格常数、弹性常数和结合能，与实验值对比。
- 误区四：将 ReaxFF 用于非反应过程。如果模拟中不涉及化学键断裂，使用 ReaxFF 只会浪费计算资源，且其数值噪声可能大于 EAM。
本章系统梳理了五类核心势函数的物理图像与适用边界。在实际工程中，势函数的选择往往需要在精度、效率与可移植性之间做出权衡。建议学习者在完成本章学习后，针对自己研究的材料体系，建立一个小型“势函数验证基准集”，包含晶格常数、弹性常数、缺陷能、热膨胀系数等关键指标，以此作为评估任何新势函数的“试金石”。