机器学习势函数与主动学习：解析粗糙铜-水界面原子尺度催化活性位点

1. 项目概述：当机器学习“遇见”粗糙的铜-水界面

在电化学催化，尤其是二氧化碳还原反应（CO2RR）的研究中，铜催化剂表面的微观结构一直是个“黑箱”。教科书里光滑平整的（100）、（111）晶面模型，在真实的反应环境中往往不复存在。实验观察告诉我们，在反应电位下，铜表面会动态重构，形成从埃米级到纳米级不等的粗糙结构，布满台阶、扭折和角落原子。这些“不完美”的位点，恰恰被认为是催化活性提升的关键。然而，想从原子层面看清水分子在这些复杂粗糙表面上的吸附、排列和反应行为，传统计算模拟方法显得力不从心。密度泛函理论（DFT）虽然精度高，但算力消耗巨大，通常只能处理几百个原子的理想模型，面对数万乃至数十万原子的真实粗糙界面，只能望洋兴叹。

这就是机器学习势函数（MLIP）大显身手的地方。它像一位技艺高超的“翻译官”，通过学习大量高精度的DFT计算结果，构建出一个既能保持量子力学精度、又能进行大规模分子动力学（MD）模拟的“代理模型”。我们这次的工作，核心就是构建一个专门针对“粗糙铜-水界面”的高精度MLIP，并利用它进行大规模模拟，目的就是掀开这个“黑箱”的盖子，看看在真实的、粗糙的铜表面上，水分子到底在干什么，哪些位点才是真正的“活性中心”。这不仅仅是算力的提升，更是一种研究范式的转变——从研究理想的“模型系统”转向模拟更接近实验条件的“真实系统”。

2. 核心思路与方案设计：如何为“粗糙”定制一个势函数

构建一个可靠的MLIP，最关键的挑战在于训练数据的质量与代表性。如果你的训练数据只包含光滑表面的构型，那么训练出的模型在预测粗糙表面的原子相互作用时，其可靠性必然存疑。这就好比只学过平原驾驶的自动驾驶模型，突然被扔到了崎岖的山路上。因此，我们的核心思路是设计一个“主动学习”工作流，让模型自己告诉我们它在哪里“没把握”，然后我们针对性地为这些“知识盲区”补充DFT数据。

2.1 粗糙铜表面的生成：从“理想”到“现实”

在模拟水之前，我们首先需要创建一系列具有不同粗糙度的铜表面。我们采用了两种物理上合理的“造景”方法，而不是简单地随机扰动原子位置。

第一种是“纳米颗粒熔融-再结晶法”。想象一下，我们在一个平整的铜基底上随机撒上一些大小不一的铜纳米颗粒（平均直径10-30 Å）。然后，在分子动力学模拟中，我们只加热基底的上半部分和这些颗粒，让它们熔化并流动，而下半部分保持低温结晶态作为“地基”。随后，我们快速淬火（降温），熔融的部分会重新结晶，但过程受初始颗粒分布和温度场影响，自然形成凹凸不平的粗糙表面。通过调整纳米颗粒的大小、最高熔融温度和淬火速率，我们可以控制最终表面的粗糙度。

第二种是“压痕诱导变形法”。我们从一个平整的铜表面开始，加热其顶层至高温（1500 K），使其软化。然后，引入一系列虚拟的“压头”（通过纯排斥的伦纳德-琼斯势函数模拟），将这些压头压入软化的表层，像用手指在黏土上按压一样，将表层原子推开，形成凹陷。保持压头位置的同时淬火，让变形的表层再结晶固定，从而形成粗糙结构。通过改变压头的分布、下压深度和半径，我们可以创造出各种形态的粗糙度。

我们最终生成了46个不同粗糙度的表面，其均方根粗糙度（Rq）范围从约2 Å到18 Å，覆盖了实验观测到的典型范围。这为我们后续研究提供了一个具有统计意义的表面模型库。

2.2 主动学习工作流：让模型“主动提问”

有了粗糙表面，接下来的核心任务是为铜-水界面收集训练数据。我们采用了一个迭代式的主动学习循环，其精妙之处在于“空间分辨不确定性”的利用。

1. 初代模拟与不确定性评估：我们使用一个初步的MLIP（基于已有的光滑界面数据库训练）对选定的粗糙铜-水界面进行大规模MD模拟。模拟过程中，我们使用一个“委员会模型”（由多个略有差异的MLIP组成）。对于一个给定的原子构型，如果委员会中所有模型的预测结果高度一致，说明模型对这个构型很“有把握”；如果预测结果分歧很大，则说明此处存在高“不确定性”，可能是模型未曾学习过的全新局部环境。

