一种适用于机器学习势能面构造的物质结构描述方法与流程

本发明属于计算化学和物理技术领域，具体涉及一种适用于机器学习高精度势能面构造的物质结构描述方法，用于数值上区分物质体系原子所处周围环境。

背景技术：

物质结构预测和化学反应路径搜索是当代物理和化学计算模拟研究的核心任务，对于理解和预测材料的热力学和动力学性质具有重要作用。由于物质体系的复杂性，目前可行的计算模拟方案一般依赖于基于量子力学的电子结构计算，建立体系几何结构和能量的关系，即势能面。然而，基于精确电子结构计算的势函数代价高昂，很难进行高效势能面搜索，从而不能有效筛选出关键结构和预测未知反应。因此，如何构建高精度势能面是当前科学研究的热点。

过去20年，通过大数据拟合精确电子结构计算结果，机器学习技术已经被证明可以用于构建精确的势能面。与电子结构计算方法一样，机器学习方法，也需要利用体系结构的信息，如笛卡尔坐标、成键矩阵(j.phys.chem.lett.,2015,6,2326；j.phys.chem.lett.,2017,8,2689)等，来预测体系的相关能量和其他性质。笛卡尔坐标等传统结构信息，虽然可以准确的定义原子位置，但是也具有两个明显局限。第一，坐标结构信息与原子个数一一对应，一般只能用于预测具有相同原子数体系的势能面；第二，坐标结构信息与坐标体系、原子排列的顺序等相关，不能准确反映结构内在对称性信息。这导致传统的结构输入信息不能用于产生通用性强、能处理众多原子体系的机器学习势能面。目前，已有尝试采用不同的结构输入信息，如结构指纹函数(j.chem.inf.model.,2010,50,742)、图像卷积(arxiv:1509.09292；acscent.sci.,2016,2,725)等，这些函数虽然是旋转不变，但是它们一般是对整体结构进行操作，局部结构区分度低，也不能解决从小体系到大体系原子数变化带来的势能面通用性问题。

近期，belher等课题组(phys.rev.lett.,2007,98,146401；j.chemphys.,2011,134,074106)提出了一类新的特征函数，与之前的区别是，这些特征函数不是描述整体结构特征，而是单个原子的周围环境。通过采用截断函数，中心原子的周围结构信息可以由在截断半径内的其他邻居原子与中心原子的键长、键角来确定。具体是通过建立基于这些键长、键角的特征函数，作为机器学习的输入参数，拟合得到单个原子的能量。这种构建势能面的方法，将体系能量有效拆分为单个原子能量的简单加和，使得机器学习势能面得以通用于不同原子数体系。但是，目前这些以原子为中心的结构特征函数，形式相对比较简单(一般是高斯函数，三角函数)，涉及参数多且复杂，缺乏三体以上的结构描述，不能有效形成统一的结构描述标准，不利于通用机器学习势能面的大规模构建和推广。

技术实现要素：

本发明的目的是为提供一种适用于机器学习高精度势能面构造的物质结构描述方法。

本发明首先构造一种简便、普适、标准化的结构特征函数。这些结构特征函数用于构造机器学习通用型势能面，可以广泛用于全局结构搜索，自由能势能面采样，反应预测等分子模拟。该方法克服了以往经验势函数精度不足，及量子力学势能面效率不高的缺点，具有很强的可移植性和通用性，能对不同材料体系均取得良好预测能力。

本发明提供系列标准化的结构特征复杂函数，它们基于邻居原子与中心原子的键长、键角，由截断函数、径向函数、角向函数和球谐函数四类基本函数组合而成。其中，截断函数主要由双曲正切函数构成，径向函数由幂指数和截断函数构成，角向函数由三角函数构成，球谐函数是拉普拉斯方程的球坐标系形式解的角度部分。这些函数的有效组合，可以高分辨区分原子周围结构环境，具体形式如下：

a)截断函数：

i，j＝1,2,…n，为原子序号且i≠j(以下内容相同定义)，tanh为双曲正切函数，rij为i，j原子间的距离，rc为截断半径，选取范围为1.0到10埃。

b)径向函数：

这里，为i，j原子间距离的幂函数，n为可调参数，通常取值范围-3～24。附图1显示了，在截断半径为6埃下，径向函数随次方数变化的曲线。显然，通过调节不同的次方数，可以对距离中心原子不同远近的邻居的信息进行径向取样，最终达到对不同径向位置侧重取样的目的。

c)角向函数：

a¹(θijk)＝[(1+λcosθijk)/2]^ζ(sin²θijk)^ξ

