在应变较小的情况下，应变后体系的总能 按应变张量 进可按泰勒级数展开为：
(110)
其中
应变前的基矢
是应变前体系的总能，
之间的关系为：
(111)
是应变前原胞的体积。另外，应变后基矢
与
其中 为单位矩阵。
因此选取特定的应变 ，计算出在一组不同幅度时应
变前后体系总能的变化( )，再根据总能的变化－应变幅度对应的一组数据点，进行二次函数拟合得到二次项系数。即可得到晶体的某个弹性常数或弹性常数的组合。对不同的晶系的晶体，因为对称性的关系，它独立的弹性常数是确定的。比如对六角晶系的晶体，它独立的弹性常数为：
,
,
和
。
,
六角晶体的弹性应变能与应变关系
对六角晶系的晶体，其原胞基矢可以取为：
(112)
其中 和 是晶体的晶格常数。
可以施加这样的应变 来计算
[#!WangSQ:Abiec:2003!#]：
(113)
施加应变
来得到
：
(114)
对于
，施加应变
来得到：
(115)
对于
，施加应变
来得到：
(116)
最后，施加应变
得到
、
、
和
的组合：
(117)
由此可见，通过施加五个特定的应变，选取一系列幅度 的应变，得到 数据点，再分别按上面五个关系式对相应的 进行拟合得到二次项系数，最后联立方程得到求出六角晶系晶体的独立弹性常数。
具体计算步骤
在计算时，有几点要特别注意的：a)，原胞内原子在应变是否驰豫；b)，k点网格大小是否足够，因为在应变后，原胞的对称性会发生变化，即使同样的k点网格，在简约布里渊区产生的k点数目是不同的。因此，k点网格大小要取得足够，以保证弹性常数计算的精确性；c)，应变幅度 要取得适中，如果太小的话，得到的应变能（应变前后体系的变化）很小，在计算弹性常数时，会引起计算误差。
下面以计算六角AlN（纤锌矿结构）为例，在计算过程考虑了应变后原子的位置需要驰豫，具体步骤如下：
?
先对六角AlN体材料的晶格参数（晶格常数和原子位置）优化好，这样得到未应变时的POSCAR，并把它拷贝成下面defvector.f需要用到的一个的输入文件OLDPOS，并对OLDPOS做一个处理。因为OLDPOS在格式上有特殊要求：
a
、在OLDPOS的第一行，在title的字符串之后，至少空一格再加上OLDPOS中原子的种类数目，比如AlN中有两类原子，title写为AlN，那么就OLDPOS的第一行就是\。 b
、OLDPOS的格式与POSCAR的类似，但是它最好是分数坐标来写出原子的位置。 以六角AlN为例，这个OLDPOS的内容如下：
towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0 AlN 2 3.11553
1.000000 0.000000 0.000000 -0.500000 0.866025 0.000000 0.000000 0.000000 1.605000 2 2 Direct
0.00000000 0.00000000 0.00000000 0.333333333 0.666666667 0.50000000 0.00000000 0.00000000 0.381483673 0.333333333 0.666666667 0.881483673
towidth 0.8pt height 0.25cm depth 0cmwidth 1.025height 0.8pt depth 0cmwidth 0.8pt height 0.25cm depth 0cm
?
对特定的应变，在下面的defvector.f中\部分，把特定的应变通过给strain(i)矩阵赋值。比如对上面提到的
用来计算
，
那么就对defvect.f中的\改写成如下的形式：
towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0 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
towidth 0.8pt height 0.25cm depth 0cmwidth 1.025height 0.8pt depth 0cmwidth 0.8pt height 0.25cm depth 0cm
其中整个defvector.f是用来得到某个应变后，新的POSCAR。应变的类型按\the strain\的部分来定义，而 应变的幅度需在程序编译后，运行编译得到的模块时输入。defvect.f的内容如下：
towidth 0.8pt height 0cm depth 0.25cmwidth 1.025height 0cm depth 0.8ptwidth 0.8pt height 0cm depth 0.25cm 1.0
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 OLDPOS C >defvector.x C >type defvector.x > create new POSCAR in file fort.3
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
program defvector
real*8 privect,strvect,delta,strten,strain,pos, alat dimension
privect(3,3),strvect(3,3),strten(3,3),strain(6) dimension pos(50,3)
character*10 bravlat, title, direct integer i,j,k,ntype, natomi, nn dimension natomi(10)
C%%%%%%%%% Read the undeformed primitive vector and atomic postion %%%%%%%
open(7,file='OLDPOS')
C%% In first line of OLDPOS, please add the number C%% of the type of atoms after the title
read(7,*) title, ntype read(7,*) alat do i=1,3
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)=0.0 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 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 strain%%%%%
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
百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库VASP几个计算实例(3)在线全文阅读。
相关推荐：