2. 高不确定性构型的提取与“包装”：我们识别出不确定性最高的那些原子（通常都位于复杂的界面处）。直接对这些包含数万原子的快照进行DFT计算是不可行的。因此，我们需要进行“手术式”提取：以高不确定性原子为中心，截取一个包含其周围一定截断半径（如4 Å）内所有原子的“局部环境”小盒子。但这个小盒子边界处被生硬切断，会产生不合理的悬挂键和真空层，无法直接用于周期性边界条件的DFT计算。

3. 非晶矩阵嵌入与两步退火：为了解决边界问题，我们采用了改进的“非晶矩阵嵌入”方法。将提取出的“局部环境”原子放入一个稍大的计算盒子后，我们执行一个两步退火流程：

   • 第一步：铜层重构。保持“核心关注区域”（即我们想研究的局部环境）和水分子固定，将盒子边缘的铜原子加热到较高温度（如600 K），允许它们移动、重组，从而消除不合理的真空间隙，形成一个连续（尽管可能是非晶态）的铜基质。
   • 第二步：水层弛豫。保持核心区域和铜原子固定，让水在室温下平衡，填充剩余空间，调整水分子取向以形成合理的氢键网络。 这个过程就像为珍贵的局部环境样本（一块奇石）精心制作一个稳固的底座（铜）和填充物（水），使其能够独立、稳定地接受DFT“鉴定”，同时保证核心样本的原始状态不被破坏。

4. DFT计算与数据库更新：对“包装”好的小体系进行DFT计算，获得精确的能量和力。将这些新的数据加入训练数据库。

5. 迭代优化：用更新后的数据库重新训练MLIP，然后用新模型再次进行大规模MD模拟，评估不确定性，并开始下一轮采样。如图4所示，经过几轮迭代后，模型在整个粗糙界面上的平均不确定性显著下降，说明其“知识盲区”被有效填补。

通过这个工作流，我们最终采样了238个高不确定性的界面构型，并补充了38个中等不确定性的构型作为测试集。这些数据共同构成了一个能充分描述粗糙铜-水界面复杂性的高质量数据库。

2.3 势函数模型选型：在精度与效率间权衡

有了数据库，下一步是选择合适的MLIP架构进行最终训练。我们主要对比了两种主流框架：原子簇展开（ACE）和图原子簇展开（GRACE）。

• GRACE 2-Layer：精度最高，但计算成本也最高，且受限于实现，通常只能处理小于1万原子的体系。我们的目标体系包含约20万个原子，因此它不适合作为生产模拟的势函数。
• ACE：计算效率极高，比GRACE快一个数量级以上，非常适合大规模模拟。但在我们的测试中，某些ACE模型（特别是截断半径较小的）在描述阶梯状表面（如(322), (433)）的水密度分布时，与高精度参考结果存在细微偏差，例如第一个化学吸附峰的形状再现不够准确。
• GRACE 1-Layer：在精度和效率之间取得了最佳平衡。虽然计算成本高于ACE，但仍在可接受范围内（能处理20万原子体系），并且其对阶梯表面水结构的描述与GRACE 2-Layer参考结果吻合得非常好。

考虑到我们的核心目标是精确解析粗糙表面（富含台阶、边缘）上的水分子行为，我们最终选择了截断半径为6 Å的GRACE 1-Layer模型作为最终模拟的势函数。这个选择背后是研究需求的直接体现：对于界面模拟，结构细节的准确性往往比绝对能量误差的微小差异更为重要。

3. 粗糙界面结构解析：无监督学习“看见”隐藏的图案

使用训练好的GRACE 1-Layer势函数，我们对一个具有代表性粗糙度（Rq = 4 Å）的铜-水界面进行了大规模分子动力学模拟。传统的分析方法是观察沿表面法线方向的平均水密度分布曲线。如图5c所示，由于表面凹凸不平，水密度峰出现在从最低到最高的多个铜原子层之间，这与光滑表面的单一尖锐峰截然不同。第一个氧峰出现在约2.9 Å处，对应于化学吸附水，这与之前光滑表面的研究一致。然而，这种全局平均的描述过于笼统，无法揭示粗糙表面上千变万化的局部环境。

为了深入挖掘，我们采用了无监督机器学习方法对界面铜原子的局部环境进行分类。具体来说，我们使用球贝塞尔描述符对每个表层铜原子及其周围环境进行数学编码，将高维的结构信息转换为一个固定长度的特征向量。然后，利用UMAP降维技术将这些高维向量投影到二维平面上，相似的局部环境会在投影图上聚集在一起。

3.1 发现“隐藏”的活性位点集群

