CN106599408A - 一种基于改进BHMC算法的FenPtm合金团簇结构优化方法 - Google Patents

Abstract

一种基于改进BHMC算法的Fe_nPt_m合金团簇结构优化方法，包括以下步骤：1)初始化种群，得到原子总数为N的合金团簇在各种不同比例情况下的初始结构；2)对初始化的结构进行最速下降局部优化，得到初始后化的局部最优结构；3)对每个局部最优结构的每个原子进行扰动处理，并进行最速下降局部优化，保存更低能量结构；4)对能量最低的结构进行遗传局部优化处理，得到当前状态下最优结构的元素最优排布；5)不断重复扰动算子操作和当前最优结构下的遗传局部优化，直到搜索到最终的最优结构。本发明对计算机硬件开销要求较低，能有效处理合金团簇出现的同分异构体问题，避免陷入局部最优，从而能够快速有效地实现对Fe_nPt_m合金团簇的结构寻优。

Description

一种基于改进BHMC算法的FenPtm合金团簇结构优化方法

技术领域

本发明涉及一种基于改进BHMC(Basin-Hopping Monte Carlo)算法的Fe_nPt_m合金团簇结构优化方法。

背景技术

二元合金纳米团簇可以充分利用两种金属之间的电子结构与晶体结构的协同效应来实现材料的双功能特性，因而在化学和物理方面表现出优于单金属团簇的性质，在催化、生物医学、新材料等领域具有广泛的应用。在金属团簇中，铁(Fe)团簇作为过渡金属，不仅具有磁性，而且与其他贵金属团簇一样，具有高反应性和催化活性，因此在信息存储、能源化工和生物医学等领域中具有广阔的应用前景。铁与铂(Pt)、钯(Pd)和金(Au)等贵金属结合形成二元或多元合金，不仅可以显著提高铁团簇的催化活性，还能实现团簇的多功能特性。

团簇的结构决定了它们的物理或化学性质，因此研究合金纳米团簇的稳定结构是研究合金团簇其他性质的基础。另外，团簇的结构依赖于团簇尺寸、原子比例和原子有序化分布，因此可以通过改变团簇的这些特征来改变团簇的物理与化学性质。目前，对于合金团簇结构稳定性研究，最常用的方法是基于密度泛函的计算方法和基于原子间经验势函数的全局优化算法。基于密度泛函理论的计算方法计算精度高，但对于原子数较多且存在大量同分异构体的合金团簇的结构研究，计算量相当大，难以获得团簇结构与其尺寸、元素比例的依赖关系。而利用经验势函数来描述原子间相互作用，将合金团簇结构稳定性的研究转换成一个全局最优化问题，使得很多用于求解全局优化问题的方法如蒙特卡洛(MonteCarlo)，盆地跳变算法(Basin-Hopping)，模拟退火算法和进化算法等可以用来研究团簇的结构优化。然而，随着原子数的增长，势能面上的局部最小值也呈指数增长，并且二元合金团簇将会出现大量的同分异构体，因此，探索一种有效的全局优化算法对于合金团簇结构优化十分重要。

纳米团簇中原子间的相互作用通常采用经验势函数来描述，研究Fe-Pt合金团簇稳定结构的实质就是寻找势能面上的能量最小值。能量越低的团簇，其结构越稳定。为研究Fe-Pt合金团簇的结构稳定性和结构特征，首先需要计算基于势能函数的能量最小值，从而获得最稳定的几何结构。寻找势能最小值是一个全局优化问题，而合金团簇的结构寻优问题可以理解为结构优化和元素原子排布优化的组合优化问题，因此求解该问题的算法必须同时具备较好的结构优化特性和元素原子排布优化特征两个方面的能力。如图1所示，结构优化是指获取某个原子数下合金团簇能量值最低的结构(图1a)，元素原子排布优化是指在某个结构下，获得合金元素在该结构内不同位置分布下的最低能量值结构(图1b)。由于在所有随机优化算法中，算法的收敛速度和搜索能力在很大程度上取决于在当前解的基础上探索新解的策略。例如传统的BHMC算法，主要利用简单的原子位移来产生新的构型，随着原子数目的增加，算法的有效性将会变差。因此，有研究者在传统BHMC算法的基础上，改进了扰动策略，引入表面削弱和位置互换等，这种类型的改进，能够很好地处理二元合金团簇结构优化问题，但不能很好地解决二元合金团簇元素原子排布优化的问题。

