在晶体固体中寻找低周期子结构：一种新方法

《ACS Applied Materials & Interfaces》：Search for Low-Periodic Substructures in Crystalline Solids: A Novel Approach

　　高分辨率图像
下载MS PowerPoint幻灯片

新2D材料的发现对于推进电子技术和量子技术的发展至关重要。由于大多数2D材料都来源于层状体结构，因此识别可剥离的晶体并估计分离单层所需的能量是关键步骤。为了解决这个问题，我们开发了一种基于Voronoi划分构建晶体图、通过键价理论估计相互作用强度以及在追踪周期性变化的同时迭代去除弱链接的稳健且计算成本较低的方法。我们通过文献和从头算（ab initio）结果验证了该方法，证明它可以可靠地识别层并提供适合作为筛选参数的层间键合能的近似估计值。随后，我们将该方法应用于分析48,504个预选的实验晶体结构，发现了694种以前未报道的2D材料，这些材料属于530种不同的结构原型。最后，我们使用从头算模拟提供了分离层的结构和电子特性的概述。

**引言**
低维（1D和2D）或更准确地称为低周期性（1-周期性和2-周期性）材料是一些最有前途的候选材料，适用于超越硅的电子、光电子和量子计算应用。近年来，它们被广泛认可的重要性激发了科学界发现和表征具有潜在突破性特性的新型低周期性材料的热情。在各种发现策略中，一种特别有效的方法是对3-周期性（体相）晶体结构进行系统分析，以识别可以分离并进一步利用的嵌入式低周期性子结构。(1?4)

当考虑像石墨或过渡金属硫族化合物这样的已知层状材料时，即使通过直接观察体结构也可以很容易地识别出相应的层状子结构。然而，在更大范围内进行系统发现需要更严格和自动化的方法。为了应对这一挑战，研究人员开发了许多自动子结构识别方法。这些方法根据其基本原理大致分为几类。第一类方法依赖于计算体晶体的物理性质，以指示低周期性子结构的存在。这些方法包括：(i) 分析力常数(5)，反映原子间相互作用的强度；在选定的力常数阈值下获得的子结构的周期性决定了晶体的周期性；(ii) 杨氏模量分析(6)，它能够识别经常反映子结构周期性的弹性性质的各向异性。第二类方法包括技术，如(iii) 实验结构和理论优化结构的直接比较(7)，这可以通过结构松弛模式的特性揭示子结构之间的弱相互作用。虽然这些方法在识别低周期性子结构方面取得了显著成功，但它们存在一些实际限制。特别是，它们严重依赖于密度泛函理论（DFT）计算，因此需要大量的计算资源，这限制了它们在大规模结构数据库中的广泛应用。

第三类重要且广泛使用的方法是基于简单原子间距离的算法，有几种不同的实现方式。(1,8,9) 这些算法的基本概念相当直接：算法在晶体结构中搜索强键网络，并计算识别出的基元的周期性。如果识别出的强键网络的周期性低于3-p，则可以认为晶体包含相应的层(2-p)或链(1-p)子结构。在实践中，基于原子间距离的算法使用范德华（RvdW）或共价（Rcov）半径系统来定性地估计键强度。在这个框架内，如果原子间距离rij小于相关半径R与一个小容差Δ之和，则认为存在键；否则认为不存在键，从而可以构建强键网络。然而，这些方法也有一些缺点，包括需要适当选择参考半径系统和可调的Δ参数，这些参数会显著影响结果。此外，这些方法不能直接估计分离识别出的层所需的能量，因此通常需要通过DFT计算来补充。因此，目前明显缺乏一种低成本、广泛适用的方法，能够提供可剥离性的定量估计。

在晶体化学领域，键价（BV）理论长期以来一直是一种有价值和成熟的方法，用于晶体结构分析。(10) 如今，BV理论的应用包括验证结构解、估计氧化态、评估潜在的晶格应变，以及分析和预测晶体结构内的离子迁移路径。(11) 这一理论框架提供了几个优势，使其特别适合我们的目的。其核心是，键价理论允许估算键级，即反映形成键的原子对之间共享的电子对数量的量。如方程1所示，键级或键价（vij）的计算需要的输入数据非常少：只需原子间距离rij和表列参数R0和b，后者假定对于不同的化合物中的同一原子对ij是通用的。

??????=exp[(??????0???????)/??????]

BV方法对于我们的工作的一个特别相关的优势是众所周知的键级与许多化学键的键强度之间的相关性，其中较大的键级意味着更强的键合。最近的大量研究(12)展示了键级与电子密度拓扑的各种描述符之间的明确关系，包括键临界点的电子密度ρc、电子密度的拉普拉斯算子，以及与键伸缩频率平方之间的强相关性。对于大多数原子对，ρc(ij)和vij之间的关系几乎是线性的，而这些线的斜率似乎与相互作用类型一致。这些关系为键价提供了成为许多类型化学键的键强度的定量估计器的理论基础。因此，键级使得能够创建一组新的描述符，可以系统地表征低周期性子结构内部和之间的键合强度分布。

