当前位置:
文档之家› 采用VASP如何计算晶体的弹性常数
采用VASP如何计算晶体的弹性常数
SYSTEM = AlN ENCUT = 400 ISTART = 0 ICHARG = 2
5
3
ISMEAR = 0; SIGMA = 0.2 NSW = 60; IBRION = 2 EDIFF = 1E-5 EDIFFG = -1E-2 ISIF = 2 POTIM = 0.2 PREC = Accurate LWAVE = .FALSE.
4
3
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 enddo C%%%%%%%% Write the new vector under strain%%%%%%%%%%%% do i=1,3 write(*,100)(strvect(i,j),j=1,3) enddo 100 format(3f20.15) C%%%%%%%%% Create the POSCAR for total energy calculation %%%%%%%%%%%%%%5 write(3,’(A10)’) title 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) write(3,’(A6)’) Direct do i=1, nn write(3,100) (pos(i,j),j=1,3) enddo C%%%%%%% end
ij
¦
Î
[1, 2]
σ
ùþǑñ Ǒñ
ǫ=
e1 1 2 e6 1 2 e5
1 2 e6
e2
1 2 e5 1 2 e4
1 2 e4
e3
.
(1)
σi =
¦
1 V
∂E (V, ǫj ) ∂ǫi
,
ǫ=0
(2)
Cij =
1 V
∂ 2 E (V, ǫk ) ∂ǫi ∂ǫj 1
,
ǫ=0
(3)
2
ÆÇ
ù
•
ª
 ”Define the strain”Í
ù
³
¿
defvector.f
”Define the strain”
e = (δ, δ, 0, 0, 0, 0)
» ñ
Ç
Ö ù
C11 + C12
Ó
Ä
strain(i)
defvect.f
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
¿
¢ ìõ
E (V0 , 0) − → R
Æ
ÓÍE(V, ǫ)Í
6
ǫ V0 2
6
Í Î
ÑǑñ
(4) − →′ R
E (V, {ǫi }) = E (V0 , 0) + V0
σi ei +
i=1
Cij ei ej + . . .
i,j =1
ð Æ Ñ ³ Ǒñ
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
(9)
ñ
æ
∆E = C44 δ 2 V0 e = (δ, δ, δ, 0, 0, 0) C11 C12 C13 C33 ∆E = (C11 + C12 + 2C13 + C33 /2)δ 2 V0
(10)
ñ
(11)
Äæ Í ³ Æ ¦
ù
∆E ∼ δ
ü è ¯ δ ½ á
2
∆E ∼ δ
à
¾ Ç Æ
¿defvector.f ¦ ntype ùþ Ǒ10 ¦ î” write(3,’(10I4)’) (natomi(i), i=1,ntype)” Á Æ natomi(10) ”10I4” ¿defvector.f èª ù defvector.f g77 -o defector.x defector.f) defvector.x • Ç VASPÇ KPOINTS POTCAR ½¦ èò Ç INCAR.relax ËÆ ñ
¸©¿Ç
SYSTEM = AlN ENCUT = 400 ISTART = 0 ICHARG = 2 ISMEAR = -5 EDIFF = 1E-5 PREC = Accurate LWAVE = .FALSE.
Þ
´¨ ½ÓÍ ½
INCAR.static
Ǒ
− →′ − → R = R • (I + ǫ)
î è Æ
Æ ½ Æ
§2
Æ ð Æ æ
Æ ¦ì
Æ Æ ¦
a1
Ǒñ
Í
1 2a 1 2a
³
0 0 . c (6)
− → R = a2 = a3
√ 3 a 2 √ 3 − 2 a
0
0
a
c
e = (δ, δ, 0, 0, 0, 0)
defvector.f
ð
¸
3
POSCAR
ºÍ”Define the
3
strain”
ËÆ ñ
ùþ
Á¿ Ë °½
î
defvect.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 >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)
VASP
Ç
2006
Òñ õ 8 Ò
21
Æ ¦
Cij
(zfhou@)
Æ é Æ î
ý Ǒ ¼« ÄVASPÇ Æ ¦ C ËÆ æ » ììì
ij
§1 §2 §3
¦ Î Æ Æ ¦ ÆÇ ù
õ
[1, 2] 1
Í
§1
³
2 3
¦ ª Æ © ǫ Ǒ ¿ ¢ ìõ Æ ËÍ ¢ ¿ ³ § ù ¦ C ð ³ VoigtÁÈñ xx → 1, yy → 2, zz → 3, yz → 4, xz → 5 xy → 6 ǫù á þǑñ
ÓÍ V ð
0
¦ì Æ ¸©
(5)
ª ü ù ǫ = e = (e , e , e , e , e , e ) Ç ¿è ÓÍ Þ(△E = E(V, ǫ) − E(V, 0)) ¾ ÓÍ Þ ø á Æ ¸ ¦ ù¦ Æ Ǒ ³ ¦ ð´ù ³ ¦ ǑñC , C , C , C C
I
1 2 3 4 5 6 11 12 13 33 44
Ç ñ
C11 + C12 [3]
ñ
(7)
æ
33
∆E = (C11 + C12 )δ 2 V0 e = (0, 0, 0, 0, 0, δ ) C11 − C12 ∆E 1 = (C11 − C12 )δ 2 V0 4 e = (0, 0, δ, 0, 0, 0)