平衡晶格常数及体积模量
硅的晶格常数和体弹模量的计算
Project #1硅的晶格常数和体弹模量的计算一、平衡晶格常数和内聚能自然条件下硅为金刚石结构(dc )。
计算模拟时，我们可以假定它为各种结构，f cc, bcc, sc, dc. 可以预测，模拟的dc 结构的硅的体系能量最低，也即最稳定。
下面我们将运用LAMMPS 来对硅的各种结构进行模拟。
定义晶格能量为Φ， 数密度为 ρ：potE N Φ=N Vρ= 其中E pot 为势能，N 为体系总原子数，V 为体系的体积。
选取 Stillinger-Weber (SW)，以下面命令执行 lammps 运算：其中，lmp_serial 为 lammps 命令；”<” 符号为读取符；in.Silicon 为输入文件，里面包含运算所需要的各种数据和命令；-log 指定输出文件的名称。
可以看到屏幕上显示出lammps 运行的信息。
这个计算量很小，所以很快就结束。
接下来以如下命令来查看计算得到的数据：grep 是linux 中一个很重要的命令，用来搜索文本，读取匹配的行并打印出来。
这里是搜索 dc.log 文件，将 @ 开头的行打印出来。
如下：晶格参数为5.4305埃,数密度为0.0499540303，每个原子的能量为-4.336599609eV.下面具体来看刚才给的输入文件，in.Silicon .dc.log 文件中有原子总数的信息，每个金刚石晶胞中有8个原子，383216⨯=，所以是216个原子。
如下给出各种结构下的体系的原子数：晶体结构类型晶胞中的原子数 总原子数 简单立方SC1 27 体心立方BCC2 54 面心立方FCC4 108 金刚石DC 8 216表1.不同晶体结构中的原子数下图是计算模拟得出的各种结构下的数密度与每个原子能量的关系图。
横坐标为数密度， 以金刚石为例，ρ= 8/5.4315^3=0.049926，也即我们直接通过 grep 命令得到的第二项值；纵坐标为每个原子的能量，为第三项值。
TiN多型体高压相变的第一性原理计算-物理学报
ⓒ2011中国物理学会Chinese Physical Society http ://TiN 多型体高压相变的第一性原理计算∗顾 雄 高尚鹏†(复旦大学材料科学系,上海 200433)(2010年8月30日收到;2010年9月3日收到修改稿) 基于密度泛函理论框架下的赝势平面波方法,计算了B 1(氯化钠结构)、B 2(氯化铯结构)、B 3(闪锌矿结构)、B k (六方氮化硼结构)、B h (碳化钨结构)和B 81(砷化镍结构)6种TiN 多型体的晶体结构、体积弹性模量和相对稳定性.计算指出,不存在B 4(纤锌矿)结构的TiN.通过不同外压下的晶格弛豫得到每种结构的焓,发现外压低于345GPa 时TiN 最稳定的结构是B 1结构,在345GPa 以上B 2结构的TiN 最稳定.TiN 从B 1结构到B 2结构的相变伴随着体积和弹性模量的变化.关键词:氮化钛,赝势,高压相变,密度泛函理论PACS :71.15.Mb,64.70.K -,62.20.de,61.50.Ah∗国家自然科学基金(批准号:10804018)和教育部留学回国人员科研启动基金资助的课题.†通讯联系人.E⁃mail:gaosp@1.引言TiN 是一种新型多功能材料,具有高强度、高硬度、耐高温、耐酸碱侵蚀、耐磨损等一系列优点,在机械工业中被用做材料表面的耐磨、耐腐蚀、耐火材料涂层和切削工具、力学部件上的硬质防磨损覆层[1—4].此外TiN 薄膜还在微电子业被用来做扩散势垒,在汽车和玻璃制造业充当反射材料,在珠宝业用作表面的黄金色涂[3,4].TiN 的热力学稳定相为氯化纳结构,由于TiN 被广泛应用做耐磨材料和切削工具,瞬时可能承受很大的应力,所以对TiN 晶体的各种可能存在的结构在高压环境下的相对稳定性和性能的研究具有实际的意义.基于密度泛函理论的第一性原理计算已广泛用于材料高压相变和高压下的材料性能研究[5—17].在计算中可以很方便地设定不同的加压条件,如实验室中金刚石压腔可以实现的从MPa 到102GPa 范围的各种压强,或者大星体内部以及极端化学条件下的TPa(103GPa)量级的极高压强[5].TiN 的结构、电子结构、力学及表面等性质的密度泛函理论计算已有较多报道[18—20],但是高压下TiN 结构相变的第一性原理计算研究还很少,我们只检索到最近Wang 等人根据总能⁃体积曲线推断出TiN 由NaCl 结构到纤锌矿结构转变的压强为-21.0GPa,即不能通过加正的外压发生[21].TiN 的高压相变需要进一步的系统研究,来澄清在高压下是否存在比NaCl 结构更稳定的相,以及各亚稳相在不同压强下的相对稳定关系.本文采用赝势平面波方法,研究多种TiN 异构体的结构和性能以及在不同压强下的相对稳定性.主要研究的多型体包括B 1(NaCl 结构)、B 2(CsCl 结构)、B 3(闪锌矿结构)、B k (六方氮化硼结构,h⁃BN)、B h (WC 结构)以及B 81(NiAs 结构).为了表述简洁,下文中很多地方采用Strukturbericht 符号(即B 1,B 2等)来指代各种结构.2.计算方法 计算采用密度泛函理论框架下的赝势平面波方法,交换关联能采用PBE 形式的GGA 近似[22].布里渊区采用Monkhorst⁃Pack 撒点[23],最大近邻K 点间距为0.043Å-1.选择BFGS (Broyden Fletcher Goldfarb Shanno)方法[24]进行几何优化来确定零压和不同外压条件下的平衡结构和性质.本工作中所加外压方式均为流体静压强.在给定外压下,各相的相对稳定性可由焓H =E +PV 来确定,焓最小的结构最稳定.所有计算工作均采用CASTEP 程序[25]完成.在赝势平面波方法用于高压研究时,赝势构造过程中所采用的芯区截断半径r c的值是一个需要特别注意的参量.由于所采用的压强可能很高,作为理论预测研究也不必受限于目前技术手段可能达到的压强范围,在高压下晶体的晶格常数变小,近邻原子距离变小,而赝势方法的基本概念则要求最相邻原子如Ti和N的r c之和要等于或小于原子间距,也就是近邻原子的芯区不能重叠.程序赝势库所提供的赝势的r c取值一般是保证适用于常规压强下的,所以对高压下的适用性要进行测试.在所计算的各种结构的TiN中,B3TiN在同等压强下近邻Ti原子和N原子的间距最小.在200GPa时计算得到的近邻的Ti原子与N原子的距离r0为1.727Å,小于CASTEP程序赝势库里的Ti原子与N原子超软赝势[26]的芯区半径之和,r c(Ti)+ r c(N)=1.799Å.这就意味着在此压强条件下采用此芯区半径的赝势不再适用,必须修改r c(Ti)与r c (N),使研究压强下的不同结构TiN的最近邻原子距离r0都大于或等于r c(Ti)+r c(N).将r c(Ti)修改为0.900Å,r c(N)修改为0.423Å,即r c(Ti)+r c (N)=1.323Å,采用新的芯区半径产生一组超软赝势.表1的测试结果表明直到1900GPa,六种结构的TiN均满足r0>r c(Ti)+r c(N).取更小的芯区截断半径会使赝势的转移性更好,但是也会使赝势更“硬”,也就是平面波截断能更高,经收敛条件测试后,我们取平面波截断能为550eV.理论上讲,取两种不同的芯区半径所构造的赝势在其适用范围内所得到的结果应一致,图1(a)给出分别采用两组赝势B1TiN和B3TiN在100GPa范围内的各压强下 的总能差值和焓差值,两组赝势给出的结果差别在0.014eV以内,图1(b)给出分别采用两组赝势计算的晶格常数,差别在0.02Å以内,计算结果基本一致和预期相符.在本论文以下计算中,我们将统一采取芯区截断半径较小的一组超软赝势.图1 B1(氯化纳结构)和B3(闪锌矿结构)TiN在100GPa外压范围内的总能差及焓差(a)和晶格常数(b).叉号代表采用CASTEP程序超软赝势库里的赝势(Ti和N的芯区截断半径分别为0.952Å和0.847Å)的计算结果,三角符号代表采用新的较小的芯区截断半径(Ti和N分别为0.900Å和0.423Å)产生的赝势的计算结果表1 6种结构的TiN在0GPa和1900GPa下近邻原子的距离(单位Å)B1B2B3k h1 1900GPa1.5351.6361.3391.5281.5611.6133.计算结果和讨论3.1.TiN异构体的结构和性能我们采用7种常见的晶体结构,即B1(NaCl结构)、B2(CsCl结构)、B3(闪锌矿结构)、B4(纤锌矿结构)、B k(h⁃BN结构)、B h(WC结构)以及B81(NiAs结构),作为初始模型来进行晶格弛豫计算得到不同结构TiN的平衡结构参数和体积弹性模量.计算中我们发现以纤锌矿结构作为起始结构进行几何优化后,内部参数u(对于理想纤锌矿结构为3/ 8)变为1/2.也就是说Ti原子和N原子将在同一平面上,每一层的Ti,N原子相间排列成六角环状网络,相邻的原子层有沿c轴方向的Ti—N键,这是类似于六方氮化硼的结构(即B k).其他6种结构的TiN的结构参数和弹性模量在表2中给出.实验报道的B1TiN晶格常数为4.238Å[27], Stampfl等使用FPLAPW⁃GGA计算的数值为4.26Å[18],Marlo等采用赝势平面波方法和PBE近似计算的数值为4.238Å[19],我们的理论计算得出晶格常数为4.248Å,与实验和其他理论计算结果的差别符合一般认可的GGA近似计算精度.其他结构的TiN作为可能存在的亚稳相还没有相应的实验结构参数报道,我们的理论结果可供将来的实验和计算工作参考.B2TiN的密度最大,而B3TiN的密度最小.具有高硬度是TiN获得广泛应用的一个重要原因,体积弹性模量是与硬度相关的一个重要参量,零压下B1TiN的体积弹性模量是最大的,而B3TiN的体积弹性模量是最小的.由图2给出B1和B2TiN的体积弹性模量随压强的变化关系可以看出,B2TiN的体积弹性模量在低压下略小于B1 TiN,随着压强的增大与B1TiN的差别变得更小,在约993GPa以上大于B1TiN.表2 6种结构TiN的结构参数和体弹模量结构空间群压强/GPa晶格常数/Åa(b)c密度/g·cm3体弹模量/GPaB1(NaCl)Fm⁃3m04.2484.2485.301271 (225)504.0724.0726.0874653003.6793.6798.254128410003.2973.29711.4673198 B2(CsCl)Pm⁃3m02.6482.6485.533247 (221)502.5212.5216.4144401002.4422.4427.0546153002.2682.2688.812124710002.0282.02812.3203200 B3(闪锌矿)F⁃43m04.6134.6134.186191 (216)504.3474.3475.003373 B k(h⁃BN)P63/mmc03.5364.2094.510213(194)503.3334.0175.3174011003.2143.8975.893569 B h(WC)P⁃6m202.9142.7085.158247(187)502.7582.6065.9884341002.6592.5436.59760210002.0982.27211.8682940 B81(NiAs)P63/mmc02.9825.2135.118254 (194)502.8275.0185.9164471002.7304.8986.50062310001.9685.09312.03326053.2.TiN各异构体外压下的相对稳定性在有外界给定压强的条件下,焓最低的结构最稳定.在压强变化时,异构体的相对稳定性可能会有变化.图3给出了6种结构的TiN在不同压强下的焓,我们计算的压强范围从0到1900GPa,为了使图线能更清楚的表达,各压强下的焓均取B1结构的焓作参考值.图4给出了6种结构TiN的总能和体积随压强的变化关系,从而可以分析内能和体积变化分别对焓的相对大小改变的贡献.从图3中可以看出,在无外压和压强较小的情况下,B1TiN焓最小所以是最稳定的结构,这和TiN 最常见的相是B1一致.在无外压时,焓值最大的是B2TiN,但是随着压强的增大,约36GPa以后,B2图2 B 1与B 2结构的TiN体积弹性模量随压强的变化图3 6种结构的TiN 的焓随压强的变化,以B 1结构的焓作参考值.为便于比较,图示焓的数据均对应于一个TiN 单元(两个原子)TiN 的焓将比B 3结构的小,压强大于72GPa 时,B 2TiN 的焓将比B k 结构的小,压强大于162GPa 时,B 2TiN 焓先将比B h 结构的小,接着在165GPa 将比B 81结构的小.压强增加到345GPa 以后,B 2结构的TiN 焓最小,所以B 2TiN 结构最稳定.在345GPa 左右,TiN 的B 1结构和B 2结构具有不同的体积(见图4(b)),含有两个原子的原胞体积分别为12.00Å3和11.23Å3,所以TiN 从B 1结构到B 2结构的相变伴随着体积的变化.由于压强P =-d E /d V 为能量⁃体积(E⁃V 曲线)的负斜率,在某个压强P 下两种结构的焓相等即P =-(E 2-E 1)/(V 2-V 1),所以材料在高压下的结构相变除了通过观察图3所示的焓⁃压强曲线图4 6种结构的TiN 的总能(a)和体积(b)随压强的变化.为便于比较,图示的能量和体积数据均对应于一个TiN 单元(两个原子)外,利用两种结构的E⁃V 曲线也可以判断转变压强.图5是B 1TiN 和B 2TiN 的E⁃V 曲线,共切线斜率的负值即转变压强约为2.184eV /Å3,即约350GPa.图5 B 1(NaCl 结构)TiN 和B 2(CsCl 结构)TiN 的能量⁃体积曲线.共切线(虚线)给出转变压强有实验报道在0—45GPa 的压强范围内,TiN样品的体积随着压强的增加曲线变化比较平滑,TiN 没有发生相变[28,29].我们的计算结果表明在0到345GPa范围内,B1结构都是最稳定的,与该实验结果一致.只有在345GPa以上压强,B2TiN才是最稳定的.由于从B1结构到B2结构的相变可能存在的势垒,在实验上要观察到此相变应该在345GPa 以上某个更高的压强处.CaO,SrO和CdO在高压下也会发生B1结构到B2结构的相变,实验上观测到的转变压强分别为60(±2)GPa(CaO)[30],36(±4)GPa(SrO)[31]和90.6GPa(CdO)[32],均小于TiN 的B1到B2结构的转变压强.量子蒙特卡罗计算预测MgO从B1结构到B2结构的转变压强要接近600GPa[33].FeO随着压强的增大则从B1结构转变为一种菱形结构,在高温高压条件下(70GPa以上, 1600K)可发生B1结构到B81结构的转变[34].B3和B k结构的TiN在较低压强下比B2TiN 焓更小,而与B h结构以及B81结构较为接近,随着压强的增大,在36GPa以后B3TiN焓变为最大,在72GPa以后B k TiN的焓仅次于B3TiN而大于其他4种结构的TiN.从图4(b)可以看出,B3结构和B k 结构的TiN的体积一直是最大的两种,内能随压强的迅速增加是导致它们焓超过B2TiN的主要原因.由图3可以明显看到,在较高压强下B3TiN和B k TiN的焓明显高于其他四种相,因此属于很难在高压下存在的相.B h与B81结构的TiN的焓在研究压强范围内十分接近,在外压较小时,略大于B1TiN.随着压强增大,在大约1010GPa时,B h TiN的焓将比B1结构的小,在大约1260GPa时,B81TiN的焓将比B1结构的小.6种结构的相对稳定性排序由无外压时的B1>B k>B81>B3>B h>B2变为高压下(1260 GPa以上)的B2>B h≈B81>B1>B k>B3. 4.结论利用密度泛函理论框架下的赝势平面波方法,研究了TiN的B1(NaCl结构)、B2(CsCl结构)、B3(闪锌矿结构)、B k(h-BN结构)、B h(WC 结构)和B81(NiAs结构)6种异构体的平衡结构和在不同压强下的相对稳定性.为了保证在极高压下近邻原子的距离大于或等于原子赝势芯区半径之和,在计算中我们使用了特别构造的具有较小芯区半径的超软赝势.计算结果给出了六种异构体相对稳定性随压强的变化,在高压下B3结构和B k结构的TiN都是焓明显较其他四种结构高也就是说相对较不稳定,而B81结构和B h结构的焓在低压下比B1结构的焓略大,随着压强的增大分别在1010GPa和1260GPa的极高压强下焓会小于B1结构,但会大于高压下的B2结构,即不存在焓最小的压强范围.B1TiN直至345GPa都具有最小的焓,也就是结构最稳定,在345GPa以上, B2TiN焓最小.从B1结构到B2结构的相变是体积不连续相变.B2TiN的弹性模量在低压下略小于同等压强下的B1TiN,而在约993GPa以上略大于B1TiN.我们的计算结果表明TiN在承受较大应力的切削等应用中,当等效压强小于345GPa 的范围内,应不需担心结构相变问题.本论文对于各种TiN亚稳相在不同压强下的结构和性能的研究结果为将来的高压理论和实验工作提供有益的参考.[1]Schwarz K1987Crit.Rev.Solid State Mater.Sci.13211[2]Rickerby D S,Burnett P J1987Surf.Coat.Technol.33191[3]Valvoda V1996Surf.Coat.Technol.8061[4]Yan P X,Wu Z G,Xu J W,Zhang Y J,Li X,Zhang W W2004J.Synth.Cryst.33974(in Chinese)[闫鹏勋、吴志国、徐建伟、张玉娟、李 鑫、张伟伟2004人工晶体学报33974][5]Pickard C J,Needs R J2010Nat.Mater.9624[6]Silas P,Yates J R,Haynes P D2008Phys.Rev.B78174101[7]Yu R,Zhan Q,De Jonghe L C2007Angew Chem.Int.Edit461136[8]Alptekin S,Durandurdu M2009Solid State Commun.149345[9]Ma X G,Liang P,Miao L,Bie S W,Zhang C K,Xu L,JiangJ J2009Phys.Status Solidi B2462132[10]Guan P F,Wang C Y,Yu T2008Chin.Phys.B173040[11]Zhu J,Yu J X,Wang Y J,Chen X R,Jing F Q2008Chin.Phys.B172216[12]Feng H J,Liu F M2009Chin.Phys.B181574[13]Tan L N,Hu C E,Yu B R,Chen X R2007Chin.Phys.163772[14]Chen D,Chen J D,Zhao L H,Wang C L,Yu B H,Shi D H2009Chin.Phys.B180738[15]Ji Z H,Zeng X H,Hu Y J,Tan M Q2008Acta Phys.Sin.573753(in Chinese)[季正华、曾祥华、胡永金、谭明秋2008物理学报573753][16]Ji G F,Zhang Y L,Cui H L,Li X F,Zhao F,Meng C M,SongZ F2009Acta Phys.Sin.584103(in Chinese)[姬广福、张艳丽、崔红玲、李晓凤、赵 峰、孟川民、宋振飞2009物理学报584103][17]LüM Y,Chen Z W,Li L X,Liu R P,Wang W K2006ActaPhys.Sin.553576(in Chinese)[吕梦雅、陈洲文、李立新、刘日平、王文魁2006物理学报553576][18]Stampfl C,Mannstadt W,Asahi R,Freeman A J2001Phys.Rev.B63155106[19]Marlo M,Milman V2000Phys.Rev.B622899[20]Liu L M,Wang S Q,Ye H Q2005J.Phys.:Condens.Matter.175335[21]Wang A J,Shang S L,Du Y,Kong Y,Zhang L J,Chen L,Zhao D D,Liu Z K2010Comput.Mater.Sci.48705 [22]Perdew J P,Burke K,Ernzerhof M1996Phys.Rev.Lett.773865[23]Monkhorst H J,Pack J D1976Phys.Rev.B135188[24]Pfrommer B G,Cote M,Louie S G,Cohen M L1997J.Comput.Phys.131133[25]Clark S J,Segall M D,Pickard C J,Hasnip P J,Probert M J,Refson K,Payne M C2005Z.Kristallogr.220567[26]Vanderbilt D1990Phys.Rev.B417892[27]Schöenberg N1954Acta Chem.Scand.8213[28]Guo Z T,Peng F,Chen H H2010J.Southwest Univ.National.(Nat.Sci.Edit.)36145[郭振堂、彭 放、陈海花2010西南民族大学学报(自然科学版)36145][29]Chen H H,Peng F,Mao H K,Shen G Y,Liermann H P,Li Z,Shu J F2010J.Appl.Phys.107113503[30]Jeanloz R,Ahrens T,Mao H K,Bell P M1979Science206829[31]Sato Y,Jeanloz R1981J.Geophys.Res.86B11773[32]Liu H,Mao H,Maddury M,Ding Y,Meng Y,Häusermann D2004Phys.Rev.B70094114[33]AlfèD,Alfredsson M,Brodholt J,Gillan M J,Towler M D,Needs R J2005Phys.Rev.B72014114[34]Murakami M,Hirose K,Ono S,Tsuchiya T,Isshiki M,Watanuki T2004Phys.Earth.Planet.Inter.146273Ab initio calculation of pressure⁃induced phasetransition of TiN polytypes∗Gu Xiong Gao Shang⁃Peng†(Department of Materials Science,Fudan University,Shanghai 200433,China)(Received30August2010;revised manuscript received3September2010)AbstractBased on a plane wave pseudopotential method within the framework of density functional theory,equilibrium structure,bulk modulus,and relative stability were calculated for6kinds of TiN polytypes including B1(NaCl structure),B2(CsCl structure),B3(zincblende structure),B k(hexagonal B N structure),B h(WC structure)and B81 (NiAs structure).Theoretical calculation also showed that TiN can not exist in B4(wurtizite)structure.Through geometry optimization under hydrostatic pressure,the enthalpy of each TiN phase at different pressures was obtained.It was found that TiN with B1structure is the most stable phase at pressure lower than about345GPa,whereas B2TiN is the most stable at pressure above345GPa.Volume discontinuity and bulk modulus change can be observed during the transition from B1to B2phase.Keywords:TiN,pseudopotential,high pressure phase transition,density functional theoryPACS:71.15.Mb,64.70.K-,62.20.de,61.50.Ah∗Project supported by the National Natural Science Foundation of China(Grant No.10804018)and the Scientific Research Foundation for the Returned Overseas Chinese Scholars,State Education Ministry.†Corresponding author.E⁃mail:gaosp@。
ABINIT实战手册
ABINIT-5.8.4软件的Windows版本实战守则（by H.J.Zhao）前言第一性原理计算发表文章的要点1. 你算的是什么，2. 你为什么要算这个东西，3. 你用了怎样的算法，4. 你算出来了什么，5. 你怎样解释你算出来的东西，有什么结论。
第一章ABINIT软件的运行/user1/11542/archives/2005/363031.shtml以BaTiO3为例，需建立输入文件*.in，在abinit网站下载Ba、Ti、O的赝势文件，之后执行abinis程序，依次输入：1. bto.in2. bto.out3. btoi4. btoo5. bto.temp6. 把三个赝势文件拖入程序执行框中运行即可！Abinit使用篇简介abinis的输入文件分为三类：定义了所有输入文件名的文件（比如in.files，这个文件就是用来告诉abinit哪些文件是abinit读入参数的文件名，以及输出主要结果的文件的名称），定义了计算的控制参数的文件（比如取名为INP)，赝势的文件。
下面举例如何做好输入文件计算金刚石结构的Si的状态方程，由此也得到Si的晶格常数和体弹性模量。
大致的步骤是在Si的晶格常数的实验值附近取11个数据点，也就是说取11个晶格常数或体积，然后计算在这些晶格常数下的总能。
在计算得到总能后，采用状态方程拟合得到状态方程、平衡态时的体积（或晶格常数）和体弹性模量。
本例子中采用的是LDA-HGH赝势。
赝势的文件名为：14si.4.hgh。
in.files的内容为（紫色标示）：#####BeginINP #设置关键词的文件名为INPOUT #主要的输出文件为OUT，该文件将被写入计算最重要的结果siisisi14si.4.hgh #赝势的文件名######ENDINP文件的内容为：# Crystalline cubic Si#ndtset 11 #说明下面将有11组数据acell: 3*9.8112 # 晶格常数a=b=c，将从9.8112.... a.u.开始增加acell+ 3*0.09 #晶格常数将以0.09 a.u.的间隔进行增加#Ground state calculationkptopt 1 #在k点网格取样时根据对称性来取样，并由下面的# ngkpt和kptrlatt, 或者nshiftk和shiftk来确定k点的数目iscf 5 #采用CG方法对能量进行优化，用在基态计算中。
2-1(黄昆-固体物理)-教案
第二章固体的结合本章的主要内容：阐明原子是依靠怎样的相互作用而结合成为固体的。
教学重点：晶体结合的基本类型。
教学难点：晶体结合的物理本质。
教学时数：3学时。
讲授方式：PPT文档。
§ 2.1 离子性结合1. 教学目的和要求: 通过讲解使学生理解并掌握离子性结合成为晶体的本质。
2．教学重点：离子晶体结合的特点及性质。
3．教学难点：离子晶体结合的基本物理量：内能、晶格常数、体变模量和结合能的计算。
4．讲授时间：45分钟。
5．讲授方式：PPT文档。
6．作业：2.1，2.2，2.3。
一．离子晶体概述元素周期表中第I族碱金属元素（Li、Na、K、Rb、Cs）与第VII族的卤素元素（F、Cl、Br、I）化合物（如NaCl，CsCl，晶体结构如图XCH001_009_01和XCH001_010所示）所组成的晶体是典型的离子晶体，半导体材料如CdS、ZnS等亦可以看成是离子晶体。
二．离子晶体结合的特点以CsCl 为例，在凝聚成固体时，Cs 原子失去价电子，Cl 获得了电子，形成离子键。
以离子为结合单元，正负离子的电子分布高度局域在离子实的附近，形成稳定的球对称性的电子壳层结构；Xe ICs Kr BrRb Ar ClK Ne FNa ⇒⇒⇒⇒-+-+-+-+,,,（1）离子晶体的模型：可以把正、负离子作为一个刚球来处理；离子晶体的结合力：正、负离子之间靠库仑吸引力作用而相互靠近，当靠近到一定程度时，由于泡利不相容原理，两个离子的闭合壳层的电子云的交迭会产生强大的排斥力。
当排斥力和吸引力相互平衡时，形成稳定的离子晶体；（2）一种离子的最近邻离子为异性离子；（3）离子晶体的配位数最多只能是8（例如CsCl 晶体）；（4）由于离子晶体结合的稳定性导致了它的导电性能差、熔点高、硬度高和膨胀系数小； （5）大多数离子晶体对可见光是透明的，在远红外区有一特征吸收峰。
举例：氯化钠型(NaCl 、KCl 、AgBr 、PbS 、MgO)(配位数6)氯化铯型(CsCl 、 TlBr 、 TlI)（配位数8）离子结合成分较大的半导体材料ZnS 等（配位数4）三. 离子晶体结合的性质1）系统内能的计算晶体内能为所有离子之间的相互吸引库仑能和重叠排斥能之和。
分子动力学模拟用贵金属势函数的应用与发展
分子动力学模拟用贵金属势函数的应用与发展夏璐;陈永泰;王塞北;魏宽;李爱坤;慕阳;任县枬;陈松;陆建生;谢明;潘勇;胡洁琼;杨有才;张吉明;王松【摘要】势函数是用来描述介观层次系统中分子或原子间相互作用的一种数学模型，是在介观层次上进行材料显微结构分子动力学模拟的关键。
针对贵金属及其合金总结了分子动力学模拟研究中常用势函数的枑型、数学形式和基本参数，介绍了各种势函数在实际研究中的具体应用情况，分析讨论了贵金属势函数的构建和应用中急需解决的问题、以及势函数研究的发展趋势。
%Potential function is a kind of mathematical model used for describing the interactions between molecules or atoms in system at the mesoscopic scale. It’s crucial in studying the micro-structure of materials at the mesoscopic scale in molecular dynamic simulation. Some common interatomic potentials used in precious metals and their alloys are given in the present paper, including their specific functions, the parameters, the characteristics as well as practical applications. Finally, issues which need to work out were discussed in the process of building the potential models, and the tendency of future development of interatomic potentials of noble metals is put forward.【期刊名称】《贵金属》【年(卷),期】2013(000)004【总页数】9页(P82-90)【关键词】金属材料;贵金属;势函数;分子动力学;模拟;介观【作者】夏璐;陈永泰;王塞北;魏宽;李爱坤;慕阳;任县枬;陈松;陆建生;谢明;潘勇;胡洁琼;杨有才;张吉明;王松【作者单位】昆明理工大学材料科学与工程学院昆明 650093; 昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明理工大学材料科学与工程学院昆明 650093; 昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明理工大学材料科学与工程学院昆明 650093;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明 650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明650106;昆明贵金属研究所稀贵金属综合利用新技术国家重点实验室昆明650106【正文语种】中文【中图分类】TG146.3势函数是描述体系粒子总能量E与粒子坐标(r)的数学函数形式。
第二章 晶体的结合
固体材料是由大量的原子(或离子)组成约1 mol / cm 3原子(或离子)结合成晶体时,外层电子作重新分布,外层电子的不同分布产生了不同类型的结合力.Na Cl ＋＝NaCl 离子键共价键金属键结合力类型决定了晶体的微观结构乃至宏观物理性质.本章主要介绍不同结合类型中原子间相互作用与晶体内能、晶体的微观结构和宏观物理性质之间的联系.共价键结合（金刚石）－－原子间束缚非常强，导电性差金属键结合（金属Cu ）－－对电子束缚较弱，良导体——原子的电子分布核外电子的分布遵从泡利不相容原理、能量最低原理和洪特规则.能量最低原理电子在原子核外排布时，要尽可能使电子的能量最低1s、2s、2p、3s、3p、4s、3d、4p、4d……泡利不相容原理每一轨道中只能客纳自旋相反的两个电子.洪特规则电子在原子核外排布时，将尽可能分占不同的轨道，且自旋平行——原子的电子分布_----电离能_----电离能_----电离能_----电子亲和能_----原子电负性_----原子电负性2．Pauling鲍林提出的电负性计算方法（较通用）：_----原子电负性•横向•电离能•亲和能•电负性按结合力的性质和特点，晶体可分为5种类型：离子晶体（离子结合）共价晶体（共价结合）金属晶体（金属结合）氢键晶体（氢键结合）如何理解各种晶体呢？离子晶体:正离子与负离子的吸引力就是库仑力.共价结合:靠近的两个电负性大的原子各出一个电子形成电子共享的形式.金属结合:原子实依靠原子实与电子云间的库仑力紧紧地吸引着.氢键结合:氢先与电负性大的原子形成共价结合后, 氢核与负电中心不在重合, 迫使它通过库仑力再与另一个电负性大的原子结合.分子结合:电偶极矩把原本分离的原子结合成了晶体. 电偶极矩的作用力实际就是库仑力.可见, 所有晶体结合类型都与库仑力有关.原子间相互作用势能----结合力的共性吸引力排斥力库仑引力库仑斥力泡利原理引起(1)吸引力和排斥力都是原子间距离r 的函数.注:(2)吸引力是长程力,排斥力短程力.(3)当r =r 0时, 原子间合力为零, 原子处于平衡.类比于弹簧振子()()⎟⎠⎞⎜⎝⎛−−=−=++11n m r nB r mA dr r du r f 为什么排斥力是短程力？()()()B A r u r u r u +−=+=()()⎜⎛−−=−=nB mA r du r f设晶体中第i个原子与第j个原子之间的相互作用势能u(r)为ij()()∑∑∑==NNNr u r u U 1晶体的结合能：()()∑=N r u Nr u晶格常数由于晶格具有周期性，设临近两原子间距R，则晶体体积可写成体弹性模量单位压强引起的体积的相对变化率。
基于第一性原理的晶体材料力学性能计算方法
基于第一性原理的晶体材料力学性能计算方法孙金芳; 姚星星; 倪程鹏; 耿杰【期刊名称】《《西安工程大学学报》》【年(卷),期】2019(033)005【总页数】7页(P555-561)【关键词】第一性原理; 晶体材料; 力学性能; 弹性性能; 强度性能; 晶格缺陷【作者】孙金芳; 姚星星; 倪程鹏; 耿杰【作者单位】安徽信息工程学院通识教育与外国语学院安徽芜湖 241000【正文语种】中文【中图分类】TP393.410 引言晶体材料的主要性能为力学性能,目前获取晶体材料性能最好的方法是试错,该方法可以从周期表内找出多种化学元素组合,分析材料蕴含成分公式,效果比较显著,但花费时间很长,成本很高[1]。
所以寻找有效的力学性能计算方法十分重要。
伴随着计算机与建模方法的高速发展,不同材料的建模方法与电子结构的原理研究得到了很大进步,通过求解化学与晶体的自洽方程和产生的电子波函数,可以计算出晶体材料力学性能,计算过程大多数需要依赖电子,因此需要用原子核与固体来约束,从而实现电子结构预测[2]。
目前,相关学者已经对晶体材料力学性能计算方法进行了研究并取得了一些成果。
文献[3]提出了一种基于密度泛函理论的第一性原理平面波赝势法,在理想状态下,运用密度泛函理论分析空位因素对B2-NiAl力学性能的影响作用,实验结果表明该方法能够计算出晶体材料的力学性能,但是计算结果与实际值偏差较大；文献[4]运用第一性原理对B2型TiSi合金的缺陷结构进行了计算,包括弹性性质和晶格常数,根据计算结果可知,B2型TiSi合金结构的韧性不高,虽然得出了比较可靠的结论,但是该方法计算时间较长,实用性不高。
根据上述分析可知,虽然目前对于如何更好地计算出晶体材料的力学性能已经有了一些研究成果,但是已有方法仍然存在一些缺陷。
运用传统计算方法对晶体材料的力学性能进行计算时得到的结果并不理想,主要是因为计算时电子、原子、介质等各种因素的影响,除此之外,传统方法具有晶体结构敏感、晶体缺陷(晶体的空缺、间隙、错位、错层、结界)等问题。
BCC结构TiCrTaV合金的电子结构和力学性能的第一性原理计算
第４０卷 第１期原 子 与 分 子 物 理 学 报Ｖｏｌ．４０Ｎｏ．１２０２３年２月ＪＯＵＲＮＡＬＯＦＡＴＯＭＩＣＡＮＤＭＯＬＥＣＵＬＡＲＰＨＹＳＩＣＳＦｅｂ．２０２３!"#$"%&'"()*+",-.-/,0.1.2;..23<556ＢＣＣ结构ＴｉＣｒＴａＶ合金的电子结构和力学性能的第一性原理计算茹金豆１，２，马 瑞１，２，万明攀１，谢 泉２（１ 贵州大学材料与冶金学院，贵阳５５００２５；２ 贵州大学新型光电子材料与技术研究所，贵阳５５００２５）摘 要：利用基于密度泛函理论的第一性原理方法计算ＴｉＣｒＴａＶ多组元合金中两种ＢＣＣ结构的结构稳定性、力学性能、德拜温度、电子结构和布居分析 生成焓和内聚能结果表明ＢＣＣ１的结构稳定性更好，更容易形成 弹性常数和模量表明ＢＣＣ１的强度和韧性更强，ＢＣＣ２的抗剪切能力和刚度更好，两种结构均具有弹性各向异性 德拜温度和Ｇｒüｎｅｉｓｅｎ参数结果表明ＢＣＣ２的键合强度和热稳定性更好 电子结构和布居分析表明两种结构均包含共价键和金属键 Ｔａ原子形成的共价键强度更大，金属键仅存在于Ｔｉ、Ｃｒ和Ｖ原子之间 元素成键后Ｔｉ和Ｖ原子失去电子，Ｃｒ和Ｔａ原子得到电子关键词：多组元合金；第一性原理；力学性能；电子结构；布居分析中图分类号：Ｏ４８２ 文献标识码：Ａ ＤＯＩ：2."274889:"2... ./;0"-.-/".2;..2Ｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓｃａｌｃｕｌａｔｉｏｎｓｏｎｅｌｅｃｔｒｏｎｉｃｓｔｒｕｃｔｕｒｅｓａｎｄｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆＴｉＣｒＴａＶａｌｌｏｙｓｗｉｔｈＢＣＣｓｔｒｕｃｔｕｒｅｓ，２，ＭＡＲｕｉ１，２，ＷＡＮＭｉｎｇ Ｐａｎ１，ＸＩＥＱｕａｎ２ＲＵＪｉｎ Ｄｏｕ１（１ ＣｏｌｌｅｇｅｏｆＭａｔｅｒｉａｌｓａｎｄＭｅｔａｌｌｕｒｇｙ，ＧｕｉｚｈｏｕＵｎｉｖｅｒｓｉｔｙ，Ｇｕｉｙａｎｇ５５００２５，Ｃｈｉｎａ；２ ＩｎｓｔｉｔｕｔｅｏｆＡｄｖａｎｃｅｄＯｐｔｏｅｌｅｃｔｒｏｎｉｃＭａｔｅｒｉａｌｓａｎｄＴｅｃｈｎｏｌｏｇｙ，ＧｕｉｚｈｏｕＵｎｉｖｅｒｓｉｔｙ，Ｇｕｉｙａｎｇ５５００２５，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｔｈｅｓｔｒｕｃｔｕｒａｌｓｔａｂｉｌｉｔｉｅｓ，ｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓ，Ｄｅｂｙｅｔｅｍｐｅｒａｔｕｒｅｓ，ｅｌｅｃｔｒｏｎｉｃｓｔｒｕｃｔｕｒｅｓａｎｄｐｏｐｕ ｌａｔｉｏｎａｎａｌｙｓｉｓｏｆｔｈｅＢＣＣ１ａｎｄＢＢＣ２ｓｔｒｕｃｔｕｒｅｓｉｎＴｉＣｒＴａＶｍｕｌｔｉ－ｃｏｍｐｏｎｅｎｔａｌｌｏｙｓａｒｅｃａｌｃｕｌａｔｅｄｕｓｉｎｇｔｈｅｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓｍｅｔｈｏｄｂａｓｅｄｏｎｄｅｎｓｉｔｙｆｕｎｃｔｉｏｎａｌｔｈｅｏｒｙ ＴｈｅｆｏｒｍａｔｉｏｎｅｎｔｈａｌｐｙａｎｄｃｏｈｅｓｉｏｎｅｎｅｒｇｙｉｎｄｉｃａｔｅｔｈａｔＢＣＣ１ｓｔｒｕｃｔｕｒｅｈａｓｂｅｔｔｅｒｓｔｒｕｃｔｕｒａｌｓｔａｂｉｌｉｔｙａｎｄｓｔｒｏｎｇｅｒｆｏｒｍａｔｉｏｎａｂｉｌｉｔｙ Ｔｈｅｅｌａｓｔｉｃｃｏｎｓｔａｎｔｓａｎｄｍｏｄｕ ｌｕｓｉｎｄｉｃａｔｅｔｈａｔｔｈｅｓｔｒｅｎｇｔｈａｎｄｔｏｕｇｈｎｅｓｓｏｆＢＣＣ１ａｒｅｂｅｔｔｅｒ，ｗｈｉｌｅｔｈｅｓｈｅａｒｒｅｓｉｓｔａｎｃｅａｎｄｓｔｉｆｆｎｅｓｓｏｆＢＣＣ２ａｒｅｂｅｔｔｅｒ Ｂｏｔｈｓｔｒｕｃｔｕｒｅｓｈａｖｅｅｌａｓｔｉｃａｎｉｓｏｔｒｏｐｙ ＤｅｂｙｅｔｅｍｐｅｒａｔｕｒｅａｎｄＧｒüｎｅｉｓｅｎｐａｒａｍｅｔｅｒｓｓｈｏｗｔｈａｔＢＣＣ２ｈａｓｂｅｔｔｅｒｂｏｎｄｉｎｇｓｔｒｅｎｇｔｈａｎｄｔｈｅｒｍａｌｓｔａｂｉｌｉｔｙ Ｅｌｅｃｔｒｏｎｉｃｓｔｒｕｃｔｕｒｅａｎｄｐｏｐｕｌａｔｉｏｎａｎａｌｙｓｉｓｓｈｏｗｔｈａｔｂｏｔｈｐｈａ ｓｅｓｃｏｎｔａｉｎｃｏｖａｌｅｎｔａｎｄｍｅｔａｌｌｉｃｂｏｎｄｓ ＴｈｅｃｏｖａｌｅｎｔｂｏｎｄｓｆｏｒｍｅｄｂｙＴａａｔｏｍｓａｒｅｓｔｒｏｎｇｅｒ，ａｎｄｔｈｅｍｅｔａｌｌｉｃｂｏｎｄｓｅｘｉｓｔｏｎｌｙｂｅｔｗｅｅｎＴｉ，ＣｒａｎｄＶａｔｏｍｓ Ｗｈｅｎｔｈｅｅｌｅｍｅｎｔｓｂｏｎｄ，ＴｉａｎｄＶｌｏｓｅｅｌｅｃｔｒｏｎｓａｎｄＣｒａｎｄＴａｏｂｔａｉｎｅｌｅｃｔｒｏｎｓＫｅｙｗｏｒｄｓ：Ｍｕｌｔｉ－ｃｏｍｐｏｎｅｎｔａｌｌｏｙ；Ｆｉｒｓｔｐｒｉｎｃｉｐｌｅｓ；Ｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓ；Ｅｌｅｃｔｒｏｎｉｃｓｔｒｕｃｔｕｒｅｓ；Ｐｏｐｕｌａｔｉｏｎａｎａｌｙｓｉｓ收稿日期：２０２２ ０２ ０７基金项目：国家自然科学基金（１２０６５００５）作者简介：茹金豆（１９９５—），女，河北省石家庄人，在读硕士，主要研究方向：高熵合金第一性原理计算．通讯作者：马瑞．Ｅ ｍａｉｌ：ｒｍａ１＠ｇｚｕ．ｅｄｕ．ｃｎ第４０卷原 子 与 分 子 物 理 学 报 第１期１ 引 言核聚变作为一种新型能源引起世界范围内广泛关注 核聚变反应堆内部环境十分恶劣，要求材料具有非常优异的性能，能够承受极高的温度和应力梯度、辐照损伤、高浓度的嬗变产物和等离子体暴露 多组元合金（又称高熵合金）是一种由四种、五种或者更多元素以近似等原子比组成的新型合金 这些新型的合金具有高强度和硬度［１］、优异的延展性和损伤耐受性［２］、良好的耐磨性［３］、优异的热稳定性［４］、优异的抗氧化性［５］和耐腐蚀性［６］，因此具有很好的应用潜力由难熔金属元素组成的难熔多组元合金由于具有优异的高温综合力学性能和良好的抗辐照性能，因此人们将其作为先进核反应堆的候选材料．如ＡｌＮｂ０ ５ＴｉＺｒＭｏ０ ５Ｔａ０ ５难熔多组元合金室温屈服强度为２０００ＭＰａ，８００℃屈服强度为１５９７ＭＰａ，室温和高温的屈服强度均优于Ｎｉ基合金（Ｉｎｃｏｎｅ１７１８），并具有良好的塑性［７］ Ｓｅｎｋｏｖ等人研究发现ＮｂＴｉＺｒ基难熔多组元合金具有较高的塑性和屈服强度，如ＮｂＴｉＺｒＶＭｏ室温屈服强度为１５１０ＭＰａ，１０００℃时表现出良好的压缩塑性［８］由Ｙａｏ等人设计的难熔多组元合金ＮｂＴｉＭｏＶＴａ，具有优异的抗氧化性能，同时很高的屈服强度和延展性［９］ 由Ｚｈｏｕ等人设计的Ｈｆ０ ５Ｎｂ０ ５Ｔａ０ ５ＴｉＺｒ难熔多组元合金与Ｎｉ基高熵合金（ＡｌｘＣｏＣｒＦｅＮ）相比具有更优异的抗腐蚀性能［１０］ Ｎａｍ等人研究发现ＮｂＭｏＴａＷ难熔高熵合金高温下表现出优异的抗氧化性能和较高的纳米压痕硬度［１１］关于难熔多组元合金在核反应堆方面的研究最早可以追溯到２０１３年 Ｎａｇａｓｅ等人研究发现具有单相ＢＣＣ结构的ＨｆＮｂＺｒ合金在２ＭｅＶ电子辐照后具有较高的结构稳定性（２５℃为１０ｄｐａ和－１７０℃为５０ｄｐａ），表现出很高的抗辐照性能［１２］ 随后Ｅｇａｍｉ等人对ＨｆＮｂＺｒ合金的抗辐照性能进行理论分析，并得出结论高熵合金由于具有“自愈效应”从而表现出较高的抗辐照性能［１３］从此，高熵合金在核反应堆中应用引起了人们的兴趣 ２０１９年Ｌｕ等人设计了单相ＢＣＣ结构的Ｔｉ２ＨｆＺｒＶ０ ５Ｍｏ０ ２难熔多组元合金，在Ｈｅ离子辐照后几乎没有辐照硬化，辐照后晶格常数降低，表现出良好的抗辐照性能［１４］ ＴｉＮｂＴａＶ难熔多组元合金具有单相ＢＣＣ结构，在２ＭｅＶ的Ｖ＋辐照后显示出０ ６６ＧＰａ（８％）的辐照硬化［１５］ ＴｉＺｒＮｂＨｆＴａ难熔多组元合金在５ＭｅＶ的Ｈｅ离子辐照后，提高了屈服和极限拉伸强度，延展性没有损失，具有优异的抗损伤性能［１６］ 上述合金虽然具有良好的抗辐照性能，但都含有高激活元素Ｎｂ和Ｍｏ，不利于其在核反应堆中应用 由Ｋａｒｅｅｒ等人设计的ＴｉＣｒＴａＶ低活化难熔多组元合金，在Ｖ离子辐照后，几乎没有辐照硬化，具有优异的抗辐照性能［１５］ 该合金含有多个结构相，与其他合金相比具有更复杂的内部环境，同时该合金具有较低的热中子吸收截面，因此十分具有研究价值目前对该合金的研究主要集中于合金的显微组织结构、压痕硬度、模量以及抗辐照性能，缺乏电子层面的深入研究 本文采用准随机近似（ＳＱＳ）方法建立了合金中两种ＢＣＣ（ＢＣＣ１Ｔｉ２３Ｃｒ２５Ｔａ２３Ｖ２９和ＢＣＣ２Ｔｉ２７Ｃｒ２７Ｔａ１８Ｖ２８）的结构模型，利用第一性原理计算，系统计算了两种结构的结构稳定性、力学性能、德拜温度、电子结构和布居分析，从原子层面分析了他们的力学性能和原子间成键特性，为低活化核用多组元合金的研究提供一定的理论依据２ 计算方法和模型基于第一性原理密度泛函理论，采用ＣＡＳＴＥＰ软件包对ＴｉＣｒＴａＶ合金两种ＢＣＣ结构进行计算 利用准随机近似（ＳＱＳ）方法搭建３×３×４的超胞来模拟原子化学无序结构 计算结构模型如图１所示 交换关联泛函采用广义梯度近似（ＧＧＡ）下的ＰＢＥ泛函，采用超软赝势处理电子－离子关系 计算平面波函数的截断能为４５０ｅＶ，布里渊区积分采用２×２×２Ｍｏｎｋｈｏｒｓｔ－Ｐａｃｋ形式的高对称特殊Ｋ点方法 自洽计算（ＳＣＦ）采用Ｐｕｌａｙ密度混合法，自洽计算的误差为１ ０×１０－６ｅＶ／ａｔｏｍ，体系总能量低于１ ０×１０－５ｅＶ／ａｔｏｍ，每个原子上的最大力小于０ ０３ｅＶ／?，公差偏移小于０ ００１? 对于所含原子其赝势分别采用Ｔｉ３ｄ２４ｓ２，Ｃｒ３ｄ５４ｓ１，Ｔａ５ｄ３６ｓ２，Ｖ３ｄ３４ｓ２３ 结果与讨论３ １ 晶格常数及结构稳定性为研究不同ＢＣＣ结构相的形成能力和稳定性，计算了生成焓（Ｈｆ）和结合能（Ｅｃｏｈ） 计算公式如下：Ｈｆ＝１Ｎｉ（Ｅｔｏｔ－∑ｎｉＥｉｓｏｌｉｄ）（１）第４０卷茹金豆，等：ＢＣＣ结构ＴｉＣｒＴａＶ合金的电子结构和力学性能的第一性原理计算第１期图１ 采用ＳＱＳ方法建立的ＢＣＣ结构计算模型：（ａ）ＢＣＣ１，（ｂ）ＢＣＣ２Ｆｉｇ １ ＢＣＣｓｔｒｕｃｔｕｒｅｍｏｄｅｌｓｅｓｔａｂｌｉｓｈｅｄｂｙＳＱＳｍｅｔｈｏｄ：（ａ）ＢＣＣ１，（ｂ）ＢＣＣ２Ｅｃｏｈ＝１Ｎｉ（Ｅｔｏｔ－∑ｎｉＥｉａｔｏｍ）（２）式中，Ｅｔｏｔ为总能量，Ｅｉｓｏｌｉｄ和Ｅｉａｔｏｍ分别是单晶完全弛豫后单个原子的能量和第ｉ个组分的孤立原子的能量 Ｎｉ和ｎｉ分别为超胞中的全部原子数和第ｉ个组分的摩尔分数表１为晶格优化后两种ＢＣＣ结构的晶格常数、结构总能、生成焓和结合能 晶格常数的计算值与Ａ Ｋａｒｅｅｒ的实验值相接近，表明计算参数较为合理 生成焓越小，反应越容易进行 结合能越低，晶体越稳定［１７］ 两种结构的生成焓和结合能都是负值，表明两种结构相的生成过程均为放热反应，而且两种结构都是稳定的 其中ＢＣＣ１结构相的生成焓和结合能更低，表明在相同反应条件下ＢＣＣ１结构更容易生成且更稳定 表１ 两种ＢＣＣ结构的晶格常数（ａ，?）、结构总能（Ｅｔｏｔ，ｅＶ／ａｔｏｍ）、生成焓（Ｈｆ，ｅＶ／ａｔｏｍ）和结合能（Ｅｃｏｈ，ｅＶ／ａｔｏｍ）Ｔａｂｌｅ１ Ｌａｔｔｉｃｅｃｏｎｓｔａｎｔｓ（ａ，?），ｔｏｔａｌｅｎｅｒｇｉｅｓ（Ｅｔｏｔ，ｅＶ／ａｔｏｍ），ｆｏｒｍａｔｉｏｎｅｎｔｈａｌｐｉｅｓ（Ｈｆ，ｅＶ／ａｔｏｍ）ａｎｄｃｏｈｅｓｉｖｅｅｎｅｒｇｉｅｓ（Ｅｃｏｈ，ｅＶ／ａｔｏｍ）ｆｏｒｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓＳｔｒｕｃｔｕｒｅａＰｒｅｓｅｎｔａＥｘｐ［１５］ＥｔｏｔＨｆＥｃｏｈＢＣＣ１３ ０９３３ １１０－３２１５－１７ １６６－２５ ９３０ＢＣＣ２３ ０９１３ １１０－３３６５－５ ０６９－１３ ７４５３ ２ 力学性能表２所示为两种ＢＣＣ结构的弹性常数 通常体心立方晶体有三个独立的弹性常数 这里的ＳＱＳ模型生成了９个弹性常数 因此，取弹性常数的平均值以减少ＳＱＳ模型对称性低、尺寸大的影响 弹性常数平均值的计算结果如表３所示ｃ１１＝ｃ′１１＋ｃ′２２＋ｃ′３３３，ｃ１２＝ｃ′１２＋ｃ′１３＋ｃ′２３３，ｃ４４＝ｃ′４４＋ｃ′５５＋ｃ′６６３（３）表２ 两种ＢＣＣ结构的弹性常数（ｃｉｊ，ＧＰａ）Ｔａｂｌｅ２ Ｅｌａｓｔｉｃｃｏｎｓｔａｎｔｓ（ｃｉｊ，ＧＰａ）ｏｆｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓＳｔｒｕｃｔｕｒｅｃ′１１ｃ′２２ｃ′３３ｃ′１２ｃ′１３ｃ′２３ｃ′４４ｃ′５５ｃ′６６ＢＣＣ１２４５２３８２３６１３８１３７１３９３３３０１８ＢＣＣ２２４５２４２２４６１３１１３３１３５２９３２２５表３ 两种ＢＣＣ结构的多晶弹性常数（ｃｉｊ，ＧＰａ）、柯西压力（ｃ１２－ｃ４４，ＧＰａ）、体积模量（Ｂ，ＧＰａ）、剪切模量（Ｇ，ＧＰａ）、Ｇ／Ｂ比、杨氏模量（Ｅ，ＧＰａ）、泊松比（ν，ＧＰａ）、各向异性分数因子（ＡＵ）和硬度（ＨＶ，ＧＰａ）Ｔａｂｌｅ３ Ｅｌａｓｔｉｃｃｏｎｓｔａｎｔｓ（ｃｉｊ，ＧＰａ），Ｃａｕｃｈｙｐｒｅｓｓｕｒｅｓ（ｃ１２－ｃ４４，ＧＰａ），ｂｕｌｋｍｏｄｕｌｉ（Ｂ，ＧＰａ），ｓｈｅａｒｍｏｄｕｌｉ（Ｇ，ＧＰａ），Ｇ／Ｂｒａｔｉｏｓ，Ｙｏｕｎｇ’ｓｍｏｄｕｌｉ（Ｅ，ＧＰａ），Ｐｏｉｓｓｏｎ＇ｓｒａｔｉｏｓ（ν），ａｎｉｓｏｔｒｏｐｙｆｒａｃｔｉｏｎｆａｃｔｏｒｓ（ＡＵ）ａｎｄｈａｒｄｎｅｓｓ（ＨＶ，ＧＰａ）ｏｆｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓＳｔｒｕｃｔｕｒｅｃ１１ｃ１２ｃ４４ｃ１２－ｃ４４ＢＧＧ／ＢＥνＡＵＨＶＢＣＣ１２４０１３８２７１０５１７１３４０ ２０９５０ ４１０ ９０１ ７５ＢＣＣ２２４６１３３２９１０４１７０３７０ ２２１０５０ ６５０ ６３２ １４ 通常，将Ｂｏｒｎ规则作为确定晶体力学稳定性的决定性条件：ｃ１１＞ｃ１２，ｃ４４＞０，（ｃ１１＋２ｃ１２）＞０（４）计算结果表明，这两个ＢＣＣ结构均具有力学稳定性 ＢＣＣ２结构的ｃ１１较大，表明其在Ｘ轴上具有较强的抗压缩能力［１８］ 同时ＢＣＣ２的ｃ４４值较大，也表明其在（１００）面抗单斜剪切能力较强由Ｖｏｉｇｔ－Ｒｅｕｓｓ－Ｈｉｌｌ（ＶＲＨ）计算了ＴｉＣｒＴａＶ合金两种ＢＣＣ结构的多晶弹性模量 计算公式如下：Ｂ＝ＢＶ＋ＢＲ２（５）Ｇ＝ＧＶ＋ＧＲ２（６）Ｅ＝９ＢＧ３Ｂ＋Ｇ（７）第４０卷原 子 与 分 子 物 理 学 报 第１期ν＝３Ｂ－Ｅ６Ｂ（８）式中，ＢＶ（ＧＶ）和ＢＲ（ＧＲ）分别表示体积模量和剪切模量的Ｖｏｉｇｔ和Ｒｅｕｓｓ估计值体积模量Ｂ，反应材料的抗压缩变形的能力，ＢＣＣ１的体积模量较大，说明其抗压性较强，强度更高 剪切模量Ｇ反应材料的抗剪切变形能力．杨氏模量Ｅ反应材料的刚度，杨氏模量越大，材料的硬度越好［１９］ ＢＣＣ２结构的Ｇ和Ｅ均大于ＢＣＣ１结构，表明ＢＣＣ２结构具有更高的抗剪切变形能力和刚度 柯西压力（ｃ１２－ｃ４４）可以评价材料中的金属键特性和共价键特性［２０］ （ｃ１２－ｃ４４）＞０表明材料具有韧性表现为金属键特性 Ｇ／Ｂ比和泊松比ν反映了材料的延性，当Ｇ／Ｂ＜０ ５７，ν＞０ ２６时，材料具有延展性，否则为脆性［２１］ 由计算结果可知两种结构都具有延展性 （ｃ１２－ｃ４４）和Ｇ／Ｂ的计算结果表明，ＢＣＣ１具有更好的金属键特性和延展性，而ν的计算结果相反 它们之间的差异可能是由于弹性各向异性的影响，导致多晶弹性模量的不确定性，从而造成变化趋势的不同ＨＶ＝０ ９２（Ｇ／Ｂ）１ １３７Ｇ０ ７０８（９）由式（９）可知，Ｇ／Ｂ比越大，硬度越大，ＢＣＣ２结构具有更高的硬度为了研究两种结构的弹性各向异性，分别计算了他们的各向异性因子ＡＵ 其表达式如下：ＡＵ＝５ＧＶＧＲ＋ＢＶＢＲ－６（１０）ＡＵ＝０说明晶体是各向同性的，ＡＵ对０的偏离程度越大，晶体的各向异性程度越大 弹性各向异性因子的计算结果如表３ 两种结构的ＡＵ均大于０，说明均存在弹性各向异性 但ＢＣＣ１结构的各向异性因子更大，表明ＢＣＣ１结构具有较强的弹性各向异性为了更直观地理解弹性各向异性，我们计算了弹性各向异性三维图形 其表达式如下：Ｂｉｊｋ＝１／［（ｓ１１＋ｓ１２＋ｓ１３）ｌ２ｉ１＋（ｓ１２＋ｓ２２＋ｓ２３）ｌ２ｊ２＋（ｓ１３＋ｓ２３＋ｓ３３）ｌ２ｋ３］（１１）Ｅｉｊｋ＝１／［ｓ１１－２（ｓ１１－ｓ２２－１／２ｓ４４）（ｌ２ｉ１ｌ２ｊ２＋ｌ２ｊ２ｌ２ｋ３＋ｌ２ｉ１ｌ２ｋ３）］（１２）Ｇｉｊｋ＝１／［ｓ４４＋４（ｓ１１－ｓ１２－ｓ４４／２）（ｌ２ｉ１ｌ２ｊ２＋ｌ２ｊ２ｌ２ｋ３＋ｌ２ｉ１ｌ２ｋ３）］（１３）其中Ｂｉｊｋ为晶体在＜ｉｊｋ＞方向上的体积模量 Ｅｉｊｋ为晶体在＜ｉｊｋ＞方向上的杨氏模量 Ｇｉｊｋ为晶体沿＜ｉｊｋ＞方向的剪切模量 ｓｉｊ为弹性柔度系数 ｌｉ１，ｌｊ２，ｌｋ３为球坐标下的方向余弦图２为三种弹性模量的三维表面结构 Ｂ的表面接近球形，表明Ｂ是各向同性的 Ｅ和Ｇ的表面偏离球形，表明Ｅ和Ｇ存在各向异性，且Ｅ的各向异性程度大于Ｇ 杨氏模量Ｅ在［１００］方向取最大值，在［１１１］方向取最小值 但剪切模量Ｇ，情况相反，原因是对于立方晶体，Ｅ的各向异性取决于（ｌ１２ｌ２２＋ｌ２２ｌ３２＋ｌ３２ｌ１２） （ｓ１１－ｓ１２－ｓ４４／２）＜０，Ｅ在［１００］方向取最大值，在［１１１］方向取最小值 如图所示ＢＣＣ１结构的体积模量Ｂ更大，ＢＣＣ２结构的杨氏模量Ｅ和剪切模量Ｇ更大，这与上述多晶弹性模量计算结果一致图２ 弹性模量的三维表面结构Ｆｉｇ ２ Ｔｈｒｅｅ－ｄｉｍｅｎｓｉｏｎａｌｓｕｒｆａｃｅｃｏｎｓｔｒｕｃｔｉｏｎｓｏｆｅｌａｓｔｉｃｍｏｄｕｌｕｓ３ ３ 德拜温度Ｄｅｂｙｅ温度θＤ与平均声速ｖｍ的关系如下：θＤ＝ｈｋＢ３ｎ４πＮＡρ()[]Ｍ１／３νｍ（１４）νｍ＝１３２ν３ｔ＋１ν３()[]ｌ－１／３（１５）νｔ＝Ｇ槡ρ（１６）νｌ＝３Ｂ＋４Ｇ３槡ρ（１７）其中ｈ为普朗克常数，ｋＢ为玻尔兹曼常数，ＮＡ为阿第４０卷茹金豆，等：ＢＣＣ结构ＴｉＣｒＴａＶ合金的电子结构和力学性能的第一性原理计算第１期伏伽德罗常数，Ｍ为分子量，ρ为密度，ｎ为单胞内原子数，ｖｔ和ｖｌ为横向和纵向声速Ｇｒüｎｅｉｓｅｎ参数γａ计算公式如下：γａ＝３２３ν２ｌ－４ν２ｔν２ｌ＋２ν２()ｔ（１８）声速、德拜温度和γａ的计算结果如表４所示 两种结构的横向声速ｖｔ远远小于纵向声速ｖｌ，原因是其体积模量Ｂ总是大于剪切模量Ｇ ＢＣＣ２的德拜温度和声速总是大于ＢＣＣ１，表明ＢＣＣ２结构的键合强度较强，热稳定性更好［２２］ Ｇｒüｎｅｉｓｅｎ参数γａ表示温度或压力对晶体体积的影响 ＢＣＣ２结构的γａ较小，表明ＢＣＣ２结构具有更好的热力学稳定性．表４ 晶胞密度（ρ，ｇ／ｃｍ３），横向声速（ｖｔ，ｍ／ｓ），纵向声速（ｖｌ，ｍ／ｓ），平均声速（ｖｍ，ｍ／ｓ），德拜温度（θＤ，Ｋ）和Ｇｒüｎｅｉｓｅｎ参数（γａ）Ｔａｂｌｅ４ Ｃｅｌｌｄｅｎｓｉｔｉｅｓ（ρ，ｇ／ｃｍ３），ｔｒａｎｓｖｅｒｓｅｓｏｕｎｄｖｅｌｏｃｉ ｔｉｅｓ（ｖｔ，ｍ／ｓ），ｌｏｎｇｉｔｕｄｉｎａｌｓｏｕｎｄｖｅｌｏｃｉｔｉｅｓ（ｖｌ，ｍ／ｓ），ａｖｅｒａｇｅｓｏｕｎｄｖｅｌｏｃｉｔｉｅｓ（ｖｍ，ｍ／ｓ），Ｄｅｂｙｅｔｅｍｐｅｒａｔｕｒｅｓ（θＤ，Ｋ）ａｎｄＧｒüｎｅｉｓｅｎｐａｒａｍｅｔｅｒｓ（γａ）ＳｔｒｕｃｔｕｒｅρｖｔｖｌｖｍθＤγａＢＣＣ１８ ７３５１９７３４９７７１３１５１５８２ ７１ＢＣＣ２８ ２２７２１３４５１７１１４１６１８８２ ５９３ ４ 电子性质图３为ＢＣＣ结构的态密度（ＤＯＳ）图，图中的黑色虚线为费米能级 计算结果表明两种ＢＣＣ结构的ＤＯＳ曲线图十分相似，故取其中一个ＤＯＳ图可以有效分析两个ＢＣＣ结构中原子成键情况 如图所示，费米能级处ＤＯＳ曲线大于０，表明两种ＢＣＣ结构都具有金属键特性［２３］ 在－１０ｅＶ和５ｅＶ之间，Ｔｉ、Ｃｒ、Ｖ、Ｔａ四种原子的电子轨道明显重叠，表明两种ＢＣＣ结构中存在严重的电子杂化形成共价键 在能量较低的区域有几个又尖又小的峰，这些峰主要由Ｔｉ、Ｃｒ和Ｖ的ｓ轨道和ｐ轨道贡献，没有电子杂化，没有Ｔａ原子参与，表明金属键只存在于Ｔｉ、Ｃｒ和Ｖ原子间，而Ｔａ原子只生成共价键［２４］ 图４为电荷密度图 如图所示，Ｔａ原子和Ｔｉ、Ｃｒ、Ｖ原子之间有比较强烈的电子云重叠，表明共价键的存在［２５］ 同时，Ｔｉ、Ｃｒ和Ｖ原子周围电子云呈圆形分布，他们之间的电子云略有重叠，表明他们之间同时存在共价键和金属键［２４］ ３ ５ 布居分析表５为两种ＢＣＣ结构的Ｍｕｌｌｉｋｅｎ原子布居平图３ ＢＣＣ结构的态密度图Ｆｉｇ ３ ＴｈｅｄｅｎｓｉｔｉｅｓｏｆｓｔａｔｅｓｏｆＢＣＣｓｔｒｕｃｔｕｒｅ图４ 两种ＢＣＣ结构的电荷密度图（ａ）ＢＣＣ１；（ｂ）ＢＣＣ２Ｆｉｇ ４ ＣｈａｒｇｅｄｅｎｓｉｔｙｄｉａｇｒａｍｓｏｆｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓ（ａ）ＢＣＣ１；（ｂ）ＢＣＣ２均值 Ｔｉ和Ｖ原子的Ｍｕｌｌｉｋｅｎ布居均为正值，表示元素成键后失去电子，而且Ｔｉ的电荷值更大，表明其失电子能力更大［２６］ Ｃｒ和Ｔａ原子的Ｍｕｌｌｉｋｅｎ布居均为负值，表示得到电子，而且Ｔａ的负电荷值更小，表明其得电子能力更强 表６为两种结构中原子的重叠布居和键长 原子间重叠布居正值代表共价键的形成，负值代表反键（金属键）的形成［２７］ 计算结果表明两种结构中同时存在共价键和金属键 Ｔａ原子只形成共价键，而且其共价键强度最高，导致塑性和延展性降低 Ｃｒ、Ｖ和Ｔｉ原子之间同时存在共价键和金属键，金属键的形成增强原子间斥力，有利于提高塑性和延展性４ 结 论基于密度泛函理论的第一性原理方法研究了ＴｉＣｒＴａＶ合金中两种ＢＣＣ结构的形成能力、弹性第４０卷原 子 与 分 子 物 理 学 报 第１期表５ 两种ＢＣＣ结构的Ｍｕｌｌｉｋｅｎ原子布居平均值Ｔａｂｌｅ５ ＡｖｅｒａｇｅＭｕｌｌｉｋｅｎｐｏｐｕｌａｔｉｏｎｓｏｆｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓＳｔｒｕｃｔｕｒｅＡｔｏｍｓｐｄｆＴｏｔａｌＣｈａｒｇｅ／ｅＶＢＣＣ１Ｔｉ２ ２３６ ５５２ ６８０１１ ４６０ ５４Ｃｒ２ ３８６ ５９５ １００１４ ０７－０ ０７Ｔａ０ ６０１ ５７３ ８８０６ ０５－１ ０５Ｖ２ ２８６ ５３３ ７６０１２ ５７０ ４３ＢＣＣ２Ｔｉ２ １９６ ５６２ ６８０１１ ４３０ ５７Ｃｒ２ ３７６ ７０５ １００１４ １６－０ １６Ｔａ０ ６１１ ６７３ ８９０６ １７－１ １７Ｖ２ ２７６ ５８３ ７７０１２ ６２０ ３８表６ 两种ＢＣＣ结构的原子的重叠布居和键长Ｔａｂｌｅ６ ＡｔｏｍｉｃＯｖｅｒｌａｐＰｏｐｕｌａｔｉｏｎｓ（ｅ）ａｎｄＢｏｎｄＬｅｎｇｔｈｓ（）ｏｆｔｈｅｔｗｏＢＣＣｓｔｒｕｃｔｕｒｅｓＳｔｒｕｃｔｕｒｅＢＣＣ１ＢＣＣ２ｂｏｎｄＯｖｅｒｌａｐＰｏｐｕｌａｔｉｏｎＢｏｎｄｌｅｎｇｔｈｂｏｎｄＯｖｅｒｌａｐＰｏｐｕｌａｔｉｏｎＢｏｎｄｌｅｎｇｔｈＴａ６－Ｔａ１３０ ７１２ ８４１Ｔａ７－Ｔａ８０ ６５２ ８５０Ｔａ１２－Ｃｒ１２０ ３９２ ７２１Ｔａ１０－Ｃｒ１６０ ４３２ ７３１Ｔａ１０－Ｖ１２０ ３０２ ７５Ｔａ１１－Ｖ１０ ３０２ ７１８Ｔａ１４－Ｔｉ８０ ２３２ ８２Ｔａ７－Ｔｉ７０ ２４２ ８２０Ｃｒ１－Ｃｒ１８０ １９２ ６８Ｃｒ１９－Ｃｒ２００ ２２ ６２２Ｖ６－Ｖ１００ １８２ ７８Ｖ７－Ｖ８０ １９２ ４８９Ｖ６－Ｃｒ２０ １７２ ６１Ｖ１４－Ｃｒ２００ １８２ ６９３Ｔｉ５－Ｖ９０ １６２ ７６５Ｔｉ１３－Ｖ９０ １５２ ７１６Ｔｉ１－Ｔｉ２０ １６２ ７９５Ｔｉ１１－Ｔｉ１２０ １３２ ７９２Ｔｉ７－Ｃｒ１１０ １２２ ６６４Ｔｉ１１－Ｃｒ４０ １２２ ５９５Ｃｒ３－Ｃｒ８－０ ０１２ ９５１Ｃｒ２－Ｃｒ７－０ ０４２ ９９３Ｖ１６－Ｃｒ１１－０ ０２２ ８２０Ｖ１１－Ｃｒ１４－０ ０２２ ８３３Ｔｉ１－Ｃｒ６－０ ０３２ ８３７Ｔｉ１８－Ｃｒ２０－０ ０１２ ９５０性能、德拜温度、电子结构和布居分析 计算的Ｈｆ和Ｅｃｏｈ表明，ＢＣＣ１结构更稳定，更容易形成两种结构均具有金属性能和延展性 ＢＣＣ１结构的Ｂ较大，强度更高 ＢＣＣ２结构的Ｅ，Ｇ和ＨＶ更大，表明其具有更好的抗剪切能力和刚度 两种结构都存在不同程度的弹性各向异性 其中Ｅ的各向异性程度最大，而Ｂ为各向同性 Ｅ在［１００］方向取最大值，［１１１］方向取最小值 Ｄｅｂｙｅ温度和Ｇｒüｎｅｉｓｅｎ参数γａ表明ＢＣＣ２相具有较强的键合作用和更好的热力学稳定性 电子结构和布居分析结果表明两种结构中同时存在共价键和金属键 Ｔａ原子只形成共价键，而且其强度最大 金属键仅存在于Ｔｉ、Ｃｒ和Ｖ原子之间 元素成键后Ｃｒ和Ｔａ原子得到电子，Ｔｉ和Ｖ原子失去电子，而且Ｔａ原子和Ｔｉ原子的得失电子能力更强 参考文献：［１］ ＤａｓＳ，ＲｏｂｉＰＳ ＡｎｏｖｅｌｒｅｆｒａｃｔｏｒｙＷＭｏＶＣｒＴａｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙｐｏｓｓｅｓｓｉｎｇｆｉｎｅｃｏｍｂｉｎａｔｉｏｎｏｆｃｏｍｐｒｅｓｓｉｖｅｓｔｒｅｓｓ－ｓｔｒａｉｎａｎｄｈｉｇｈｈａｒｄｎｅｓｓｐｒｏｐｅｒｔｉｅｓ［Ｊ］．Ａｄｖ ＰｏｗｄｅｒＴｅｃｈｎｏｌ ，２０２０，３１：４６１９ ［２］ ＣｈｅｎＹ，ＸｕＺ，ＷａｎｇＭ，ｅｔａｌ Ａｓｉｎｇｌｅ－ｐｈａｓｅＶ０ ５Ｎｂ０ ５ＺｒＴｉｒｅｆｒａｃｔｏｒｙｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙｗｉｔｈｏｕｔ ｓｔａｎｄｉｎｇｔｅｎｓｉｌｅｐｒｏｐｅｒｔｉｅｓ［Ｊ］ Ｍａｔｅｒ Ｓｃｉ Ｅｎｇ Ａ，２０２０，７９２：１３９７７４［３］ ＪｏｓｅｐｈＪ，ＨａｇｈｄａｄｉＮ，ＳｈａｍｌａｙｅＫ，ｅｔａｌ ＴｈｅｓｌｉｄｉｎｇｗｅａｒｂｅｈａｖｉｏｕｒｏｆＣｏＣｒＦｅＭｎＮｉａｎｄＡｌｘＣｏＣｒＦｅＮｉｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓａｔｅｌｅｖａｔｅｄｔｅｍｐｅｒａｔｕｒｅｓ［Ｊ］Ｗｅａｒ，２０１９，４２８：３２［４］ ＣｅｍｉｎＦ，ＪｉｍｅｎｅｚＭＪＭ，ＬｅｉｄｅｎｓＬＭ，ｅｔａｌ Ａｔｈｅｒ ｍｏｄｙｎａｍｉｃｓｔｕｄｙｏｎｐｈａｓｅｆｏｒｍａｔｉｏｎａｎｄｔｈｅｒｍａｌｓｔａｂｉｌｉｔｙｏｆＡｌＳｉＴａＴｉＺｒｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙｔｈｉｎｆｉｌｍｓ［Ｊ］Ｊ Ａｌｌｏｙ Ｃｏｍｐ ，２０２０，８３８：１５５５８０［５］ ＬｉＬＣ，ＬｉＭＸ，ＬｉｕＭ，ｅｔａｌ Ｅｎｈａｎｃｅｄｏｘｉｄａｔｉｏｎｒｅ ｓｉｓｔａｎｃｅｏｆＭｏＴａＴｉＣｒＡｌｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓｂｙｒｅｍｏｖａｌｏｆＡｌ［Ｊ］ Ｓｃｉ ＣｈｉｎａＭａｔｅｒ ，２０２１，６４：２２３ ［６］ ＹｕＹ，ＷａｎｇＪ，ＹａｎｇＪ，ｅｔａｌ Ｃｏｒｒｏｓｉｖｅａｎｄｔｒｉｂｏｌｏｇｉ ｃａｌｂｅｈａｖｉｏｒｓｏｆＡｌＣｏＣｒＦｅＮｉ－Ｍｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓｕｎｄｅｒ９０ｗｔ ％Ｈ２Ｏ２ｓｏｌｕｔｉｏｎ［Ｊ］ Ｔｒｉｂｏｌ Ｉｎｔ ，２０１９，１３１：２４［７］ ＭｉｒａｃｌｅＤＢ，ＳｅｎｋｏｖＯＮ Ａｃｒｉｔｉｃａｌｒｅｖｉｅｗｏｆｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓａｎｄｒｅｌａｔｅｄｃｏｎｃｅｐｔｓ［Ｊ］ＡｃｔａＭａｔｅｒ，２０１７，１２２：４４８［８］ ＳｅｎｋｏｖＯＮ，ＲａｏＳ，ＣｈａｐｕｔＫＪ，ｅｔａｌ ＣｏｍｐｏｓｉｔｉｏｎａｌｅｆｆｅｃｔｏｎｍｉｃｒｏｓｔｒｕｃｔｕｒｅａｎｄｐｒｏｐｅｒｔｉｅｓｏｆＮｂＴｉＺｒ－ｂａｓｅｄｃｏｍｐｌｅｘｃｏｎｃｅｎｔｒａｔｅｄａｌｌｏｙｓ［Ｊ］ ＡｃｔａＭａｔｅｒ ，２０１８，１５１：２０１ ［９］ ＹａｏＨＷ，ＱｉａｏＪＷ，ＨａｗｋＪＡ，ｅｔａｌ Ｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆｒｅｆｒａｃｔｏｒｙｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙｓ：Ｅｘｐｅｒｉｍｅｎｔｓａｎｄｍｏｄｅｌｉｎｇ［Ｊ］ Ｊ Ａｌｌｏｙ Ｃｏｍｐ ，２０１７，６９６：１１３９［１０］ ＺｈｏｕＱ，ＳｈｅｉｋｈＳ，ＯｕＰ，ｅｔａｌ ＣｏｒｒｏｓｉｏｎｂｅｈａｖｉｏｒｏｆＨｆ０ ５Ｎｂ０ ５Ｔａ０ ５Ｔｉ１ ５Ｚｒｒｅｆｒａｃｔｏｒｙｈｉｇｈ－ｅｎｔｒｏｐｙｉｎａｑｕｅｏｕｓｃｈｌｏｒｉｄｅｓｏｌｕｔｉｏｎｓ［Ｊ］ Ｅｌｅｃｔｒｏｃｈｅｍ Ｃｏｍｍｕｎ ，２０１９，９８：６３ ［１１］ ＮａｍＳ，ＳｏｎｇＹ，ＳｏｎＭ，ｅｔａｌ Ｈｉｇｈ－ｔｅｍｐｅｒａｔｕｒｅ第４０卷茹金豆，等：ＢＣＣ结构ＴｉＣｒＴａＶ合金的电子结构和力学性能的第一性原理计算第１期ｍｅｃｈａｎｉｃａｌａｎｄｔｈｅｒｍｏｃｈｅｍｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆＮｂＭｏ ＴａＷｒｅｆｒａｃｔｏｒｙｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙｃｏａｔｉｎｇｓｐｒｏｄｕｃｅｄｖｉａＮａｎｏＰａｒｔｉｃｌｅＤｅｐｏｓｉｔｉｏｎＳｙｓｔｅｍ（ＮＰＤＳ）［Ｊ］ Ｉｎｔ Ｊ Ｒｅｆｒａｃｔ Ｍｅ Ｈ ，２０２１，９９：１０５５９４ ［１２］ ＮａｇａｓｅＴ，ＡｎａｄａＳ，ＲａｃｋＰＤ，ｅｔａｌ ＭｅＶｅｌｅｃｔｒｏｎ－ｉｒｒａｄｉａｔｉｏｎ－ｉｎｄｕｃｅｄｓｔｒｕｃｔｕｒａｌｃｈａｎｇｅｉｎｔｈｅｂｃｃｐｈａｓｅｏｆＺｒ Ｈｆ Ｎｂａｌｌｏｙｗｉｔｈａｎａｐｐｒｏｘｉｍａｔｅｌｙｅｑｕｉａｔｏｍｉｃｒａｔｉｏ［Ｊ］ Ｉｎｔｅｒｍｅｔａｌｌｉｃｓ，２０１３，３８：７０［１３］ ＥｇａｍｉＴ，ＧｕｏＷ，ＲａｃｋＰ，ｅｔａｌ Ｉｒｒａｄｉａｔｉｏｎｒｅｓｉｓｔ ａｎｃｅｏｆｍｕｌｔｉｃｏｍｐｏｎｅｎｔａｌｌｏｙｓ［Ｊ］ Ｍｅｔａｌｌ Ｍａｔｅｒ Ｔｒａｎｓ Ａ，２０１４，４５：１８０ ［１４］ ＬｕＹ，ＨｕａｎｇＨ，ＧａｏＸ，ｅｔａｌ Ａｐｒｏｍｉｓｉｎｇｎｅｗｃｌａｓｓｏｆｉｒｒａｄｉａｔｉｏｎｔｏｌｅｒａｎｔｍａｔｅｒｉａｌｓ：Ｔｉ２ＺｒＨｆＶ０ ５Ｍｏ０ ２ｈｉｇｈ－ｅｎｔｒｏｐｙａｌｌｏｙ［Ｊ］ Ｊ Ｍａｔｅｒ Ｓｃｉ Ｔｅｃｈｎｏｌ ，２０１９，３５：３６９ ［１５］ ＫａｒｅｅｒＡ，ＷａｉｔｅＪＣ，ＬｉＢ，ｅｔａｌ Ｓｈｏｒｔｃｏｍｍｕｎｉｃａ ｔｉｏｎ：‘Ｌｏｗａｃｔｉｖａｔｉｏｎ，ｒｅｆｒａｃｔｏｒｙ，ｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓｆｏｒｎｕｃｌｅａｒａｐｐｌｉｃａｔｉｏｎｓ’［Ｊ］ Ｊ Ｎｕｃｌ Ｍａｔｅｒ ，２０１９，５２６：１５１７４４ ［１６］ ＭｏｓｃｈｅｔｔｉＭ，ＸｕＡ，ＳｃｈｕｈＢ，ｅｔａｌ Ｏｎｔｈｅｒｏｏｍ－ｔｅｍｐｅｒａｔｕｒｅｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆａｎｉｏｎ－ｉｒｒａｄｉａｔｅｄＴｉＺｒＮｂＨｆＴａｒｅｆｒａｃｔｏｒｙｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙ［Ｊ］ ＪＯＭ，２０２０，７２：１３０ ［１７］ ＨｕＹＬ，ＢａｉＬＨ，ＴｏｎｇＹＧ，ｅｔａｌ Ｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｃａｌｃｕｌａｔｉｏｎｉｎｖｅｓｔｉｇａｔｉｏｎｏｆＮｂＭｏＴａＷｂａｓｅｄｒｅｆｒａｃｔｏｒｙｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓ［Ｊ］ Ｊ Ａｌｌｏｙ Ｃｏｍｐ ，２０２０，８２７：１５３９６３ ［１８］ ＭａＬ，ＤｕａｎＹ，ＬｉＲ Ｐｈａｓｅｓｔａｂｉｌｉｔｙ，ａｎｉｓｏｔｒｏｐｉｃｅｌａｓｔｉｃｐｒｏｐｅｒｔｉｅｓａｎｄｅｌｅｃｔｒｏｎｉｃｓｔｒｕｃｔｕｒｅｓｏｆＣ１５－ｔｙｐｅＬａｖｅｓｐｈａｓｅｓＺｒＭ２（Ｍ＝Ｃｒ，ＭｏａｎｄＷ）ｆｒｏｍｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓｃａｌｃｕｌａｔｉｏｎｓ［Ｊ］ Ｐｈｉｌｏｓ Ｍａｇ ，２０１７，９７：２４０６［１９］ ＷｕＨ，ＨｕａｎｇＳ，ＺｈｕＣ，ｅｔａｌＩｎｆｌｕｅｎｃｅｏｆＣｒｃｏｎｔｅｎｔｏｎｔｈｅｍｉｃｒｏｓｔｒｕｃｔｕｒｅａｎｄｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆＣｒｘＦｅＮｉＣｕｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓ［Ｊ］Ｐｒｏｇ Ｎａｔ Ｓｃｉ－Ｍａｔｅｒ ，２０２０，３０：２３９ ［２０］ ＬｉｕＸ，ＳｈａＧ，ＷｕＱ，ｅｔａｌ Ｐｈａｓｅｓｔａｂｉｌｉｔｙｏｆａｎｈｉｇｈ－ｅｎｔｒｏｐｙＡｌ－Ｃｒ－Ｆｅ－Ｎｉ－Ｖａｌｌｏｙｗｉｔｈｅｘｃｅｐｔｉｏｎａｌｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓ：ｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓａｎｄＡＰＴｉｎｖｅｓｔｉｇａｔｉｏｎｓ［Ｊ］ Ｃｏｍｐ Ｍａｔｅｒ Ｓｃｉ ，２０１９，１７０：１０９１６１［２１］ ＧａｏＹ，ＱｉａｏＬ，ＷｕＤ，ｅｔａｌ ＦｉｒｓｔｐｒｉｎｃｉｐｌｅｃａｌｃｕｌａｔｉｏｎｏｆｔｈｅｅｆｆｅｃｔｏｆＣｒ，ＴｉｃｏｎｔｅｎｔｏｎｔｈｅｐｒｏｐｅｒｔｉｅｓｏｆＶＭｏＮｂＴａＷＭｘ（Ｍ＝Ｃｒ，Ｔｉ）ｒｅｆｒａｃｔｏｒｙｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙ［Ｊ］ Ｖａｃｕｕｍ，２０２０，１７９：１０９４５９［２２］ ＺｈｏｕＸ，ＬｉＧ，ＧｕｏＦ，ｅｔａｌ ＦｉｒｓｔｐｒｉｎｃｉｐｌｅｓｉｎｖｅｓｔｉｇａｔｉｏｎｏｆｔｈｅＳｔｒｕｃｔｕｒａｌ，Ｍｅｃｈａｎｉｃａｌ，ＥｌｅｃｔｒｏｎｉｃＰｒｏｐ ｅｒｔｉｅｓａｎｄＤｅｂｙｅｔｅｍｐｅｒａｔｕｒｅｏｆＷ－Ｃｏａｌｌｏｙｓ［Ｊ］ ＦｕｓｉｏｎＥｎｇ Ｄｅｓ ，２０２１，１７０：１１２７１５ ［２３］ ＺｈａｎｇＪ，ＸｕＱ，ＨｕＹ，ｅｔａｌ ＩｎｔｅｒｆａｃｉａｌｂｏｎｄｉｎｇｍｅｃｈａｎｉｓｍａｎｄａｄｈｅｓｉｖｅｔｒａｎｓｆｅｒｏｆｂｒａｚｅｄｄｉａｍｏｎｄｗｉｔｈＮｉ－ｂａｓｅｄｆｉｌｌｅｒａｌｌｏｙ：Ｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓａｎｄｅｘｐｅｒｉｍｅｎｔａｌｐｅｒｓｐｅｃｔｉｖｅ［Ｊ］ Ｃａｒｂｏｎ，２０１９，１５３：１０４ ［２４］ ＢａｉＬ，ＨｕＹ，ＬｉａｎｇＸ，ｅｔａｌ ＴｉｔａｎｉｕｍａｌｌｏｙｉｎｇｅｎｈａｎｃｅｍｅｎｔｏｆｍｅｃｈａｎｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆＮｂＴａＭｏＷｒｅｆｒａｃｔｏｒｙｈｉｇｈ－Ｅｎｔｒｏｐｙａｌｌｏｙ：Ｆｉｒｓｔ－ｐｒｉｎｃｉｐｌｅｓａｎｄｅｘ ｐｅｒｉｍｅｎｔｓｐｅｒｓｐｅｃｔｉｖｅ［Ｊ］ Ｊ Ａｌｌｏｙ Ｃｏｍｐ ，２０２１，８５７：１５７５４２ ［２５］ ＳｈｏｕＨ，ＸｉｅＲ，ＰｅｎｇＭ，ｅｔａｌ Ｓｔａｂｉｌｉｔｙａｎｄｅｌｅｃｔｒｏｎ ｉｃｓｔｒｕｃｔｕｒｅｓｏｆｔｈｅＴｉＺｎｉｎｔｅｒｍｅｔａｌｌｉｃｃｏｍｐｏｕｎｄｓ：ＡＤＦＴｃａｌｃｕｌａｔｉｏｎ［Ｊ］ ＰｈｙｓｉｃａＢ，２０１９，５６０：４１［２６］ ＺｈｏｕＹ，ＬｉＹ，ＷａｎｇＷ，ｅｔａｌ Ｅｆｆｅｃｔｏｆｉｎｔｅｒｓｔｉｔｉａｌｎｉ ｔｒｏｇｅｎｉｎＦｅ１８Ｃｒ６Ｍｎ８ａｕｓｔｅｎｉｔｉｃａｌｌｏｙｓｆｒｏｍｄｅｎｓｉｔｙｆｕｎｃｔｉｏｎａｌｔｈｅｏｒｙ［Ｊ］ Ｊ Ｍａｇｎ Ｍａｇｎ Ｍａｔｅｒ ，２０１８，４６３：５７［２７］ ＴｏｎｇＹ，ＢａｉＬ，ＬｉａｎｇＸ，ｅｔａｌ ＩｎｆｌｕｅｎｃｅｏｆａｌｌｏｙｉｎｇｅｌｅｍｅｎｔｓｏｎｍｅｃｈａｎｉｃａｌａｎｄｅｌｅｃｔｒｏｎｉｃｐｒｏｐｅｒｔｉｅｓｏｆＮｂＭｏＴａＷＸ（Ｘ＝Ｃｒ，Ｚｒ，Ｖ，ＨｆａｎｄＲｅ）ｒｅｆｒａｃｔｏｒｙｈｉｇｈｅｎｔｒｏｐｙａｌｌｏｙｓ［Ｊ］ Ｉｎｔｅｒｍｅｔａｌｌｉｃｓ，２０２０，１２６：１０６９２８。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
平衡晶格常数及体弹模量的模拟计算（0K）
一、实验原理
1.1 平衡晶格常数
通过分子动力学模拟，在给定条件下，计算晶体结构最稳定，也就是体系能量最小时所对应的晶格间距，即为平衡晶格常数。
1.2 体弹模量
在弹性变形范围内,物体的体应力与相应体应变之比的绝对值称为体弹模量。
表达式为
式中，P 为体应力或物体受到的各向均匀的压强，dV/V为体积的相对变化。
对于立方晶胞，总能量可以表示为ε=ME，E 为单个原子的结合能，M 为单位晶胞内的原子数。
晶胞体积可以表示为V=a^3，那么压强P为
故体弹模量可以表示为
根据实验第一部分算出的平衡晶格常数，以及能量与晶格间距的函数关系，可以求得对应晶格类型的体积模量。
二、拟合方法
2.1 多项式拟合
使用下式对计算数据进行拟合，计算系数a、b、c。
平衡晶格常数即为-b/2a，二阶导数即为2a。
2.2 Birch-Murnaghan方程拟合
Birch-Murnaghan方程如下
通过这种方法可以直接拟合得出平衡晶格常数及体弹模量。
三、操作步骤
3.1 步骤及解释
$ cp -r share/1_lattice ～ &复制文件夹
$ cd 1_lattice &依次进入包含某一元素运行文件的文件夹中$ cd Cu (or Al, Si, Fe, Mg)
$ gedit ttice &编辑运行文件
$ lmp < ttice &使用lammps运行文件
$ A.i686 a0.cfg &使用ayomeye观察晶体结构
$ gnuplot plot.2nd.gnu (plot.bm.gnu) &拟合数据并作图
3.2 实际步骤（以Cu为例）
[user022@cluster ~]$ cp -r share/1_lattice/ ~
[user022@cluster ~]$ ls
1_lattice 2_point bin Desktop share
[user022@cluster ~]$ cd 1_lattice/
[user022@cluster 1_lattice]$ ls
Al Cu Fe Mg Si
[user022@cluster 1_lattice]$ cd Cu/
[user022@cluster Cu]$ ls
ttice jin_copper_lammps.setfl plot.2nd.gnu plot.bm.gnu [user022@cluster Cu]$ lmp < ttice
[user022@cluster Cu]$ A.i686 a0.cfg
[user022@cluster Cu]$ gnuplot plot.2nd.gnu
[user022@cluster Cu]$ gnuplot plot.bm.gnu
四、模拟数据
4.1 Mg
4.1.1 多项式拟合
Final set of parameters Asymptotic Standard Error ======================= ========================== a = 2.27117 +/- 0.002085 (0.09179%) b = -14.464 +/- 0.01328 (0.0918%) c = 21.4997 +/- 0.02114 (0.09833%) 4.1.2 Birch-Murnaghan 方程拟合
4.2 Al
HCP Lattice:
E_0 = -1.52868631023185 eV
a_0 = 3.18431542679217 Angstrom V_0 = 22.76161517799 Angstrom**3 B_0 = 36.0301849559129 GPa B_0'= -0.761464207502543
4.2.1 多项式拟合
Final set of parameters Asymptotic Standard Error ======================= ==========================
a = 2.20808 +/- 0.005565 (0.252%)
b = -17.8654 +/- 0.04503 (0.252%)
c = 32.7261 +/- 0.09108 (0.2783%)
4.2.2 Birch-Murnaghan方程拟合
FCC Lattice:
E_0 = -3.41065714040381 eV
a_0 = 4.04527130437683 Angstrom
V_0 = 16.5494273213047 Angstrom**3
B_0 = 77.7803912766984 GPa
B_0'= 7.11079881556912
4.3.1 多项式拟合
Final set of parameters Asymptotic Standard Error ======================= ========================== a = 1.93485 +/- 0.001141 (0.05897%) b = -21.0163 +/- 0.01239 (0.05896%) c = 52.7331 +/- 0.03365 (0.06381%) 4.3.2 Birch-Murnaghan 方程拟合
diamond Lattice:
E_0 = -4.33660000718975 eV
a_0 = 5.43095170306466 Angstrom V_0 = 20.0234005455525 Angstrom**3 B_0 = 101.425444944596 GPa B_0'= 2.85073370817905
4.4.1 多项式拟合
Final set of parameters Asymptotic Standard Error ======================= ========================== a = 7.12859 +/- 0.001914 (0.02686%) b = -40.7092 +/- 0.01093 (0.02686%) c = 53.9969 +/- 0.01561 (0.02891%) 4.4.2 Birch-Murnaghan 方程拟合
BCC Lattice:
E_0 = -4.1224351934112 eV
a_0 = 2.85532720281661 Angstrom V_0 = 11.6395892035167 Angstrom**3 B_0 = 177.84840115414 GPa B_0'= 1.41894857732246
4.5.1 多项式拟合
Final set of parameters Asymptotic Standard Error ======================= ========================== a = 3.48337 +/- 0.005239 (0.1504%) b = -25.1497 +/- 0.03783 (0.1504%) c = 41.9049 +/- 0.06829 (0.163%) 4.5.2 Birch-Murnaghan 方程拟合
FCC Lattice:
E_0 = -3.49001356518869 eV
a_0 = 3.6098665239098 Angstrom
V_0 = 11.7601656929196 Angstrom**3 B_0 = 137.631357881733 GPa B_0'= 4.22407949095044
表 1 多项式拟合结果及实验数值
表 2 bm拟合数值及实验数值
表 3 晶格常数模拟值与实验值的相对误差
表 4 体弹模量模拟值与实验值的相对误差
五、结论
对上述几种元素晶体晶格常数的模拟计算，两种拟合方法得出的结论与实验结果均符合得很好，两种方法产生的相对误差值也很接近。
但是计算体积模量时，计算公式放大了误差。
采用多项式拟合方法时，若采用4阶或5阶拟合，平衡晶格常数的计算效果会更好。