键合强度各向异性描述符已在几项研究中进行了探索，并通过成功应用于分子晶体结构和机械性质的表征而证明了它们的实用性。例如，已经提出了对分子间相互作用维度的全面分析，并将其应用于硫族二氰化物系列的比较研究(13)、对氢键连接的胍结构的详细描述(14)，以及2-溴马拉醛晶体中键合模式的阐明(15)。这种计算方法使得可以将总晶格能量分解，并将晶格能量分数分配给晶体的低周期性片段（层、链、二聚体等），从而从能量角度提供了对晶体结构主要结构的宝贵见解。另一个交互各向异性描述符的示例是X参数，它是专门为识别分子晶体中的解理面而开发的。(16)

在这项工作中，我们开发了一种计算工具，解决了当前子结构搜索算法的局限性，并将其应用于扩展可剥离晶体的范围。具体来说，我们创建了一个开源的CrystalSubstructureSearcher代码，该代码能够通过利用键价理论来高效识别体晶体结构中的低周期性子结构。这种方法结合了基于简单原子间距离的方法的计算效率与键价理论的定量描述能力。它能够快速筛选大量晶体结构，并分析超出基于DFT方法范围的系统。此外，它还提供了有价值的晶体化学数据，包括子结构电荷、子结构鲁棒性的代理指标以及表征键合强度各向异性的描述符，这些数据对于后续过滤非常有用。我们将我们的方法与文献结果、其他低周期性材料数据库以及专门的从头算模拟进行了基准测试，然后使用它研究了来自无机晶体结构数据库（ICSD）(17)和Pierson’s Crystal Data (PCD) (18)数据库的48,504个实验结构。在这些结构中，我们发现了694种以前未报道的、坚固且易于剥离的材料，后来我们使用DFT模拟对其进行了优化和表征，排除了含有稀土的化合物。

**结果和讨论**
**基础算法**
CrystalSubstructureSearcher是一种用于分析晶体结构并识别其中片段或低周期性子结构的代码。接收到包含晶体结构的cif文件作为输入后，CrystalSubstructureSearcher通过构建晶体结构的图表示来启动分析。晶体结构中的原子间连接性是使用pymatgen中的VoronoiNN类实现的Voronoi算法计算的。(19) 构建的图称为结构图（SG），它由代表原子的节点和代表原子间键的边组成。计算了一组原子间键的特征，然后将它们作为属性分配给每条边。选定的原子间接触特征集旨在作为键强度的代理指标。它包括原子间距离（R）；对应于给定原子间接触的Voronoi–Dirichlet多面体的面积（A）和solid角（SA）；(20) 根据参考文献(21)计算的穿透指数（PI）；以及使用键价参数R0计算的键的键价（BV）。前98种化学元素的4851个原子对的所有R0值都是通过基于874个实验R0值的ML回归模型估计的，如方法部分所述。对于A、SA、PI和BV描述符，较大的值对应较强的键，而对于R，相反，较小的值对应原子间的强相互作用。Solid角已经被用作分子堆积的多级拓扑描述的键强度标准。(22) 用户从集合中选择一个属性作为结构图中的边权重。

实现代码的主要任务是提取低周期性子结构。通过在SG中迭代移除权重小于阈值的边集来识别SG中的子结构。每次迭代的阈值逐渐增加，如果使用A、SA、PI和BV作为权重，则从最小值开始；如果使用R作为边权重，则从最大的原子间距离值开始破坏边。在每次迭代中，破坏了一组特定的键后，会检查编辑后的结构图的维度。如果在迭代过程中SG组分的最大周期性发生了变化，则在该步骤获得的图将被保存下来，以便进一步处理。这样的SG快照对应于具有降低周期性的识别出的子结构，它可以由几个不同的SG组分组成。例如，可能会出现两个分别对应于单元格中两个不同但可能等效的层的图。之后还需要进行一次额外的SG编辑步骤，因为在之前的迭代中，可能会破坏组分内部的接触（组分内接触）以及组分间的接触。因此，需要进行一个辅助的键集恢复程序。恢复的键集是那些恢复后不会导致图组分最大周期性变化的键集。总体而言，所描述的算法属于贪婪算法类，它在每一步都通过做出局部最优选择来解决其任务。在我们的案例中，每次迭代我们都会破坏最弱的（由选定的边权重反映的）相互作用。虽然它不能保证产生全局最优解，但在子图搜索的情况下，这是唯一可行的方法。表示算法的整体流程图显示在支持信息中的方案S1中。图1示意性地表示了对Ca2Sb结构的分析。

