用VASP计算硅的能带-(改进稿)

用VASP计算晶体硅的电子能带结构本文主要是用VASP进行了晶体硅的电子能带结构的计算。硅的晶体结构如图1所示，它是两个嵌套在一起的FCC布拉菲晶格，相对位置为(a/4,a/4,a/4)。在计算中，采用了FCC的原胞，每个原胞里面有两个硅原子。图1：硅的晶体结构首先优化了晶格参数，求出了能量最低所对应的晶格常数a。然后进行静态计算，得到了自洽的电荷密度。再由静态计算得到的自洽的电荷密度计算了硅的能带结构和电子态密度。（一）优化晶格常数只需要优化一个晶格常数a。首先准备四个文件：INCAR,POSCAR,POTCAR,KPOINTS，内容分别如下INCARPOSCARPOTCARKPOINTSsystem=Si-fccISTART=0ICHARGE=2ENCUT=300EDIFF=1E-5IBRION=2POTIM=0.25NSW=100EDIFFG=-1E-2ISMEAR=0SIGMA=0.1PREC=AccurateISIF=2Si-fcc5.10.00.50.50.50.00.50.50.50.02Direct0.00.00.00.250.250.25#采用PAW_PBE赝势auto0Monkhorst999000运行VASP，读取OUTCAR文件，可以得到a=5.1的情况下对应的能量。然后改变a的值，再运行VASP，得到相应晶格常数对应的能量。a的取值范围为5.1到5.8，每隔0.1取一个点。最终得到晶格常数与能量的对应值如下表所示晶格常数a(Å)Energy(eV)5.1-10.255601655.2-10.548330925.3-10.732725825.4-10.824812165.5-10.839103125.6-10.788285235.7-10.683956305.8-10.53611370将以上数据用MATLAB进行二次多项式拟合，求出拟合二次多项式的极小值点为5.50Å，即为优化的晶格常数，与实验值5.43Å比较接近。后面的计算均基于此优化的晶格常数。（二）静态计算自洽的电荷密度修改输入文件，在INCAR中定义NSW=0、LCHARG=T，POSCAR中的晶格常数定义为优化的晶格常数5.5，其它参数不变。运行VASP，即得到了自洽的面电荷密度，将其保存下来，以便后面计算能带结构和电子态密度之用。（三）计算能带结构选取6个特殊K点，特殊K点间的分割点数分别为20、20、20、10、20。从自洽的电荷密度计算得到的OUTCAR文件中可以找到倒格子基矢和费米能级。将以上信息写入syml文件，然后将syml文件输入到程序gk.x，产生K点，得到KPOINTS文件。另外设置INCAR中的ISTART=1，ICHAGE=11，NSW=0。POSCAR中的晶格常数定义为优化的晶格常数5.5。其他参数不变。运行VASP，计算完后得到本征值文件EIGENVAL。输入文件EIGENVAL和syml到程序pbnf.x，可得到输出文件bnd.dat和highk。将bnd.dat和highk中的数据导入到Origin，即可画出如图2所示的硅的电子能带结构图。导带底的能量为0.4045eV，价带顶的能量为-0.1422eV，故带隙为0.5647eV，大大低于实验值1.12eV。由于导带底和价带顶位于不同的K点处，故该带隙为间接带隙。图2：硅的电子能带结构图（四）计算电子态密度从POTCAR文件中找到硅的RWIGS=1.312，在INCAR中定义ISTART=1，ICHARGE=11，NSW=0，RWIGS=1.312,POSCAR中的晶格常数定义为优化的晶格常数5.5，其它参数不变。运行VASP，计算完后得到包含了态密度值的DOSCAR文件。采用split_dos对态密度文件DOSCAR进行分割，得到总态密度DOS0以及两个硅原子的分波态密度DOS1、DOS2。将DOS0的数据导入到Origin作图，即得到了总态密度图，如图3所示。图3：总态密度图将DOS1的数据导入到Origin作图，即得到了第一个原子的分波态密度图，如图4所示。黑色的线对应s态电子的分波态密度，红色的线对应p态电子的分波态密度。第二个原子的分波态密度与第一个原子相同。图4：第一个原子的分波态密度图由能带图可知，导带的能量范围大致在0—2.5eV，由上面的分波态密度图可知，此能量范围内既有s电子又有p电子两者比例接近1:1。价带的能量范围大致在-5—0eV，由上面的分波态密度图可知，此能量范围内主要是p电子。查看POSCAR文件，考察各个K点处导带和价带的组分，得到的精确结果与上面的定性分析一致，即价带基本都是p电子，导带既有s电子又有p电子。前面的计算所选用的参数主要是参照侯柱峰的《VASP软件包的使用入门指南》以及师兄的经验数据，并不一定是最优的选择。完成上面的作业后，杨吉辉师兄又建议我做k点数目和k-mesh的优化、切断动能的优化以及SIGMA的优化。（五）k点数目和k-mesh的优化选取不同的k点网格划分，分别进行计算，得到相应的约化k点数和能量值，如下表所示K点网格约化K点数Energy(eV)2*2*22-10.681606083*3*34-10.241025904*4*410-10.838020465*5*510-10.778370316*6*628-10.840599227*7*720-10.831470428*8*8160-10.840706689*9*935-10.8391031210*10*10110-10.84070717以每种网格一个基矢方向上的分割点数N和相应的能量值作图，得到如图5所示的Energy-N关系图。由图可见，在N=6时能量就基本收敛了，故可以选择6*6*6的网格。图5：能量与一个基矢方向上的分割点数N的关系图以约化K点数与能量值作图，如图6所示。由图可见，在K点数为10时能量即基本收敛了。图6：约化K点数和能量的关系图（六）切断动能ENCUT的优化选取不同的切断动能ENCUT，分别进行计算，得到相应的能量值，所得数据如下表所示ENCUT(eV)Energy(eV)60-10.5718490070-10.4814626080-10.49859273100-10.65940511150-10.75226690200-10.80172279250-10.82617381300-10.83910312350-10.84293686作出能量和切断动能ENCUT的关系图，如图7所示。由图可见，在ENCUT=250eV时，能量基本收敛了，故ENCUT可减小到250eV。图7：能量和切断动能ENCUT的关系图（七）SIGMA的优化选取不同的SIGMA，分别进行计算，得到相应的entropyT*S值EENTRO，所得数据如下表所示SIGMAEENTRO0.02-0.000000000.05-0.000000000.1-0.000000010.2-0.000094740.3-0.001082060.4-0.004222810.6-0.020220280.8-0.05656195作出EENTRO和SIGMA的关系图，如图8所示。由图可见，在SIGMA=0.2时，EENTRO就已经非常接近于0了，故SIGMA的值可优化为0.2。图8：EENTRO和SIGMA的关系图