a²(δijkl)＝[(1+λcosδijkl)/2]^ζ(sin²δijkl)^ξ

这里，cos为余弦三角函数，sin为正弦三角函数，i，j，k，l＝1,2,…n，为原子序号且i,j,k,l两两不等(以下内容相同定义)；θijk为以i原子为中心，i，j，k原子之间的夹角；δijkl为以i,j原子为中心，i,j,k,l原子之间的二面角；λ，ζ，ξ为可调参数，且λ为+1或者-1，ζ，ξ取范围0～24的整数。

举例来说，附图2显示了a¹(θijk)函数在不同λ，ζ，ξ取值下，角度从0度变化到180度函数值的变化。图2(a)说明，当固定ξ＝0，ζ为不同整数时，λ为+1和-1分别着重对锐角和钝角取样。ζ指数取值为4，配合λ为+1和-1时，可以对整个角度范围都有适当的取样，且随着ζ指数的增加，对应取值范围逐渐变小，偏离90度。图2(b)说明，通过固定ζ指数，而调节ξ，可以达到对不同角度侧重取样的目的。

d)球谐函数

球谐函数是拉普拉斯方程的球坐标系形式解的角度部分，球谐函数可以从直角坐标系转化为球坐标系：

其中，表示向量的x分量，下标y，z的含义相同。在球坐标系中，球谐函数可以表达为：

这里的ylm称为l和m的球谐函数，i是虚数单位，是伴随勒让德多项式，用方程定义为：

pl(x)是l阶勒让德多项式，可用罗德里格公式表示为：

本发明中，根据以上径向函数、球谐函数，对于体系内任意原子(以i标记)，构造具有如下形式的两种双体特征函数

为球谐函数；nbonds为加和项所有原子对的个数；为j原子指向i原子构成的空间向量；l为可调参数，通常取值范围是2、4、6。

举例来说，图3显示了常见原子的不同配位环境，包括平面正多边形(线性、平面三角形、正方形、五边形及六边形)和三维正多面体(正四、六、八、十二及二十面体)配位，其中粉色原子为中心原子，蓝色原子为配位原子，且所有原子间的键长均为1.0埃。图4显示了这些不同配位环境下，其中心原子的双体特征函数的数值，其中，可调参数l＝2，4，6，双体函数的径向部分截断函数的半径为1.8埃，幂函数的次方数为2。如图所示，结合不同的特征函数值可以达到很好区分中心原子的配位环境。

本发明中，根据以上径向函数、角向函数、球谐函数，对于体系内任意原子(以i标记)，构造具有如下形式的三种三体特征函数

。

本发明中，根据以上径向、角向函数，对于体系内任意原子(以i标记)，构造具有如下形式的四体特征函数

本发明中，所有的特征函数：都可以根据需要规定函数中原子所属于的元素种类。对于指定元素，在计算特征函数时，若涉及的原子不属于指定元素，则其对特征函数贡献为零。

本发明中，在给定数据集和给定网络参数个数条件下，所使用的结构特征函数的个数以及组合方式决定了机器学习势能面精度。对于任意原子(以i标记)，使用下列结构特征函数的组合以及个数作为一种标准组合，作为训练高维机器学习势能面的输入层，可以使机器学习势能面精度达到能量误差小于15毫电子福特每原子(mev/atom)的常见精度(误差以在数据集中，机器学习势能面和量子力学势能面的单原子能量标准均方差表示，以下能量精度、误差相同含义)。该标准组合为：

单元素体系

(一)24个双体函数,包含16个和8个

(二)16个三体函数，包含4个8个和4个

(三)2个四体函数，包含2个

实施例1中，为pt体系，该标准组合的格式和参数取值见附表1-6。

以此类推，多元素体系需要不同元素组合以形成特征函数，标准组合的个数如下：

双元素体系

(一)48个双体函数,包含32个和16个

(二)48个三体函数，包含12个24个和12个

(三)8个四体函数，包含8个

实施例2中，为金属氧化物氧化锰mnox体系，该标准组合的格式和参数取值见附表7-12。

n个元素体系。

(一)24×n个双体函数,包含16×n个和8×n个

(二)16×mn个三体函数，包含4×mn个8×mn个和4×mn个

(三)2×pn个四体函数，包含2×pn个

其中

综上，本发明提供的适用于机器学习势能面构造的物质结构描述方法，是通过构建一系列数学函数来描述物质结构特征，用于构建机器学习复杂体系势能面；具体来说，利用物质体系的原子笛卡尔坐标qi(i＝1,2,…n为原子序号,n为体系的原子数)，根据给定的截断半径构造截断函数和径向函数,再结合角向函数和球谐函数，形成六类双体、三体、四体的结构特征函数：这些结构特征函数具有旋转不变、原子交换不变、一阶和二阶导数连续性质；这些结构特征函数的组合作为训练高维机器学习势能面的输入层，用于构建高维机器学习势能面。

