CN108733973B - 一种自动高效dftb排斥势拟合方法 - Google Patents

Abstract

本发明提供一种自动高效DFTB排斥势拟合方法，包括以下步骤:步骤(1):优化包含原子对的分子；根据排斥势作用区间,进行原子对间距扫描；步骤(2):对扫描结果进行高精度能量计算；步骤(3):能量差分计算力和能量的二阶导数；步骤(4):将结构、能量、力和能量二阶导数存入数据库；步骤(5):得到未包含排斥势的能量、力和能量二阶导数,获得排斥势目标值；步骤(6):将分离的排斥势目标值转换成样条曲线；步骤(7):将多条样条曲线在平衡位置点间的平均位置拼接,得到全新的训练数据;步骤（8）：用必须样本和选择样本构建训练集；步骤（9）：使用奇异值分解拟合排斥势。本发明克服了键类型不平衡拟合问题。

Description

一种自动高效DFTB排斥势拟合方法

技术领域

本发明属于生物信息技术,涉及功能材料的设计,尤其涉及一种自动高效DFTB排斥势拟合方法。

背景技术

DFTB是一种快速计算的高精度量子化学方法,做到了计算规模和精度的有效平衡,在生物蛋白模拟和功能材料设计等领域有广泛的应用。排斥势是DFTB参数的重要组成部分,其可靠性严重影响计算模拟的准确性。

目前的DFTB参数化方案通过手动建立训练集,采用大量训练的样本点进行样条曲线拟合,来确定排斥势曲线，由于同一个排斥势需要描述两个原子对的多种成键形式,如果使用大量训练的样本点进行拟合,极易造成无法准确覆盖不同成键形式,如,C...O相互作用有C-O/C＝O两种成键方式,平衡位置明显不同,直接样条拟合难以平衡两种键长。

样条拟合的方式一次只能拟合一条原子对曲线,拟合顺序具有明显的先后关系,不仅对训练集的复杂度要求极高,而且会造成拟合排斥势的兼容匹配问题。比如要开发O...O相互作用,必须先准备H...H和O...H排斥势,才可以进行O...O的开发,但是开发的O...O相互作用不确定是否适合于除了这三种作用的体系中.会造成开发参数流程的反复迭代,开发速度较慢，并且拟合涉及大量的训练集和测试集的建立,参数调整,测试校正,需要进行大量的人工手动干预。

发明内容

为了克服现有技术的不足,本发明提供了一种基于差分进化算法的自动高效DFTB排斥势能拟合方法,包括以下几个步骤:

步骤(1):根据目标原子对,根据原子自身的信息生成氢饱和的分子体系,使用Gaussian09等量子化学软件中的精度较高的密度泛函理论(DFT)方法，如B3LYP/cc-pVTZ和PBE0/cc-pVTZ方法优化分子；根据排斥势的作用区间,使用相同精度的DFT方法对原子对进行间距扫描；

步骤(2):使用Gaussian09等软件中高精度方法，如：G3B3、CCSD/CBS方法对间距扫描进行高精度能量计算；

步骤(3):使用高精度方法获得的能量差分计算力和能量的二阶导数；

步骤(4):将结构、高精度方法计算得到的能量,力和能量的二阶导数保存入数据库；

步骤(5):将结构使用DFTB+等可以实现DFTB计算的程序得到未包含排斥势的能量,力和能量的二阶导数,与高精度方法计算结果做差,得到排斥势拟合的目标值；

步骤(6):根据原子间距和排斥势的目标值，进行样条插值,获得多条样条曲线,将分离的能量样本转换成为连续的能量曲线；

步骤(7):将平衡位置点间的平均位置为拼接点,并根据拟合终点的位置,将终点能量确定为0.0,从拟合终点开始倒序,在分立点处将多条样条曲线进行拼接,得到全新的训练能量数据；

步骤(8):选择离目标排斥势起点和终点最近的两个点,以及平衡位置点作为必须样本点,选择离拟合起点较近2-3个样本点为补充样本点,构建训练集池；

