我印象很深的一次翻车经历用ATOMSK生成铝的堆垛层错模型在OVITO里看ABC序列错位得清清楚楚扔进LAMMPS一算层错能出来310 mJ/m²文献值才160左右。一开始怀疑势函数不对换了三个EAM势依然偏高最后才发现模型本身带了两个层错面能量差忘了除以二。这种“每一步都好像没问题连起来结果就离谱”的事在面缺陷模拟里实在太常见。这篇就把我从ATOMSK建模、LAMMPS算层错能到OVITO可视化验证的完整流程以及所有踩过的坑一起整理出来适合正在做FCC金属层错能计算、被数值对不上文献折磨的人。1. 先搞清楚模型里要造的是什么堆垛层错的几何与能量1.1 面缺陷的本质ABCABC里的秩序破坏FCC晶体沿111方向的密排面正常堆叠顺序是ABCABCABC……无限重复。所谓堆垛层错就是某个环节的堆叠顺序出现了异常。最简单的理解方式把每一层原子看成一层“地板”正常FCC是三种位置循环铺一旦某两层之间的相对位置错了就破坏了循环形成一个二维的面缺陷。按照缺陷形态层错分两种。内禀层错ISF相当于在某个位置抽掉了一层堆垛序列中会出现局部的ABAB特征外禀层错ESF则是插入了一层形成ABA的局部特征。几何上二者都表现为夹在FCC区域中间的一小段HCP密排结构这也是后面OVITO里用CNA分析时层错面上的原子会显示成HCP颜色、而不是FCC颜色的根本原因。层错不是凭空存在的它总是和位错联系在一起。Shockley不全位错的扫过区域内部就是一片内禀层错扩展位错的两条Shockley部分位错之间也夹着一片层错带。这也是为什么材料学家这么关心层错能——它决定了扩展位错宽度、交滑移难易、孪生倾向对塑性变形机制影响很大。1.2 为什么用模拟算层错能能量差公式的来历实验测层错能的方法不少比如透射电镜下测量扩展位错节点宽度、层错四面体尺寸等但误差范围动辄百分之三四十。原子模拟的好处是可以直接对比“有层错”和“无层错”两个体系的总能量差再除以层错面积得到单位面积的层错能。标准公式是[ \gamma_{SF} \frac{E_{SF} - E_{perfect}}{n \cdot A} ]其中E_SF和E_perfect分别是层错模型和完美晶体模型的总能量A是层错面面积n是模型中包含的层错面数量。这个n极其关键我后面会反复强调。很多人照文献公式抄作业抄了个除以A结果自己模型里因为周期性边界条件出现了两个层错面算出来正好是文献值的两倍。层错面积A怎么取如果你的模型z轴沿[111]方向xy平面就是(111)密排面层错面是和xy平面平行的完整平面A就等于盒子的Lx乘以Ly。单位方面LAMMPS的metal单位制下能量是eV、长度是Å算出来的层错能单位是eV/Ų。1 eV/Ų约等于16.02 J/m²换算成常用单位就是1 eV/Ų 16021.8 mJ/m²。这个换算关系每次都要写对否则结果和文献差好几个数量级。1.3 建模的核心目标两个体系除了层错没有任何区别这个原则是整个建模工作的纲领。完美晶体模型和层错模型必须具备完全一致的原子数、盒子尺寸、原子类型和势函数参数唯一区别是后半部分晶体是否沿特定方向平移了层错矢量。任何多余的删除原子、修改盒子尺寸、改变原子数量都会让能量差里混入别的贡献层错能的计算就失去意义。所以建模型时我不推荐“建完层错再删掉边界原子做真空层”这种操作。表面能、镜像作用、原子缺失带来的弛豫都会淹没层错面的能量信号。正确做法是保持三维周期性边界条件让两个模型在完全相同的边界条件下做能量对比误差项自然抵消。2. ATOMSK建模从命令写法到几何原理2.1 方向选择为什么非要把[111]摆到盒子某个轴ATOMSK建层错的思路是先生成一个完美FCC超胞然后把某一层以上的原子整体平移一个层错矢量。要做这件事首先必须保证密排面(111)和盒子的某个坐标面平行这样“原子层”才能和“z坐标轴”对应起来。我常用的命令是atomsk --create fcc 4.05 Al orient [1-10] [11-2] [111] -duplicate 20 20 20 Al_perfect.cfg这条命令生成一个晶格常数为4.05Å的FCC铝超胞盒子的x、y、z轴分别对应[1-10]、[11-2]、[111]三个晶向duplicate参数表示沿三个方向重复的周期数。如果把FCC结构想象成一摞盘子[111]方向就是穿过盘心的垂直方向xy平面正好是盘面——也就是密排面。这里必须提一个极易踩的坑orient后面的三个晶向要满足右手定则。很多人写orient [111] [1-10] [11-2]看着没问题实际上[111]叉乘[1-10]得到的是[-1-12]和[11-2]正好反方向ATOMSK会报错或者生成一个镜面对称的盒子。我在新旧版本ATOMSK上都遇到过这个问题换个顺序或者把[11-2]改成[-1-12]通常就能解决。模型建好以后应该把盒子坐标输出到屏幕上确认三轴的夹角是90度再开始duplicate。2.2 duplicate参数与盒子尺寸层错面积和镜像距离的平衡duplicate的数量不是随便定的。x和y方向的重复次数决定了层错面积A。比如铝的晶格常数4.05Å沿[1-10]方向一个周期长度是4.05乘以根号2再除以2大约2.86Å重复20次就是约57Å沿[11-2]方向一个周期大约是4.05乘以根号6除以2约4.96Å重复20次约99Å。这样层错面积大概有5600Ų足够让能量差的信号远高于数值涨落。z方向重复次数决定了两个相邻层错镜像之间的距离。由于周期边界条件一个层错模型里实际有两个层错面如果z方向太短两个层错面靠得太近彼此的弹性交互作用会污染能量。我的经验是z方向至少要有20层以上的(111)原子层也就是重复次数在15到20以上比较好。模型总原子数通常在几万个到十几万个LAMMPS跑起来毫无压力没必要为了省计算量把盒子做太小。2.3 制造层错平移半块晶体而不是删掉一层ATOMSK新版本虽然有专门的堆垛层错生成选项但不同版本的参数格式有差异而且很多人用的还是老版本。我最推荐的通用做法是用-select配合-shift手动平移原理透明、版本兼容性好出了问题也知道去哪查。操作分三步。第一步先查一下生成盒子沿z方向的总长度找一个位于原子层间隙处、接近盒子一半高度的z坐标作为切割面。第二步选中该平面以上的所有原子。第三步把选中的原子整体平移层错矢量R。对FCC内禀层错R通常取Shockley不全位错对应的位移矢量也就是a/6112类方向。当盒子x、y、z方向是[1-10]、[11-2]、[111]时这个位移矢量正好落在xy平面内。对铝来说a4.05ÅR的模长是a/6乘以根号6算出来约1.6536Å。完整命令如下# 生成完美晶体 atomsk --create fcc 4.05 Al orient [1-10] [11-2] [111] -duplicate 20 20 20 Al_perfect.cfg # 平移z坐标20.5以上的原子注意这个20.5要按实际盒子高度调整 atomsk Al_perfect.cfg -select above 20.5 -shift above 0.0 1.6536 0.0 -wrap Al_sf.cfg # 转成LAMMPS data文件 atomsk Al_sf.cfg lmp注意这里-shift above后面的位移向量是0、1.6536、0因为层错矢量恰好沿着y轴方向。如果你用的盒子方向和我不同向量分量要重新投影。-wrap选项会把越出盒子的原子折回盒子内部这是必要的否则后续LAMMPS读入时会因为坐标越界报错。还有一个关键点这种平移方案不会删除任何原子也不会改变原子层之间的垂直距离所以不需要再手动清理重叠原子。平移后上半部分和下半部分在切割面处重新叠放形成层错同时顶层原子被wrap到底层区域在盒子顶底边界处也会形成一个等价层错。2.4 模型自检用“数”而不是用“看”来验证建模完成后别急着进LAMMPS先做几个最基础的检查。第一原子总数是否和完美晶体一致。理论上平移不改变原子数但如果你的ATOMSK版本在-wrap后悄悄合并了重复原子总数就会变化这时候能量差公式里的体积项会出问题。第二检查层错面对应的原子层数是否正确。我习惯用atomsk自带的-a选项或者直接看输出文件统计每个z坐标附近的原子数量。FCC的(111)面是三层一个周期如果平移后某两层的面内位置完全重合而其他层不重合说明层错矢量选错了产生的可能不是单一内禀层错。第三用OVITO快速看一眼。载入Al_sf.cfg按z坐标给原子着色你会看到沿z轴方向原子层颜色周期性变化唯独在某个z位置附近上下层的水平偏移规律被破坏——那里就是层错面。这个检查不用做CNA肉眼足够确认方向对不对、层错面是否平行于xy平面。如果层错面是斜的说明盒子方向设置有问题回去查orient的三个晶向。3. 进入LAMMPS格式、边界与弛豫的设计3.1 从cfg到data你以为转完格式就完了ATOMSK转LAMMPS data文件很简单一条atomsk Al_sf.cfg lmp就搞定。但生成data文件不代表万事大吉有几个细节必须打开文件确认。首先是原子类型。纯铝只有一个原子类型data文件里atom_types应该是1对应的Masses段必须给出质量。如果你用的势函数文件里不带质量而data文件里又没写LAMMPS读到pair_coeff时会直接报错。其次是晶格方向信息data文件里会有一个xlo xhi ylo yhi zlo zhi的盒子范围务必确认z方向盒长和ATOMSK输出的盒子一致否则后面算层错面积时A取错结果全错。还有一个容易忽视的坑LAMMPS读入data文件后原子编号顺序和ATOMSK里的顺序可能不一致。这是正常的但如果你后面想用group按原子序号选层错附近的原子就需要小心。我更推荐按region选择因为在层错建模中层错面的位置和盒子坐标的关系是明确的用region加group按坐标范围选择比依赖编号靠谱得多。LAMMPS的安装也顺带提一句。无论你用源码编译还是conda装好的二进制版本确认mpirun能正常起并行任务再开始建模。层错能计算本身很小几万原子单核也能跑但如果你后面要做大体系的扩展位错模拟MPI并行是必须的。编译LAMMPS时记得把需要的势函数相关包、额外包都打开否则读到某个pair_style时才发现没编译就尴尬了。3.2 周期性边界和“双面层错”效应LAMMPS的in文件里边界条件默认是boundary p p p也就是三个方向都是周期性。层错能计算要求三个方向都保持周期性这样才能消除表面能干扰。很多人会想既然周期性会在z方向产生两个层错面那我是不是把z方向设成ffixed或者mshrink-wrapped就能只留一个面这个想法我在早期也动过。实际做了以后发现固定边界会在z方向引入两个自由表面表面的弛豫和重构带来的能量变化比层错能本身还大算出来的数值完全不可用。正确的做法是接受两个层错面的事实在能量差公式中除以2倍的层错面积。这不是妥协而是周期性条件下最干净的处理方案。还有一点要提醒如果你用OVITO打开data文件默认只显示一个盒子可能只看到一条层错打开周期性副本功能Show Periodic Images后才能看到上下镜像处隐藏的第二条层错。所以“OVITO里只看到一条层错”是正常的能量计算时该除以2还是要除以2。3.3 先让原子松弛到局部能量最低再采能量静态计算直接run 0得到的是建模样式下的瞬时能量。ATOMSK生成的完美晶体虽然几何上是理想FCC但因为你用-wrap做了坐标折叠原子位置可能出现微小的坐标不匹配直接run 0的能量会包含这部分人为涨落。我建议在计算能量前先做一次能量最小化让体系弛豫到局部能量最低状态。对纯单质晶体几十步共轭梯度就能收敛成本极低。minimize后完美晶体和层错模型都处于各自的局域势能极小值点能量差才反映真实的层错能。但这里也有一个两难如果让体系充分弛豫到全局最小某些势函数下层错面可能会迁移、湮灭尤其当层错能很低、体系又很小时层错面甚至会自动消失。这说明该势函数对这个体系的描述不够好或者模型尺寸不足。遇到这种情况可以改用固定部分原子、只让界面附近原子弛豫的受限优化方案或者换势函数验证。4. LAMMPS计算层错能从in文件到单位换算4.1 势函数是整个计算里最大的“自由参数”层错能对势函数极其敏感。同样是铝用不同EAM势算出来的内禀层错能可以从120 mJ/m²到250 mJ/m²差异接近一倍。这不代表模拟不可信而是说明不同势函数拟合数据集不同对偏离平衡构型的缺陷描述能力不同。选势函数时先查文献里这个势函数是否拟合过层错能或广义层错能GSFE。很多EAM势的拟合目标主要是晶格常数、弹性常数、空位形成能、内聚能等层错能只是“顺带”预测的物理量。如果你要研究层错相关的机制尽量选择在论文中明确报告了层错能、孪生能的势函数。以铝为例Mishin等人发展的EAM势在层错能、表面能方面校准得较好是面缺陷模拟的常用选择而某些早期EAM势的层错能明显偏高用它们算SF你能得到“铝层错能特别高”的错误结论。另外势函数文件的格式必须和pair_style匹配。我用的是eam/alloy对应势文件通常是Al99.eam.alloy这类格式。用错eam/fs格式或者文件路径放错LAMMPS会直接报错或者给出一堆乱数值。4.2 in文件骨架可以直接套用下面这个in文件我习惯命名成in.sfe完美晶体和层错模型共用只改read_data读取的文件名units metal atom_style atomic boundary p p p read_data data.sf pair_style eam/alloy pair_coeff * * Al99.eam.alloy Al neighbor 2.0 bin neigh_modify delay 10 every 1 check yes minimize 1.0e-10 1.0e-10 2000 5000 run 0 variable pe equal pe print PE ${pe}几个值得注意的细节。neighbor宽度2.0是经验值必须大于势函数截断半径和模拟盒子最短边长之半的差否则部分原子对会被漏算。minimize的前两个参数是能量和力的收敛容差建议至少到1e-10量级因为层错能是几十上百个meV级别的能量差容差太粗会让你无法区分真实信号和数值噪声。run 0不推进时间积分只是让LAMMPS完成一次势能计算输出体系总能。完美晶体跑一遍层错模型跑一遍分别记录PE。注意两个模型必须用完全相同的in文件逻辑、相同的势函数参数这是能量差法能成功的前提。4.3 从eV到mJ/m²一步步算给你看假设一个层错模型盒子沿x方向[1-10]的长度Lx57.2Å沿y方向[11-2]的长度Ly99.2Å层错面积A57.2×99.25674Ų。完美晶体minimize后体系总能PE_perfect-4000.00eV层错模型PE_sf-3999.70eV数值仅为示意。能量差ΔE0.30eV。周期性条件下模型含两个层错面n2所以每层错面能量贡献为0.15eV。层错能γ0.15/56742.64e-5 eV/Ų。乘以16021.8换算成mJ/m²结果是0.42... 等等这个结果太小说明我的示意数值取得太离谱实际Al层错能应该在150 mJ/m²左右反推ΔE大概是0.3eV、面积几百Ų。为了不误导读者我应该把示意数值调合理一点。实际上铝晶格常数4.05Åx方向[1-10]一个周期长度a/√22.864Å20个周期57.3Åy方向[11-2]一个周期长度a√6/24.961Å20个周期99.2Å面积约5685Ų。如果文献值Al的ISF约150mJ/m²换算成eV/Ų就是150/16021.80.00936eV/Ų。单个层错面对应的能量是0.00936×568553.2eV两个面加起来ΔE≈106.4eV。这个数量级就对了。示意计算时要拿准量级否则读者按你的流程算完发现差2个数量级又得排半天错。4.4 ESF外禀层错的额外说明外禀层错的建模和ISF略有不同需要在正常堆垛中插入一层原子不能简单通过平移半块晶体来实现。ATOMSK建模外禀层错时通常先选择要插入的原子面复制一份再整体平移然后用-remove-doubles清理重叠原子。这个操作比较繁琐我建议第一次做的时候先用小盒子跑通流程确认原子数和结构无误后再放大到正式尺寸。外禀层错能一般大约是内禀层错能的1.5到2倍这个关系可以作为验证模型的粗略参考。5. OVITO可视化验证别停留在“看起来很对”5.1 CNA的根本逻辑和截断半径陷阱OVITO里最常用来看层错的是Common Neighbor Analysis。CNA会分析每个原子周围12个近邻的拓扑连接方式然后把环境分为FCC、HCP、BCC、ICOSAHEDRAL和OTHER几类。在FCC完美晶体里所有原子都被识别为FCC在层错面上由于上下堆垛变成ABAB原子环境具备HCP特征于是被识别为HCP。所以一张图上哪个位置出现了一条HCP颜色的“带”那里就是层错面。CNA有个参数叫截断半径这个值设置不当会得出完全错误的结论。截断半径太小最近邻统计不全FCC原子会被误判为OTHER截断半径太大第二近邻也被算进来FCC和HCP就没法区分了。我的经验是对FCC金属截断半径取最近邻距离的1.05到1.2倍比较合适。以铝为例最近邻距离是a/√22.86ÅOVITO的CNA默认截断半径通常是0.0需要手动填3.0左右。判据是在完美晶体上做CNA所有原子都显示为FCC稍微调大截断半径如果开始出现大量OTHER原子说明取值偏大了回到出现过拟合之前的节点。5.2 在OVITO中同时确认层错面位置、数量和方向打开OVITO后我推荐的检查流程是载入data文件或dump文件先按粒子类型着色确认原子无重叠、盒子边界整齐。添加Common Neighbor Analysis把FCC设为灰色、HCP设为亮色其余结构设为浅色或隐藏层错面会以一条亮色平面凸显出来。点击“Create Bond”之外在显示设置里打开周期性副本设成3×3×1或1×1×3你会看到层错面周期性重复出现。数一数在一个盒子范围内层错面到底有几个。ISF模型应该看到两条亮色平面并且相互平行。使用Slice工具沿z方向取一个厚度为0.5Å的薄片配合CNA着色直接读堆垛顺序。切面的x轴和y轴分别是[1-10]和[11-2]可以确认层错面确实横贯整个盒子的xy方向。如果切面上看到的层错线不是一条贯穿盒子的直线而是歪歪扭扭或者中断多半是模型构建时层错矢量方向不对或者盒子方向的选取让层错面与xy平面不平行。此时返回ATOMSK检查orient和shift方向不要继续往LAMMPS流程走否则后面所有计算结果都是错的。5.3 一个常见的误判你把HCP带当成了层错初学者经常在OVITO里看到大片HCP颜色以为生成了大片层错。实际上当模型中含有孪晶界面、晶界、表面时这些界面处的原子CNA也会显示为HCP或OTHER。更常见的场景是势函数弛豫之后层错面附近的原子不会全都整齐地保持HCP拓扑部分原子由于畸变会变成OTHER只有层错面中央一小圈干净的区域显示为HCP。所以别看到一小条HCP带就认为模型做对了一定要配合坐标切片确认堆垛序列确实在某一层发生了变化、其他地方完全正常。6. 全流程避坑实录六个高频翻车点下面这张表是我这些年摸爬滚打的浓缩版先给结论再逐项拆解。坑典型症状解决方向层错面数量没数清层错能高出一倍周期性镜像会产生两个面能量差除以2A层错矢量方向或符号错CNA看不到规则层错面核对R的分量手动算位移投影CNA截断半径乱填FCC显示成HCP或OTHER取最近邻距离的1.05-1.2倍势函数没校准层错能结果与文献差30%以上查势函数拟合文献换势验证盒子z太小层错能依赖盒子尺寸增加z方向重复次数做尺寸收敛测试弛豫过头层错消失minimize后HCP带消失用受限弛豫或增大盒子重新测试6.1 除以2A还是除以A周期性镜像的隐藏层错这是层错能计算中最多人翻车的坑没有之一。用平移法建层错模型时内部界面有一层错盒子顶底周期性边界处又会出现一层错。OVITO里默认显示下你可能只看到一条亮色带但周期性物理图像里就是两条。文首提到我那次310 mJ/m²的错误结果就是只除了一个A。解决办法很简单建模后用OVITO开周期性副本数清楚层错面数量再决定公式里n取1还是2。6.2 层错矢量方向搞错能量对不上但OVITO可能有不出来层错矢量R只有正确方向才能产生稳定的内禀层错。如果你把符号写反相当于把上半部分沿反方向平移很多情况下生成的是另一条Shockley部分位错扫过的等价构型层错能数值不变但层错面的位置或者堆垛特征和预期不同如果分量算错可能生成外禀层错或者更复杂的两层错结构。我的排查方法是把R的三个分量按盒子基矢投影手算一遍再对比ATOMSK命令里的-shift数值确认投影后的位移量和理论值一致。6.3 同是EAM层错能差出一倍这件事势函数对层错能的影响远超多数人的预期。不同EAM势对FCC金属层错能的预测可以差30%到100%这不是代码问题是势函数本身的性质。做计算前先看势函数论文里有没有报告SFE和GSFE曲线。没有的话至少要和实验值或者其他势的计算结果做交叉验证。我习惯在正式跑大模型前用同一个势函数算一次完美晶体和一次最小层错模型先确认层错能落在合理区间——这个步骤十分钟就能完成能替你省下后面整整一天的纠结。6.4 盒子z方向不够大层错能和自己的镜像在“打架”如果z方向原子层数太少两个层错面之间的距离只有十几埃它们之间的弹性交互作用会让层错能偏高或偏低并且结果会明显依赖盒子尺寸。判断方法是做一个尺寸收敛测试把z方向重复次数从10、20、30逐渐增大算三次层错能看数值是否趋于稳定。一般到20层以上变化就很小了但不同材料、不同势函数表现不一样以你的实际收敛曲线为准。6.5 弛豫后层错面“消失”了什么minimize完成后发现OVITO里原本清晰的HCP带不见了这一般有两种情况。第一种势函数对这个层错构型本身不稳定体系自动翻回完美FCC说明该势函数描述的层错能极低甚至为负不适合用来研究这个缺陷。第二种原子的局域弛豫让层错面附近结构发生微小畸变CNA的截断半径无法捕获正确的拓扑分类需要微调CNA参数再看不代表层错真的消失了。区分二者的方法是看minimize前后总能量变化率如果层错模型和完美晶体模型最终能量几乎一样那就是真消失了。6.6 ATOMSK版本差异导致的命令格式问题不同版本的ATOMSK对-select、-shift、-stacking-fault的参数解析并不完全一样。老版本里-shift之后的向量顺序可能是笛卡尔坐标新版本可能加入了更多选项。我的建议是第一次在某个版本的ATOMSK上建模先用最小盒子和纯FCC铜或铝跑通流程用OVITO确认层错面正确后再扩大尺寸正式建模不要一上来就追求大模型。就我个人的使用体会而言ATOMSK这种脚本化的建模工具最大的价值在于可重复性命令写清楚模型就能被任意复核。层错能计算本身不复杂难点全在“模型到底建对没有”这件事上。如果你刚开始接触先用最小体系把“完美晶体层错模型”的能量差跑通再逐步扩大盒子、换势函数、做收敛性测试。后续如果想进阶还可以在这个模型基础上做广义层错能曲线、构建扩展位错、研究层错与溶质原子的交互建模思路都是一脉相承的。