3.3 电声耦合

发表于 2023-05-09 更新于 2024-06-15 分类于 DFT学习

MC-电声耦合

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

PHON_LMC = .TRUE.
IBRION = 6

全MC采样

标签PHON_NSTRUCT设置由于MC抽样而生成的结构的数量。应该监视观测值相对于这个数量的收敛。
标签TEBEG=0也是需要的，用于选择运行抽样的温度。
另外，PHON_LBOSE可以设置为.TRUE.或.FALSE. (默认PHON_LBOSE=.TRUE.)，分别选择玻色-爱因斯坦统计或麦克斯韦-玻尔兹曼统计。

0K的一个INCAR文件如下:

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毫电子伏特的精度已经足够时。

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

给出了一个温度范围为0-700K（步长为100K）的示例INCAR文件：

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.

输出

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

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

one-shot例子（金刚石带隙重整化）

输入文件

INCAR：

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 修改为

 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0 
SIGMA = 0.1  #去掉PHON_标签

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

三、提取ZPR

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

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

#!/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)}'

输出如下：

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：

 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 修改为

 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0 
SIGMA = 0.1  #去掉PHON_标签

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

#!/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

!/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

#! /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数据处理。
脚本内容：

#!/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获得温度-体积关系。
脚本内容：

#!/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

输出如下：

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.sh 和extract_temp.sh.获得数据画图。

通过添加体积效应，与实验获得了更好的一致性。在这个教程中，我们只使用了一个4x4x4的单元格，因为更大的单元格已经相当耗时，但对于收敛的5x5x5单元格，两条曲线都应该与实验相比略有变差。由于在本例中使用的PBE无法充分描述电子交换和相关性，因此预计实验与理论之间会有差异。为了获得真正的优异一致性，需要使用GW近似方法。
严格来说，将体积效应添加到电子-声子相互作用的正确方法是首先为每个温度更改体积，然后为该温度计算电子-声子相互作用。在这个教程和参考文献[5]中，我们是反过来做的。因此，电子-声子相互作用只需要计算一次。在参考文献[5]中，我们观察到这两种方法给出了非常相似的结果。