VASP固定基矢优化结构方法

本文章为原创,版权归作者刘锦程所有,文章转载请先取得作者的同意,非常欢迎转发文章链接!严禁以任何方式挪用本文内容,用于以盈利为目的各种活动。

在对称性允许的条件下,VASP的晶胞优化(ISIF = 3)是允许在9个自由度上自由弛豫的,如果想要固定其中几个自由度需要用重新编译过的vasp。

重写constr_cell_relax.F,直接覆盖./src/的文件即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
!-----------------------------------------------------------------------
!
! At present, VASP does not allow to relax the cellshape selectively
! i.e. for instance only cell relaxation in x direction.
! To be more precisse, this behaviour can not be achived via the INCAR
! or POSCAR file.
! However, it is possible to set selected components of the stress tensor
! to zero.
! The most conveninent position to do this is the routines
! CONSTR_CELL_RELAX (constraint cell relaxation).
! FCELL contains the forces on the basis vectors.
! These forces are used to modify the basis vectors according
! to the following equations:
!
! A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
! DO J=1,3
! DO I=1,3
! DO K=1,3
! A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
! ENDDO
! ENDDO
! ENDDO
! where A holds the basis vectors (in cartesian coordinates).
!
!-----------------------------------------------------------------------

SUBROUTINE CONSTR_CELL_RELAX(FCELL)
USE prec
REAL(q) FCELL(3,3)

! just one simple example
! relaxation in x directions only
! SAVE=FCELL(1,1)
! FCELL=0 ! F90 style: set the whole array to zero
! FCELL(1,1)=SAVE
! relaxation in z direction only
! SAVE=FCELL(3,3)
! FCELL=0 ! F90 style: set the whole array to zero
! FCELL(3,3)=SAVE

LOGICAL FILFLG
INTEGER ICELL(3,3)
INQUIRE(FILE='OPTCELL',EXIST=FILFLG)
IF (FILFLG) THEN
OPEN(67,FILE='OPTCELL',FORM='FORMATTED',STATUS='OLD')
DO J=1,3
READ(67,"(3I1)") (ICELL(I,J),I=1,3)
ENDDO
CLOSE(67)
DO J=1,3
DO I=1,3
IF (ICELL(I,J)==0) FCELL(I,J)=0.0
ENDDO
ENDDO
ENDIF

RETURN
END SUBROUTINE

比如优化单层MoS2,POSCAR如下:

1
2
3
4
5
6
7
8
9
10
11
MoS2
1.0
3.1659998894 0.0000000000 0.0000000000
-1.5829999447 2.7418363326 0.0000000000
0.0000000000 0.0000000000 18.4099998474
S Mo
2 1
Direct
0.000000000 0.000000000 0.413899988
0.000000000 0.000000000 0.586099982
0.666666687 0.333333343 0.500000000

VASP优化的时候,用OPTCELL就限制真空层就行了。OPTCELL如下:

1
2
3
100
110
000

ISIF = 3 优化之后CONTCAR如下:

1
2
3
4
5
6
7
8
9
10
11
MoS2                                    
1.00000000000000
3.1822561751580456 0.0000000000000000 0.0000000000000000
-1.5911280875790228 2.7559146890376600 0.0000000000000000
0.0000000000000000 0.0000000000000000 18.4099998474000017
S Mo
2 1
Direct
-0.0000000000000000 -0.0000000000000000 0.4150537162890733
-0.0000000000000000 0.0000000000000000 0.5849462537109242
0.6666666870000029 0.3333333429999996 0.5000000000000000

由此可见 a = b = 3.1659998894,变成了 a = b = 3.1822561751580456。而 c = 18.4099998474没有变。

如果不固定基矢优化,得到的CONTCAR如下:

1
2
3
4
5
6
7
8
9
10
11
MoS2                                    
1.00000000000000
3.1822639795633445 0.0000000000000000 0.0000000000000000
-1.5911319897816723 2.7559214478509104 0.0000000000000000
0.0000000000000000 -0.0000000000000000 18.2494484388264482
S Mo
2 1
Direct
-0.0000000000000000 -0.0000000000000000 0.4143158624764718
-0.0000000000000000 -0.0000000000000000 0.5856841075235258
0.6666666870000029 0.3333333429999996 0.5000000000000000

