离子相关性解释了在扩散受限的固态合成反应中的动力学选择性

《Nature Materials》：Ion correlations explain kinetic selectivity in diffusion-limited solid-state synthesis reactions

　　摘要：建立新型无机材料的可行固态合成路径仍然是材料科学中的一个主要挑战。以往的路径设计方法采用成对反应的方式，并利用第一性原理数据来分析热力学特性，但缺乏动力学信息，这限制了它们的有效性。这种差距导致了前体选择和预测的不理想，尤其是在形成能量相似的竞争性相的反应中，离子扩散是一个关键因素。在这里，我们通过将机器学习得出的传输特性融入到“类液体”产物层的热力学细胞反应模型中，展示了一个无机合成框架。在以竞争性多态性著称的Ba–Ti–O系统中，我们获得了随时间和温度变化的BaO:TiO2比例下的相形成准确预测。我们发现，扩散-热力学相互作用决定了相组成，其中跨离子传输系数对于预测受扩散限制的选择性至关重要。这项工作通过整合固态反应动力学、第一性原理热力学和空间反应性，跨越了不同的长度和时间尺度。

由于对改进能源技术的新型无机材料的需求，强大的数据驱动材料设计应运而生。然而，由于缺乏通用的固态反应机制模型，通过无机合成实现计算机设计的材料仍然是一个挑战，这落后于有机合成。目前最先进的固态合成原子级建模是根据高通量密度泛函理论（DFT）计算数据库（如Materials Project (MP)）中的批量热力学性质来描述反应行为的。显著的例子包括反应网络（识别热力学上有利的路径）和基于实验反馈迭代改进合成配方的主动学习算法。然而，仅基于热力学的预测可能不准确，特别是在形成能量相似的竞争性相的系统中。在这种情况下，关键成分的有限传输可能会阻止全局稳定产物的形成，从而阻碍热力学平衡的达成。

先前的尝试通过使用经验速率表达式来理解这些效应，这些表达式是根据反应物的转化程度来拟合“有效速率常数”，因此是事后确定的，并不能预测固态反应产物。对这些反应的计算研究主要依赖于相场方法，这需要解决耦合的多物理场偏微分方程。以往的工作集中在简化场景上，包括无中间体的单步反应、烧结驱动的形态演变以及在熔体中生长预定义的成核相。尽管这些研究提供了有价值的信息，但它们依赖于实验或建模输入，因此需要特征明确的系统来解析扩散和成核动力学。

一般来说，固态反应动力学可以分为成核限制和扩散限制两种情况。虽然前者可以洞察形成的第一相，这在例如薄膜合成中很有用，但这些信息并不能确定粉末反应中产物的整体分布，后者是通过前体成分向反应区的扩散控制传递来进行的。我们假设这类系统的合成演变可以描述为在时间依赖的离子流限制下，通过具有与候选成核相同化学计量的缺陷类液体界面来优化局部能量。简单来说，我们提出，在成对粉末反应中，晶体核的生长受其成分通过具有相同组成的类液体界面的离子传输的调控。

为了展示我们的模型，我们选择了具有特别竞争性多态性的Ba–Ti–O化学空间（图2a）。技术上相关的Ba–Ti–O三元化合物（例如BaTiO3、Ba2Ti9O20、BaTi5O11和BaTi2O5）通常是通过以不同比例混合BaCO3和TiO2前体来合成的。铁电体BaTiO3的形成已经被广泛研究，通常通过在1,000–1,300?°C下加热混合粉末来实现。尽管BaTiO3是主要产物，但在TiO2/BaO界面处，Ba2TiO4在能量上更受青睐（图2a）。杂质的存在取决于合成条件：在1,050?K以下Ba2TiO4占主导地位，而在更高温度（约1,200?K）下会增加BaTiO3的产率，但会促进BaTi2O5等次要相的形成。

图2：Ba–Ti–O化学空间中的热力学和动力学景观。图中蓝色方块（除Ba3TiO5外）是在固态反应中实验观察到的相，红色方块是MP中的额外相。该图的稳定性外壳形状（即相及其相对形成能量）在本研究考虑的所有温度下保持不变。

a. 在600?°C下，使用MP中的数据计算得到的BaO和TiO2的热力学稳定性边界。有限温度效应通过机器学习估计器来考虑振动对熵的贡献。a，BaO和TiO2在600?°C下的热力学稳定性边界。b. 计算得到的有效扩散速率常数（KD），它是衡量所有离子通过反应产物相在1,000–1,750 K温度范围内BaO–TiO2反应界面的平均通量的指标，按Ba:Ti比例递增排序。KD在Ba:Ti比例中等的相处达到峰值，在Ba:Ti?>?1的相中显著下降。KD的误差条（可见时）是使用五个MD轨迹的标准偏差计算得出的。

