“分而治之”方法构建复杂分子体系高精度势能面及应用_宋国梁.pdf

文章编号：（）收稿日期：基金项目：国家自然科学基金（，）作者简介：宋国梁（），博士，；李振华，教授，通信作者，：“分而治之”方法构建复杂分子体系高精度势能面及应用宋国梁，李振华（复旦大学 化学系，上海 ）摘要：“分而治之”是处理复杂分子势能面的一种快速简洁方案。本文介绍了近年来本课题组在这一领域的方法研究进展及相关应用研究。在新的组装势能面算法和数据库支持下，我们将复杂大分子势能面计算量降低到（）次方的标度，并且计算精度接近 甚至 （）的水平。采用该方法，我们对两个不同的应用场景进行了测试。结果表明，高性能的组装势能面可以完成以往难以实现的大规模的分子结构搜索、稳定性分析以及新型分子结构设计等关键性任务。关键词：分而治之；复杂分子体系；势能函数；势能面中图分类号：文献标志码：从头算方法可以提供媲美实验精度的分子结构、能量、性质以及反应通道的预测，不过其高昂的计算代价（）使得其只能应用于中小分子体系。对于复杂的大分子体系，通常采用分子力场或者半经验方法进行计算，这些方法在速度上可以满足模拟要求但精度上仍缺乏普适性。严格的从头算方法在大体系中通常应用于 方法中，用于解决那些有明确反应中心的体系的大规模模拟问题。“分而治之”（，）方法是另外一种处理方案。该方案中，大体系将被切割成小的区域，每个区域及其周围的虚拟环境进行独立计算，最后再进行整体拼装。这一方案主要针对那些反应中心不是特别明确或者根本无法区分反应中心的分子体系。年以来，类方法的很多细分方法如 ，等都在复杂的生物和材料领域取得了相当多的进展。从理论上看，方法将大体系切分为小的子体系（图），这必然导致部分价键被切断，从而引发波函数在边界处的变异并进而影响电子密度和能量。为减小这种影响，最佳切分位置尽量选择在相对局域化图大分子体系切分成小的子体系及能量计算公式 第 卷第期 年月复 旦 学 报（自然科学版）（）DOI:10.15943/ki.fdxb-jns.20230210.001的单键处。因此选取高效的切分策略是 类方法需要面对的第一个重要问题。进一步，在切断处，通常需要采用补 或者补充基团的方式使得子体系计算时的边界尽量接近原本大体系的化学环境。不过两者间始终存在一些关键参数如电荷、能量、力的差异。如何将这一差异最小化是另外一个关键点。最后，在化学反应过程中，由于各种子体系间存在原子的迁移，键连关系同样在不断变化，如何确保势能面的连续性同样是一个严峻的问题。本文第、和小节分别介绍了本课题组发展的基于能量的 方法、子体系区域划分策略、多体展开式断键恢复方法以及多体区域平滑过渡方案。第和小节介绍了数据库和分子生成的一些技术细节。第小节讨论了本方法的速度和精度。最后通过个应用例子测试了本方法的有效性。基于能量分解的“分而治之”方法（）一个大分子体系的势能可以表达为一个低精度方法的势能 叠加一个高精度势能和低精度势能 之差。该差可以利用多体展开，进一步近似表达为各个子体系的高精度势能和低精度势能差值的和：（）。（）对能量差势能面进行子体系切分比直接对总能量进行切分具有更好的可移植性。这是因为低精度方法也是薛定谔方程的一个近似解，残余的误差主要来自于对交换相关作用描述得不够精确。而后者通常认为是局域的，具有很好的可切分特性。其次，在断键处理上，虽然依然需要采用补充基团并进行二体或者三体逐步校正，但低精度方法已经提供了接近正确的价键和电子密度描述，因此断键校正也相对容易。最后，高低精度方法能量差的变化幅度比总能量变化幅度要小个数量级以上，这在提高势能面拟合精度时也是一个巨大的优势。本文中所提及的势能面若非特别标示均是指高低精度方法之差的势能面。子体系区域划分方法对于一个稳定的复杂有机大分子体系，可以采用价键方法对该分子每个价键进行归类。然后再根据价键特征确定最优的截断方案。通常认为，单键是较好的电子局域化键。切分位置恰好处于单键时不会大范围破坏子体系中其它离域区域的电子结构特征，同时通过补 或其它基团也更加容易恢复截断处价键的化学环境。在可截断的单键的判断过程中要注意两个问题：首先是避免选择可以共振形成双键的单键，例如苯环的根共轭键；其次，还要尽量避开两个共轭双键之间的单键。为此，我们首先需要对分子进行自动的共振结构分析。共振结构是 年鲍林（）创立的一种分子结构理论。当一个分子、离子或自由基的结构不能用路易斯结构式正确地描述时，可以用多个路易斯式表示，这些路易斯式称为共振结构（），又称极限式或正则结构。共振结构根据八隅体规则书写，共振结构式之间只允许键和电子的移动，而不允许原子核位置的改变。所有的共振结构式必须具有相同数目的未成对电子。在共振分析中，首先根据所有键长是否小于平衡键长（均以单键平衡键长为准）的 倍来判断两个原子是否相连，这样我们就可以获得整个分子的拓扑连接图。然后根据每个原子的连接数确定该原子的杂化类型，再根据两个原子的杂化类型确定其连接的价键类型。该过程需要反复递归执行，并且由于过程中会对不确定的双键引入双键假设，因此会形成双键认定分歧以致形成多种不同的共振结构。首先根据平衡键长（均以单键键长为准）倍判断两个原子间是否连接成键。再针对不同原子周围的连接数目确定价键类型：）单键、双键和三键的直接确定：原子周围只有一个连接，则该键是单键；原子周围有个连接，则这两个键是单键；类似的，的根单键和的根单键都可以通过该方式确定下来。其他确定的双键如=和三键等等也可以在这一步确认下来。部分特殊的复杂结构也可以在这一步预先判定，例如，等。第期宋国梁等：“分而治之”方法构建复杂分子体系高精度势能面及应用）单键、双键和三键的推断确定：在上一步完全确定的键的基础上，推断其他部分可以确定的键：如连接的原子已经确认两个键是单键，剩余的一个连接对象是 原子，且该 原子仅一个连接，则该 连接是双键。通过类似推断可以确定体系中的所有可推断确定的键。图苯乙烯分子中苯环的共振结构 ）对于依然没有确定连接，采用树结构遍历其可能的键的类型：即：对于每个不确定连接产生一个树的分支，每个分支假定该键为一种键的类型，然后重复的推断过程。如果依然有不确定连接，则继续产生分支假定并继续的推断，直到所有连接都已确认则得到一种共振价键类型。遍历树的所有分支则可以获得所有共振结构。以苯乙烯为例，在前两步中只能确定乙烯基上所有价键，但无法确认苯环上的键类型。此时需要对 进行分支假定，两个分支分别为单键和双键，在两个分支上分别继续进行）步骤的价键推断，可以获得图中两个共振结构。程序中预设了平衡键长（表），对于其他元素，针对具体体系需要建立简单模型优化后确定相应的键长。表 的预设平衡键长 键类型键类型键类型单键 单键=双键 单键 单键=双键 单键 单键=双键 单键=双键 叁键 单键=双键 叁键 单键=双键 叁键 遍历并记录所有可能的共振结构，并给予每个共振结构相应的权重：（），（），（）。（）其中，是该共振结构中某个键的键长；是平衡键长（详见）；是平滑过渡函数，确保其值在从 到 时从 平滑过渡到 。这也表示分子只有开始脱离平衡键长至少 倍后才会发生势能面的平滑过渡行为，在平衡结构附近（）依然会采用原本的势能面。这个平滑过渡的权函数的存在可以很好的处理整个势能面的平滑问题。当存在多种共振结构时（），需要对每种共振结构的权重进行归一化处理。归一化处理后可以方便后面多种不同类型势能面拼接时进行简单加权处理：（）。（）共振分析后每个键在每种共振结构下都有明确的价键类型，可以根据数据库中预存的平衡键长确定其相对键长。当键长超过 倍平衡键长时就要开启该共振结构逐步平滑过渡到的过程。参见公式（）和（）。本文中引入的共振分析主要是对 倍以外的势能面区域起作用，因此平衡键长的少许差别不严重影响平滑过渡的历程。当然，对平滑过渡区域的能量是有影响的，尤其是大体系切分后反应过渡态恰好在切分位置附近时，这将使得过渡态结构附近能量的误差相对较大。复 旦 学 报（自然科学版）第 卷表部分键的类别及相应的切割优先度 类别键类别键 =经过共振分析，我们精确确定了每个连接的键的类型，这时可以开始对大分子进行切割。由于切分方法有很多种，我们对每个连接引入了一个切分时的优先度，越小越优先切割，这样我们在切分时可以按照一定顺序进行，并获得相对最优的切分方案。我们将经过共振分析的所有价键划分为 类。类键是纯粹的单键（至少一个端点是 杂化原子），类键是单键但其两端原子有 杂化原子，类键是单键双键共振键（即在部分共振结构中是单键，部分共振结构中是双键），类键是确定性的双键，类键是确定性的三键及以三键以上的键。每个键还定义了类别内的排序值（见表）。整体类别排序按照优先度逐步升高。、和类键由于有共轭键的存在因此尽量不要切断，如果必须切断，其优先度排序也非常靠后。最优先的是类键。每种键的优先度值通常依据键的相对极化强弱由人工预先指定，无极性键通常排在最前面，详见表。确定了每个键是否可截断以及截断优先度后我们就可以将大分子进行切分处理。以 为例（图），切分中要遵循以下规则：）首先根据分子图的拓扑结构，收缩所有的叶子节点。）剩余的图结构再进行单连接截断（单连接：图的两个节点间仅有一条连接通道）。图示例大分子体系（）的一种切分为多个子体系的方案 （）重复以上两步，直至最终剩余结构的每个节点至少包含两条边，每条边都不是单连接边。这种结构被称之为最简拓扑分子结构。最简拓扑结构通常是复杂网状连接分子，将这些结构从整体分子结构中切分出来，剩余的子体系通常是链状分子，我们需要对其进行进一步切分，使得每个切分出的子体系的非 原子数目小于等于。）最后还需要进行一次合并，使得那些在原分子中是直接相连且两个体系非 原子数目小于等于第期宋国梁等：“分而治之”方法构建复杂分子体系高精度势能面及应用的子体系重新合并。或者两个子体系合并后恰好是势函数数据库中存储过的体系。）以上过程中切分位置仅限于 类键，如果无法找到切分方案，则依次放松限制，允许对 类键进行截断，直到找到切分方法后停止。经过以上处理，通常得到切分方案数目依然很多，我们采用如下方案确定最优切分方案：）优先寻找切分出的子体系数目最少的方案。）切分数目相同时，优先选择最小的子体系的原子数最接近的方案。）依然无法分辨时，优先选择的累计值最低的切分方案。断键复原切断的价键通常要补充 原子或者其它补充基团。=例如对于双键的切断通常要补充。不过这一补偿仅仅能够使得子体系（单体）的薛定谔方程的求解相对接近其在整体分子中的真实值，我们还需要进行多体校正使得分子总能量和电荷能有效逼近整体计算值。根据我们之前的 方法，通过二体和三体电荷和能量外推，就能使能量充分逼近整体计算值。在 方法中，共振结构展开到三体的能量和电荷通过如下公式得到：，（），（），（），（）。（）其中，是单体单独计算的能量。更精确的总能量通过如下公式外推得到：（）。（）这里，是归一化电荷净值，是电荷能量校正系数，详细算法可以见文献 。在本工作中用高精度和低精度能量差来替换，而使用低精度半经验方法（如 ，）的全体系电荷来作为参照基准。要注意的是，为了确保长程效应不丢失，全部二体都需要被计算，而三体体系则仅需计算那些有直接断键的相连三体。在一些分子中，部分二体和三体的体系可能会比较大，采用高精度计算方法计算这些三体的开销较大，这会成为整体计算中的瓶颈。这时也可以使用较低精度的方法例如用 方法计算这些二体和三体的多体校正值。由于这部分的影响主要是断键的外推补偿，因此采用低精度方法代替而引入的相对误差不算太大。平滑过渡在大分子体系化学反应的过程中，部分原子可能从一个子体系逐步迁移到另外一个子体系。在这一过程中，键连关系在不断发生变化中间的断键结构包含着向多种键连关系转移的可能。因此，每一步结构的变化都需要对整个结构重新进行共振结构分析。在分析中会产生新的共振结构（）或者删除已有的共振结构（），并更新每个现有共振结构的权重。通过这些平滑过渡的共振权重我们可以将不同的子体系及其二体三体校正势能面加权组合成最终完全平滑的整体势能面。这种多体平滑过渡我们已经从数学的角度进行了证明和测试。最终整体的组装势能面表达式如下：。（）为总能量（高低精度能量差，见公式（）；表示每个共振结构的总能量（同样是高低记得能量差），由公式（）给出；为公式（）中计算出的权重；表示不同的共振结构模式。复 旦 学 报（自然科学版）第 卷小分子势能函数数据库构建在我们的方案中，子体系以及它们之间的二体三体作用依然需要用高精度的 或 方法进行计算。为了加速这些子体系的计算，我们建立了小分子势能函数数据库 ，该势能函数数据库包含了一些常见小分子的 （）和半经验方法 之间的能量差的解析函数（如 （）、（）等）。该小分子势函数数据库主要包含个非 原子以下（仅限）的小分子，部分超过个非 原子的部分常见分子，以及 原子以下的最简拓扑结构。在实际应用中，对于某一大分子，我们首先按第部分的方法将该分子拆分成小分子，然后在数据库中搜索相应小分子的势能函数，并把它们按公式（）（）计算出该大分子的总能量。如果该小分子的势函数不在数据库中，则其能量由调用量化计算程序计算获得。这就是我们的组装势能面策略。这一策略使得绝大多数的体系的计算中需要进行高精度 的计算量降低到原有的不到。分子生成分子生成技术是新分子设计中不可缺少的部分。本工作中通常采用两种方法。第种是在常见分子骨架的活性位点上进行官能团替换，通过程序自动组装生成新的分子结构。对新分子结构组装整体势能面，在该势能面上采用 算法 研究该分子的最优构型，异构化和解离通道，获得最低反应通道以确认该分子的稳定性和最可能的异构化通道。第种是有明确的目标局部结构，我们对数据库中全部分子搜索该目标结构。进一步可以搜索该目标结构的相似结构（对目标结构进行等电等键替换形成的结构，例如目标结构的可以被等电替换为）。搜索出键连相似的结构后再替换回目标结构，排除重复后再用组装势能面和 进一步研究最优构型并进行稳定性分析。组装势能面的精度和速度我们的构造势能面仍然需要对大分子整体采用一次半经验方法（如 或者 ）进行计算，这部分计算目前占主要计算量，计算标度是体系大小的（）。子体系的部分由于主要采用数据库预存并匹配的 势函数 或者 势函数，因此计算开销相对较小（注：少部分不常见共振结构仍需要实时计算 能量和力）。整体性能要比直接采用 方法平均提高两个数量级。对于平衡结构（）的计算，本方法的 （）可以达到 的水平，因此采用构造势能面计算局域最低点（搜索构象）基本没有问题。但对于过渡态，测试集的 在 左右。这对于初步筛选是足够精确的，在得到初筛结构和能量后，还需要对 范围内的过渡态进行进一步精确的 方法或者 （）方法确认最终各条反应通道的过渡态能量精确值。应用案例本方法针对大分子体系实现了高精度、高速度势能面的构建，这对于相当多的应用都有重要意义。下面我们用两个应用实例来阐述本方法在各种不同方向上的应用。第一个是新型含能材料的设计。在这个例子中我们采用组装势能面对上千种待选分子的稳定性进行评估，筛选出稳定且高能的含能材料分子。稳定性的评估需要对每个分子进行全部解离和异构化反应通道的探索。对于这样的任务，无论采用 还是其他反应力场都是非常困难的。前者速度过于缓慢，后者由于误差高达 以上，导致稳定性的评估极不可靠。第二个例子是针对 活化的均相催化剂设计。在该案例中，我们利用现有数据库中的官能团构建目标分子结构。由于它们存在=这样的特殊结构，未进行特别修正的反应力场同样无法很好地应对。我们构建的势能面在核心结构基本不变的条件下可以实现外围基团的可靠替换，构造出的势能面具有很好的一致性。第期宋国梁等：“分而治之”方法构建复杂分子体系高精度势能面及应用新型高能化合物（）的设计含能材料的发展多年来遭遇瓶颈，主要原因在于高能和稳定性是一个巨大矛盾，能量越高的化合物稳定性通常越差 。我们对富氮五元环进行了研究，以双富氮五元环为骨架，遍历生成了其全部，及衍生物组合。然后采用我们的组装势能面技术对每个化合物搜索最优结构，计算其爆轰参数，并找出其全部解离和异构化通道。在确认化合物基本稳定且能量密度基本符合要求后，再用 方法精确计算其爆轰参数，最低解离及异构化通道的过渡态能量，以及最小的键解离能。在近千个化合物中我们发现了约 种爆速超过 且键解离能 和 超过 的 候选化合物。图显示了其中一种，该化合物预测的爆速超过 。从最稳结构出发，在步以内总共扫描到 条解离及异构化通道，其中能量最低的通道为中心 键断裂，过渡态势垒 为 。最容易解离的 的 为 。该化合物在 年已被合成出离子盐。图高通量搜索出的一种 化合物及其预测的爆轰参数 体系氮气（）活化催化剂设计活化一直是合成氨过程中的核心难题，在均相催化体系中鲜有成功案例。我们发现 体系具有类似过渡金属的活化 、的可能性。其活化的一种关键特征结构为（）（）。因此我们以该特征结构及其相似的结构（例如：，=）为骨架遍历了数据库中近 万个基本分子结构，快速筛选出近 种候选化合物，并对每种候选化合物进行解离和异构化通道搜索以确认稳定性。最终我们发现几种能够稳定吸附并活化 的 化合物。图中列出了与已知已合成化合物（）较为接近的化合物。该化合物 上的各种基团形成的位阻效应和堆叠效应恰好合适，使得 可以稳定地结合在两个 原子中间。该化合物也具有足够的稳定性，从而可以被实验合成。周围的取代基很好地保护了 路易斯酸中心，使得其他具有 的基团（如自身以及溶剂）无法接近中心而被活化，从而造成化合物的分解；同时该化合物也因位阻无法形成自身的二聚体结构，从而优先与 结合。这些有效吸附并活化 的化合物的各种解离和异构化通道最后都采用 方法进行了精确自由能计算。图高通量搜索出的 化合物用于吸附并活化 总结针对复杂大分子体系，我们发展了一种组装势能面方法。该方法将能量表达为低精度方法势能量加上高低精度方法能量差的方式。然后针对能量差，采用“分而治之”策略，将大分子体系的能量差势能面复 旦 学 报（自然科学版）第 卷分解成小的子体系的组合。采用断键多体校正方法解决断键复原；引入平滑过渡的共振结构权函数解决势能面的平滑过渡；进一步使用数据库预存小分子势能函数加速势能面计算速度；最终我们的组装势能面方法已经能将复杂大分子势能面计算量降低到（）次方的标度，同时保持计算精度接近 的水平。采用该方法，我们对两个不同的应用场景进行了测试。结果表明，高性能的组装势能面可以加速完成以往难以实现的大规模的分子结构搜索、稳定性（解离和异构化反应通道搜索）分析以及新型分子结构设计验证等关键性任务。参考文献：，（）：，（）（）？，（）：，（）：，（）：，（）：，（）：，（）：，：，：，：，：，：，：，（）：，（）：，（）：，（）：，（）：，（）：，第期宋国梁等：“分而治之”方法构建复杂分子体系高精度势能面及应用 ，（）：，（）：，（）：，（）：，：，（）：，（）：，（）：，（）：，（）：，（）：，：，：，：，（）：，（）：，：，（）：，：，（）：，：，：，（）：，（）：，：，（）：，（）：（下转第 页）复 旦 学 报（自然科学版）第 卷 ，（，）：，“”，：；（上接第 页）“”，（，）：“”（），（）（），：；复 旦 学 报（自然科学版）第 卷