发明内容

为此，本发明提出了一种基于改进BHMC算法的Fe_nPt_m合金团簇结构优化方法，以提高对铁(Fe)铂(Pt)二元合金元素原子排布优化的能力。本发明的具体技术方案是：

一种基于改进BHMC算法的Fe_nPt_m合金团簇结构优化方法，包括以下步骤：

S1、初始化种群，得到原子总数为N的Fe_nPt_m合金团簇在各种不同比例情况下的M个初始结构，其中，N＝n+m，M为种群中个体数；

S2、对初始化的每个结构进行最速下降局部优化，得到初始化后的局部最优结构；

S3、对每个局部最优结构进行原子移位扰动和球切扰动处理，并进行最速下降局部优化，对比处理前后两个结构的能量，保存较低能量结构，在完成种群中M个个体扰动及最速下降局部优化处理后，可以得到M个能量较低的结构；

S4、对S3得到的每个能量较低结构进行遗传局部优化处理，得到当前状态下最优结构的原子最优排布，其中，所述遗传局部优化处理包括：

S41、对所述步骤S3得到的能量较低结构，令每个原子坐标保持不变，即固定每个原子位置；

S42、利用遗传局部种群初始化方法，通过改变每个位置上的原子类型(Fe或Pt)来产生在能量较低结构下的随机原子排布初始种群；

S43、利用遗传算法，分别对原子排布进行交叉、变异和选择操作，以实现原子排布相对寻优，其中，交叉采用交叉率为0.4的两点交叉，变异采用概率为0.1的单点变异，选择采用轮盘赌方式；

S44、通过不断重复S43，找到在固定原子位置下的原子最优排布；

S5、对步骤S4得到的每个最优结构进行最速下降局部优化处理，保存更低能量结构为全局最优结构，并与上一代得到的全局最优结构进行能量对比，如果前后两次能量差在预定精度范围内，表示能量不变，在有预定代数能量不变时或达到最大进化代数时，则寻优结束，否则回到S3。

进一步的，所述最速下降局部优化处理的具体步骤是:

1)首先利用Gupta势能函数计算能量E，其中，原子总数为n的团簇的Gupta势能函数为，

其中r_ij表示原子i与原子j间的距离，r_ij ⁽⁰⁾表示两类原子间的特征长度，一般以块体材料的第一近邻距离表示，A_ij是衡量原子间的排斥强度的量，ξ_ij是有效跳跃积分，通常只与原子类别有关；

2)再对Gupta势能函数(1)求导，计算得到梯度F；

3)利用迭代公式X_i+1＝X_i-ρ_iF更新原子坐标位置，其中ρ_i是梯度方向上的步长，X_i是第i步的坐标变量，X_i+1是第i+1步的坐标变量；

4)将X_i＝X_i+1重复步骤3)不断迭代，直到满足前后两次迭代能量差ΔE<ε，其中ε为1e-6。

进一步的，所述步骤S1的所述初始结构种群是通过引入相同原子个数的Fe和Pt单金属的最稳定结构作为初始构型来产生的。

进一步的，所述步骤S3中的所述原子移位扰动是对原子坐标进行随机移位，得到新的原子坐标位置的结构，从而产生新的团簇结构；所述球切扰动是对已知的父代和母代两个团簇个体进行球域范围内外的交叉组合，从而得到新的子代团簇结构。

优选的，N的取值范围为5～24。

本发明相对于现有技术，引入了单金属稳定结构初始化、球切扰动和遗传局部优化等对传统BHMC算法进行改进，能有效处理合金团簇出现的同分异构体问题，避免陷入局部最优，从而更快速有效的寻找合金团簇最低能量稳定结构。

附图说明

图1a是合金团簇的结构优化的示意图；

图1b是合金团簇的元素原子排布优化的示意图；

图2是根据本发明实施例的Fe_nPt_m合金团簇的结构优化的流程图；

图3示出了N＝13时Fe原子与Pt原子各种比例情况下的初始最优结构；

图4是球切扰动处理过程的示意图；

图5示出了N＝24时不同比例情况下的Fe-Pt合金团簇稳定结构优化结果；

