3.3 电声耦合

MC-电声耦合

需要vasp6.0之后的版本,INCAR开启

PHON_LMC = .TRUE.
IBRION = 6

全MC采样

  • 标签PHON_NSTRUCT设置由于MC抽样而生成的结构的数量。应该监视观测值相对于这个数量的收敛。

  • 标签TEBEG=0也是需要的,用于选择运行抽样的温度。

  • 另外,PHON_LBOSE可以设置为.TRUE.或.FALSE. (默认PHON_LBOSE=.TRUE.),分别选择玻色-爱因斯坦统计或麦克斯韦-玻尔兹曼统计。

0K的一个INCAR文件如下:

1
2
3
4
5
6
7
8
System = DEFAULT
PREC = Accurate
ISMEAR = 0; SIGMA = 0.1;
IBRION = 6

PHON_LMC = .TRUE.
PHON_NSTRUCT = 100
TEBEG = 0.0

随后生成Wycoff positions畸变的许多POSCAR,输出为POSCAR.TEBEG.NUMBER

ZG configuration即时采样(one-shot)

这种方法仅使用单一的扭曲结构,因此比完全的MC抽样快上几个数量级,同时保持了与收敛的超晶胞尺寸的完全MC抽样非常接近的准确性。例如,我们展示了对于带隙的零点重整化,准确度在ZG配置和完全MC抽样之间在5毫电子伏特之内。因此,我们建议优先使用这种方法,当很难实现超晶胞尺寸的收敛或者5毫电子伏特的精度已经足够时。

  1. 输入
    要选择ZG配置,必须在INCAR文件中设置PHON_NSTRUCT=0。不同温度的数量和温度列表(以K为单位)必须使用标签PHON_NTLISTPHON_TLIST,在INCAR文件中分别提供。一个示例如下:
    1
    2
    PHON_NTLIST = 4
    PHON_TLIST = 0.0 100.0 200.0 350.0

给出了一个温度范围为0-700K(步长为100K)的示例INCAR文件:

1
2
3
4
5
6
7
8
System = DEFAULT
PREC = Accurate
ISMEAR = 0; SIGMA = 0.1;
IBRION = 6
PHON_NTLIST = 8
PHON_TLIST = 0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0
PHON_NSTRUCT = 0
PHON_LMC = .TRUE.
  1. 输出

类似于MC抽样,ZG配置方法会产生多个具有不同扭曲Wycoff位置但不变的布里渊矩阵的POSCAR文件。这些文件被标记为

POSCAR.TEMP
其中TEMP遍历由PHON_TLIST定义的所有温度。

one-shot例子(金刚石带隙重整化)

输入文件

INCAR:

1
2
3
4
5
6
7
8
9
10
11
general:
System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0
SIGMA = 0.1;
IBRION = 6 #获得动力学矩阵声子计算
PHON_LMC = .TRUE. #开启电声耦合计算
PHON_NSTRUCT = 0 #选择one-shot方法
PHON_NTLIST = 1 #温度点个数
PHON_TLIST = 0.0 #对应温度

准备优化后的结构,需要扩胞,以保证声子准确。KPOINTS、POTCAR根据POSCAR调整。

计算

一、获得结构

按照超胞结构和上面的INCAR提交任务,结束后获得OUTCAR。
执行: cp OUTCAR OUTCAR.init
新的畸变后的结构输出为POSCAR.T=0.

二、计算特定位移畸变结构的电子能级

cp POSCAR.T=0. POSCAR

INCAR 修改为

1
2
3
4
5
 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0
SIGMA = 0.1 #去掉PHON_标签

执行计算,完成后cp OUTCAR OUTCAR.T=0.

三、 提取ZPR

具有特殊位移的带隙减去没有特殊位移的带隙,就可以得到零点带隙重整的值。带隙保存在OUTCAR文件里。

执行例子中的脚本:bash extract_zpr.sh,
脚本内容:

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
27
28
29
30
31
32
33
34
35
36
#!/bin/bash 

i="OUTCAR.T=0."
j="OUTCAR.init"

