接下来进行的这一步，和之前 Vasp wiki 上面的教程 Fcc Si 优化晶格常数 基本过程是一样的。
首先，准备好 INCAR, KPOINTS, POTCAR 文件，POTCAR 与 bulk 结构的一样。
KPOINTS
k-points
0
Monkhorst Pack
8 8 1
0 0 0
INCAR
SYSTEM = CrI3 single layer
ISIF = 2
NSW = 1000
EDIFF = 1E-5
EDIFFG = -5E-3
PREC = High
IBRION = 2
ISMEAR = 0
ISPIN = 2
ENCUT = 600
注意 INCAR 这里 ISIF = 2，而三维计算的时候我们是让 ISIF=3，差别在于是否可以改变元胞的形状和体积。
修改晶格常数计算
接下来要做的主要思路是：由 CrI3 bulk 结构切得到的单层的POSCAR (上一篇博客中提到)，写脚本，目的是每次修改 POSCAR中的晶格常数值，然后进行计算，得到 total free energy，最后能量最低的晶格常数值对应最佳晶格常数。
下面介绍具体操作。
下图中的蓝色菱形线条框起来的就是单层CrI3原胞，绿色的原子表示的是 Cr，红色和蓝色的原子是 I，分别在上一层和下一层。
图中 a1, a2 就是晶格基矢，设晶格常数为 a（即图中棱形的边长），则 a1=(a,0,0), a2 = (-a/2, sqrt(3)*a/2, 0) ，a3 = (0, 0, c)。
以下的 python2 脚本实现了根据晶格常数得到 POSCAR 的功能
import sys
import math
a = float(sys.argv[1]) # lattice constant
fout = open('POSCAR', 'w')
fout.write(
'''CrI3 monolayer
1.0
%.10f 0.0000000000 0.0000000000
%.10f %.10f 0.0000000000
0.0000000000 0.0000000000 23.2178993225
I Cr
6 2
Direct
0.318349988 0.333710033 0.138589990
0.666289984 0.984630016 0.138589990
0.015369989 0.681649979 0.138589990
0.014990008 0.332949996 0.000000000
0.667049973 0.682030038 0.000000000
0.317970037 0.985010075 0.000000000
0.000000000 0.000000000 0.068870003
0.333330026 0.666670051 0.069729999
'''%(a, -a/2, math.sqrt(3)*a/2))
fout.close()
还需要修改在服务器提交计算任务的脚本 sub_vasp，在末尾 # running program 前后加上一段
54 rm WAVECAR SUMMARY.fcc
55 for i in 5.6 5.8 6.0 6.2 6.4 6.8 7.0 7.2; do
56 python pos.py $i
57
58 # running program
59 $OPEN_MPI $IB_FLAG -np $NCPUS -machinefile .NODES_to_RUN.${job id} $VASP_EXEC
60
61 # E=`awk '/F=/ {print $0}' OSZICAR` ; echo $i $E >>SUMMARY.fcc
62 mkdir $i
63 mv CHG CONTCAR EIGENVAL OSZICAR PCDAT vasprun.xml WAVECAR CHGCAR DOSCAR IBZKPT OUTCAR POSCAR XDATCAR $i
64
65 done
也就是说，在提交任务的脚本中，我们选择尝试晶格常数从 5.6~7.2，每隔0.2取一个值计算，并且调用了 pos.py。计算完成后我们又建立了以晶格常数命名的文件夹，把刚刚的计算结果移到新的文件夹，以免被后面的计算结果所覆盖。
之所以取这些晶格常数尝试，是因为之前那篇文章给出的计算结果是 7.0008（然而我看走眼了，看成了 6.051)。大概每一个晶格常数的需要算3个多小时，所以算这么多还是挺花时间的，大概一两天。（遇到服务器关机，计算节点存储满了什么的，则需要更多的时间 = =）
根据计算结果画图
算完以后，可以取更密的值继续计算。然后再写一个 python3 脚本，从不同的OSZICAR中提取每次计算的能量
sum.py
a = [5.6, 5.8, 6.0, 6.2, 6.4, 6.8, 6.9, 7.0, 7.1, 7.2, 7.4]
with open("SUMMARY.txt", "wt") as f:
for i in a:
with open("%.1f/OSZICAR"%i) as fin:
for line in fin:
pass
print(i,line[:-1], file=f)
得到 SUMMARY.txt 如下，第一列是晶格常数
5.6 10 F= -.25231765E+02 E0= -.25185740E+02 d E =-.829857E-03 mag= 5.3969
5.8 12 F= -.27302282E+02 E0= -.27267145E+02 d E =-.256012E-03 mag= 6.0472
6.0 11 F= -.28870757E+02 E0= -.28850240E+02 d E =-.123916E-04 mag= 6.0440
6.2 10 F= -.30009049E+02 E0= -.30004278E+02 d E =0.422491E-05 mag= 6.0047
6.4 13 F= -.30774062E+02 E0= -.30773558E+02 d E =-.483616E-05 mag= 6.0001
6.8 9 F= -.31488252E+02 E0= -.31488244E+02 d E =-.170946E-03 mag= 6.0000
6.9 9 F= -.31544414E+02 E0= -.31544408E+02 d E =-.308623E-03 mag= 6.0000
7.0 11 F= -.31562794E+02 E0= -.31562788E+02 d E =-.100152E-03 mag= 6.0000
7.1 14 F= -.31547196E+02 E0= -.31547189E+02 d E =-.146476E-04 mag= 6.0000
7.2 10 F= -.31500496E+02 E0= -.31500486E+02 d E =-.181810E-03 mag= 6.0000
7.4 11 F= -.31322697E+02 E0= -.31322675E+02 d E =-.108923E-02 mag= 6.0000
再写一个作图脚本 plot.py
import matplotlib.pyplot as plt
a = []
energy = []
with open("SUMMARY.txt", "r") as fin:
for line in fin:
mylist = line.split()
a.append(mylist[0])
energy.append(mylist[5])
plt.plot(a, energy, ".-")
plt.ylabel('energy')
plt.xlabel('a')
plt.show()
#plt.savefig('energy.png')
就能得到所需的晶格常数、自由能量关系曲线图
放大一点
所以，最低能量对应的最佳晶格常数就是7.0，和我们上一篇博客中的 POSCAR完全一样 (⊙﹏⊙)，和论文中用PBE functional 得到的 7.008 也是非常的一致。
