扩散分子动力学模拟：从势函数选择到统计平均的完整方案

扩散过程的分子动力学模拟不能只算一个扩散系数就收工。真正的分析层次需要向下追溯到原子尺度的扩散机制——是空位跳跃、间隙迁移、还是协同换位——然后向上关联到宏观的可测量扩散率。

势函数的选取是扩散模拟中最影响结果质量的决定性因素。对于金属体系，EAM（嵌入原子势）对扩散相关性质的描述优于Lennard-Jones类对势——后者缺少多体效应，低估了鞍点能量，导致跳跃频率被高估。打个比方，用EAM算出的空位迁移势垒一般是0.5~1.0eV，用LJ对势可能跌到0.2~0.3eV——这不是数值误差，是物理模型层面就错了。对于离子晶体和氧化物，Born-Mayer-Huggins势或者Buckingham势外加长程库仑作用是目前的主流方案，但库仑项的长程求和需要用PPPM或Ewald方法，直接用截断半径会引入实质性的人工梯度。

温度区间的设计需要覆盖物理上有意义的范围。通常从略低于熔点往下降——空位扩散在熔点附近是主导机制，在低温下逐渐被晶界扩散和表面扩散取代。每降50~100K跑一组模拟，做出Arrhenius图（lnD vs 1/T），活化能从斜率中提取——这条斜线的线性度本身就是在验证你的模拟是否捕捉到了单一机制主导的扩散行为。如果Arrhenius图出现了弯折，那说明在中温区间可能有机制切换。

有一个实操层面的温度梯度策略。从高温往低温做的好处是——高温体系的原子迁移快、统计收敛容易——先跑出高温端几个温度点的可靠数据，然后用这些数据外推低温端的扩散系数量级，以此判断低温模拟需要多长的采集时间窗口。直接按低温开跑，在没有数量级概念的情况下很容易低估所需模拟时长——比如一个活化能1.0eV的过程，温度从1000K降到500K，扩散系数大约降了5个数量级，这意味着低温下的MSD曲线爬到能拟合线性段所需的模拟时长也需要相应增加5个数量级。

统计层面，扩散系数的标准误差来自两个独立的方向：一是时间序列的统计波动——均方位移曲线本身的线性拟合误差；二是系综平均的不足——单个模拟盒子里空位数量有限，稀缺陷浓度的统计代表性天然不够。前一个问题靠加长模拟时间和多次重复来解决，后一个问题需要有意识地构造多个独立初始构型，每个构型里空位的微环境不同。经验是5~8个独立构型搭配每条轨迹不少于2ns的采集时间，得到的扩散系数和活化能的可重复性基本能通过同行评审的审视。
晶界扩散是一个在纳米晶材料中无法回避的补充话题。当晶粒尺寸降到纳米级别时，晶界体积分数急剧增大，晶界扩散对总扩散率的贡献在Arrhenius图上表现为低温区的激活能显著低于体扩散激活能。做晶界扩散的MD模拟需要构造一个双晶模型——两个不同取向的晶粒夹一个晶界——然后分别追踪晶界区域内和晶粒内部的原子的MSD。晶界区域的原子识别可以通过共近邻分析（CNA）或者原子体积（Voronoi体积）来实现——晶界原子和体相原子在这两个度量上有明显的数值分界线。晶界扩散系数通常比体扩散高一到三个数量级，这个差异在纳米晶材料的使用温度下是决定性的。

https://www.keyanxueshu.com/category/md/