吴夏1)2) 魏征1)2)

1)(安庆师范大学化学化工学院,安庆 246011)

2)(安庆师范大学,光电磁功能材料安徽省重点实验室,安庆 246011)

1 引 言

在纳米团簇的基础研究及应用领域,三元金属团簇因其特殊的光、电、磁和催化等性质成为近期研究热点[1−4].针对单一金属及二元合金团簇已经开展了广泛的实验与理论研究[5−7].与二元合金团簇和单一金属团簇相比,三元金属团簇拥有了更高的催化活性和选择性[8−10].例如,与Pt-Ru二元团簇相比,Pt-Ru-Co纳米粒子对于甲醇氧化反应表现出更强的催化性能[9].对于乙醇氧化反应,与Pd/C催化剂相比,Cu1Au0.5Pd1.5/C具有更高的活性和稳定性[8].然而,三元金属团簇的性质易受到尺寸和原子序列变化的影响.对于实验合成的Au23Cu1Pd(SC12H25)18团簇,Pd原子的存在增强了团簇的稳定性[1].理论计算方面使用量子修正的Sutton-Chen(SC)势对二十四面体Cu-Au-Pt团簇的稳定结构进行了系统研究,发现Cu或Au的比例影响到表面偏聚性能[3].另一团簇结构优化结果显示75原子Ag-Pd-Pt团簇的主要构型是十面体/密堆积的混合结构,而相应的Au-Pd-Pt团簇为Marks十面体结构,且在结构的外层Au原子比Ag原子形成更为聚集的结构[10].

在分子模拟领域,优化二元及更多元异种类型组成的团簇稳定结构是一项艰巨的任务.对于单一原子类型A团簇,确定其稳定几何结构的全局优化属于非确定性多项式难题.已经广泛应用的优化算法包括遗传算法(GA)[11],BH(basin-hopping)算法[12]、动态格点搜索(DLS)算法[13]、自适应免疫优化算法(AIOA)[14]等.而在二元AmBn团簇中,除了单一原子类型中存在的几何构型同分异构体,还存在由不同类型原子相对位置差异引起的homotop同分异构体.任一几何构型均存在(m+n)!/(m!n!)个homotop异构体,如原子数为10的A5B5团簇homotop异构体数目为252.优化这类团簇时常采用GA,算法等BH基本框架,并普遍使用原子交换的方法解决homotop问题[15−17].而在三元团簇AmBnCl结构中任一几何构型的homotop数目为(m+n+l)!/(m!n!l!),如原子总数为10、组成为A5B2C3的homotop数目为2520.可见,从优化单一原子类型团簇到二元团簇,homotop问题使得优化过程变得尤为困难.因此,对于存在更多homotop异构体的三元团簇结构优化来说,亟待发展更为高效的全局优化算法.

Northby[18]提出了静态建模的思想,在优化单一原子类型Lennard-Jones(LJ)团簇结构时设计了二十面体的格点,并使用GA进行格点搜索得到了LJ147团簇最稳定结构.后来,采用相同的思想将团簇可优化的尺寸扩大到1610原子[19].受该思想启发,为了减小搜索空间,提出了一种内核构建(IC)的办法,即在产生初始结构时构建一个特定构型的内核,如Ino十面体、Mackay二十面体等.Shao教授课题组[20]提出了基于IC的动态格点搜索(DLSc)算法成功优化到LJ670团簇.对于原子数目小于200的团簇优化,DLSc算法是无偏的,即最稳定结构不受初始构型影响.此外,DLSc算法也成功用于优化Ag120—310团簇[21].

本文提出了基于IC的自适应免疫优化算法(称为AIOA-IC算法)来优化三元Cu-Au-Pd团簇最稳定结构.为了测试AIOA-IC算法的效率,将其运用于原子数为60的Ag-Pd-Pt团簇结构优化.结果发现,最稳定结构的能量均低于文献报道的值.使用改进的算法分别优化了原子数为38及55的三元Cu-Au-Pd团簇稳定结构,并分析了结构生长以及原子分布规律.此外,研究发现147原子的Cu12Au93Pd42团簇具有完整的核层二十面体结构.

2 计算原理与方法

2.1 Gupta势函数

