基于典型镧系金属-有机框架的高熵结构能量稳定性的计算分析

《The Journal of Physical Chemistry C》：Computational Analysis of the Energetic Stability of High-Entropy Structures of a Prototypical Lanthanide-Based Metal–Organic Framework

高熵金属有机框架（HE-MOFs）的稳定性研究通过DFT和机器学习势（MLIPs）方法结合，分析原型La系元素MOF的高熵版本。计算表明，仅少数异金属结构位于凸包上，说明高熵MOF的合成需平衡热力学稳定性和熵效应，且MLIPs（如CHGNet）能高效筛选候选结构，但需结合DFT验证。

　　高熵材料因其复杂的成分结构而受到广泛关注，通常由五种或更多元素以接近等原子比的形式组成。将这一概念应用于金属离子与金属有机框架（MOFs）的结合，为探索新型高熵MOFs提供了可能。尽管高熵材料的组成策略可以利用构型熵来增强热力学稳定性，但其复杂的成分空间和多样的相变特性也带来了显著的分析挑战。本文通过计算方法对高熵版本的代表性镧系MOF的多种复杂性进行了系统研究，利用密度泛函理论（DFT）和机器学习原子间势（MLIP）方法计算了这些异金属MOFs的混合能。通过MLIP方法，可以系统地探索包含最多五种不同金属的热力学稳定MOF结构的凸包。这种方法在计算效率方面具有优势，特别是在处理复杂的成分空间时，能够在不牺牲预测精度的前提下显著减少计算成本。

MOFs作为一类具有高度适应性和巨大表面积的材料，因其可调节的孔隙结构而受到广泛关注。这类晶体材料由金属离子或簇与有机连接体结合，形成一维、二维或三维的结构。随着新型材料设计和功能化策略的发展，一类被称为混合金属MOFs（MM-MOFs）或异金属MOFs的子类逐渐出现。通过结合不同金属，MM-MOFs可以展现出比其同金属对应物更广泛的化学和物理特性。尽管目前大多数异金属MOFs仅涉及两种金属物种，将三种或更多金属离子引入单一框架结构可以创造出具有更复杂成分的MOFs。这些材料的灵感来源于高熵合金，它们最早在2004年被引入。Yeh等人提出，多种元素的存在可以增加混合时的构型熵，甚至能够克服正的形成焓，这一现象在金属合金和纳米颗粒中得到了广泛研究。基于这些材料，MOFs中的多种金属物种也被称为“高熵MOFs”（HE-MOFs）。

在金属合金中，构型熵的增加可以促进固溶体的形成，但研究发现，仅仅最大化构型熵并不足以确保热力学稳定性。因此，有必要研究焓和熵的相对贡献，以全面评估高熵材料的热力学稳定性。DFT作为一种广泛使用的量子力学方法，已被广泛应用于MOFs的研究。在本文中，我们使用DFT计算了混合能，并与机器学习原子间势方法进行了对比。由于高熵材料的多维成分空间，使用DFT进行计算筛选存在挑战。例如，Chen等人研究了658,008个候选的单相等摩尔五元合金，使用40种金属元素的组合，最终只识别出30,201个（约5%）作为可能的单相等摩尔合金。这种高计算成本促使开发其他高通量筛选方法。近年来，机器学习原子间势和相关方法的发展为解决这一挑战提供了可能。通用的MLIP可以作为DFT的替代方法，用于探索高熵材料的成分空间。这些势函数通常基于大量的无机化合物数据集进行训练，尽管在泛化方面仍存在挑战，但它们仍然是研究复杂问题的有力工具。例如，Vazquez等人使用图神经网络（GNN）模型和簇扩展方法，研究了（Cr,Hf,Mo,Ta,Ti,Zr）B2和难熔高熵合金AlHfNbTaTiZr的短程有序性。