图6示出了原子总数N＝5-24，Fe:Pt在比例为(a)2:8,(b)5:5和(c)8:2时Fe-Pt合金团簇的最低能量构型。

具体实施方式

为进一步说明各实施例，本发明提供有附图。这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理。配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点。图中的组件并未按比例绘制，而类似的组件符号通常用来表示类似的组件。

现结合附图和具体实施方式对本发明进一步说明。参照图2，具体描述Fe_nPt_m合金团簇的结构优化步骤。改进BHMC算法主要包括如下步骤：

1、初始化种群

应用改进的BHMC算法进行Fe_nPt_m合金团簇的结构优化，涉及两个初始化过程：一是种群的初始化，得到整个算法流程的初始种群；二是遗传局部种群的初始化，得到较好结构下的随机元素排布的初始种群。

1)种群初始化：对于某个合金结构，通过在连续空间里初始原子的配置即在随机的球形容器{R}内，找到N个原子粒子的基态构型{R}＝r₁,r₂,r₃,...,r_N，其中

R_c＝(1/M)Σ_αR_c ^α，R_c ^α是元素α对应的共价半径^[20]，种群的大小会随着原子个数发生变化。在改进算法中，通过引入相同原子个数的Fe和Pt单金属的最稳定结构作为初始构型，从而避免完全随机方式无法产生所需的初始稳定结构。

通过上述初始化过程得到原子总数为N(N＝n+m)的Fe_nPt_m合金团簇在各种不同比例情况下的M个初始结构，M为种群中个体数。

2)遗传局部种群初始化：对原子总数为N的Fe_nPt_m二元合金团簇，在结合扰动算子得到较好结构的基础上，对原子位置进行相对固定，通过改变原子种类随机生成一个个体P，个体P的表示方式如下：

P＝[p₁，p₂，...，p_N]. (5)

其中，每个p_i的取值为0或者1，分别代表Fe原子或Pt原子，i表示第i个原子。p_i的取值还必须满足一个条件：

p₁+p₂+...+p_N＝n. (6)

该初始化过程仅用于遗传局部优化步骤，用于产生遗传局部优化的初始种群。

2.最速下降局部优化处理

通过对上述已经初始化的结构进行最速下降局部优化，可以得到初始后化的局部最优结构。最速下降局部优化处理的具体步骤是：

i)首先利用Gupta势能函数计算能量E，其中，原子总数为n的团簇的Gupta势能函数为，

其中r_ij表示原子i与原子j间的距离，r_ij ⁽⁰⁾表示两类原子间的特征长度，一般以块体材料的第一近邻距离表示，A_ij是衡量原子间的排斥强度的量，ξ_ij是有效跳跃积分，通常只与原子类别有关；

ii)再对Gupta势能函数(1)求导，计算得到梯度F；

iii)利用迭代公式X_i+1＝X_i-ρ_iF更新原子坐标位置，其中ρ_i是梯度方向上的步长，X_i是第i步的坐标变量，X_i+1是第i+1步的坐标变量；

iv)将X_i＝X_i+1重复步骤iii)不断迭代，直到满足前后两次迭代能量差ΔE<ε，其中ε为1e-6。

不断重复上述操作，可分别得到某个尺寸某个比例下基于Fe单金、Pt单金和原子数目较小的Fe-Pt合金团簇的初始最稳定结构。以原子总数N＝13为例，图3给出了Fe原子与Pt原子各种比例情况下的初始最优结构。图3中，深色表示Fe原子，白色表示Pt原子，n表示Fe原子数，m表示Pt原子数。

由图3可以发现，对于Fe_nPt_13-n，无论是单金属还是合金的最优结构，都十分类似。我们比较了其他原子数较小的单金属结构和合金团簇稳定结构，也得到了类似的结论。这说明在单金属最优结构的基础上再进行合金结构优化是合理的，并能保证获得比随机初始化更优的结构。

3.扰动处理

通过对每个局部最优结构的每个原子进行扰动处理，并进行最速下降局部优化，对比处理前后两个结构的能量，保存较低能量结构，在完成种群中M个个体扰动及最速下降局部优化处理后，可以得到M个能量较低的结构。在利用改进的BHMC算法进行Fe_nPt_m合金团簇能量优化的过程中，为避免陷入局部最优，我们引入了两种扰动方式：一种是原子移位扰动，即对原子坐标进行随机的位移，得到新的原子坐标位置的结构，从而产生新的团簇结构；另一种是球切扰动(参见图4)，即对已知的父代和母代两个团簇个体进行球域范围内外的交叉组合，从而得到新的子代团簇结构。前者通过小扰动可以实现局部寻优，后者通过增加大的扰动可以帮助跳出小的局部最优，利用这两种扰动方式进行团簇结构更新，不断搜索新的团簇结构，达到团簇结构寻优的目的。

