VASP固定基矢优化结构方法

啦啦黑还黑

本文章为原创，转发自http://blog.wangruixing.cn/2019/05/05/constr/，版权归作者Jin-Cheng Liu所有，文章转载请先取得作者的同意，非常欢迎转发文章链接）！

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

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

1. !-----------------------------------------------------------------------
   
2. !
   
3. ! At present, VASP does not allow to relax the cellshape selectively
   
4. ! i.e. for instance only cell relaxation in x direction.
   
5. ! To be more precisse, this behaviour can not be achived via the INCAR
   
6. ! or POSCAR file.
   
7. ! However, it is possible to set selected components of the stress tensor
   
8. ! to zero.
   
9. ! The most conveninent position to do this is the routines
   
10. ! CONSTR_CELL_RELAX  (constraint cell relaxation).
    
11. ! FCELL contains the forces on the basis vectors.
    
12. ! These forces are used to modify the basis vectors according
    
13. ! to the following equations:
    
14. !
    
15. !      A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
    
16. !      DO J=1,3
    
17. !      DO I=1,3
    
18. !      DO K=1,3
    
19. !        A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
    
20. !      ENDDO
    
21. !      ENDDO
    
22. !      ENDDO
    
23. ! where A holds the basis vectors (in cartesian coordinates).
    
24. !
    
25. !-----------------------------------------------------------------------
    
26. ​
    
27.       SUBROUTINE CONSTR_CELL_RELAX(FCELL)
    
28.       USE prec
    
29.       REAL(q) FCELL(3,3)
    
30. ​
    
31. !     just one simple example
    
32. !     relaxation in x directions only
    
33. !      SAVE=FCELL(1,1)
    
34. !      FCELL=0   ! F90 style: set the whole array to zero
    
35. !      FCELL(1,1)=SAVE
    
36. !     relaxation in z direction only
    
37. !      SAVE=FCELL(3,3)
    
38. !      FCELL=0   ! F90 style: set the whole array to zero
    
39. !      FCELL(3,3)=SAVE
    
40. ​
    
41.       LOGICAL FILFLG
    
42.       INTEGER ICELL(3,3)
    
43.       INQUIRE(FILE='OPTCELL',EXIST=FILFLG)
    
44.       IF (FILFLG) THEN
    
45.          OPEN(67,FILE='OPTCELL',FORM='FORMATTED',STATUS='OLD')
    
46.          DO J=1,3
    
47.             READ(67,"(3I1)") (ICELL(I,J),I=1,3)
    
48.          ENDDO
    
49.          CLOSE(67)
    
50.          DO J=1,3
    
51.          DO I=1,3
    
52.             IF (ICELL(I,J)==0) FCELL(I,J)=0.0
    
53.          ENDDO
    
54.          ENDDO
    
55.       ENDIF
    
56.    
    
57.       RETURN
    
58.       END SUBROUTINE
    
59. 
    

复制代码

比如优化单层MoS2，POSCAR如下：

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

复制代码

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

1. 100
   
2. 110
   
3. 000

复制代码

ISIF = 3 优化之后CONTCAR如下：

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

复制代码

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

如果不固定基矢优化，得到的CONTCAR如下：

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

复制代码

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

再举个例子 Theta-Al2O3的优化，POSCAR如下：

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

复制代码

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

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

复制代码

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

1. 100
   
2. 110
   
3. 111

复制代码

最后得到的CONTCAR如下：

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. A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
   
2. DO J=1,3
   
3. DO I=1,3
   
4. DO K=1,3
   
5.    A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
   
6. ENDDO
   
7. ENDDO
   
8. ENDDO

复制代码

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

1. 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

飞翔的猪

老师您好，我想请问一下：
1、这个脚本对于晶系有没有要求？

2、我是三维和二维材料都有涉及，如果在src文件夹下同时存在原来的constr_cell_relax.F和重新编译的constr_cell_relax.F.fixed C，在调用vasp计算的时候该如何处理呢？

3、“VASP优化的时候，用OPTCELL就限制真空层就行了”这句话不是很理解？是要单独写一个文件还是？

dingniu2

你好，关于修改constr_cell_relax.F文件重新编译vasp计算二维材料有两个问题请教：
1. 这种方法是否只适用于正交的晶胞，对非正交是否适用？
2. 这种方法在用vasp计算二维材料方面是否已被大家广泛采用？
非常感谢！

huangtf

采用楼主的方法，优化晶胞，结果γ角增大了一点。请问各位老师和同学，编译和计算方面有问题吗？这个结果能用于扩胞和后续研究吗？  谢谢！

luchaolu

dingniu2 发表于 2020-6-25 17:55
你好，关于修改constr_cell_relax.F文件重新编译vasp计算二维材料有两个问题请教：
1. 这种方法是否只适用 ...

正交适用的，请问有关OPTCELL怎么设置，以及怎么操作使其在VASP中生效，还是只是简单的出现在计算目录下就可以了

physics_xw

huangtf 发表于 2022-1-2 17:14
采用楼主的方法，优化晶胞，结果γ角增大了一点。请问各位老师和同学，编译和计算方面有问题吗？这个结果能 ...

你试试110，110，000看看，我之前也遇到过这种情况，应该是冗余约束条件导致的。另一种可能性是VASP计算时对称性的识别出了问题，和ISYM以及SYMPREC相关参数有关，可以借助VASPKIT在计算之前检查一下。
其实我认为你这个结果用起来应该也没有太大问题，如果实在不放心可以按我上面提到的检查一下。我之前出现过γ偏离120度，也出现过a和b不相等的情况，按上面的方法处理过后基本就没问题了。

huangtf

physics_xw 发表于 2022-1-7 16:48
你试试110，110，000看看，我之前也遇到过这种情况，应该是冗余约束条件导致的。另一种可能性是VASP计算 ...

老师，您好！根据您的讲解，我注释掉ISYM = 0，采用默认值（ISYM = 2），其余保持不变；重新优化晶胞，结果γ角仍为120度。谢谢老师！

huangtf

请问各位老师，使用固定基矢的版本（vasp_z）优化好二维（单层）晶胞。对于后面的单点、band、DOS、phonon和AIMD等计算任务，都需要改用没有固定基矢的标准版本（vsap_std）吗？  谢谢！

少广

重新编译过的vasp ，这里仅仅是把按作者说的文件覆盖就行了，还是需要重新make？

shameless

huangtf 发表于 2022-4-13 16:47
请问各位老师，使用固定基矢的版本（vasp_z）优化好二维（单层）晶胞。对于后面的单点、band、DOS、phonon ...

是的