苗　竹　张　海　杨海瑞

(清华大学热能工程系，热科学与动力工程教育部重点实验室，北京100084)



升温烧结过程中TiO2纳米颗粒的原子分类分析(I)：表面原子识别

苗竹张海*杨海瑞

(清华大学热能工程系，热科学与动力工程教育部重点实验室，北京100084)

开发表面原子识别模型，对单个TiO2颗粒升温烧结过程中表面原子进行分类研究。模型对颗粒空间立方体网格化，利用标准球形颗粒体积积分法确定最佳网格尺寸为0.3 nm。通过表面网格识别实现表面原子分类，采用近邻网格中外部网格数量(Next)作为准则数判断目标网格是否为表面网格，确定最佳Next=9。基于LAMMPS软件模拟了半径为0.75 nm颗粒的升温过程，发现系统能量弛豫速度明显高于结构弛豫速度；利用表面识别模型分类分析原子特性，表面原子平均位移大于内部原子，且表面O原子迁移活性高于Ti原子；表面原子配位数低于内部原子，佐证表面结构规律性较差。研究结果为深入分析纳米材料活性位等结构分布奠定基础。

TiO2；纳米颗粒；分子动力学；表面识别；表面原子

www.whxb.pku.edu.cn

1　引言

二氧化钛(TiO2)是一种典型n型半导体金属氧化物1,2，其锐钛矿型晶相具有较高活性3,4，在气体传感器5－7、催化剂8,9和染料敏化太阳能电池10,11等方面得到广泛应用。近年来，人们发现火焰合成法制备TiO2纳米颗粒不仅经济，而且具有粒度分布可控、产物纯度高的特点12,13；同时，火焰合成法可以直接成膜制备气体传感器，具有明显优势14,15。

TiO2纳米颗粒的表面特性是影响传感器气敏感应性能的主要因素。气体分子与TiO2颗粒表面发生相互作用，颗粒表面相比内部区域存在更多结构缺陷，充当活性位用于气敏感应，颗粒表面特性分析对气体传感器性能研究有重要作用。分子动力学是研究纳米颗粒结构的有效方法16,17，可以直观确定颗粒结构在不同情况下的变化。Lu等18基于分子动力学方法研究了TiO2纳米颗粒之间的相互作用，直观观察纳米颗粒团聚融合过程中的结构变化；Mushnoori19和Tillotson20等分别通过分子动力学模拟TiO2表面与不同反应物作用过程，确定了反应物与表面活性位作用机制，Zhang等21利用分子动力学方法研究了纳米颗粒之间偶极相互作用对颗粒碰撞-烧结过程的影响。然而，现有研究主要基于纳米颗粒结构变化的直观观察，缺少对纳米颗粒表面结构的针对性分析。基于表面原子分类，可以进一步确定活性位等特征结构分布，对表征和分析纳米颗粒结构特性有重要作用。

颗粒表面结构分析的基础是表面原子识别，为此本文基于分子动力学原理建立了表面原子识别模型，对TiO2纳米颗粒表面原子和内部原子进行分类分析。空间网格化后以表面网格识别实现原子分类，通过近邻网格种类分析判断目标网格是否为表面网格。利用开源软件LAMMPS22,23模拟半径0.75 nm颗粒升温过程，通过表面原子识别模型对比分析升温过程中表面原子和内部原子迁移规律，确定颗粒表面特性。

2　模拟方法

本文计算模型为锐钛矿型TiO2，晶胞参数a= 0.37845 nm、c=0.95143 nm、Z=4，属于四方晶系24。构建半径0.75 nm的球形TiO2颗粒模型如图1所示，随机删除若干外侧原子保证颗粒整体电中性。图中粉色原子为钛(Ti)原子，绿色原子为氧(O)原子，颗粒中共有原子150个，其中Ti原子50个，O原子100个。

模拟通过开源软件LAMMPS完成，将TiO2纳米颗粒模型导入后模拟升温过程。模拟在正则系综(NVT)中进行，原子间相互作用势选用Matsui-Akaogi(MA)势24，该势函数被认为是刻画金属氧化物晶格结构的良好选择25。势函数形式为：

