VASP计算力学常数1.pdf

.VASP 计算 -力学常数 摘要 本文主要介绍了用 VASP对弹性模量、剪切模量、体积模量以及泊松比等力学常数计算，首先介绍了计算所需的相关基础知识，然后详细的阐述了理论的推导过程和对结果的处理方法，并介绍了 VASP所需文件和生成的文件，最后提供了计算的一个例子和其程序流程图。目录 一、基础知识 .1 二、VASP 计算时解析推导.3 三、VASP 计算.9 四、有待继续研究的地方.10 五、参考文献.10 六、附录（一）程序流程图.11 七、附录（二）-一个例子，TaN.12 一、基础知识12 这部分主要介绍了进行 VASP 计算时所需要的概念的解释，其主要部分来自弹性力学，详细的介绍可阅读参考文献。1、应力与应变 a、应力：某描述单位面积上一点的内力称为应力。单位：帕斯卡(Pa)，由于这个单位很小，通常使用 MPa或 GPa。0limAFA 反应的是材料在横截面A 上的内力的合力F 的强弱程度。b、应变：描述一点处变形的程度的力学量是该点的应变。量纲为 1。.0limxsx 反应的是在外力作用下材料形变量s与其原长x之间的比值。c、正应力()与正应变()：沿截面法线方向。d、切应力()与切应变()：沿截面切向方向。2、胡克定律（Hookes law）：在弹性限度内，物体的形变与引起形变的外力成正比。a、表达式：Fkx 其中 F-物体受力，k-弹性系数，x-形变量 b、材料力学表达式：;EG 其中 E-弹性模量（杨氏模量），G-切变模量，其量纲都是 GPa c、广义胡克定律：111121314151612212223242526233132333435363141424344454612515253545556236162636465663ssssssssssssssssssssssssssssssssssss d、体积胡克定律：mB 其中m-三个主应力的平均这值，-体积改变量，B-体积弹性模量 经过推导计算可以得到体积模量与弹性模量和泊松比之间的关系 3(1 2)EB 3、泊松比()：横向正应变与轴向正应变之比的绝对值。由于横向正应变与轴向正应变的变化是相反的，所以去掉绝对值要加负号。4、Voigt标记：用向量表示对称矩阵 111652211624123456221154322eeeeeeeeeeeeeee 5、张量 零阶张量就是标量，有 30=1 个量 一阶张量是矢量，有 31=3 个量.二阶张量两个相关的矢量，有 32=9 个量，如：应力张量，应变张量 四阶张量，有 34=81，如：弹性常数 二、VASP 计算时解析推导345 这部分主要对 VASP 计算过程的理论推导，并且介绍对计算结果的处理方法。这部分推导只限于结构为各向同性的正六面体，如需对其他结构进行计算这部分也列出了不同结构的弹性常数结构。1、忽略：a、忽略温度变化对体系总能的影响。b、在小变形的条件下，忽略切应力()对正应变()的影响 2、对胡克定律变形 上一部分介绍的胡克定律的标准形式，将每个方向单独的应力应变关系及泊松比带入矩阵，并且就可得到如下矩阵形式：111122133111122133000000000000000000000000EEEEEEEEEGGG 这样就得出了材料在 x、y、z 三个方向上的应变与各应力之间的的关系。由于 VASP 计算需要的是应变与能量的关系，所以需要将上式变成用应变来表示应变的形式，只需将矩阵求逆即可得到。(1)(1 2)(1)(1 2)(1)(1 2)(1)1(1)2(1 2)(1)(1 2)(1)(1 2)(1)(1)3(1 2)(1)(1 2)(1)(1 2)(1)123000000000000000000000000EEEEEEEEEGGG 123123 用i统一表示正应变和剪切应变，用i统一表示正应力与剪切应力，用 Cij表示其中的系数，这也是胡克定律的另一种标准形式。.111121314151612212223242526233132333435363441424344454645515253545556566162636465666cccccccccccccccccccccccccccccccccccc 其中 Cij就是我们要求的弹性常数。将上面两个矩阵进行系数对比，不难看出：112233121321233132445566(1)(1 2)(1)(1 2)(1)EcccEcccccccccG 其他都为零 而且由于 E、G、存在如下关系：2(1)EG 所以，实际上独立的变量只有 C11，C12 *胡克定律最终变形为：111212111211122212121133111124421111255211112662000000000000()000000()000000()ccccccccccccccc 3、力学常数的表达：a、剪切模量 G：11121()2(1)2EGGcc b、体积模量 B：11121(2)3(1 2)3EBBcc c、弹性模量 E：.2(1)123(1 2)93EGEGEEGBGE 93BGEBG d、泊松比：3262BGBG 4、力学常数的求解：1）系统总能 EWV 其中 W-内能密度，V-系统体积 123456()iiiWf 2(,)1kijijijE VWCV 2）在应变较小的情况下，应变后体系的总能 E(V，)按应变张量进可按泰勒级数展开为:660001,1(,)(,0)2ii iij ijii jVE VE VVC 对上面偏微分方程的求解后的到应变后体系的总能的变化量为：660112ij ijijVEC 3）应变后基矢与应变前的基矢之间的关系为:()I 其中I为单位矩阵，为应变的张量矩阵。4）对的选取：a、求剪切模量 G：求解剪切模量时，要求应变前与应变后的体积不变。每个晶胞的体积可以由基矢求得。1112366()V 如果要求体积不变，就是要求有如下关系：.()I 由此我们可以得出这样的一个应变的矩阵：2000000(1+)1 用 Voigt标记该矩阵：2(1+)1000 这就是在 VASP 计算剪切模量时的程序中所需的应变的形式，将其代入前面介绍的体系总能变化的式子，带入过程如下：21232456(2);=(1+)1;(1)0 代入 660112ij ijijVEC 得：1 1 1 11 2 121 3 1 302 1 2 12 2 222 3 233 1 3 13 2 323 3 331 11 21 221 21 11 221 21 2221120000(2)(1)(2)(1)(2)(2)(1)(1)(2)ECCCVCCCCCCCCCCCCCCC 222111242(2)(1)(1)(2)(2)1 11 1 4(1)(1)CC 由于0.05，所以22;11 ，代入上式，得：.1112022111202663()6ECCVECCGV 其中 G 就是我们要求的剪切模量，现在我们找出了体系能量变化与剪切模量和应变值之间的关系，当我们取多个不同的值，通过 VASP 计算，就可得到相应的体系能量变化的量，然后可以拟合出一条E 的二次曲线，得出的二次项系数 A0（乘出结果后单位不是 GPa，需将结果乘以 160.2）006AGV 剪切模量即为所求。b、体积模量 B：求体积模量的过程与前面求剪切模量的过程类似，我们在三个方向上都取相同的应变，就可以求得体积模量。000000 用 Voigt标记 000 代入上面能量变化的式子，过程与上面的类似 22392211120(2)EccBV 取多个不同的值，通过 VASP 计算后可得到相应的体系能量变化值，然后拟合出一条E 的二次曲线，得出的二次项系数 A1 1029ABV 体积模量即为所求。c、弹性模量 E 与泊松比 93BGEBG 3262BGBG 将前面数据代入即可得。.5、其他晶体类型的弹性常数 1）单斜晶系 11121316212223263132333644455455616263660000000000000000cccccccccccccccccccc 2）正交晶系 111213212223313233445566000000000000000000000000cccccccccccc 3）三角晶系 111213141521221314151313331414444515154414111245140000000000002ccccccccccccccccccccccccc 4）六角晶系 1112131211131313334444111122000000000000000000000000()ccccccccccccc 5）立方晶系 111212121112121211444444000000000000000000000000cccccccccccc .三、VASP 计算6 1、需要准备的文件 defvector.f OLDPOS KPOINTS POTCAR optimize a、defvector.f 这个文件经过编译后，是被 optimize 文件调用的子程序。它使用 FORTRAN语言编写的。这个文件主要进行的内容是，对应变的 Voigt标记的形式进行定义，变换基矢，以及生成 VASP 计算所需的 POSCAR 文件的数据。对于 defvector.f 需要进行编译，生成 defvector.x 文件才能被调用，在终端中输入：ifort o defvector.x defvector.f b、OLDPOS 文件中是defvector.f所需的最初的原子排布结构，其形式与POSCAR 的相似，只是在第一行多一个原子种类数。c、KPOINTS 对倒空间 K 点的选择 d、POTCAR 计算中，所需元素的赝势文件 如果需要计算不止一个元素则需要把元素的 POTCAR 合并 cat POTCAR_Ta POTCAR_N POTCAR e、optimize 这个文件是计算时最重要的文件，是它对计算进行控制。它的内容包括计算时所需的不同的应变值，以及驰豫和计算能量时的两个 INCAR 文件。它的功能是，(1)对不同的应变值进行循环的 VASP 计算，(2)调用子程序，并生成计算所需的 POSCAR 文件，(3)生成 INCAR 文件，并进行切换，(4)控制 VASP 软件开始计算，(5)提取计算结果数据。VASP 计算对两个 INCAR 文件有以下的参数调整要求：驰豫时：IBRION=2 ISIF=2 能量计算时：ISMEAR=-5 2、计算 在终端中输入 bash optimize 或./optimize，回车后，如果没有错误 VASP 就开始计算直至得到最终结果。3、VASP 计算后得到的有关文件，以及对数据的处理 SUMMARY文件记录了系统能量 E 和相应的应变，将能量与当=0 时的能量想减，得到E，然后拟合出一条E 的二次曲线。二次曲线的二次项系数就是所需的 Ai OUTCAR 文件中记录了很多信息，其中 volume of cell是指晶胞的体积 V0。代入下式即可得到力学常数。006AGV 1029ABV 93BGEBG 3262BGBG.四、有待继续研究的地方 1、P3 页最下边，得出的胡克定律中14411122()ccc，而张亮论文中得出的数值两者并不相等，原因有待继续研究。2、不同类型的晶格所对应的应变的形式没有求出，有待确定。3、对于非各项同性的立方晶体，弹性常数 Cij与各力学常数的关系没有得到，有待解决。五、参考文献 1 刘鸿文.材料力学(M).第四版.北京:高等教育出版社,2004.1:238-239 2 程昌钧.弹性力学(M）.第一版.1995.11：87-106 3 单耀祖.材料力学(M).第二版.北京:高等教育出版社,2004.8:256-257 4 侯柱锋.采用 VASP 如何计算晶体的弹性常数 Cij 5 朱晓玲.超硬材料 PtN_2 和 ReB_2 力学性质的第一性原理计算(D).四川师范大学。2209.4 6 张亮.Ti-Si-N 类超硬薄膜的结构和成形机理的第一原理计算(D).内蒙古科技大学.2009.3 .六、附录（一）程序流程图 Optimize文件 defvector.f 文件 开始./optimize 定义不同的应变值 对不同循环应变 调用子程序 defvector.x 调用 OLDPOS 并提取基矢 调用应变值 进行基矢的应变处理 生成 fort.3 文件存放处理后的数据 用 fort.3 生成驰豫的 POSCAR 生成驰豫的 INCAR 进行 VASP 计算 用 CONTCAR 生成计算能量的 POSCAR 生成计算能量的 INCAR 进行 VASP 计算 提取能量值 是否还有新的应变 生成 SUMMARY 结束 Y N.七、附录（二）-一个例子，TaN 1、defvector.f C%!注释行 C this simple program to get the primitive vectors after C$delta$strain,in order to calculate the independent C elastic constants of solids.C usage:C Please first prepare the undeformed POSCAR in C OLDPOS C defvector.x C type defvector.x create new POSCAR in file fort.3 C%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 postion C%open(7,file=OLDPOS)!打来OLDPOS文件，标明单元号位 7 供后面调用，并默认为顺序打开 C%In first line of OLDPOS,please add the number C%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)!读取原子总数行数据并赋值 enddo C%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.0 C%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.0 C%Transform the primitive vector to the new vector under Cstrain%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 enddo C%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 calculation C%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、OLDPOS TaN 2 4.258 2.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 2.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 2.0000000000000000 4 4 Direct 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.00 0.25 0.00 0.25 0.25 0.25 0.25 0.25 0.00 0.00 0.00 0.25 0.00 0.00 0.00 0.25 3、KPOINTS tan 0 Monkhorst-Pack 7 7 7 0.0.0.4、optimize#!/bin/bash for 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.x cp fort.3 POSCAR#将defvector.x产生的fort.3 变成POSCAR备用#cat INCAR INCAR SUMMARY#将应变与能量赋给SUMMARY done#完成循环 cat SUMMARY#生成SUMMARY 5、SUMMARY-0.03 1 F=-.70020775E+03 E0=-.70020775E+03 d E=0.000000E+00-0.02 1 F=-.70301675E+03 E0=-.70301675E+03 d E=0.000000E+00-0.01 1 F=-.70490289E+03 E0=-.70490289E+03 d E=0.000000E+00 0.00 1 F=-.70556982E+03 E0=-.70556982E+03 d E=0.000000E+00 0.01 1 F=-.70483837E+03 E0=-.70483837E+03 d E=0.000000E+00 0.02 1 F=-.70242837E+03 E0=-.70242837E+03 d E=0.000000E+00 0.03 1 F=-.69823921E+03 E0=-.69823921E+03 d E=0.000000E+00 6、A0的拟合.应变量 总能 E/ev 能量差 E/ev -0.03-700.20777 5.36204-0.02-703.01670 2.55311-0.01-704.90289 0.66692 0.00-705.56981 0.00000 0.01-704.83837 0.73144 0.02-702.42838 3.14143 0.03-698.23923 7.33058 y=7055.7x2+25.524x+0.0042-1.00.01.02.03.04.05.06.07.08.0-0.04-0.03-0.02-0.01 0.000.010.020.030.04应变量能量差E/ev 二次曲线方程得到的二次项系数 A0=7055.7。7、计算结果 C11/GPa C12/GPa C44/GPa G/GPa B/GPa PW91 709.3 157.1 75.8 276.1 341.2