**图1**
图1. Ca2Sb结构（ICSD 154）分析的第一次迭代，导致识别出2-周期性子结构。单元格沿a轴显示。在此示例中，键价被选为边权重。每次迭代中断裂的SG的边用红色突出显示，它们的特性在下面的表中显示。最后一步恢复的边用绿色显示。BVS总和是SG组件中保留边的键合价之和。高分辨率图像下载MS PowerPoint幻灯片
当编辑后的SG中子结构的最大周期性为零时，算法终止。换句话说，我们在所有剩余的边都对应于0周期性分子片段中的内部键时停止。因此，整个过程产生了一个给定晶体结构的所有低周期性子结构的列表。接下来，CrystalSubstructureSearcher为每个子结构以及初始结构图计算一组键合价总和（BVS）描述符，这些描述符旨在表征晶体结构及所有识别出的子结构中的键强度分布。计算出的BVS描述符的完整列表以及其他关于子结构及其相互之间关系的潜在有用数据显示在表S1中。
在计算出包含每个识别出的子结构的总BVS数据的BVS_x_periodicity描述符后，我们可以为晶体固体生成另一个重要的定量描述符，该描述符反映了键合强度的各向异性。我们将其称为内在固体周期性。它用一个连续的值来量化晶体结构在多大程度上是3周期性的、2周期性的或1周期性的，这是基于不同周期性的子结构之间的键合强度，与传统的离散周期性分类相反。对于Ca2Sb晶体结构，只识别出2周期性的中间子结构；因此，BVS_x_periodicity描述符中收集的子结构的总键合价总和如下：3周期性：36.556 BV单位，2周期性：30.069 BV单位，0周期性：2.412 BV单位。从这些数据中，我们可以评估将层状Ca2Sb子结构连接成块状Ca2Sb固体的接触强度为36.556–30.069 = 6.487 BV单位。类似地，将Ca原子和CaSb哑铃状结构组合成Ca2Sb层的相互作用强度为30.069–2.412 = 27.657 BV单位。通过对构成给定周期性子结构的接触作用的BVS进行归一化，我们得到了所需的内在固体周期性描述符。对于Ca2Sb晶体结构，最大的内在周期性出现在2周期性子结构中。因此，大部分键合发生在形成层的构建单元（Ca原子和CaSb分子片段）之间（2-p：0.81），而不是在形成块状晶体的层之间（3-p：0.19）。
除了子结构识别和BVS描述符计算外，CrystalSubstructureSearcher还能估计子结构的电荷。这一过程涉及计算SG组件之间接触点的电负性差异和键合价，从而提供了对识别出的子结构上电荷分布的宝贵估计。总之，CrystalSubstructureSearcher采用逐步算法的方法来分析晶体结构，包括低周期性子结构识别、BVS描述符和内在周期性的计算、界面特性以及子结构电荷估计。
为了对我们的方法进行基准测试和校准，我们首先在一个包含1713个理论优化的块状晶体结构的数据集上进行了2周期性子结构的搜索，这些结构最初用于创建MC2D数据库的一部分。(1,2) 该数据集包括之前被分类为具有层状结构的结构，这些结构要么容易剥离（784个结构，层间键合能不超过30 meV/?2），要么可能剥离（929个结构，键合能在30至120 meV/?2之间）。CrystalSubstructureSearcher中实现的算法仅在容易剥离子集中的两个结构和可能剥离子集中的26个结构中未能识别出2周期性子结构（见Supporting Information.xlsx文件中的no_layers_identified表格）。此外，只有对于52个结构（见Supporting Information.xlsx文件中的different_layers_identified表格），新算法识别的层与之前使用基于距离的算法识别的层不同。(1) 因此，新算法产生了与基于距离的算法一致的有意义的结果，因为在1713个结构中的1633个（95.3%）结构中，子结构是相同的。在容易剥离子集的结构中识别出的大多数2周期性子结构（775个结构）是中性电荷的，这与层间弱范德华相互作用的固体预期一致。仅在两个结构中发现了带电层，但其估计的电荷可以忽略不计，不超过0.26|qe|。相反，在可能剥离子集的30个结构（858个结构）中发现了带电的低周期性子结构。
可视化内在周期性的最佳方法是使用三元图，因为所有内在周期性值的总和为1.0。图2a显示了容易剥离子集中的晶体结构在内在周期性坐标中的位置。正如预期的那样，只有少数点接近表示内在3周期性晶体的三角形角落，其中亚硝酸汞(I)最接近内在3周期性固体（3-p：0.67，2-p：0.33）。总共，这个子集中的65%的结构具有内在2周期性特征2-p ≥ 0.5，即这些确实是内在2周期性晶体。除了2-p角落外，结构也集中在1-p/2-p的中点附近，以及在较小程度上靠近1-p角落。这些结构对应于由逐渐增大的棒状结构形成的层状结构，其中棒内的原子间键合比棒间键合更多。例如SiP结构（3-p：0.007，2-p：0.049，1-p：0.945），其1-p内在周期性最大。与图2b中绘制的可能剥离子集的比较表明，点的分布从1-p/2-p边缘向三角形中心偏移，这对于具有较弱层间键合能的结构来说是合乎预期的。总体而言，这种键合各向异性的描述符是一个有价值的工具，用于区分稳定的层状结构和容易分解成1周期性片段的棒状层状结构。
图2
图2. 用于创建MC2D数据库的数据集中晶体结构的内在周期性的表示。(a) 三元图显示了来自容易剥离子集（Ebinding小于30 meV/?2）的782个结构在内在周期性坐标中的位置。(b) 同一张图显示了来自可能剥离子集（Ebinding在30至120 meV/?2范围内）的903个结构。点的颜色根据inter_bvs_per_unit_area描述符的值进行区分。
图2a和图2b之间的整体相关性适中（Spearman等级相关性ρ = 0.83，Pearson’s r = 0.58），这已经相当好，特别是考虑到散点图（图S1b）中显示了401种不同接触类型和754种不同接触类型组合的界面键。（图S1b）然而，如果我们“放大”并考虑特定键类型之间的Ebind/A和BVS/A的相关性，我们可以发现更强的相关性。例如，在图S2中，散点图显示了同周期原子形成的同原子类 pnictogen (Pn···Pn)、chalcogen (Ch···Ch) 和 halogen (Hal···Hal) 接触具有更高的相关系数（ρ = 0.69 ÷ 0.96）。对具有丰富层间接触类型结构（至少5个结构）的分析表明，原则上不同的阈值BVS/A值应该对应于不同的接触类型，因为Ebind/A与BVS/A的依赖性不同（表S2）。不幸的是，只有在615个结构中遇到的26种接触类型（37.7%）的OLS拟合得到了统计上显著的（α = 0.05）斜率系数。大多数层间接触类型及其组合（图S1b）仅在1-3个结构中出现，这显然使得从OLS拟合中获得可靠的依赖性变得不可能。因此，我们需要设计一种方法来扩展算法的适用性，以适应具有多种界面接触类型的大型结构集合。我们提出了一种简化方法，决定通过将具有相似属性的原子分组来“聚合”这些原子。例如，碱金属和碱土金属被共同称为“电正性金属”（EPMs）。我们根据它们的相似性设立了九个粗略的原子组，如表S3所示。然后，我们将形成界面接触的原子简单地映射到这些聚合的原子组中，从而得到了新的聚合界面接触类型。例如，Cl···Cl和Br···Br接触都被映射到NM···NM（非金属）接触组，而Cl···Cl|Cl···Sr和Br···Br|Br···Ca则被映射到相应的EPM..NM|NM..NM接触组。总共获得了263种不同的聚合接触组组合，与最初的754种接触类型相比减少了近3倍。之后，我们使用至少在7个代表性结构中出现的接触组，利用statsmodel库中实现的Robust Linear Model (RLM)来建模Ebind/A与BVS/A的线性关系。（23）使用Huber损失函数的RLM模型进行拟合，因为它在接触分组后存在大异常值时能提供更稳定的回归系数估计。结果，对于覆盖1119个结构（68.5%）的29种聚合接触类型，我们获得了具有统计上显著回归系数的线性拟合（图S3），这与简单的接触类型相比几乎提高了2倍。最后，使用获得的回归系数，我们计算了每种聚合接触类型对应的BVS/A阈值，该阈值对应于30 meV/?2的键合能（见Supporting Information中的表S4）。此外，我们将BVS/A阈值提高了20%，以便接受更多接近任意选定的30 meV/?2边界线的结构。虽然这种分组策略可能会模糊同一组成员之间的某些化学差异，但我们估计相关的误差相对较小，如Supporting Information中所讨论的。
为了验证获得的Ebind/A与BVS/A依赖性，我们使用了一组新的具有层状结构的晶体结构，这些结构的晶胞足够小，可以直接用DFT计算处理（最多12个原子），并且不在用于创建MC2D数据库的结构中（使用的结构来自ICSD版本2018）。测试集包含162个计算出的键合能低于120 meV/?2的结构。键合能的计算使用了创建MC2D数据库时相同的程序，如参考文献(1和2)中所述。在这些晶体结构中，只有95个包含属于29种聚合接触类型之一的层间接触类型，从而可以使用BVS/A描述符来估计Ebinding（见Supporting Information.xlsx文件中的Ebinding_calc_vs_estimated表格）。图3中的散点图表明，对于聚合接触类型，简单的线性关系提供了足够准确的键合能估计，适用于筛选目的。值得注意的是，回归模型倾向于高估Ebinding，从而通过最小化假阳性来确保保守的筛选：难以剥离的材料将被排除，尽管一些可剥离的结构可能会被排除。平均绝对偏差仅为约10 meV/?2，证实了这种简单直接方法的有效性。
图3
图3. 使用DFT方法计算的Ebinding与通过线性Ebind/A与BVS/A关系估计的95个结构中的键合能的比较。
为了验证我们的方法是否可以直接在实验数据上可靠地使用，而无需中间且可能昂贵的优化步骤，我们还使用CrystalSubstructureSearcher分析了直接从晶体学数据库中获取的用于创建MC2D数据库的实验晶体结构，这些结构没有对块状单元几何进行任何DFT优化。Supporting Information中的图S4显示了DFT优化结构相对于实验晶体结构获得的BVS/A的相对百分比差异。总体而言，DFT优化导致BVS/A描述符增加，这在72%的结构中观察到，这是由于层间接触缩短所致。在56.1%的结构中，相对于原始结构优化的BVS每单位面积描述符的绝对百分比差异不超过25%。ICSD过滤和在大量结构集上的程序应用

