使用vasp计算α-SiO2弹性模量(个人记录)

使用vasp计算α-SiO₂弹性模量

文件下载

由于vasp wiki的example里有alpha-SiO₂的样例，故INCAR、KPOINTS、POSCAR、POTCAR文件可以从vasp wiki直接下载

样例链接：https://www.vasp.at/wiki/index.php/Alpha-SiO2

输入文件下载链接：https://www.vasp.at/wiki/images/4/4e/SiO2_NMR.tgz

输入文件选用

因为需要计算弹性模量，所以修改（添加）如下INCAR文件参数

IBRION = 6
ISIF   = 3
NFREE  = 4
NSW    = 1

ENCUT收敛性测试

选用官网样例中默认的KPOINTS.4x4x3作为KPOINTS文件，测试ENCUT收敛性

4x4x3
 0
G
 4 4 3
 0 0 0

修改INCAR中的ENCUT参数为AAA并另存为INCAR0，编写bash如下，测试ENCUT收敛性, 将不同ENCUT值时的能量写入energy文件

#!/bin/bash
for i in 200 250 300 350 400 450 500 550 600 # 候选ENCUT值
    do
        sed "s/AAA/$i/g" INCAR0 > INCAR # 文本替换
        mkdir $i
        cp vasp_std INCAR POTCAR POSCAR KPOINTS $i/
        cd $i
        ulimit -s unlimited # 关闭堆栈大小限制
        mpirun -np 8 vasp_std 2>&1 | tee log # vasp主程序运行
        echo $i
        cd ../
    done
    
echo "#encut   Energy" > energy
for i in 200 250 300 350 400 450 500 550 600
    do
        E=`grep -A 2 "free\ energy" $i/OUTCAR | tail -1 | cut -c33-44` # 字符串处理，读取能量值
        echo $i " " $E >> energy
    done

计算得到energy文件如下

#encut   Energy
200   -69.62312758
250   -71.76037373
300   -70.73853595
350   -69.96559936
400   -69.75761344
450   -69.65430805
500   -69.61260083
550   -69.62444570
600   -69.64168679

规律不明显，故选取ENCUT=400作为截断能（和官网样例相同）

KPOINTS收敛性测试

如上一步所示，选取ENCUT=400作为截断能

修改KPOINTS中的晶格参数为AAA并另存为KPOINTS0，编写bash如下，测试KPOINTS收敛性, 将不同KPOINTS值时的能量写入energy文件

#!/bin/bash
for i in "3 3 2" "4 4 3" "5 5 4" # 几个候选KPOINTS值
    do
        sed "s/AAA/$i/g" KPOINTS0 > KPOINTS # 文本替换
        mkdir "$i"
        cp vasp_std INCAR POTCAR POSCAR KPOINTS "$i"/
        cd "$i"
        ulimit -s unlimited # 关闭堆栈大小限制
        mpirun -np 8 ./vasp_std 2>&1  | tee log # 主程序运行
        echo $i
        cd ../
    done

echo "#k-points    Energy" >energy
for i in "3 3 2" "4 4 3" "5 5 4"
    do
        E=`grep -A 2 "free\ energy" "$i"/OUTCAR | tail -1 | cut -c33-44` # 字符串处理，读取能量值
        echo $i "     " $E >> energy
    done

计算得到energy文件如下

#k-points    Energy
3 3 2       -69.75780682
4 4 3       -69.75761344
5 5 4       -69.76029387

规律不明显，故选取4*4*3作为KPOINTS（也和官网样例相同）

结构弛豫计算 优化晶格

修改INCAR文件中的以下参数如下，准备做结构弛豫计算

IBRION = 6
NFREE = 4
ISIF = 3
NSF = 100
ENCUT = 520 # 400*1.3

计算后使用vimdiff对比POSCAR、CONTCAR两个文件的差异如下

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hKOnaCPG-1665304896101)(/home/ian/.config/Typora/typora-user-images/image-20221009160546639.png)]

但由于结构优化最终并未收敛，故并未采取CONTCAR优化过的晶格替代初始POSCAR进行最终计算

最终计算 计算弹性模量

修改INCAR文件中的以下参数如下，进行力学性能计算，计算弹性模量

IBRION = 6
NFREE = 4
ISIF = 3
NSW = 1
ENCUT = 400

计算完成后可在OUTCAR文件中找到弹性矩阵如下：

 TOTAL ELASTIC MODULI (kBar)
 Direction    XX          YY          ZZ          XY          YZ          ZX
 --------------------------------------------------------------------------------
 XX        1866.1602    729.1319    756.9221     -0.0000     44.6970     -0.0000
 YY         729.1319   1866.1602    756.9221      0.0000    -44.6970      0.0000
 ZZ         756.9221    756.9221   2800.3034     -0.0000      0.0000      0.0000
 XY          -0.0000      0.0000      0.0000    568.5142     -0.0000     44.6970
 YZ          44.6970    -44.6970      0.0000     -0.0000    835.5214     -0.0000
 ZX          -0.0000      0.0000      0.0000     44.6970     -0.0000    835.5214
 --------------------------------------------------------------------------------

之后便可以使用力学性质在线作图工具（链接：http://progs.coudert.name/elate）读取数据并作图