Contents
  1. 1. 修改晶格常数计算
  2. 2. 根据计算结果画图

接下来进行的这一步,和之前 Vasp wiki 上面的教程 Fcc Si 优化晶格常数 基本过程是一样的。
首先,准备好 INCAR, KPOINTS, POTCAR 文件,POTCAR 与 bulk 结构的一样。

KPOINTS

1
2
3
4
5
k-points
0
Monkhorst Pack
8 8 1
0 0 0

INCAR

1
2
3
4
5
6
7
8
9
10
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,差别在于是否可以改变元胞的形状和体积。

ISIF的设置

修改晶格常数计算

接下来要做的主要思路是:由 CrI3 bulk 结构切得到的单层的POSCAR (上一篇博客中提到),写脚本,目的是每次修改 POSCAR中的晶格常数值,然后进行计算,得到 total free energy,最后能量最低的晶格常数值对应最佳晶格常数。

下面介绍具体操作。

下图中的蓝色菱形线条框起来的就是单层CrI3原胞,绿色的原子表示的是 Cr,红色和蓝色的原子是 I,分别在上一层和下一层。

图片 Zhang, Wei-Bing, et al. Journal of Materials Chemistry C 3.48 (2015): 12457-12468

图中 a1, a2 就是晶格基矢,设晶格常数为 a(即图中棱形的边长),则 a1=(a,0,0), a2 = (-a/2, sqrt(3)*a/2, 0) ,a3 = (0, 0, c)

以下的 python2 脚本实现了根据晶格常数得到 POSCAR 的功能

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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 前后加上一段

1
2
3
4
5
6
7
8
9
10
11
12
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

1
2
3
4
5
6
7
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 如下,第一列是晶格常数

1
2
3
4
5
6
7
8
9
10
11
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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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 也是非常的一致。

Contents
  1. 1. 修改晶格常数计算
  2. 2. 根据计算结果画图