其中Uij为i粒子与j粒子之间的相互作用势，rij是i粒子与j粒子之间的距离，qi和qj分别是i粒子和j粒子的等效电荷，Ci和Cj分别是i粒子和j粒子之间的范德华力系数，f是一个标准力参数，其值大小为41.84 kJ·nm－1·mol－1，Ai和Aj分别为i粒子和j粒子的斥力半径，Bi和Bj是与泡利斥力有关的参数。势函数具体参数如表1所示。

系统边界条件选择非周期性(m,m,m)，计算时间步长为1 fs，每隔几纳秒输出中间结果，对系统势能、均方位移(MSD,mean-squared displacement)和配位数变化过程进行分析，其中MSD指颗粒中所有原子距离初始位置位移平方的平均值，用于描述系统结构弛豫过程；配位数指该原子周围晶格排布中配位原子个数，用于表征原子周围晶格排布周期性。

单颗粒升温模拟中，颗粒先在10 K下稳定1 ns，然后经过1 ns升温至设定温度，随后稳定弛豫约40 ns。在373－2373 K范围内选择19个计算点进行稳态计算，每隔几个纳秒输出计算结果用于过程分析。模拟完成后，分析系统势能和MSD变化规律，同时通过表面原子识别模型对表面原子进行识别和特性分析。

3　表面原子识别模型

表面原子识别模型旨在将表面原子和内部原子进行分类分析，以便对比分析两种原子的特性。由于颗粒弛豫之后原子坐标分布无规则，且不连续，直接处理原子坐标难度相对较大；为此，模型通过空间网格化后以网格为处理对象，根据网格内是否有原子将网格分类为内部网格和外部网格，再通过近邻网格种类确认目标网格的所处位置，实现表面网格识别，表面网格中所包含的原子即为表面原子。

3.1空间网格化

采用立方体网格对纳米颗粒所在空间进行网格化，理论上网格尺寸越小，原子类型识别的准确性越高，然而由于颗粒中原子分布不连续，网格尺寸过小时颗粒内部会出现无效空网格，干扰表面网格识别，如图2所示，因此网格尺寸存在最佳优化值。

图1　TiO2纳米颗粒模型Fig.1　TiO2nanoparticle model

表1　TiO2模拟的MA势函数参数表Table 1　Parameters of MApotential function for TiO2simulation

网格尺寸优化使用标准球形颗粒体积积分法。以半径2 nm颗粒为对象，在0.2－1.0 nm范围内选择九个值进行网格化，积分内部网格体积得到颗粒计算体积。如图3所示为颗粒计算体积随网格尺寸的变化关系。随网格尺寸增大，颗粒计算体积逐渐增大；相比颗粒体积理论值33.5 nm3，网格尺寸在0.3 nm时颗粒计算体积最为合理，即为最佳值。当网格尺寸小于0.3 nm时，内部无效空网格导致颗粒计算体积小于理论值；随网格增大，边界网格增大会导致更多外部体积计入颗粒内，导致颗粒计算体积随网格尺寸增大而增大。

综上，通过立方体网格实现纳米颗粒空间网格化，并利用标准球形颗粒体积积分法确定网格最佳尺寸为0.3 nm，为表面原子识别模型确定最佳网格化方案。

图2　不同尺寸网格效果示意图Fig.2　Schematic of the meshing results with different mesh sizes

图3　半径为2 nm的颗粒计算体积与网格尺寸关系Fig.3　Relation between calculated particle volume and mesh size with particle radius of 2 nm

3.2表面原子识别

空间网格化后，通过表面网格识别实现表面原子和内部原子分类，原理如图4所示，选定某目标网格，统计其近邻网格中外部网格数量(Next)，Next越大，该目标网格位于颗粒表面的可能性越高。选择不同若干Next值分别进行表面原子识别，对比识别效果，确定最佳Next值。