4.遗传局部优化

合金团簇的结构优化比单金属团簇要复杂得多。在应用扰动算子时，不但要保证团簇的总原子数不变，还要保证团簇中每种元素的原子个数不变。此外，即使它们的几何构型极其相近，但两种原子的位置不同，就可产生出多种不同能量的团簇结构。特别是在两种合金原子数越接近1:1比例的情况下，这样的团簇就越多，这大大增加了寻找稳态结构的难度^]。对于原子总数为N的合金团簇Fe_nPt_m的单个同分异构体而言，其homotop异构构型的个数是

例如，就Fe₁₀Pt₁₀而言，它存在184756个homotop异构构型。

在基本的BHMC算法中，通常采用表面削弱和结构重组等扰动算子，以增强其跳出局部最优的能力。但这种扰动力度很小，很难保证在存在大量同分异构体的情况下跳出局部最优。为解决此类问题，我们在已有的BHMC算法基础上引入了遗传局部优化算子，以在较短的时间范围内找到每单个同分异构体的最优的homotop异构构型，从而大大节约在最优homotop异构构型上的求解时间。

遗传局部优化的具体操作如下：1)对经过扰动处理后得到的良好个体(即，结构优化后的最优结构)，令每个原子坐标保持不变，即固定每个原子位置；2)利用遗传局部种群初始化方法，通过改变每个位置上的原子类型(Fe或Pt)来产生在能量较低结构下的随机原子排布初始种群；3)利用遗传算法，分别对原子元素的分布进行交叉、变异和选择操作，以实现相对寻优，其中，交叉采用交叉率为0.4的两点交叉，变异采用概率为0.1的单点变异，选择采用轮盘赌方式；4)通过不断的迭代(重复3))找到在固定原子位置下的最优同构构型(即，最优元素原子分布)。

5.选择比较

对遗传局部优化得到的每个最优结构进行最速下降局部优化处理，保存更低能量结构为全局最优结构，并与上一代得到的全局最优结构进行能量对比，如果前后两次能量差在预定精度(10^-3)范围内，表示能量不变，在有预定代数(例如，250代)能量不变时或达到最大进化代数(例如，2000代)时，则寻优结束，否则回到扰动处理。

结果验证

基于本发明提出的改进BHMC算法，我们对N≤24的Fe_nPt_m合金团簇进行了结构寻优。图4示出了原子总数N＝24时不同比例情况下的Fe-Pt合金团簇稳定结构优化结果，深色表示Fe原子，白色表示Pt原子，其中n:m，n表示Fe原子数，m表示Pt原子数。从图4中可以看出合金团簇的稳定结构具有较好的对称性。同时，部分不同原子比例合金团簇之间有很强的结构相似性，如Fe₅Pt₁₉与Fe₁₂Pt₁₂。但也有部分不同比例合金团簇之间结构相似性较差，如Fe₁₉Pt₅与Fe₁₂Pt₁₂。

图5示出了原子总数N＝5-24，在Fe:Pt比例为(a)2:8,(b)5:5和(c)8:2时Fe-Pt合金团簇的最低能量构型，深色表示Fe原子，白色表示Pt原子，其中N(n:m)，N表示总原子数、n表示Fe原子数，m表示Pt原子数。由图6可以看出，Fe-Pt合金团簇稳定结构随总原子个数N的增加而呈现出结合成环的趋势，特别是在总原子数为N＝13时首次形成类二十面体的稳定结构。此外，在Fe原子个数占优的情况下，Fe-Pt合金团簇更容易成环，形成对称性更高的结构。这是因为在Fe原子数较多时，Pt原子会更加接近合金团簇稳定结构的质心，这样更加有益于形成稳定结构。

综上所述，本发明引入了单金属稳定结构初始化、球切扰动和遗传局部优化等对传统BHMC算法进行改进，能有效处理合金团簇出现的同分异构体问题，避免陷入局部最优，并且对计算机硬件开销要求较低，能够快速有效地实现对Fe_nPt_m合金团簇的结构寻优。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。