步骤(9):使用奇异值分解方案,拟合得到初始排斥势；再进一步使用初始排斥势的拟合方法,将初始排斥势转化为排斥势的标准形式，使用标准的排斥势进行测试验证。

本发明采用以上技术方案，其优点在于，使用平衡键长点和少量样本点,并基于差分进化算法的确定分立点,联合拟合多条排斥势曲线。

由于为保证拟合效果,除了正常的划分键类型区域的分立点,还需要插入1-2个训练样本,以及相应的分立点.而且分立点的选择在一定程度上影响拟合效果.因此,需要优选分立点.这里使用差分进化算法进行优选训练集合和分立点，包括以下步骤：

步骤A：使用起点-平衡点,平衡点-平衡点,平衡点-终点的平均位置初始化分立点,初始化分离点.这里,平衡点-平衡点间的至少包含一个分立点；如图5所示。

步骤B：根据样本和分立点密度,从训练集池中随机抽取1-2个补充样本,保持分离点区间至少包含一个样本,并对分立点进行不越过必须样本的平移扰动,以此产生15-50个个体；

步骤C：每个个体进行二次拟合，即：原始排斥势,目标排斥势的拟合,并验证评估误差：第一次,使用随机或演化产生的分离点线性拟合得到初始排斥势；均匀划分分立点的划分，取得初始排斥势上的能量,力和能量的二阶导数，并线性拟合得到目标排斥势；使用拟合得到的排斥势进行计算,将相对能量和结构相似性作为打分函数,对每次的拟合进行评估；

步骤D：使用差分进化算法产生分离点扰动和样本扰动；对每一代拟合按照打分评估的值进行排序,并将分立点和补充样本点作为一个矢量.将最优个体矢量和其它个体矢量做差,生成微扰方向矢量；对其它矢量进行多个微扰方向的扰动,产生新的分立点和补充样本；如果新的分立点越界,则随机产生一个满足条件的分立点；如果新的补充样本不再训练集内,则在训练集中的补充样本中随机抽选一个作为新的补充样本；

步骤E：迭代步骤C-D,直到误差收敛。

本发明带来了如下效果：

1.由于拟合过程中,重点使用了不同键类型的平衡位置,能量和力,仅使用少数样本点辅助拟合,平衡位置的信息重点保留,且权重近似相等,克服了键类型不平衡拟合问题。

2由于使用的线性方程组的方式拟合,可以将多条排斥势的拟合平时放进一个方程组里,实现多条排斥势能的拟合,解决了参数的兼容问题。

3使用差分进化算法和数据库,实现排斥势的自动拟合。

附图说明

图1为本发明自动高效DFTB排斥势能拟合方法的流程图。

图2是本发明中的一种实施例，两种C...O键类型扫描得到的C-O原子间距VS排斥势能量样本,其中黑点为H₂C＝O(灰色为其平衡位置),圈为H₃C-HO(虚线圈为其平衡位置)。

图3是本发明中的一种实施例，拟合得到两种C...O键类型对应的样条曲线。

图4本发明中的一种实施例，将两条C...O键样条曲线拼接成一条曲线,平衡两种键类型，选择离目标排斥势起点和终点最近的两个点,以及平衡位置点作为必须样本点,选择离拟合起点较近2-3个样本点为补充样本点,构建训练集池。

图5是本发明中的一种实施例，选择起点,终点和平衡点为必须样本点,选择平衡点附近的2-3个样本为补充样本。

图6本发明中的一种实施例，分立点的初始化。

图7图5中根据必须样本点个数将拟合区间划分为4个区间,随机向第三个区间内添加一个平衡位置附近的样本,并对分立点进行平移。

图8是本发明中的一种实施例，使用线性方程组拟合得到初始排斥势。

图9是本发明中的一种实施例，均匀划分分立点的划分。

具体实施方式

下面结合附图，对本发明的较优的实施例作进一步的详细说明：

实施例1

如图1所示，一种基于差分进化算法的排斥势能拟合方法,包括以下步骤:

步骤(1):根据目标原子对,根据原子自身的信息生成氢饱和的分子体系,使用Gaussian09中的B3LYP/cc-pVTZ方法优化分子对；根据排斥势的作用区间,使用B3LYP/cc-pVTZ方法对原子对进行间距扫描；

