VASP计算-力学常数1

VASP 计算
----------力学常数
摘要
本文主要介绍了用VASP 对弹性模量、剪切模量、体积模量以及泊松比等力学常数计算，首先介绍了计算所需的相关基础知识，然后详细的阐述了理论的推导过程和对结果的处理方法，并介绍了VASP 所需文件和生成的文件，最后提供了计算的一个例子和其程序流程图。

目录
一、 基础知识 .................................................................................................................... 1 二、 VASP 计算时解析推导 .............................................................................................. 3 三、 VASP 计算 .................................................................................................................. 9 四、 有待继续研究的地方 .............................................................................................. 10 五、 参考文献 .................................................................................................................. 10 六、 附录（一）程序流程图 .......................................................................................... 11 七、
附录（二）------一个例子，TaN (12)
一、 基础知识[1][2]
这部分主要介绍了进行VASP 计算时所需要的概念的解释，其主要部分来自弹性力学，详细的介绍可阅读参考文献。

1、 应力与应变
a 、 应力：某描述单位面积上一点的内力称为应力。

单位：帕斯卡(Pa)，由于
这个单位很小，通常使用MPa 或GPa 。

0lim
A F A ∆→∆∆ 反应的是材料在横截面△A 上的内力的合力△F 的强弱程度。

b 、 应变：描述一点处变形的程度的力学量是该点的应变。

量纲为1。

0lim
x s x ∆→∆∆ 反应的是在外力作用下材料形变量△s 与其原长△x 之间的比值。

c 、 正应力(σ)与正应变(ε)：沿截面法线方向。

d 、 切应力(τ)与切应变(γ)：沿截面切向方向。

2、 胡克定律（Hooke's law ）：在弹性限度内，物体的形变与引起形变的外力成正比。

a 、 表达式：
F kx =-
其中F---物体受力，k---弹性系数，x---形变量 b 、 材料力学表达式：
;
E G σετγ==
其中E---弹性模量（杨氏模量），G---切变模量，其量纲都是GPa c 、 广义胡克定律：
1111213141516122122
232425262331323334353631414243444546125152535455562361
62
63
64
65
663s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s εσεσεσγτγτγτ⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
d 、 体积胡克定律：
m B σθ=
其中σm ---三个主应力的平均这值，θ---体积改变量，B---体积弹性模量
经过推导计算可以得到体积模量与弹性模量和泊松比之间的关系
3(12)E B μ=
-
3、 泊松比(μ)：横向正应变与轴向正应变之比的绝对值。

由于横向正应变与轴向正应变
的变化是相反的，所以去掉绝对值要加负号。

εεμεε''=
=-
4、 Voigt 标记：用向量表示对称矩阵
()
111
6522
1
16
2
412345
622115
4
322
e e e e e e e e e e e e e e e εε⎛⎫ ⎪=⇔= ⎪ ⎪⎝⎭
5、 张量
零阶张量就是标量，有30=1个量 一阶张量是矢量，有31=3个量
二阶张量两个相关的矢量，有32=9个量，如：应力张量，应变张量 四阶张量，有34=81，如：弹性常数
二、 VASP 计算时解析推导[3][4][5]
这部分主要对VASP 计算过程的理论推导，并且介绍对计算结果的处理方法。

这部分推导只限于结构为各向同性的正六面体，如需对其他结构进行计算这部分也列出了不同结构的弹性常数结构。

1、 忽略：
a 、 忽略温度变化对体系总能的影响。

b 、 在小变形的条件下，忽略切应力(τ)对正应变(ε)的影响 2、 对胡克定律变形
上一部分介绍的胡克定律的标准形式，将每个方向单独的应力应变关系及泊松比带入矩阵，并且就可得到如下矩阵形式：
1
111221*********
330000000
000000
000
0000
00
E E
E E E
E
E
E E
G
G
G μ
μμμμ
μεσεσεσγτγτγτ⎡⎤--⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎣⎦⎣⎦⎣⎦ 这样就得出了材料在x 、y 、z 三个方向上的应变与各应力之间的的关系。

由于VASP 计算需要的是应变与能量的关系，所以需要将上式变成用应变来表示应变的形式，只需将矩阵求逆即可得到。