目前关于高熵MOFs的形成研究仍处于初步阶段，实验和模拟中缺乏相关数据。因此，我们对高熵MOFs的合成或设计中涉及的关键因素仍存在认知空白。本文使用DFT和MLIP方法分析了代表性镧系MOF的高熵变种的构型景观。我们的研究结果为寻找稳定高熵MOFs提供了某些物理原理，并为这一材料类别未来的发展方向提供了建议。我们首先选择了五种不同的金属物种，从镧系元素中生成了异金属MOFs的结构，并使用DFT和MLIP方法（M3GNet/CHGNet）进行了柔性晶格松弛。对于MLIP方法，我们使用了Broyden-Fletcher-Goldfarb-Shanno（BFGS）算法进行结构优化。为了保持一致性，DFT和MLIP方法均使用了相似的力收敛标准。尽管每个MOF的单位晶胞中包含的金属种类不同，但为了简化，我们仅通过单位晶胞中的总金属组成来表示整个结构。例如，一个包含等原子比Nd和Tb的MOF结构被表示为Nd3Tb3。

在进行比较时，由于不同方法中参考能量点的任意性，绝对能量值的直接比较并不具有意义。然而，我们仍然可以研究单个计算框架内不同结构的能量趋势。例如，通过将Tb6-MOF的总能量作为参考点，我们计算了不同结构的相对能量，如表S1所示。对于DFT计算，能量值的范围可以达到±2 eV/单位晶胞，但MLIP计算中某些金属（如Eu、Gd和Yb）的能量值出现了显著的异常。这可能与我们在计算中使用的赝势与材料项目数据集中训练MLIP所使用的赝势存在差异有关。因此，我们选择了Ce、Nd、Sm、Tb和Er这五种金属作为进一步研究的对象。我们还评估了CHGNet在预测高熵MOFs的混合能时的准确性，发现其平均绝对误差（MAE）为10.1 kJ/mol·单位晶胞，表明在不考虑构型熵的情况下，CHGNet的预测能力仍有局限。然而，这一误差并未显著影响整体趋势的判断。

为了进一步分析构型熵的影响，我们在300 K的条件下加入了理想构型熵项，–nR（xi × ln xi），其中n是不同金属物种的数量，xi是第i种金属的摩尔分数。通过这种方式，我们计算了不同金属比例下的混合能，并将其投射到一个二维正五边形上，如图4a所示。加入构型熵后，稳定结构的数量从97增加到207，其中78个异金属结构也位于凸包上。这表明构型熵对高熵MOFs的稳定性有轻微的促进作用。然而，即使在考虑构型熵后，仍有49个（约0.6%）可能的4金属结构和5个（约0.3%）可能的5金属结构位于凸包上，说明这些结构在整体结构数据库中仍占极小比例。因此，随机选择的结构不太可能位于凸包上，而是处于亚稳态。

我们还评估了构型熵在900 K时的影响。由于900 K远高于MOF分解的温度，因此这种高温条件为构型熵效应提供了上限。温度对稳定结构数量的影响并不显著，整体分布与300 K时的情况相似，如图S1和S2所示。考虑到我们结构数据库中大量结构处于亚稳态，我们计算了它们的能量与凸包之间的距离（Ehull），即垂直距离。对于稳定结构，Ehull值为零；而对于亚稳态结构，Ehull值大于零。我们评估了不同金属组合下的Ehull值，并将其与DFT计算进行了对比，发现这些值在所有情况下均低于35 kJ/mol·单位晶胞。2-3金属MOFs的Ehull值多集中在较低范围，表明它们与凸包的接近程度较高。相比之下，4-5金属MOFs的Ehull值更多分布在较高范围内，这表明这些结构在热力学平衡下可能更容易分解为含有更少金属物种的更稳定子簇。我们还发现，加入构型熵后，Ehull值的分布变化不大。