先前的研究，特别是最近的工作，基于第一性原理热力学数据分析了固态合成反应中的选择性。参考文献5提出了一个热力学阈值，即每个原子60?meV以上，可以可靠地预测初始产物。相比之下，对于多个具有相似驱动力的竞争性相的系统，初始产物通常由动力学决定。对于Ba–Ti–O系统，形成驱动力最高的产物（Ba2TiO4）与BaTiO3和BaTi2O5之间的形成能量差分别为每个原子约51?meV和46?meV。因此，我们预计该系统的合成结果将受到扩散通量和反应能量之间温度依赖性相互作用的影响。实际上，实验上通常首先观察到的产物是Ba2TiO4；然而，随后会发现驱动力较低的产物，如BaTi2O5和BaTiO3。以Ba–Ti–O作为一个严格的测试案例，我们在这里提出了一个通用的、预测性的固态合成结果框架，该框架将严格计算的离子传输和第一性原理热力学与细胞自动机（ReactCA）相结合，这是一种离散的计算模型，其中网格单元根据局部邻居相互作用进行演化（图3）。如图所示，预测结果与四个经过仔细表征的实验在时间和温度上的结果非常吻合。

图3：用于模拟固态粉末合成反应的ReactCA仿真框架概述。a，通过成对界面反应进行反应的一个步骤的进展。b，该框架的主要阶段示意图：(i) 用户指定要研究的合成配方（包括要考虑的化学系统、前体及其用量、反应气氛和加热曲线），根据这些信息从MP获取形成能量，并使用机器学习估计器计算振动自由能、熔点和扩散速率常数，然后对反应进行枚举和评分；(ii) 生成粒子的随机初始分布（体素），并使用分配给每个反应的评分重复应用演化规则来模拟反应；(iii) 将仿真步骤串联成一个轨迹，分析该轨迹以确定反应过程中的相分数和所遵循的路径。

在Ba–Ti–O相空间中，我们考虑了两种反应物（这里是BaO和TiO2）以及Ba–Ti–O空间中的九种可能的固态反应产物（图2a，蓝色方块）。我们使用BaO而不是BaCO3，因为已知实验前体BaCO3在约1,100?K时会分解为BaO（参考文献17），在任何三元反应发生之前。我们使用Onsager分析从界面上的化学势差计算成分离子的通量，以及Ba2+、Ti4+和O2?的传输系数。这些分析基于通过机器学习原子间势（AIMD）数据为Ba–Ti–O系统中九种可能的类液体非晶类似物生成的5-ns分子动力学（MD）轨迹。

图2b显示了在1,000–1,750?K温度范围内，BaO/TiO2界面处九种可能产物组成的KD。有趣的是，在1,000?K以上，Ti富集相中的KD比Ba富集相中的KD高一个数量级以上。在Ti富集相中，KD随着温度每升高250?K增加一个数量级，并在1,750?K时大多数相达到一个相似的值。而在Ba富集相（Ba:Ti比例?>?1）中，KD仅随着温度升高750?K增加一个数量级。

为了将产品的温度和时间依赖性演变与其它地方报告的合成实验进行比较，我们使用了最近开发的细胞自动机仿真框架ReactCA。ReactCA旨在通过模拟基于可定制的局部规则演化的三维网格单元来模拟固态反应，这些规则结合了热力学和动力学输入，通过反应速率的代理（“评分”函数）来实现。在这里，我们扩展了ReactCA框架，允许评分函数依赖于瞬时生长率，该生长率是液体产物相在温度T下的有效离子扩散常数（KD）、修正的热力学驱动力（\(\frac{\Delta {G}^{* }}{{k}_{{\rm{B}}}T}\)和Tammann规则的启发式函数的函数。在Tammann温度以下，反应速率较低但仍然可能。在其以上，由于扩散和热力学的共同作用，速率随温度增加而增加；然而，在高温下，扩散率的饱和使平衡偏向于热力学控制的结果。我们使用与实验相同的加热曲线（图4和5）模拟了从1:5到1:1的前体Ba:Ti化学计量的反应。

图4：使用ReactCA进行的反应模拟，考虑了扩散通量和反应热力学与仅基于热力学的反应。a，反应1（BaO?+?TiO2）(a) 和反应2（BaO?+?5TiO2）(b)。灰色背景的深浅区分了不同的反应选择性机制：浅灰色区域表示活化控制机制，深灰色区域表示动力学控制机制，深色区域表示热力学控制机制。