(1)
(12)(1)(12)(1)(12)(1)1(1)2(12)(1)(12)(1)(12)(1)(1)3(12)(1)(12)(1)
(12)(1)
1230
000000000000000000
00E E E E E E E E E G G G μμ
μ
μμμμμμμμ
μ
μμμμμμμμ
μ
μμμμμμσσστττ--+-+-+--+-+-+--+-+-+⎡⎤⎡⎤⎢⎢⎥⎢⎢⎥⎢⎢⎥⎢=⎢⎥⎢⎢⎥⎢⎢⎥⎢⎢⎥⎢⎢⎥⎣⎦⎣123123εεεγγγ⎡⎤⎥⎢⎥
⎥⎢⎥
⎥
⎢⎥
⎥⎢⎥⎥⎢⎥⎥⎢⎥
⎥⎢⎥
⎥⎢⎥⎣⎦⎦ 用εi 统一表示正应变和剪切应变，用σi 统一表示正应力与剪切应力，用C ij 表
示其中的系数，这也是胡克定律的另一种标准形式。

111121314151612212223
242526233132333435363441424344454645515253545556566162
63
64
65
666c c c c c c c
c c c c c c c c c c c c c c c c c c c c c c c c c c c c c σεσεσεσεσεσε⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
其中C ij 就是我们要求的弹性常数。

将上面两个矩阵进行系数对比，不难看出：
112233121321233132445566(1)(12)(1)
(12)(1)
E c c c E c c c c c c c c c G
μμμμ
μμ-===
-+======-+===
其他都为零
而且由于E 、G 、μ存在如下关系：
2(1)E G μ=
+
所以，实际上独立的变量只有C 11，C 12
*胡克定律最终变形为：
11121211121112221212113311112442
11112552
1
1112662000
000000
000()0000
00()0
()c c c c c c c c c c c c c c c σεσεσεσεσεσε⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎢
⎥⎢⎥⎣⎦⎣⎦⎣⎦ 3、 力学常数的表达：
a 、 剪切模量G ：
11121
()
2(1)2E G G c c μ=
⇒=-+
b 、 体积模量B ：
11121
(2)
3(12)3E B B c c μ=
⇒=+-
c 、 弹性模量E ：
2(1)1
23(12)93E G E G E EG
B G E
μμμ=
+=-==--
93BG E B G =
+
d 、 泊松比μ：
3262B G B G μ-=
+
4、 力学常数的求解：
1） 系统总能
E WV =
其中W---内能密度，V---系统体积
123456()
i i i
W
f σεεεεεεε∂=
=∂
2(,)1k ij i j i j
E V W C V εεεεε∂∂==
∂∂∂∂
2） 在应变较小的情况下，应变后体系的总能E(V ，ε)按应变张量进可按泰勒级数展
开为:
6
6
001
,1
(,{})(,{0})2i i i ij i
j
i i j V
E V E V V C εσεεε
===+++⋅⋅⋅
∑∑
对上面偏微分方程的求解后的到应变后体系的总能的变化量为：
66
11
2
ij i
j
i j V
E C εε
==∆=∑∑
3） 应变后基矢与应变前的基矢之间的关系为: ()I ααε'
=+ 其中I 为单位矩
阵，ε为应变的张量矩阵。

4） 对ε的选取：
a 、 求剪切模量G ：
求解剪切模量时，要求应变前与应变后的体积不变。

每个晶胞的体积可以由基矢求得。

11
12366()V αααα
=⨯⋅=
如果要求体积不变，就是要求有如下关系：
()I ααεα
'=+=
由此我们可以得出这样的一个应变的矩阵：
2000000(1+)1δεδδ-⎛⎫
⎪
= ⎪
⎪-⎝⎭
用Voigt 标记该矩阵：
()
2(1+)1000εδδδ-=-
这就是在VASP 计算剪切模量时的程序中所需的应变的形式，将其代入
前面介绍的体系总能变化的式子，带入过程如下：
21232
456(2);;=(1+)1;
(1)0
δδ
εδεδεδδεεε-+==-=-+===
代入
66
11
2
ij i
j
i j V E C εε
==∆=∑∑
得：
1111121213130
2121222223233131323233331112122
1211122
121222
1120000
(2)[](1)(2)[](1)
(2)(2)[][](1)(1)(2)[E
C C C V C C C C C C C C C C C C C C C εεεεεεεεεεεεεεεεεεδδ
δδδδδδδδ
δδδδδδδδδδδδδδδ∆=+++++++++++++=++-
+++++-++++-+-++++-22
2111242
(2)][](1)(1)
(2)(2)
[11][114](1)(1)C C δδδδδδδδδδδδδ+-++++=++++-++
由于0.05δ<，所以22;11δδ+≈+≈，代入上式，得：
11120
2211120
2663()6E
C C V E
C C G V δδδδδδ∆=-∆=-=
其中G 就是我们要求的剪切模量，现在我们找出了体系能量变化与剪切模量和应变值之间的关系，当我们取多个不同的δ值，通过VASP 计算，就可得到相应的体系能量变化的量，然后可以拟合出一条E δ∆≈的二次曲线，得出的二次项系数A 0（乘出结果后单位不是GPa ，需将结果乘以160.2）
06A G V =
剪切模量即为所求。

