下载文档原格式
团簇物理研究王丽丽摘要：团簇物理在发展过程中，从原子核物理、凝聚态物理和原子分子物理等学科引入了许多概念和方法，构成团簇研究的中心议题，逐渐形成一门介于原子分子物理与凝聚态物理之间的交叉学科。
文章对团簇物理作了简单介绍，从团簇的概念、发展史、研究范畴到它的性质、研究方法。
作为一个初学者，利用gussian03程序包对二、三小原子团簇的结构进行了计算，算出了键长与键角并分析了它们的结构。
关键词：团簇，gussian03，团簇结构，密度泛函理论（DFT）模型STUDY ON CLUSTER PHYSICSAbstract: In the development of claster physics many concepts and methods have been introduced from atomic and molecular physics ,nuclear physics and matter physics , forming an interdisciplinary field between atomic and molecular physics on the one hand and condensed matter physics on the other .Key words: claster ,gaussian03, structure of claster, DFT1、引言团簇的研究开始于上个世纪七十年代，到了八十年代有了较大的发展。
由于团簇的知识构成的特殊性，即它是从原子分子物理、凝聚态物理、表面科学、量子化学、材料科学，甚至核物理学引入了概念和方法，构成其知识框架。
所以团簇物理是一门交叉学科，它的研究需要掌握原子分子、量子化学、凝聚态、电子计算机技术等一系列的知识。
在前一阶段的调研中我阅读了王广厚的《团簇物理学》，这本书较全面的介绍了团簇的合成、结构、性质等，还有四川大学毛华平的博士学位论文《金、铜、钇小团簇的几何结构、势能函数、能级分布和电子特性研究》等一些文章涉及到团簇的性质的研究，对于结构的计算只是简单介绍而没有具体阐述。
MaterialsStudio软件一、Materials Studio软件的主要应用领域包括：金属材料研究无机非金属材料研纳米材料研究高分子及其复合材料研究表界面研究化学反应研究含能材料研究生物、医药研究在晶体结构、形貌研究中的应用QSAR 的应用Perl 语言的应用Accelrys(美国)公司是世界领先的计算科学公司，是一系列用于科学数据的挖掘、整合、分析、模建与模拟、管理和提交交互式报告的智能软件的开发者，是目前全球范围内唯一能够提供分子模拟、材料设计、化学信息学和生物信息学全面解决方案和相关服务的软件供应商，所提供的全面解决方案和科技服务满足了当今全球领先的研究和开发机构的要求。
Materials Studio多尺度分子模拟平台是Accelrys公司（美国）在材料设计领域的核心产品。
它融合多种模拟方法，整合多达23 个功能模块，实现从电子结构解析到宏观性能预测的全尺度科学研究。
在国内拥有近400家用户，分布在石油、化工、环境、能源、制药、电子、食品、航空航天和汽车等工业领域和教育科研部门；相关的研究工作在Nature、Science等各类权威期刊上发表论文过万篇。
Materials Studio分子模拟软件采用了先进的模拟计算思想和方法，如量子力学（QM）、线性标度量子力学(Linear Scaling QM)、分子力学(MM)、分子动力学(MD)、蒙特卡洛(MC)、介观动力学(MesoDyn)和耗散粒子动力学（DPD）、统计方法QSAR(Quantitative Structure - Activity Relationship )等多种先进算法和X射线衍射分析等仪器分析方法；同时产品提供了界面友好的的模拟环境，研究者能方便地建立三维结构模型，并对各种小分子、纳米团簇、晶体、非晶体以及高分子材料的性质及相关过程进行深入的研究，得到切实可靠的数据。
Materials Studio分子模拟软件支持Windows和Linux操作平台，用户可以自由定制、购买自己的模拟方法和模块，以满足特定领域研究需求。
fe团簇的穆斯堡尔谱
穆斯堡尔谱是一种研究原子核与电子之间相互作用的谱学技术。
对于fe团簇的穆斯堡尔谱，可以提供关于其原子核的信息，例如原子核的电荷状态、磁性以及与周围电子的相互作用等。
穆斯堡尔谱的实验通常通过测量样品辐射出的γ射线频率和能量分布来获得。
对于fe团簇的穆斯堡尔谱，可以观察到不同的峰，每个峰对应于不同的穆斯堡尔参数。
穆斯堡尔参数可以提供关于样品中铁原子的信息。
其中，最重要的参数是穆斯堡尔位移，它表示了γ射线相对于无磁场时的频率偏移。
穆斯堡尔位移可以反映铁原子的电子云密度以及与周围电子的相互作用。
正的穆斯堡尔位移通常表示电子云密度较大，而负的穆斯堡尔位移则表示电子云密度较小。
另一个重要的参数是穆斯堡尔宽度，它表示了穆斯堡尔峰的宽度。
穆斯堡尔宽度可以提供关于铁原子的磁性信息。
对于磁性样品，穆斯堡尔宽度通常较大，而对于非磁性样品，穆斯堡尔宽度较小。
此外，穆斯堡尔谱还可以提供关于铁原子与周围晶格的相互作用以及样品中可能存在的其他相的信息。
需要注意的是，fe团簇的穆斯堡尔谱可能会受到多种因素的影响，例如样品的制备方法、温度、压力等。
因此，在
进行穆斯堡尔谱实验和解析时，需要综合考虑这些因素，以获得准确的结果。
Materials Studio的模块Materials Studio是一个全尺度材料模拟平台。
平台以可视化视窗界面Materials visualizer为核心，在其上共整合了24个功能模块，囊括了量子力学方法、经典模拟方法、介观模拟方法、有限元模拟等各种常见分子模拟方法，以及晶体结构解析、晶体形貌预测、定量构效关系分析等实用工具，实现了从电子结构解析到宏观性能预测的跨尺度研究。
Materials visualizerMaterials visualizer是Materials Studio的图形化界面，也是整个平台的核心。
Materials visualizer的功能包括：●搭建、调整各类三维可视的结构模型，包括晶体、小分子、聚合物、纳米材料、团簇、表界面以及各种缺陷结构；●提供模块参数设置、结果分析的视窗界面；提供结构文件、参数文件以及结果文件的管理界面；提供计算进程的监控界面；●对模拟结果进行各种分析，可与结构模型相结合进行数据的二维、三维显示，可以给出数据的图表，可以对特定的结果进行动画演示或给出矢量图；Materials visualizer的特性包括：●支持多种结构、图形、文本文件格式的输入和输出；●支持不同功能模块间结构数据的共享；●提供Perl语言环境，以及脚本编译工具；●提供不规则多面体表面积、体积的计算工具。
量子力学方法量子力学方法（Quantum Mechanics）是一种能够对材料体系电子结构特点进行解析的方法，精度高且几乎不依赖于任何经验参数，因此被广泛应用在各类材料的模拟研究中。
半经验量子力学方法（Semi-empirical Quantum Mechanics）同样能够对材料体系电子结构特点进行解析，但是包含有更多的经验参数以及数学、物理近似，因此，计算效率相比于纯粹的量子力学更高，但是精度会略低。
量子力学以及半经验量子力学方法均以定态薛定谔方程为核心，计算原子核满足特定排列、堆积时，核外电子的空间、能量分布，并由此进一步得到体系的电学性质、磁学性质、光学性质、热力学性质、力学性质，所能研究的材料体系类型包括：各类晶体材料及可能的各种缺陷结构，各种维度的纳米材料，各种分子及团簇材料。
量子力学方法：是一种能够对材料体系电子结构特点进行解析的方法，精度高且几乎不依赖于任何经验参数，因此被广泛应用在各类材料的模拟研究中。
CASTEP（平面波赝势方法）Cambridge Sequential Total Energy Package)是由剑桥凝聚态理论研究组开发的一款基于密度泛函理论的先进量子力学程序。
程序采用平面波函数描述价电子状态，利用赝势替代内层电子，因此也被称为平面波赝势方法。
研究对象主要有半导体、陶瓷、金属、分子筛等各类晶体材料，以及掺杂、位错、界面、表面等各种缺陷结构。
CASTEP的主要功能：能量计算--吸附热，缺陷形成能，内聚能，表面能等结构优化--优化原子坐标和晶胞参数，支持原子分数坐标、晶胞参数、键长、键角、二面角限定，支持外加应力过渡态--过渡态搜索（Synchronous Transit方法）电子结构解析--能带结构，电子态密度（局域态密度、分波态密度），电荷密度差分电荷密度，电子局域函数，电子轨道，扫描隧道显微镜(STM)图像模拟，共价键级，静电荷（Mulliken、Hirshfeld），静电势，功函数，自旋磁矩，费米面介电性质 --波恩有效电荷，静态介电常数张量，极化率张量力学性质计算 --弹性力常数张量，体模量，剪切模量，杨氏模量，泊松比热力学性质计算 --声子态密度、色散谱（linear response的方法扩展至金属体系6.0）；熵，焓，自由能，零点能，德拜温度，等容热容随温度的变化曲线光学性质计算 --红外光谱，拉曼光谱5.0（计算指定频率范围的拉曼活性模强度6.1），核磁共振谱a（化学位移、电场梯度张量），电子能量损失谱4.4（旋轨耦合效应5.5），X射线吸收（发射）谱4.4（旋轨耦合效应5.5），光频介电常数虚（实）部，吸收系数，折射率，光导率虚（实）部，能量损失函数动力学计算（结果分析采用Forcite plus的分析工具，具体内容参考Forcite plus介绍）---支持NVE，NVT，NPT以及NPH等系综，以及多种控温控压函数D Mol3（原子轨道线性组合方法）；DMol3是由Bernard Delley教授发布的一款基于密度泛函理论的先进量子力学程序，它采用原子轨道线性组合的方法描述体系的电子状态，因此也被称为原子轨道线性组合方法。
前言团簇是指几个乃至上千个原子、分子或离子通过一定的物理或化学结合力组成的相对稳定的微观或亚微观聚集体，团簇是介于原子、分子与宏观固体物质之间的物质结构的新层次，是各种物质由原子分子向大块物质转变的过渡状态。
它的很多新的物理、化学现象日益引起人们的兴趣，对团簇的研究是新材料设计以及研究凝聚态物质性质和特征的一项重要课题[1-9]，在大块凝聚态物质的特征、化学吸附以及固体表面反应、燃烧机理、晶体生长、催化等领域的研究中被广泛应用[2-16]。
团簇科学起源于19世纪中叶人们对烟雾、云以及碰撞等现象的研究和后来对成核现象的研究。
目前团簇科学研究的几个主要方向是：（1）研究团簇的组成及电子构型的规律、幻数和几何结构、稳定性的规律；（2）研究团簇的成核和形成过程及机制，研究团簇的制备方法、尤其是获取尺寸均一与可控的团簇束流；（3）研究金属、半导体及非金属和各种化合物团簇的光、电、磁、力学、化学等性质，它们与结构和尺寸的关系，及向大块物质转变的关节点；（4）研究团簇材料的合成和性质；（5）探索新的理论，不仅能解释现有团簇的效应和现象，而且能解释和预知团簇的结构，模拟团簇动力学性质，指导实验；（6）发展新的方法对团簇表面进行修饰和控制。
团簇的研究方法有很多，既有基础实验方法又有理论计算方法。
早期实验方面利用蒸发和热解的方法形成团簇，利用溅射或喷雾的方法可以从表面或凝聚相中获得团簇离子，离子的缔合和生长形成团簇，超声分子束的方法产生团簇。
现在主要通过直接激光气化法来形成团簇；利用团簇的丰度分布来推测团簇结构，偏离丰度分布的结构为“幻数”结构；实验上还采用离子色谱、团簇的光解离或碰撞解离以及团簇与气态小分子的吸附反应等方法来探究团簇的结构以及增长过程。
理论计算方法，基于经验势的，基于半经验方法的，基于量子计算方法的(不带任何经验参数)。
经验势方法是用一组含有参数的解析表达式来描述原子间的相互作用。
近年来提出的势能函数有：嵌入原子模型(EAM)势，Gupta势等。
第44卷第3期2021年3月V ol.44，No.3March2021核技术NUCLEAR TECHNIQUES相互作用势对原子核α团簇结构的影响狄宣平1,2方德清1,3何万兵3马余刚1,31（中国科学院上海应用物理研究所上海201800）2（中国科学院大学北京100049）3（复旦大学现代物理研究所教育部核物理与离子束应用重点实验室上海200433）摘要本文使用扩展的量子分子动力学模型（Extension of quantum molecular dynamics，EQMD），通过调整泡利势参数，使模型计算出的原子核平均结合能和均方根半径更符合实验值和经验公式。
在EQMD模型中加入动量相关势，根据无限大核物质的状态方程推导得到不同的核相互作用势参数，研究不同核势深度对4≤A≤40的α共轭核（N=Z，偶偶核）内的α团簇结构形成的影响。
结果表明：核势越深越容易在α共轭核中产生α团簇结构；同时，核势深度的变化也会影响12C和16O原子核内α团簇的构型。
关键词核相互作用势，α团簇，EQMD模型中图分类号O571.2，O571.25DOI:10.11889/j.0253-3219.2021.hjs.44.030501Effect of nuclear potential on formation ofαcluster in nucleiDI Xuanping1,2FANG Deqing1,3HE Wanbing3MA Yugang1,31(Shanghai Institute of Applied Physics,Chinese Academy of Sciences,Shanghai201800,China)2(University of Chinese Academy of Sciences,Beijing100049,China)3(Key Laboratory of Nuclear Physics and Ion-beam Application(MOE),Institute of Modern Physics,Fudan University,Shanghai200433,China)Abstract[Background]The cluster structure in light nuclei has been a hot topic in the nuclear structure research.Cluster phenomenon is related to the effective nuclear interaction,but the detailed mechanism of clusters in nuclei is still unclear.[Purpose]This study aims to investigate the effect of nuclear potential on formation ofαcluster in nuclei.[Methods]By adjusting the Pauli potential parameters in the extension of quantum molecular dynamics model(EQMD),the average binding energy and root mean square radius of the nucleus calculated by the model are in better agreement with the experimental values and empirical formulas.Momentum-dependent nuclear interactions (MDI)is added to the EQMD model.According to the state equation of infinite nuclear matter,different nuclear interaction potential parameters are derived.[Results]The deeper the Skyrme potential,the easier theαcluster structure will be produced.[Conclusions]The nuclear interaction potential has a great influence on the formation of αcluster structure in nuclei.Key words Skyrme potential,αcluster,EQMD model国家自然科学基金(No.11925502,No.11935001,No.11961141003,No.11927901,No.11705031)、国家重点研发计划(No.2018YFA0404404)、上海市自然科学基金(No.19ZR1403100)和中国科学院战略性先导科技专项(B类)(No.XDB34030100)资助第一作者：狄宣平，女，1995年出生，2014年毕业于山东师范大学，现为硕士研究生，研究领域为原子核物理通信作者：方德清，E-mail：；何万兵，E-mail：收稿日期：2020-12-15，修回日期：2021-01-29Supported by National Natural Science Foundation of China(No.11925502,No.11935001,No.11961141003,No.11927901,No.11705031),the National Key R&D Program of China(No.2018YFA0404404),Shanghai Development Foundation for Science and Technology(No.19ZR1403100) and the Strategic Priority Research Program of Chinese Academy of Sciences(No.XDB34030100)First author:DI Xuanping,female,born in1995,graduated from Shandong Normal University in2014,master student,focusing on nuclear physics Corresponding author:FANG Deqing,E-mail:HE Wanbing,E-mail:Received date:2020-12-15,revised date:2021-01-29狄宣平等：相互作用势对原子核α团簇结构的影响结团现象在自然界中普遍存在，原子核内的团簇结构一直都是原子核物理的重要研究方向之一。
原子间相互作用势对δ'相在过渡区沉淀影响的计算机研究王永欣;陈铮;李晓玲;刘兵;马良;赵宇宏【期刊名称】《自然科学进展》【年(卷),期】2003(013)002【摘要】基于离散格点形式的微扩散方程(Langevin方程)和非平衡自由能函数,编制了引入原子间相互作用能变化的Al3Li(δ')相沉淀原子层面计算机模拟程序.该程序包容亚稳区到失稳区的全部温度、成分范围,孕育期至粗化的全过程,可以处理与时间相关的过程问题.开展了不同原子间相互作用势下原子图像、序参数的计算机模拟,进而探讨了最近邻原子相互作用能(W1)对Al3Li相沉淀的影响机制.探明过渡区合金,随W1的增大,A13Li相沉淀的孕育期缩短,形核率增加,合金有序化速度和原子簇聚速度加快,在所研究的时间范围内达到的长程序参数和成分偏离序参数最大值增大.W1的增大促进Al3Li相以失稳分解形式析出.【总页数】7页(P179-185)【作者】王永欣;陈铮;李晓玲;刘兵;马良;赵宇宏【作者单位】西北工业大学材料科学与工程学院,西安,710072;西北工业大学材料科学与工程学院,西安,710072;西北工业大学材料科学与工程学院,西安,710072;西北工业大学材料科学与工程学院,西安,710072;西北工业大学材料科学与工程学院,西安,710072;西北工业大学材料科学与工程学院,西安,710072【正文语种】中文【中图分类】TG1【相关文献】1.相场法研究化学组成对Ni-Al-Ti合金γ′相早期沉淀过程的影响 [J], 孙远洋;赵宇宏;侯华;郑晓娟;郭慧俊2.原子间相互作用势对中Al浓度Ni75AlxV25-x合金沉淀序列的影响 [J],3.三元体系γ'相和θ相沉淀早期计算机研究 [J], 赵宇宏;陈铮;王永欣;唐丽英;卢艳丽4.应变能对Ni_(75)Cr_(16.4)Al_(8.6)合金沉淀过程影响的微观相场研究 [J], 霍进良;王永欣;陈铮;张静;赵彦;张利鹏5.共格沉淀相粗化行为的计算机模拟研究 [J], 赵宇宏;陈铮;王永欣;刘兵;马良因版权原因，仅展示原文概要，查看原文内容请购买。
Fe团簇的自旋-轨道耦合效应曾加莉西南大学物理科学与技术学院，重庆400715摘要：最近实验研究发现，在小尺寸铁团簇中，互补的轨道磁矩被强烈的壳弱化。
（Niemeyer等人，PRL，108，057201（2012年））通过考虑强大的电子相关效应，自旋 - 轨道耦合，以及非线性原子间的磁化，我们给出了自旋和轨道磁矩的理论解释以及中性和阳离子铁团簇（n=2-20）电子特性的理论解释。
使用了GGA + U这种方法，同时强调其对磁矩的作用。
我们发现，在不考虑库仑相互作用U时，自旋（轨道）力矩的平均值在2.69到3.50µB/atom（0.04到0.08µB/atom）之间；考虑库仑相互作用U时，磁场值在2.75到3.80µB/atom（0.10和0.30µB/atom）之间，这与实验测量结果非常吻合。
我们的研究结果证实在小尺寸铁团簇中，自旋磁矩较少被壳弱化，而轨道磁矩却被强烈地壳弱化。
通过GGA和GGA+ U函数优化Fe团簇后，总是得到共线的的磁基态。
当团簇尺寸发生改变时，几何演化表现出，为了形成大尺寸Fe团簇，二十面体形态将会与六角棱柱形态进行竞争。
此外，计算所得的电离势，电子亲和势，分裂能和极化趋势普遍与相应实验观察结果相同。
关键词：铁团簇；自旋磁矩；自旋轨道耦合；The spin and orbital moment of Fe n (n=2-20) clustersZENG JialiSchool of Physical Science and Technology, Southwest China University, Chongqing, 400715.China Abstract:Complementary to the recent experimental finding that the orbital magnetic moment is strongly quenched in small Fe clusters (Niemeyer et al., Phys. Rev. Lett. 108, 057201 (2012), we provide the theoretical understanding of the spin and orbital moments as well as the electronic properties of neutral and cation Fe n clusters (n=2-20) by taking into account the effects of strong electronic correlation, spin-orbit coupling, and noncollinearity of inter-atomic magnetization. The GGA+U method is used and its effluence on the magnetic moment is emphasized. We find that without inclusion of the Coulomb interaction U, the spin (orbital) moments have an average value between 2.69 and 3.50 µB/atom (0.04 and 0.08 µB/atom). With inclusion of U, the magnetic value is between 2.75 and 3.80 µB/atom (0.10 and 0.30 µB/atom), which provide an excellent agreement with the experimental measurements. Our results confirm that the spin moments are less quenched while the orbital moments are strongly quenched in small Fe clusters. Both GGA andGGA+Ufunctionals always yield collinear magnetic ground-state solutions for the fully relaxed Fe structures. Geometrical evolution, as a function of cluster size, illustrates that the icosahedral morphology competes with the hexagonal-antiprism morphology for large Fe clusters. In addition, the calculated trends of ionization potentials, electron affinities, fragment energies, and polarizabilities generally agree with respective experimental observations.Key word:Fe clusters；spin moments；spin-orbit coupling；1I.引言最近，科学家们对研究磁性过渡金属（TM）团簇，表现出了相当大的兴趣，因为这些团簇，有可能呈现出显著的磁各向异性能，这对于传统的信息存储是至关重要的。
α⁃Fe 和γ⁃Fe 长程F ⁃S 势的分子动力学模拟程江伟张先明吴永全∗王秀丽郑少波蒋国昌(上海大学材料科学与工程学院,上海市钢铁冶金重点实验室,上海200072)摘要：根据Sutton 等推导面心立方金属长程F ⁃S 势函数的方法推导得到α⁃Fe 和γ⁃Fe 的最优参数分别为,ε=0.2453,a =0.28664nm,n =7,m =4,c =7.7525和ε=0.0006,a =0.36467nm,n =15,m =4,c =1104.7351.用所推导的势参数对常压下不同温度时α⁃Fe 和γ⁃Fe 的性质进行了分子动力学(MD)模拟,结果表明,计算得到的Fe 的微观结构(径向分布函数,配位数,原子的配位状态)和宏观物性(线性膨胀,密度)都能与实验结果相吻合,说明所推导的长程F ⁃S 势函数参数适用于α⁃Fe 和γ⁃Fe 的MD 模拟.关键词：长程F ⁃S 势;α⁃Fe;γ⁃Fe;微观结构中图分类号：O641;TF01MD Simulation of α⁃Fe and γ⁃Fe with Long ⁃range F ⁃S PotentialCHENG Jiang ⁃WeiZHANG Xian ⁃Ming WU Yong ⁃Quan ∗WANG Xiu ⁃LiZHENG Shao ⁃BoJIANG Guo ⁃Chang(Shanghai Enhanced Laboratory of Ferrometallurgy,School of Material Science and Engineering,Shanghai University,Shanghai 200072,P.R.China )Abstract ：Several sets of long ⁃range F ⁃S potential parameters for α⁃Fe and γ⁃Fe were deduced according to the method introduced by Sutton et al .,and the optima were determined as follows:ε=0.2453,a =0.28664nm,n =7,m =4,c =7.7525for α⁃Fe and ε=0.0006,a =0.36467nm,n =15,m =4,c =1104.7351for γ⁃Fe.Accordingly,the isothermal ⁃isobaric MD simulations were carried out with the set of optimal parameters for α⁃Fe and γ⁃Fe at different temperatures in barometric condition.The most ⁃scale agreements between simulations and experiments both for microstructures and macroscopic properties strongly validate the application of this set of parameters to the MD simulation of α⁃Fe and γ⁃Fe,and further to the study of interfacial dynamics between metal and oxide.Key Word ：Long ⁃range F ⁃S potential;α⁃Fe;γ⁃Fe;Microscopic structure[Note]物理化学学报(Wuli Huaxue Xuebao )Acta Phys.鄄Chim.Sin .,2007,23(5)：779-785May Received:August 30,2006;Revised:January 5,2007;Published on Web:April 21,2007.∗Corresponding author.Email:yqwu@ or yqwu@;Tel:+8621⁃56332144.国家自然科学基金重点基金(钢铁联合基金)(50334050),国家自然科学基金(50504010),上海市自然科学基金(04ZR14054)资助项目ⒸEditorial office of Acta Physico ⁃Chimica Sinica在早期对金属体系的研究中,对势因其简单性而获得了普遍的应用,但它并不能体现金属体系常见而又重要的多体电子效应.为此,Daw 和Baskes [1]于1983年提出了基于量子力学理论的半经验势函数方法研究金属中的氢脆问题,发展出原子嵌入势(embedded atom method,EAM)[2],同年Finnis 等[3]提出过渡金属的简单经验多体势,简称为F ⁃S 势.这两种将多体势引入到金属性质理论计算的方法代表了原子间作用力经验描述的显著进步及当时用经验或半经验方法描述金属凝聚的计算机模拟的最高水平.基于紧束缚态密度的二阶动量近似,这两种势考虑了局域电子密度对原子间作用力的影响,相对于对势的优点就在于能从物理意义上更为准确地描述金属的特性.在此基础上,长程F ⁃S 势[4]综合了F ⁃S 势和L ⁃J 势的优点(F ⁃S 势和EAM 势的作用范围仅扩展到原子的第三近邻,长程F ⁃S 势具有较大的作用范围,这为我们在后续的工作中模拟界面特性提供了方便,也是我们选择长程F ⁃S 的重要原因),以779Acta Phys.鄄Chim.Sin.,2007Vol.23期实现原子团簇间的机械性能的模拟,取得了很好的结果[4-7].Fe及其合金在异相时的物理特性,及Fe与氧化物之间的界面结构及行为,一直是材料学界关注的焦点之一,除了实验手段的革新之外,也发展了多种有效的理论计算方法,其中包括经验或半经验的分子力学及分子动力学方法.当前,对于具体考察金属Fe与氧化物(比如硅酸盐、铝酸盐等)界面结构及行为,尤其是液态Fe与纳米氧化物颗粒之间的界面行为的研究还非常少见,而这又是当前冶金界的研究难点之一.在具体考察Fe与氧化物界面之前,首先必须分别对金属Fe和氧化物的体材料微观结构和宏观属性有完整的认识和描述.据此,吴永全等人[8-12]以分子动力学及Raman光谱为手段对氧化物体材料的结构及性质进行了研究并得到了很好的结果.到目前为止,对于包括Fe在内的过渡金属的分子动力学模拟,其成果多集中在采用F⁃S势对同一晶相内点缺陷、断裂、表面驰豫、空位形成能、内聚性、相变、表面问题、磁性能和合金特性等的研究上[1-5,13-29],而针对Fe的体材料微观结构(尤其是γ相)及其在相变过程中的变化的研究报道还相对较少[30-34].本文工作期望获得适用于体材料Fe的描述的长程F⁃S势函数参数,从而对不同相时金属Fe的微观结构和宏观物性进行模拟研究.作者根据长程F⁃S势函数参数的推导原理推导出一套用于描述α⁃Fe 和γ⁃Fe的势函数参数,并将之应用于纯铁在升温相变过程中的微观结构变化及相应的体积相关量变化的分子动力学描述中.1长程F鄄S势函数参数确定方法在长程F⁃S势中总能量由体现局域原子和电子密度对总能量贡献的多体势和对势组成.基于紧束缚理论中二阶动量近似的观点,Finnis和Sinclair认为多体相互作用的函数具有电子密度平方根的形式,通常F⁃S势的作用范围仅扩展到第三紧邻而没有1/r6这一范德华尾项,但该势在物理意义上能更为准确地描述金属的性质;而L⁃J势虽然在描述金属原子间的相互作用上存在不足,但是却能很好地描述长程相互作用.Sutton等[4]综合考虑了以上两点,把长程F⁃S势函数总能量写成如下形式:E r=ε［12∑i≠j∑V(r ij)-c j∑ρ1/2i](1)V(r ij)=(a/r ij)n(2)ρi=i≠j∑(a/r ij)m(3) V(r ij)代表排斥项,多体项代表引力项,r ij是原子i和j 的间距.c是正的无量纲参数,ε是具有能量量纲的参数,a是具有长度量纲的参数,m和n是正整数.在长程F⁃S势中需要确定ε、a、n、m和c这五个参数.在长程F⁃S势中,Sutton等仅给出了面心立方金属的势函数参数的推导方法,我们据此推导了面心立方Fe的F⁃S势参数,同时借鉴推导面心立方Fe (γ⁃Fe)势参数的方法外延推导了体心立方Fe(α⁃Fe)的F⁃S势函数参数.正如Sutton等推导面心立方金属时考虑了晶体的结构因素一样,在推导体心立方金属的势参数时我们同样考虑了α⁃Fe的结构因素(包括原胞构型和晶格常数,这一点是通过晶格和S f n来体现的).推导方程见公式(4)-(7).S f n=j∑(a f/r j)n S f m=j∑(a f/r j)m(4) c=nS f n/[m(S f m)1/2](5)ε=2mE f c/[2n-m)S f n](6) 18Ωf B f/E f c=nm(7)γ⁃Fe的势函数参数由方程(4)-(7)[4]确定.对γ⁃Fe而言,方程(4)中的a f为面心立方Fe的晶胞常数, r j为晶胞中两个原子的间距;S f n、S f m是Sutton等在长程F⁃S势中定义的晶格和,n、m为待定参数.可以看出它与晶体的结构相关,是反映晶体构型的参数.方程(5)和(6)中c和ε即长程F⁃S势中的参数.方程(7)中Ωf为面心立方Fe的原子体积,B f为面心立方Fe 的体变模量.由方程(4)可以求出对应于不同的n与m时晶胞的晶格和S f n与S f m,由方程(5)可以求出对应于不同的n 和m的c值,然后由方程(6)求出对应于不同n和m的ε值.由此可以看出,n和m是两个独立变量,在方程(7)的约束条件下,它们可以任意取整数,然后结合三个外来实验参数———晶格常数a f,内聚能E f c和体变模量B f,从而形成多套势函数参数.然后对这些势函数参数分别进行对比验证,获得的最优的一套势函数参数就是我们需要的势函数参数.根据上述原则,在确定势函数参数过程中采用的外来实验参数列于表1.最终得到α⁃Fe和γ⁃Fe的最优势函数参数列于表2.2模拟结果与势函数参数的验证对由686个Fe原子组成的立方体系进行了模780No.5程江伟等：α⁃Fe 和γ⁃Fe 长程F ⁃S 势的分子动力学模拟拟,以验证我们推导得到的势函数参数的可行性.采用NPT 系综和Berendsen 方法进行压力和温度控制,立方晶胞的初始边长为2.0065nm,截断半径为1.0000nm,每次模拟时间40ps,时间步长为1fs.为使体系充分达到平衡,趋衡时间为35ps.对所有模拟体系都做了体系能量和温度随时间的变化图,如图1(体系的目标温度为1667K)所示,35ps 的模拟时间足以保证体系达到平衡.据此,后续的分析都是模拟时间在35ps 以后进行的.2.10K 时势参数结构稳定性的检验采用所推导的两套势参数对0K 时的Fe 进行模拟.采用α⁃Fe 和γ⁃Fe 的势参数得到体系的能量曲线如图2所示,可以看出,常压0K 时以α⁃Fe 的势参数对Fe 进行模拟时得到的体系能量低于以γ⁃Fe 的势参数模拟时体系能量.说明常压0K 时体心立方结构的α⁃Fe 比面心立方结构的γ⁃Fe 的结构稳定,这与低温时Fe 以α⁃Fe 的形式存在相一致.图3和图4分别为0K 与293K 时α⁃Fe 径向分布函数和配位数函数图.图3(a)与图4(a)相比,0K 时Fe 原子的径向分布函数的峰尤为尖锐,晶格的周期性也更为明显.在图3(b)中可以看到0K 时Fe 的第一和第二配位数分别为8.0和6.0.采用α⁃Fe 和γ⁃Fe 的势参数对0K 时Fe 的模拟结果说明了所选取的势参数的结构稳定性,在检验了势参数的结构稳定性之后,采用这些势参数对Fe 的物性进行了模拟.2.2α⁃Fe 和γ⁃Fe 的微观结构图5是1185K 时γ⁃Fe 的径向分布函数和配位数函数图,表3列出了293K 时α⁃Fe 和1185K 时γ⁃Fe 的径向分布函数的峰位和配位数情况.图4(a)中293K 时α⁃Fe 的径向分布函数很好地反映了该温度下α⁃Fe 的长程有序性,图4(b)标示了293K 时α⁃Fe 的第一配位数及第一、第二配位数之和.表3反映了该温度下α⁃Fe 各峰位和第一配位数以及第一、第二配位数之和的实验值与计算值的比较(峰位的实验值由该温度下α⁃Fe 晶格常数的实验值[35]推出,配位数实验值由体心立方晶型推出),所列举的前四峰位都能与实验结果相吻合,这在过渡金属α⁃Fe 结构计算的相关文献的报道中是很少见a /nmE c /eV B /GPa γ⁃Fe 0.36467[35] 4.26[36]132[37]α⁃Fe0.28664[35]4.28[38]173[3]表1确定势函数参数时采用的实验参数Table 1Experimental data for determiningpotential parametersε/eVa /nm n m c α⁃Fe 0.24530.28664747.7525γ⁃Fe0.00060.364671541104.7351表2α⁃Fe 和γ⁃Fe 的势函数参数Table 2Potential parameters of α⁃Fe and γ⁃Fe图11667K 时体系能量(a)和温度(b)随时间的演化Fig.1The system energy (a)and temperature (b)evolutions with time at 1667K图20K 时分别由α⁃Fe 和γ⁃Fe 的势参数得到的体系的能量曲线Fig.2The system energies obtained by the potential parameters of α⁃Fe and γ⁃Fe at 0K781Acta Phys.鄄Chim.Sin.,2007Vol.23到的.事实上,Shimomura 等人[23]发现,用EAM 方法计算体心立方金属时的结果不如计算面心立方金属时的结果与实验值吻合得好.而本文采用F ⁃S 势模拟的体心立方Fe 的结果却很理想.此外,计算得到293K 时α⁃Fe 的第一配位数为8.69,与理论上体心立方Fe 的配位数8.0相差0.69,同时第一、第二配位数之和为14.0,与体心立方Fe 的第一、第二配位数之和为14.0相吻合.第一配位数与理论值产生的偏差很可能是以下原因引起的:293K 时α⁃Fe 的第一近邻原子和第二近邻原子之间的间距很小(仅为0.0384nm [35]),使得第一配位原子和第二配位原子不易区分开来,并且随着温度的升高就会更难分辨(相应地在图30K 时α⁃Fe 的径向分布函数g (r )(a)和配位数函数n (r )(b)Fig.3The RDF g (r )(a)and coordination numberfunction n (r )(b)of α⁃Fe at 0K 图4293K 时α⁃Fe 的径向分布函数g (r )(a)和配位数函数n (r )(b)Fig.4The RDF g (r )(a)and coordination numberfunction n (r )(b)of α⁃Fe at 293K图51185K 时γ⁃Fe 的径向分布函数g (r )(a)和配位数函数n (r )(b)Fig.5The RDF g (r )(a)and coordination numberfunction n (r )(b)of γ⁃Fe at 1185KI p /nm1st CN 1st and 2nd CNs 1st2nd3rd4thα⁃Fe calculated 0.25330.28840.41260.48828.6914.0experimental [35]0.24820.28660.40540.49658.014.0γ⁃Fecalculated0.25210.36120.44180.498712.1experimental [39]0.25780.36470.44670.515612.0表3293K 时α⁃Fe 和1185K 时γ⁃Fe 的径向分布函数的峰位(I p )和配位数Table 3The peaks (I p )of RDFs and coordinate number (CN)of α⁃Fe at 293K and γ⁃Fe at 1185K782No.5程江伟等：α⁃Fe 和γ⁃Fe 长程F ⁃S 势的分子动力学模拟径向分布函数的演化上表现为一个宽化的强峰,这一点将在α⁃Fe 的原子局部配位状态的分析中给出明确的说明);此外,体系的大小、时间步长、所采用的系综和控温控压方法都会对模拟结果产生影响,比如体系越大越能反映真实的情况,而体系太小则会带来较大的体积效应,模拟结果与真实情况相差较大;模拟时的时间步长越短,也越能接近真实情况,但是要花费较大的机时.而时间步长过大,则会使模拟结果偏离真实值.在模拟过程中我们采用的是NPT 系综,立方晶胞的初始边长为2.0065nm,采用Berendsen 方法进行压力和温度控制.每次模拟时间40ps,时间步长为1fs.客观地说,模拟的原胞体积有偏小之嫌,而时间步长为1fs,则有偏大之嫌.这些都可能是导致模拟结果中配位数接近9.0的原因.表3中反映了1185K 时γ⁃Fe 各峰位实验值[39](各峰的实验值由1185K 时γ⁃Fe 的晶格常数推出)与计算值的比较.从表3可以看到计算所得各峰位值与实验值基本吻合.与α⁃Fe 相比,γ⁃Fe 的各峰尖锐程度下降,除了第一峰外其余的峰都较弱而且较宽.随着温度升高,有强峰变弱、弱峰消失的趋势,这是温度升高引起原子运动加剧的结果.2.3α⁃Fe 和γ⁃Fe 的热膨胀为验证势参数对宏观物性模拟的适用性,对α⁃Fe 的线性膨胀和γ⁃Fe 的密度随温度变化的关系进行考察,其中参比值L 0为293K 时α⁃Fe 晶胞的轴向长度,L 为各温度点下α⁃Fe 晶胞的轴向长度.图6(a)中比较了α⁃Fe 的线性膨胀的计算值与实验值[40],1185K 时计算值与实验值相差达到最大,相对误差为13.8%,误差较大;图6(b)比较了γ⁃Fe 密度随温度变化时的计算值与实验值[40,41],可以看出,计算值与两组实验值都能很好地吻合,最大绝对误差约为0.05g ·cm -3.2.4α⁃Fe 和γ⁃Fe 结构的变化与温度的关系图7是不同温度时Fe 的径向分布函数,从图7可以看出,293K 时Fe 具有很好的长程有序性,各峰都清晰可见;随着温度的升高各个峰不断地展宽,造成两个或者多个峰的融合.当温度为1185K 时第一峰和第二峰已完全演化成一个峰,此时径向分布函数发生了显著的变化,Fe 从α相转变为γ相.在α相内,由于第一配位原子和第二配位原子间距很小,随着温度升高更难以分辨,当温度从293K 升高到相变点1185K 时,在径向分布函数上第一峰和第二峰融合在一起,相应地,此时的第一配位数不单包图6不同温度时α⁃Fe 的线性膨胀(a)和γ⁃Fe 的密度(b)Fig.6The linear expansions of α⁃Fe(a)and densities of γ⁃Fe(b)at differenttemperatures图7Fe 在不同温度时的径向分布函数g (r )(a)和配位数CN(b)Fig.7The RDF(a)and coordinate number N (b)of Fe at differenttemperatures783Acta Phys.鄄Chim.Sin.,2007Vol.23括Fe的第一配位原子(理论值为8.0),也包括了Fe的第二配位原子(理论值为6.0),其数值也从8.7增加到14.0,这一点在α⁃Fe原子的局部配位状态中也将给出进一步的说明;在γ相内,随着温度的升高,径向分布函数的各峰也不断弱化,但径向分布函数没有出现像α相内那么明显的峰的消失与融合,其配位数变化也没有较大变化,随着温度从1185K升高至1667K,γ⁃Fe的配位数从12.1增加至12.5.图7 (b)显示了温度从293K升高到1185K的过程中α⁃Fe和γ⁃Fe配位数的变化情况,在1185K时晶体结构的转变导致在配位数上也有一个急剧的变化. 2.5α⁃Fe和γ⁃Fe原子的局部结构图8(a)-(g)中所有构型均来自于模拟得到的平衡构像,对应各图中心原子的配位数分别为8、9、10、11、12、13和14,反映了温度从293K升高至1185K时α⁃Fe原子局部配位状态的变化;图8中(h)为γ⁃Fe原子典型的局部配位状态,其中心原子的配位数为12.从图8(a)-(g)可以看出,随着温度的升高,原子热运动的增强,体心立方构型的α⁃Fe中心原子先与邻近的中心原子键合,从而其配位数不断增大.当中心原子的配位数分别为9、10、11、12和13时,由于配位原子间键长、键角等的不对称性而产生的相互挤压和排斥,导致这种状态的配位单元结构稳定性较差.随着配位数增加至14,配位原子与中心原子间的键长再次趋于均匀化,此时的配位单元结构处于相对稳定的状态.可以直观地看出,当中心原子配位数为14时,配位原子的确包括了第一近邻的8个顶点原子和邻近的6个中心原子,这与之前对配位数变化的分析相吻合;同时,面心立方构型的γ⁃Fe的原子局部配位状态也与其理论值12相吻合,并可以直观地看出中心原子与配位原子的空间位置关系.需要说明的是,为便于清楚地显示出中心原子的配位状态而又不影响考察中心原子的局部结构,适当地去除了部分非中心原子之间的连键.最后需要指出的是,在势参数的选取过程中,虽然我们参考了Sutton等的做法,但是不同的模拟需要使得我们在势参数的选择上也有所不同.对Sutton 等来说,在确定长程F⁃S势参数时,能否准确地模拟金属原子团簇的机械作用是势参数优劣的首要标准.我们则侧重于不同晶相时Fe的微观结构及相应的宏观体积物性的模拟,为使进行计算机模拟时势能作用范围不是太大,Sutton等建议m≥6,但在我们选取的所有参数中,m≥6时并不能得到合理的结果,正如Sutton等在文献[4]中提到的那样,他们对不同面心立方金属进行模拟时也发现,相对较小的m时模拟结果与实验值更吻合.3结论本文按照Sutton等推导面心立方金属长程F⁃S 势函数的方法,获得了γ⁃Fe的势参数,并借鉴该方法推导出体心立方α⁃Fe的长程F⁃S势函数参数,得到一套α⁃Fe和γ⁃Fe的最优参数;对该套势函数参数的模拟结果给出了详细的分析和验证,并得到如下结论:按照本文采用的方法推导出的长程F⁃S势函数参数在对α⁃Fe和γ⁃Fe的模拟上有很好的适用性,在微观结构(径向分布函数、配位数函数及原子局部配位状态)和体积相关量(线性膨胀和密度)等性质上都能与实验值相吻合;从长程F⁃S势对α⁃Fe的模拟结果来看,长程F⁃S势不单适用于对面心立方铁的模拟,也适用于对体心立方铁的模拟.本文所获图8不同配位状态时中心α⁃Fe原子的局部结构Fig.8The local structure of central atoms with various coordinates inα⁃Fe(a)-(f):coordination number ofα⁃Fe atoms vary from8to14;(h):typical coordination status of14atom inγ⁃Fe784No.5程江伟等：α⁃Fe和γ⁃Fe长程F⁃S势的分子动力学模拟得的F⁃S势函数的最优参数,使得我们可以进一步从微观和宏观两个层面对Fe进行更深入的研究,并将之应用于金属与氧化物界面的研究中.References1Daw,M.S.;Baskes,M.I.Phys.Rev.Lett.,1983,50(17):12852Daw,M.S.;Baskes,M.I.Phys.Rev.B,1984,29(12):64433(a)Finnis,M.W.;Sinclair,J.E.Philos.Mag.A,1984,50(1):45;(b)Finnis,M.W.;Sinclair,J.E.Philos.Mag.A,1986,53:161(Erratum)4Sutton,A.P.;Chen,J.Philos.Mag.Lett.,1990,61(3):1395Ackland,G.J.;Vitek,V.Phys.Rev.B,1990,41(15):103246Pethica,J.B.;Sutton,A.P.J.Vac.Sci.Technol.A,1988,59:321 7Smith,J.R.;Bozzolo,G.;Banerjea,A.;Ferrante,J.Phys.Rev.Lett.,1989,63:12698Wu,Y.Q.;Jiang,G.C.;You,J.L.;Hou,H.Y.;Chen,H.J.Chem.Phys.,2004,121(16):78839Wu,Y.Q.;Jiang,G.C.;You,J.L.;Hou,H.Y.;Chen,H.Acta Phys.Sin.,2005,54(2):961[吴永全,蒋国昌,尤静林,侯怀宇,陈辉.物理学报,2005,54(2):961]10Wu,Y.Q.;Jiang,G.C.;You,J.L.;Hou,H.Y.;Chen,H.Chin.Phys.Lett.,2002,19(12):188011Wu,Y.Q.;You,J.L.;Jiang,G.C.;Chen,H.Chin.J.Inorg.Chem., 2004,20(2):133[吴永全,尤静林,蒋国昌,陈辉.无机化学学报,2004,20(2):133]12You,J.L.;Jiang,G.C.;Hou,H.Y.;Chen,H.;Wu,Y.Q.;Xu,K.D.J.Raman Spectrosc.,2005,36:23713Liu,R.S.;Qi,D.W.;Wang,S.Phys.Rev.B,1992,45:45114Li,H.;Bian,X.F.;Zhang,J.X.Matter.Sci.Eng.A,1999,271: 11615Zhou,G.R.;Wu,Y.S.;Zhang,C.J.;Zhao,F.Acta Phys.⁃Chim.Sin.,2003,19(1):13[周国荣,吴佑实,张川江,赵芳.物理化学学报,2003,19(1):13]16Zhang,T.;Gu,T.K.;Qi,Y.H.Acta Phys.⁃Chim.Sin.,2005,21(2):173[张弢,谷廷坤,齐元华.物理化学学报,2005,21(2):173]17Liu,X.;Meng,C.G.;Liu,C.H.Acta Phys.⁃Chim.Sin.,2004,20(3):280[刘新,孟长功,刘长厚.物理化学学报,2004,20(3):280]18Ackland,G.J.;Thetford,R.Philos.Mag.A,1987,56:1519Daw,M.S.;Foiles,S.M.Phys.Rev.B,1986,35:212820Daw,M.S.;Foiles,S.M.Phys.Rev.Lett.,1987,59:275621Ackland,G.J.;Tichy,G.;Vitek,V.;Finnis,M.W.Philos.Mag.A, 1987,56:73522Li,H.;Wang,G.H.;Ding,F.;Wang,J.L.;Shen,W.F.Phys.Lett.A.,2001,280:32523Shimomura,Y.;Sugio,K.;Kogure,Y.;Doyama,put.Mater.Sci.,1999,14:3624Li,X.H.;H,J.F.J.Solid State Chem.,2003,176:23425Foiles,S.M.Surf.Sci.,1987,191:32926Ackland,G.J.;Bacon,D.J.;Calder,A.F.;Harry,T.Philos.Mag.A,1997,75:71327Yang,Z.;Johnson,R.A.Modeling Simul.Mater.Sci.Eng.,1993, 1:70728Phythian,W.J.;Foreman,A.J.E.;Stoller,R.E.;Bacon,D.J.;Calder,A.F.J.Nucl.Mater.,1995,223:24529Besson,R.;Morillo,J.Phys.Rev.B,1997,55:19330Colakoglu,K.;Ugur,G.;Cakmak,M.;Tutuncu,H.M.Turkish Journal of Physics,1999,23:47931Bacon,D.J.;Diaz de la Rubia,T.J.Nucl.Mater.,1994,216:275 32Bacon,D.J.;Calder,A.F.;Kapinos,V.G.;Wooding,S.J.Nucl.Instrum.Meth.B,1995,102:3733Belonoshko,A.B.;Ahuja,R.;Johansson,B.Phys.Rev.Lett., 2000,84:363834Koci,L.;Belonoshko,A.B.;Ahuja,R.Phys.Rev.B,2006,73: 22411335Massalski,T.B.;Okamoto,H.;Subramanian,P.R.;Kacprzak,L.Binary alloy phase diagrams.2nd ed.Metals Park,Ohio:ASMInternational,1990(1):T1036Rosato,V.Acta.Metall.,1989,37:275937Zarestky,J.;Stassis,C.Phys.Rev.B,1987,35:4500.38Kittle,C.Introduction to solid state physics.5th.New York: Wiley,1976:7439Stern,M.J.Electrochem.Soc.,1955,102(12):66340ASM,Metals handbook.Volume2.9th ed.Trans.Ma,J.R.;Chen,L.M.Beijing:Machinery Industry Press,1994:946[金属手册.第2卷.第9版.马九荣,陈立敏译.北京:机械工业出版社,1994:946]41Basinski,Z.S.;Hume⁃Rothery,W.;Sutton,A.L.Proceedings of the Royal Society of London,Series A:Mathematical and Physical Sciences,1995,229(1179):459785。