反应1（图4a）代表了典型的1:1 BaO到TiO2的合成反应，包括加热到1,500?K后进行短暂的退火过程，结果主要是BaTiO3。在我们的模拟中，低温下首先形成的产物是Ba2TiO4和BaTiO3，这与实验结果非常吻合，实验表明在1,000到1,050?K之间的反应中Ba2TiO4是主要杂质（高达27?mol%）。在更高温度下，我们发现Ba2TiO4被消耗，而BaTi2O5成为主要杂质，这一结果与参考文献15,16中的数据一致。随着温度的继续升高，BaTi2O5被大量消耗，但仍作为主要杂质存在（<5?mol%）。实际上，在参考文献21中，当反应在1,300?K下进行超过100?分钟后，反应完全进行，得到了几乎纯度的BaTiO3。如果不包括扩散率，进行相同的模拟会得到不同的结果。仅基于热力学数据和Tammann规则的模拟结果显示，在1,500?K退火后，产物混合物包括BaTiO3（约70?mol%）、Ba2TiO4（约15?mol%）和BaTi2O5（约20?mol%）。特别是，BaTiO3仅占产品的70%，这与实验观察到的几乎纯相的BaTiO3形成了鲜明对比。反应2（图4b）遵循了之前的实验22，其中BaO和TiO2以1:5的比例在一系列温度（1,250–1,500 K）下反应。在低温（<1,300 K）下，我们的模拟表明最初会形成少量的BaTiO3（<10 mol%）和BaTi2O5（约10 mol%），然后它们转化为BaTi5O11、Ba2Ti9O20和中间体BaTi4O9。BaTi4O9在1,500 K时达到约15 mol%的最大值，随后被高温产物BaTi5O11和Ba2Ti9O20以及剩余的前体TiO2消耗掉。有趣的是，在参考文献22的最终产物分布中只发现了BaTi5O11、BaTi4O9和未反应的金红石（TiO2）。由于他们没有指定退火或反应时间，我们推测他们的实验对应于我们模拟的早期阶段，在Ba2Ti9O20大量形成之前。如果在更高温度下继续反应更长时间，将会促进Ba2Ti9O20的形成，这也与参考文献23中的实验结果一致，在那里Ba2Ti9O20的形成是以牺牲BaTi4O9和TiO2为代价的。在参考文献18中，BaO和TiO2以2:9的比例（反应3）反应，通过快速将温度升高到约1,400 K，然后进行3小时的退火处理。样品随后在大约1,700 K下烧结6小时，得到了含有BaTi5O11和BaTi4O9杂质的Ba2Ti9O20产物。特别是，Ba2Ti9O20主要是以牺牲BaTi5O11和BaTi4O9为代价形成的。在参考文献18中，他们还观察到了这三个相的形成速率：BaTi5O11 > BaTi4O9 > Ba2Ti9O20。在我们的模拟（图5a）中，最初形成少量的BaTiO3和BaTi2O5，然后迅速转化为Ba2Ti9O20、BaTi5O11和BaTi4O9。特别是，我们的结果成功地再现了实验观察到的初始形成速率趋势，即BaTi5O11 > BaTi4O9 > Ba2Ti9O20。此外，与实验结果定性一致的是，在1,400 K的退火过程中BaTi5O11和BaTi4O9的生长停滞，而在1,400 K以上则持续以正速率生成Ba2Ti9O20。在参考文献18中，当反应在1,400 K以下进行时，还检测到了微量的Ba4Ti13O30和一些未反应的金红石（TiO2），这些也在我们的模拟中观察到。

图5：使用ReactCA进行的反应模拟，考虑了不同的BaO和TiO2前体比例，同时考虑了扩散通量和反应热力学以及仅考虑热力学的影响。

参考文献24提供了一种通过将1:2的BaO:TiO2前体在三个连续的温度步骤（约1,200 K、1,300 K和1,400 K）下退火来合成纯相BaTi2O5的配方，我们通过反应4（图5b）进行了模拟。在第一步之后，他们的样品主要含有BaTiO3和一定量的BaTi2O5。在第二步退火之后，BaTiO3的含量减少，而BaTi2O5的比例增加。在第三步退火之后，观察到了接近纯相的BaTi2O5，BaTiO3和Ba6Ti17O40作为杂质相出现。在1,400 K以上更长的加热时间内，BaTiO3和Ba6Ti17O40（1,473 K以上的平衡相）的相分数增加。我们的模拟显示了极好的一致性；我们展示了在第一步退火过程中BaTiO3迅速成为主要相，而在第二步退火过程中被BaTi2O5取代，最终得到了含有微量（<5%）Ba6Ti17O40和BaTiO3杂质的纯相BaTi2O5。最后，在研究空气中通过传统固态反应形成BaTiO3的过程中，参考文献15提到在反应的早期阶段同时形成了BaTi4O9和Ba4Ti13O30相。他们将这些杂质相的存在归因于动力学效应，即尽管这些相的形成热力学驱动力不大，但由于它们的快速扩散动力学可以成核和生长。我们的模拟与实验观察结果非常吻合，显示了这些相与BaTiO3的同时形成，其出现取决于整个系统的组成。富钛的系统主要产生微量的BaTi4O9（反应4），而组成更平衡的系统产生微量的Ba4Ti13O30（反应1）。

