能带计算大概是VASP用得最多的功能之一但我在各个群里看到最多的求助帖也是“能带计算跑出来是乱的”“带隙是0”“K路径又断了”。VASP做能带说难不难说简单里面不少坑确实是文档里不会明说的。这篇文章我按自己平时实际跑能带的流程来写从原理到输入文件再到脚本提取数据一次性把全流程捋清楚尤其适合刚接触DFT、被能带计算绕晕的新手。先说结论VASP能带计算的标准流程是结构优化、自洽计算、非自洽计算、数据提取四步前两步保证体系处于正确的基态第三步在固定自洽电荷密度条件下沿高对称K路径求解本征值第四步把本征值映射成能带图。每一环的设置都会影响最终结果的物理意义下面逐个拆开讲。1. 先弄明白能带计算到底在算什么十个人里有八个第一次跑能带都是拿到流程模板就开算算出图来再对着图发呆。这样也不是不行但一旦结果异常你连该怀疑哪个环节都不知道。所以先把物理图像建立起来。1.1 为什么必须是“三步走”很多人不理解为什么算能带不能像算单点能那样丢一个POSCAR进去直接得出能带图原因是能带图本身就是“能量随K点的变化关系”而K点只是布里渊区里的坐标点体系的实际电荷密度分布、势场分布并不取决于你沿哪条K路径去看它。换句话说电荷密度是K网格的函数不是K路径的函数。于是这里有个计算顺序问题如果你想沿一条一维K路径即高对称线求解各个K点上的本征值那么这些K点上的波函数是依赖于体系实际电荷密度的。而电荷密度需要在一个足够密、能完整覆盖布里渊区的K网格上自洽求解出来。这就是分步的核心逻辑结构优化和自洽计算都使用均匀分布的自动K网格让体系达到基态得到准确的电荷密度和势场非自洽计算读入第一步得到的电荷密度在K路径上逐点求解本征值因为此时不需要重新自洽速度会快得多把路径上每个K点对应的本征值按能量从低到高整理出来就是所谓的“能带”。我自己刚学的时候犯过错误在非自洽这一步也开了自洽结果跑了几十步才收敛K路径上的K点又少出来的带子断断续续。后来理解了流程就是“均匀网格算密度路径网格算本征值”再没出过这类问题。1.2 能带图里藏着的物理信息能带图横轴是倒空间里的K点路径通常是沿布里渊区的高对称点连线从Γ点出发经过X、M、R等点再回到Γ纵轴是能量一般以费米能级为原点E−EF。这条线上每一条曲线就是一个能带反映的是电子在周期性势场中某个本征态的能量随动量的色散关系。从能带图里我们能提取出大量物理信息这也是VASP能带计算最核心的价值所在带隙价带顶与导带底之间的能量差决定材料是金属、半导体还是绝缘体直接带隙与间接带隙价带顶和导带底是否位于同一个K点。如果不在同一个K点就是间接带隙半导体光电应用场景完全不同有效质量能带在带边附近的曲率曲率越大有效质量越小电子的输运性质就越好能带宽度某条能带在路径上的能量展宽反映电子态的离域程度能带简并与分裂轨道杂化、晶体场效应、自旋轨道耦合都会在能带上留下痕迹。我把这一步放在最前面讲是因为整个实操过程的参数选择几乎都是为了能清晰可靠地得到这些信息。理解了“能带图是结果不是目标”你才会自然地理解后面定量设置K点数、能带条数、费米能级归零这些操作。2. 输入文件逐个拆解四件套不是随便填的VASP计算的基本输入文件是INCAR、POSCAR、POTCAR、KPOINTS能带计算也不例外。但能带计算里这四件套每个都有专门的讲究尤其是INCAR和KPOINTS稍微填错一点结果就可能变成一团乱线。2.1 INCAR 参数的“能带专用”配置先给一个典型的能带计算INCAR以非自洽计算为例SYSTEM Band calculation ENCUT 520 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 PREC Accurate LREAL .FALSE. LORBIT 11 ICHARG 11 LWAVE .FALSE. LCHARG .FALSE.逐个说参数含义。结构优化时常用ISMEAR1或Methfessel-Paxton但能带计算的非自洽步骤里体系通常是半导体或绝缘体推荐用ISMEAR0高斯展宽SIGMA取0.05电子伏特即可。金属体系如果遇到能带穿越费米面用ISMEAR1配合更小的SIGMA也是常見做法但非自洽计算里我一般都会用ISMEAR0因为此时已经不需要自洽迭代展宽函数只影响占据数的统计方式。ICHARG是关键参数。非自洽计算必须设置ICHARG11表示读取已有的CHGCAR文件并且在计算过程中保持电荷密度固定不变不再进行自洽迭代。这一步配合前面自洽计算生成的CHGCAR文件是能带计算的“开关”。很多人忘了加ICHARG11系统会默认从初始电荷密度开始重新自洽K路径网格又稀疏最后得到的本征值根本没意义。LREALFALSE保证计算在倒空间中进行精度更高。LORBIT11用于输出每个能带在原子轨道上的投影成分后面可以用这个数据来做“轨道分辨能带图”或者“fat band”。如果你只想要纯能带图LORBIT可以不开但开了也不亏后续画projected band会省掉重新算一遍的麻烦。2.2 KPOINTS能带计算的灵魂所在能带计算的KPOINTS要分两步分别设置。自洽步骤用自动K网格格式是这样的Automatic mesh 0 Gamma 8 8 8 0 0 0这是以Gamma点为中心的8×8×8均匀网格。具体网格密度和体系的晶格常数、原子数目有关一般原则是晶格常数越大需要的K点越少原子数越多K点也可以适当减少。做收敛性测试时把总能量或带隙作为收敛判据网格加密到能量差小于1毫电子伏/原子就可以。非自洽的能带计算KPOINTS则必须用Line-mode格式如下Band path 5 Line-mode 0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.5 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.0第二行的5表示沿5条高对称线计算第四行开始的每一行是一段K路径的起点和终点。VASP会自动在两段之间插入K点K点总数目由每条线上的点数决定。这里关键要弄懂的是每行两个K点的坐标必须落在倒空间高对称点上不同晶系的布里渊区高对称点坐标表可以在Materials Project或者一些VASP教程中找到。VASP在Line-mode下每条线默认插入10个K点可以通过在INCAR里加“NEDOS”间接控制但更直接的是在KPOINTS第五行之后如果写“10”就表示每条线10个点。我一般设20路径长一点的话会设到40K点太疏能带会显得很粗糙太密则浪费计算资源。2.3 POSCAR与POTCAR的准备工作POSCAR是体系结构的描述文件包括晶格常数、晶格矢量、原子种类和坐标。做能带计算时POSCAR必须与结构优化结束后保持一致不能有任何漂移。我在实际项目里遇到过一次诡异的问题非自洽计算时用了自洽计算前未完全优化的POSCAR结果能带在某个高对称点附近出现两条带子交叉后来发现是结构优化中途停止导致的原子位移没收敛。建议你在自洽计算完成后直接复制CONTCAR为POSCAR这样能保证结构完全一致。POTCAR则用与自洽计算完全相同的文件按POSCAR中的原子顺序拼接。VASP对POTCAR顺序很敏感顺序不对会直接报错或者悄悄算错有些版本甚至不报错就那么算下去结果全是错的这类错误特别坑。还有个小细节做能带计算时POSCAR里的原子种类和数量必须和POTCAR里的元素一一对应。我曾经在一次合金体系计算里POSCAR里写的是Fe、Ni两种元素POTCAR拼接时手滑把Ni的位置复制成了Fe的VASP居然也能跑完出来的能带全是铁的铁磁性特征查了半天才定位到是POTCAR出错。3. 完整实操从优化结构到画出能带图下面以闪锌矿结构的砷化镓GaAs为例走一遍全流程。GaAs是直接带隙半导体的经典代表也是DFT能带计算里最常见的入门体系拿来演示最有说服力。3.1 结构优化别在化简之前偷懒结构优化的对象是POSCAR初始结构可以从Materials Project网站或者实验晶体结构数据库里拿。GaAs的闪锌矿结构空间群是F-43m晶格常数约5.65埃原胞包含2个原子Ga和As各一个。初始POSCAR可以写成GaAs zincblende 1.0 0.0 2.825 2.825 2.825 0.0 2.825 2.825 2.825 0.0 Ga As 1 1 direct 0.00 0.00 0.00 0.25 0.25 0.25结构优化的INCAR建议这样设置SYSTEM GaAs optimization ENCUT 520 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 EDIFFG -0.01 IBRION 2 ISIF 3 NSW 100 PREC Accurate LREAL .FALSE.IBRION2表示使用共轭梯度法进行离子弛豫ISIF3表示同时优化原子位置、晶胞体积和晶胞形状。如果是做晶格常数比较确定的块体材料ISIF3是首选但是如果是表面体系、分子体系或者有固定约束的体系就得把ISIF改成其他值或者固定某些原子。EDIFFG-0.01表示离子弛豫的收敛标准是原子受力小于0.01电子伏特每埃。注意负号表示力的单位正值表示能量差。很多人把EDIFFG设成1E-4结果离子步极其缓慢跑几天都完不成其实对于大多数体系-0.01的力收敛标准已经足够精细了。结构优化跑完后先看OUTCAR末尾的“reached required accuracy - stopping structural energy minimisation”这句话确认真的收敛了然后复制CONTCAR为POSCAR进入下一步。3.2 自洽计算拿到可靠的电荷密度把上一步优化好的POSCAR保持不变KPOINTS换成均匀自动网格INCAR改为自洽计算配置SYSTEM GaAs SCF ENCUT 520 ISMEAR 0 SIGMA 0.05 EDIFF 1E-6 PREC Accurate LREAL .FALSE. LORBIT 11 LWAVE .TRUE. LCHARG .TRUE.这一步的关键是LWAVE和LCHARG都为.TRUE.确保波函数和电荷密度文件都能写出来。VASP默认就会产生CHGCAR和WAVECAR但显式写出来更保险也方便你检查它们的存在。自洽迭代的收敛判据是EDIFF1E-6电子伏特即两次电子步之间的能量差小于这个值就停止自洽循环。对于大部分半导体体系这个精度足够了。自洽计算通常几十步就能收敛。跑完后打开OUTCAR搜索“free energy TOTEN”看最后几个数值是否不再明显变化。如果能量在某个值附近来回振荡不收敛多半是ISMEAR和SIGMA设置不对金属体系用了ISMEAR0且SIGMA太小就会振荡半导体体系反而没问题。金属体系自洽推荐ISMEAR1和SIGMA0.1配合使用半导体系用ISMEAR0就好。自洽完成后重要的事情是检查CHGCAR文件是否生成并且文件大小合理。CHGCAR动辄几百兆一切正常后就可以进入非自洽步骤。3.3 非自洽计算沿着K路径取能量本征值现在把POSCAR保持不变KPOINTS换成前面写的Line-mode路径INCAR改成开头那段“能带专用配置”。加上ICHARG11让VASP读入已有的CHGCAR并固定电荷密度。运行VASP时一定要确保当前目录下有CHGCAR文件。我记得有一次跑能带程序卡在自洽循环里迭代了一百多步都不收敛一看是忘了复制CHGCARICHARG11根本读不到文件于是VASP就从初始电荷密度开始重新自洽了K路径上网格又稀疏自然不收敛。这种低级错误特别容易犯建议养成固定套路自洽计算完成后立即写个脚本把CHGCAR复制到能带计算的目录里。非自洽计算会比自洽快很多因为它不需要迭代求解电荷密度。如果这一步跑得很慢检查是不是INCAR里忘了加ICHARG11或者CHGCAR没有被正确读取。VASP在OUTCAR里会明确显示“charge density read from file”搜一下就知道有没有读到。3.4 提取能带数据的实战操作非自洽计算跑完后能带信息存放在EIGENVAL文件里。EIGENVAL的格式很简单前几行是体系信息然后每个K点一个数据块块内每一行对应一条能带的本征值。手动提取能带是体力活工具最简单的是用VASPKIT或者一行Python脚本搞定。VASPKIT提取能带特别方便安装好之后在能带计算目录下输入vaspkit然后在交互菜单里选择“1”进入能带相关功能再选“1”提取能带数据VASPKIT会自动生成BAND.dat文件包含三条数据列第一列是K点在路径上的累积距离第二列开始是各条能带的能量值。然后用gnuplot或者Python画图即可。如果你不习惯交互式工具也可以直接用Python解析EIGENVAL。下面这个脚本是我自己惯用的简化版适合直接复制使用import numpy as np def read_eigenval(filename): with open(filename, r) as f: content f.readlines() nk int(content[1].split()[0]) nbands int(content[1].split()[1]) kpoints [] eigenvalues [] idx 7 for i in range(nk): parts content[idx].split() kpoints.append([float(parts[0]), float(parts[1]), float(parts[2])]) idx 1 band_block [] for j in range(nbands): band_block.append(float(content[idx].split()[1])) idx 1 eigenvalues.append(band_block) return np.array(kpoints), np.array(eigenvalues) kpoints, eigenvalues read_eigenval(EIGENVAL) # 每个K点的累积距离用于横轴 kdist np.zeros(len(kpoints)) for i in range(1, len(kpoints)): kdist[i] kdist[i-1] np.linalg.norm(np.array(kpoints[i]) - np.array(kpoints[i-1])) # 假设费米能级已知从OUTCAR里读 efermi -0.0 # 替换成你从OUTCAR读到的费米能级 eigenvalues - efermi画图部分自己用matplotlib处理横纵轴就行。重点是把OUTCAR里的费米能级读出来常见做法是搜索“E-fermi :”字段然后用每个本征值减去费米能级这样能带图的零点就落在费米能级上。很多人画出来的能带图整体偏移了几个电子伏特往往就是忘了这一补正。4. 能带图怎么看结果解读与常见误区跑完能带之后真正的重头戏是看图说话。能带图一眼看去就是一堆线但怎么从里面准确读出物理量怎么发现可能算错了这些都需要经验。4.1 带隙值怎么读能带图的纵轴能量以费米能级为零点后半导体和绝缘体的能带图会在某个能量区间出现一条水平空白带——这个空白就是带隙。价带顶的能量记为EVBM价带最大值导带底的能量记为ECBM导带最小值带隙就是二者之差。以GaAs为例用PBE泛函算出来的带隙大约在0.5电子伏特左右而实验值是1.42电子伏特。这不算算错而是PBE泛函系统性低估带隙的已知问题背后原因是PBE没有正确描述电子自相互作用。想要得到接近实验值的带隙就要用杂化泛函HSE或者GW方法这部分后面单说。读带隙时还有两个容易犯的错误一是把费米能级附近的部分占据带误认为导带底。自洽计算用的是有限温度展宽费米能级附近可能有一些部分占据的态在能带图上表现为几条能带穿过费米面。对半导体来说你要读的带隙是两条能带之间的最大能量间隔而不是所有能带和费米能级的交点。二是搞混直接带隙和间接带隙。GaAs的价带顶在Γ点导带底也在Γ点这是直接带隙半导体带隙值直接看Γ点附近就能带间距。但硅的价带顶在Γ点导带底在[100]方向靠近X点的地方两个极值不在同一个K点此时带隙是“间接带隙”不能只盯着Γ点看。4.2 高对称字母与路径对应关系的坑横轴上的每一个字母都对应布里渊区里的一个高对称点。简单立方、面心立方、体心立方、六方晶系各自的高对称点名称和坐标都不同找错坐标是能带计算里最常见的低级错误。我见过不少同学做六方体系比如石墨烯或者六方氮化硼KPOINTS里的路径还按立方晶系习惯写Γ-X-W-K结果能带图完全看不出狄拉克锥。石墨烯的布里渊区高对称路径应该是Γ-M-K-Γ其中K点是狄拉克锥所在位置。路径选错物理图像整个就没了。路径选择有一个简单可靠的参考来源Materials Project网站的材料详情页上会给出这个材料的标准能带路径直接抄过来用就行。也可以查一些专门整理高对称点坐标的文献比如Setyawan和Curtarolo在2010年发表的“High-throughput electronic band structure calculations”论文里面有所有布拉伐格子的标准K路径表非常推荐收藏。4.3 能带图常见的几个大坑第一个坑是费米能级归零处理不当。非自洽计算读自洽计算的CHGCAR但非自洽的OUTCAR里也有费米能级信息。自洽计算和非自洽计算由于K点取样不同读取的费米能级可能有细微差别建议统一从自洽计算的OUTCAR里读费米能级保持一致。第二个坑是能带出现“断带”。如果你在能带图的某个高对称点附近看到线突然断掉大概率是KPOINTS里相邻两段高对称路径的终点和起点坐标不连续。比如第一段终点是M点0.5, 0.5, 0.0第二段起点却写成了K点0.333, 0.333, 0.0线段之间出现跳跃图就断了。仔细核对路径坐标能避免这个问题。第三个坑是为了省时间把K点设得太少。每条高对称线只放5个或10个K点能带就是几条折线连成的锯齿看不出真实的色散关系。我的经验是每段路径至少20个K点关键的Γ点附近如果有带边极值那段路径加密到40个点也不为过。第四个坑是投影能带和非投影能带混用。你在INCAR里设置了LORBIT11提取数据时能同时得到total能带和project能带两者加和应该等于总能带。画图时如果不小心把project数据当total来处理能带数目会对不上图形看着也怪。5. 进阶方向HSE、贝里曲率与缺陷态能带能带计算的基本盘掌握之后很多实际课题会逼着你往上走一步。PBE能带的带隙低估问题怎么解决拓扑材料怎么判断缺陷体系怎么做能带这些都是能带计算的延伸场景这里交代一些我的经验。5.1 杂化泛函HSE修正带隙PBE带隙和实验值偏差是DFT的老大难问题。对大部分半导体PBE会低估带隙30%~50%而杂化泛函HSE06通过把一部分精确交换项混入PBE交换项能显著改善带隙预测精度。HSE计算能带的流程和PBE类似但有两个重要区别一是自洽计算时就要用HSE泛函。设置方法是在INCAR里加LHFCALC .TRUE. AEXX 0.25 HFSCREEN 0.2 PRECFOCK FastAEXX0.25是HSE06混入25%精确交换的标准值HFSCREEN0.2是屏蔽参数。二是计算量剧增HSE自洽计算的时间一般是PBE的几十倍。对于超过几十个原子的体系HSE自洽要么等很久要么直接算不动。一个折中的方案是用PBE做结构优化和自洽计算然后用“HSE single shot”模式即用PBE的电荷密度作为初值只算几步HSE得到修正后的本征值。做法是在PBE自洽完成后用上述HSE参数再跑一次非自洽计算同时保留ICHARG11。这个方案的计算量比完整HSE自洽小很多带隙修正效果也能大致接近。但注意它不是一个完全自洽的HSE结果严谨来说还是得做HSE自洽具体取舍看课题需要。5.2 贝里曲率与拓扑性质计算热词里有人搜“vasp贝里曲率计算”这确实是当前凝聚态计算的热门方向。VASP从5.4版本开始内置了计算贝里曲率和反常霍尔电导的功能路径是先用PBE或者HSE做自洽计算然后利用WannierTools或Wannier90拟合瓦尼尔函数再用紧束缚模型计算贝里曲率。VASP里跑瓦尼尔拟合需要在INCAR里设置相关参数并且KPOINTS要选一个均匀网格。贝里曲率计算的流程比普通能带复杂很多但只要基础能带流程熟练转换过去并不难。核心是先得到准确的能带结构在这个前提下再做瓦尼尔插值得到的拓扑性质才有意义。5.3 缺陷体系的能带计算注意事项热词里还有“vasp缺陷浓度”这也是能带计算的延伸应用。做缺陷体系能带时你通常要构建一个超胞比如3×3×3的原胞扩展然后移除一个原子或者替换一个原子。此时能带图会变得非常密集因为超胞的布里渊区缩小了能带发生了折叠。缺陷体系的能带解读要特别小心缺陷态会在带隙中间出现几条孤立的能带它们往往比较平直色散很小对应局域化的缺陷态。如果你关心缺陷态能级的位置和占据情况最稳妥的做法是同时画出总态密度和投影态密度看看这些带隙中间的能带有多少成分来自缺陷原子附近的轨道。单独看能带图很容易漏掉或者误解这些信息。另外超胞的尺寸直接影响缺陷态之间的相互作用。超胞太小时周期性镜像缺陷之间会有耦合导致缺陷能级展宽看起来像连续带而不像分立态。一般建议至少用4×4×4的超胞具体收敛性要测试。6. 常见问题排查速查表最后整理一个能带计算的高频问题速查表这些都是我在实际项目里踩过或者帮别人排过的雷按“问题→可能原因→解决方案”的格式整理成表格方便大家直接对照现象可能原因解决方案非自洽计算迭代超多步不收敛忘了ICHARG11或CHGCAR没复制过来设置ICHARG11确认目录下有CHGCAR能带图在某个K点突然断开KPOINTS相邻路径的端点不连续核对高对称坐标是否连续能带图全是锯齿状没有平滑色散每条K路径的K点数太少Line-mode下每条线增加到20~40个点带隙明显偏低PBE已知局限泛函本身的缺陷用HSE或GW修正带隙能带整体能量偏移没有把费米能级归零从OUTCAR读E-fermi所有本征值减去它能带图里有额外交叉线结构未完全优化或高对称路径不对确认使用结构优化后的CONTCAR核对K路径金属体系自洽不收敛ISMEAR0配合过小SIGMA改用ISMEAR1或MP展宽SIGMA取0.1能带条数和预期不符自旋极化打开或关闭设置错误检查ISPIN参数磁性体系要ISPIN2提取数据时EIGENVAL解析失败文件版本不同格式偏移用VASPKIT提取避免手写解析脚本出错投影能带图数据不完整LORBIT没设置或没开能带计算前在INCAR里加LORBIT11几个额外的独家避坑tips都是常规教程不写的经验算完结构优化后先用同样的K网格做一次自洽不要直接跳到Line-mode。这样能确保电荷密度确实是收敛结构的电荷密度而不是零号步电荷密度。能带计算前先看一眼OUTCAR里的费米能级和自洽计算结果对比如果差别超过0.1电子伏特多半是结构变了或者CHGCAR没读对别急着提取数据。测试K点收敛性的时候把带隙值也作为收敛指标之一只盯总能量有可能错过带隙还没有完全稳定的时候。计算表面体系能带时POSCAR里必须加足够厚的真空层一般至少15埃并检查能带图里有没有出现虚假的平带那是真空层中spurious态的标志。我个人在实际操作中最深的一点体会是能带计算的问题十有八九不是VASP算错了而是输入设置和数据处理对不对的问题。流程熟练之后一次能带计算从准备好结构到出图通常半小时内就能搞定。但前提是你真的理解每一步在做什么而不是机械地套模板。希望这篇全流程梳理能帮你少走一些弯路把时间花在解读物理而不是debug上。