在开发出算法、将其实现到CrystalSubstructureSearcher程序中，并在具有已知层间结合能量的数据集上验证了其性能后，我们继续将其应用于更大的结构集，以展示其识别具有新型低周期子结构的晶体结构的能力。来自ICSD（17）和PCD（18）数据库的晶体结构数据经过了一系列过滤步骤，旨在确保只选择最高质量的结构。最终，使用CrystalSubstructureSearcher代码分析了48,504种无机化合物。在对这些数据集的分析中，有1390个结构（占2.9%）被中止，主要原因包括原子间距离过短、Voronoi多面体构建过程中的错误以及其他结构异常和特殊性。在成功完成分析的47,114个晶体结构中，有27,350个结构被识别出2周期子结构（见图4a）。图4b中的三元图显示了这些结构的固有周期性分布。同时存在1周期和2周期子结构的结构集中在三角形的中心附近；而仅包含1周期子结构的结构则聚集在三角形的3周期顶点附近，这意味着它们大多数具有更强的3周期特性。相反，仅包含2周期子结构的结构集中在2周期/3周期三角形边的2周期顶点附近。总体来看，之前讨论的易于剥离的结构子集位于三元图的稀疏区域（参见图2a）。

图4
(a) 文氏图显示了含有低周期子结构的晶体结构数量。发现所有周期子结构的结构占分析结构的34%，仅含2周期子结构的结构占25%，仅含1周期子结构的结构占20%，其余21%的结构没有低周期子结构。
(b) 三元图显示了从ICSD中筛选出的结构在固有周期性坐标上的分布。注意，9468个没有低周期子结构的结构位于三角形的3周期顶点。由于结构特性，有955个结构的固有周期性无法计算。