通过将严格计算的类液体反应相间的离子扩散系数整合到有效的扩散率常数（KD）中，我们的方法改进了基于温度和时间的固态合成路径和结果的先验预测。图4和图5与经过仔细表征的实验结果非常吻合，准确地捕捉了产品和中间体作为温度和时间的函数。特别是，当仅使用热力学反应性和Tamman规则进行相同模拟时，相同的框架通常能预测正确的主要产物，但无法准确捕捉中间产物和主要杂质，这些杂质要么缺失，要么被错误识别。实验和从头算预测之间的强一致性支持了这样的假设：通过类液体相间的传输与目标相相匹配，很好地近似了实验条件。利用这一见解，我们分析了不同类液体组成对扩散率的影响，预计这种影响会表现出强烈的离子相关性25,26。尽管在之前关于通过浓缩离子液体传输的研究中已经详细研究了这种效应27,28，但在固态反应中间体的讨论中还缺乏这种性质的分析。为了展示相关离子运动的影响，我们绘制了通过阳离子（Ba2+和Ti4+）的自扩散部分归一化的计算出的不同离子Onsager传输系数。图6a显示了归一化的不同Ba离子和Ti离子传输系数，与每种组成中类液体相间结构中每个原子的Ba–O、Ti–O、Ba–Ba和Ti–Ti配位分数进行了比较。数据揭示了Ba2+和Ti4+之间的相关运动程度与局部配位环境之间的强烈趋势。随着系统中Ba:Ti比例的增加，Ti4+的环境被O2?主导，这导致TiOx簇的比例增加。这反过来又导致每个Ba原子的Ba–Ba配位数增加（图6b）和Ba–Ba相关运动的程度增加。这一观察与计算出的不同离子相关性一致：在富Ba相中，Ba2+离子显示出负相关性，表明在化学势梯度下存在相互排斥；而Ti4+离子显示出较小的正相关性，表明它们有较小的聚集倾向。我们还注意到，与Ba2+相比，Ti4+的尺寸较小，以及钛在玻璃态结构中形成网络的已知倾向29，这与改变系统Ba:Ti比例时传输的变化一致。因此，由于在富Ba相间中存在交叉相关运动，Ba2+的净移动速度减慢，其移动性降低。在富Ti相中也出现了类似但较弱的效果，导致扩散率常数在中等Ba:Ti比例时达到峰值（图2b）。为了展示Ba和Ti离子之间的相关运动如何影响合成过程中的选择性相形成，我们比较了包含和不包含交叉离子效应的模拟，即仅使用离子移动性的自扩散估计（补充图1）。尽管结果在定性上与图4和图5中的数据相似，但在初始阶段涉及更多富Ba相时存在关键差异。例如，反应1没有形成实验验证的Ba4Ti13O30中间体15。在反应3中，加入交叉离子通量减少了初始富Ba中间体的数量，这与早期的观察结果一致18。此外，在反应3中，最终产物分布中仅观察到Ba2Ti9O20和BaTi4O918，当包含交叉传输项时，这一趋势在模拟中得到了更准确的捕捉。当仅考虑自传输项时，模拟高估了Ba2Ti9O20的数量，低估了BaTi5O11的数量，反之亦然。我们得出结论，尽管自扩散项捕获了大多数动力学效应，但结合交叉离子通量对于准确捕捉相对相形成至关重要。在这些涉及多个阴离子基团的反应（如离子交换或交换反应）中，这些交叉效应预计会更加显著30,31。

图6：在1,500 K下考虑的类液体系统中的相关效应分析。