b 、 体积模量B ：
求体积模量的过程与前面求剪切模量的过程类似，我们在三个方向上都取相同的应变，就可以求得体积模量。

000000δεδδ⎛⎫
⎪
= ⎪
⎪⎝⎭
用Voigt 标记
()
000εδδδ=
代入上面能量变化的式子，过程与上面的类似
223
92211120
(2)E c c B V δδ∆=+=
取多个不同的δ值，通过VASP 计算后可得到相应的体系能量变化值，然后拟合出一条E δ∆≈的二次曲线，得出的二次项系数A 1
1029A B V =
体积模量即为所求。

c 、 弹性模量E 与泊松比ε
93BG E B G =
+ 3262B G B G μ-=
+ 将前面数据代入即可得。

5、 其他晶体类型的弹性常数
1） 单斜晶系
11121316212223
263132333644455455616263
66000000000000000
c c c c c
c c c c c c c c c c c c c c c ⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎣
⎦
2） 正交晶系
1112
13212223
31323344556600000000000000000000
00
c c c c
c c c c c c c c ⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎣⎦
3） 三角晶系
11
12131415
21221314151313331414444515154414111245
14
0000000000
2c c c c c c c c c c c c c c c c c c c c c c c c c ⎡⎤⎢⎥--⎢⎥⎢⎥
⎢⎥--⎢⎥⎢⎥
-⎢⎥
-⎢⎥-⎢⎥⎣
⎦ 4） 六角晶系
11121312111313133344441
11122000
000000
0000000
000
()c c c c c c c c c c c c c ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎢⎥⎢⎥-⎢⎥⎣⎦
5） 立方晶系
11
121212111212121144444400000000000000000000
c c c c c c c c c c c c ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
三、 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 δ∆≈的二次曲线。

二次曲线的二次项系数就是所需的A i
OUTCAR 文件中记录了很多信息，其中volume of cell 是指晶胞的体积V 0。

代入下式即可得到力学常数。

006A G V =
1
029A
B V =
93BG E B G =
+ 3262B G
B G μ-=
+
四、有待继续研究的地方
1、P3页最下边，得出的胡克定律中
1
441112
2
()
c c c
=-
，而张亮论文中得出的数值两者
并不相等，原因有待继续研究。

2、不同类型的晶格所对应的应变的形式没有求出，有待确定。

3、对于非各项同性的立方晶体，弹性常数C ij与各力学常数的关系没有得到，有待解决。

五、参考文献
[1]刘鸿文.材料力学Ⅰ(M).第四版.北京:高等教育出版社,2004.1:238-239
[2] 程昌钧.弹性力学(M）.第一版.1995.11：87-106
[3] 单耀祖. 材料力学Ⅰ(M).第二版. 北京:高等教育出版社,2004.8:256-257
[4] 侯柱锋.采用VASP如何计算晶体的弹性常数C ij
[5] 朱晓玲. 超硬材料PtN_2和ReB_2力学性质的第一性原理计算(D).四川师范大学。

2209.4
[6]张亮. Ti-Si-N类超硬薄膜的结构和成形机理的第一原理计算(D).内蒙古科技大
学.2009.3
六、附录（一）程序流程图
七、附录（二）------一个例子，TaN
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
tan
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，并复制!!中间的文字SYSTEM = TaN
ENCUT = 400
ISTART = 0
ICHARG = 2
ISMEAR = 0; SIGMA = 0.2
NSW = 60; IBRION = 2
EDIFF = 1E-5
EDIFFG = -1E-2
ISIF = 2
POTIM = 0.2
PREC = Accurate
LWAVE = .FALSE.
LCHARG =.FALSE.
!
echo "delta = $i" ; vasp #显示应变值，开始驰豫VASP计算cp CONTCAR pos.$i
cp CONTCAR POSCAR #将驰豫完的数据生成计算能量的POSCAR cat > INCAR <<! #生成计算能量的INCAR SYSTEM = AlN
ENCUT = 400
ISTART = 0
ICHARG = 2
ISMEAR = -5
EDIFF = 1E-5
PREC = Accurate
LWAVE = .FALSE.
LCHARG =.FALSE.
!
echo "delta = $i "; vasp
E=`grep "TOTEN" out.$i | tail -1 | awk' {printf "%12.6f \n", $5 }'` #提取能量值echo $i $E >>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
7、计算结果
C11/GPa C12/GPa C44/GPa G/GPa B/GPa PW91 709.3 157.1 75.8 276.1 341.2。