通过BVS/A描述符为29种聚合接触类型确定的阈值，我们仅评估了8343个结构的可剥离性（在27,350个结构中），其中2135个结构的单位面积BVS低于30 meV/?2的阈值，因此我们预计这些结构含有弱结合的层状子结构。此外，我们还过滤掉了含有带电子结构的结构，最终得到了一组1989个潜在易于剥离的晶体结构，这些结构的电荷估计值<0.25|qe|（这个值接近之前在MC2D易于剥离子集中识别的层的最大绝对电荷）。最后，我们将识别出的子结构列表与现有数据库进行了比较，以评估我们结果的新颖性（以下称为CSS结果）。为了比较，我们使用了截至2025年6月的最新版本的MC2D（1,2）（https://mc2d.materialscloud.org）、2DMatpedia（3）（http://www.2dmatpedia.org）和C2DB（4）（https://c2db.fysik.dtu.dk）数据库，这些数据库包含了基于距离的算法在ICSD、COD和Materials Project数据库中的大量实验和理论晶体结构上的应用结果。重要的是，除了通过自上而下的方法理论剥离的层外，2DMatpedia和C2DB数据库还包含了通过自下而上的方法创建的2周期子结构，该方法系统地对使用自上而下方法识别的层进行了元素替换。为了避免因不同数据库中使用的结构优化程序差异而导致的错误识别新结构的情况，我们选择了更保守的基于组成的比较方法。我们比较了这三个数据库中通过自上而下方法识别的2周期子结构的组成与CrystalSubstructureSearcher识别的组成。如图5a中的文氏图所示，CrystalSubstructureSearcher识别出了929种独特的层组成（存在于963个晶体结构中），而这些组成在其它2D材料数据库中并不存在。还需要再进行一次过滤，以得到最终的可靠层列表。正如参考文献（6）中提到的，“……如果层状晶体的面内强度不足以抵抗层间作用力，那么在机械剥离过程中该材料的层面将会破裂，因此这种层状晶体就不是可剥离的晶体”。为了解决这个问题并从集合中剔除脆弱层，我们设计了一个额外的基于BV的描述符，称为max_intracomponent_bond_strength。这个描述符反映了构成2周期子结构的1周期或0周期子结构之间的键合强度。图5b给出了该描述符应用的一些例子。从W5I11晶体结构中提取的层由通过强W–I键连接的钨簇构成（0.57 BV单位），正是这种键赋予了该层坚固性。相反，从Cu(SCN)2中提取的2周期子结构由通过弱Cu–S键连接的链组成（0.05 BV单位）。因此，这样的层在剥离过程中很可能会破裂。为了确定max_intracomponent_bond_strength描述符的阈值，我们检查了MC2D中1633种易于剥离材料中的相应分布。由于数据集中99%的层的最大层内键合强度大于0.1 BV单位，并且在这个值附近的分布有明显的倾斜变化，因此将这个值设定为阈值。MC2D和CSS结果数据集中max_intracomponent_bond_strength描述符的分布报告在支持信息中。