此外，通过分析仅基于热力学数据的模拟与结合动力学信息的模拟之间的差异，我们可以识别出Ba–Ti–O系统中的三种选择性温度区域：激活控制区、动力学控制区和热力学控制区。在低温（对于Ba–Ti–O系统低于约1,100 K）下，所有产品的扩散动力学都很慢，与前体比例相当且具有显著反应能量ΔGrxn的相被优先形成（激活控制区）。在中等温度（1,100 K < T < 1,700 K）下，扩散动力学有利于某些组成在可能相中的形成（动力学控制区）。在高温（T > 1,700 K）下，接近系统的熔点时，所有离子物种的动力学都很快，最接近热力学设定组成的相混合物形成热力学平衡。这一见解与现有的固态反应性理论一致1,32,33,34。例如，BaTiO3和BaTi2O5在壳层中富Ti的一侧表现出较大的负形成能量和相对较快的动力学（图2），倾向于在较低温度下形成（反应2和3）。将温度升高到1,700 K以上，对应于热力学区域，促进了BaTi5O11和Ba2Ti9O20的形成，这是给定前体比例下热力学上最有利的产品混合物。最后，我们评估了固态合成中结构相似性的影响。几项先前的工作5,18,35认为，前体、中间体和产物之间的结构相似性或“模板作用”促进了特定相的形成。实际上，在参考文献18中，假设Ba2Ti9O20仅在1,373 K以上形成，这是因为中间体BaTi5O11通过模板作用降低了界面和应变能量。在这里，我们再现了参考文献18中的产物分布，没有纳入任何结构相似性效应。具体来说，我们发现反应2和3的早期阶段（图4和5）中三种三元相的相对形成速率相似（BaTi5O11 > BaTi4O9 > Ba2Ti9O20），并且重要的是，Ba2Ti9O20仅在1,373 K以上热力学上稳定，如BaO–TiO2相图所示36。其形成的温和驱动力（图2）解释了它在参考文献18和我们的模拟中仅在1,373 K以上温度出现的原因，以及在其他人使用的较短反应时间中未出现的原因16。为了总体分析Ba–Ti–O系统中的结构效应，我们使用了参考文献35中计算出的成核率。他们的指标显示，BaTiO3和BaTi2O5在低温到中等温度下表现出最低的成核障碍，应该占主导地位（补充表1）。然而，在两个实验和模拟中，Ba2TiO4通常作为低温中间体形成，这无法用这种方法解释。因此，尽管结构相似性为进行拓扑化学的合成反应提供了一个有价值的概念框架，但我们发现，当反应物发生显著的结构转变和/或完全分解时，所提出的方法能更好地解释产物结果。我们预期我们的方法将非常适合研究这种热力学-动力学相互作用在更多受扩散限制的固态合成反应中的表现，并且能够以高吞吐量执行此类模拟，以寻找现有和未实现材料的更好配方，这些材料可能由于动力学原因而无法形成。

为了捕捉密集的类液体或非晶相中的离子动力学效应，我们采用了参考文献38中的框架，以获得自相关和相关的扩散传输系数，形式为Onsager传输矩阵。通过使用基于原子簇展开（ACE）的机器学习原子间势能39，并在75个来自非晶系统的AIMD轨迹上进行训练，得到了传输系数的可靠估计。所得的动力学数据以及来自MP的热力学数据被实现到ReactCA模拟框架19中，该框架允许描述在规定的固态反应过程中空间和时间的相演变。下面，我们将更详细地描述模拟框架的组成部分。

我们使用atomate2 Python包中实现的MPMorph40,41工作流程来生成用于训练ACE的非晶配置。首先使用packmol生成所需组成的随机非晶结构。MPMorph工作流程首先将此结构的体积调整为初始体积的0.8倍和1.2倍，然后进行4皮秒的NVT AIMD运行以确定给定温度下的状态方程。然后根据状态方程选择一个试验体积，再进行4皮秒的NVT AIMD运行以检查能量是否收敛。如果能量收敛，则进行20皮秒的AIMD“生产”运行以在该体积下平衡结构。如果能量不收敛，则继续迭代调整体积，直到达到能量收敛。一旦找到一个收敛的体积，就使用它进行20皮秒的生产运行。该工作流程使用jobflow和fireworks Python包在高性能计算资源上运行。

根据BaTiO3、Ba2TiO4、BaTi2O5和BaTi5O11的组成，在1000、1250和1500 K下分别进行了实验，并且每100飞秒采样一个帧，以确保采样配置的多样性。这些采样帧与MP数据一起用于在Ba–Ti–O化学空间中训练ACE势能，使用pacemaker42,43 Python包。ACE势能每个元素用500个基函数参数化，考虑5 ?范围内的邻域列表。遵循参考文献44的方法，使用更高阶的Finnis–Sinclair嵌入函数。通过基于幂次的“阶梯拟合”方法43进行训练，其中对力的损失贡献赋予更高的权重（99%），共进行了2000次迭代。一旦训练收敛，就通过将力的损失权重降低到90%再进行另外2000次迭代，并保持势能的形状和复杂性不变。