步骤(2):使用Gaussian09中G3B3方法对间距扫描进行高精度能量计算；

步骤(3):使用G3B3方法获得的能量差分计算力和能量的二阶导数；

步骤(4):将结构、G3B3方法计算得到的能量,力和能量的二阶导数保存入数据库；

步骤(5):将结构使用DFTB+等可执行DFTB计算的程序得到未包含排斥势的能量,力和能量的二阶导数,与G3B3结果做差,得到排斥势的相应值。

步骤(6)：如图2所示，两种C...O键类型扫描得到的C-O原子间距VS排斥势能量样本,其中黑点为H₂C＝O(灰色为其平衡位置),圈为H₃C-HO(虚线圈为其平衡位置)。根据多个键类型的扫描结果,进行样条插值,获得多条样条曲线,将分离的能量样本转换成为连续的能量曲线,方便后续能量曲线的合并,以及多种键类型的拟合；

步骤(7)：如图3所示，拟合得到两种C...O键类型对应的样条曲线；将平衡位置点间的平均位置为拼接点,并根据拟合终点的位置,将终点能量确定为0.0,从拟合终点开始倒序,在分立点处将多条样条曲线进行拼接,得到全新的训练能量数据,此处初步处理了多种键类型问题；如图4所示，将两条C...O键样条曲线拼接成一条曲线,平衡两种键类型。

步骤(9):选择离目标排斥势起点和终点最近的两个点,以及平衡位置点作为必须样本点,选择离拟合起点较近2-3个样本点为补充样本点,构建训练集池；

如图5所示：选择起点,终点和平衡点为必须样本点,选择平衡点附近的2-3个样本为补充样本

定义初始排斥势曲线

的形式为多段二阶连续可导4次导多项式，其中第i段的表达式为：

其中，A和B为原子对中原子对应的元素，k为次数，i为原子间距r对应分立区间序号，s_ABik为第i个区间第k次多项式的系数，r_ABi为第i个区间的起点位置。因此，整个体系的排斥能(E^rep)可以表达为：

即，根据体系中的原子a和b的元素E_a和E_b选择元素对应的排斥势多项式，并根据原子间距r_ab,确定所属的分立区间i，然后计算得到排斥能；

类似的，原子a的排斥力表达式为：

其中

根据初始排斥势的表达,可以将拟合表达为线性方程组的形式:M·s＝b,这里的M根据多段多项式的形式、分立点上的二阶连续可导信息和样本点信息确定的矩阵，s为全部待拟合的参数确定的矢量，b为二阶连续可导和样本点信息对应的能量，导数，二阶导数的数值。

根据排斥势的连续可导(cont),和样本点拟合数值(fit),将方程整理成如下形式:

这里的M^cont是仅包含分立点上的二阶连续可导信息的矩阵，而M^fit为仅包含样本点信息的矩阵。b^fit为仅包含样本点能量、导数，二阶导数的数值的矢量。由于矩阵变换和连续可导，b^cont为0矢量。进一步，将s矢量按照精确拟合分解为由样本点所在的区间确定的必须精确拟合部分(expl)，以及其余部分(other)：

M^cont,others^cont,other+M^cont,expls^cont,expl＝0

M^fit,others^fit,other+M^fit,expls^fit,expl＝b^fit

将上两式变换并合并得到下式：

M^last·s^expl＝b^fit

使用奇异值分解方案,拟合得到初始参数；由于初始排斥势形式与最终的形式存在一定差异,我们进一步使用初始排斥势的拟合方法,将初始排斥势转化为排斥势的标准形式(二次拟合)，使用标准的排斥势进行测试验证.

由于为保证拟合效果,除了正常的划分键类型区域的分立点,还需要插入1-2个训练样本,以及相应的分立点.而且分力点的选择在一定程度上影响拟合效果.因此,需要优选分力点.这里我们使用差分进化算法进行优选训练集合和分立点，包括以下步骤：

