VASP 固定基矢优化结构方法

chempeng

一直以来在本站学习颇多，感谢 sob 与大佬们的答疑解难。此贴为经验整理总结帖，希望对萌新们（我算一个~）有所帮助，如有疏漏还请各位批评斧正

一般对于三维材料的结构晶格参数优化有两种方法（教程参见 Learn VASP The Hard Way）：

1. Birch-Murnaghan 状态方程拟合 http://www.bigbrosci.com/newsitem/278011633

2. VASP 计算中调节 ISIF 参数 http://www.bigbrosci.com/newsitem/278019758

而上述方法对于优化表面，二维或者一维体系不适用，优化此类材料时，须固定一个或两个晶胞方向不优化。一般可通过晶格参数扫描再拟合函数实现，引侯柱峰老师对该类问题解答回帖。

如果二维参数里有一个参数（也就是a,b相等）,那通过固定 c，改 a 或 b 的值，对每一个 a 值下的结构进行 ISIF=2 的计算，得到总能，然后根据一系列的 总能-晶格常数 a 的值进行一元二次函数的拟合，得到最小值的位置。
如果是含两个常数（即 a,b 不等），通过地，固定 c，改变a和b地值，对每一组 a 和 b 值下的结构进行 ISIF=2 的计算，得到总能，然后根据一系列的 总能-晶格常数 a 和 b 的值，进行二元二次函数的拟合得到最小值的位置。

这样的方法普适性很高，但却是繁琐。于是有了在特定条件下，重新编译 vasp 使得该优化简化的方法，即当材料为正交晶系或单斜晶系时，可通过修改 constr_cell_relax.F 编译软件实现固定基矢进行结构优化。

对于 VASP 5.4.1 之后的版本，软件的安装编译已经方便了很多，本文以 intel 编译器 + mkl + impi 编译 VASP 5.4.1 为例进行介绍。VASP 的安装编译无需 root 权限，普通用户在家目录即可安装。

1. 查看编译器，mkl 等安装及环境变量设置正确

1. export LD_LIBRARY_PATH=/share/intel/composer_xe_2015.6.233/mkl/lib/intel64/:$LD_LIBRARY_PATH
   
2. export LD_LIBRARY_PATH=/share/intel/composer_xe_2015.6.233/compiler/lib/intel64:$LD_LIBRARY_PATH
   
3. export PATH=/share/intel/composer_xe_2015.6.233/bin/intel64:$PATH
   
4. export PATH=/share/intel/impi/5.0.3.049/intel64/bin:$PATH

复制代码

2. 修改 constr_cell_relax.F

以固定 c 方向为例，constr_cell_relax.F （src 目录下）修改后如下，显示部分修改内容

1.       SUBROUTINE CONSTR_CELL_RELAX(FCELL)
   
2.       USE prec
   
3.       REAL(q) FCELL(3,3)
   
4.       DO I=1,3
   
5.         FCELL(3,I)=0
   
6.         FCELL(I,3)=0
   
7.       ENDDO
   
8. !     just one simple example
   
9. !     relaxation in x directions only
   
10. !      SAVE=FCELL(1,1)
    
11. !      FCELL=0   ! F90 style: set the whole array to zero
    
12. !      FCELL(1,1)=SAVE
    
13. !     relaxation in z direction only
    
14. !      SAVE=FCELL(3,3)
    
15. !      FCELL=0   ! F90 style: set the whole array to zero
    
16. !      FCELL(3,3)=SAVE
    
17. 
    
18.       RETURN
    
19.       END SUBROUTINE

复制代码

对于其他方向晶格基矢的修改同理：对于a方向基矢，将 FCELL(3,I) 和 FCELL(I,3) 分别改为 FCELL(1,I) 和 FCELL(I,1) ；对于b方向基矢，则分别改为 FCELL(2,I) 和 FCELL(I,2) ；固定两个基矢则应该同时将两个方向对应的矩阵元设置为0。或参考 VASP并行可执行软件包，可对晶胞参数进行部分优化(自行搜索) 可实现添加参数文件实现不同基矢方向固定。（未做测试）

3. 编译

建议重新解压源码编译软件

1. [yourname@localhost ~]$ tar -zxvf vasp.5.4.1.05Feb16.tar.gz
   
2. [yourname@localhost ~]$ tar -zxvf vasp.5.lib_.tar_2.gz
   
3. [yourname@localhost ~]$ cd vasp.5.4.1
   
4. [yourname@localhost ~]$ cp arch/makefile.include.linux_intel makefile.include
   
5. [yourname@localhost ~]$ make std

复制代码

普通用户在家目录安装编译前在 makefile.include 文件中添加 $MKLROOT，位置如下

1. MKLROOT    =/share/intel/composer_xe_2015.6.233/mkl
   
2. MKL_PATH   = $(MKLROOT)/lib/intel64

复制代码

编译结束后，进入 bin 目录，可将 vasp_std 重命名为 vasp_fixZ ，将软件移入系统 bin 目录下或将该目录添加至环境变量。

1. [yourname@localhost ~]$ echo 'export PATH=/home/yourname/vasp.5.4.1/bin:$PATH' >> ~/.bashrc
   
2. [yourname@localhost ~]$ source ~/.bashrc

复制代码

完成编译后，在 INCAR 文件中设置 ISIF=3，便可实现固定 Z 轴的结构优化。