图5
(a) 文氏图显示了四个来源中层的组成的交集。注意，只考虑了从晶体结构中获得的层——即所谓的自上而下的方法。每个数据集中独特层组成的总数如下：CSS结果——1683种组成，2DMatpedia——2504种组成；MC2D——2522种组成；C2DB——671种组成。此外，还展示了一组694个可靠层。
(b) new_prospective_exfoliable子集中的一些结构示例，这些结构具有更明显的2周期特性（W5I11，ICSD 244435，oP64）和主要为1周期的晶体结构（Cu(SCN)2，ICSD 131336，aP7），在其识别的2周期子结构内部有弱的Cu..S键合。固有周期性以饼图形式显示。

通过过滤其他数据库中不存在的929种独特层组成，我们得到了694个可靠层的最终名单（见Supporting Information.xlsx文件中的new_prospective_exfoliable_shortlist表），这些层存在于712个晶体结构中。这些可靠层分布在530种不同的结构原型中，与之前的研究相比，可能的结构模板数量增加了近40%，从而为基于元素替换的未来理论探索打开了重要领域。这些新识别出的材料大部分来源于之前因计算成本过高而被排除在外的大单元胞块结构。另一大一部分来源于最近才在实验数据库中沉积的结构。有趣的是，对这份名单的分析揭示了一些之前算法未能捕捉到的层状结构。例如，层状化合物Bi3Te2S（ICSD 107587，hP12）、Gd2Te5（ICSD 636468，oS28）以及不太明显的层状化合物如EuI2（ICSD 260561，oP12）、Cs3Bi2I9（ICSD 410726，hP28）或Cs2C4O4（ICSD 154357，mS20）并未在上述任何2D材料数据库中收录。

该算法有效性和普遍性的最终证明是它成功识别了一组最近出现的非范德华2D材料，这些材料在两篇近期出版物中被报道。（24,25）在这些研究中，计算了一组松弛结构的剥离能量，主要包括混合金属氧化物以及NaMnCl3、Al2S3和In2S3——所有这些都属于刚玉结构类型。总的来说，研究中提到的36种化合物中有30种出现在之前引用的经过过滤的ICSD子集中，并使用CrystalSubstructureSearcher算法进行了分析（见Supporting Information.xlsx文件中的non_vdW_2D表）。在这些化合物中，有27种成功识别出了中性2周期子结构。对于NaSbO3和KSbO3，算法检测到了带有碱金属阳离子的带负电的SbO3层。而对于Rh2O3，由于两组Rh–O距离的差异可以忽略不计，未能检测到层状子结构。

作为最后一步，为了提供提取层的性质概览，我们使用从头算方法优化并表征了550个提取的可靠层，排除了那些由于f轨道局部化而难以在普通DFT中正确描述的含稀土化合物。计算细节在方法部分有详细报告。由于这些层的尺寸比母体块状材料小，且维度降低，这种表征在计算上的成本比完整地从头算研究母体材料低一个数量级。在图6中，我们报告了缩放原子坐标的均方根位移以及提取层与理论优化层之间的2D单元胞面变化。这两个量可以作为提取层稳定性的定性估计，因为对于容易剥离的材料，预期在其内部原子排列和胞体在块状和孤立层之间不会有显著变化。有趣的是，对于成功优化的材料，所得到的值与MC2D数据库中的材料观察结果相当，进一步表明了提取层的高质量。

图6
(a) 缩放原子坐标在提取层和理论优化层之间的均方根位移分布。
(b) 提取层和优化层之间的单元胞面变化分布。

最后，我们对优化层的电子性质进行了基本表征，计算了沿着适当2D高对称路径的能带结构。完整的能带结构报告在支持信息中，而半导体和绝缘化合物的带隙分布报告在图7中。与参考文献（2）中获得的结果相比，分布峰值向更高带隙区域偏移（在DFT水平约为2 eV），这一区域特别适合光化学、透明和半透明电子学以及可能在室温下实现稳定激子的器件。值得注意的是，在高带隙区域（在DFT水平高于6 eV）也存在几个候选材料，这对于高性能电介质材料可能非常重要。

图7
图7. 在DFT水平下，绝缘和半导体2D材料的基本带隙分布。

