静电筛选与机器学习势函数加速：高通量预测材料分裂空位缺陷

1. 项目概述：从“点”到“裂”的缺陷认知革命

在材料研究的微观世界里，点缺陷——尤其是空位缺陷——长久以来被我们视为一个简单的“缺失”。想象一下，一个完美的晶体格子，其中一个原子被移走了，留下一个空洞，这就是我们熟知的简单空位。这个模型简洁、直观，是理解掺杂、离子电导、甚至催化活性的基础。然而，就像平静的湖面下可能暗藏复杂的涡流，真实的材料世界远比这个简化模型要复杂。在某些特定的化学环境和晶体结构中，这个“空洞”并不安分，它会自发地发生重构，一个空位“分裂”成两个空位，并伴随一个间隙原子的产生，形成所谓的“分裂空位”复合体。这不仅仅是几何构型的变化，它深刻地改变了局部的电荷分布、应力场，进而对材料的电子结构、离子迁移势垒乃至宏观性能产生颠覆性影响。

传统上，发现这样的分裂空位依赖于研究者的物理直觉和大量的试错式第一性原理计算，效率低下且极易遗漏。我们这次的工作，正是要打破这个瓶颈。我们发展并验证了一套结合了静电学预筛选与机器学习势函数加速弛豫的高通量计算框架。其核心思路非常清晰：首先，利用基于氧化态和库仑相互作用的静电学模型，从海量的可能分裂构型中快速、廉价地筛选出能量可能较低的候选者；然后，对这些“种子选手”使用高精度的机器学习势函数进行快速的几何结构优化，最终用少数几次昂贵的密度泛函理论计算进行确认。这套方法就像一位经验丰富的侦探，先用快速的线索（静电学）缩小嫌疑犯范围，再用精准的技侦手段（机器学习弛豫）进行排查，最后才动用终极的DNA鉴定（DFT），从而在材料项目的庞大数据库中，系统性地搜寻那些隐藏的、能量更优的分裂空位缺陷。我们的目标，是为材料缺陷研究提供一套可扩展、自动化的“勘探”工具，将分裂空位从偶然的发现，转变为可预测、可设计的对象。

2. 核心原理：为什么空位会“分裂”？

要理解我们方法背后的逻辑，首先得深入探究分裂空位形成的物理根源。这绝非一个随机的过程，而是晶格能量最小化驱动下的必然结果，主要受两大因素支配：局域键合重构与长程静电相互作用。

2.1 晶格弛豫与局域键合重构

当一个原子被移除形成空位后，周围的原子不再处于受力平衡状态。它们会朝着空位中心移动，以降低系统的总能量，这个过程称为晶格弛豫。对于某些特定的元素和晶体结构，这种弛豫可能不是对称的、均匀的收缩。例如，当一个高价阳离子（如Ga³⁺、Sb⁵⁺）缺失时，留下的高正电荷空位会对邻近的阴离子（如O²⁻）产生极强的库仑排斥。为了缓解这种强烈的静电排斥，邻近的一个阴离子可能会被“推离”其原有的晶格位置，进入一个间隙位点，从而在原始空位两侧形成两个电荷密度较低的空位。这就完成了从V_X（一个X空位）到[V_X - X_i - V_X]（两个X空位夹一个间隙原子X_i）的转变。

这个过程可以类比为紧绷的弹簧网：当你剪断一根弹簧（移除原子），周围的弹簧会收缩。但如果某个连接点特别脆弱（静电排斥强），它可能会断裂并重新连接，形成两个新的、更稳定的节点（分裂空位），而不是简单地拉紧。在我们的工作中，通过分析数千个DFT弛豫后的结构，我们发现分裂空位构型中原子的最大位移，相对于其体相位置，呈现出与简单空位显著不同的分布模式（如原文图S1所示），这为从几何上识别分裂空位提供了直观依据。

2.2 静电相互作用的主导角色

静电相互作用是驱动这一分裂过程的远程力。一个带电的空位缺陷会在晶体中产生一个长程的静电势场。在离子性较强的材料中，如大多数金属氧化物和氮化物，这一效应尤为显著。分裂空位构型本质上是对原始点电荷缺陷产生的强静电场的响应。通过将缺陷电荷分布分散到两个空位和一个间隙原子上，系统可以更有效地屏蔽缺陷电荷，降低整体的静电能。