目前,Gupta势函数[22]、Murrel-Mottram势函数[23]、SC势函数[3]以及嵌入原子势[24]常被用于描述金属团簇原子间相互作用以确定团簇的稳定结构及热力学性质.其中,Gupta势函数已被用于描述单一金属团簇、二元合金团簇、三元甚至于四元金属团簇原子间的相互作用[25].本研究采用Gupta势函数描述Cu-Au-Pd团簇中同种及异种类型原子间相互作用.此势函数的形式可以分解为斥力项Vr(i)和引力项Vm(i),对于总原子数为N的Gupta势函数VN表述为

其中,rij代表原子i和j之间的距离,r(0)ij是最近邻距离,Aij是斥力项系数,ξij代表原子i和j之间的有效跳跃积分参数,Pij与qij分别描述为对斥力与引力作用的权重.Cu-Au-Pd团簇的Gupta势函数采用实验值拟合的参数引自文献[26,27],列于表1中.

表1 三元Cu-Au-Pd团簇Gupta势函数参数Table 1.Gupta potential parameters of trimetallic Cu-Au-Pd clusters.

2.2 IC

在优化大尺寸单一原子类型团簇结构时,采用固定内核的方法极大地提高了DLS的效率.固定内核是指在产生初始结构时使用十面体(Dh)、二十面体(Ih)、面心立方(fcc)等结构作为其内核,而多余的原子则随机地产生在内核的外层.实验证明通过构建内核初始结构变得更加规则,减小了搜索空间,提高了方法效率.研究结果已证实,对于优化原子数为200以内的团簇结构时,内核的构建并不会改变优化算法的无偏性.

根据文献[28]报道的单一原子类型及多元团簇常见的结构,选用十面体、二十面体、fcc、六折叠(six-fold)、Leary四面体(LT)结构作为初始结构的内核,如图1所示.十面体与二十面体同属于五折叠(f i ve-fold)结构.图1中使用不同颜色圆球来区分Ih55,six-fold38,six-fold55,LT34和LT98结构的内、外层原子.针对待优化原子数目来选择使用相应尺寸的内核.若优化的原子数目小于55时,从Dh13,Ih13,fcc19,LT34,six-fold38中随机选择作为内核.当优化原子数目小于98且大于等于55时,被选用的内核包括Dh55,Ih55,fcc44,LT34和six-fold51.对于原子数大于等于98的团簇,选择Dh55,Ih55,fcc85,LT98和six-fold51作为内核.当构建好内核以后,随机地在内核的外层产生余下原子的坐标,构成初始结构.

图1 (网刊彩色)典型的十面体(Dh)、二十面体(Ih)、面心立方体(fcc)、六折叠(six-fold)和Leary四面体(LT)结构Fig.1.(color online)Typical structures of decahedron(Dh),icosahedron(Ih),face centered cubic(fcc),six-fold and Leary tetrahedron(LT).

2.3 AIOA-IC算法

受生物免疫系统理论的启发,AIOA算法模拟了抗体的产生、克隆与变异等过程.AIOA算法基本步骤包括免疫克隆和变异操作.用于团簇结构优化的AIOA算法首先随机地产生初始结构及初始基因库.然后,采用免疫选择程序从初始基因库中选择一定数量的个体,再对这些个体采取免疫操作以产生新个体.最后,更新基因库操作以获取能量更低的结构.不断重复选择、变异和更新基因库的操作,通过多次迭代以搜索团簇最稳定结构.AIOA算法及其改进算法已经成功应用于LJ200单一原子类型团簇[29]、Cu-Ag和Cu-Au等二元团簇[30]、及Ar-Kr-Xe[31]和Cu-Ag-Au等[32]三元团簇的结构优化.AIOA-IC算法采用AIOA算法的基本框架,其具体流程如下.1)为了更大程度地保持算法的无偏性,除了采取上述构建内核的方式产生初始结构,还保留随机创建初始结构的方法.随机地产生Nlib个初始结构,并使用限制内存的类牛顿算法进行局部优化操作[33].这些结构构成了初代(it=0)的初始基因库.2)免疫克隆操作.计算基因库中第it代个体的适应值(依据Gupta势函数计算的势能量)与个体密度的比值来确定个体的选择可能性概率.个体密度则由被选择的个体与基因库中所有个体的相似性来判断.并通过免疫选择程序从基因库中挑选出Npop个体.具体计算参考AIOA[14].3)变异操作.针对多元团簇存在的几何构型同分异构与homotop同分异构问题分别采取不同的变异方法.为了得到更优的几何构型,50%的概率用来选择个体中能量高的原子并将其随机地移动到团簇表面.另外的50%概率则采用原子交换的方式用于解决homotop难题[30].原子交换方式是通过随机地选择一对异种类型的原子再对换其原子类型来实现的.4)更新基因库.对比变异获得的新结构与基因库中的所有结构,采用基于连接表的相似性检测程序[29]更新基因库以获得能量更低的结构.拥有更低能量及更少相似性的结构则保留在基因库中.5)此时,it=it+1.若达到既定的循环次数,输出能量最低的结构即为全局最优结构,否则,程序转到步骤2.为充分搜寻三元合金团簇对应的势能面,AIOA-IC算法运行100次.

