LAMMPS 这个坑我算是踩出经验了。以前做高熵合金的组织演化最怕的就是别人问你“退火之后结构到底变成什么样了”。实验上要做长时间扩散退火、然后再做 EBSD 或 TEM成本高、周期长而且很多中间状态根本抓不住。后来我把整套分析挪到 LAMMPS 上用分子动力学做退火模拟至少能先把“成分-温度-时间”这三个变量对结构的影响摸一遍再回头指导实验。这一讲就用 FCC 结构的高熵合金 CoCrCuFeNi 做例子把模型怎么搭、退火流程怎么设计、结果怎么分析讲清楚适合已经有 LAMMPS 基础、想真正上手做合金退火模拟的同学参考。这篇内容我尽量按实际操盘的顺序来讲先想清楚这个模拟要回答什么问题再决定用多大的盒子、选什么势函数然后才是退火工艺怎么在 LAMMPS 里实现。最后我会把常见的报错、结果异常和排查思路一并列出来省得你在网上翻半天找不到答案。1. 案例预设为什么选 CoCrCuFeNi 做退火模拟1.1 高熵合金和 CoCrCuFeNi 的选择理由高熵合金(HEA)这个概念已经火了很多年核心就是打破传统“一种主元微量合金元素”的思路让多种元素以接近等原子比混合。CoCrCuFeNi 是面心立方(FCC)高熵合金里非常典型的一个体系五种元素的原子半径差不算大容易形成稳定的 FCC 固溶体实验上也有很多文献可以直接对照。但对分子动力学模拟来说这个体系有个很有意思的点Cu 和 Fe、Co、Ni 之间混溶度并不好。也就是说虽然初始结构可以做成均匀的 FCC 固溶体但在退火过程中Cu 往往会偏聚、析出甚至向表面或晶界方向迁移。这种“名义上是单相高熵合金退火后却出现局部成分起伏”的行为正是很多实验文章里讨论的焦点。用 MD 来实现它既不需要搭复杂的界面模型也不需要预设缺陷只要给一个随机固溶体初始结构再加上合适的退火工艺就能看到元素再分布的趋势。另外一点CoCrCuFeNi 的 FCC 晶格常数约在 3.553.62 Å 之间具体数值取决于成分和温度。这个范围对 EAM 势函数非常友好直接建盒子、能量最小化就能跑起来。对刚接触高熵合金模拟的同学来说这是一个上手的绝佳体系难度适中但又能引出很多深层的物理问题。1.2 退火模拟在分子动力学里对应什么实验上的退火是把材料加热到某一温度保温再缓慢冷却目的是消除缺陷、均匀化成分、调控析出相。MD 里的退火本质上也一样只不过时间和空间尺度被压缩了很多。我们通常没法像实验那样退火几小时只能把过程压缩到纳秒量级用“更高温度 更高冷却速率”来换取足够明显的原子迁移。所以做 MD 退火模拟之前心里一定要清楚我们关注的是退火过程中的物理趋势而不是精确复现某一炉次的实验曲线。比如升温到哪个温度、保温多久、降温速率多快这些参数都有“模拟味”需要通过多组对照实验找出规律。这个思路我会在第 3 章展开。1.3 模拟整体流程速览这次模拟我把它拆成 4 个阶段每一步都有明确的目的阶段一建立随机的 FCC 固溶体模型超胞规模取 10×10×10 或 12×12×12对应 4000 到 6912 个原子。阶段二用共轭梯度法做能量最小化消除初始结构中可能存在的原子重叠再在目标温度下用 NPT 系综做预平衡。阶段三实施退火工艺包括升温、保温和降温三个环节。升温到 10001200 K 观察元素扩散保温一定时间后再逐步降温。阶段四分析输出数据比如径向分布函数(RDF)、原子团簇分析(CNA)、均方位移(MSD)判断退火后到底发生的是偏聚还是保持无序固溶体。计算成本方面4000 个原子的 EAM 模拟在普通工作站上用 8 核并行跑几百皮秒大概几小时就能完成。如果是更大尺寸、更长退火时间建议先想清楚统计需求再决定要不要上 GPU 或更大集群。这个后面会细说。2. 模型构建从晶格常数到随机固溶体2.1 晶格常数与超胞尺寸的确定在 LAMMPS 里建模第一步永远是确定晶格常数。对 CoCrCuFeNi 这类多主元合金有两种常见做法第一种经验加权平均。把各元素的纯 FCC 晶格常数按原子比例加权。Co 约 3.54 Å、Cr 在 FCC 结构下约 3.65 Å、Cu 约 3.61 Å、Fe 的 FCC 相约 3.64 Å、Ni 约 3.52 Å等摩尔平均下来大约在 3.583.60 Å 附近。这个方法最快但只能作为初值。第二种先用 NPT 预平衡“自动找”平衡体积。初始结构可以用上述估算值建立一个盒子然后在目标温度、零压下跑一段 NPT让盒子体积自己收敛到平衡值。EAM 势函数给出的平衡晶格常数和实验值会有一定偏差所以最终应该以势函数自身的预测为准。这也是为什么我强调势函数选择很重要后面专门讲。超胞尺寸方面10×10×10 的 FCC 盒子有 4000 个原子做 RDF、MSD、局部结构统计已经能看出明显趋势但如果要跟踪第二相析出形貌、晶界迁移这类空间关联较强的现象建议至少 20×20×20也就是 32000 个原子。退火模拟本身受有限尺寸影响比较大盒子太小析出的团簇会和周期性镜像相互作用导致结果失真。2.2 随机固溶体原子占位的实现这一步是最容易卡住新手的环节。LAMMPS 内置的lattice、create_atoms命令虽然能很方便地建立一个 FCC 单质晶体但并不能直接把不同元素随机放在各个格点上。你需要先产生一个“原子类型随机分配”的初始结构通常有三种实现路径路径一直接用 LAMMPS 命令建好单质 FCC 盒子再用文本处理脚本按原子 ID 随机改写 type 字段。做法是# 生成只有一种原子的 FCC 超胞 lattice fcc 3.58 region box block 0 10 0 10 0 10 create_box 1 box create_atoms 1 box之后用 Python 脚本读取生成的 data 文件把所有原子的 type 按等摩尔比随机替换成 15再写回 data 文件。这个方法完全可控适合熟悉脚本的同学。路径二用 ASE、pymatgen 等原子模拟库生成结构后再导出成 LAMMPS data 文件。这部分代码量很小ASES 里可以先把晶格建好再用random_replace之类的方式把 5 种元素随机填充上去。我个人更推荐这条路径因为生成的结构可以直接可视化检查一遍避免进 LAMMPS 之后才发现原子位置有问题。路径三先建一个包含了多种元素的超胞再用 LAMMPS 的set命令重新指定原子类型。实际操作中这种方式的灵活性不好因为你很难精确控制每个原子的类型分布。我一般不推荐除非你的占位规则特别简单。不管用哪种路径最后要确认 data 文件里每个原子都带有正确的新质量五种元素的比例严格等于 1:1:1:1:1。否则后面跑出来的成分偏析结果就没有意义了。2.3 势函数选用EAM 合金势高熵合金模拟能否可信90% 取决于势函数选得合不合适。对 FCC 金属和合金来说EAM嵌入原子法势是默认选择因为它能描述金属键和原子间的多体相互作用计算量又可控。针对 CoCrCuFeNi 这个体系一般能找到 Zhou 等人发展的 EAM 合金势它包含 Co、Cr、Cu、Fe、Ni 五种元素的交叉作用参数在 LAMMPS potentials 目录或者 NIST 势函数库里都能找到对应文件。选用势函数时有几个坑要特别注意。第一确认势文件里确实包含了你要用的所有元素并且元素顺序和 LAMMPS 里的原子类型一一对应。比如pair_style eam/alloy pair_coeff * * /path/to/CoCrCuFeNi.eam.alloy Co Cr Cu Fe Ni这行命令的意思是data 文件里的原子类型 1 对应 Co类型 2 对应 Cr以此类推。顺序一旦写错结果必然是错的而且很难从肉眼上发现异常。第二确认势函数预测的晶格常数、弹性常数和实验值偏差不大。你可以先建一个单质 FCC 晶体跑一个 0 K 能量最小化看晶格常数是否落在合理范围。如果纯 Co 的最小化结果比实验值偏大超过 2%那么这个势函数用于合金的可靠性就要打问号。第三如果现有 EAM 势精度不够可以考虑用机器学习势如神经网络势、深度学习势来替代。这类势函数对局部化学环境的描述更精细但需要大量的 DFT 训练数据普通课题组不一定具备条件。我的建议是先用 EAM 跑通流程确认物理趋势没有问题再考虑是否升级势函数精度。2.4 单位、边界条件和初始设置LAMMPS 常用的单位组合是metal和real。我这次用的是metal单位长度单位 Å能量单位 eV时间单位 ps温度单位为 K。EAM 势的参数一般都是按 metal 单位标定的所以除非你很清楚自己在做什么否则不要随便混用。边界条件默认用周期性边界boundary p p p这个对块体退火模拟是合理的。盒子尺寸的初始值按晶格常数 3.58 Å 来设置10 个晶格常数就是 35.8 Å可以直接用units metal boundary p p p atom_style atomic lattice fcc 3.58 region box block 0 10 0 10 0 10 create_box 5 box create_atoms 1 box注意这里我只是建了一个单质 FCC 盒子随后需要把类型替换成五种元素见 2.2。如果你已经用外部工具生成了 data 文件则不需要执行上述lattice命令改为read_data hea_CoCrCuFeNi.data还要设置每个原子类型的质量mass 1 58.933 mass 2 51.996 mass 3 63.546 mass 4 55.845 mass 5 58.693最后在正式跑之前用write_data输出一个中间文件或者直接跑一段 0 K 能量最小化检查一下。这样如果后面出问题可以迅速定位是势函数、初始结构还是工艺参数的问题。3. 退火工艺设计升温、保温和降温怎么设置3.1 预热与能量最小化别让原子重叠毁掉你的盒子读取初始结构后第一件事永远是能量最小化。LAMMPS 里常用共轭梯度法min_style cg minimize 1.0e-6 1.0e-8 1000 10000第一个参数是能量容忍度第二个是力容忍度后面是最大迭代次数。如果最小化过程中能量不收敛或者出现 NaN90% 的原因是随机占位时两个原子距离太近也就是产生了物理上不允许的原子重叠。遇到这种情况我一般会先检查 data 文件里最近邻原子间距再用write_dump导出结构到 VMD 里看一眼。最小化完成后还需要给原子赋一个随机的初始速度。这个过程要在目标温度下设置velocity all create 300 12345 mom yes rot yes12345是随机数种子可以随意指定。mom yes和rot yes会把系统整体的平动和转动动能扣除防止温度计算出现虚假偏移。3.2 高温平衡给扩散一个足够强的驱动力实验退火通常是在合金熔点以下 0.60.8 倍的温度进行目的是让原子扩散能够发生但又不至于让样品熔化。MD 模拟里我们往往会把温度取得更高一些因为在纳秒时间尺度内原子的迁移距离实在有限温度不够高就什么都看不到。对 CoCrCuFeNi 这个体系经验熔点大约在 1600 K 以上。我常用 10001200 K 作为退火温度。这个温度下FCC 结构保持稳定但原子已有足够的动能进行短程扩散。如果想更清楚地观察 Cu 的偏聚倾向可以先把系统在 1200 K 保温一段时间让成分起伏发展起来再降温“冻结”住结构。高温平衡用什么系综我建议先用 NPT恒定原子数、压力、温度跑 50100 ps让盒子体积在零压下充分弛豫得到一个真实的平衡密度。然后再切到 NVT恒定原子数、体积、温度进行正式退火因为体积固定后温度控制更稳定也方便后续做结构分析。如果你在退火过程中一直开着 NPT盒子可能会随着温度变化不断收缩或膨胀导致局部压力波动这对判断相变点会产生干扰。具体命令参考# 300 K 到 1200 K 升温100 ps 线性升到 fix anneal1 all nvt temp 300 1200 $(100*dt) run 100000 # 1200 K 保温 200 ps fix anneal2 all nvt temp 1200 1200 $(100*dt) run 200000这里的 $(100*dt) 是温度阻尼参数单位是时间。阻尼参数的经验值是每 1 ps 左右一个量级我通常设在 0.1 ps 到 0.5 ps 之间太小会导致温度波动大太大会让温度控制响应迟缓。3.3 阶梯降温和连续降温两种退火策略的比较保温结束后进入降温阶段。LAMMPS 里fix nvt本身就支持从Tstart到Tstop的线性温度斜坡也就是你告诉它起始温度和终止温度它会在这段时间内均匀地升降温度。比如# 从 1200 K 经过 200 ps 降到 300 K fix cool all nvt temp 1200 300 1000 run 200000这里的1000仍然是阻尼参数单位是时间步的倍数需要换算成时间。我习惯先写成固定时间步长下的绝对值方便后期调整。不过连续线性降温在实际模拟中有一个问题温度变化太快时系统的相变可能来不及响应导致你看到的是“过冷液体”而非真正的退火组织。更稳妥的办法是阶梯降温把温度范围切成几个区间每个区间保温一段让结构充分松弛后再继续降# 第一段1200 K 降到 1000 K保温 100 ps fix cool1 all nvt temp 1200 1000 500 run 50000 fix hold1 all nvt temp 1000 1000 500 run 100000 # 第二段1000 K 降到 800 K保温 100 ps fix cool2 all nvt temp 1000 800 500 run 50000 fix hold2 all nvt temp 800 800 500 run 100000 # 第三段800 K 降到 300 K fix cool3 all nvt temp 800 300 500 run 100000阶梯降温的好处是每一步都给了原子足够的时间去弛豫降温过程中的瞬态效应更小。缺点是计算时间会明显增加。我在实际操作中会先用连续降温跑一版快速结果确认整体趋势没问题再用阶梯降温做精细模拟。3.4 时间步长和系综选择的实战建议对这个体系2 fs 的时间步长在低温下是安全的但在高温1200 K 以上或者快速升温阶段我建议改成 1 fs。原因是温度越高原子热运动越剧烈如果时间步长太大能量积分会出现明显漂移最终可能导致 LAMMPS 报出Bad global thermo energy之类的错误。你可以在一个输入脚本里直接设置timestep 0.001 # 单位 ps即 1 fs另外温度控制和积分器之间也存在搭配问题。LAMMPS 的fix nvt和fix npt内置了 Nose-Hoover 恒温器如果你同时还想用fix langevin做局部温度控制两个恒温器叠加会导致真正有效的温度状态很混乱我不建议同时使用。退火流程中统一用一个恒温器就好别贪多。最后提醒一个很容易被忽略的点NVT 下盒子体积固定温度变化时压力会随之起伏。如果在降温后期你发现压力出现很大的负值或正值说明这个温度下平衡体积已经和初始体积差得太多应该在中途改用 NPT 再平衡一段而不是硬扛着让系统处在高静水压状态。4. 结果分析怎么知道退火之后到底发生了什么4.1 热力学曲线的判读模拟结束后先看热力学输出。LAMMPS 的thermo命令每 N 步输出一行时间、温度、势能、总能量、压力等信息我会固定设置thermo 1000 thermo_style custom step temp pe ke etotal press vol升温阶段温度应该平滑上升势能随温度升高而增加保温阶段温度在目标值附近波动能量趋于平稳降温阶段能量和体积应逐渐下降。如果温度曲线出现明显的平台甚至倒挂说明系统正在发生相变比如局部熔化或者析出这时候需要结合结构分析来确认。另一个关键参数是体积或者密度。如果用 NPT 做平衡能直接看到体积随温度的变化热膨胀曲线如果出现明显的转折往往对应相变点。这一步可以帮你确定体系的熔点或者固溶线对后续退火工艺选择很有指导意义。4.2 径向分布函数RDF和结构因子RDF 是判断短程序最直观的工具。FCC 晶体的 RDF 第二峰是一个典型的双峰劈裂结构液体或非晶态则只有一到两个宽峰没有长程序特征。LAMMPS 里计算 RDFcompute rdf all rdf 100 fix avgrdf all ave/time 100 10 1000 c_rdf[*] file rdf.dat mode vector跑到平衡后rdf.dat会给出不同原子对之间的距离分布。我一般同时输出总 RDF 和分元素对的 RDF比如 Cu-Cu 的 RDF 和 Fe-Fe 的 RDF。这样就能看出退火后哪些元素更容易聚在一起。如果 Cu-Cu 的第一峰明显比随机固溶体状态更高、更尖锐这说明 Cu 已经发生了偏聚。这个信息实验上不容易拿到但对理解高熵合金的稳定性非常重要。4.3 局部结构分析CNA 和 PTMRDF 是全局平均意义上的结构信息想看看每个原子所处的局部环境到底是不是 FCC就要用 CNA公共近邻分析或者 PTM多面体模板匹配。LAMMPS 里的 CNA 命令compute cna all cna/atom 3.2 fix avgcna all ave/atom 10 10 100 c_cna[1] file cna.out这里的 3.2 是近邻截断距离单位是 Å需要根据第一近邻距离来调整。具体截断值可以通过 RDF 第一峰谷底位置来估计一般取峰谷最低点附近。CNA 会输出每个原子属于 FCC、HCP、BCC 还是其他非晶态结构进一步统计即可得到退火后 FCC 比例是否提高有没有新的 HCP 堆垛层错出现是否存在 BCC 相的局部临界核晶界或无序区的占比。如果你考虑析出也可用compute cluster/atom结合成对截断来判断不同元素的团簇尺寸分布。4.4 MSD 和扩散行为辅助判断退火过程中元素扩散是核心物理过程最好顺手计算一下均方位移MSD看每种元素的迁移能力。LAMMPS 里compute msdco all msd fix avmsd all ave/time 100 10 1000 c_msdco[4] file msd.dat这里的 c_msdco[4] 是总均方位移或者你可以按元素分组分别计算 Co、Cr、Cu、Fe、Ni 的 MSD比较它们的扩散快慢。一般情况下Cu 在高熵合金里扩散较快这和它容易偏聚的现象是相互印证的。通过 MSD 曲线的线性段斜率可以估算扩散系数D MSD / 6t。不过要注意如果退火过程中系统发生相变或者成分偏析扩散系数会在不同时间段有明显变化这时单纯用一个 D 值已经没有太大意义最好分段分析。5. 常见问题与排查技巧5.1 原子重叠和“missing atom”报错新手最常遇到的错误是 LAMMPS 在最小化阶段报出原子缺失或者 bonds 错误。这个体系用的是原子类势函数没有 bond所以最常见的表现是能量无法收敛、力分量出现 NaN。排查步骤检查 data 文件里最近邻原子距离是否小于势函数的截止半径通常 EAM 的截止在 5 Å 左右如果两个原子间距小于 1.5 Å基本就是重叠。用write_dump导出结构在 VMD 里查看是否有原子位置重叠。如果确认是随机占位导致的重叠可以先把盒子稍微扩大一点或者用高对称占位算法重新生成初始结构。5.2 势函数与结构不匹配导致的异常熔化我踩过最深的坑是换了势函数后体系的熔点比实验值偏高了 300 多 K。这意味着我按实验退火温度设计的模拟在势函数里其实还处于很低的同源温度原子扩散很慢结果什么都看不到。所以正式做退火前先用 NPT 扫一遍温度看看势函数预测的熔点大概在哪。方法很简单从 300 K 升温到 2000 K每 100 K 为一个节点看体积和 RDF 是否发生突变。根据这个结果再决定退火的具体温区比自己拍脑袋定温度可靠得多。5.3 降温速率敏感性和相变滞后MD 的降温速率通常都在 1 K/ps 以上比实验快了十几个数量级因此观察到的结构转变温度一般低于实验值。降温越快过冷度越大析出或者有序化的临界温度就越低。如果你要做严格定量分析建议至少跑三组不同的降温速率看结果是否收敛。实际操作中我还喜欢用“快速淬火短时保温”的组合先快速从高温降到目标温度再在该温度下保温足够长时间观察结构是否继续演变。这样既节省了均匀降温的时间又能捕捉到亚稳相的长大过程。5.4 有限尺寸效应与统计取样最后要说的坑是盒子尺寸。4000 个原子对于元素偏聚的趋势分析基本够用但如果你想统计析出相颗粒的尺寸分布或者观察晶界上的偏聚就很有可能被周期边界坑到。析出相一旦长大到盒子尺寸的 1/3 以上它的长大行为就会受到镜像相互作用的影响结果不可信。遇到这种问题我通常会把超胞扩大到 15×15×15 或者更大再做一次同样的退火流程看关键指标比如 Cu 团簇尺寸、FCC 比例是否发生显著变化。如果两次结果一致说明有限尺寸效应可以接受如果不一致那只能老老实实用更大的盒子。我在实际跑 CoCrCuFeNi 退火时还有一个体会这个体系的结构演化非常依赖初始随机占位的种子。换句话说不同随机种子生成的同一成分固溶体退火后局部成分分布可能差异很大。建议同一种子体系至少重复 3 次独立退火统计平均之后再下结论这样才能把“随机涨落”和“真实物理过程”区分开。模拟这条路没有太多捷径参数扫得越多、重复次数越多结果才越经得起推敲。