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