3 结果与讨论

3.1 AIOA-IC算法效率

以三元团簇结构优化为例来测试AIOA-IC全局优化算法的效率.目前,多元团簇中仅有二元LJ团簇模型可被用于测试二元团簇结构优化算法的效率,而三元团簇的计算模拟尚缺少多种优化算法相互检验以建立标准模型.本研究以文献[10]报道的Ag20PdnPt40−n(n=5,10,15,20,25,30,35)团簇结构作为优化对象,采用与文献[10]一致的优化算法参数以及相同的算法运行次数.

图2(a)显示了Ag20PdnPt40−n团簇文献[10]报道的局部极优,图2(b)显示了本文得到的势能量最优结构.从团簇势能量(E)对比结果可以看出AIOA-IC算法优化得到的结构能量均低于文献值.比如,Ag20Pd10Pt30团簇的势能量值差异为0.32 eV.观察最优结构可知Ag20Pd5Pt35,Ag20Pd10Pt30,Ag20Pd15Pt25,Ag20Pd20Pt20和Ag20Pd25Pt15团簇采取规则的五折叠饼状结构类型,而Ag20Pd30Pt10与Ag20Pd35Pt5团簇呈现出基于55原子完整二十面体生长的构型.然而,对于文献[10]中报道的局部最优结构,Ag20Pd5Pt35为无定型结构,Ag20Pd25Pt35为由三个19原子双二十面体面面相连组成的环状结构,Ag20Pd30Pt10为六折叠饼状结构.研究显示作为主要构型的五折叠饼状结构其内核为13原子的十面体.可见十面体内核的构建易于形成更为稳定的五折叠结构.因此,AIOA-IC算法显著地提高了AIOA算法的效率.

图2 (网刊彩色)(a)文献[10]报道及(b)本文得到的能量更低的Ag20PdnPt40−n团簇结构Fig.2.(color online)The structures and energies of the investigated Ag20PdnPt40−n(n=5,10,15,20,25,30 and 35)clusters reported(a)from references and(b)in this work.

3.2 38原子Cu-Au-Pd最稳定结构

本文使用AIOA-IC方法优化得到了38原子Cu6AunPd32−n,CunAu6Pd32−n及CunAu32−nPd6(n=1—31)团簇最稳定结构.图3显示了随着n的变化,Cu-Au-Pd团簇势能量的变化规律.图中绘制了部分团簇的结构图形,使用不同颜色圆球分别代表Cu,Au和Pd原子.

如图3(a)所示,Cu6AunPd32−n团簇结构可被分为五折叠结构和六折叠结构.当Au原子比例低时,即n值小,团簇主要是五折叠结构;而n＞14时,拥有六折叠结构.图3(b)所示的CunAu6Pd32−n团簇包括五折叠结构和完整的截角八面体(TO)结构.当n=3—14时为五折叠结构,其他均为TO结构.TO构型常出现在38原子单一原子类型及二元团簇中[16].在图3(c)所示的CunAu32−nPd6团簇中, 当n=1—3和13—16时形成五折叠结构,当n=4—12时为六折叠结构,而当n=16—31时形成TO结构.

