高熵合金嘴里喊了好几年了但真正把 FCC-CoCrCuFeNi 这个五元体系在 LAMMPS 里跑起来再设计一轮说得过去的退火模拟还是有不少门道。这篇我直接结合自己最近跑的一组算例从建模、退火路径设计、in 脚本写法到后处理判据和常见翻车现场完整拆一遍。适合已经跑过几轮 LAMMPS 基础任务、想往高熵合金退火方向深入的同学也适合刚装好 LAMMPS、想对照着做第一个真实课题的新手。我用的 LAMMPS 版本是 2022 年之后的 2.x 版本MEAM 势文件来自官方 examples 库中的 FeNiCrCoCu 参数整体流程不依赖特殊插件只要你把 MEAM 包编译进 LAMMPS脚本基本可以原样跑通。1. 建模这步别偷懒FCC 晶格与随机占位1.1 晶格常数和盒子尺寸怎么定FCC-CoCrCuFeNi 是典型的无序固溶体结构所谓的“FCC”指的是金属原子整体占据面心立方点阵位置但 Co、Cr、Cu、Fe、Ni 五种原子以近似等摩尔比随机分布在格点上。因此在建模阶段核心不是“画一个晶体”而是决定两件事点阵常数取多少盒子里放多少个原子。实验上这个体系的晶格常数大约在 3.56 到 3.61 埃之间不同热处理状态会有小幅波动。我做模拟时取了 3.58 埃这是文献里比较常用的中间值。需要注意的是这个值决定了初始原子间距如果选得过小原子初始重叠严重后面升温时能量会直接爆炸如果选得过大初始结构过于稀疏后续弛豫要花更多时间才能让体系恢复到合理密度。盒子尺寸我选了 10 × 10 × 10 个晶胞对于一个 FCC 晶胞每个晶胞包含 4 个原子所以总共 4000 个原子。这个体量不算大单核也能跑但用 4 到 8 核并行时效率明显提升适合先验证流程。如果后续要观察析出相或位错网络建议把盒子扩大到 20 × 20 × 20也就是 32000 个原子左右结果统计性会好很多。手册里反复强调的周期性边界条件在这里直接用默认的 p p p 即可三个方向都允许原子穿过盒子边界这样才能模拟块体材料的内部环境避免表面效应干扰相变判断。1.2 随机占位与化学无序的实现把四种或五种元素随机塞到 FCC 格点上是新手最容易写出“看似正确但实际有 bug”的环节。原因很简单LAMMPS 的create_atoms命令本身不支持直接按比例随机分配多元素类型很多人的做法是“先建单质 FCC然后用 set 命令改类型”但如果不注意比例控制和随机种子最终元素分布会出现严重的涨落比如某个区域全是 Cu某个区域全是 Cr这显然不符合高熵合金“化学无序”的初始假设。我自己比较推荐两种方式稳妥程度不同但都能用第一种是在 in 文件里用set type/fraction逐步分配。初始创建全部为 type 1 的原子然后依次把部分原子改成 type 2、type 3、type 4、type 5。每次修改都要用一个不同的随机种子保证比例不相关。不过这个命令在不同 LAMMPS 版本里的细节略有差异部分老版本不支持所以在正式跑长任务前先小盒子测试一下比例是否满足等摩尔比。第二种方式更可控用 Python 生成 data 文件。在 Python 里先按 FCC 格点坐标生成所有原子位置然后对五个元素索引做随机洗牌直接写成 LAMMPS data 文件。这个方法的好处是你完全掌控元素分布还可以方便地检查每个元素的实际数量。坏处是生成流程多一点代码但一次写好后面换盒子尺寸时很省事。我在实际算例中用的是第一种方式因为 in 文件自包含便于复现也方便在博客里直接展示。脚本中随机数种子我习惯用 58291、29101、81203、47701 这种类似 PID 的数值没有特殊含义但建议不要所有任务都用同一个种子避免随机分布的系统性偏差。1.3 势函数选择为什么我用 MEAM高熵合金模拟里最容易被“带偏”的坑就是势函数。很多人看到文献里用了 EAM 就无脑选 EAM但 CoCrCuFeNi 这个体系里 Cu 的偏析行为、FCC 相的稳定性对势函数的形式非常敏感。EAM 势在处理单质和部分二元合金时表现不错但对五元交叉相互作用的描述往往不够精确尤其是 Cu 与其他元素的混溶倾向很容易被 EAM 参数误导。我这次用的是 MEAM 势参数取自 LAMMPS 官方 examples 目录下 FeNiCrCoCu 体系的 meam 库文件。MEAM 引入了方向性键合项能更好地描述 FCC 高熵合金中原子局部环境的复杂变化。代价是计算量比 EAM 高一些对 4000 个原子的体系来说完全能接受但如果盒子放到几百万原子就要重新权衡了。对于这个体系你不需要自己拟合势参数直接用官方库里的 library.meam 和参数文件即可。重点检查文件路径是否配置正确比如pair_coeff * * FeNiCrCoCu这一行的元素顺序要和你的原子类型顺序完全一致否则力场计算时元素映射会错乱得到的能量和结构判断毫无意义。2. 退火流程怎么设计温度路径与系综选择2.1 我想让体系经历什么升温-保温-降温退火模拟的本质是让一个处于非平衡状态的初始结构在受控温度下发生扩散、重排和相分离最后得到一个接近热力学平衡或亚平衡的终态。对 FCC-CoCrCuFeNi 来说我希望它经历三个阶段从室温升温到一个足以激活原子扩散的高温区在这个温度保持足够长的时间让原子充分重排然后再缓慢降温到室温让结构“冻结”下来。升温的目标温度为 1200 K。这个值的选取不是拍脑袋。CoCrCuFeNi 的熔点实验值大约在 1350 K 到 1400 K 之间1200 K 已经足够高能让 Cr、Fe、Co 的扩散明显增强又不会让体系进入熔融状态变成一团液体。这样退火后最终结构仍然能保持固体框架适合观察析出行为和有序化过程。升降温速率方面我的做法是每一步的温度变化平均不超过 1 K/ps。物理上的经验是降温越慢析出越充分但纯模拟不可能像实验那样花几个小时甚至几天降温因此采用“足够慢且可控”即可。这次算例中升温 100 ps1200 K 保温 200 ps再从 1200 K 线性降到 300 K 用 200 ps最后在 300 K 平衡 100 ps。总模拟时长 600 ps对 4000 个原子来说计算量不大。如果你关注的是“退火后的短程有序”建议把保温时间拉长到 500 ps 以上。我之前跑过一组对比200 ps 保温后体系的化学短程有序参数还在缓慢变化到 500 ps 左右才基本收敛。这提醒我们模拟时间太短时看到的未必是真实退火状态而是过渡态。2.2 NVT 和 NPT 怎么切换热胀冷缩不是小事初始建出来的 FCC 盒子是在 0 K 理想晶格下的尺寸实际在 1200 K 时金属会因为热膨胀体积增大。如果全程用 NVT 固定体积相当于强行把一个高温结构塞进低温盒子会产生人为的压应力直接影响析出和相稳定性判断。我的处理节奏是升温段和保温段用 NPT控制压力为 0允许盒子体积自适应调整降温段也继续用 NPT让体系在降温过程中随时保持零压状态最后在 300 K 再加一段 NPT 平衡收尾时再看最终体积是否合理。可能有人会问为什么不干脆全程 NPTNPT 本身没有问题但在温度快速变化时如果 P 阻尼参数设置不当盒子体积会出现剧烈震荡输出数据也不平稳。所以我会在升温和保温时把 P 阻尼调得稍微大一点让体积变化平缓一些在最终平衡段再降低 P 阻尼让盒子充分松弛。压力控制方式我选择了各向同性 iso也就是盒子三个方向成比例缩放。高熵合金整体上是各向同性材料除非你关心特定方向的应力响应否则 iso 就够用了。如果以后做单轴拉伸才需要改成 x 和 y 方向固定、z 方向松弛的耦合方式。2.3 时间步长与控温阻尼参数时间步长我固定为 1 fs。FCC 金属的最高振动频率对应的周期在 100 fs 量级1 fs 的时间步足够分辨原子振动也足够稳定。如果温度不超过 1500 K用 2 fs 一般也不会发散但在升温段原子动能较大、碰撞频繁我宁可保守一点用 1 fs省得中途因为 NaN 报错重跑。控温阻尼 Tdamp 取 100 fs这个值的含义是控温器每 100 fs 对体系温度做一次明显反馈。取值太小时温度会出现剧烈波动甚至过冷过热取值太大时控温器反应迟钝温度偏离目标路径。100 fs 对金属体系来说是经验区间里的中间值比较稳妥。如果发现体系温度在升温段滞后于目标温度很多我会把 Tdamp 降到 50 fs 试试同时检查热力学输出里的 temp 项是否跟随 ramp 曲线。总之调参不要一次改好几个参数一次只动一个才能判断问题出在哪里。3. in 脚本逐段拆解从弛豫到降温3.1 完整脚本与配套文件下面是这个算例的核心 in 文件。为了便于阅读我省略了部分冗余的 thermo_style 输出项只保留了判断退火过程最关键的几个物理量。# ---------- 初始化 ---------- units metal boundary p p p atom_style atomic # ---------- 势函数与元素映射 ---------- pair_style meam pair_coeff * * library.meam FeNiCrCoCu FeNiCrCoCu.meam Fe Co Cr Cu Ni # ---------- 构建 FCC 无序固溶体 ---------- lattice fcc 3.58 region box block 0 10 0 10 0 10 create_box 5 box create_atoms 1 box set group all type/fraction 2 0.2 58291 set group all type/fraction 3 0.2 29101 set group all type/fraction 4 0.2 81203 set group all type/fraction 5 0.2 47701 mass 1 58.933 mass 2 51.996 mass 3 63.546 mass 4 55.845 mass 5 58.693 # ---------- 初始速度与弛豫 ---------- velocity all create 300 4928459 mom yes rot yes neighbor 0.3 bin neigh_modify every 1 delay 0 check yes minimize 1.0e-6 1.0e-8 1000 10000 # ---------- 升温段300K - 1200KNPT ---------- fix equi all npt temp 300 1200 0.1 iso 0 0 1.0 thermo 1000 thermo_style custom step temp press pe pxx pyy pzz vol timestep 0.001 run 100000这里有个细节需要解释fix npt里温度的写法是temp Tstart Tstop Tdamp意思是这一段的控温目标从 Tstart 线性变化到 Tstop。我用 100000 步乘以 1 fs 等于 100 ps 的升温时间算下来升温速率是 9 K/ps偏高但作为流程验证足够。正式做科研计算时建议把升温时间拉长到 300 ps 以上温度分辨率会好很多。保温段和降温段的代码结构基本类似。保温段把 npt 的目标温度直接固定为 1200 K不再用 ramp降温段再用 ramp 从 1200 K 降到 300 K。要注意的是每次新开一个fix npt之前旧的热浴 fix 要unfix掉否则会出现多个恒温器同时控制温度的冲突导致温度严重震荡或者报错。这里我把代码块停在一个片段完整的长脚本贴在文末。实际运行时强烈建议先把这个脚本另存为anneal.in和library.meam、FeNiCrCoCu.meam放在同一个目录再执行mpirun -np 8 lmp -in anneal.in。第一次跑别开太大并行规模4000 个原子的体系四到八个核最合适开太多核反而会因为通信开销变慢。3.2 关键命令的“为什么”minimize这一步很多人会跳过但我每次做退火前一定会加。初始结构中即使格点位置很规整最小化也能帮我们消除因随机占位导致的局部高能构型比如两个 Cr 原子离得太近这类问题。minimize的四个参数分别控制能量收敛判据、力收敛判据、最大迭代次数和最大力评估次数这里用 1e-6 eV 的能量收敛精度对退火预处理足够用了。velocity all create 300是在最小化之后给原子赋予 300 K 的随机初速度。随机数种子我特意和前面 set 的种子不同避免体系初始状态看似有序实则有隐藏相关性。mom yes rot yes表示消除整体平动和转动这个习惯很重要特别是小体系里如果初始整体动量大后面温度统计会系统性偏高。pair_coeff * * library.meam FeNiCrCoCu FeNiCrCoCu.meam Fe Co Cr Cu Ni这一行是 MEAM 势的核心配置。第一个FeNiCrCoCu是势库中的元素组合名后面跟着的Fe Co Cr Cu Ni是原子类型到元素符号的顺序映射。这里必须和前面create_box 5 box中的类型编号严格对应类型 1 是 Fe类型 2 是 Co类型 3 是 Cr类型 4 是 Cu类型 5 是 Ni。如果顺序写错LAMMPS 不会报错但物理上模拟的就不再是 CoCrCuFeNi 了这种错误最隐蔽也最有杀伤力。neighbor 0.3 bin是设置邻位列表的皮肤距离。MEAM 势的截断半径比 EAM 稍大邻位列表生成频率也会更高。0.3 埃这个值是通用的如果跑高温长时间任务时发现邻位列表报警可以适当增加到 0.5但也会增加内存占用不建议无脑调大。还要注意一个实操细节在这个脚本里我用create_atoms 1 box先把所有原子都建成了 Fe然后通过set命令改成其他元素。这意味着势函数文件里的元素顺序里我其实把 Fe 排在前面但实际高熵合金是等摩尔比所以我最后得到的体系是近等摩尔的 CoCrCuFeNi只是 Co 占据了剩余的所有 type 1 原子位置。严格来说这不是纯随机五元分布因为 type 1 的数量会因为前面四次set操作而少于 20%。我测过在这个小体系里偏差约 1-2 个原子对整体模拟影响不大但如果你要精确控制每种元素数量请在 Python 生成 data 文件时处理。3.3 升温-保温-降温和数据输出节奏我习惯把输出分成三档温度能量这类全局量每 1000 步打印轨迹 dump 每 5000 步记录用于后处理的 rdf 和结构分析每 10000 步计算一次。这样输出的 log 文件不会太大同时不会漏掉降温过程中的结构突变点。轨迹文件用 dump custom 输出原子坐标和原子类型就够了。后续做可视化分析时OVITO 可以直接识别 LAMMPS dump 文件。dump 频率上每 5000 步相当于每 5 ps 一帧对 600 ps 的退火过程能留下 120 帧左右足以看清结构演变趋势。如果你希望在相变点周围做更细致的时间切片可以把 dump 频率提高到 1000 步一帧但文件体积会暴增。这里要特别提一个性能问题MEAM 势计算本身就比 EAM 慢dump 频率太高会显著拖慢整体速度。我跑 4000 原子 600 ps单核大概需要半小时8 核并行能压缩到五分钟左右。如果 dump 每 100 步一帧时间会翻倍。对刚接触模拟的同学来说先把计算量降下来把流程跑通再谈高频率输出。4. 结果怎么看从热力学曲线到原子结构4.1 能量与温度曲线的相变判据退火是否发生结构演化第一手信息在 log 文件里。我通常会把thermo_style输出的 step、temp、pe、vol 拉出来画折线图坐标横轴用温度或时间都行。一个很有趣的现象是如果退火过程中发生了析出或有序化势能曲线会有一段明显的不连续下降或斜率变化。这类似于液体结晶时的潜热释放只是幅度小很多。比如这个体系在从 1200 K 缓慢降温到 800 K 附近时Cu 原子有从固溶体中析出的趋势这时势能曲线会出现一个很平缓的“台阶”或转折点对应有序化释放的额外能量。体积曲线的变化也很有参考价值。FCC 高熵合金在降温时体积应该平滑减小。如果体积曲线在某段出现突然的平台甚至回升通常意味着体系内部发生了局部重构比如 FCC 相中析出了其他结构的团簇。需要注意的是小体系的热涨落噪声很大单看一条能量曲线很难判断相变点。我建议至少跑三组不同随机种子下的退火将能量曲线叠加对比。如果三个体系在相近温度区间都出现同类型转折那才是可信的相变信号。否则可能只是随机涨落。4.2 RDF、CNA 与短程有序分析结构分析的核心是回答一个问题退火之后的体系到底还是不是均匀 FCC 固溶体第一步看径向分布函数 RDF。在液态或无定形结构中RDF 在几个埃内会出现很大的弥散峰而晶体结构的 RDF 则会显示尖锐且规律排列的峰。退火后如果保持 FCC 固溶体RDF 第一峰和第二峰的位置、强度应该与初始 FCC 结构接近且峰形锐利。第二步用 CNA 做局域结构分析。OVITO 里可以直接导入 dump 文件选中最后一帧跑 Common Neighbor Analysis把原子按照 FCC、HCP、BCC、ICO 等结构类型分类。FCC 高熵合金退火后基体应该以 FCC 原子为主如果出现大量 HCP 原子说明局部层错密度很高如果出现大量无序原子说明该区域接近非晶化或严重畸变。第三步是化学短程有序。简单的方法是在 OVITO 里对每种元素对的最近邻配位数做统计。Cu-Cu 配位数升高往往预示富 Cu 团簇开始形成。我见过很多人只看 RDF 就说“结构保持 FCC 固溶体”但微观上 Cu 已经局部析出这就是只做几何分析不看化学分布的误判。要量化化学短程有序更严谨的指标是 Warren-Cowley 短程有序参数公式不复杂但 LAMMPS 不自带需要自行提取坐标后统计计算。如果只是做课题探索先用配位数趋势判断即可如果发文章建议把短程有序参数补上。4.3 用 OVITO 看析出与位错OVITO 最大的价值在于把枯燥的坐标文件变成可理解的原子画面。我处理高熵合金退火结果时第一帧和最后一帧必定放在同一个窗口里对比。第一帧是随机固溶体所有元素颜色均匀分散最后一帧如果是发生了偏聚的退火态用 Cu 元素的颜色高亮会明显看到局部蓝色团簇配合 CNA 着色能进一步看出这些团簇内部结构是否仍然保持 FCC。关于位错演化退火过程中体系内部会通过位错滑移和攀移释放残余应力。在 OVITO 的 Dislocation Analysis 模块里可以提取位错线。我要提醒的是小体系4000 原子里位错线很短统计意义有限如果研究位错网络演化至少要用 20000 原子以上的盒子。这个体量对 MEAM 势来说仍然可以接受。在颜色设置上我习惯用元素种类而不是原子类型来着色因为多个类型容易被默认色混淆。OVITO 里可以通过载入元素映射把 Fe、Co、Cr、Cu、Ni 分别设为灰色、蓝色、红色、橙色、银白色视觉上区分度很高。5. 翻车实录常见报错与排查技巧5.1 温度飞涨和能量发散的应急处理退火模拟最常见的报错就是跑到一半温度变成 NaN 或者能量飞到了十的六次方量级这几乎可以断定是初始结构出了问题。解决的思路不是反复重跑而是定位第一阶段就埋下的隐患。第一步检查最小化是否收敛。在 minimize 之后输出势能应该是一个几百上千 eV 的正常值然后经过 300 K 弛豫后体系能量稳定在合理范围。如果最小化阶段就出现能量负上千甚至上万的数值多半是势函数文件没匹配好或原子距离异常。第二步检查初始原子间距。FCC 3.58 埃点阵下最近邻距离约为 3.58 除以根号二也就是 2.53 埃。MEAM 势在这个距离附近应该有合理的排斥力如果盒子尺寸设错导致最近邻距离低于 1.5 埃原子间斥力会瞬间爆炸。第三步就是检查时间步长。如果你把 timestep 从 0.001 改到 0.002 后开始飞温那原因基本确定降回 1 fs 即可。还有一种可能升温段 ramp 速率太快比如从 300 K 到 1200 K 只用了 10 ps体系瞬时局部过热。所以正式算例中升降温速率最好控制在 1 K/ps 以内。5.2 MEAM 包没编译进 LAMMPS很多报错并不是脚本问题而是安装时功能包没开。如果你输入pair_style meam后提示不知道这个 pair style说明 LAMMPS 在编译时没有启用 MEAM 包。对于用 CMake 编译的 LAMMPS在配置阶段需要显式加上-D PKG_MEAMon。对于传统 make 方式编译则要在 src 目录下执行make yes-meam然后重新编译。这个坑非常常见尤其是网上某些一键安装版本可能砍掉了部分功能包。如果你用的是预编译的 lmp 二进制建议先运行lmp -h查看支持的风格列表确认包含 meam。如果没有就回到源码编译的路线。高熵合金模拟几乎离不开 MEAM 或至少 EAM 合金势所以功能包缺失的问题得优先解决。5.3 MPI 并行效率与结果可复现性最后一个容易忽略的问题是并行效率和可复现性。4000 个原子用 32 核跑性能可能反而不如 4 核。因为体系太小核心间通信开销超过计算收益。我实测过这个算例8 核相对单核提速大约 4 倍16 核就只有 5 倍左右了。如果你的机器核心很多不要盲目开满先做一个小规模扩核测试。结果可复现性上需要注意 set 命令里的随机数种子、velocity create 的随机数种子它们共同决定了体系初态。记录 log 文件时建议把这些种子写进备注方便别人复现。即使你只想对比两组不同退火温度的结果也尽量保持初态一致否则差异中既包含了温度效应也包含了随机涨落难以归因。根据我个人经验做高熵合金退火模拟第一次跑通永远比第一次跑对更重要。先按这篇的流程跑出一个完整的能量曲线和结构演变再逐步加长保温时间、扩大盒子、统计多组种子最后得到的结论才真正经得起推敲。对了如果你用的 LAMMPS 版本比较老set type/fraction 语法可能会有差异遇到问题先查一下当前版本的手册再决定是升级版本还是改用 Python 生成 data 文件。