Claims (5)

1.一种基于改进BHMC算法的Fe_nPt_m合金团簇结构优化方法，其特征在于，包括以下步骤：

S1、初始化种群，得到原子总数为N的Fe_nPt_m合金团簇在各种不同比例情况下的M个初始结构，其中，N＝n+m，M为种群中个体数；

S2、对初始化的每个结构进行最速下降局部优化，得到初始化后的局部最优结构；

S3、对每个局部最优结构进行原子移位扰动和球切扰动处理，并进行最速下降局部优化，对比处理前后两个结构的能量，保存较低能量结构，在完成种群中M个个体扰动及最速下降局部优化处理后，可以得到M个能量较低的结构；

S4、对步骤S3得到的每个能量较低结构进行遗传局部优化处理，得到当前状态下最优结构的原子最优排布，其中，所述遗传局部优化处理包括：

S41、对步骤S3得到的能量较低结构，令每个原子坐标保持不变，即固定每个原子位置；

S42、利用遗传局部种群初始化方法，通过改变每个位置上的Fe或Pt原子类型来产生在能量较低结构下的随机原子排布初始种群；

S43、利用遗传算法，分别对原子排布进行交叉、变异和选择操作，以实现原子排布相对寻优，其中，交叉采用交叉率为0.4的两点交叉，变异采用概率为0.1的单点变异，选择采用轮盘赌方式；

S44、通过不断重复S43，找到在固定原子位置下的原子最优排布；

S5、对步骤S4得到的每个最优结构进行最速下降局部优化处理，保存更低能量结构为全局最优结构，并与上一代得到的全局最优结构进行能量对比，如果前后两次能量差在预定精度范围内，表示能量不变，在有预定代数能量不变时或达到最大进化代数时，则寻优结束，否则回到S3。

2.如权利要求1所述的方法，其特征在于，所述最速下降局部优化处理的具体步骤是:

1)首先利用Gupta势能函数计算能量E，其中，原子总数为n的团簇的Gupta势能函数为，

$V_n=\underset{i=1}{\overset{n}{\&Sigma;}}\left(V^r\right(i)-V^m(i\left)\right)---\left(1\right)$

$V^r\left(i\right)=\underset{j=1(j\&NotEqual;i)}{\overset{n}{\&Sigma;}}{A}_{ij}\exp (-p_{ij}(\frac{r_{ij}}{{r_{ij}}^{\left(0\right)}}-1)),---\left(2\right)$

$V^m\left(i\right)={\&lsqb;\underset{j=1(j\&NotEqual;i)}{\overset{n}{\&Sigma;}}{{\mathrm{\&xi;}}_{ij}}^2\exp \left(-2q_{ij}\left(\frac{r_{ij}}{{r_{ij}}^{\left(0\right)}}-1\right)\right)\&rsqb;}^{1/2}.---\left(3\right)$

其中r_ij表示原子i与原子j间的距离，r_ij ⁽⁰⁾表示两类原子间的特征长度，一般以块体材料的第一近邻距离表示，A_ij是衡量原子间的排斥强度的量，ξ_ij是有效跳跃积分，通常只与原子类别有关；

2)再对Gupta势能函数(1)求导，计算得到梯度F；

3)利用迭代公式X_i+1＝X_i-ρ_iF更新原子坐标位置，其中ρ_i是梯度方向上的步长，X_i是第i步的坐标变量，X_i+1是第i+1步的坐标变量；

4)将X_i＝X_i+1重复步骤3)不断迭代，直到满足前后两次迭代能量差ΔE<ε，其中ε为1e-6。

3.如权利要求1所述的方法，其特征在于，所述步骤S1的M个初始结构是通过引入相同原子个数的Fe和Pt单金属的最稳定结构作为初始构型来产生的。

4.如权利要求1所述的方法，其特征在于，所述步骤S3中的所述原子移位扰动是对原子坐标进行随机移位，得到新的原子坐标位置的结构，从而产生新的团簇结构；所述球切扰动是对已知的父代和母代两个团簇个体进行球域范围内外的交叉组合，从而得到新的子代团簇结构。

5.如权利要求1所述的方法，其特征在于，N的取值范围为5～24。