结合案例图分析，我把这个项目的完整过程、参数设置、结果解读和踩坑经验整理如下。
客户要的是纯金属的熔化温度和熔化焓，用于铸造工艺模拟的材料数据库。他提供的是 hcp 结构的纯金属（图上给出 T_m=933 K、ΔH_m=0.111 eV/atom），要求用分子动力学给出与实验可比的数据。
我采用的是两相共存法（two-phase coexistence method）：先构建一半固态、一半液态的界面体系，然后在 NPT 系综下缓慢升温，观察固-液界面移动。软件用 LAMMPS，势函数用 EAM（嵌入原子法）势，适合金属体系。
参数设置：体系约 8000 原子，x/y 方向周期边界，z 方向有固液界面；NPT 系综，Nose-Hoover 热浴和压力控制；温度从 700 K 升到 1100 K，步长 10 K，每个温度平衡 200 ps；共做了 8 次独立模拟，每次用不同随机速度种子，以估计统计误差。
第一张案例图是单个模拟的势能随温度变化曲线。蓝色线是固态 hcp 相，红色线是液态。两条线的交点对应熔点 T_m=933 K，与实验值（如果体系是 Mg，实验 T_m=923 K；偏差约 1%）吻合得很好。从交点处的能量跳跃可以读出熔化焓 ΔH_m=0.111 eV/atom。
第二张案例图是 8 次独立 MD 运行得到的 ΔH_m 统计分布。每个蓝点是一次独立模拟的结果，误差棒代表该次模拟的块平均标准差；红色虚线是平均值 0.111 eV/atom（≈10.7 kJ/mol，与图上实验值一致），绿色虚线标注实验值 10.7 kJ/mol。8 次结果分布在 0.105–0.120 eV/atom 之间，平均值与实验值几乎重合，说明两相共存法对这个体系是可靠的。
第一，体系尺寸要足够大。我第一次用 2000 原子试算，固液界面波动很大，T_m 偏差超过 50 K；提到 8000 原子后，结果稳定在 933±5 K。
第二，升温速率不能快。每 10 K 平衡 200 ps 是经验下限，如果缩短到 50 ps，界面移动还没达到稳态，会得到虚假的过热现象。
第三，势函数决定一切。EAM 势对金属熔化描述通常不错，但不同文献参数给出的 T_m 能差几十 K。我对比了三个版本的 EAM 势，选了与实验最接近的那个。
下次再做金属熔化模拟，我会按这个 checklist：1) 选合适的 EAM 或 MEAM 势，先用块体性质（晶格常数、弹性常数、空位能）做验证；2) 构建固液两相体系时，让界面尽量平整且取向明确；3) 每个温度点充分平衡并做 3–8 次独立运行；4) 同时输出 T_m 和 ΔH_m，并和实验值对比。这样做出来的数据可以直接进材料数据库。
计算工程师团队：500+ ｜ 成功计算案例：5W+ ｜ 1W+ 助力顶刊发布经验 需要模拟代算服务可以联系我们，可免费评估项目的计算费用和免费出计算方案。