homo1=`awk '/NELECT/ {print $3/2}' $i`
homo2=`awk '/NELECT/ {print $3/2-1}' $i`
homo3=`awk '/NELECT/ {print $3/2-2}' $i`
lumo1=`awk '/NELECT/ {print $3/2+var+1}' $i`
lumo2=`awk '/NELECT/ {print $3/2+var+2}' $i`
lumo3=`awk '/NELECT/ {print $3/2+var+3}' $i`
lumo4=`awk '/NELECT/ {print $3/2+var+4}' $i`
lumo5=`awk '/NELECT/ {print $3/2+var+5}' $i`
lumo6=`awk '/NELECT/ {print $3/2+var+6}' $i`
e1a=`grep " $homo1 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1b=`grep " $homo2 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1c=`grep " $homo3 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e2a=`grep " $lumo1 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2b=`grep " $lumo2 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2c=`grep " $lumo3 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2d=`grep " $lumo4 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2e=`grep " $lumo5 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2f=`grep " $lumo6 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`

homo_ref=`awk '/NELECT/ {print $3/2}' $j`
lumo_ref=`awk '/NELECT/ {print $3/2+var+1}' $j`

h_ref=`grep " $homo_ref " $j | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
l_ref=`grep " $lumo_ref " $j | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`

echo "The band gap (in eV) without zero-point vibrations is:"
echo "$h_ref $l_ref" |awk '{print ($2-$1)}'
echo "The band gap (in eV) including zero-point vibrations is:"
echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}'
echo "The zero-point renormalization of the band gap (in eV) is:"
echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $h_ref $l_ref" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-($11-$10)}'

输出如下:

1
2
3
4
5
6
The band gap (in eV) without zero-point vibrations is:
4.4049
The band gap (in eV) including zero-point vibrations is:
4.05102
The zero-point renormalization of the band gap (in eV) is:
-0.353883

精度

精度取决于超胞大小。超胞越大越准确。

温度依赖的带隙计算

输入文件

INCAR:

1
2
3
4
5
6
7
8
9
10
 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0
SIGMA = 0.1;
IBRION = 6 #获得动力学矩阵声子计算
PHON_LMC = .TRUE. #开启电声耦合计算
PHON_NSTRUCT = 0 #选择one-shot方法
PHON_NTLIST = 8 #温度个数
PHON_TLIST = 0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0

准备优化后的结构,需要扩胞,以保证声子准确。KPOINTS、POTCAR根据POSCAR调整。

计算

一、 获取特定位移的结构

按照超胞结构和上面的INCAR提交任务,结束后获得OUTCAR。
执行: cp OUTCAR OUTCAR.init
新的畸变后的结构输出为POSCAR.T=0. ~ POSCAR.T=700.

二、计算特定位移畸变结构的电子能级
cp POSCAR.T=0. POSCAR

INCAR 修改为

1
2
3
4
5
 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0
SIGMA = 0.1 #去掉PHON_标签

对每个温度下的结构执行计算。可通过脚本run_temperature.sh批量提交。
脚本内容:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#!/bin/bash

# please enter your executable path here
vasp_exec=./vasp_gam

#please enter the number of processors used for VASP here
np=8


for i in 0 100 200 300 400 500 600 700
do
cp POSCAR.T\=$i. POSCAR
mpirun -np $np $vasp_exec
mv OUTCAR OUTCAR.T\=$i
done

三、 分析OUTCAR

利用脚本extract_temp.sh

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
!/bin/bash 

if [ -f gap_vs_temp.dat ]
then
rm gap_vs_temp.dat
fi

touch gap_vs_temp.dat
counter=0

for temp in 0 100 200 300 400 500 600 700
do
i="OUTCAR.T=$temp"

homo1=`awk '/NELECT/ {print $3/2}' $i`
homo2=`awk '/NELECT/ {print $3/2-1}' $i`
homo3=`awk '/NELECT/ {print $3/2-2}' $i`
lumo1=`awk '/NELECT/ {print $3/2+var+1}' $i`
lumo2=`awk '/NELECT/ {print $3/2+var+2}' $i`
lumo3=`awk '/NELECT/ {print $3/2+var+3}' $i`
lumo4=`awk '/NELECT/ {print $3/2+var+4}' $i`
lumo5=`awk '/NELECT/ {print $3/2+var+5}' $i`
lumo6=`awk '/NELECT/ {print $3/2+var+6}' $i`
e1a=`grep "^ $homo1 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1b=`grep "^ $homo2 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1c=`grep "^ $homo3 " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e2a=`grep "^ $lumo1 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2b=`grep "^ $lumo2 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2c=`grep "^ $lumo3 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2d=`grep "^ $lumo4 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2e=`grep "^ $lumo5 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2f=`grep "^ $lumo6 " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`

if [ $temp -eq "0" ]
then
ref=`echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}'`
fi

echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $temp $ref" |awk '{print $10,(($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-$11}' >> gap_vs_temp.dat
done

运行脚本后生成gap_vs_temp.dat文件,画图即可。

gnuplot画图:
gnuplot -e "set terminal jpeg;set xlabel 'T (K)'; set ylabel 'band gap'; set style data lines; plot 'gap_vs_temp.dat', 'C_exp_points_offset0.dat' w circles, 'C_exp_fit_offset0.dat'" > gap_vs_temp.jpg

包含体积效应的温度依赖带隙

仅考虑温度效应会低估带隙,要考虑体积效应。

一、准简谐近似方法获得体积变化

三个脚本:

1
2
3
quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh
quasi_harm_4x4x4_diamond_make_energy_vs_volume_plots.sh
quasi_harm_4x4x4_diamond_obtain_fitting.sh

对超胞结构创建不同体积的结构进行体积-能量计算:

运行脚本:quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#! /bin/bash

# please enter your executable path here
vasp_exec=./vasp_gam

#please enter the number of processors used for VASP here
np=8

cp INCAR.qh INCAR

for i in 6.13521592 6.27789536 6.4205748 6.56325424 6.70593368 6.84861312 6.99129256 7.133972 7.27665144 7.41933088 7.56201032 7.70468976 7.8473692 7.99004864 8.13272808
do
sed "s/7.13397200/${i}/g" POSCAR.4x4x4 > POSCAR_$i
cp POSCAR_$i POSCAR
mpirun -np 8 $vasp_exec

mv OUTCAR OUTCAR_$i
done

该脚本创建15个POSCAR,体积每步变化2%,运行脚本quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh提交计算。

二、计算完成后运行quasi_harm_4x4x4_diamond_make_energy_vs_volume_plots.sh数据处理。
脚本内容:

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#!/bin/bash

