计算化学公社

标题: VASP 固定基矢优化结构方法 [打印本页]
作者
Author: chempeng    时间: 2018-3-21 11:21
标题: VASP 固定基矢优化结构方法
本帖最后由 chempeng 于 2018-3-21 19:41 编辑

一直以来在本站学习颇多，感谢 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

  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 基矢矩阵的更新流程

作者
Author: Superconduct-BB    时间: 2018-11-13 09:38
学长你好，请问你帖子里发的优化二维材料用侯老师说的方法对a=b改变晶格常数用ISLF=2优化是否要用到脚本，如果手动改是不是太麻烦了。
作者
Author: 啦啦黑还黑    时间: 2018-11-13 22:11
楼主的方法只能限制Z方向，可以简单写一个3*3的矩阵文件OPTCELL限制任意的晶格矢量。
      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
作者
Author: chempeng    时间: 2018-11-14 11:37
Superconduct-BB 发表于 2018-11-13 09:38
学长你好，请问你帖子里发的优化二维材料用侯老师说的方法对a=b改变晶格常数用ISLF=2优化是否要用到脚本， ...

是的，一个 for 循环脚本就可以了。脚本可以参考候博的 VASP 软件包的使用入门指南
作者
Author: chempeng    时间: 2018-11-14 11:48
啦啦黑还黑 发表于 2018-11-13 22:11
楼主的方法只能限制Z方向，可以简单写一个3*3的矩阵文件OPTCELL限制任意的晶格矢量。
      INTEGER ICELL ...


补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下
  1. $ cat OPTCELL
  2. 100
  3. 110
  4. 000
复制代码
作者
Author: Superconduct-BB    时间: 2018-11-14 13:56
chempeng 发表于 2018-11-14 11:37
是的，一个 for 循环脚本就可以了。脚本可以参考候博的 VASP 软件包的使用入门指南

与侯老师的脚本有所不同。应该是直接把变量i×到a,b基失里吧，毕竟c不用变。
作者
Author: Superconduct-BB    时间: 2018-11-14 21:16
chempeng 发表于 2018-11-14 11:37
是的，一个 for 循环脚本就可以了。脚本可以参考候博的 VASP 软件包的使用入门指南

学长能否发下你用来循环优化的脚本我刚才试着写了个，发现似乎把变量成到POSCAR的基失里并没有效果。
作者
Author: chempeng    时间: 2018-11-16 16:53
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
复制代码
作者
Author: Superconduct-BB    时间: 2018-11-17 19:27
chempeng 发表于 2018-11-16 16:53
参考修改下

多谢
作者
Author: Xpeinny    时间: 2019-3-26 10:24
啦啦黑还黑 发表于 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
作者
Author: Xpeinny    时间: 2019-3-26 12:14
请问一下，您测试过三斜体系吗？其中gamma角度和c的大小会受影响发生改变吗？我测试出来c轴大小不变，但是gamma角度会发生变化...您有什么办法固定gamma角度不变吗？其中OPTCELL设置为110,110,000
作者
Author: 科研菜菜    时间: 2020-8-20 16:08
chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

OPTCELL如何只是固定z轴，是不是只需要在第一行输入110就好
作者
Author: luchaolu    时间: 2022-1-4 14:13
chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

请问这个OPTCELL文件怎么设置呢？？怎么能够让脚本将设置的OPTCELL读进去并且按固定轴优化呢
作者
Author: luchaolu    时间: 2022-1-4 14:15
Xpeinny 发表于 2019-3-26 12:14
请问一下，您测试过三斜体系吗？其中gamma角度和c的大小会受影响发生改变吗？我测试出来c轴大小不变，但是g ...

OPTCELL固定c轴这样对吗？不是应该为100 010 000吗
作者
Author: shameless    时间: 2022-6-18 11:19
luchaolu 发表于 2022-1-4 14:15
OPTCELL固定c轴这样对吗？不是应该为100 010 000吗

你好 固定两个轴的OPTCELL 应该怎么写
作者
Author: fjy    时间: 2022-6-30 17:01
shameless 发表于 2022-6-18 11:19
你好 固定两个轴的OPTCELL 应该怎么写

把要固定的两个轴都写成0就好了
作者
Author: helpcal    时间: 2023-12-29 05:04
请问直接isif=3固定c轴优化和调整ab轴画E-a曲线得到的晶格常数相差约0.002，这时候应该取哪一个值更准确。

另外E-a曲线方法是不是取能量最低点然后对其进行isif-2的优化？这个方法里我对近邻几个a也优化原子位置，结果发现别的a驰豫完能量更低一点点。
作者
Author: 汪杰    时间: 2024-10-29 14:43
啦啦黑还黑 发表于 2018-11-13 22:11
楼主的方法只能限制Z方向，可以简单写一个3*3的矩阵文件OPTCELL限制任意的晶格矢量。
      INTEGER ICELL ...

请问一下，这样改之后，OPTCELL应该怎么设置呢？能不能给一个具体的例子
作者
Author: 汪杰    时间: 2024-10-29 14:47
chempeng 发表于 2018-11-14 11:48
补充提示下，计算前编写输入文件 OPTCELL ，对应格矢，1为优化，0为固定，如下

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3