修改 constr_cell_relax.F 资料来源科大 @用 在遇见大师兄群内的答疑，另附有 VASP 基矢矩阵的更新流程

汪杰

chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

请问一下，最上面的100是代表x方向的基矢吗？还是竖着来看，110才是x方向的基矢呢？如果只是固定z方向，为什么不是下面这样：
100
010
000
能否给一下具体原理，感激不尽！

汪杰

啦啦黑还黑 发表于 2018-11-13 22:11
楼主的方法只能限制Z方向，可以简单写一个3*3的矩阵文件OPTCELL限制任意的晶格矢量。
      INTEGER ICELL ...

请问一下，这样改之后，OPTCELL应该怎么设置呢？能不能给一个具体的例子

helpcal

请问直接isif=3固定c轴优化和调整ab轴画E-a曲线得到的晶格常数相差约0.002，这时候应该取哪一个值更准确。

另外E-a曲线方法是不是取能量最低点然后对其进行isif-2的优化？这个方法里我对近邻几个a也优化原子位置，结果发现别的a驰豫完能量更低一点点。

fjy

shameless 发表于 2022-6-18 11:19
你好 固定两个轴的OPTCELL 应该怎么写

把要固定的两个轴都写成0就好了

shameless

luchaolu 发表于 2022-1-4 14:15
OPTCELL固定c轴这样对吗？不是应该为100 010 000吗

你好 固定两个轴的OPTCELL 应该怎么写

luchaolu

Xpeinny 发表于 2019-3-26 12:14
请问一下，您测试过三斜体系吗？其中gamma角度和c的大小会受影响发生改变吗？我测试出来c轴大小不变，但是g ...

OPTCELL固定c轴这样对吗？不是应该为100 010 000吗

luchaolu

chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

请问这个OPTCELL文件怎么设置呢？？怎么能够让脚本将设置的OPTCELL读进去并且按固定轴优化呢

科研菜菜

chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

OPTCELL如何只是固定z轴，是不是只需要在第一行输入110就好

Xpeinny

请问一下，您测试过三斜体系吗？其中gamma角度和c的大小会受影响发生改变吗？我测试出来c轴大小不变，但是gamma角度会发生变化...您有什么办法固定gamma角度不变吗？其中OPTCELL设置为110,110,000

Xpeinny

啦啦黑还黑 发表于 2018-11-13 22:11
楼主的方法只能限制Z方向，可以简单写一个3*3的矩阵文件OPTCELL限制任意的晶格矢量。
      INTEGER ICELL ...

请问一下，您这种矩阵写法的话，是不是不会限制晶胞类型啊？比我下面这种写法的优势在哪里呢？
     SUBROUTINE CONSTR_CELL_RELAX(FCELL)
      USE PREC
      REAL(Q) FCELL(3,3), SAVE(3)
      LOGICAL FILFLG
      INTEGER ICELL(3)
      INQUIRE(FILE='OPTCELL',EXIST=FILFLG)
      IF (FILFLG) THEN
        OPEN(67,FILE='OPTCELL',FORM='FORMATTED',STATUS='OLD')
          READ(67,"(3I1)") (ICELL(I),I=1,3)
        CLOSE(67)
        DO I=1,3
          SAVE(I)=FCELL(I,I)
        ENDDO
        FCELL=0.0d0
        DO I=1,3
          IF (ICELL(I)==1) FCELL(I,I)=SAVE(I)
        ENDDO
      ENDIF

Superconduct-BB

chempeng 发表于 2018-11-16 16:53
参考修改下

多谢

chempeng

Superconduct-BB 发表于 2018-11-14 21:16
学长能否发下你用来循环优化的脚本我刚才试着写了个，发现似乎把变量成到POSCAR的基失里并没有效果。

参考修改下

1. #!/bin/sh
   
2. rm WAVECAR
   
3. for i in   1.00 1.02 1.04 1.06 1.08 1.10 1.12
   
4. do
   
5. rm POSCAR
   
6. j=$(echo "$i*3" |bc)
   
7. cat > POSCAR <<!
   
8. 1-1
   
9. 1.0
   
10.         $j         0.0000000000         0.0000000000
    
11.         0.0000000000        $j         0.0000000000
    
12.         0.0000000000         0.0000000000         20
    
13. xxx
    
14. xxx
    
15. xxx
    
16. !
    
17. echo "c = $i Ang"
    
18. mpirun -np 8 vasp > out </dev/null
    
19. cp CONTCAR CONTCAR-$i
    
20. E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
    
21. V=`grep "volume" OUTCAR | tail -1 | awk '{printf "%12.4f \n", $5 }'`
    
22. cat OUTCAR | awk '$1=="in"&&$2=="kB" {print $3,$4,$5}' >> stress
    
23. echo $i $V $E  >> energy-strain
    
24. done
    

复制代码

Superconduct-BB

chempeng 发表于 2018-11-14 11:37
是的，一个 for 循环脚本就可以了。脚本可以参考候博的 VASP 软件包的使用入门指南

学长能否发下你用来循环优化的脚本我刚才试着写了个，发现似乎把变量成到POSCAR的基失里并没有效果。

Superconduct-BB

chempeng 发表于 2018-11-14 11:37
是的，一个 for 循环脚本就可以了。脚本可以参考候博的 VASP 软件包的使用入门指南

与侯老师的脚本有所不同。应该是直接把变量i×到a,b基失里吧，毕竟c不用变。