《Journal of Chemical Information and Modeling》：Rapid Generation of Transition-State Conformer Ensembles via Constrained Distance Geometry
编辑推荐：
本综述介绍了一种名为racerTS的新方法，该方法利用经过修改的距离几何算法，能够快速生成过渡态（TS）构象集合。与现有方法（如CREST和GOAT）相比，racerTS在保持相当的构象空间覆盖率和更高过渡态结构有效性的同时，显著降低了计算成本（最高可达500倍），为高通量计算化学筛选和构建高质量的机器学习（ML）催化剂数据集提供了高效工具，有望推动可持续催化剂的发现。
在计算化学中，准确建模化学反应需要充分考虑过渡态（TS）的构象集合，这在计算催化剂设计中起着关键作用。虽然CREST和GOAT是当前主流的TS构象集合生成方法，但其相关的计算成本仍然是计算化学流程中的一个主要瓶颈，尤其是在为催化剂设计生成大型机器学习数据集时。为此，本研究提出了racerTS（RApid Conformer Ensembles with RDKit for Transition States），这是一种基于约束距离几何学的高效TS构象集合生成方法。
引言
开发新型、更可持续化学反应的催化剂是一个繁琐的过程，可以通过计算来加速。包括活化能（E?）和反应选择性在内的多种反应性质都由过渡态决定，因此需要对它们的结构和能量进行精确建模。而精确建模的前提是计算目标化合物的构象集合。Laplaza等人最近的研究展示了过渡态构象集合的不同处理方式甚至可能逆转预测的反应选择性，凸显了其在TS建模中的重要性。同时，在使用机器学习模型预测反应性质时，也越来越需要考虑过渡态的构象集合。例如，Zhu等人指出，在预测Ru催化的烯酰胺不对称氢化反应的对映体过量时，考虑TS构象集合中的信息与仅考虑最低能量构象相比，可以降低20%的预测误差。这些研究清楚地表明，从TS构象集合中获取信息对于实现更准确的计算预测至关重要，而这是将计算工具成功整合到催化剂开发流程中以发现可持续化学过程所必需的。
对于基态结构，构象集合已越来越多地被包含在广泛用于机器学习的可用数据集中。对于基态结构的构象集合计算，有多种开源程序可用，例如CREST、全局优化算法（GOAT）、RDKit、Balloon、机器学习方法以及商业软件等。这些程序在准确性、构象空间覆盖率和计算成本方面提供了不同的权衡。其中，CREST在需要精确构象集合时常用，但其基于元动力学的物理模拟方法导致其计算成本较高，在许多情况下，构象集合生成是计算化学流程中的瓶颈。基于距离几何算法的RDKit ETKDG方法是基于物理模拟的快速替代方案，它通过从距离边界矩阵中随机采样原子间距离，显著降低了计算成本。此外，ETKDG算法中还利用额外的化学知识和实验数据对构象进行细化，生成的构象通常通过使用快速力场方法进行几何优化来后处理。
与基态相比，可用于过渡态的构象生成方法较少。最近，机器学习已被应用于过渡态生成任务，包括Kim等人提出的TSDiff模型。通过对相同的反应物和产物重复采样模型，可以生成多个TS构象。其他模型还包括Heid等人的模型以及Duan、Liu、Du等人提出的React-OT模型，这些模型可以从不同的反应物和产物构象对开始生成TS构象集合。尽管这些进展很有前景，但所有机器学习方法都固有地受限于其训练数据。目前，大多数TS结构生成模型都是在原子数不超过11个重原子的数据集上训练的，这限制了它们的应用范围。大多数数据集也只包含有机结构。目前，只有Kraka等人考虑了涉及过渡金属的反应。量子化学和化学信息学基础的构象生成器目前具有更广泛的应用范围，并将在为生成模型生成更多样化和更准确的训练数据方面发挥关键作用。缺乏此类方法的原因在于，为过渡态定义一个具有合适原子类型的有效分子拓扑结构仍然具有挑战性，这是所谓过渡态力场领域的一个活跃研究领域。特别是反应中心，即反应中涉及键断裂和形成的原子和键，对构象生成提出了重大挑战。由于TS的原子间距离和角度通常在基态结构中观察不到，力场对这些键/角的准确性较低。它们也没有被整合到基于知识的算法中，且机器学习生成器在训练期间没有见过它们。TS构象集合生成器通过采样具有约束原子位置、距离和/或角度的构象来应对这一挑战，这些约束通常取自用户提供的TS模板。
方法
racerTS的设计旨在与CREST和GOAT类似地工作，用户需要提供一个TS猜测的3D结构（即xyz坐标文件）、总电荷以及反应中心原子的索引列表。用户还可以选择性地提供反应物或产物的SMILES字符串。基于此输入，racerTS生成一个包含输入过渡态多个构象的xyz坐标文件。
工作流程首先使用RDKit的rdDetermineBonds模块，根据用户提供的3D结构和总电荷（以及可能的SMILES）生成一个RDKit Mol对象。该对象包含键类型和形式电荷信息，这些信息要么直接从3D结构推断得出，要么取自用户提供的产物或反应物的SMILES字符串。为反应中心分配原子和键类型通常是容易出错的，因为反应中心不显示典型的键和价态。对于这些原子和键，程序接受任何有效的RDKit分配，因为在构象生成和精修过程中，这些原子和键是受到约束的，所以正确的标记并非严格必需。如果无法分配键序，racerTS默认仅推断连接性，然后进行清理。
接下来，使用获得的Mol对象生成初始TS构象。与RDKit用于基态构象生成的方法类似，racerTS依赖高效的ETKDG方法进行快速TS构象生成。为了获得合理的构象，ETKDG在距离几何方法中考虑了实验扭转角和几何启发式规则。对于过渡态，为了保持关键的反应中心，racerTS约束了相应原子的位置。此外，在开发过程中发现，约束反应中心原子及其所有邻居原子之间的距离可以得到更可靠的TS构象。反应中心及其直接邻居的统一在下文中被称为冻结原子。随后，通过重复从受约束的距离边界矩阵中随机采样距离，然后按照ETKDG方法进行启发式结构细化，从而生成多个构象。
初始生成后，获得的构象需要进一步的结构优化。默认情况下，所有生成构象的能量都使用MMFF94力场的RDKit实现进行最小化，该力场针对有机分子进行了参数化。一种失败模式是无法识别的原子类型，这可能特别与反应中心（例如，分配给反应中心的自由基）和有机金属反应有关。作为备用方案，则使用RDKit的UFF实现进行优化，所有无法识别的原子类型都被分配默认参数。
最后一步是对构象集合进行修剪以去除重复项，同时考虑能量和结构因素。与CREST和GOAT使用的结构准则类似，如果某个构象与某个较低能量构象的RMSD < 0.125 ?，则将其移除。与CREST不同的是，racerTS默认仅考虑重原子进行RMSD计算，侧重于移除高度相似的构象。为了仅考虑相关能量窗口内的构象，CREST会修剪所有能量比最低能量构象高6 kcal/mol的构象。由于力场衍生能量的不准确性，并为了避免移除相关构象，racerTS默认应用了基于力场能量的更保守的20 kcal/mol能量窗口。在此步骤中未被移除的所有优化后的构象将连接成一个单一的多结构xyz文件作为输出。
性能基准测试
为了评估racerTS生成的构象集合的质量，我们将其与最先进的方法CREST和GOAT进行了比较。我们评估了理想构象生成器的五个标准。为了评估racerTS在高通量计算化学工作流中的可用性，我们从文献中选取了nrxn= 20个反应的过渡态来生成TS构象集合，其中包括具有挑战性的有机金属配合物。这20个反应涵盖了多种化学反应类型，确保了能够准确展现racerTS的优势和劣势。在正文中，我们比较了五种不同的方法：使用3D结构作为输入的racerTS、使用3D结构和反应物SMILES作为输入的racerTSSMILES、使用3D结构作为输入且后处理级别降低的racerTS(no xTB)、在GFN2-xTB//GFN-FF级别运行并约束冻结原子位置和距离的CREST、以及在GFN2-xTB级别运行的GOAT。为了增加racerTS和racerTSSMILES方法中构象能量排序的准确性，所有构象都使用xtb程序在GFN-FF力场下进行后优化，并精确固定冻结原子，最终能量使用GFN2-xTB进行单点计算评估。后优化后，再次以相对于最低能量构象6 kcal/mol的能量窗口修剪高能量构象，并以重原子RMSD < 0.125 ?的阈值修剪潜在的重复构象。racerTS(no xTB)则不进行此类后优化。
为了比较生成的构象集合，我们对每个集合应用了一个常规的高通量计算化学工作流程。作为该工作流程的第一步，我们为每种方法提取五个最低能量的构象。此外，我们还使用marc程序从每个集合中提取最多10个代表性构象。在整个过程中，提取的构象将分别称为低能量构象和marc构象。然后，使用ORCA 6.0.1，在r2SCAN-3c理论水平上，对每个提取的构象进行DFT计算。在此步骤中，首先在固定反应原子之间距离的情况下进行约束几何优化，然后是过渡态优化，最后进行快速反应坐标计算以确保过渡态连接了所需的反应物和产物。
评估指标
我们根据上述理想构象生成器的五个属性来分析所描述的TS构象生成方法的性能。
计算成本：为了分析方法i的计算成本，我们考虑每种方法的挂钟时间，并将其相对于racerTS(no xTB)的时间进行归一化。总体而言，我们发现使用racerTS相比CREST可以提供4倍的加速，相比GOAT则高达500倍。这种加速是可以预期的，这归功于距离几何算法公认的计算效率。与GOAT的巨大速度差异也可以理解，因为在撰写本文时，GOAT只能在GFN2-xTB级别运行，这比GFN2-xTB//GFN-FF慢得多，但也产生更高质量的结构。相比之下，由于CREST、racerTS和racerTSSMILES都在GFN2-xTB//GFN-FF级别生成结构，因此racerTS和racerTSSMILES的效率提升证明了距离几何方法在采样速度上优于元动力学。值得注意的是，racerTS/racerTSSMILES与racerTS(no xTB)之间的时间差异完全归因于racerTS/racerTSSMILES内部的GFN2-xTB//GFN-FF后优化，这大约占总运行时间的四分之三。虽然racerTS(no xTB)内部较低的理论水平导致所需的挂钟时间显著减少，但构象的结构和质量较差。后续对低能量区域准确性的分析表明，需要后优化来提供更准确的构象集合。尽管如此，如果构象排序对于应用来说不是决定性的（例如，当整个集合无论如何都要在更高水平上进行优化时），racerTS提供的加速可能更加显著。因此，我们的结果表明，racerTS实现了所需的显著加速，与最先进的方法相比，计算成本降低了数个数量级。
穷举性：为了研究每种方法在探索构象空间方面的穷举性，我们首先查看每个反应生成的构象数量。总体趋势是，CREST产生的构象数量最多，在20个基准测试反应中有16个如此。这一趋势与早期的研究一致。我们假设这一观察结果可以部分解释为冻结原子固定较不严格（与racerTS和GOAT相比，CREST不允许精确固定，只允许约束）以及RMSD修剪较不严格（CREST在RMSD修剪中还考虑氢原子，而racerTS不考虑）。相反，我们发现没有一种方法能始终产生最少的构象，其中racerTS有11个反应如此，racerTSSMILES有4个，GOAT有6个。比较racerTS和racerTS(no xTB)产生的构象数量，可以清楚地看到，GFN2-xTB//GFN-FF后优化以及随后更严格的能量过滤移除了许多高能量构象。平均而言，只有41%的来自racerTS(no xTB)的构象保留在后优化的racerTS集合中。
虽然构象数量可以指示所探索的构象空间，但该指标并不能直接反映多样性，因为即使在RMSD修剪后仍可能发现许多高度相似的构象。为了研究构象空间的覆盖情况，我们计算了每种研究方法的空间探索指标。结果显示，GOAT在构象空间探索方面提供了最佳能力，平均F1得分为0.93，而所有三种racerTS变体和CREST表现相似，四种方法的平均F1分数在0.74到0.80之间。F1分数提供了方法精度和召回率之间的权衡，只有当两个组成部分都很高时，F1分数才高。对每种方法的这些指标进行更仔细的研究发现，除了racerTS(no xTB)外，所有方法的精度都接近1，这表明这些方法产生的几乎所有构象都属于完整集合。唯一的例外是racerTS(no xTB)，其精度显著较低，为0.84。这一结果表明，用racerTS(no xTB)生成的一些构象不属于完整集合，该完整集合是通过考虑优化和第二次更严格的能量修剪获得的。因此，这一结果是预期的，因为racerTS(no xTB)没有经过具有更严格能量阈值的第二次修剪。未通过此修剪移除的高能量构象通常不属于组合集合的任何聚类，从而降低了观察到的精度。
对于大多数方法，精度接近1，因此召回率是限制空间探索的因素。事实上，对于所有racerTS变体以及CREST，这一指标介于0.65和0.72之间，显著低于1。同时，GOAT优越的空间探索性能可以用其0.82的优越召回率来解释，远远优于其他方法。CREST较低的召回率起初令人惊讶，因为CREST在绝大多数反应中产生的构象数量最多。然而，这一观察表明，基于元动力学的CREST产生了许多相似的构象，这些构象由于反应中心受到扰动而未被修剪。这一观察再次表明，构象数量不一定与所探索的构象空间相关。相比之下，racerTS变体较低的召回率可以解释为它们产生的构象数量较少，导致组合集合中的一些聚类未被填充。特别是，racerTS在11/20的反应中构象数量最少，召回率也最低，为0.65。在racerTS内部，初始生成的构象数量可以通过传递参数conf_factor来轻松调整。初始生成更多的构象会导致发现更多的聚类，从而也提高F1分数。在racerTS的开发过程中，我们发现30是conf_factor的良好默认值，因为它与最先进的方法CREST的性能相匹配。然而，在生成所有相关构象至关重要的应用中，可以轻松修改racerTS以实现更高的空间探索，尽管计算成本会增加。
为了进一步研究每种方法产生的构象空间，我们研究了构象集合中构象与输入3D几何结构之间RMSD的分布。将该分布与为组合集合计算的RMSD分布进行比较。与空间探索指标的结果一致，GOAT最完整地覆盖了空间，因为观察到的RMSD分布与组合集合的分布最为相似。这由平均空间分布得分0.19所证明。再次确认了空间探索指标的结果，racerTS变体显示出与CREST相似的JS散度。这一结果再次证实了构象数量不一定与探索的构象空间相关；由于CREST产生许多高度相似的构象，其RMSD分布与组合集合的分布不同，在组合集合中，这些构象大多已被修剪。
总体而言，我们的结果表明，GOAT是穷举采样构象空间最熟练的方法。虽然racerTS的表现与CREST相当，但我们也展示了如何轻松修改racerTS的设置以更穷举地采样构象空间，尽管计算成本会增加。
有效性：除了穷举探索构象空间外，构象集合生成器还应提供有效的TS结构。为了评估这一点，我们测量了在通过上述优化流程后，所选构象中有多少能在DFT水平上提供有效的过渡态。比较所有五种方法，我们发现所有racerTS变体产生有效TS的比例都高于CREST和GOAT。所有方法都以较高的成功率产生有效的TS，成功率在82%到90%之间。这些失败率在高通量计算流程中是典型的，通常会遇到10%到50%的错误率。racerTS产生更多有效TS的观察结果也体现在至少识别出一个有效TS的反应数量上。所有racerTS变体对所有20个测试反应都产生了至少一个有效TS，而CREST和GOAT分别只对19个和17个反应如此。CREST未能为SN2糖反应产生有效的TS构象，而GOAT则分别未能为氢硼化反应、NHC反应和Ru烯烃复分解反应产生有效TS。就GOAT而言，氢硼化反应和NHC反应失败是由于计算流程中的几何收敛问题，而Ru烯烃复分解反应识别出的TS在进行QRC后导致了不同反应的反应物和产物。使用CREST作为构象生成器时，唯一没有有效TS的反应是SN2糖反应，其中识别出的TS在进行QRC后导致相同的反应物和产物。与CREST和GOAT的失败相反，racerTS能够为每个测试反应产生有效的TS，超越了最先进的方法。racerTS在所有情况下生成有效TS的能力证明了其在广泛多样的应用中的适用性，包括高度复杂、柔性和有机金属体系。
低能量构象的准确性：虽然由于计算效率的考虑，在整个计算化学中越来越普遍地考虑整个构象集合，但活化能通常仅通过考虑过渡态的最低能量构象来计算。由于低能量构象由于具有最高的玻尔兹曼权重而对整体能量贡献最大，因此正确识别低能量构象对于每个构象集合生成器都高度相关。
我们首先确定任何方法在DFT水平上找到的总体最低能量构象。然后，我们评估每种方法提取的低能量构象是否提供了一个与最低能量构象相同的构象。我们考虑了前1和前5准确性，即最低能量构象和五个最低能量构象是否包含一个与总体最低能量构象相同的构象。总体而言，识别出最低能量构象的比例出人意料地低，racerTS显示出最高的前1准确性，为25%。其他方法分别在三个或四个反应中正确识别出最低能量构象。对于racerTS，正确识别出前1低能量构象的反应是环化、环氧化物重排、周环反应、炔丙基化和SNAr反应。考虑到前5准确性适度地提高了大多数方法识别最低能量构象的比例，最高达到35%。我们假设这种相对较低的检索率是由于在半经验水平或力场水平上构象排序不准确造成的。对于5个反应，总体最低能量构象仅在marc构象中发现，而不在lowe构象中，即被DFT排序为低能量，但未被半经验方法排序。对于其他75%的反应，识别出的最低能量构象分布在各个构象生成器中，racerTS的表现与最先进水平相当。考虑到所有方法的高构象空间穷举性，可以假设实际上生成了更多在DFT水平上能量较低的构象，但未被半经验方法如此排序。这种对关键构象的忽略在其他工作中也被识别出来。这一假设得到了racerTS(no xTB)较低准确性的支持，但只有通过在DFT水平上优化每个构象才能确认，这超出了本研究的范围。因此，最大的准确性提升可以通过在DFT水平或其他更新的半经验方法或机器学习势上优化和评估每个构象来实现。我们的结果证实了racerTS作为一种高效的构象筛选器，其准确性达到或超过了当前的最先进水平。
为了进一步探索构象生成器对于低能量构象的准确性，我们研究了计算活化能时的偏差。首先，我们研究了当选择每种方法识别出的最低能量构象时，与总体最低能量构象相比的活化能偏差。其次，我们还研究了每种方法选定的marc构象的活化能与总体marc集合的活化能偏差。比较构象集合生成器之间的ΔElowe?，我们看到所有racerTS变体都能实现准确的活化能计算。对于这些方法中的每一种，ΔElowe?的中位数显著低于化学精度。对于除racerTS(no xTB)外的所有方法，ΔElowe?< 0.18 kcal/mol，展示了racerTS对于活化能计算的适用性。racerTS(no xTB)相对较高的中位偏差是由于通过MMFF或UFF进行能量排序，而其他构象生成器是在GFN2-xTB//GFN-FF级别排序的。不准确的排序导致能量较高的构象被估计为低能量，从而导致较高的ΔElowe?。尽管如此，所有方法的ΔElowe?中位数都远低于化学精度，这表明每种方法在“平均”用例中都适合寻找低能量构象，尽管个别感兴趣的反应的偏差可能高于平均值。图中更宽的分布显示，我们发现racerTS变体在特殊情况下更容易产生高能量异常值。作为显著的异常值（ΔElowe?> 5 kcal/mol），我们在多个racerTS变体中观察到了酰胺甲基化、环氧化、肽、SN2糖和Pd氧化加成反应。对这些反应进行更深入的研究发现，导致高能量TS构象的两个主要因素是：相对于酰胺键的不利旋转异构体导致TS的空间位阻冲突，以及倾向于为某些六元环生成船式构象。增加采样构象的数量可能在一定程度上缓解这些问题，但会增加计算成本。CREST在环氧化反应中也显示出一个异常值，其中稠环的构象空间探索过于刚性，因此未能识别出最低构象。
虽然基于最低能量构象计算活化能仍然被广泛采用，但更复杂的方法（如marc）会考虑代表玻尔兹曼分布的多个构象，以获得更准确的活化能。为了模拟这种用例，我们计算了基于marc选定构象的各种构象生成器的活化能。我们的观察表明，对于所有构象生成器，ΔEmarc?的中位偏差≤ 0.1 kcal/mol，这再次证明了racerTS在日常用例中的实用性。与对ΔElowe?的研究类似，我们观察到当考虑ΔEmarc?时，racerTS变体也更容易出现异常值。在racerTS变体中频繁出现的高偏差异常值与ΔElowe?相同，原因如上所述。CREST的异常值与上面讨论的一致。
总体而言，我们观察到最高的异常值出现在racerTS(no xTB)中，这很可能是由于MMFF/UFF与GFN2-xTB//GFN-FF相比能量排序性能较差。因此，在容易出现上述异常值的情况下，需要谨慎应用racerTS。尽管如此，我们证明了racerTS生成的构象集合能够在大多数用例中实现准确的活化能估计。这一观察结果展示了racerTS在计算化学工作流程中作为构象集合生成器的实用性，既可用于基于最低能量构象的活化能计算，也可用于更复杂的选择策略。结合我们之前的分析，结果表明，racerTS以远低于最先进方法的计算成本，提供了相当准确的过渡态构象集合、广泛的构象空间覆盖率和较高的有效性。
结论与展望
本文介绍了racerTS，一种用于快速生成过渡态构象集合的方法。racerTS使用了基于距离几何学的构象生成方法的修改版本，实现了约束距离采样，这是保持过渡态中反应中心原子所必需的。其唯一的依赖是流行的化学信息学库RDKit，这使得它可以轻松集成到现有的化学信息学和计算化学工作流程中。racerTS的使用灵感来源于最先进的过渡态构象集合生成器CREST和GOAT，它们只需要一个过渡态的3D几何结构作为输入，以及关于哪些原子属于反应中心的信息。通过为20个涉及四个有机金属反应的反应生成TS构象集合，我们将racerTS与C