步骤A：使用起点-平衡点,平衡点-平衡点,平衡点-终点的平均位置初始化分立点,初始化分离点.这里,平衡点-平衡点间的至少包含一个分立点.如图6所示。

步骤B：根据样本和分立点密度,从训练集池中随机抽取1-2个补充样本,保持分离点区间至少包含一个样本,并对分立点进行不越过必须样本的平移扰动,以此产生20个个体；

步骤C：每个个体进行二次拟合，即：原始排斥势,目标排斥势的拟合,并验证评估误差；第一次,使用随机或演化产生的分离点线性拟合得到初始排斥势.如图7所示为图6中根据必须样本点个数将拟合区间划分为4个区间,随机向第三个区间内添加一个平衡位置附近的样本,并对分立点进行平移；进一步将初始排斥势能转化为样条曲线,进一步根据起点和终点的距离,以根据实际情形，选择0.02-0.05Bohr为间隔将拟合区间进行均匀划分,将样条数据转换为线性方程组,求解线性方程,得到目标排斥势曲线.如图8所示。

如图9所示，均匀划分分立点的划分：使用拟合得到的排斥势进行计算,将相对能量和结构相似性作为打分函数,对每次的拟合进行评估,能量误差越小,结果误差越接近，拟合效果越好。

步骤D：使用差分进化算法产生分离点扰动和样本扰动；对每一代20次(20个个体)拟合按照打分评估的值进行排序,并将分立点和补充样本点作为一个矢量.将最优个体矢量和其它个体矢量做差,生成微扰方向矢量；对其它矢量进行多个微扰方向的扰动,产生新的分立点和补充样本；如果新的分立点越界,则随机产生一个满足条件的分立点；如果新的补充样本不再训练集内,则在训练集中的补充样本中随机抽选一个作为新的补充样本。

步骤五：迭代步骤C-D,直到误差收敛。

实施例2

由于拟合过程中,重点使用了不同键类型的平衡位置、能量和力,仅使用少数样本点辅助拟合,平衡位置的信息重点保留,且权重近似相等,克服了键类型不平衡拟合问题；如：C—O需要拟合2.0到4.6Bohr区间内排斥能,按照0.05的均匀分立区间,至少需要52个训练样本.如前图2所示,平衡位置处的排斥势相关参数较小,起点处参数较大,造成拟合权重从起点到终点递减拟合结果倾向获得更近的平衡位置；当不同键类型时,会造成排斥势相关参数更大的键类型相应更可靠.C-O中,C＝O键会描述相对更准.如直接使用样条拟合C＝O键和C-O键的距离为1.18和1.25A,而实验C＝O键和C-O键的距离为1.23A和1.43A。

由于化学、材料和生物等研究中,平衡位置附近才是研究者关心的重点.因此,在拟合C-O过程中,实际上仅使用2个平衡点位置的参数,为了得到更准确的拟合效果,并使用1-2个平衡点附近的样本作为补充,平衡了多个键类型的拟合权重.如单次拟合后,C＝O键和C-O键的距离可达1.22和1.42A.明显优于样条直接拟合。

实施例3

由于使用的线性方程组的方式拟合,可以将多条排斥势的拟合平时放进一个方程组里,实现多条排斥势能的拟合,解决了参数的兼容问题；举例,对于C-O的样条拟合，首先要准备好H-H,C-C,O-O,O-H,C-H的全部参数；如果缺少任意参数,则必须准备好相应的参数,才可继续.而且,参数间存在耦合兼容问题.如,H-H间的作用会造成C-O训练样本值的计算,影响C-O的拟合.要得到兼容的参数,必须对多套参数进行反复的调整,确定最终的排斥势参数,效率低下；该方法使用线性方程组进行拟合,可以将H-H的参数和C-O参数同时放进一个线性方程组中进行求解,一次得到对两种原子对兼容性好的参数.

本发明使用差分进化算法和数据库,实现排斥势的自动拟合。样条拟合可以通过反复的测试,以及改变相关点的拟合权重调整,将平衡键长调整至合理,手动操作频繁,且效果通常不理想；本发明所采用的拟合在解决键类型平衡和参数兼容性的问题基础上,通过差分进化算法调整分立点,实现排斥势的微调,进一步将拟合的效果提升。综合以上结果,本发明实现了排斥势的自动拟合。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (8)