图5为四种Next选值情况下表面原子和内部原子识别效果图，统计得到表面原子与内部原子数之比与Next关系如图6所示。随Next增大，表面原子与内部原子数比逐渐减小，变化速度先减小后增大；在Next=9时，曲线变化率最小，表明此时表面原子识别效果最稳定，即最佳Next。对比半径为1、2和3 nm颗粒识别结果基本一致，进一步确定表面原子识别最佳的Next=9。最佳Next也对应与目标网格一个侧面的近邻网格数吻合，具有一定的理论依据。

综上，通过表面网格识别实现表面原子和内部原子分类，开发表面原子识别模型。通过近邻网格种类分析法，以近邻网格中外部网格数量Next为准则数判断目标网格是否位于颗粒表面，并确定最佳Next=9。

图4　表面网格识别方法示意图Fig.4　Schematic of surface mesh identification method

图5　不同Next情况下颗粒表面原子识别效果图Fig.5　Surface atoms identification result with different Nextvalues

4　结果分析

4.1颗粒升温过程

图6　不同半径(r)颗粒识别出表面原子与内部原子数之比与Next关系Fig.6　Relation between ratio of surface atom number to interior atom number and Nextfor particles with different radii(r)

通过LAMMPS模拟颗粒升温过程，各升温过程中温度变化曲线如图7所示。图8为半径0.75 nm颗粒升温过程系统势能随时间变化曲线，由图可知，系统势能随温度升高迅速增大，温度稳定后系统势能也随之稳定；温度越高，稳态势能越大。颗粒MSD随时间变化曲线如图9所示，对比势能变化曲线可看出，MSD变化速度明显慢于势能变化速度；温度稳定后，系统MSD仍缓慢增大而后趋于稳定。结果表明：系统势能弛豫速度明显快于系统结构弛豫速度，不能仅靠能量稳定判断系统是否达到稳态。当温度高于1173 K时，不同温度稳定MSD趋于一致，表明颗粒在高于此温度时会发生熔化现象，液相中原子运动性一致，对应一致的MSD结果。

4.2表面识别及特性分析

通过表面原子识别模型将颗粒中表面原子与内部原子分类，其中表面原子将Ti原子与O原子分类处理。统计得到内部原子、表面原子、表面Ti原子和表面O原子分别在573、973、1173、1573和2073 K情况下的原子平均位移，如图10所示。由图可知，所有原子平均位移均随温度升高而增大；在同一温度下，表面原子平均位移大于内部原子，且表面中O原子平均位移大于Ti原子，当温度较高颗粒熔化时Ti原子与O原子的位移相差不大；表明TiO2纳米颗粒表面原子运动强度高于内部原子，且表面O原子迁移活性强于Ti原子。对比表面原子和内部原子配位数如图11所示，随温度升高两类原子配位数均降低；各温度情况下，表面原子配位数低于内部原子配位数。

图7　不同升温过程温度变化曲线Fig.7　Temperature curves of different heating processes

图8　不同升温过程中势能随时间变化曲线Fig.8　Variation of potential energy with time in different heating processes

图9　不同升温过程中颗粒均方位移(MSD)随时间变化曲线Fig.9　Variation of particle mean-squared displacement(MSD)with time in different heating processes

图10　不同温度下各类原子平均位移对比Fig.10　Averaged displacement comparison of different kinds of atoms under different temperatures

图11　不同温度下表面原子与内部原子配位数对比Fig.11　Coordination number comparison of interior atoms and surface atoms under different temperatures

5　结论

基于分子动力学开发纳米颗粒表面原子识别模型，对颗粒中原子进行分类分析。采用立方体网格对颗粒空间网格化，利用标准球形颗粒体积积分法确定最佳网格尺寸为0.3 nm；通过近邻网格种类分析识别表面网格，以近邻外部网格数Next为准则数判断表面网格，并确定最佳Next=9。模拟半径0.75 nm纳米颗粒升温过程，分析系统势能和MSD变化规律，发现系统能量弛豫速度明显高于结构弛豫速度，不能通过能量稳定判断系统稳态。通过表面原子识别模型对比分析表面原子与内部原子特性，其中表面原子平均位移大于内部原子，且表面O原子迁移活性高于表面Ti原子；内部原子配位数整体高于表面原子。综上，表面原子识别模型效果稳定，能够准确分辨不同温度下表面原子和内部原子性质差异，可以用于进一步分析表面特征结构；该模型可以推广应用与其他晶体颗粒结构，但最佳网格尺寸与该晶体的晶格参数有关。

