利用VASP进行弹性模量的计算笔记

本文是基于《基于能量-应变关系计算弹性常数—VASPKIT v1.00新功能》、《基于应力-应变关系计算弹性常数—VASPKIT v1.2.0新功能》，并在此基础上进行整理记录，在此对原分享者表示感谢。

1 准备优化彻底的POSCAR文件，注意通常采用标准的惯用原胞计算弹性常数，如果不确信POSCAR文件中是否是标准的惯用原胞，可以用vaspkit-603/604生成标准结构；

补充：在Materials Project中下载所需材料的CIF文件，利用VESTA导出为XX.VASP文件，重命名为POSCAR

2 结构优化（弛豫）

INCAR 获取：vaspkit → 101 → LR

KPOINTS 获取： vaspkit →102→ 2→ 0.03 （由于计算弹性常数对K-mesh要求很高，因此对于半导体（金属体）体系，生成K点的精度应不小于0.03(0.02) * 2π Å-1）

2.1 如果弛豫发现报错

最好重新走一遍结构优化，或者INCAR中加上【ISYM=0】忽略对称性（但是增加计算量）

3 弹性模量计算文件准备

POSCAR : 使用结构优化后的新坐标， cp CONTCAR POSCAR

KPOINTS : vaspkit →102 →2 →0.03

INCAR : vaspkit →101 →DC , 可以除了PREC参数，也可以手动设置ENCUT=1.5*ENMAX

Global Parameters
ISTART = 0        
LREAL = F      
PREC   = High （截断能设置默认值1.5-2倍）
LWAVE = F      
LCHARG = F    
ADDGRID= .TRUE.    
Electronic Relaxation
ISMEAR = 0          
SIGMA = 0.05      
NELM   = 40          
NELMIN = 4          
EDIFF = 1E-08      
Ionic Relaxation
NELMIN = 6          
NSW   = 100        
IBRION = 2          
ISIF   = 2   （切记选择2,如果选择3会把施加应变后原胞重新优化成平衡原胞）
EDIFFG = -1E-02

4 准备VPKIT.in文件，内容如下

1                   ! 设置1将产生计算弹性常数的输入文件，2则计算弹性常数
3D                   ! 2D为二维体系，3D为三维体系
7                   ! 7个应变
-0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 ! 应变变化范围

运行vaspkit 200（或201），分别对应【用应力-应变法计算弹性常数】、【能量-应变法计算弹性常数】

根据POSCAR产生对应张量元文件夹和计算文件

5 准备批量提交vasp作业的脚本（txml.sh）

#这条代码运行后，会在每个张量元的应变计算后，把超算脚本vasp.job复制到下一文件夹，继续张量元的计算
#!/bin/bash
root_path=`pwd`
for cij in `ls -F | grep /$`
do
cd ${root_path}/$cij
for s in strain_*
do
cd ${root_path}/$cij/$s
echo `pwd`
    cp ../../vasp.job .            
    ./vasp.job
# 上面两行vasp.job 替换为 超算运算的脚本文件名
# Add here your vasp_submit_job_script     把这一行替换为提交运算的命令
done
done

上述脚本要授予权限，否则如下图

利用chmod命令：chmod 777 txml.sh
6 计算结束后，再次修改VPKIT.in文件中第一行为2（后处理），然后再次运行vaspkit并选择200/201，即可得到弹性模量矩阵