为了进一步验证凸包计算的可靠性，我们使用DFT对CHGNet预测的稳定结构进行了评估。由于这些结构的混合能值为负，我们还随机选择了具有正混合能值的结构，并与DFT结果进行对比。我们将CHGNet的混合能值分为三个类别：–15到–5 kJ/mol·单位晶胞（类别1）、–5到10 kJ/mol·单位晶胞（类别2）、以及10到25 kJ/mol·单位晶胞（类别3）。我们对每个类别中的8个结构进行了DFT计算，总共有11个结构的混合能值小于零并位于凸包上。这些24个结构的混合能值与DFT结果之间的关系如图7所示。对于位于凸包上的结构，混合能值与DFT结果之间表现出半定量的一致性，其平均绝对误差（MAE）为5 kJ/mol·单位晶胞，并与之前的研究结果相符。然而，在类别3中，DFT混合能值与力场预测值之间的变化更大，MAE达到18.5 kJ/mol·单位晶胞。当将这些24个结构的数据与前文提到的26个结构合并后，总的MAE为8.3 kJ/mol·单位晶胞。

我们还使用DFT计算了这些稳定结构的带隙，发现它们的带隙范围在2.77–2.81 eV之间，表明金属的变化并未显著改变这些MOFs的带隙特性。此外，我们还探讨了单金属替换对凸包分析的影响。我们选择了Pr作为Ce的替代，并对Pr–Nd–Sm–Tb–Er MOFs进行了额外的凸包计算。我们重复了图2中的对比图，将Pr基结构替换为Ce，并进行了结构优化。最终的MAE为5.0 kJ/mol·单位晶胞，相关结果如图S4所示。我们还评估了Pr–Nd–Sm–Tb–Er MOFs的CHGNet混合能值，发现其分布范围在–15到25 kJ/mol·单位晶胞之间，如图S5所示。经过全面的凸包分析，我们发现总共有113个MOF组成和结构位于凸包上，其中没有5金属组成。在加入理想构型熵后，稳定结构的数量增加到204，其中包含5个5金属MOF结构。这些结构的组成及其金属比例分布如图S6和S7所示。亚稳态结构的Ehull值分布如图S8所示。最终，我们比较了CHGNet和DFT在部分稳定结构上的混合能值，如图S9所示。位于凸包上的结构的MAE为2.6 kJ/mol·单位晶胞，而24个结构的总体MAE为6.7 kJ/mol·单位晶胞。当将这些24个结构的数据与图S4中的26个结构合并后，总的MAE为5.9 kJ/mol·单位晶胞。

尽管存在一些定量差异，但Pr–Nd–Sm–Tb–Er高熵MOFs的凸包分析趋势与Ce–Nd–Sm–Tb–Er高熵MOFs的分析趋势一致，这表明我们的结论在其他具有相似相对能量的镧系元素中具有可迁移性。我们发现，在将Ce替换为Pr后，MAE值有所降低，与表S1的结果一致。

此外，我们还探讨了金属短程有序性对材料能量稳定性的影响。即使具有相同的组成，不同金属物种的构型排列也会导致混合能值的显著差异。例如，对于组成CeNd2SmTbEr的MOF，混合能值的变化范围为–16到+16 kJ/mol·单位晶胞。为了研究这些短程有序性效应，我们考虑了更大的超胞结构。由于在研究MOF时，单位晶胞中的原子位置数量有限，难以捕捉短程有序性现象。因此，通过考虑更大的超胞结构，可以更有效地研究这些效应。我们进行了系列计算，以探讨这些效应在原型高熵MOF中的表现。

为了研究与单位晶胞层面的组成变化相关的混合能变化，我们算法性地探索了包含多个单位晶胞的超胞结构的特性。具体步骤如下：
1. 我们从结构优化的单位晶胞中生成一个包含多个单位晶胞的超胞。
2. 在每次迭代中，我们通过交换两种不同金属物种的位置生成新的结构，并对所有可能的金属组合进行单位置交换，生成一组新的结构。
3. 我们对所有这些新结构进行结构优化，并计算它们的混合能值。为简化计算，我们仅优化原子位置，保持体积不变。
4. 在每次迭代后，我们确定当前步骤中的最低混合能值，并使用对应的结构进行下一步迭代。在下一步迭代中，我们重复这一过程以进一步生成新的构型。