本发明中，高维机器学习势能面中机器学习包含如下形式：

神经网络(包括：前馈网络，卷积网络)；

支持向量机；

决策树和随机森林。

机器学习所需要的势能面大数据来源于基于量子力学电子结构计算的全局势能面扫描，主要用到中国发明专利zl201210421939.6的随机表面行走(ssw)算法对复杂势能面进行高效采样。每一个结构数据均包含体系总能量、原子的力、体系的应力。

机器学习的训练过程，需要拟合每个结构数据的总能量、原子的力、体系的应力。通过最小化目标损失函数，训练得到机器学习势能面。

按本发明，得到的材料势能面可以直接用于具有不同元素、或众多原子数的复杂体系(比如原子数大于1000的体系)计算模拟。具体模拟内容包括，全局势能面结构搜索、分子动力学模拟、反应通道搜索和反应机理研究等。

相比于传统的势能面构建方法，本发明提出的依靠新型结构特征函数的机器学习势能面构造方法，具有如下优点：1.适用的材料体系范围广，精度高，可以用于不同原子数、不同组分的材料全局势能面的预测；2.自动化程度高，对科研人员的专业背景要求低，一般只需要根据标准组合模式，就可以用于构建高精度的全局势能面，得到的势能面可以直接用于不同要求的分子模拟。

附图说明

图1：在截断半径为6埃下，径向函数随次方数变化的曲线，为了比较不同次方数的形状，所有的曲线均相对于4.5ⁿ，n与内的可调参数一致。

图2：不同可调参数取值下，角向函数a¹(θijk)的取值范围变化。

图3：对于常见原子的配位环境，对应双体特征函数的变化。

图4:对于常见原子的配位环境，心原子的双体特征函数的数值，其中，可调参数l＝2，4，6，双体函数的径向部分截断函数的半径为1.8埃，幂函数的次方数为2。

具体实施方式

下面通过实施例来详述本发明，但本发明的内容并不局限于此。

实施例1：

金属铂pt全局势能面的构建。运用第一性原理密度泛函势能面及势能面随机搜索算法进行全局数据集采样，共得到26063个结构，包括体相、层状、团簇结构。数据集的结构能量范围为6电子福特每原子(能量参考点为全局势能面最稳定结构)，力的范围为10电子福特每埃，张力的范围为260兆帕斯卡。使用本发明中的标准特征函数集作为输入信息(共42个特征函数，其中双体特征函数24个，三体特征函数16个，四体体征函数2个)，采用前馈神经网络训练该数据集，得到的全局势能面精度为能量误差9.9毫电子福特每原子，力误差0.11电子福特每埃，张力误差1.3兆帕斯卡。具体标准特征函数集参数如表1-6所示。

实施例2：

金属氧化物氧化锰mnox(含不同价态mn)全局势能面的构建。运用第一性原理密度泛函势能面及势能面随机搜索算法进行全局数据集采样，共得到102134个结构，包括体相、层状、团簇结构。数据集的结构能量范围为3.2电子福特每原子(能量参考点为全局势能面最稳定结构)，力的范围为40电子福特每埃，张力的范围为84兆帕斯卡。使用本发明中的标准特征函数集作为输入信息(每个元素含共104个特征函数，其中双体特征函数48个，三体特征函数48个，四体体征函数8个)，采用前馈神经网络训练该数据集，得到的全局势能面精度为能量误差8.5毫电子福特每原子，力误差0.15电子福特每埃，张力误差1.3兆帕斯卡。具体标准特征函数集参数如表7-12所示。

表1：pt的特征函数个数及其可调参数

表2：pt的特征函数个数及其可调参数

表3：pt的特征函数个数及其可调参数，其中，n，m分别为径向部分对应rij和rik的幂指数参数

表4：pt的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和rjk的幂指数参数

表5：pt的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和rjk的幂指数参数

表6：pt的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和ril的幂指数参数

表7：mno2的特征函数个数及其可调参数

表8：mno2的特征函数个数及其可调参数

表9：mno2的特征函数个数及其可调参数

表10：mno2的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和rjk的幂指数参数

表11：mno2的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和rjk的幂指数参数

表12：mno2的特征函数个数及其可调参数，其中，n，m，p分别为径向部分对应rij，rik和ril的幂指数参数