遗传算法优化矩张量势：提升材料模拟效率与精度

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^m，m为缩并索引数）。

Shapeev在提出MTP时，想到了一个聪明的办法：将每个标量函数B_α的缩并过程表示为一棵计算树。这棵树的叶子节点是最基础的矩张量M，中间节点是缩并过程中产生的中间张量B_η，根节点就是最终的标量B_α。计算时，从叶子节点（矩张量）开始，按照树的结构，两两结合生成父节点（中间张量），最终得到根节点。

核心洞察：不同的标量函数B_α和B_β，它们的计算树可能共享相同的中间子树（即相同的中间张量B_η）。例如，B_α和B_β可能都需要先计算(M1·M2)这个中间结果。传统方法为每个标量函数独立构建计算树，会重复计算这些共享的中间张量。

因此，优化的目标就变成了：为整个基组（所有需要计算的标量函数集合）寻找一组全局最优的计算树，使得需要存储和计算的独立中间张量分量总数、以及执行缩并运算的“乘法-加法”规则总数最小化。这本质上是一个复杂的组合优化问题。传统MTP采用了一种基于“最小化维度之和”和“最大化张量复用”的启发式贪婪算法来寻找这组树。但这种方法在面临大量选择时，容易陷入局部最优，并且当多个分解路径的“代价”相同时，无法做出最优决策。

2.3 遗传算法：为“计算树”注入进化力量

面对这个搜索空间巨大（例如，对于level-28的基组，可能性超过10^2195种）的优化问题，我们引入了遗传算法（GA）。GA是一种受生物进化启发的全局优化算法，特别适合解决这类组合爆炸问题。

我们将整个优化过程设计如下：

1. 基因编码：每个标量函数B_α都有多个可能的计算树（分解方式）构成一个“树集合”{T_i^α}。我们为每个B_α分配一个“基因”，这个基因的值k表示我们选择其树集合中的第k棵树。那么，对于包含n个标量函数的整个基组，一个完整的分解方案（个体）就可以用一个基因序列[k1, k2, ..., kn]来表示。

2. 种群初始化：随机生成一定数量（如40个）的个体，构成初始种群。为了加速收敛，我们并非完全随机，而是给每棵树一个初始“评分”（例如，选择维度之和最小的树、能复用已有中间张量的树会得分更高），优先选择高得分的树来构建初始个体。

3. 适应度函数：这是进化的“指挥棒”。我们定义了一个成本函数F(R) = m + 3t，其中R代表一个分解方案，m是独立张量分量的总数，t是计算规则的总数。给t乘以权重3，是因为在实际计算中，执行这些规则（乘加运算）的开销通常比存储张量更大。适应度函数f(R)被设计为与F(R)负相关，且采用类似费米-狄拉克分布的形式，确保种群多样性，避免过早收敛于次优解。

4. 选择、交叉与变异：

   • 选择：根据个体的适应度，采用“轮盘赌”等方式，选择优秀的个体进入下一代。
   • 交叉：随机选取两个父代个体，交换它们的一部分基因（即交换一部分标量函数所选择的计算树），产生新的子代个体。这模拟了基因重组。
   • 变异：以一定概率随机改变某个个体中某个基因的值（即为某个标量函数随机更换一棵计算树）。这为种群引入了新的可能性，避免陷入局部最优。

5. 迭代与收敛：重复进行选择、交叉、变异，并计算新一代的适应度，如此迭代数万乃至数十万代。直到达到预设的最大迭代次数，或最优解在连续多代中不再改善。

6. 冲突解决与后优化：在进化过程中，一个关键挑战是基因间的“冲突”。例如，标量函数A和B可能都需要用到同一个中间张量X，但A选择的计算树要求X由(M1·M2)生成，而B选择的树却要求X由(M3·M4)生成，这就产生了矛盾。我们的算法在生成计算规则时，会强制统一这些冲突，确保一致性。在GA搜索结束后，我们还会执行一个快速的贪婪算法进行局部微调，从低阶到高阶张量逐一尝试替换为更优的计算树，确保得到最终的最优方案。

通过这套流程，遗传算法能够在浩如烟海的分解方案中，高效地寻找到那个能最大程度复用中间结果、最小化计算量的全局优解。图2清晰地对比了传统方案与我们优化后方案的工作流差异：传统方案是“各自为政”的独立计算，而我们的方案则是经过全局规划后的“协同流水线”。

3. 基组设计策略：在精度与效率间寻找新平衡

在优化了“怎么算”之后，下一个问题是“算什么”，即如何选择构成势函数的标量基函数B_α的集合（基组）。基组决定了模型的表达能力和复杂度。传统MTP使用一个称为“水平（level）”的参数lev(M_μ,ν) = 4μ + ν + 2来筛选基函数，只保留总水平低于某个阈值d的基函数。

然而，我们通过分析发现，计算成本主要来自两方面：