bandshift=0
gwrun=-1
dgbd=-1
val=-1
con=-1
gwldadiff=-1
test=-1
while [[ $# -gt 0 ]]
do
key="$1"
case $key in
esac
shift
done

if [ -f "helpscript.perl" ]; then
rm helpscript.perl
fi


cat > helpscript.perl <<EOF
#!/bin/perl

use strict;
use warnings;

my \$zahler=0;
my @entropy=0;
my \$ezp=0;
my \$fhelmholtz=0;
my \$uenergy;
my \$kboltzmann=8.6173303*10**(-5.0);
my \$ntemp=8;
my \$tmax=700;
my \$tmin=0;
my \$tstep=100;
for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
{
\$entropy[\$itemp]=0;
}
while (<>)
{
chomp;
\$_=~s/^/ /;
my @help=split(/[\t,\s]+/);
\$zahler=\$zahler+1;
if (\$zahler == 1) {\$uenergy=\$help[1];}
else
{
my \$homega=\$help[2]/1000;
\$ezp=\$ezp+\$homega*0.5;
for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
{
my \$temp=(\$itemp-1)*\$tstep+\$tmin;
my \$kbt=\$kboltzmann*\$temp;
if (\$temp < 0.0000001)
{
}
else
{
\$entropy[\$itemp]=\$entropy[\$itemp]-\$kboltzmann*log(1-exp(-\$homega/\$kbt));
\$entropy[\$itemp]=\$entropy[\$itemp]+\$kboltzmann*\$homega/\$kbt*(1/(exp(\$homega/\$kbt)-1));
}
}
}
last if eof;
}

\$ezp=\$ezp; #/ \$zahler;

for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
{
my \$temp=(\$itemp-1)*\$tstep+\$tmin;

\$entropy[\$itemp]=\$entropy[\$itemp]; # / \$zahler;
\$fhelmholtz=\$uenergy+\$ezp-\$temp*\$entropy[\$itemp];
printf ("%15.8e %15.8e\n",\$temp,\$fhelmholtz);
}
printf ("#NOTEMP %15.8e\n",\$uenergy);
EOF

rm OUTTEMP*
first=0
for i in OUTCAR_*; do

echo "Starting $i"
v2=${i##OUTCAR_}
if [ -f "helpfile.help" ]; then
rm helpfile.help
fi
touch helpfile.help
cp OUTCAR_$v2 OUTCAR
awk '/free energy/' OUTCAR | tail -n 1| awk '{print $5}' >> helpfile.help
awk '/[0-9]* f .* THz/ {print $1,$10}' OUTCAR >> helpfile.help
awk '/[0-9]* f.i.*THz/ {print $1,$9}' OUTCAR >> helpfile.help
volume=`awk '/volume of cell/ {print $5}' OUTCAR | tail -n 1`

perl helpscript.perl helpfile.help > hhhhelp.txt

runcount=0
while read line; do
runcount=$((runcount+1))
if [[ $first -eq 0 ]]; then
echo $line | awk -v var="$volume" '{print var,$2,$1}' > OUTTEMP_$runcount
else
echo $line | awk -v var="$volume" '{print var,$2}' >> OUTTEMP_$runcount
fi
done < ./hhhhelp.txt

first=1
done

for i in OUTTEMP_*; do
v2=${i##OUTTEMP_}
mv OUTTEMP_$v2 OUTHELP
sort -n -k 1 OUTHELP > OUTTEMP_$v2
done

rm helpfile.help
rm helpscript.perl
rm hhhhelp.txt
rm OUTHELP

自由能曲线保存在OUTTEMP_*.

运行脚本quasi_harm_4x4x4_diamond_obtain_fitting.sh获得温度-体积关系。
脚本内容:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#!/bin/bash

for i in OUTTEMP_*
do
cp $i OUTTEMP.current
#extract temperature
temp=`head -n 1 OUTTEMP.current|awk '{print $3}'`
#do fitting
gnuplot -e "E(V)=E0+9.0/8.0*B0*V0*((V0/V)**(2.0/3.0)-1)**2 + 9.0/16.0*B0*\
(B0P-4)*V0*((V0/V)**(2.0/3.0)-1.0)**3.0 + R*((V0/V)**(2.0/3.0)-1.0)**4.0;\
B0P = 1;B0 = 1;V0 = 720;E0 = -1150;R = -1.0;fit E(x) 'OUTTEMP.current' u 1:2 via B0P,B0,V0,E0,R" &> suppress_output
#extract volume from fit
a=`grep "V0" fit.log|grep "=" |tail -n 1|awk '{print ($3/2.0)**(1.0/3.0)}'`
#print temperature and volume to
echo "temperature: $temp, a_latt: $a"
done

rm suppress_output
rm OUTTEMP.current

输出如下:

1
2
3
4
5
6
7
8
9
temperature: 0.00000000e+00, a_latt: 7.18012
temperature: 1.00000000e+02, a_latt: 7.18014
temperature: 2.00000000e+02, a_latt: 7.18064
temperature: 3.00000000e+02, a_latt: 7.18267
temperature: 4.00000000e+02, a_latt: 7.18613
temperature: 5.00000000e+02, a_latt: 7.19037
temperature: 6.00000000e+02, a_latt: 7.19493
temperature: 7.00000000e+02, a_latt: 7.19959
temperature: #NOTEMP, a_latt: 7.15218

三、 最后将QHA获得的晶胞参数,替换掉POSCAR.T=* files中的晶胞参数。再次运行run_temperature.shextract_temp.sh.获得数据画图。

  • 通过添加体积效应,与实验获得了更好的一致性。在这个教程中,我们只使用了一个4x4x4的单元格,因为更大的单元格已经相当耗时,但对于收敛的5x5x5单元格,两条曲线都应该与实验相比略有变差。由于在本例中使用的PBE无法充分描述电子交换和相关性,因此预计实验与理论之间会有差异。为了获得真正的优异一致性,需要使用GW近似方法

  • 严格来说,将体积效应添加到电子-声子相互作用的正确方法是首先为每个温度更改体积,然后为该温度计算电子-声子相互作用。在这个教程和参考文献[5]中,我们是反过来做的。因此,电子-声子相互作用只需要计算一次。在参考文献[5]中,我们观察到这两种方法给出了非常相似的结果。