1. 项目概述当材料模拟遇上“基因进化”在计算材料科学的世界里我们这些做模拟的每天都在和两个“冤家”打交道一个是精度极高但贵得离谱的量子力学方法比如密度泛函理论DFT另一个是速度快如闪电但精度时常“掉链子”的半经验势函数比如嵌入原子法EAM。前者能告诉你电子在原子间如何“跳舞”但算个几百个原子的体系就得等上几天几夜后者能让你模拟上百万个原子跑上纳秒尺度但预测复杂缺陷的性质时结果可能和实验差之千里。这中间的鸿沟就是材料设计中的“精度-效率”困局。近年来机器学习原子间势MLIP的出现像一座精心设计的桥梁试图连接这两端。它的核心思想很直观用大量高精度的DFT计算结果作为“标准答案”训练一个机器学习模型让它学会根据原子的位置和种类预测整个体系的能量、原子受力和应力。这样在后续的分子动力学模拟中我们就能用这个训练好的、计算成本极低的“代理模型”来代替昂贵的DFT实现既准又快的模拟。在众多MLIP模型中矩张量势MTP因其数学上的优雅和计算上的高效受到了广泛关注。它用一组旋转协变的“矩张量”来描述一个原子周围邻居的几何环境形象点说就是给每个原子的“邻里关系”拍了一张多维度的、包含距离和方向信息的“快照”。最终的总能量就是这些张量经过一系列特定规则收缩组合后的线性加权和。然而问题就出在这个“组合”过程上。传统的MTP在如何高效地组合收缩这些高维张量上采用的是一种启发式策略并非最优解。这就好比你要组装一个复杂的乐高模型说明书只给了你大概的步骤但没告诉你哪些通用零件可以提前批量拼好导致你在拼装过程中反复进行重复劳动效率低下。随着模型复杂度基组水平提升这种低效的收缩操作会消耗大量计算资源成为性能瓶颈。本文要分享的就是我们针对这个瓶颈的一次“外科手术式”优化。我们引入了一个在优化领域久经沙场的“老兵”——遗传算法GA来为MTP的矩张量收缩过程寻找最优的“装配流水线”。通过这套基于遗传算法的优化方案我们成功地将MTP模型在评估时的计算速度提升了近一个数量级并且在描述Ni-Al合金中各种棘手缺陷时精度还有所提升。如果你正在为大规模材料模拟的计算成本发愁或者对如何将智能优化算法与物理模型结合感兴趣那么接下来的内容或许能给你带来一些新的思路。2. 核心原理拆解MTP的“张量乐高”与遗传算法优化要理解我们的优化工作首先得弄明白MTP是怎么工作的以及传统方法的瓶颈在哪。这就像修车你得先知道发动机的构造才能找到提升动力的方法。2.1 矩张量势MTP的核心从原子环境到标量能量MTP的核心在于如何数学化地描述一个原子i的局部环境。它定义了一个称为矩张量描述符的量M_μ,ν(n_i) Σ_j [ f_μ,ν(|r_ij|, z_i, z_j) * f_c(|r_ij|) * (r_ij ⊗ ... ⊗ r_ij) ]这个公式看起来复杂但我们可以分部分理解Σ_j对中心原子i的所有邻居原子j求和。f_μ,ν(|r_ij|, z_i, z_j)这是一个径向基函数比如切比雪夫多项式它负责处理原子间距离|r_ij|和原子种类z的信息。μ是这个径向函数的索引。f_c(|r_ij|)截断函数保证在某个截断半径外相互作用平滑衰减至零这是所有经验势的常见做法以保证计算的局域性。(r_ij ⊗ ... ⊗ r_ij)这是最关键的张量部分。r_ij是从原子i指向原子j的向量。符号⊗表示克罗内克积一种张量积。ν表示这个向量被自乘了多少次也就是张量的阶数rank。物理图像你可以把ν0的矩张量想象成这个原子邻居的总“质量”标量ν1的矩张量一个向量描述了这个邻居团的“质心”位置ν2的矩张量一个矩阵则类似于“转动惯量”张量描述了质量分布的形状。更高阶的张量则包含了更精细的几何不对称性信息。MTP通过组合不同μ径向和ν角度的矩张量就能以多项式展开的形式逼近任意光滑的原子间相互作用势。单个矩张量M_μ,ν本身并不是旋转不变的旋转坐标系它的分量会变因此不能直接作为能量的贡献。MTP的巧妙之处在于它通过将多个不同阶的矩张量按特定规则进行缩并即对某些共同的索引进行求和最终生成一个旋转不变的标量B_α。这个缩并规则由一个对称索引矩阵α来定义。整个模型的原子能量U_i就是许多这样的标量基函数B_θ的线性组合U_i Σ_θ ξ_θ B_θ其中ξ_θ是通过拟合DFT数据得到的线性系数。2.2 传统收缩的瓶颈与“计算树”思想现在到了关键部分如何高效计算成千上万个这样的标量基函数B_α每个B_α都是多个矩张量缩并的结果。直接暴力计算所有可能的缩并索引计算量会随着缩并索引的数量指数增长3^mm为缩并索引数。Shapeev在提出MTP时想到了一个聪明的办法将每个标量函数B_α的缩并过程表示为一棵计算树。这棵树的叶子节点是最基础的矩张量M中间节点是缩并过程中产生的中间张量B_η根节点就是最终的标量B_α。计算时从叶子节点矩张量开始按照树的结构两两结合生成父节点中间张量最终得到根节点。核心洞察不同的标量函数B_α和B_β它们的计算树可能共享相同的中间子树即相同的中间张量B_η。例如B_α和B_β可能都需要先计算(M1·M2)这个中间结果。传统方法为每个标量函数独立构建计算树会重复计算这些共享的中间张量。因此优化的目标就变成了为整个基组所有需要计算的标量函数集合寻找一组全局最优的计算树使得需要存储和计算的独立中间张量分量总数、以及执行缩并运算的“乘法-加法”规则总数最小化。这本质上是一个复杂的组合优化问题。传统MTP采用了一种基于“最小化维度之和”和“最大化张量复用”的启发式贪婪算法来寻找这组树。但这种方法在面临大量选择时容易陷入局部最优并且当多个分解路径的“代价”相同时无法做出最优决策。2.3 遗传算法为“计算树”注入进化力量面对这个搜索空间巨大例如对于level-28的基组可能性超过10^2195种的优化问题我们引入了遗传算法GA。GA是一种受生物进化启发的全局优化算法特别适合解决这类组合爆炸问题。我们将整个优化过程设计如下基因编码每个标量函数B_α都有多个可能的计算树分解方式构成一个“树集合”{T_i^α}。我们为每个B_α分配一个“基因”这个基因的值k表示我们选择其树集合中的第k棵树。那么对于包含n个标量函数的整个基组一个完整的分解方案个体就可以用一个基因序列[k1, k2, ..., kn]来表示。种群初始化随机生成一定数量如40个的个体构成初始种群。为了加速收敛我们并非完全随机而是给每棵树一个初始“评分”例如选择维度之和最小的树、能复用已有中间张量的树会得分更高优先选择高得分的树来构建初始个体。适应度函数这是进化的“指挥棒”。我们定义了一个成本函数F(R) m 3t其中R代表一个分解方案m是独立张量分量的总数t是计算规则的总数。给t乘以权重3是因为在实际计算中执行这些规则乘加运算的开销通常比存储张量更大。适应度函数f(R)被设计为与F(R)负相关且采用类似费米-狄拉克分布的形式确保种群多样性避免过早收敛于次优解。选择、交叉与变异选择根据个体的适应度采用“轮盘赌”等方式选择优秀的个体进入下一代。交叉随机选取两个父代个体交换它们的一部分基因即交换一部分标量函数所选择的计算树产生新的子代个体。这模拟了基因重组。变异以一定概率随机改变某个个体中某个基因的值即为某个标量函数随机更换一棵计算树。这为种群引入了新的可能性避免陷入局部最优。迭代与收敛重复进行选择、交叉、变异并计算新一代的适应度如此迭代数万乃至数十万代。直到达到预设的最大迭代次数或最优解在连续多代中不再改善。冲突解决与后优化在进化过程中一个关键挑战是基因间的“冲突”。例如标量函数A和B可能都需要用到同一个中间张量X但A选择的计算树要求X由(M1·M2)生成而B选择的树却要求X由(M3·M4)生成这就产生了矛盾。我们的算法在生成计算规则时会强制统一这些冲突确保一致性。在GA搜索结束后我们还会执行一个快速的贪婪算法进行局部微调从低阶到高阶张量逐一尝试替换为更优的计算树确保得到最终的最优方案。通过这套流程遗传算法能够在浩如烟海的分解方案中高效地寻找到那个能最大程度复用中间结果、最小化计算量的全局优解。图2清晰地对比了传统方案与我们优化后方案的工作流差异传统方案是“各自为政”的独立计算而我们的方案则是经过全局规划后的“协同流水线”。3. 基组设计策略在精度与效率间寻找新平衡在优化了“怎么算”之后下一个问题是“算什么”即如何选择构成势函数的标量基函数B_α的集合基组。基组决定了模型的表达能力和复杂度。传统MTP使用一个称为“水平level”的参数lev(M_μ,ν) 4μ ν 2来筛选基函数只保留总水平低于某个阈值d的基函数。然而我们通过分析发现计算成本主要来自两方面计算矩张量分量本身对于一个ν阶张量其独立分量数以ν^2规模增长。计算每个分量都需要遍历中心原子的所有邻居而邻居数量与截断半径rc的立方成正比。这部分计算无法通过优化缩并过程来减少。执行缩并规则这正是我们上一节用遗传算法优化的部分。因此传统的水平定义4μ ν 2对μ径向复杂度的惩罚因子4过高导致为了达到高多项式精度不得不引入高阶ν大的矩张量而这会急剧增加第一部分的计算量。我们的策略是重新定义水平lev(M_μ,ν) 2μ ν 1。我们降低了对径向复杂度μ的惩罚从4降到2但严格限制了张量的阶数ν在我们的工作中通常ν ≤ 3。这意味着我们使用更多低阶ν小的矩张量但通过让更多这样的张量参与缩并即考虑更高阶的多体相互作用如五体、七体来提升模型的表达能力。这样做的优势非常明显低阶张量的分量数少计算成本低。虽然增加参与缩并的张量数量会增加第二部分缩并的复杂度但这部分恰恰是我们用遗传算法重点优化的部分。实测表明这种“低阶张量高阶相互作用优化缩并”的策略能在不损失精度甚至有所提升的前提下带来数量级的计算加速。例如我们定义的“2653”基组包含2653个线性参数考虑至五体相互作用ν_max4其精度与传统的level-28基组相当但计算速度却快了近11倍。4. 实战构建用于Ni-Al合金缺陷模拟的高效MTP势理论再好也需要实战检验。我们选择Ni-Al合金体系作为试金石因为其内部的γ/γ’相结构以及各类缺陷空位、反位、层错等对力学性能至关重要也是检验势函数准确性的经典难题。4.1 训练集构建主动学习确保“见过世面”一个MLIP的可靠性极大程度上取决于其训练集是否足够全面和具有代表性。我们采用了一种两阶段主动学习策略来构建高质量训练集第一阶段基于贝叶斯误差的在线学习。我们使用VASP软件进行从头算分子动力学AIMD模拟。在模拟过程中集成一个基于高斯过程回归GAP的在线学习模块。这个模块会实时评估当前原子构型对于当前GAP模型是否属于“未知领域”即预测不确定性高。一旦发现这类构型就暂停MD调用DFT计算其精确的能量、力和应力并将其加入训练集然后更新GAP模型。这个过程能自动探索相空间捕获动力学过程中出现的稀有但重要的构型。第二阶段基于D-最优性准则的主动采样。用第一阶段得到的训练集初步拟合一个MTP势。然后我们用这个MTP势进行更长时间的MD模拟200 ps生成大量新的候选构型。接着利用Shapeev提出的广义D-最优性准则从这些候选构型中筛选出那些能最大程度提升模型泛化能力即最大程度降低参数不确定性的构型进行DFT计算并加入训练集。随后重新拟合MTP。如此循环直到没有新的构型能显著扩展模型的认知边界。最终我们得到了一个包含8450个构型的训练集涵盖了fcc-Ni、L12-Ni3Al、相界面、各类堆垛层错GSF以及它们与空位、空位团、反位缺陷的组合构型。同时我们构建了一个独立的包含815个构型的验证集用于评估模型的泛化能力。通过主成分分析PCA可以看到验证集很好地覆盖了训练集所占据的构型空间。4.2 势函数拟合两步法求精在确定了优化的2653基组后我们使用完整的训练集进行最终拟合。拟合过程分为两步以确保找到全局最优的线性系数ξ_θ初始化线性系数首先从完整的2653基组中提取一个最小的子集例如只包含最重要的基函数用这个子集去拟合训练数据。拟合得到的势函数其径向基函数是确定的。然后我们固定这些径向基函数在完整的2653基组上以子集拟合的结果作为初始点进行一次线性最小二乘优化得到一组较好的线性系数初值。这一步的目的是避免直接在高维参数空间中随机初始化可能导致的收敛缓慢或陷入局部极小。非线性精修以上一步得到的系数为起点采用L-BFGS算法进行最多5000次迭代的非线性优化进一步调整所有参数包括线性系数和径向基函数的参数。L-BFGS迭代完成后再执行一次最终的线性最小二乘拟合对线性系数进行微调。在拟合过程中能量、力和应力的权重分别设置为1、0.01和0.005。这个权重的设置很有讲究力的数据量远大于能量每个构型有3N个力分量适当降低其权重可以防止力误差主导优化过程而应力数据则更少权重也最低。4.3 效率与精度评估优化带来的显著提升我们对优化前后的不同基组进行了全面的基准测试测试在一个2048原子的Ni3Al超胞上进行10万步MD模拟。结果令人振奋参见表I对于简单基组如level-18优化将中间张量数量、独立分量数和计算规则数减少了约8%-19%。但由于此时计算成本主要由计算矩张量分量本身主导MD模拟速度仅提升了约4%。对于复杂基组如level-28优化的效果极为显著。中间张量数量减少18.2%独立分量数减少41.3%计算规则数减少43.9%。这直接转化为71%的MD模拟加速时间从4418秒减少到2580秒。这是因为在复杂基组下缩并操作的成本占据了主导我们的优化直击要害。我们的新基组2653 vs 传统基组level-28这是最关键的对比。两者具有相近的线性参数数量~2500个但我们的2653基组由于采用了低阶张量设计其MD模拟速度是level-28基组的11倍。更令人惊喜的是在验证集上2653基组的能量、力、应力的预测误差均低于level-28基组。这证明了我们“降低张量阶数、增加多体阶数、优化收缩”策略的成功它在降低计算复杂度的同时反而获得了更好的拟合能力。表II展示了训练阶段的成本。在L-BFGS迭代过程中优化带来的加速比与MD模拟中观察到的趋势基本一致。对于线性拟合步骤由于计算成本主要来自求解大型线性方程组优化带来的加速不明显但对于小基组level-18仍有提升。5. 性能验证Ni-Al合金缺陷模拟的实战表现一个势函数好不好最终要看它解决实际科学问题的能力。我们将优化后的MTP势基于2653基组应用于Ni-Al合金的一系列缺陷性质预测并与DFT结果以及经典的EAM势进行对比。5.1 基础性质晶格、弹性与声子首先我们检验了模型对完美晶体基本性质的再现能力晶格常数与能量-体积曲线MTP势精确地复现了Ni和Ni3Al的平衡晶格常数以及它们在压缩和膨胀时的能量变化曲线E-V曲线与DFT结果高度吻合。弹性常数对于立方晶体关键的C11, C12, C44等弹性常数MTP的预测误差在几个GPa以内远优于EAM势这对于模拟材料的力学响应至关重要。声子色散谱我们通过冻结声子法计算了Ni和Ni3Al的声子谱。MTP势不仅准确预测了声子支的整体走势甚至复现了DFT计算中一些细微的“扭结”kinks这表明它很好地捕捉到了原子间相互作用的细节。而EAM势通常在高频部分偏差较大。5.2 点缺陷空位与反位缺陷点缺陷是合金中最基本的缺陷类型。空位形成能我们计算了Ni和Ni3Al中单个空位的形成能。MTP的结果与DFT值几乎一致而EAM势通常系统性偏低或偏高。空位团结合能这是更严峻的考验。我们计算了从二空位到小空位团如四面体空位团的结合能。MTP势准确地描述了空位之间的相互作用随距离和构型的变化趋势成功预测了DFT计算中发现的某些特定构型的相对稳定性。这对于研究辐照或淬火过程中空位的聚集行为至关重要。反位缺陷形成能在Ni3Al中一个Al原子占据Ni的位置Ni位反位或一个Ni原子占据Al的位置Al位反位其形成能不同。MTP势精确区分了这两种能量与DFT相符而EAM势往往难以准确描述这种化学有序性。5.3 面缺陷堆垛层错与相界面面缺陷直接影响材料的塑性变形机制。广义堆垛层错能GSFE曲线我们计算了Ni和Ni3Al中不同滑移面如{111}面上的GSFE曲线。这条曲线决定了位错分解的宽度和层错能是计算材料强度的关键输入。MTP势给出的GSFE曲线与DFT计算的重合度极高特别是内禀层错能ISF和非内禀层错能USF的数值。相比之下不同EAM势给出的结果分散且常与DFT有较大偏差。γ/γ’相界面能Ni基高温合金的强度源于γ相基体中弥散分布的γ’Ni3Al强化相。两相之间的界面能是控制组织形貌的关键参数。MTP势计算得到的界面能与DFT参考值吻合良好。5.4 动力学性质空位扩散最后我们通过分子动力学模拟研究了空位在Ni中的扩散行为。我们统计了空位的均方位移随时间的变化计算了扩散系数。MTP势模拟得到的扩散系数与激活能与实验值及基于DFT的过渡态理论计算结果处于合理范围内。这表明优化后的MTP势不仅能描述静态性质还能可靠地用于研究有限温度下的动力学过程。6. 实操心得与避坑指南基于这个项目的完整开发和应用流程我总结了一些对于从事MLIP开发或应用的研究者可能极具价值的实操经验。6.1 遗传算法参数调优耐心与平衡种群大小与迭代次数种群大小不宜过小如小于20否则多样性不足容易早熟收敛也不宜过大如大于100否则每次迭代的计算开销剧增。我们发现在40-60之间是个不错的起点。迭代次数需要足够多对于复杂基组如461310万次迭代是必要的。可以观察最优个体适应度随代数的变化曲线当曲线进入漫长的平台期时即可考虑停止。适应度函数设计成本函数F(R) m w * t中的权重w是关键。我们通过 profiling 发现在目标硬件上执行一条“乘法-加法”规则对应t的成本大约是访问一个张量分量对应m成本的3倍左右因此设置了w3。这个权重需要根据具体的代码实现和硬件架构进行微调。如果代码对内存访问做了极致优化w值可能降低如果计算规则部分存在瓶颈w值可能需要提高。交叉与变异概率过高的交叉概率会导致优良基因模式被破坏过高变异概率则会使搜索退化为随机游走。我们采用自适应策略初期使用较高的交叉概率如0.8和中等变异概率如0.1以快速探索空间后期降低交叉概率如0.5略微提高变异概率如0.15以加强局部搜索和跳出局部最优的能力。6.2 训练集构建质量重于数量主动学习是必须的而非可选的单纯从完美晶体、弹性变形等常规构型采样得到的训练集其势函数在模拟缺陷、表面、熔化等“非典型”过程时极易外推失败extrapolation。在线学习on-the-fly和基于D-最优性准则的采样是确保模型“见过世面”的最有效手段。不要吝啬在主动学习循环中花费的DFT计算资源这能从根本上避免后续大规模模拟得到错误结果。警惕“数据分布不均”在主动学习过程中可能会在某个区域如熔融状态采集到过多构型。需要在构建最终训练集时进行适当的下采样或者为不同区域的构型设置不同的权重防止模型过度拟合某一类构型而忽略其他。必须包含应力数据对于固体材料尤其是要计算弹性常数或模拟力学变形训练集中必须包含DFT计算的应力张量信息。这对正确描述体积相关的能量变化至关重要。我们的经验是应力数据的权重可以设得比力数据低一个数量级。6.3 基组选择不是越复杂越好从简开始逐步增加不要一开始就使用参数最多的基组。从一个中等水平的基组如我们的“2653”开始训练和测试。如果验证误差已经满足要求且进一步增加基组复杂度如到“4613”带来的精度提升微乎其微但计算成本显著增加那么就应该停止。模型的“奥卡姆剃刀”原则同样适用。关注截断半径截断半径rc需要仔细选择。太小会丢失重要的长程相互作用对于某些合金或半导体很重要太大会急剧增加计算量邻居数 ~rc^3。通常rc需要至少包含第二近邻或第三近邻原子。可以通过绘制径向分布函数RDF或测试不同rc下对结合能、弹性常数的影响来确定。验证“能量-力-应力”的一致性一个健康的势函数其预测的能量、力和应力应该是自洽的。一个简单的检查方法是对一个平衡构型施加一个微小应变用势函数计算能量变化和应力同时用有限差分法通过位移原子计算力再推导应力来验证。两者应该基本一致。6.4 应用模拟中的注意事项外推预警所有MLIP都内置了外推预警机制通常基于训练集构型在描述符空间的距离。在运行MD时务必监控这个“外推指标”extrapolation grade。一旦有原子构型的指标超过阈值比如我们设为5模拟就不可信了。此时应该停止模拟将该构型加入训练集重新拟合势函数。温度与相稳定性用MLIP进行高温MD模拟前最好先测试一下它在高温下对相稳定性的预测。例如对于Ni-Al合金可以跑一个升温模拟看γ和γ’相是否能在正确的温度范围内保持稳定以及熔点是否合理。我们的MTP势在预测Ni和Ni3Al熔点方面表现良好这是许多EAM势的短板。缺陷模拟的有限尺寸效应在计算缺陷形成能时超胞尺寸必须足够大以消除周期性镜像缺陷之间的相互作用。一个实用的方法是逐步增大超胞尺寸计算缺陷形成能直到其变化小于目标精度如0.01 eV。对于空位、层错等通常需要包含数百个原子的超胞。通过将遗传算法这一全局优化工具深度融入MTP势的构建流程我们不仅获得了一个在计算速度上提升近一个数量级的高效势函数更重要的是建立了一套可复用的、系统化的势函数开发和优化方法论。这项工作表明在计算材料科学中对底层数学模型和算法进行“精雕细琢”的优化与采集更多数据、使用更复杂的网络结构同样重要甚至能带来事半功倍的效果。对于Ni-Al合金乃至更复杂的多元合金体系这种高效的模拟工具将为在原子尺度上深入理解缺陷演化、相变动力学和力学性能打开一扇新的大门。
遗传算法优化矩张量势：提升材料模拟效率与精度
发布时间：2026/5/24 16:19:49
1. 项目概述当材料模拟遇上“基因进化”在计算材料科学的世界里我们这些做模拟的每天都在和两个“冤家”打交道一个是精度极高但贵得离谱的量子力学方法比如密度泛函理论DFT另一个是速度快如闪电但精度时常“掉链子”的半经验势函数比如嵌入原子法EAM。前者能告诉你电子在原子间如何“跳舞”但算个几百个原子的体系就得等上几天几夜后者能让你模拟上百万个原子跑上纳秒尺度但预测复杂缺陷的性质时结果可能和实验差之千里。这中间的鸿沟就是材料设计中的“精度-效率”困局。近年来机器学习原子间势MLIP的出现像一座精心设计的桥梁试图连接这两端。它的核心思想很直观用大量高精度的DFT计算结果作为“标准答案”训练一个机器学习模型让它学会根据原子的位置和种类预测整个体系的能量、原子受力和应力。这样在后续的分子动力学模拟中我们就能用这个训练好的、计算成本极低的“代理模型”来代替昂贵的DFT实现既准又快的模拟。在众多MLIP模型中矩张量势MTP因其数学上的优雅和计算上的高效受到了广泛关注。它用一组旋转协变的“矩张量”来描述一个原子周围邻居的几何环境形象点说就是给每个原子的“邻里关系”拍了一张多维度的、包含距离和方向信息的“快照”。最终的总能量就是这些张量经过一系列特定规则收缩组合后的线性加权和。然而问题就出在这个“组合”过程上。传统的MTP在如何高效地组合收缩这些高维张量上采用的是一种启发式策略并非最优解。这就好比你要组装一个复杂的乐高模型说明书只给了你大概的步骤但没告诉你哪些通用零件可以提前批量拼好导致你在拼装过程中反复进行重复劳动效率低下。随着模型复杂度基组水平提升这种低效的收缩操作会消耗大量计算资源成为性能瓶颈。本文要分享的就是我们针对这个瓶颈的一次“外科手术式”优化。我们引入了一个在优化领域久经沙场的“老兵”——遗传算法GA来为MTP的矩张量收缩过程寻找最优的“装配流水线”。通过这套基于遗传算法的优化方案我们成功地将MTP模型在评估时的计算速度提升了近一个数量级并且在描述Ni-Al合金中各种棘手缺陷时精度还有所提升。如果你正在为大规模材料模拟的计算成本发愁或者对如何将智能优化算法与物理模型结合感兴趣那么接下来的内容或许能给你带来一些新的思路。2. 核心原理拆解MTP的“张量乐高”与遗传算法优化要理解我们的优化工作首先得弄明白MTP是怎么工作的以及传统方法的瓶颈在哪。这就像修车你得先知道发动机的构造才能找到提升动力的方法。2.1 矩张量势MTP的核心从原子环境到标量能量MTP的核心在于如何数学化地描述一个原子i的局部环境。它定义了一个称为矩张量描述符的量M_μ,ν(n_i) Σ_j [ f_μ,ν(|r_ij|, z_i, z_j) * f_c(|r_ij|) * (r_ij ⊗ ... ⊗ r_ij) ]这个公式看起来复杂但我们可以分部分理解Σ_j对中心原子i的所有邻居原子j求和。f_μ,ν(|r_ij|, z_i, z_j)这是一个径向基函数比如切比雪夫多项式它负责处理原子间距离|r_ij|和原子种类z的信息。μ是这个径向函数的索引。f_c(|r_ij|)截断函数保证在某个截断半径外相互作用平滑衰减至零这是所有经验势的常见做法以保证计算的局域性。(r_ij ⊗ ... ⊗ r_ij)这是最关键的张量部分。r_ij是从原子i指向原子j的向量。符号⊗表示克罗内克积一种张量积。ν表示这个向量被自乘了多少次也就是张量的阶数rank。物理图像你可以把ν0的矩张量想象成这个原子邻居的总“质量”标量ν1的矩张量一个向量描述了这个邻居团的“质心”位置ν2的矩张量一个矩阵则类似于“转动惯量”张量描述了质量分布的形状。更高阶的张量则包含了更精细的几何不对称性信息。MTP通过组合不同μ径向和ν角度的矩张量就能以多项式展开的形式逼近任意光滑的原子间相互作用势。单个矩张量M_μ,ν本身并不是旋转不变的旋转坐标系它的分量会变因此不能直接作为能量的贡献。MTP的巧妙之处在于它通过将多个不同阶的矩张量按特定规则进行缩并即对某些共同的索引进行求和最终生成一个旋转不变的标量B_α。这个缩并规则由一个对称索引矩阵α来定义。整个模型的原子能量U_i就是许多这样的标量基函数B_θ的线性组合U_i Σ_θ ξ_θ B_θ其中ξ_θ是通过拟合DFT数据得到的线性系数。2.2 传统收缩的瓶颈与“计算树”思想现在到了关键部分如何高效计算成千上万个这样的标量基函数B_α每个B_α都是多个矩张量缩并的结果。直接暴力计算所有可能的缩并索引计算量会随着缩并索引的数量指数增长3^mm为缩并索引数。Shapeev在提出MTP时想到了一个聪明的办法将每个标量函数B_α的缩并过程表示为一棵计算树。这棵树的叶子节点是最基础的矩张量M中间节点是缩并过程中产生的中间张量B_η根节点就是最终的标量B_α。计算时从叶子节点矩张量开始按照树的结构两两结合生成父节点中间张量最终得到根节点。核心洞察不同的标量函数B_α和B_β它们的计算树可能共享相同的中间子树即相同的中间张量B_η。例如B_α和B_β可能都需要先计算(M1·M2)这个中间结果。传统方法为每个标量函数独立构建计算树会重复计算这些共享的中间张量。因此优化的目标就变成了为整个基组所有需要计算的标量函数集合寻找一组全局最优的计算树使得需要存储和计算的独立中间张量分量总数、以及执行缩并运算的“乘法-加法”规则总数最小化。这本质上是一个复杂的组合优化问题。传统MTP采用了一种基于“最小化维度之和”和“最大化张量复用”的启发式贪婪算法来寻找这组树。但这种方法在面临大量选择时容易陷入局部最优并且当多个分解路径的“代价”相同时无法做出最优决策。2.3 遗传算法为“计算树”注入进化力量面对这个搜索空间巨大例如对于level-28的基组可能性超过10^2195种的优化问题我们引入了遗传算法GA。GA是一种受生物进化启发的全局优化算法特别适合解决这类组合爆炸问题。我们将整个优化过程设计如下基因编码每个标量函数B_α都有多个可能的计算树分解方式构成一个“树集合”{T_i^α}。我们为每个B_α分配一个“基因”这个基因的值k表示我们选择其树集合中的第k棵树。那么对于包含n个标量函数的整个基组一个完整的分解方案个体就可以用一个基因序列[k1, k2, ..., kn]来表示。种群初始化随机生成一定数量如40个的个体构成初始种群。为了加速收敛我们并非完全随机而是给每棵树一个初始“评分”例如选择维度之和最小的树、能复用已有中间张量的树会得分更高优先选择高得分的树来构建初始个体。适应度函数这是进化的“指挥棒”。我们定义了一个成本函数F(R) m 3t其中R代表一个分解方案m是独立张量分量的总数t是计算规则的总数。给t乘以权重3是因为在实际计算中执行这些规则乘加运算的开销通常比存储张量更大。适应度函数f(R)被设计为与F(R)负相关且采用类似费米-狄拉克分布的形式确保种群多样性避免过早收敛于次优解。选择、交叉与变异选择根据个体的适应度采用“轮盘赌”等方式选择优秀的个体进入下一代。交叉随机选取两个父代个体交换它们的一部分基因即交换一部分标量函数所选择的计算树产生新的子代个体。这模拟了基因重组。变异以一定概率随机改变某个个体中某个基因的值即为某个标量函数随机更换一棵计算树。这为种群引入了新的可能性避免陷入局部最优。迭代与收敛重复进行选择、交叉、变异并计算新一代的适应度如此迭代数万乃至数十万代。直到达到预设的最大迭代次数或最优解在连续多代中不再改善。冲突解决与后优化在进化过程中一个关键挑战是基因间的“冲突”。例如标量函数A和B可能都需要用到同一个中间张量X但A选择的计算树要求X由(M1·M2)生成而B选择的树却要求X由(M3·M4)生成这就产生了矛盾。我们的算法在生成计算规则时会强制统一这些冲突确保一致性。在GA搜索结束后我们还会执行一个快速的贪婪算法进行局部微调从低阶到高阶张量逐一尝试替换为更优的计算树确保得到最终的最优方案。通过这套流程遗传算法能够在浩如烟海的分解方案中高效地寻找到那个能最大程度复用中间结果、最小化计算量的全局优解。图2清晰地对比了传统方案与我们优化后方案的工作流差异传统方案是“各自为政”的独立计算而我们的方案则是经过全局规划后的“协同流水线”。3. 基组设计策略在精度与效率间寻找新平衡在优化了“怎么算”之后下一个问题是“算什么”即如何选择构成势函数的标量基函数B_α的集合基组。基组决定了模型的表达能力和复杂度。传统MTP使用一个称为“水平level”的参数lev(M_μ,ν) 4μ ν 2来筛选基函数只保留总水平低于某个阈值d的基函数。然而我们通过分析发现计算成本主要来自两方面计算矩张量分量本身对于一个ν阶张量其独立分量数以ν^2规模增长。计算每个分量都需要遍历中心原子的所有邻居而邻居数量与截断半径rc的立方成正比。这部分计算无法通过优化缩并过程来减少。执行缩并规则这正是我们上一节用遗传算法优化的部分。因此传统的水平定义4μ ν 2对μ径向复杂度的惩罚因子4过高导致为了达到高多项式精度不得不引入高阶ν大的矩张量而这会急剧增加第一部分的计算量。我们的策略是重新定义水平lev(M_μ,ν) 2μ ν 1。我们降低了对径向复杂度μ的惩罚从4降到2但严格限制了张量的阶数ν在我们的工作中通常ν ≤ 3。这意味着我们使用更多低阶ν小的矩张量但通过让更多这样的张量参与缩并即考虑更高阶的多体相互作用如五体、七体来提升模型的表达能力。这样做的优势非常明显低阶张量的分量数少计算成本低。虽然增加参与缩并的张量数量会增加第二部分缩并的复杂度但这部分恰恰是我们用遗传算法重点优化的部分。实测表明这种“低阶张量高阶相互作用优化缩并”的策略能在不损失精度甚至有所提升的前提下带来数量级的计算加速。例如我们定义的“2653”基组包含2653个线性参数考虑至五体相互作用ν_max4其精度与传统的level-28基组相当但计算速度却快了近11倍。4. 实战构建用于Ni-Al合金缺陷模拟的高效MTP势理论再好也需要实战检验。我们选择Ni-Al合金体系作为试金石因为其内部的γ/γ’相结构以及各类缺陷空位、反位、层错等对力学性能至关重要也是检验势函数准确性的经典难题。4.1 训练集构建主动学习确保“见过世面”一个MLIP的可靠性极大程度上取决于其训练集是否足够全面和具有代表性。我们采用了一种两阶段主动学习策略来构建高质量训练集第一阶段基于贝叶斯误差的在线学习。我们使用VASP软件进行从头算分子动力学AIMD模拟。在模拟过程中集成一个基于高斯过程回归GAP的在线学习模块。这个模块会实时评估当前原子构型对于当前GAP模型是否属于“未知领域”即预测不确定性高。一旦发现这类构型就暂停MD调用DFT计算其精确的能量、力和应力并将其加入训练集然后更新GAP模型。这个过程能自动探索相空间捕获动力学过程中出现的稀有但重要的构型。第二阶段基于D-最优性准则的主动采样。用第一阶段得到的训练集初步拟合一个MTP势。然后我们用这个MTP势进行更长时间的MD模拟200 ps生成大量新的候选构型。接着利用Shapeev提出的广义D-最优性准则从这些候选构型中筛选出那些能最大程度提升模型泛化能力即最大程度降低参数不确定性的构型进行DFT计算并加入训练集。随后重新拟合MTP。如此循环直到没有新的构型能显著扩展模型的认知边界。最终我们得到了一个包含8450个构型的训练集涵盖了fcc-Ni、L12-Ni3Al、相界面、各类堆垛层错GSF以及它们与空位、空位团、反位缺陷的组合构型。同时我们构建了一个独立的包含815个构型的验证集用于评估模型的泛化能力。通过主成分分析PCA可以看到验证集很好地覆盖了训练集所占据的构型空间。4.2 势函数拟合两步法求精在确定了优化的2653基组后我们使用完整的训练集进行最终拟合。拟合过程分为两步以确保找到全局最优的线性系数ξ_θ初始化线性系数首先从完整的2653基组中提取一个最小的子集例如只包含最重要的基函数用这个子集去拟合训练数据。拟合得到的势函数其径向基函数是确定的。然后我们固定这些径向基函数在完整的2653基组上以子集拟合的结果作为初始点进行一次线性最小二乘优化得到一组较好的线性系数初值。这一步的目的是避免直接在高维参数空间中随机初始化可能导致的收敛缓慢或陷入局部极小。非线性精修以上一步得到的系数为起点采用L-BFGS算法进行最多5000次迭代的非线性优化进一步调整所有参数包括线性系数和径向基函数的参数。L-BFGS迭代完成后再执行一次最终的线性最小二乘拟合对线性系数进行微调。在拟合过程中能量、力和应力的权重分别设置为1、0.01和0.005。这个权重的设置很有讲究力的数据量远大于能量每个构型有3N个力分量适当降低其权重可以防止力误差主导优化过程而应力数据则更少权重也最低。4.3 效率与精度评估优化带来的显著提升我们对优化前后的不同基组进行了全面的基准测试测试在一个2048原子的Ni3Al超胞上进行10万步MD模拟。结果令人振奋参见表I对于简单基组如level-18优化将中间张量数量、独立分量数和计算规则数减少了约8%-19%。但由于此时计算成本主要由计算矩张量分量本身主导MD模拟速度仅提升了约4%。对于复杂基组如level-28优化的效果极为显著。中间张量数量减少18.2%独立分量数减少41.3%计算规则数减少43.9%。这直接转化为71%的MD模拟加速时间从4418秒减少到2580秒。这是因为在复杂基组下缩并操作的成本占据了主导我们的优化直击要害。我们的新基组2653 vs 传统基组level-28这是最关键的对比。两者具有相近的线性参数数量~2500个但我们的2653基组由于采用了低阶张量设计其MD模拟速度是level-28基组的11倍。更令人惊喜的是在验证集上2653基组的能量、力、应力的预测误差均低于level-28基组。这证明了我们“降低张量阶数、增加多体阶数、优化收缩”策略的成功它在降低计算复杂度的同时反而获得了更好的拟合能力。表II展示了训练阶段的成本。在L-BFGS迭代过程中优化带来的加速比与MD模拟中观察到的趋势基本一致。对于线性拟合步骤由于计算成本主要来自求解大型线性方程组优化带来的加速不明显但对于小基组level-18仍有提升。5. 性能验证Ni-Al合金缺陷模拟的实战表现一个势函数好不好最终要看它解决实际科学问题的能力。我们将优化后的MTP势基于2653基组应用于Ni-Al合金的一系列缺陷性质预测并与DFT结果以及经典的EAM势进行对比。5.1 基础性质晶格、弹性与声子首先我们检验了模型对完美晶体基本性质的再现能力晶格常数与能量-体积曲线MTP势精确地复现了Ni和Ni3Al的平衡晶格常数以及它们在压缩和膨胀时的能量变化曲线E-V曲线与DFT结果高度吻合。弹性常数对于立方晶体关键的C11, C12, C44等弹性常数MTP的预测误差在几个GPa以内远优于EAM势这对于模拟材料的力学响应至关重要。声子色散谱我们通过冻结声子法计算了Ni和Ni3Al的声子谱。MTP势不仅准确预测了声子支的整体走势甚至复现了DFT计算中一些细微的“扭结”kinks这表明它很好地捕捉到了原子间相互作用的细节。而EAM势通常在高频部分偏差较大。5.2 点缺陷空位与反位缺陷点缺陷是合金中最基本的缺陷类型。空位形成能我们计算了Ni和Ni3Al中单个空位的形成能。MTP的结果与DFT值几乎一致而EAM势通常系统性偏低或偏高。空位团结合能这是更严峻的考验。我们计算了从二空位到小空位团如四面体空位团的结合能。MTP势准确地描述了空位之间的相互作用随距离和构型的变化趋势成功预测了DFT计算中发现的某些特定构型的相对稳定性。这对于研究辐照或淬火过程中空位的聚集行为至关重要。反位缺陷形成能在Ni3Al中一个Al原子占据Ni的位置Ni位反位或一个Ni原子占据Al的位置Al位反位其形成能不同。MTP势精确区分了这两种能量与DFT相符而EAM势往往难以准确描述这种化学有序性。5.3 面缺陷堆垛层错与相界面面缺陷直接影响材料的塑性变形机制。广义堆垛层错能GSFE曲线我们计算了Ni和Ni3Al中不同滑移面如{111}面上的GSFE曲线。这条曲线决定了位错分解的宽度和层错能是计算材料强度的关键输入。MTP势给出的GSFE曲线与DFT计算的重合度极高特别是内禀层错能ISF和非内禀层错能USF的数值。相比之下不同EAM势给出的结果分散且常与DFT有较大偏差。γ/γ’相界面能Ni基高温合金的强度源于γ相基体中弥散分布的γ’Ni3Al强化相。两相之间的界面能是控制组织形貌的关键参数。MTP势计算得到的界面能与DFT参考值吻合良好。5.4 动力学性质空位扩散最后我们通过分子动力学模拟研究了空位在Ni中的扩散行为。我们统计了空位的均方位移随时间的变化计算了扩散系数。MTP势模拟得到的扩散系数与激活能与实验值及基于DFT的过渡态理论计算结果处于合理范围内。这表明优化后的MTP势不仅能描述静态性质还能可靠地用于研究有限温度下的动力学过程。6. 实操心得与避坑指南基于这个项目的完整开发和应用流程我总结了一些对于从事MLIP开发或应用的研究者可能极具价值的实操经验。6.1 遗传算法参数调优耐心与平衡种群大小与迭代次数种群大小不宜过小如小于20否则多样性不足容易早熟收敛也不宜过大如大于100否则每次迭代的计算开销剧增。我们发现在40-60之间是个不错的起点。迭代次数需要足够多对于复杂基组如461310万次迭代是必要的。可以观察最优个体适应度随代数的变化曲线当曲线进入漫长的平台期时即可考虑停止。适应度函数设计成本函数F(R) m w * t中的权重w是关键。我们通过 profiling 发现在目标硬件上执行一条“乘法-加法”规则对应t的成本大约是访问一个张量分量对应m成本的3倍左右因此设置了w3。这个权重需要根据具体的代码实现和硬件架构进行微调。如果代码对内存访问做了极致优化w值可能降低如果计算规则部分存在瓶颈w值可能需要提高。交叉与变异概率过高的交叉概率会导致优良基因模式被破坏过高变异概率则会使搜索退化为随机游走。我们采用自适应策略初期使用较高的交叉概率如0.8和中等变异概率如0.1以快速探索空间后期降低交叉概率如0.5略微提高变异概率如0.15以加强局部搜索和跳出局部最优的能力。6.2 训练集构建质量重于数量主动学习是必须的而非可选的单纯从完美晶体、弹性变形等常规构型采样得到的训练集其势函数在模拟缺陷、表面、熔化等“非典型”过程时极易外推失败extrapolation。在线学习on-the-fly和基于D-最优性准则的采样是确保模型“见过世面”的最有效手段。不要吝啬在主动学习循环中花费的DFT计算资源这能从根本上避免后续大规模模拟得到错误结果。警惕“数据分布不均”在主动学习过程中可能会在某个区域如熔融状态采集到过多构型。需要在构建最终训练集时进行适当的下采样或者为不同区域的构型设置不同的权重防止模型过度拟合某一类构型而忽略其他。必须包含应力数据对于固体材料尤其是要计算弹性常数或模拟力学变形训练集中必须包含DFT计算的应力张量信息。这对正确描述体积相关的能量变化至关重要。我们的经验是应力数据的权重可以设得比力数据低一个数量级。6.3 基组选择不是越复杂越好从简开始逐步增加不要一开始就使用参数最多的基组。从一个中等水平的基组如我们的“2653”开始训练和测试。如果验证误差已经满足要求且进一步增加基组复杂度如到“4613”带来的精度提升微乎其微但计算成本显著增加那么就应该停止。模型的“奥卡姆剃刀”原则同样适用。关注截断半径截断半径rc需要仔细选择。太小会丢失重要的长程相互作用对于某些合金或半导体很重要太大会急剧增加计算量邻居数 ~rc^3。通常rc需要至少包含第二近邻或第三近邻原子。可以通过绘制径向分布函数RDF或测试不同rc下对结合能、弹性常数的影响来确定。验证“能量-力-应力”的一致性一个健康的势函数其预测的能量、力和应力应该是自洽的。一个简单的检查方法是对一个平衡构型施加一个微小应变用势函数计算能量变化和应力同时用有限差分法通过位移原子计算力再推导应力来验证。两者应该基本一致。6.4 应用模拟中的注意事项外推预警所有MLIP都内置了外推预警机制通常基于训练集构型在描述符空间的距离。在运行MD时务必监控这个“外推指标”extrapolation grade。一旦有原子构型的指标超过阈值比如我们设为5模拟就不可信了。此时应该停止模拟将该构型加入训练集重新拟合势函数。温度与相稳定性用MLIP进行高温MD模拟前最好先测试一下它在高温下对相稳定性的预测。例如对于Ni-Al合金可以跑一个升温模拟看γ和γ’相是否能在正确的温度范围内保持稳定以及熔点是否合理。我们的MTP势在预测Ni和Ni3Al熔点方面表现良好这是许多EAM势的短板。缺陷模拟的有限尺寸效应在计算缺陷形成能时超胞尺寸必须足够大以消除周期性镜像缺陷之间的相互作用。一个实用的方法是逐步增大超胞尺寸计算缺陷形成能直到其变化小于目标精度如0.01 eV。对于空位、层错等通常需要包含数百个原子的超胞。通过将遗传算法这一全局优化工具深度融入MTP势的构建流程我们不仅获得了一个在计算速度上提升近一个数量级的高效势函数更重要的是建立了一套可复用的、系统化的势函数开发和优化方法论。这项工作表明在计算材料科学中对底层数学模型和算法进行“精雕细琢”的优化与采集更多数据、使用更复杂的网络结构同样重要甚至能带来事半功倍的效果。对于Ni-Al合金乃至更复杂的多元合金体系这种高效的模拟工具将为在原子尺度上深入理解缺陷演化、相变动力学和力学性能打开一扇新的大门。