我们使用外推等级和高温MD方法45进行了主动学习。我们使用packmol为BaTi5O11、Ba2Ti9O20、BaTi4O9、Ba4Ti13O30、Ba6Ti27O40、BaTi2O2、BaTiO3、Ba2TiO4和Ba3TiO5的组成生成随机堆积的结构，作为基于ACE的NVT MD模拟的输入，模拟时间为100皮秒（100,000步），温度为2000 K。我们使用外推等级来过滤掉与训练数据差异较大的结构，该等级量化了通过ACE-MD采样的测试配置与训练数据之间的偏差。在高温MD模拟期间，如果模型遇到的结构的外推等级超过5，则将其纳入主动学习。使用D-最优性标准45确定收集的结构子集，以确保主动学习数据的多样性，并进行DFT单点计算以计算DFT能量和力。然后用主动学习数据对训练得到的势能进行额外1000次的重新训练，这个主动学习循环重复四次。总共，包括在主动学习期间获得的数据，训练数据集包含了4,628个结构，共包含492,216个原子。

需要强调的是，这里的目标并不是开发一个能够处理Ba–Ti–O化学空间中广泛任务的势能（如晶体松弛、性质计算等），这已经在更简单的系统中成功实现。具体来说，训练得到的势能的目标是在类非晶液体配置上获得与DFT精度相似的长时间NVT MD轨迹。出于这个目的，关于训练模型质量的基准测试在补充部分3中进行了讨论。

在线性不可逆热力学的框架中，由于驱动力在离子i中引起的扩散通量可以表示为
$$
{J}_{i}=-\mathop{\sum }\limits_{j}{L}_{ij}\nabla {\widetilde{\mu }}_{j}.
$$
(1)
这个表达式考虑了驱动力(\(\{\nabla {\widetilde{\mu }}_{j}\}\))对系统中所有离子的影响，以及离子i通过Onsager传输矩阵(L = [Lij])所经历的通量。然后使用Green–Kubo关系的微分形式从MD轨迹计算传输系数Lij：
$$
{L}_{ij}=\frac{1}{6{k}_{{\rm{B}}}TV}\mathop{{\lim}}\limits_{t\to \infty }\frac{\mathrm{d}}{{\mathrm{d}}t} < \mathop{\sum }\limits_{\alpha }\left[{{\bf{r}}}_{\alpha }^{{\rm{i}}}{(\rm{t})}-{{\bf{r}}}_{\alpha }^{{\rm{i}}}(0)\right]\cdot \mathop{\sum }\limits_{\beta }\left[{{\bf{r}}}_{\beta }^{{\rm{j}}}(\rm{t})-{{\bf{r}}}_{\beta }^{\,{\rm{j}}}(0)\right] > .
$$
(2)
这里V和T分别是系统的体积和温度；\({{\bf{r}}}_{\alpha }^{i}(t)-{{\bf{r}}}_{\alpha }^{i}(0)\)是时间t时物种i的第α个粒子的位移。自传输系数\({L}_{ii}^{{\rm{self}}}\)可以使用类似的关系计算：
$$
{L}_{ii}^{\mathrm{self}}=\frac{1}{6{k}_{{\rm{B}}}TV}\mathop{{\text{lim}}}\limits_{t\to \infty }\frac{{\rm{d}}}{{\rm{d}}t} < \mathop{\sum }\limits_{\alpha }{\left[{{\bf{r}}}_{{\rm{\alpha }}}^{{\rm{i}}}({\rm{t}})-{{\bf{r}}}_{{\rm{\alpha }}}^{{\rm{i}}}(0)\right]}^{2} > .
$$
离子i的自传输系数与Nernst–Einstein扩散率估计（也称为自扩散系数(Di)通过以下关系相关：
$$
{L}_{ii}^{{\rm{self}}}=\frac{{{D}}_{i}{c}_{i}}{{k}_{{\rm{B}}}T},
$$
(4)
其中ci是物种i的浓度。这个自项精确描述了稀薄系统的离子迁移率，其中扩散是理想的。实际上，Lij是时间相关函数\(< {\Sigma }_{\alpha }[{{\bf{r}}}_{\alpha }^{i}(t)-{{\bf{r}}}_{\alpha }^{i}(0)]\cdot {\Sigma }_{\beta }[{{\bf{r}}}_{\beta }^{j}(t)-{{\bf{r}}}_{\beta }^{j}(0)] >\)的线性拟合（扩散区域）的斜率。因此，在时间相关函数与时间的图中至少使用200皮秒（200,000步）的线性区域来计算所有传输系数。关于如何拟合传输系数以及在使用这一理论于固态反应时的其他注意事项的更多细节在补充部分4中提供。

