VASP计算力学常数1Word下载.docx

VASP计算力学常数1Word下载.docx

《VASP计算力学常数1Word下载.docx》由会员分享，可在线阅读，更多相关《VASP计算力学常数1Word下载.docx（18页珍藏版）》请在冰豆网上搜索。

沿截面切向方向。

2、胡克定律（Hooke'

slaw）：

在弹性限度内，物体的形变与引起形变的外力成正比。

a、表达式：

其中F---物体受力，k---弹性系数，x---形变量

b、材料力学表达式：

其中E---弹性模量（杨氏模量），G---切变模量，其量纲都是GPa

c、广义胡克定律：

d、体积胡克定律：

其中σm---三个主应力的平均这值，θ---体积改变量，B---体积弹性模量

经过推导计算可以得到体积模量与弹性模量和泊松比之间的关系

3、泊松比（μ）：

横向正应变与轴向正应变之比的绝对值。

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

4、Voigt标记：

用向量表示对称矩阵

5、张量

零阶张量就是标量，有30=1个量

一阶张量是矢量，有31=3个量

二阶张量两个相关的矢量，有32=9个量，如：

应力张量，应变张量

四阶张量，有34=81，如：

弹性常数

二、VASP计算时解析推导[3][4][5]

这部分主要对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）

剪切模量即为所求。

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

用Voigt标记

代入上面能量变化的式子，过程与上面的类似

取多个不同的δ值，通过VASP计算后可得到相应的体系能量变化值，然后拟合出一条的二次曲线，得出的二次项系数A1

体积模量即为所求。

c、弹性模量E与泊松比ε

将前面数据代入即可得。

5、其他晶体类型的弹性常数

1）单斜晶系

2）正交晶系

3）三角晶系

4）六角晶系

5）立方晶系

三、VASP计算[6]

1、需要准备的文件

defvector.fOLDPOSKPOINTSPOTCARoptimize

a、defvector.f

这个文件经过编译后，是被optimize文件调用的子程序。

它使用FORTRAN语言编写的。

这个文件主要进行的内容是，对应变的Voigt标记的形式进行定义，变换基矢，以及生成VASP计算所需的POSCAR文件的数据。

对于defvector.f需要进行编译，生成defvector.x文件才能被调用，在终端中输入：

ifort–odefvector.xdefvector.f

b、OLDPOS

文件中是defvector.f所需的最初的原子排布结构，其形式与POSCAR的相似，只是在第一行多一个原子种类数。

c、KPOINTS

对倒空间K点的选择

d、POTCAR

计算中，所需元素的赝势文件

如果需要计算不止一个元素则需要把元素的POTCAR合并

catPOTCAR_TaPOTCAR_N>

POTCAR

e、optimize

这个文件是计算时最重要的文件，是它对计算进行控制。

它的内容包括计算时所需的不同的应变值，以及驰豫和计算能量时的两个INCAR文件。

它的功能是，

（1）对不同的应变值进行循环的VASP计算，

（2）调用子程序，并生成计算所需的POSCAR文件，（3）生成INCAR文件，并进行切换，（4）控制VASP软件开始计算，（5）提取计算结果数据。

VASP计算对两个INCAR文件有以下的参数调整要求：

驰豫时：

IBRION=2ISIF=2

能量计算时：

ISMEAR=-5

2、计算

在终端中输入bashoptimize或./optimize，回车后，如果没有错误VASP就开始计算直至得到最终结果。

3、VASP计算后得到的有关文件，以及对数据的处理

SUMMARY文件记录了系统能量E和相应的应变δ，将能量与当δ=0时的能量想减，得到ΔE，然后拟合出一条的二次曲线。

二次曲线的二次项系数就是所需的Ai

OUTCAR文件中记录了很多信息，其中volumeofcell是指晶胞的体积V0。

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

四、有待继续研究的地方

1、P3页最下边，得出的胡克定律中，而张亮论文中得出的数值两者并不相等，原因有待继续研究。

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文件

七、附录

（二）------一个例子，TaN

1、defvector.f

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!

注释行

C>

thissimpleprogramtogettheprimitivevectorsafter

C$\delta$strain,inordertocalculatetheindependent

Celasticconstantsofsolids.

Cusage:

CPleasefirstpreparetheundeformedPOSCARin

COLDPOS

defvector.x

typedefvector.x>

createnewPOSCARinfilefort.3

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

programdefvector!

程序名为defvector

real*8privect,strvect,delta,strten,strain,pos,alat!

定义8位实变量

dimensionprivect（3,3）,strvect（3,3）,strten（3,3）,strain（6）

dimensionpos（50,3）!

定义相应规格的数组

character*10bravlat,title,direct!

定义长度为10的字符串

integeri,j,k,ntype,natomi,nn!

定义整型变量

dimensionnatomi（10）

C%%%%%%%%%Readtheundeformedprimitivevectorandatomicpostion

C%%%%%%%

open（7,file='

OLDPOS'

）

!

打来OLDPOS文件，标明单元号位7供后面调用，并默认为顺序打开

C%%InfirstlineofOLDPOS,pleaseaddthenumber

C%%ofthetypeofatomsafterthetitle

read（7,*）title,ntype!

读取单元号为7的第一行赋给两个变量

read（7,*）alat!

读取第二行并赋值

doi=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!

清空变量

doi=1,ntype

nn=nn+natomi（i）!

计算原子总数

enddo

read（7,*）direct!

字符串赋值

doi=1,nn

re