在这项工作中，我们开发、基准测试并应用了一种灵活且计算效率高的算法，用于识别块状晶体中的低周期子结构并评估其剥离的难易程度。该算法实现在开源的CrystalSubstructureSearcher代码中，通过构建一个用反映键合强度的描述符加权的晶体结构图表示来运行。通过迭代移除最弱的链接来识别低周期子结构，同时监控结果子结构的维度。该算法在MC2D数据库的结构数据集上进行了验证，这些数据集中已经使用基于距离的算法识别出了2周期子结构，并且有基于DFT估计的层间结合能量。这种比较还使我们能够设计一种有效的策略，将DFT得到的结合能量与单位面积的BVS描述符相关联，以便快速近似估计剥离能量，而无需进行直接的DFT计算。我们进一步将这种方法应用于从ICSD和PCD数据库中筛选出的48,504个条目，其中包括在之前的研究中因计算成本过高而被放弃的大型结构。这使我们能够识别出694种坚固且易于剥离的二维材料，这些材料的成分在现有的二维材料数据库（MC2D、C2DB和2DMatpedia）中并不存在，从而展示了显著扩展已知二维材料领域的潜力，尤其是在那些对基于密度泛函理论（DFT）的方法具有挑战性的大单元胞系统中。此外，我们提出了一种晶体固有周期性的概念以及相应的描述符，用于表征晶体中的键合各向异性。该描述符定量反映了晶体结构的1-p/2-p/3-p特性。作为最后一步，我们优化并计算了新识别出的大部分二维材料的能带结构，提供了它们电子特性的大致概述，这有助于指导未来的应用导向筛选。作为一个副产品，我们获得了前98种元素之间所有可能键的R0 BV参数的估计值。预测的R0值与测试集中的真实R0值进行了比较，回归器的性能表现出色，如RMSE（0.034 ?）、R2（0.965）和MAE（0.026 ?）指标所示。这套新的R0 BV参数可以用于键价和映射方法，以便在固体电解质的晶体结构中可视化移动离子的扩散路径。

（26）利用机器学习回归模型扩展键价参数集
BV参数R0和b通常是通过使用实验晶体结构数据来拟合的，以确保给定原子周围的键价和等于其价态，这通常是针对结构中每个原子手动设置的。这意味着晶体化学中的键价参数拟合是基于共享电子的形式值。虽然可以拟合这两个BV参数，但通常只拟合R0参数，同时将b参数限制在0.37的固定值。然而，并非所有键的BV参数都是可获得的。目前可用的键价参数的全面列表来自专门的IUCr网页。（27）然而，在这个来源中，只有889个键的BV参数是可用的，且b = 0.37；因此，对于在晶体学数据库分析晶体结构时可能遇到的不太常见的原子对或金属对，实测数据是有限的。为了解决这一挑战，我们决定采用机器学习（ML）方法来估计前98种化学元素的所有4851个原子对的R0参数。我们的方法遵循了参考文献（28）中提出的类似策略，其中未知的R0参数被近似为相应元素大小参数r与校正项（Δij）的总和，该校正项包含了与原子电负性（EN）相关的第二个c参数（公式2）。研究中的b参数被固定为0.37，这是估计BV参数时的一个常见选择，尽管有研究表明优化b参数也可能有益，并能更好地建模键序（参考文献10,29）。
??????0=????+?????Δ????=????+?????????????(???????√????????√)2(????????+????????)
R0ij=ri+rj?Δij=ri+rj?rirj(ci?cj)2(ciri+cjrj)
（2）

??????0=??cov??+??cov???Δ????，其中Δ????=??(EN,??eff,period,vdw_crust_V_portion)
R0ij=rcovi+rcovj?Δij，?
其中Δij=f(EN,Zeff,period,vdw_crust_V_portion)（3）
同样，我们将r大小参数固定为元素的共价半径（30），并将b参数值固定为0.37，以便训练回归模型使用从mendeleev python库（31）中获取的一组原子属性作为独立变量来估计校正因子Δ（公式3）。原则上，可以使用多目标回归模型来获得更完善的R0和b BV参数集，但这将是未来研究的主题。在这里，我们专注于测试晶体结构分析算法以及扩展其中使用的BV理论在ICSD中所有化合物的适用性。
在获取的IUCr数据中，R0参数还取决于形成键的原子的氧化态和自旋态，尽管在大多数情况下这种依赖性是边缘的，这也得到了早期研究的支持（32）。尽管如此，对于过渡金属与O/N/C原子之间的键以及Tl–X键，观察到了相对较大的R0值范围。我们在训练ML模型时决定忽略这种依赖性，因为这样可以在没有关于元素氧化数和自旋态的额外数据的情况下分析晶体结构。这很重要，因为许多晶体结构缺乏这类信息，或者可能难以为结构中的原子分配特定的氧化态，例如铜硫族化合物（33）。因此，对于具有多个R0值的原子对，目标R0值是通过平均实测R0参数获得的。总体而言，这种方法为包含90种元素的874个键生成了一组实测R0参数数据。原始数据库中标记为“未检查”的参数被排除，因为它们不可靠。支持信息中的图S5显示了每种不同键类型的实测R0参数存在的键的比例。对于有足够实测数据（我们设置了该组中所有可能原子对的20%）的键类型，模型对R0的估计是可靠的。对于实测数据少于20%的原子对的键类型，模型对R0的估计被认为是不可靠的。记录了每个分析结构的具有不可靠ML估计R0参数的键的比例。
在收集、处理和准备BV参数数据之后，我们创建了用于后续模型训练的信息丰富的原子对特征表示。选定的8个原子属性集被用来描述形成键的每个原子。该集合包括五种不同的EN尺度：EN_Rahm-Hoffmann（34）、EN_Allred-Rochow（35）、EN_Nagle（36）、EN_Martynov-Batsanov（37）、有效核电荷Zeff、周期数以及vdw_crust_V_portion，后者计算为1?(??cov/??vdw)31?(rcov/rvdw)3，反映了原子体积中被vdW外壳占据的部分（21）。使用了最新的元素vdW半径集合（38）。由于每对原子都生成了两组原子描述符，特征的组合顺序很重要。为了解决排序问题，根据原子在Allred-Rochow EN尺度上的EN对每个原子对进行了排序；EN较小的原子的描述符在结果合并向量中位于前面。

