接下來(lái)進(jìn)行的這一步，和之前 Vasp wiki 上面的教程 Fcc Si 優(yōu)化晶格常數(shù) 基本過(guò)程是一樣的。
首先，準(zhǔn)備好 INCAR, KPOINTS, POTCAR 文件，POTCAR 與 bulk 結(jié)構(gòu)的一樣顷蟆。
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，而三維計(jì)算的時(shí)候我們是讓 ISIF=3缎患，差別在于是否可以改變?cè)男螤詈腕w積慕的。
修改晶格常數(shù)計(jì)算
接下來(lái)要做的主要思路是：由 CrI3 bulk 結(jié)構(gòu)切得到的單層的POSCAR (上一篇博客中提到)，寫(xiě)腳本挤渔，目的是每次修改 POSCAR中的晶格常數(shù)值肮街，然后進(jìn)行計(jì)算，得到 total free energy判导，最后能量最低的晶格常數(shù)值對(duì)應(yīng)最佳晶格常數(shù)嫉父。
下面介紹具體操作。
下圖中的藍(lán)色菱形線(xiàn)條框起來(lái)的就是單層CrI3原胞眼刃，綠色的原子表示的是 Cr绕辖，紅色和藍(lán)色的原子是 I，分別在上一層和下一層擂红。
圖中 a1, a2 就是晶格基矢仪际，設(shè)晶格常數(shù)為 a（即圖中棱形的邊長(zhǎng)），則 a1=(a,0,0), a2 = (-a/2, sqrt(3)*a/2, 0) 昵骤，a3 = (0, 0, c)树碱。
以下的 python2 腳本實(shí)現(xiàn)了根據(jù)晶格常數(shù)得到 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()
還需要修改在服務(wù)器提交計(jì)算任務(wù)的腳本 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
也就是說(shuō)变秦，在提交任務(wù)的腳本中成榜，我們選擇嘗試晶格常數(shù)從 5.6~7.2，每隔0.2取一個(gè)值計(jì)算蹦玫，并且調(diào)用了 pos.py赎婚。計(jì)算完成后我們又建立了以晶格常數(shù)命名的文件夾刘绣，把剛剛的計(jì)算結(jié)果移到新的文件夾，以免被后面的計(jì)算結(jié)果所覆蓋挣输。
之所以取這些晶格常數(shù)嘗試纬凤，是因?yàn)橹澳瞧恼陆o出的計(jì)算結(jié)果是 7.0008（然而我看走眼了，看成了 6.051)歧焦。大概每一個(gè)晶格常數(shù)的需要算3個(gè)多小時(shí)移斩，所以算這么多還是挺花時(shí)間的肚医，大概一兩天绢馍。（遇到服務(wù)器關(guān)機(jī)，計(jì)算節(jié)點(diǎn)存儲(chǔ)滿(mǎn)了什么的肠套，則需要更多的時(shí)間 = =）
根據(jù)計(jì)算結(jié)果畫(huà)圖
算完以后舰涌，可以取更密的值繼續(xù)計(jì)算。然后再寫(xiě)一個(gè) python3 腳本你稚，從不同的OSZICAR中提取每次計(jì)算的能量
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 如下瓷耙，第一列是晶格常數(shù)
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
再寫(xiě)一個(gè)作圖腳本 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')
就能得到所需的晶格常數(shù)、自由能量關(guān)系曲線(xiàn)圖
放大一點(diǎn)
所以刁赖，最低能量對(duì)應(yīng)的最佳晶格常數(shù)就是7.0搁痛，和我們上一篇博客中的 POSCAR完全一樣 (⊙﹏⊙)，和論文中用PBE functional 得到的 7.008 也是非常的一致宇弛。