我们方法的第一步——静电学筛选——正是基于这一原理。我们为晶体中的每个原子分配了整数氧化态（利用pymatgen和doped中的算法），然后将候选的分裂空位构型（即预设两个空位和一个间隙原子的位置）建模为一系列点电荷。通过计算这些点电荷构型在介电介质中的静电能，并与简单点空位的静电能进行比较，我们可以快速估算出哪些分裂构型在静电意义上更有利。这是一个计算量极小的步骤，但能有效过滤掉绝大多数能量上不利的构型，将候选数量降低几个数量级。

2.3 机器学习势函数的桥梁作用

然而，静电模型是高度简化的，它忽略了键合的方向性、短程排斥以及电子结构的细节。因此，通过静电筛选的构型，必须经过更精确的原子级弛豫来验证。传统上，这完全依赖于DFT计算，成本高昂。我们引入机器学习势函数（本文中主要使用MACE-mp）作为“加速器”。

机器学习势函数通过在大量DFT数据上训练，能够以接近DFT的精度预测原子间的相互作用力和能量，但计算速度比DFT快成千上万倍。在我们的流程中，所有通过静电初筛的候选结构，都会使用MACE-mp势函数进行完整的几何优化。这一步可以快速淘汰那些在考虑原子细节后仍然不稳定的构型，只将最有希望的少数候选者送入最终的DFT弛豫和能量计算环节。这种“静电粗筛 + ML精炼 + DFT确认”的三级漏斗式流程，在保证结果可靠性的前提下，实现了计算效率的极大提升。

3. 方法实现：从理论到代码的完整工作流

我们的方法不是一个孤立的算法，而是一个集成在开源生态中的自动化工作流。下面我将拆解其中的关键步骤和实现细节。

3.1 分裂空位的几何识别算法

如何判断一个弛豫后的缺陷结构是简单空位、分裂空位还是其他复杂构型？我们采用了基于doped包中站点匹配算法的几何判据。该算法的逻辑清晰而严谨：

1. 输入：一个弛豫后的缺陷超胞结构，以及对应的完美晶体（体相）结构。
2. 匹配：尝试将缺陷结构中的每一个原子位置与体相结构中的原子位置进行匹配。允许存在一定的距离容差（通常设置为体相键长的50%）。这是一个寻找“谁是谁”的过程。
3. 分类：
   • 简单空位：如果体相结构中有且仅有1个站点在缺陷结构中找不到匹配项，并且缺陷结构中的所有原子都能在体相中找到匹配项。这对应V_X。
   • 分裂空位：如果体相结构中有2个站点在缺陷结构中找不到匹配项（对应两个空位V_X），并且缺陷结构中有1个站点在体相中找不到匹配项（对应一个间隙原子X_i）。这完美对应[V_X - X_i - V_X]模型。
   • 非平凡空位：所有其他不满足以上两种情况的结构，可能涉及更复杂的重构或多种缺陷的复合。

这个基于计数的几何规则，通过doped高效实现，成为了我们自动化分类的基石。在实际操作中，距离容差的选择需要谨慎。50%的键长是一个经验值，对于大多数氧化物和氮化物适用，但对于某些键长差异很大或弛豫特别剧烈的体系，可能需要微调。

3.2 静电筛选模型的构建���计算

静电筛选是我们高通量流程的“第一道滤网”。其核心是计算点电荷模型的静电能。

1. 氧化态赋值：这是静电计算的前提。我们利用pymatgen的BVAnalyzer（键价分析器）和ICSD数据库的先验知识，为晶体中的每个原子分配整数氧化态。对于约15万种材料项目数据库中的化合物，我们成功为其中约11万种确定了氧化态。
2. 缺陷电荷模型：对于带电量为q的缺陷，我们将其建模为位于缺陷位置的一组点电荷。对于分裂空位，这涉及到两个带+q/2（如果是阳离子空位）的空位点和一个带-q的间隙原子点（具体符号取决于缺陷类型和氧化态）。电荷的精确分配基于我们确定的氧化态。
3. 静电能计算：在连续介质近似下，计算这些点电荷在具有特定介电常数（各向同性或各向异性）的晶体环境中的相互作用能。我们通常使用Ewald求和或类似的方法来处理周期性边界条件。计算出的静电能E_elec是一个相对值，用于比较不同分裂构型之间的相对稳定性。

注意：静电模型是快速筛选工具，而非精确能量预测。它忽略了原子弛豫、电子局域化、键合等效应。因此，其绝对能量值意义不大，但其相对排序对于识别低能候选结构非常有效。我们的测试表明，静电能量排名靠前的构型，有很大概率在DFT弛豫后仍然是低能态。