由此可见,不仅ab 优化了,c也被优化了。对于二维材料我们并不希望在真空层的方向优化。

再举个例子 Theta-Al2O3的优化,POSCAR如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
al4 o6
1.0
6.1415114403 0.0000000000 0.0000000000
-5.4373541295 2.8554058978 0.0000000000
-1.3346698477 0.3291362294 5.5013695377
Al O
4 6
Direct
0.909878987 0.090120998 0.204389006
0.090121022 0.909879037 0.795611026
0.657850062 0.342150010 0.316776984
0.342150002 0.657850009 0.683223004
0.839793974 0.160205982 0.891047027
0.160205983 0.839793908 0.108952995
0.504978997 0.495020970 0.742557983
0.495020993 0.504979010 0.257441984
0.173605024 0.826394978 0.567412014
0.826394930 0.173604994 0.432588018

如果不用OPTCELL进行晶胞优化,在保证整体的空间群对称性条件下,晶胞的整体会发生旋转,最后得到的CONTCAR如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
al4 o6                                  
1.00000000000000
6.0834761338287953 0.0022219821590104 0.0007877232894261
-5.3849398048848336 2.8303904510334665 -0.0007877232893012
-1.3212365137134128 0.3258235023583434 5.4612696048070495
Al O
4 6
Direct
0.9093415914296588 0.0906583935703435 0.2039066722878985
0.0906584175703384 0.9093416414296629 0.7960933597121042
0.6582648301748026 0.3417352418251964 0.3167705605521709
0.3417352338251957 0.6582647771748018 0.6832294274478280
0.8376321892727857 0.1623677667272178 0.8911128886483282
0.1623677677272144 0.8376321232727802 0.1088871333516737
0.5047059001915972 0.4952940668084034 0.7439206754098269
0.4952940898084018 0.5047059131915947 0.2560792915901740
0.1729719163094664 0.8270280856905338 0.5656625638453888
0.8270280376905298 0.1729718863094710 0.4343374681546067

可见在y和z分量为0的 a也发生了变化,这并不意味着晶胞的对称性发生了变化,只是晶胞整体发生了转动。如果不想让晶胞整体转动只需要固定3个上对角矩阵元或者3个下对角矩阵元,比如:

1
2
3
100
110
111

最后得到的CONTCAR如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
al4 o6                                  
1.00000000000000
6.0833317459767864 0.0000000000000000 0.0000000000000000
-5.3846770309770191 2.8328089483277861 0.0000000000000000
-1.3211733148948146 0.3263555668976587 5.4614651686107925
Al O
4 6
Direct
0.9093549922684255 0.0906438462883539 0.2039382615101327
0.0906450167315719 0.9093561887116525 0.7960617704898700
0.6582203631760277 0.3417794562164434 0.3168160759831926
0.3417797008239705 0.6582205627835545 0.6831839120168064
0.8377447190014768 0.1622536363223279 0.8911091000542825
0.1622552379985234 0.8377462536776701 0.1088909219457195
0.5047198237549158 0.4952780523628080 0.7438557622658220
0.4952801662450831 0.5047219276371903 0.2561442047341789
0.1729935617079242 0.8270063674354625 0.5657116781331480
0.8270063922920718 0.1729936045645422 0.4342883538668472

可见,除了OPTCELL中为0的位置,其他的6个自由度都被优化了。

修改constr_cell_relax.F实现固定基矢优化结构的原理

FCELL contains the forces on the basis vectors. These forces are used to modify the basis vectors according to the following equations:

FCELL 是个3 * 3的矩阵,包含了晶格矢量的受力,当ISIF = 3时,晶格矢量会按照如下的方法更新:

1
2
3
4
5
6
7
8
A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style 
DO J=1,3
DO I=1,3
DO K=1,3
A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
ENDDO
ENDDO
ENDDO

当通过OPTCELL文件把晶格矢量在某个方向的受力改成0以后,

IF (ICELL(I,J)==0) FCELL(I,J)=0.0

在结构更新的时候,该方向就不会改变。

参考资料:

固定晶格部分无缺大师兄贡献:https://zhuanlan.zhihu.com/p/34750103

http://bbs.keinsci.com/thread-9494-1-1.html


转载请注明来源,欢迎对文章中的引用来源进行考证,欢迎指出任何有错误或不够清晰的表达。