分析结果令人振奋（图6）。我们将粗糙界面上的表层铜原子局部环境划分成了19个不同的簇（Cluster）。

• 来自“教科书”的熟悉面孔：一部分簇（如簇8, 9）对应于光滑的(111)晶面环境（有/无化学吸附水），另一部分（如簇3, 4）则与之前研究过的阶梯表面（如(211)）上的边缘原子环境高度相似。这说明我们的粗糙表面包含了所有已知模型表面的结构基元。
• 粗糙表面独有的“新角色”：更重要的是，我们发现了多个在任何理想模型表面中都不存在的簇。这主要包括三类：
  1. 低配位角落原子（簇1, 2, 6）：这些原子位于多个台阶边缘的交汇处，平面内的铜原子邻居数极少（只有2-3个），是粗糙表面特有的、配位严重不足的位点。
  2. 堆垛层错相关环境（簇13, 14, 16）：由表面制备过程中的晶体结构缺陷导致。
  3. 表面空位底部的原子（簇19）。

3.2 化学吸附行为的原子级洞察

无监督分类不仅区分了几何结构，还敏锐地捕捉到了化学状态的差异。例如，簇3和簇4在UMAP投影上位置很近，说明它们的局部铜原子排列非常相似（都是阶梯边缘）。但通过计算径向分布函数（RDF）和配位数分析（图7），我们发现了一个关键区别：簇3中的所有铜原子上都稳定地化学吸附着一个水分子（氧配位数≈1），而簇4的原子上则没有化学吸附水（氧配位数≈0）。这种同构不同性的区分，完全是由数据驱动发现的，无需我们预先定义“是否有水吸附”这个标签。

最有趣的发现来自于那些独特的“角落原子”簇（簇1和簇2）。在所有被归类到簇1和簇2的角落原子上，我们都观察到了持续、稳定的水分子化学吸附。而且，在粗糙表面上，我们没有发现任何具有类似低配位几何结构但却没有水吸附的铜原子。这表明，对于这类极低配位的角落位点，水分子化学吸附可能是一种热力学上非常有利的必然行为。

进一步分析RDF发现，簇1（配位更低）的第一个氧峰比簇2更尖锐、更强，意味着水分子与这些位点的结合更强。这种增强的吸附作用，很可能使得第二个水分子更容易在3 Å范围内接近吸附位点，这在RDF的次近邻峰特征中有所体现。这些强烈吸附水分子的低配位角落原子，因其极高的反应活性，极有可能是二氧化碳还原等电化学反应的关键活性位点。

3.3 回望主动学习：它学到了什么？

我们还可以将主动学习过程中采样的238个高不确定性构型中的中心原子，映射到上述19个簇的UMAP图谱中（图8）。结果显示，绝大多数采样点都落在了簇1、2、3这三个区域。这完美验证了主动学习策略的有效性：

• 簇1和簇2（角落原子）在初始训练数据库中完全不存在，模型对其一无所知，因此不确定性最高，被优先采样。
• 簇3（有化学吸附水的阶梯边缘）被大量采样则有些出乎意料，因为类似的阶梯边缘环境在初始数据库中本应存在。这可能意味着，粗糙表面上的阶梯边缘其局部应变或电子结构发生了细微变化，导致水吸附行为与理想模型有差异，从而触发了模型的不确定性。

而那些同样为粗糙表面独有的、但位于堆垛层错或空位底部的原子（如簇13, 14, 19），则很少被采样。这是因为这些位点更多地暴露给体相铜，化学环境相对简单，且没有发生强烈的水吸附，因此模型对其预测的不确定性较低。

4. 实操要点、挑战与经验总结

这个项目涉及从表面生成、势函数开发到模拟分析的完整链条，每一步都有需要特别注意的“坑”。

4.1 势函数训练与验证的陷阱

1. 测试集构建的独立性：绝对不能从主动学习的采样轨迹中随机抽取构型作为测试集！因为模型在训练过程中已经“见过”这些轨迹的整体特征。正确的做法是：a) 保留一部分初始数据；b) 从后续模拟的新轨迹（用训练好的模型跑出来的）中抽取；c) 像我们一样，专门采样一些“次高不确定性”的构型作为测试集。这能更真实地反映模型的泛化能力。
2. 误差指标的局限性：不能只看能量和力的均方根误差（RMSE）。一个在测试集上RMSE很低的模型，可能在模拟中产生非物理的相变或结构。必须进行“物理验证”。对我们而言，就是对比不同模型在标准模型表面（如(100), (111), (211)）上跑出的水密度分布、吸附构型等，看其是否与高精度参考（如GRACE 2-Layer或有限的DFT-MD）一致。这是选择最终生产模型的黄金标准。
3. 主动学习收敛的判断：不确定性下降曲线会逐渐平缓。需要设定一个阈值（如平均力不确定性低于0.02 eV/Å），并结合检查新采样构型的结构类型是否开始重复。当新增数据不再带来新的结构类型时，即可认为收敛。

