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