3.3 机器学习势函数加速弛豫的实战细节

我们选择MACE-mp作为机器学习势函数，因为它在大规模材料数据库上进行了预训练，对氧化物和氮化物体系具有较好的泛化能力。

1. 模型选择与精度权衡：MACE-mp提供了不同大小的模型（small, medium, large）。我们的测试（如原文图S7-S10）表明，small和large模型在预测精度上表现相近，而medium模型稍差。鉴于计算速度的显著差异（small模型比large模型快数倍），我们在生产性筛选中统一使用small模型和32位浮点精度。这能在保证能量和力预测平均绝对偏差约1-2 meV/atom的前提下，最大化计算吞吐量。
2. 优化器选择：我们对比了ASE（原子模拟环境）中多种几何优化算法，包括FIRE、BFGS和GOQN等。实测发现，GOQN（Good Old Quasi-Newton）算法在稳定性和收敛速度上取得了最佳平衡。像GPMin（高斯过程最小化器）这类算法，虽然在某些问题上高效，但对于我们涉及上百个原子的大超胞，其内存需求过高，容易导致计算崩溃。
3. 工作流集成：我们使用doped的流程来生成缺陷超胞，并自动准备静电筛选。通过自定义脚本，将筛选出的候选结构传递给ASE，调用MACE-mp势和GOQN优化器进行弛豫。弛豫的收敛标准通常设置为力小于0.01 eV/Å。整个过程可以通过Python脚本进行批量化管理，实现无人值守的自动化运行。

3.4 有限尺寸修正的考量

在超胞方法中，计算带电缺陷的形成能必须考虑有限尺寸修正E_corr，以消除周期性镜像电荷之间的虚假相互作用。对于分裂空位这类多电荷中心的缺陷，修正项的处理比简单点缺陷更复杂。

1. 修正项的影响：E_corr与电荷平方q²成正比，与介电常数ε成反比，与超胞尺寸L近似成反比。对于相同电荷态q的分裂空位和点空位，虽然q相同，但由于电荷分布不同（分裂空位的电荷更分散），它们的E_corr也可能不同。
2. 我们的处理方式：在关键的对比计算中（如表1中的测试集），我们采用了Kumagai & Oba发展的eFNV修正方案。该方案能原生处理各向异性介电屏蔽和不同方向静电势平移的平均。对于分裂空位，我们将其质心作为缺陷位置用于修正计算。
3. 影响评估：我们的分析表明，对于大多数体系，点空位与基态分裂空位之间的有限尺寸修正能差异在10-65 meV之间，约占修正总量的5%，远小于两者之间的超胞能量差（通常>0.5 eV）。然而，对于某些亚稳态分裂空位，特别是当V_X-X_i距离较大、构型对称性较低时，修正能差异可能达到~0.1-0.3 eV。这意味着，在能量差异很小的临界情况下，必须谨慎评估有限尺寸修正的影响，有时甚至需要在更大的超胞中重新计算以得到可靠结论。

4. 大规模筛选结果与材料学启示

我们将这套方法应用于Kumagai等人的金属氧化物数据集（~600种阳离子空位）和材料项目数据库的扩展筛选，获得了系统性的发现。

4.1 分裂空位的普遍性与元素偏好

我们的ML加速筛选在材料项目数据库中预测了大量可能存在低能分裂空位的化合物。通过统计分析（如原文图S11的热图），我们发现了清晰的元素化学趋势：

• 阳离子特性：高氧化态、中等离子半径的阳离子更容易形成分裂空位。例如，Sb⁵⁺, V⁵⁺, W⁶⁺, Mo⁶⁺, As⁵⁺等元素出现的频率和预测的能量降低幅度都显著偏高。这是因为高电荷产生强局域电场，驱动晶格重构；而离子半径适中，则为原子迁移提供了可能的空间。
• 晶体结构环境：非中心对称、具有柔性多面体连接（如层状、链状结构）的晶体，比分立四面体或八面体紧密堆积的结构更有利于分裂。例如，在C2/c结构的Sb₂O₅和Ga₂O₃中，我们都观察到了显著的分裂空位稳定化现象。
• 宿主化合物化学：含有高电负性阴离子（如O²⁻）和可极化阳离子的化合物是分裂空位的“温床”。静电驱动是核心，因此离子性越强，效应往往越明显。

4.2 亚稳态分裂空位的多样性