4.2 大规模模拟与分析的技巧

1. 表面粗糙度的量化：均方根粗糙度（Rq）是一个直观的全局指标，但在原子尺度，它可能掩盖局部特征。我们同时使用了局部曲率分析和配位数分布来辅助表征。例如，统计表面铜原子配位数小于6（体相铜为12）的比例，能直接反映表面“不饱和”位点的密度。
2. 无监督分类的参数敏感性：UMAP的降维结果受其参数（如n_neighbors,min_dist）影响。我们的经验是：n_neighbors不宜过小（避免过度碎片化），通常设置为数据量的平方根量级；min_dist可以设得较小（如0.1），以便在投影图上更好地区分密集的簇。关键是要多次尝试不同参数，观察聚类结构的稳定性，并结合化学直觉（如手动检查不同簇的代表性原子环境）进行验证。
3. 水吸附的判断标准：仅凭原子间距离（如Cu-O < 2.5 Å）判断化学吸附有时不可靠，因为热涨落会导致距离波动。更稳健的方法是结合配位数分析（积分RDF第一峰）和动力学稳定性：追踪一个水分子中的氧原子与特定铜原子的距离随时间的变化，如果在一段相当长的模拟时间内（如几十皮秒）保持近距离且振动频率高，则可判定为强化学吸附。

4.3 计算资源与工作流管理

1. 分层计算策略：DFT计算（用于主动学习）使用高精度设置，但只算小体系（<200原子）。MLIP训练在中等规模GPU集群上进行。最终的大规模MD模拟（20万原子）则在CPU集群上使用LAMMPS进行。合理分配计算资源至关重要。
2. 自动化工作流：主动学习涉及“模拟->分析不确定性->提取构型->DFT计算->重新训练”的循环。手动操作极易出错。我们使用Python脚本（结合ASE、pymatgen等库）将整个过程管道化，确保可重复性。数据库版本管理（例如用dvc）也非常重要。
3. 可视化与诊断：大量数据的分析依赖可视化。我们重度使用OVITO进行结构快照、配位数着色、轨迹动画的生成。同时，编写自定义脚本实时监控模拟过程中的能量、温度、压力以及特定结构序参量的变化，以便及时发现问题。

5. 未来展望与项目延伸

这项工作为模拟真实形貌的固液界面奠定了基础。基于此，有几个非常自然的延伸方向：

1. 引入反应物与电化学环境：当前工作只研究了纯水界面。下一步的核心是在体系中引入CO₂、CO或其他碳中间体，利用同样的无监督分类方法，监控碳物种的局部环境演化，从而在原子尺度上追踪反应路径。更进一步的挑战是模拟带电界面，即施加电极电位。这需要结合连续介质模型（如Poisson-Boltzmann）或更高级的恒电位方法，是当前领域的前沿。
2. 动态演化与时间尺度：目前的模拟是在静态（或热力学平衡）的粗糙表面上进行的。真实的电催化过程中，表面可能在反应条件下动态重构。可以尝试进行更长时间尺度的模拟，或结合蒙特卡洛方法，研究表面粗糙度、缺陷密度随电位、时间的演化。
3. 从结构到性能的桥梁：识别出潜在的活性位点（如簇1、2的角落原子）后，可以通过机器学习加速的动力学模拟或微动力学模型，定量估算不同位点的反应速率常数，最终将原子尺度结构与大尺度的催化活性（如电流密度、产物选择性）关联起来。
4. 方法论的通用化：我们发展的“主动学习+空间分辨不确定性+非晶矩阵嵌入”工作流，具有很强的通用性。它可以被应用于其他任何需要研究复杂、非周期性界面的体系，如合金催化剂表面、电极-电解质界面、生物分子-材料界面等。

构建一个可靠的MLIP并用于探索未知的复杂界面，就像训练一个探险家并为他绘制未知海域的地图。主动学习确保地图覆盖了所有危险和有趣的区域，无监督学习则帮助我们识别出地图上那些从未被标注过的、却可能埋藏着宝藏的特殊地貌。这次对粗糙铜-水界面的探索，不仅让我们看到了角落原子作为强吸附位点的清晰证据，更重要的是验证了这套方法论的强大能力——它让我们有信心去模拟和解析更多真实而复杂的材料界面。