金属原子间键能计算的方法选择与精度控制——从DFT到EAM势

金属原子间键能计算这个问题，乍看是教科书上的基础内容——两个金属原子之间的键能，DFT跑个二聚体就能出来。但实际问题里的”键能”很少是孤立的两个原子，而是合金体系中某种特定原子对在复杂化学环境下的有效键合强度。后者的计算就远不是跑个二聚体那么简单。

三年前做高熵合金力学性能预测时遇到的真正困难，不是不会算键能，而是算出来的键能怎么转化为可解释的力学性能指标。

直接DFT方法：晶体轨道哈密顿布居分析

对于有序结构的合金体系，晶体轨道哈密顿布居（COHP）分析是金属原子间键能计算最直接的手段。原理不复杂：把DFT算出来的电子波函数投影到原子轨道基组上，再拆解哈密顿矩阵元中不同原子对之间的成键/反键贡献。负的COHP积分代表成键态占据，直观反映键合强度。

实际操作上，Lobster程序是目前走这条路最成熟的选择。它可以直接读VASP的WAVECAR输出，做COHP投影。在做高熵合金的CoCrFeMnNi体系时，COHP分析给出了不同原子对之间键合强度的排序：Cr-Mn对的反键态占据明显高于Co-Fe对，这个信息和实验中观察到的Cr-Mn偏聚倾向对应得很好。

但COHP有个明显的局限：它给的是能量量纲的积分值，不同体系之间的绝对值不能直接比较。比如体心立方Fe中Fe-Fe键的COHP积分和面心立方Ni中Ni-Ni键的COHP积分，数值差得再多也不能直接说”Fe-Fe键比Ni-Ni键弱”。只有同一体系、同一基组下的相对比较才有意义。这个限制在一次跨体系的对比项目中暴露得特别明显——团队花了将近一周把COHP结果重新归一化，才把不同结构的键能放到同一个尺度上讨论。

EAM势函数路径：从DFT数据到经验势

当体系大到几百个原子、需要跑分子动力学来算弹性常数或位错行为时，DFT直接算键能就不现实了。这个时候要走嵌入原子势（EAM）的路线：用DFT算少量代表性构型的能量和力作为训练集，再拟合EAM势的参数，最后用拟合好的势函数在MD尺度上计算各类力学性能。

这个路线的精度瓶颈在拟合这一步。EAM势的核心假设是：原子i的能量等于嵌入能+对势之和，其中嵌入能只依赖于原子i所在位置的背景电子密度。这个”局域电子密度”的假设在简单金属（Cu、Al、Ni）里非常有效，但在过渡金属合金中——特别是含d电子方向性成键的体系——嵌入函数的形式不够灵活，往往低估了特定原子对之间的键合各向异性。

当时在拟合CoCrFeNi四元合金的EAM势时，团队用DFT的AFLOW数据库里已有的二元和三元数据做了初始训练，但在力学测试（比如单轴拉伸模拟）中，拟合出来的势函数预测的屈服强度比DFT基准值低了约15%。追溯原因，发现训练集里缺少足够的高应变构型——EAM势在平衡晶格参数附近拟合得很好，但在远离平衡的变形状态下就开始偏离。

这个问题的补救办法是在训练集中显式加入一系列不同应变状态的DFT构型，但每加一组构型，DFT计算量就翻倍。最后的妥协方案是：EAM势负责力学性能的整体趋势判断，关键结论再用DFT做有限构型的基准验证。严格来说这不是一个最优方案，但在当时的算力和时间约束下，这是唯一能把任务推下去的路径。

键能与力学性能的关联：从原子尺度到宏观

金属原子间键能计算最终要回答的问题，通常不是”这个键有多强”，而是”这个键的强度怎么影响材料的力学性能”。

对于高熵合金，晶格畸变导致不同原子对的键长分布变宽，键能分布也随之变宽。DFT计算可以给出每种原子对在不同局域环境下的键能分布，但这个分布怎么映射到宏观的弹性模量或屈服强度——中间需要一个连贯的多尺度桥梁。

团队的做法是：DFT算键能和局域弹性常数→EAM势做大尺度力学测试→把EAM结果中的局域应力和DFT算出来的键能分布做线性回归。在CoCrFeMnNi体系中，这个回归给出了R²≈0.85的相关性——说明键能分布确实能解释大部分局域力学响应的变化，但剩下15%的方差来自位错核心结构这类键能解释不了的贡献。这份工作后来放到了科研学术网首页的计算案例专栏，作为多尺度键能分析的参考。