References

(1)Liu,M.;Li,H.M.;Wang,W.J.Catal.Today 2016,264,236. doi:10.1016/j.cattod.2015.08.044

(2)Bet-moushoul,E.;Mansourpanah,Y.;Farhadi,K.;Tabatabaei, M.Chem.Eng.J.2016,283,29.doi:10.1016/j.cej.2015.06.124

(3)Zhang,J.;Wu,B.;Huang,L.H.;Liu,P.L.;Wang,X.Y.;Lu,Z. D.;Wu,G.L.;Zhang,E.P.;Wang,H.B.;Kong,Z.;Xi,J.H.;Ji, Z.G.J.Alloy.Compd.2016,661,441.doi:10.1016/j. jallcom.2015.11.225

(4)Xu,L.;Tang,C.Q.;Huang,Z.B.Acta Phys.-Chim.Sin.2010, 26(5),1401.[徐凌,唐超群,黄宗斌.物理化学学报,2010, 26(5),1401.]doi:10.3866/PKU.WHXB20100526

(5)Zhao,M.;Wang,J.X.;Feng,C.H.;Zou,B.;Chen,C.;Wang,Z. Y.;Wu,F.Q.;Zou,L.H.Acta Phys.-Chim.Sin.2007,23(7), 1003.[赵萌,王金兴,冯彩慧,邹博,陈骋,王竹仪,吴凤清,邹乐辉.物理化学学报,2007,23(7),1003.]doi:10.3866/ PKU.WHXB20070708

(6)Moon,H.G.;Shim,Y.S.;Jang,H.W.;Kim,J.S.;Choi,K.J.; Kang,C.Y.;Choi,J.W.;Park,H.H.;Yoon,S.J.Sens.Actuat. B-Chem.2010,149,116.doi:10.1016/j.snb.2010.06.014

(7)Chaudhari,G.N.;Bende,A.M.;Bodade,A.B.;Patil,S.S.; Sapkal,V.S.Sens.Actuat.B-Chem.2006,115,297. doi:10.1016/j.snb.2005.09.014

(8)Ao,P.;Xu,X.S.;Xu,X.X.;Li,J.H.;Yan,X.H.Acta Phys.-Chim.Sin.2014,30(5),950.[敖平,许响生,徐潇潇,李加衡,严新焕.物理化学学报,2014,30(5),950.] doi:10.3866/PKU.WHXB201403111

(9)Marin,R.P.;Ishikawa,S.;Bahruji,H.;Shaw,G.;Kondrat,S. A.;Miedziak,P.J.;Morgan,D.J.;Taylor,S.H.;Bartley,J.K.; Edwards,J.K.;Bowker,M.;Ueda,W.;Hutchings,G.J.Appl. Catal.A 2015,504,62.doi:10.1016/j.apcata.2015.02.023

(10)Gao,S.W.;Lan,Z.;Wu,W.X.;Que,L.F.;Wu,J.H.;Lin,J. M.;Huang,M.L.Acta Phys.-Chim.Sin.2014,30(3),446.[高素雯,兰章,吴晚霞,阙兰芳,吴季怀,林建明,黄妙良.物理化学学报,2014,30(3),446.]doi:10.3866/PKU. WHXB201401022

(11)Guo,W.;Wang,K.;Shen,Y.H.;Zhang,H.;Weng,T.;Ma,T.L. Acta Phys.-Chim.Sin.2013,29(1),82.[郭薇,王开,沈艺华,张贺,翁韬,马廷丽.物理化学学报,2013,29(1),82.] doi:10.3866/PKU.WHXB201211071

(12)Zhang,Y.Y.;Li,S.Q.;Deng,S.L.;Yao,Q.;Tse,S.D.J.Aeros. Sci.2012,44,71.doi:10.1063/1.3555633

(13)Zhao,B.;Uchikawa,K.;Wang,H.Proc.Combust.Inst.2007, 31,851.doi:10.1016/j.proci.2006.08.064