此外,由图3(a)可知,Cu6AunPd32−n团簇中6个Cu原子均位于结构的内层.在CunAu6Pd32−n(图3(b))与CunAu32−nPd6(图3(c))团簇中,随着n的增加,即Cu原子数目的增加,Cu原子均从内层向外层生长.因此,Cu原子趋于位于结构的内层.从图3(a)可知,Cu6AunPd32−n团簇中随着Au原子数目的增加,绝大部分Au原子都占据着结构的外层.如图3(b)所示的CunAu6Pd32−n团簇中,6个Au原子均位于最外层,但随着Cu和Pd组成的变化形成不同的分布规律.另外,从图3可以看出,Pd原子趋于分布在结构的外层.

图3 (网刊彩色)(a)Cu6AunPd32−n,(b)CunAu6Pd32−n 及(c)CunAu32−nPd6(n=1—31)团簇势能量变化及结构分布,Cu,Au和Pd原子分别用绿色、红色、白色圆球表示Fig.3. (color online)Variations of potential energies and structures of(a)Cu6AunPd32−n,(b)CunAu6Pd32−nand(c)CunAu32−nPd6(n=1–31)clusters.

3.3 55原子Cu-Au-Pd团簇结构

采用AIOA-IC方法优化了55原子三元Cu-Au-Pd团簇稳定结构.图4显示了Cu13AunPd42−n,CunAu13Pd42−n及CunAu42−nPd13(n=1—41)团簇最优结构.对于所研究的55原子Cu-Au-Pd团簇,所有构型都拥有Mackay二十面体结构.如图4(a)所示的Cu13AunPd42−n团簇中,13个Cu原子均位于团簇的内层.随着n值增加,增加的Au原子均逐渐生长在结构的外层.此外,对于所有组成,Pd原子均分布在外层.对于n=1—17,13个Cu原子构成二十面体的内核.而当n=18—41时,结构中心被1个Au原子占据,次外层为12个Cu原子,另一个Cu原子位于表面.

图4 (网刊彩色)代表性的(a)Cu13AunPd42−n,(b)CunAu13Pd42−n,(c)CunAu42−nPd13(n=1—41)团簇稳定结构Fig.4.(color online)Typical structures of(a)Cu13AunPd42−n,(b)CunAu13Pd42−nand(c)CunAu42−nPd13(n=1–41)clusters.

图4 (b)显示了CunAu13Pd42−n团簇的稳定结构.Cu1Au13Pd41团簇的1个Cu原子位于中心,13个Au原子分散在结构外层.其中,二十面体的12个顶点被12个Au原子占据,另一个Au原子也位于外层.而随着Cu原子数增加,团簇中Cu原子先占据内层,当内层均被占据之后,多余的Cu原子分布到外层.13个Au原子均分布在外层,但呈现出不同的分布方式.因此,Au原子的分布受团簇中Cu和Pd组成变化的影响.图4(c)显示了部分CunAu42−nPd13团簇稳定结构.当n=1时,内核是13原子二十面体,其中心被Cu原子占据,Cu原子被12个Pd原子包裹,另一个Pd原子则位于外层.随着Cu原子数目的增加(n值增加),增加的Cu原子先占据内层,而相应的Pd原子被排挤到外层.对于所有n值,Au原子均分布在外层.

使用序列参数(R)来描述团簇中原子分布规律,R值为不同类型原子到团簇结构中心的距离,具体形式如下:其中,xi,yi,zi为原子坐标值.序列参数R值偏大或偏小分别说明该类型原子趋于分布在团簇表面或内部,中等大小的值则说明可能形成混合模式.