CatBoost梯度提升决策树回归模型（39）被选为预测公式3中Δij校正因子的ML算法。该算法逐步将决策树添加到集成中，最小化RMSE损失函数并逐步提高其预测准确性。Δij值是通过从原子的共价半径之和中减去R0值获得的（目标Δij值的分布显示在支持信息中的图S6的直方图中）。对于键的特征表示，测试了多种技术来识别不同原子特征的最佳组合方式，包括原子特征的平均值、它们之间的差异计算以及简单连接原子特征向量。使用了5折交叉验证通过预定义的超参数集训练回归器（支持信息中的表S5）。结果表明，将原子对作为原子特征集的简单连接表示在预测校正因子Δij的任务中显示出最低的RMSE（支持信息中的图S7）。
在选择了键的特征表示之后，应用了一种简单的回归器训练策略。数据集被分为三部分：训练集包含655个条目，评估集包含120个条目，测试集包含99个条目。使用了耐心为10轮的早期停止策略来防止梯度提升回归模型的过拟合，并选择了集成中最佳数量的决策树。然后在训练集上以0.05的固定学习率训练模型，同时监控其在评估集上的性能。结果，使用了包含399棵树的最佳集成来对测试集进行预测。估计的R0值与测试集中的真实实测R0值进行了比较，使用RMSE（0.034 ?）、R2（0.965）和MAE（0.026 ?）指标来评估回归器的性能。最后，使用最佳超参数的模型在所有874个条目的数据集上重新进行了训练，并进一步用于估计98种元素中所有4851个原子对的校正因子Δij（见支持信息.xlsx文件中的ML_BV_parameters表）。模型预测的总体准确性与874个键的实测数据在支持信息中的图S8和图S9中显示。

ICSD和PCD数据集过滤
晶体结构数据来自ICSD（2023/2版本）和PCD（2014/15版本）数据库，并经过了一系列过滤步骤，这些步骤在图8的方案中简要展示。过滤后的结构根据它们的化学公式进行了分组，随后使用pymatgen中实现的StructureMatcher算法估计了结构相似性。结果，对于具有相同公式的每组化合物，识别出了一组子组，每个子组包含重复的晶体结构——同一晶体化合物的不同结构解或其特定多形体。在包含多个结构确定的组中，根据以下标准选择了单一代表结构：如果数据库条目中报告了Rf值，则选择Rf值最低的结构；如果组中所有数据库条目的Rf值都不可用，则选择发布年份最新的结构。

图8
图8. 从ICSD和PCD数据库过滤晶体结构的过程流程图。

高分辨率图像
下载MS PowerPoint幻灯片

在完成重复项去除过程后，最终数据集包含了一组精选的代表性晶体结构，每个结构代表一个独特且不同的晶体结构多形体。这些结构随后经过了进一步的精细化和验证过程。特别强调了包含CH、NH、OH、FH和SH键的结构，因为这些键的长度在X射线衍射研究中通常被低估了。执行了一个专门的归一化程序来调整这些结构中的氢原子坐标，使其与理想化的El–H键长度对齐。例如，使用了C–H的标准键长1.09 ?、N–H的1.01 ?、O–H的0.99 ?、F–H的0.95 ?和S–H的1.33 ?。如果键长小于理想值减去0.05 ?，则会进行键长重置。当检测到这种偏差时，H原子坐标会被调整以匹配通常从中子衍射研究中获得的理想化键长。

从头算计算
我们使用了Quantum ESPRESSO（40,41）软件包和SSSP PBE效率伪势库（版本1.3）（42），这是一个来自不同来源的经过测试的伪势库（43?47）。对于每种材料，波函数和电荷密度截断值是根据SSSP库为化合物中所有元素推荐的最高值选择的。对于块状材料，我们使用了两个经过适当vdW校正的泛函，即vdW-DF2泛函（48），其中包含c09交换项（49）。在优化和计算孤立2D单元的电子属性时，沿正交方向使用了20 ?的真空空间，并在开放边界条件下进行了仅涉及平面坐标的可变单元放松（50），并使用了PBE（51）泛函。使用Γ中心的Monkhorst–Pack网格对布里渊区进行采样（52），每个倒易晶格方向上的k点数最少，以确保至少0.2 ?–1的间距，并使用了Marzari–Vanderbilt–DeVita–Payne冷模糊（53）技术，冷模糊为0.02 Ry。能带结构是沿着高对称路径计算的，k点密度为0.01 ?–1。在自旋未极化的近似下，所有结构都被视为非磁性的，即使其中某些结构的原子可能呈现出磁有序的基态。所有的从头算计算都是通过AiiDA基础设施来进行的，以确保研究的完全可重复性。（54,55）