(14)Tolmachoff,E.;Memarzadeh,S.;Wang,H.J.Phys.Chem.C 2011,115,21620.doi:10.1021/jp206061h

(15)Zhao,B.;Uchikawa,K.;McCormick,J.R.;Ni,C.Y.;Chen,J. G.;Wang,H.Proc.Combust.Inst.2005,30,2569.doi:10.1016/ j.proci.2004.08.146

(16)Fereidoon,A.;Aleaghee,S.;Taraghi,I.Comp.Mater.Sci.2015, 102,220.doi:10.1016/j.commatsci.2015.02.044

(17)Gautam,S.;Cole,D.Chem.Phys.2015,458,68.doi:10.1016/j. chemphys.2015.07.012

(18)Lu,J.;Liu,D.M.;Yang,X.N.;Zhao,Y.;Liu,H.X.;Tang,H.; Cui,F.Y.Appl.Surf.Sci.2015,357,1114.doi:10.1016/j. apsusc.2015.09.142)

(19)Mushnoori,S.;Chong,L.;Dutt,M.Materials Today: Proceedings 2016,3,513.doi:1016/j.matpr.2016.01.082

(20)Tillotson,M.J.;Brett,P.M.;Bennett,R.A.;Ricardo G.C.Surf. Sci.2015,632,142.doi:10.1016/j.susc.2014.09.017

(21)Zhang,Y.Y.;Li,S.Q.;Yan,W.;Yao,Q.;Tse,S.J.Chem.Phys. 2011,124,084501.doi:10.1063/1.3555633

(22)Freitas,R.;Asta,M.;Koning,M.D.Comp.Mater.Sci.2016, 112,333.doi:10.1016/j.commatsci.2015.10.050

(23)Martin,M.G.;Thompson,A.P.Fluid Phase Equilib.2004,217, 105.doi:10.1016/j.fluid.2003.06.007

(24)Zhang,Y.Y.Dynamics of Nanoparticles in Stagnation Flames. Ph.D.Dissertation.Tsinghua University,Beijing,2013.[张易阳.基于滞止火焰合成的高温纳米颗粒动力学研究[D].北京:清华大学,2013.]

(25)Zhang,H.Z.;Chen,B.;Banfield,J.F.Phys.Chem.Chem.Phys. 2009,11(14),2553.doi:10.1039/b819623k

Atom Identification and Analysis of TiO2Nanoparticles in the Heating and Sintering Process(I):Surface Atom Identification

MIAO ZhuZHANG Hai*YANG Hai-Rui
(Key Laboratory for Thermal Science and Power Engineering of Ministry of Education,Department of Thermal Engineering, Tsinghua University,Beijing 100084,P.R.China)

In this paper,we develop a surface atom identification model to perform atom identification and analysis of the surface atoms of a single TiO2nanoparticle during the heating and sintering process.Cubic mesh was used to obtain the particle structure mesh,and the optimal mesh size of 0.3 nm was determined by volumetric integration of a spherical particle.Surface atoms were classified according to identification of surface meshes,and the number of external meshes in all of the neighbor meshes(Next)was set as a criterion to determine whether the target mesh was a surface mesh.The optimal value was Next=9.LAMMPS was used to simulate the heating process of a particle with a radius of 0.75 nm.The results show that energy relaxation is significantly faster than structure relaxation.The atom classification analysis with the developed surface atom identification model shows that displacement of the surface atoms is larger than the interior atoms,and surface O atoms are more active in migration than surface Ti atoms.The coordination number of surface atoms is lower than that of interior atoms.The present study provides fundamental information for analyzing the active structure distribution of nanoparticles.

TiO2;Nanoparticle;Molecular dynamics;Surface identification;Surface atom

March 14,2016;Revised:April 26,2016;Published on Web:April 26,2016.

O647

10.3866/PKU.WHXB201604262

*Corresponding author.Email:haizhang@tsinghua.edu.cn;Tel:+86-10-62773153.

The project was supported by the National Natural Science Foundation of China(51176095,51476088).

国家自然科学基金(51176095,51476088)资助项目

©Editorial office ofActa Physico-Chimica Sinica

［Article］