图5显示了Cu13AunPd42−n,CunAu13Pd42−n和CunAu42−nPd13团簇序列参数R值.如图5(a)所示的Cu13AunPd42−n团簇中,RCu的值小,而RAu和RPd值大.说明Cu原子位于内层,Au和Pd原子趋于分布在外层.当n=1—17时,RCu值保持不变,当n=18—41时,RCu拥有一定值,且该值略大于前者.如上所述,当n=1—17时,二十面体的内核被13个Cu原子占据.而当n=18—41时,位于内核中心的Cu原子被Au原子替换,被替换的Cu原子被放置到结构表面.图5(b)所示的CunAu13Pd42−n和图5(c)所示的CunAu42−nPd13团簇中,随着n值增加,RCu值由小逐渐变大,这是由于增加的Cu原子逐渐由内层向外层生长.这也证实了Cu原子趋于占据内层的结论.此外,对于CunAu42−nPd13团簇(图5(c)),13个Pd原子的RPd值由中等值增加到较大的值.这是由于当Au原子数比例高,Au更趋于占据最外层,而Pd趋于位于次外层,因而形成混合模式.随着Au原子比例的不断降低,Pd和Au共同位于外层.

图5 (网刊彩色)(a)Cu13AunPd42−n,(b)CunAu13-Pd42−n 与(c)CunAu42−nPd13(n=1—41)团簇序列参数变化Fig.5.(color online)Variations of the order parameter values of(a)Cu13AunPd42−n,(b)CunAu13Pd42−n and(c)CunAu42−nPd13(n=1–41)clusters.

3.4 147原子Cu-Au-Pd二十面体结构

为了更清楚地考察三元合金团簇中各种原子的分布规律,使用AIOA-IC方法优化原子数高达147原子的Cu-Au-Pd团簇稳定结构.如前所述,三元团簇结构势能面异常复杂,确定其稳定结构需要的计算非常耗时.因此,重点考察了具有代表性的团簇结构.如图6显示的Cu12Au93Pd42团簇为完整的二十面体结构.图6(b)中使用不同颜色的圆球来表示Cu,Au和Pd原子符号;图6(a)为团簇结构的直观图,外层Au原子用小圆球来表示.

众所周知,二十面体可被看成由二十个四面体构成,图6(b)显示了其中一个四面体结构单元.图6(b)中,单个Au原子代表图6(a)结构的内核原子,Cu原子构成了内层,Pd原子构成次外层,最外层则全部被Au原子所占据.图6(c)显示了图6(a)结构的中心剖面图.可以清楚地看出,Cu12Au93Pd42团簇中,中心原子为Au,12个Cu原子包裹着中心原子构成了内层,42个Pd原子占据了结构的次外层(同时也是55原子二十面体内核的外层).此外,147原子完整二十面体的最外层共有92个格点,Cu12Au93Pd42团簇的最外层完全被92个Au原子占据.从大尺寸的团簇结构可以看出,对比Au与Pd原子的分布规律,Au原子更趋于分布在最外层,Pd原子分布在次外层.

图6 (网刊彩色)(a)Cu12Au93Pd42团簇二十面体结构及能量;(b)四面体结构单元;(c)为(a)结构的截面图Fig.6.(color online)(a)The complete icosahedral structure and energy of Cu12Au93Pd42;(b)its tetrahedral unit;and(c)sectional view of the icosahedral structure.

Cu,Au和Pd的分布规律可由原子半径和表面能来解释.Cu,Au和Pd的原子半径分别为1.28,1.44和1.37 Å,表面能分别为113.9,96.8和125—131 mV·Å−2[34].Cu原子位于团簇内层的原因可由三种原子中最小的原子半径及较大的表面能来解释.Au原子半径最大,表面能最小,易于生长在外层.与上述分析结果一致.另外,Pd原子半径大,说明Pd原子易于被推到结构外层.而Pd原子表面能大,这点又不支持Pd原子位于外层.因此,理论上大尺寸Cu-Au-Pd团簇中Pd原子可能会位于内层或外层.本研究计算结果显示Pd原子更趋于分布在次外层.

4 结 论

针对全局优化算法在确定大尺寸三元合金团簇结构时优化难的问题,提出了基于IC的方法改进了自适应免疫优化算法(AIOA),称之为AIOAIC方法.首先运用该算法优化了60原子Ag-Pd-Pt稳定结构,发现了比文献[10]中报道的能量更低的结构.使用新算法优化了38原子和55原子Cu-Au-Pd团簇最稳定结构.结果显示38原子Cu-Au-Pd团簇结构类型包括五折叠、六折叠及截角八面体结构,且结构类型受到原子组成的影响.而55原子Cu-Au-Pd团簇均为二十面体结构.此外,147原子Cu12Au93Pd42团簇中,中心原子为Au,内层由12个Cu原子构成,次外层由42个Pd原子构成,最外层则被92个Au原子占据.序列参数显示大尺寸团簇中,Cu原子趋于占据内层,Pd和Au原子分别倾向于分布在次外层和最外层.