为了展示这一方法，我们选择了一个组成CeNd2SmTbEr的结构，其初始混合能值为–1.3 kJ/mol·单位晶胞，且处于凸包上方，为亚稳态结构，其Ehull值为9.7 kJ/mol·单位晶胞。我们生成了其2 × 1 × 1超胞，并进行了10次迭代。初始和最终构型如图8所示。经过多次金属位置交换后，最终构型的混合能值显著降低至–7.4 kJ/mol·单位晶胞。最终构型中的两个六核簇具有不同的组成：一个为CeNd2Sm2Tb（单单位晶胞混合能值为–6.6 kJ/mol，Ehull值为1.1 kJ/mol），另一个为CeNd2TbEr2（单单位晶胞混合能值为–12.3 kJ/mol，Ehull值为0.0 kJ/mol）。因此，初始的5金属构型分解为两个更稳定的4金属子簇。

我们还对来自Ce–Nd–Tb–Sm–Er成分空间的三个其他亚稳态结构进行了类似的计算，分别为4金属、3金属和2金属MOF结构。在图S10中，我们展示了这些结构及其混合能值在10次迭代后的变化情况。我们观察到类似的金属交换或聚集趋势，这些趋势导致混合能值显著降低。在另一个场景中，我们考察了这些结构特性对几个2 × 1 × 1超胞的影响，这些超胞被识别为在之前的单单位晶胞凸包分析中可能稳定的结构。我们选择了Ce–Nd–Tb–Sm–Er成分空间中的一个5金属和一个4金属MOF结构，并展示了它们在10次迭代后的混合能值变化，如图S11(a,b)所示。在每种情况下，最低能量结构相对于初始构型仅包含一次金属交换，混合能值的变化不超过0.1 kJ/mol·单位晶胞。我们还对位于凸包上的2金属和3金属MOF的2 × 1 × 1超胞进行了分析，发现类似的观察结果。

这些基于代表性结构的凸包分析表明，我们对单单位晶胞的计算足以估计多组分成分空间的整体凸包。这在该原型MOF中是合理的，因为六核簇的物理分离意味着簇间相互作用（如连接簇的链接体中的局部应变）对混合能的贡献可能小于直接形成簇的键能。因此，通过MLIP方法进行系统探索是可行的，尽管混合能值的变化可能较小。此外，我们还注意到，我们使用的MLIP（CHGNet）仅适用于部分镧系金属。如果要在实验研究中预测特定MOF的组成，使用MLIP生成的凸包进一步进行DFT计算以提高整体预测的精度是明智的选择。

总之，目前对高熵MOFs的研究仍处于早期阶段，其中主要挑战之一是其广泛的成分空间，使得实验试错法难以用于其合成。此外，单位晶胞层面的局部有序和无序现象也难以直接实验探测。这表明利用高通量计算方法可以更好地理解这些材料的复杂特性。我们使用DFT和MLIP方法对代表性镧系MOF的高熵变种进行了研究，尽管父晶体结构是同金属的，但我们考察了包含最多五种不同镧系金属的异金属结构。通过比较混合能值并进行凸包计算，我们识别了热力学稳定的结构。虽然这些计算无法提供关于实际MOF结晶过程中可能发生的动力学效应的信息，但我们的结果表明，许多具有多金属组成的结构中只有极小部分位于凸包上。这提示在选择预设组成以验证合成可行的高熵MOFs候选结构时，存在一定的复杂性。我们还探讨了金属有序性对材料能量稳定性的影响，发现这些效应在单单位晶胞结构上可能不显著，但在接近凸包的结构上变化较小。我们的分析局限于特定的原型MOF，但提供了关于预期趋势和复杂性的见解，而不是对这一材料或其他高熵材料做出精确的定量预测。由于可能的组成变种数量众多，我们依赖MLIP计算方法系统地探索广泛的结构。我们展示了所使用的MLIP方法（CHGNet）与DFT计算结果之间的合理一致性，但不同组成之间的能量变化可能较小。此外，我们还发现，所使用的MLIP（CHGNet）仅适用于部分镧系金属，如果要预测特定MOF的组成，使用MLIP生成的凸包进一步进行DFT计算以提高整体预测的精度是明智的选择。