1. 计算矩张量分量本身：对于一个ν阶张量，其独立分量数以ν^2规模增长。计算每个分量都需要遍历中心原子的所有邻居，而邻居数量与截断半径rc的立方成正比。这部分计算无法通过优化缩并过程来减少。
2. 执行缩并规则：这正是我们上一节用遗传算法优化的部分。

因此，传统的水平定义4μ + ν + 2对μ（径向复杂度）的惩罚因子（4）过高，导致为了达到高多项式精度，不得不引入高阶（ν大）的矩张量，而这会急剧增加第一部分的计算量。

我们的策略是重新定义水平：lev(M_μ,ν) = 2μ + ν + 1。我们降低了对径向复杂度μ的惩罚（从4降到2），但严格限制了张量的阶数ν（在我们的工作中通常ν ≤ 3）。这意味着，我们使用更多低阶（ν小）的矩张量，但通过让更多这样的张量参与缩并（即考虑更高阶的多体相互作用，如五体、七体）来提升模型的表达能力。

这样做的优势非常明显：低阶张量的分量数少，计算成本低。虽然增加参与缩并的张量数量会增加第二部分（缩并）的复杂度，但这部分恰恰是我们用遗传算法重点优化的部分。实测表明，这种“低阶张量+高阶相互作用+优化缩并”的策略，能在不损失精度甚至有所提升的前提下，带来数量级的计算加速。例如，我们定义的“2653”基组（包含2653个线性参数，考虑至五体相互作用，ν_max=4），其精度与传统的level-28基组相当，但计算速度却快了近11倍。

4. 实战：构建用于Ni-Al合金缺陷模拟的高效MTP势

理论再好，也需要实战检验。我们选择Ni-Al合金体系作为试金石，因为其内部的γ/γ’相结构以及各类缺陷（空位、反位、层错等）对力学性能至关重要，也是检验势函数准确性的经典难题。

4.1 训练集构建：主动学习确保“见过世面”

一个MLIP的可靠性，极大程度上取决于其训练集是否足够全面和具有代表性。我们采用了一种两阶段主动学习策略来构建高质量训练集：

1. 第一阶段：基于贝叶斯误差的在线学习。我们使用VASP软件进行从头算分子动力学（AIMD）模拟。在模拟过程中，集成一个基于高斯过程回归（GAP）的在线学习模块。这个模块会实时评估当前原子构型对于当前GAP模型是否属于“未知领域”（即预测不确定性高）。一旦发现这类构型，就暂停MD，调用DFT计算其精确的能量、力和应力，并将其加入训练集，然后更新GAP模型。这个过程能自动探索相空间，捕获动力学过程中出现的稀有但重要的构型。

2. 第二阶段：基于D-最优性准则的主动采样。用第一阶段得到的训练集初步拟合一个MTP势。然后，我们用这个MTP势进行更长时间的MD模拟（200 ps），生成大量新的候选构型。接着，利用Shapeev提出的广义D-最优性准则，从这些候选构型中筛选出那些能最大程度提升模型泛化能力（即最大程度降低参数不确定性）的构型，进行DFT计算并加入训练集。随后重新拟合MTP。如此循环，直到没有新的构型能显著扩展模型的认知边界。

最终，我们得到了一个包含8450个构型的训练集，涵盖了fcc-Ni、L12-Ni3Al、相界面、各类堆垛层错（GSF），以及它们与空位、空位团、反位缺陷的组合构型。同时，我们构建了一个独立的包含815个构型的验证集，用于评估模型的泛化能力。通过主成分分析（PCA）可以看到，验证集很好地覆盖了训练集所占据的构型空间。

4.2 势函数拟合：两步法求精

在确定了优化的2653基组后，我们使用完整的训练集进行最终拟合。拟合过程分为两步，以确保找到全局最优的线性系数ξ_θ：

1. 初始化线性系数：首先，从完整的2653基组中提取一个最小的子集（例如，只包含最重要的基函数），用这个子集去拟合训练数据。拟合得到的势函数，其径向基函数是确定的。然后，我们固定这些径向基函数，在完整的2653基组上，以子集拟合的结果作为初始点，进行一次线性最小二乘优化，得到一组较好的线性系数初值。这一步的目的是避免直接在高维参数空间中随机初始化可能导致的收敛缓慢或陷入局部极小。

2. 非线性精修：以上一步得到的系数为起点，采用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之间是个不错的起点。迭代次数需要足够多，对于复杂基组（如4613），10万次迭代是必要的。可以观察最优个体适应度随代数的变化曲线，当曲线进入漫长的平台期时，即可考虑停止。
• 适应度函数设计：成本函数F(R) = m + w * t中的权重w是关键。我们通过 profiling 发现，在目标硬件上，执行一条“乘法-加法”规则（对应t）的成本，大约是访问一个张量分量（对应m）成本的3倍左右，因此设置了w=3。这个权重需要根据具体的代码实现和硬件架构进行微调。如果代码对内存访问做了极致优化，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合金乃至更复杂的多元合金体系，这种高效的模拟工具将为在原子尺度上深入理解缺陷演化、相变动力学和力学性能，打开一扇新的大门。