除了能量最低的基态分裂构型，静电筛选还揭示了丰富的亚稳态分裂空位（能量比点空位低，但比基态分裂空位高）。在金属氧化物数据集中，我们发现了超过200个能量在点空位0.5 eV以内的独特亚稳态（原文表S1）。

这些亚稳态具有重要的物理意义：

• 构型空间复杂：它们对应着不同的V_X-X_i距离和相对取向，代表了缺陷在势能面上的不同局部极小点。
• 影响动力学过程：在离子迁移（如扩散）或缺陷反应过程中，这些亚稳态可能作为中间态或过渡态，从而显著影响材料的动力学性质。例如，一个迁移的离子可能会被这些亚稳态缺陷捕获一段时间。
• 对计算的要求：发现这些亚稳态需要充分的构型采样。我们的静电预筛选结合ML弛豫，能够以较低成本探索广阔的构型空间，这是传统手动或随机采样难以做到的。

4.3 与纯DFT结果的对比验证

为了验证我们混合方法的可靠性，我们在已知的测试集上进行了严格的基准测试。

1. 静电筛选的有效性：我们比较了仅基于静电能量排序的候选构型，与经过DFT弛豫后的最终能量排序。结果显示，静电能量最低的10-20个候选构型中，有很高概率（>70%）包含DFT弛豫后的最低能量分裂空位。这说明静电筛选作为“粗筛”是极其有效的。
2. ML弛豫的精度与效率：对比MACE-mp弛豫后的结构与DFT弛豫后的最终结构，我们发现：
   • 几何结构：原子位置的均方根偏差通常小于0.1 Å，对于筛��目的而言完全足够。
   • 能量排序：MACE-mp预测的相对能量顺序与DFT结果高度一致。虽然绝对能量可能存在数十meV的偏差，但用于识别“最低能量候选者”这个任务，其可靠性非常高。
   • 速度提升：MACE-mp（small模型）的单次弛豫比DFT（使用中等精度泛函）快约3个数量级。这使得对成千上万个候选结构进行弛豫成为可能。
3. 最终DFT确认：所有通过ML弛豫后看起来有希望的构型，我们都用更高精度的DFT计算（如PBEsol或HSE06）进行了最终的能量确认和电子结构分析。这确保了最终报道结果的量子力学精度。

5. 实操指南、避坑经验与扩展思考

基于大量的实际计算经验，我总结出以下关键的操作要点和常见问题解决方案，这些是在官方文档中不易找到的“实战心得”。

5.1 工作流搭建的实用步骤

1. 环境准备：建议使用Conda创建一个独立环境，安装pymatgen,doped,ase，以及MACE-mp的接口包（如mace-torch）。确保所有包的版本兼容。
2. 输入文件准备：你需要完美晶体的POSCAR文件（或CIF文件）。使用pymatgen的Structure类读入，并利用doped的DefectsGenerator生成所有对称性不等价的点空位超胞。这一步会生成初始的缺陷结构。
3. 静电预筛选脚本：

   # 伪代码逻辑 from doped.analysis import get_oxidation_states from local_electrostatic_module import calculate_split_vacancy_energy # 假设的自定义模块 structure = get_structure_from_file("POSCAR") oxidation_states = get_oxidation_states(structure) # 获取氧化态 # 为某个空位位点生成候选分裂构型 vacancy_site = get_vacancy_site(structure) candidate_split_configs = generate_split_configs(vacancy_site, structure) electrostatic_energies = [] for config in candidate_split_configs: # 将构型建模为点电荷集合 point_charges = model_as_point_charges(config, oxidation_states) # 计算静电能 (需实现或调用现有库，如考虑各向异性介电张量) e_elec = calculate_electrostatic_energy(point_charges, dielectric_tensor) electrostatic_energies.append((config, e_elec)) # 按静电能排序，选取前N个低能候选 low_energy_candidates = sorted(electrostatic_energies, key=lambda x: x[1])[:N]

4. ML弛豫批量提交：将上一步得到的候选结构写成单独的POSCAR文件。编写一个批量脚本，用ASE读取每个POSCAR，设置MACE-mp势和GOQN优化器，进行弛豫并输出最终结构和能量。
5. 后处理与分析：使用doped的站点匹配功能对弛豫后结构进行分类（简单/分裂/非平凡）。提取能量，与点空位能量比较。可视化低能分裂构型，分析键长、键角变化。

5.2 关键参数选择与常见陷阱