1.一种自动高效DFTB排斥势拟合方法，其特征在于，包括以下几个步骤:

步骤(1):优化包含原子对的分子；根据排斥势的作用区间, 对原子对进行间距扫描；

步骤(2):对间距扫描获得的结构进行高精度能量计算；

步骤(3):对获得的能量进行差分计算得到力和能量的二阶导数；

步骤(4):将结构、高精度方法计算得到的能量, 力和能量的二阶导数保存入数据库；

步骤(5):将结构计算得到未包含排斥势的能量, 力和能量的二阶导数, 与高精度方法得到的结果做差, 得到排斥势的目标值；

步骤(6):根据原子间距和排斥势能量的目标值，进行样条插值, 获得多条样条曲线,将分离的样本转换成为连续的曲线；

步骤(7):将平衡位置点间的平均位置为拼接点, 并根据拟合终点的位置, 将终点能量确定为0,从拟合终点开始倒序, 在分立点处将多条样条曲线进行拼接, 得到全新的训练能量数据；

步骤(8):选择离目标排斥势起点和终点最近的两个点, 以及平衡位置点作为必须样本点, 选择离拟合起点较近2-3个样本点为补充样本点, 构建训练集池；

步骤(9):使用奇异值分解方案, 拟合得到初始参数；再进一步使用初始排斥势的拟合方法, 将初始排斥势转化为排斥势的标准形式,即，二次拟合，使用标准的排斥势进行测试验证。

2.如权利要求1所述的方法，其特征在于，所述步骤(1)中，根据目标原子对, 根据原子自身的信息生成氢饱和的分子体系,采用Gaussian09量子化学计算软件中的B3LYP/cc-pVTZ中精度方法优化分子对；根据排斥势的作用区间,采用相同方法对原子对进行间距扫描。

3.如权利要求1所述的方法，其特征在于，所述步骤(2)中，采用Gaussian09量子化学软件中G3B3高精度计算方法对间距扫描进行高精度能量计算。

4.如权利要求1所述的方法，其特征在于，所述步骤(3)中，采用高精度方法获得的能量差分计算力和能量的二阶导数。

5.如权利要求1或4所述的方法，其特征在于，所述步骤(4)中，所述高精度采用MP2、G3B3和CCSD/CBS方法中的一种。

6.如权利要求1所述的方法，其特征在于，所述步骤(5):将结构使用DFTB+进行DFTB计算的程序得到未包含排斥势的能量, 力和能量的二阶导数, 与高精度结果做差, 得到排斥势的相应值。

7.如权利要求1所述的方法，其特征在于，所述步骤（8）中需要优选分立点，采用差分进化算法进行优选训练集合和分立点，包括以下几个步骤：

步骤A：使用起点-平衡点, 平衡点-平衡点, 平衡点-终点的平均位置初始化分立点,初始化分离点；

步骤B：根据样本和分立点密度, 从训练集池中随机抽取1-2个补充样本, 保持分离点区间至少包含一个样本, 并对分立点进行不越过必须样本的平移扰动, 以此产生个体；

步骤C：每个个体进行拟合，即：初始排斥势, 目标排斥势的拟合, 并验证评估误差；

步骤D：使用差分进化算法产生分离点扰动和样本扰动；对每一代拟合按照打分评估的值进行排序, 并将分立点和补充样本点作为一个矢量；将最优个体矢量和其它个体矢量做差, 生成微扰方向矢量；对其它矢量进行多个微扰方向的扰动, 产生新的分立点和补充样本；如果新的分立点越界, 则随机产生一个满足条件的分立点；如果新的补充样本不再训练集内, 则在训练集中的补充样本中随机抽选一个作为新的补充样本；

步骤E：迭代步骤C-D, 直到误差收敛。

8.如权利要求7所述的方法，其特征在于，所述步骤A中，平衡点-平衡点间的至少包含一个分立点。

Images