力学常数1.doc

VASP计算 -力学常数李智内蒙古科技大学2010年7月1日摘要本文主要介绍了用VASP对弹性模量、剪切模量、体积模量以及泊松比等力学常数计算，首先介绍了计算所需的相关基础知识，然后详细的阐述了理论的推导过程和对结果的处理方法，并介绍了VASP所需文件和生成的文件，最后提供了计算的一个例子和其程序流程图。目录一、基础知识1二、VASP计算时解析推导3三、VASP计算9四、有待继续研究的地方10五、参考文献10六、附录（一）程序流程图11七、附录（二）-一个例子，TaN12一、 基础知识12这部分主要介绍了进行VASP计算时所需要的概念的解释，其主要部分来自弹性力学，详细的介绍可阅读参考文献。1、 应力与应变a、 应力：某描述单位面积上一点的内力称为应力。单位：帕斯卡(Pa)，由于这个单位很小，通常使用MPa或GPa。反应的是材料在横截面A上的内力的合力F的强弱程度。b、 应变：描述一点处变形的程度的力学量是该点的应变。量纲为1。反应的是在外力作用下材料形变量s与其原长x之间的比值。c、 正应力()与正应变()：沿截面法线方向。d、 切应力()与切应变()：沿截面切向方向。2、 胡克定律（Hookes law）：在弹性限度内，物体的形变与引起形变的外力成正比。a、 表达式：其中F-物体受力，k-弹性系数，x-形变量b、 材料力学表达式：其中E-弹性模量（杨氏模量），G-切变模量，其量纲都是GPac、 广义胡克定律：d、 体积胡克定律：其中m-三个主应力的平均这值，-体积改变量，B-体积弹性模量经过推导计算可以得到体积模量与弹性模量和泊松比之间的关系3、 泊松比()：横向正应变与轴向正应变之比的绝对值。由于横向正应变与轴向正应变的变化是相反的，所以去掉绝对值要加负号。4、 Voigt标记：用向量表示对称矩阵5、 张量零阶张量就是标量，有30=1个量一阶张量是矢量，有31=3个量二阶张量两个相关的矢量，有32=9个量，如：应力张量，应变张量四阶张量，有34=81，如：弹性常数二、 VASP计算时解析推导345这部分主要对VASP计算过程的理论推导，并且介绍对计算结果的处理方法。这部分推导只限于结构为各向同性的正六面体，如需对其他结构进行计算这部分也列出了不同结构的弹性常数结构。1、 忽略：a、 忽略温度变化对体系总能的影响。b、 在小变形的条件下，忽略切应力()对正应变()的影响2、 对胡克定律变形 上一部分介绍的胡克定律的标准形式，将每个方向单独的应力应变关系及泊松比带入矩阵，并且就可得到如下矩阵形式：这样就得出了材料在x、y、z三个方向上的应变与各应力之间的的关系。由于VASP计算需要的是应变与能量的关系，所以需要将上式变成用应变来表示应变的形式，只需将矩阵求逆即可得到。用i统一表示正应变和剪切应变，用i统一表示正应力与剪切应力，用Cij表示其中的系数，这也是胡克定律的另一种标准形式。其中Cij就是我们要求的弹性常数。将上面两个矩阵进行系数对比，不难看出：其他都为零而且由于E、G、存在如下关系： 所以，实际上独立的变量只有C11，C12*胡克定律最终变形为：3、 力学常数的表达：a、 剪切模量G：b、 体积模量B：c、 弹性模量E：d、 泊松比：4、 力学常数的求解：1） 系统总能其中W-内能密度，V-系统体积2） 在应变较小的情况下，应变后体系的总能E(V，)按应变张量进可按泰勒级数展开为:对上面偏微分方程的求解后的到应变后体系的总能的变化量为：3） 应变后基矢与应变前的基矢之间的关系为: 其中为单位矩阵，为应变的张量矩阵。4） 对的选取：a、 求剪切模量G：求解剪切模量时，要求应变前与应变后的体积不变。每个晶胞的体积可以由基矢求得。如果要求体积不变，就是要求有如下关系：由此我们可以得出这样的一个应变的矩阵：用Voigt标记该矩阵：这就是在VASP计算剪切模量时的程序中所需的应变的形式，将其代入前面介绍的体系总能变化的式子，带入过程如下：代入得：由于，所以，代入上式，得：其中G就是我们要求的剪切模量，现在我们找出了体系能量变化与剪切模量和应变值之间的关系，当我们取多个不同的值，通过VASP计算，就可得到相应的体系能量变化的量，然后可以拟合出一条的二次曲线，得出的二次项系数A0（乘出结果后单位不是GPa，需将结果乘以160.2）剪切模量即为所求。b、 体积模量B：求体积模量的过程与前面求剪切模量的过程类似，我们在三个方向上都取相同的应变，就可以求得体积模量。用Voigt标记代入上面能量变化的式子，过程与上面的类似取多个不同的值，通过VASP计算后可得到相应的体系能量变化值，然后拟合出一条的二次曲线，得出的二次项系数A1体积模量即为所求。c、 弹性模量E与泊松比将前面数据代入即可得。5、 其他晶体类型的弹性常数1） 单斜晶系2） 正交晶系3） 三角晶系4） 六角晶系5） 立方晶系三、 VASP计算61、 需要准备的文件defvector.f OLDPOS KPOINTS POTCAR optimizea、 defvector.f这个文件经过编译后，是被optimize文件调用的子程序。它使用FORTRAN语言编写的。这个文件主要进行的内容是，对应变的Voigt标记的形式进行定义，变换基矢，以及生成VASP计算所需的POSCAR文件的数据。对于defvector.f需要进行编译，生成defvector.x文件才能被调用，在终端中输入：ifort o defvector.x defvector.fb、 OLDPOS 文件中是defvector.f所需的最初的原子排布结构，其形式与POSCAR的相似，只是在第一行多一个原子种类数。c、 KPOINTS 对倒空间K点的选择d、 POTCAR计算中，所需元素的赝势文件如果需要计算不止一个元素则需要把元素的POTCAR合并cat POTCAR_Ta POTCAR_N POTCARe、 optimize这个文件是计算时最重要的文件，是它对计算进行控制。它的内容包括计算时所需的不同的应变值，以及驰豫和计算能量时的两个INCAR文件。它的功能是，(1)对不同的应变值进行循环的VASP计算，(2)调用子程序，并生成计算所需的POSCAR文件，(3)生成INCAR文件，并进行切换，(4)控制VASP软件开始计算，(5)提取计算结果数据。VASP计算对两个INCAR文件有以下的参数调整要求：驰豫时：IBRION = 2 ISIF = 2能量计算时：ISMEAR = -52、 计算在终端中输入bash optimize或./ optimize，回车后，如果没有错误VASP就开始计算直至得到最终结果。3、 VASP计算后得到的有关文件，以及对数据的处理SUMMARY文件记录了系统能量E和相应的应变，将能量与当=0时的能量想减，得到E，然后拟合出一条的二次曲线。二次曲线的二次项系数就是所需的AiOUTCAR文件中记录了很多信息，其中volume of cell是指晶胞的体积V0。代入下式即可得到力学常数。 四、 有待继续研究的地方1、 P3页最下边，得出的胡克定律中，而张亮论文中得出的数值两者并不相等，原因有待继续研究。2、 不同类型的晶格所对应的应变的形式没有求出，有待确定。3、 对于非各项同性的立方晶体，弹性常数Cij与各力学常数的关系没有得到，有待解决。五、 参考文献1 刘鸿文.材料力学(M).第四版.北京:高等教育出版社,2004.1:238-239 2 程昌钧.弹性力学(M）.第一版.1995.11：87-1063 单耀祖. 材料力学(M).第二版. 北京:高等教育出版社,2004.8:256-257 4 侯柱锋.采用VASP如何计算晶体的弹性常数Cij5 朱晓玲. 超硬材料PtN_2和ReB_2力学性质的第一性原理计算(D).四川师范大学。2209.46 张亮. Ti-Si-N类超硬薄膜的结构和成形机理的第一原理计算(D).内蒙古科技大学.2009.3六、 附录（一）程序流程图 Optimize文件 defvector.f文件结束生成SUMMARYN是否还有新的应变用CONTCAR生成计算能量的POSCAR生成驰豫的INCAR进行VASP计算Y提取能量值进行VASP计算生成计算能量的INCAR用fort.3生成驰豫的POSCAR生成fort.3文件存放处理后的数据进行基矢的应变处理调用应变值调用OLDPOS并提取基矢调用子程序defvector.x定义不同的应变值开始./optimize对不同循环应变七、附录（二）-一个例子，TaN1、 defvector.fC% !注释行C this simple program to get the primitive vectors afterC $delta$ strain, in order to calculate the independentC elastic constants of solids.C usage: C! Please first prepare the undeformed POSCAR inC OLDPOSC defvector.xC type defvector.x create new POSCAR in file fort.3C% program defvector !程序名为defvector real*8 privect,strvect,delta,strten,strain,pos, alat !定义8位实变量 dimension privect(3,3),strvect(3,3),strten(3,3),strain(6) dimension pos(50,3) !定义相应规格的数组 character*10 bravlat, title, direct !定义长度为10的字符串 integer i,j,k,ntype, natomi, nn !定义整型变量 dimension natomi(10)C% Read the undeformed primitive vector and atomic postionC% open(7,file=OLDPOS)!打来OLDPOS文件，标明单元号位7供后面调用，并默认为顺序打开C% In first line of OLDPOS, please add the numberC% of the type of atoms after the title read(7,*) title, ntype !读取单元号为7的第一行赋给两个变量 read(7,*) alat !读取第二行并赋值 do i=1,3 !循环开始到第一个enddo结束 read(7,*) (privect(i,j),j=1,3) !读取并赋值给数组 write(*,*) (privect(i,j),j=1,3) !屏幕显示 enddo !循环结束 read(7,*) (natomi(i),i=1,ntype) !读取每种元素的原子个数 nn=0 !清空变量 do i =1, ntype nn=nn+natomi(i) !计算原子总数 enddo read(7,*) direct !字符串赋值 do i=1, nn read(7,*) (pos(i,j),j=1,3) !读取原子总数行数据并赋值 enddoC% Read the amti of strain % read(*,*) delta !读取外部输入并赋值C% Define the strain % strain(1)=delta !为应变的数组赋值，确定应变类型 strain(2)=delta strain(3)=0.0 strain(4)=0.0 strain(5)=0.0 strain(6)=0.0C% Define the strain tensor % strten(1,1)=strain(1)+1.0 !将Voigt标记换回矩阵，进行加上单位矩阵 strten(1,2)=0.5*strain(6) strten(1,3)=0.5*strain(5) strten(2,1)=0.5*strain(6) strten(2,2)=strain(2)+1.0 strten(2,3)=0.5*strain(4) strten(3,1)=0.5*strain(5) strten(3,2)=0.5*strain(4) strten(3,3)=strain(3)+1.0C% Transform the primitive vector to the new vector underCstrain%C strvect(i,j)=privect(i,j)*(I+strten(i,j) !注释行 do k=1,3 !矩阵相乘求的形变后的基矢 do i=1,3 strvect(i,k)=0.0 do j=1,3 strvect(i,k)=strvect(i,k)+privect(i,j)*strten(j,k) enddo enddo enddoC% Write the new vector under strain% do i=1,3 write(*,100)(strvect(i,j),j=1,3) !显示基矢，100为格式说明符标号 enddo 100 FORMAT(3F20.15)! 100:上面的格式说明，3:重复3次，F:实数型，20:字段宽度，15:小数部分宽度C% Create the POSCAR for total energy calculationC%5 write(3,(A10) title !在单元号为3的位置输出变量，A:字符型，10:宽度 write(3,(F15.10) alat !输出基矢系数 do i=1,3 write(3,100)(strvect(i,j),j=1,3) !输出基矢 enddo write(3,(10I4)(natomi(i), i=1,ntype) !输出各元素原子数，I:十进制整数型 write(3,(A6) Direct !输出字符 do i=1, nn write(3,100) (pos(i,j),j=1,3) !输出原子坐标 enddo !注：单元号除*、0、5、6有特定意义外，其他如果没有特别定义会以fort.n文件输出，n为单元号。C% end !程序结束2、 OLDPOSTaN 24.258 2.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 2.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 2.00000000000000004 4Direct0.00 0.00 0.000.00 0.25 0.250.25 0.25 0.000.25 0.00 0.250.25 0.25 0.250.25 0.00 0.000.00 0.25 0.000.00 0.00 0.253、 KPOINTStan0Monkhorst-Pack7 7 70. 0. 0.4、 optimize#!/bin/bashfor i in -0.018 -0.015 -0.012 -0.009 -0.006 -0.003 0.00 0.003 0.006 0.009 0.012 0.015 0.018 #应变所需的数值do #开始对不同的应变值进行循环echo $i | ./defvector.x #显示应变值并赋给defvector.x，开始运行defvector.xcp fort.3 POSCAR #将defvector.x产生的fort.3变成POSCAR备用#cat INCAR INCAR SUMMARY #将应变与能量赋给SUMMARYdon