• 静电筛选的电荷模型：最简单的模型是将空位视为带完整离子电荷的点电荷。但对于共价性较强的材料，这可能高估了静电效应。一个更稳健的做法是使用Bader电荷或Mulliken布居分析（从一次DFT计算中获得）来分配缺陷的有效电荷，但这会增加计算量。对于高通量筛选，从整数氧化态开始通常是可行的。
• ML势的适用性检查：MACE-mp是在宽泛的材料数据集上训练的，但对于非常特殊的局部化学环境（如极端畸变、罕见氧化态），其外推能力可能有限。强烈建议对你所研究的材料体系，随机选取几个点空位和分裂空位构型，同时进行MACE-mp和DFT弛豫，对比最终结构和能量差异。如果偏差系统性较大（>50 meV/缺陷），可能需要考虑使用专门针对该体系训练的ML势。
• 超胞尺寸敏感性：分裂空位涉及多个缺陷中心，其相互作用范围可能比点缺陷更大。我们使用的超胞（平均~100-200原子，等效立方长度~10-14 Å）对于大多数体系是足够的，但对于介电常数很小（ε < 10）或缺陷电荷很大（|q| > 2）的体系，需要格外小心。务必检查有限尺寸修正的幅度，并考虑使用更大超胞进行验证计算。
• 构型采样充分性：静电筛选依赖于预设的分裂构型。我们通常基于晶体学对称性，在空位周围一定半径内枚举可能的间隙位点和另一个空位位点。这个搜索半径至关重要。半径太小会错过稳定构型，太大则计算量激增。一个经验法则是搜索到第二或第三近邻原子壳层。对于各向异性强的晶体，可能需要沿不同晶向设置不同的搜索半径。

5.3 结果解读与物理意义挖掘

当你发现一个能量显著降低的分裂空位后，真正的材料物理分析才刚刚开始：

1. 电子结构分析：计算分裂空位和点空位的态密度（DOS）。你通常会观察到缺陷态在禁带中位置和宽度的变化。分裂空位可能引入更局域或更离域的缺陷态，从而影响光吸收、载流子捕获等性质。
2. 迁移势垒计算：分裂空位可能作为离子迁移的新通道。使用爬坡弹性带（NEB）方法，计算间隙原子X_i在两个空位V_X之间跳跃的势垒，并与体相中的本征迁移势垒比较。这可以预测分裂空位对离子电导率是促进还是抑制。
3. 与实验关联：分裂空位的稳定存在可能解释一些实验观测。例如：
   • 非对称的局域结构：在扩展X射线吸收精细结构（EXAFS）谱中，分裂空位可能导致配位壳层距离分布出现双峰或其他非对称特征。
   • 异常高的缺陷浓度：如果分裂空位的形成能低于点空位，则在热平衡条件下，其浓度会更高，这或许能解释某些材料中异常高的空位浓度测量值。
   • 退火行为：分裂空位可能在特定温度下合并或分解，这可能在差示扫描量热法（DSC）或电阻率-温度曲线上留下特征信号。

5.4 方法局限性与未来展望

没有任何方法是万能的，我们的框架也有其边界：

• 强电子关联与局域化：对于强关联电子体系（如某些过渡金属氧化物），静电模型和标准ML势的预测能力会下降。可能需要结合DFT+U或更高级的电子结构方法。
• 动力学稳定性：我们的筛选基于静态能量。一个在0K下能量低的分裂空位，在有限温度下可能由于熵效应或动力学不稳定性而变得不重要。需要进行声子谱计算或分子动力学模拟来评估其动力学稳定性。
• 带电态与费米能级：我们的分析通常针对特定的电荷态。在实际材料中，缺陷的稳定电荷态随费米能级变化。完整的缺陷形成能计算需要考虑不同电荷态，并绘制形成能随费米能级变化的图。
• 与更复杂缺陷的竞争：我们主要关注孤立的空位。在实际材料中，空位可能与杂质、其他本征缺陷（如间隙原子、反位缺陷）结合形成复合体，这些复合体可能比分裂空位更稳定。

未来的发展方向是清晰的：将这套框架与更先进的主动学习结合，让ML势在筛选中自我改进；扩展到更复杂的缺陷类型（如双空位、缺陷团簇）；以及开发直接与实验表征数据（如STEM图像、光谱）进行对比验证的工具链。分裂空位的研究，正在从一个需要运气和直觉的领域，转变为一个可计算、可预测、可设计的定量科学分支。