单层 CrI3 晶格常数优化

文章目录

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

接下来进行的这一步，和之前 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 也是非常的一致。