由于离子运动的相关性，有效扩散系数的偏差可以通过L的非对角项（即交叉离子相关性）来研究，这些相关性通过
$$
{L}_{ii}^{{\rm{distinct}}}={L}_{ii} - {L}_{ii}^{{\rm{self}}.
$$
(5)
来量化。在这项工作中，我们关注了Ba–Ba和Ti–Ti阳离子对之间的相关性效应，通过相应的不同离子相关项来研究。为了比较不同相之间的相关项，我们计算了通过相关函数的自部分归一化的不同相关性，我们称之为相关因子(f)：
$$
{f}_{i}=\frac{{L}_{ii}^{{\rm{distinct}}}}{{L}_{ii}^{{\rm{self}}}}=\frac{{L}_{ii}}{{L}_{ii}^{{\rm{self}}}}-1.
$$
相关因子的值通常在-1和+1之间，接近0的值表示类型i的离子之间没有相关性。此外，我们通过计算每个离子在5 ?截止半径内的所有几何上可行的键来分析类液体相内的局部配位环境。这种方法捕获了第一和第二配位壳层内的空间相关性或键合相互作用。包括这两个壳层是因为阳离子通常通过桥接O离子进行相互作用。这种枚举是使用pymatgen46,47中实现的CrystalNN算法来执行的。

在训练ACE势能时生成的类液体结构被用作使用ACE势能进行NVT MD模拟的输入。这些模拟在750、1000、1250、1500和1750 K的温度下进行，时间步长为1飞秒，使用Langevin恒温器，摩擦因子为0.01飞秒^-1。为了获得更好的传输系数统计数据，使用基于MPMorph流程的ACE生产的五个帧作为NVT MD运行的起始配置，从而为每种组成进行五次重复的MD运行。

遵循参考文献48的方法，作为第一阶近似，我们将方程(1)中的驱动力\(\nabla {\widetilde{\mu }}_{j}\)（即电化学势梯度）近似为物种j的体化学势梯度?μj。这些梯度是通过使用pymatgen47和MP中的所有晶相数据构建的相关化学势图获得的（见补充图5）。我们通过化学势图上对应于前体（或一般而言，反应物）的域之间的最短距离来近似物种j的化学势梯度，除以产物层的厚度(h，见补充图8）：
$$
\nabla {\widetilde{\mu }}_{j}\approx \nabla {\mu }_{j}=\frac{\min \left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\beta }\right)}{h}.
$$
(6)
这种方法确保当相处于热力学平衡时，反应驱动力和产生的通量为零。有效的“扩散率常数”(KD)定义如下：
$$
{K}_{{\rm{D}}}=\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\frac{\left|{L}_{ij}^{\gamma }{V}^{\,\gamma }\times \min \left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\,\beta }\right)\right|}{{n}_{i}^{\gamma }{N}_{{\rm{A}}}.
$$
(7)
其中Vγ是产物γ的摩尔体积，\({n}_{i}^{\gamma }\)是γ中物种i的原子的分数，NA是阿伏伽德罗常数。我们假设产物的形成采用核壳模型（推导过程见补充部分6），γ的生长由Ba2+和Ti4+在界面上的传输支持。特别是，生长速率随着产物相的厚度而变化，并且随着反应的进行而减小。对于补充图8中所示的纯扩散限制情况，速率方程可以解为
$$
1-{(1-y)}^{1/3}=\frac{\sqrt{2{K}_{{\rm{D}}}t}}{{r}_{0}},
$$
(8)
其中y (1?≥?y?≥?0)是反应的完成程度，r0是前体粉末粒子的半径，这让人联想到通常用于拟合固态反应数据的Jander方程。在这项研究中，我们假设粉末粒子混合良好且大小均匀，允许我们使用\(\frac{\sqrt{{K}_{{\rm{D}}}}}{{r}_{0}}\)项来分析固态反应的动力学可行性，其中\({r}_{{\rm{B}}{\rm{a}}{\rm{O}}}={r}_{{\rm{T}}{\rm{i}}{{\rm{O}}}_{2}}={r}_{0}\)。在我们的模拟中，参数r0统一设置为10^-1微米，这代表了固态合成中常用的粒子大小。我们强调，这种公式并不适用于所有固态反应，因为当物种表现出不同的氧化状态，或者物种之间存在大量电荷转移时，\({\widetilde{\mu }}_{j}\approx {\mu }_{j}\)的近似并不成立，这两种情况在许多涉及过渡金属的固态反应中都会发生。在本研究中，所有考虑的相中的物种都只有一种氧化状态：Ba为+2，Ti为+4，O为-2。此外，通量的现象学表达式（方程(1)）是在质心参考系下定义的，这对不可压缩系统中系统的独立材料通量总数施加了限制。这些限制在补充部分4中有所描述。

最近开发了ReactCAA元胞自动机框架（ReactCA）19，用于模拟固态反应过程中粉末混合物的演变。反应容器被离散成体素（单元格），每个体素代表一个给定相的粉末粒子。在给定温度下，使用之前工作中开发的方法3，组合列举了系统中任何其他粉末（相）之间可能发生的所有反应，并使用评分函数进行评分。系统以两种相（在这种情况下，代表两种混合均匀的前体粉末）的随机分布进行初始化。在模拟的每一步中，选择两种反应物，并根据它们的评分从列举的反应中概率性地选择一个反应。然后，根据反应的化学计量比，将反应物概率性地替换为产物。这个过程在模拟过程中重复数百万次，以准确模拟反应容器内相的演变。在最初的出版物19中，相邻粒子之间的反应速率是基于反应热力学（通过MP计算出的零温度形成能以及机器学习估计的振动吉布斯自由能37）和机器学习估计的前体材料熔点50（作为反应动力学能力的初步指标）来估算的。在这项工作中，对于反应α + β → γ，我们使用一个新的评分函数S来纳入扩散通量，具体如下：
$$\begin{array}{l}\begin{array}{l}S={\sigma }_{1}\left(\frac{{K}_{\mathrm{D}}}{{r}_{0}^{2}s}\times \frac{\Delta {G}^{\ast }}{{k}_{\mathrm{B}}T}\right)\times {\sigma }_{2}\left(\frac{T}{{T}_{{\rm{m}},\mathrm{reactant}}}\right)\\ {\sigma }_{1}(x)=\frac{1}{3}\log [1+\exp (ax)]\\ {\sigma }_{2}(x)=\frac{1}{2}\log [1+\exp (bx-c)]\\ \Delta {G}^{\ast }=1+\mathrm{erf}(-d(\Delta {G}_{{\mathrm{rxn}}}+e))\\ {K}_{\mathrm{D}}=\mathop{\sum }\limits_{i}\mathop{\sum }\limits_{j}\frac{\left|{L}_{ij}^{\gamma }{V}^{\,\gamma }\times {\mathrm{min}}\left(\,{\mu }_{j}^{\alpha }-{\mu }_{j}^{\,\beta }\right)\right|}{{n}_{i}^{\gamma }{N}_{\mathrm{A}}}\end{array}\end{array}$$ （9）
这个评分函数的形式是对ReactCA19中现有评分函数的扩展，使用了相同的常数a、b、c、d、e。因子s对所有考虑的反应都设置为相同的值，以便在CA模拟中使评分结果稳定，并确保前体之间的反应在BaCO3的分解温度以上开始（s = 5 × 10^13）。该评分函数对Tammann规则20具有“柔和”的激活作用，这是通过函数的σ2部分实现的；而对于具有负反应自由能（ΔGrxn）的反应，则通过误差函数实现“开启”行为。有关评分函数这些方面的更多细节，请参考参考文献19。在这项工作中，我们在类液态产物相（KD）中实现了显式的定量传输扩散常数，通过σ1函数实现柔和激活，该函数也是一个软加函数。这个函数还考虑了修正后的ΔGrxn和温度。因此，在Tammann温度以下（\(\frac{T}{{T}_{{\rm{m}},{\rm{r}}{\rm{e}}{\rm{a}}{\rm{c}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{t}}}}\le 0.67\)），评分结果是非零的，但仅略为正值。在Tammann温度以上，具有略微正的ΔGrxn或低离子通量的反应会被赋予非零但较小的评分。随着反应温度的升高，函数的σ2部分会增加，从而提高反应速率。同样，温度的升高也会增加扩散速率常数KD（通常是超线性方式），这通常会增强评分函数σ1部分的输出。然而，正如这项工作所展示的，随着温度的升高，KD会接近一个饱和点，这反过来会导致σ1的输出在非常高的温度下趋于平稳甚至下降。这种效应在反应温度接近或高于熔点时导致了相选择性的热力学行为。最后，评分函数与旧的评分函数保持向后兼容性，这在系统中的某些或所有相的动力学传输参数不可用时尤为重要。

DFT计算使用了Vienna ab initio模拟包51、52和Perdew–Burke–Ernzerhof53的广义梯度近似公式以及投影增强波势54、55来进行单点DFT计算和AIMD计算。此外，为了最小化计算需求，所有AIMD计算仅使用了Γ点，并且是非自旋极化的。有关所使用参数的更多细节，请参考atomate2中的MPMorphMDSet类。为了构建高温相图，使用了一个机器学习估计器来估算所有相的有限温度振动吉布斯自由能，该估计器用于估计振动对熵的贡献37。