[1]Sharma S,Kurashiqe W,Nobusada K,Neqishi Y 2015Nanoscale7 10606

[2]Zhang M,Zhang J F,Gu T,Zhang H Y,Luo Y H,Cao W 2015J.Phys.Chem.A119 3458

[3]Li T J,Sun Y,Zheng J W,Shao G F,Liu T D 2015Acta Phys.Sin.64 153601(in Chinese)[李铁军,孙跃,郑骥文,邵桂芳,刘暾东2015物理学报64 153601]

[4]Ma Z N,Jiang M,Wang L 2015Acta Phys.Sin.64 187102(in Chinese)[马振宁,蒋敏,王磊 2015物理学报64 187102]

[5]Sattler K,Mühlbach J,Recknagel E 1980Phys.Rev.Lett.45 821

[6]Ferrando R,Jellinek J,Johnston R L 2008Chem.Rev.108 845

[7]Meitzner G,Via G H,Lytle F W,Sinfelt J H 1985J.Chem.Phys.83 4793

[8]Mao H,Huang T,Yu A S 2014J.Mater.Chem.A2 16378

[9]Zhang X,Zhang F,Chan K Y 2004Catal.Commum.5 749

[10]Wu X,Liu Q M,Sun Y,Wu G H 2015RSC Adv.5 51142

[11]Deaven D M,Tit N,Morris J R,Ho K M 1996Chem.Phys.Lett.256 195

[12]Wales D J,Doye J P K 1997J.Phys.Chem.A101 5111

[13]Cai W S,Shao X G 2002J.Comput.Chem.23 427

[14]Shao X G,Cheng L J,Cai W S 2004J.Chem.Phys.120 11401

[15]Johnston R L 2003J.Chem.Soc.Dalton Trans.22 4193

[16]Doye J P K,Meyer L 2005Phys.Rev.Lett.95 063401

[17]Wu X,Liu Q M,Duan R Y,Wei Z 2016Acta Phys.Sin.65 210202(in Chinese)[吴夏,刘启满,段仁燕,魏征 2016物理学报65 210202]

[18]Northby J A 1987J.Chem.Phys.87 6166

[19]Xiang Y H,Cheng L J,Cai W S,Shao X G 2004J.Phys.Chem.A108 9516

[20]Yang X L,Cai W S,Shao X G 2007J.Comput.Chem.28 1427

[21]Shao X G,Yang X L,Cai W S 2008Chem.Phys.Lett.460 315

[22]Gupta R P 1981Phys.Rev.B23 6265

[23]Cleveland C L,Landman U,SchaaffT G,Shaf i gullin M N,Stephens P W,Whetten R L 1997Phys.Rev.Lett.79 1873

[24]Mantina M,Valero R,Truhlar D G 2009J.Chem.Phys.131 064706

[25]Wu X,Wei Z,Liu Q M,Pang T,Wu G H 2016J.Alloy Compd.687 115

[26]Darby S,Mortimer-Jones T V,Johnston R L,Roberts C 2002J.Chem.Phys.116 1536

[27]Ismail R,Johnston R L 2010Phys.Chem.Chem.Phys.12 8607

[28]Rossi G,Ferrano R,Rapallo A,Fortunelli A,Curley B C,Lloyd L D,Johnston R L 2005J.Chem.Phys.122 194309

[29]Cheng L J,Cai W S,Shao X G 2004Chem.Phys.Lett.389 309

[30]Wu X,Cai W S,Shao X G 2009J.Comput.Chem.30 1992

[31]Wu X,Sun Y,Gao Y C,Wu G H 2013J.Mol.Model.19 3119

[32]Wu X,Wu G H,Chen Y C,Qiao Y Y 2011J.Phys.Chem.A115 13316

[33]Liu D C,Nocedal J 1989Math.Program45 503

[34]Wu X,Sun Y,Wei Z,Chen T J 2017J.Alloy Compd.701 447