1. 为什么掺了个杂能带图就“变质”了你大概率遇过这种场景拿初胞算了一张干干净净的能带图然后在4×4×1超胞里替位掺杂掺进一个杂质原子再跑一遍能带出来的图密密麻麻、每条线都在“颤”。那一瞬间你可能会怀疑是自己建模建错了或者干脆怀疑VASP是不是在乱给结果——其实两者都没错。真正的问题在于超胞计算天然会把能带折叠到一起你直接画出来的超胞能带并不是我们习惯讨论的初胞色散关系。这时候就需要把超胞能带“展开”回初胞布里渊区也就是标题里那个unfold。这篇东西我主要写给两类人一类是做掺杂、缺陷、合金等体系的研究生最近被超胞能带图折磨到想弃坑另一类是已经会跑VASP基础流程但还没搞清楚unfold到底在数学上做了什么、工具又是怎么工作的进阶用户。我会用自己实际跑的替位掺杂超胞流程做主线把建模、输入文件、展开原理、工具选型和常见坑一次说清楚。相信我把这些理顺以后你再面对掺杂体系的能带图就不会觉得心虚了。1.1 周期镜像为什么替位掺杂绕不开超胞先说替位掺杂为什么要拉超胞。DFT的电子结构计算建立在周期性边界条件之上平面波基组天然假设体系在三个方向无限重复。如果你直接在初胞里把一个原子换成掺杂原子那等效于每个晶胞里都有一个杂质掺杂浓度直接就是100%。这在物理上完全不是我们想要的“稀掺杂”场景。比如做半导体里百万分之一量级的掺杂或者石墨烯里百分之几的氮替位碳就必须把原胞放大成超胞然后在超胞的几十个原子里替换其中一个。超胞尺寸的选择本质上是掺杂浓度和计算量之间的权衡。你用4×4×1的石墨烯超胞32个原子里替换1个掺杂浓度大约是3.1%用3×3×1超胞则是1/18≈5.6%。浓度越低杂质之间的周期镜像距离越大相互影响越小但这意味着计算量按原子数快速上涨。另一个隐藏问题是镜像相互作用超胞边界处的杂质和相邻超胞的杂质会“隔空对话”如果超胞太小杂质能级的色散就会被人为加宽展开之后的能带也会出现假象。我的经验是如果是定性看能带趋势2×2或3×3超胞通常够用如果要和光谱实验定量对比尽量把最小镜像间距拉到10 Å以上。1.2 能带折叠是超胞的“副作用”超胞把实空间放大代价就是倒空间缩小。原本初胞布里渊区里的一条能带在超胞布里渊区里会被折叠成好几段。打个比方你在大房间里听一个人说话声源位置清楚现在把房间隔成很多小间每个小间里听到的都是混在一起的回声。超胞能带图上那些密密麻麻的线本质上就是初胞能带被“叠”进更小布里渊区的结果。这种折叠不是VASP算错了而是周期性和倒空间映射的必然。可问题是我们做掺杂研究的最终目标通常是要回答掺杂之后原本的能带结构哪些地方变了杂质引入的能级出现在哪个位置费米面附近是不是多了一条平带这些问题都必须回到初胞布里渊区的语言里才方便讨论。直接画超胞能带图大家只会看到一团紊乱的线既没法标出对称点也没法拿去做ARPES谱的直观对照。展开unfold这一步就是把折叠产生的伪影去掉把真正属于初胞的那条色散关系“提取”出来。2. 替位掺杂建模从初胞到超胞的输入文件准备2.1 替位位点怎么选替位掺杂的意思是掺杂原子占据原晶格里某个原子的位置而不是挤进间隙。建模第一步是决定替换谁。这里有两个原则一是看化学合理性替位原子的尺寸、电负性、价电子数和被替换原子越接近形成能越低体系越稳定二是看对称性等价位置晶体里可能有好几个同种原子位点替换其中任何一个在物理学上等价所以通常只需要算一个位点。不过更严谨的做法是并行算几个候选位点的形成能公式大致是E_f E_doped - E_pure - μ_dopant μ_host这个形成能可以帮你判断替位掺杂在热力学上是否容易发生。实际工作中你会发现替位位点的选择会直接影响展开后的权重分布尤其是在低对称性掺杂位点附近能带会混入更多杂质成分展开权重普遍下降。如果你发现自己展开出来的“真”能带权重不到0.7先别急着怀疑工具回头看看是不是掺杂位点本身破坏了太多对称性。2.2 超胞构建与原子替换实操我习惯用VASPKIT处理超胞和K点用ASE做原子替换两个工具都很顺手。先用VASPKIT把初胞扩成超胞vaspkit -task 501按提示输入扩胞倍数比如4 4 1。生成的新POSCAR是4×4×1超胞里面所有原子都是原来的元素。接下来用ASE或者手动改坐标的方式把其中一个原子替换成掺杂元素from ase.io import read, write atoms read(POSCAR_super) # 确认要替换的原子序号例如超胞里所有原子中第17号 atoms[16].symbol N write(POSCAR_doped, atoms, vasp5True, directTrue)这里有个细节容易踩坑POSCAR里面的原子顺序必须和POTCAR里元素的顺序严格对应。VASP是先把所有势文件拼接成一个POTCAR然后按照POSCAR里元素的排列顺序去匹配。比如你做石墨烯掺氮POSCAR的元素行应该是“C N”那么你的POTCAR就必须先cat C的POTCAR再cat N的POTCAR顺序不能反。如果你用ASE改了原子符号ASE输出的POSCAR会把元素自动分组排序你一定要重新核对该顺序否则计算会直接报错或者张冠李戴。2.3 INCAR、KPOINTS、POTCAR的关键设置一个完整的掺杂超胞流程会经历三个阶段结构弛豫、静态自洽、能带计算。三个阶段的INCAR关键参数是不一样的我按自己常用的模板列一下。结构弛豫阶段SYSTEM doped_supercell ENCUT 500 PREC Accurate EDIFF 1E-6 IBRION 2 ISIF 2 NSW 100 ISMEAR 0 SIGMA 0.05 LREAL .FALSE. ISPIN 2 MAGMOM 31*0.0 1.0 LORBIT 11静态自洽阶段IBRION -1 ISIF 3 NSW 0 ICHARG 1 ISMEAR -5 LREAL .FALSE. LWAVE .TRUE. LCHARG .TRUE. LORBIT 11能带计算阶段ICHARG 11 ISMEAR 0 SIGMA 0.05 LORBIT 11 LWAVE .TRUE. LREAL .FALSE.先解释一下ISIF2而不是ISIF3替位掺杂超胞如果允许晶胞形状和体积自由变化很容易因为局部张力产生过大形变而且这种形变未必对应真实实验条件下的结构。所以我的习惯是先固定超胞形状只让原子位置弛豫。如果确实需要放晶格那就在后期单独做体积弛豫。LREAL.FALSE.也是专门为展开准备的。展开工具需要平面波系数来做投影如果你用了LREAL.TRUE.波函数在实空间近似处理投影权重会失真。很多人在这一步翻车展开出来的权重全都小得没法看。MAGMOM的设置要看体系。掺氮、掺硼这类体系通常带一个未配对电子开自旋极化ISPIN2以后要把掺杂原子附近的初始磁矩给定好。上面写了31个0和一个1意思是32个原子中最后一个也就是掺杂原子初始磁矩设为1.0。提示无论哪个阶段都要保证LWAVE.TRUE.展开需要保留波函数文件LORBIT11则保证PROCAR里有足够的投影信息。3. 能带展开到底在展开什么折叠机制与谱权重3.1 超胞布里渊区与初胞布里渊区的对应关系想理解展开先得搞清楚实空间超胞和倒空间超胞的对应。假设你把初胞沿a、b、c三个方向分别放大N_a、N_b、N_c倍那么超胞布里渊区就是初胞布里渊区的“缩小版”超胞BZ的体积只有初胞BZ的1/(N_a N_b N_c)。反过来初胞BZ中的一个点放到超胞BZ中会对应到一系列等价的“反折叠”点。超胞BZ里任意一个K点都能通过加上初胞的倒格矢G映射回初胞BZ中的一组k点。这个折叠的倍数不是随便定的它正好等于超胞里包含的初胞数量。比如4×4×1超胞包含了16个初胞那么超胞能带图上每一个K点的每一条能带理论上都是由初胞BZ中16个不同k点的能带折叠而来的。你看到能带图密密麻麻本质上就是这16倍的线叠在一起。3.2 展开系数与谱权重展开的数学目标是给出每个超胞本征态在初胞布洛赫态上的投影大小。常用的谱函数可以写成下面这种形式A_s(k, ε) Σ_{K,n} |ψ_{k,i}|ψ_{K,n}|² δ(ε - ε_{K,n})这里ψ_{K,n}是超胞计算得到的本征态ψ_{k,i}是与超胞K点对应的初胞布洛赫态投影系数的平方就是所谓的展开权重unfolding weight取值范围在0到1之间。权重越接近1说明这条超胞能带越“纯粹”它就是初胞能带的延续权重明显小于1说明这条能带混入了杂质态或者其他能带的成分通常表现为影子带shadow band。δ函数把能量锁定在超胞本征能量附近所以展开结果是一张三维图横轴是初胞k路径纵轴是能量颜色深浅代表权重。绝大多数后处理工具会把这个三维信息压缩成带权重的能带图权重低的线画淡一点权重高的线画实一点。这里的核心思想是超胞计算包含了完整的信息只是信息的坐标系被折叠了。展开并不是“算出”新能带而是把我们本来就知道的初胞色散成分从折叠坐标系里还原出来。3.3 权重低不代表算错很多人第一次跑完展开看到有些权重只有0.2、0.3的能带第一反应是工具坏了。其实权重低有物理原因。替位掺杂破坏了原有的平移对称性和点对称性掺杂原子附近的状态会和宿主能带杂化产生局域态和共振态。这些态投影到初胞布洛赫态上本来就不可能有接近1的权重。它们代表的物理是真实的杂质态不等于宿主能带你非要把它“展开”成初胞能带当然只能得到一部分投影。所以看展开图的时候我一般只把权重高于0.5的线当成“宿主能带的变化”把权重低但连续平缓的带当成杂质或共振信号。这种区分对后续分析特别重要能帮你在掺杂体系的能带图里分清哪些是杂质自己的态哪些是宿主态的重新杂化。4. 从自洽到展开的完整实操流程4.1 第一步结构优化结构优化是一切的前提。替位掺杂后掺杂原子和周围原子间的键长会重新调整如果直接拿初胞结构硬算力的收敛会很差展开出来的虚峰也多。我用的是上一节那个ISIF2的INCAR让原子坐标充分弛豫同时保持超胞格子不变。弛豫完成后检查CONTCAR里的原子位移。以石墨烯掺氮为例C-N键长会比原来C-C键长短一些邻近原子的位置会有几个百分点的偏移。这个偏移量看起来很小但对展开权重影响很大因为初胞和超胞的原子坐标要尽可能保持一致映射关系才清晰。注意如果你在展开阶段发现权重普遍低于0.5先回去检查结构弛豫后的晶格畸变程度这个比工具本身的bug概率高得多。4.2 第二步静态自洽计算静态计算的作用是拿到收敛的电荷密度和波函数。收敛标准要设严一点我的习惯是EDIFF1E-6因为展开权重对能量本征矢的相位极其敏感能量都不收敛权重自然飘。这一步的K点密度也要足够。超胞本身已经很大K点可以相应减少但不要牺牲到让自洽都收敛不了。4×4×1石墨烯超胞我一般用Gamma-centered的4×4×1网格3D体系则根据超胞大小做缩放保证有效K点密度和初胞12×12×12相当就好。静态计算结束后目录里必须有WAVECAR、CHGCAR、OUTCAR和PROCAR。PROCAR是VASP自动输出的但前提是你开了LORBIT11这个在前面已经强调过了。4.3 第三步沿初胞高对称路径的能带计算能带计算的主流做法是用ICHARG11读取自洽的CHGCAR然后沿着K路径做非自洽计算。这里最关键的一点是K路径必须按初胞的高对称点来写而不是按超胞自己的高对称点来写。以二维六角体系为例初胞路径通常是Γ-M-K-ΓKPOINTS文件可以这样写K-Path for Unfolded Bands line-mode c 21 0.000 0.000 0.000 ! Gamma 0.333 0.333 0.000 ! M 0.333 0.333 0.000 ! M 0.500 0.000 0.000 ! K 0.500 0.000 0.000 ! K 0.000 0.000 0.000 ! Gamma这些坐标是初胞布里渊区里的分数坐标。超胞能带计算也可以接受它们因为超胞布里渊区不过是初胞BZ的一个子区域这些点会自然落到对应的折叠位置上。关键是展开工具需要知道你的K点坐标是用初胞晶格定义的这样才能正确做折叠映射。不同工具对此的处理方式不太一样PyProcar一般直接吃KPOINTSBandUP则通常要求你额外提供初胞POSCAR。4.4 第四步执行unfold与绘图展开这一步我实际用下来最顺手的是BandUP和PyProcar两套方案。BandUP的流程是先做预处理把初胞和超胞的原子对应关系、折叠映射准备好再运行主程序。大致命令长这样BandUP_preprocess -p POSCAR_prim -s POSCAR_super bandup具体的参数选项不同版本略有差异但核心输入总是初胞POSCAR、超胞POSCAR、超胞的OUTCAR和PROCAR。BandUP的优势是处理严谨适合对权重精度有要求的体系。PyProcar则是完全脚本化的路子Python环境里一条命令就能跑完import pyprocar # 展开接口传入超胞POSCAR、初胞POSCAR、超胞PROCAR和费米能 pyprocar.unfold( modevasp, poscarPOSCAR_super, poscar_primPOSCAR_prim, procarPROCAR, fermi0.1234, # 费米能由你的OUTCAR读出 ispin2, kpath[...] # 初胞高对称路径 )跑完展开你会得到一个包含(k点、能量、权重)的数据块剩下的就是画图。我一般用matplotlib画散点图权重做透明度映射import matplotlib.pyplot as plt import numpy as np # data: 每行是 k, energy, weight data np.loadtxt(unfolded_bands.txt) k data[:, 0] e data[:, 1] w data[:, 2] plt.scatter(k, e, cw, cmaprainbow, s4, vmin0, vmax1) plt.xlabel(k path) plt.ylabel(E (eV)) plt.ylim(-4, 2) plt.colorbar(labelunfolding weight) plt.show()画出来的图里颜色鲜艳的是真正的宿主能带暗色的就是影子带。我建议把费米能放到能量零点这样看杂质能级位置最直观。5. 工具选择BandUP、PyProcar和我自己的取舍5.1 BandUP严谨但准备步骤多BandUP是专门做能带展开的老牌工具算法上处理对称性和反折叠映射最规范。它会把初胞原子映射到超胞里逐一对齐原子位置再计算投影权重逻辑很透明。如果你做的是低对称性掺杂、或者需要跟ARPES谱严格对比BandUP的置信度会更高。代价是准备流程比较繁琐除了VASP的输出你还要单独准备无掺杂初胞的POSCAR确保初胞和超胞的原子坐标一致编译环境也要花点时间调。而且它本身不自带画图展开结果还需要自己处理成散点图。5.2 PyProcar快速又省心PyProcar更适合日常跑量和快速排查。它内置了K路径生成、费米能读取、谱函数绘图全套接口一个脚本从PROCAR到能带图全搞定。我平时先用PyProcar快速看趋势如果发现某个体系展开结果异常再用BandUP复核。PyProcar的缺点是它对超胞波函数的处理做了不少近似大超胞和多原子体系下内存占用会比较夸张。此外它的文档在不同版本之间变动不小很多老教程里的函数名已经废弃了用之前最好看一眼当前版本的API。5.3 我的选择基准我整理了一个很简单的决策表每逢新体系就按这个来选场景推荐工具理由快速定性看掺杂能带趋势PyProcar配置少、出图快、脚本化方便批量处理和实验谱定量对比BandUP权重计算更严格映射逻辑更可靠对称性较低、权重异常BandUP PyProcar复核两个工具交叉验证避免工具bug误判超胞特别大、内存吃紧BandUPPyProcar在大体系下读PROCAR容易爆内存工具只是手段理解展开原理才是关键。无论用哪个工具最后都要回到谱权重去解读物理哪些是宿主能带的变化哪些是杂质态。6. 我踩过的坑与调过的参展开这个流程看起来简单实际上坑非常密集。我把自己踩过的几个典型坑和排查思路写出来希望你们能少走弯路。6.1 展开权重普遍偏低先查原子映射我第一次跑BandUP时展开权重几乎全在0.5以下所有能带都模糊不清。排查了很久发现不是工具问题而是初胞POSCAR和超胞POSCAR里的原子排列顺序没有对齐。BandUP做展开需要知道初胞里每个原子在超胞里对应哪个位置如果初胞的晶格原点、原子顺序和超胞不一致投影重叠很小权重自然低。解决方法是先把初胞和超胞都用同样的格式规范化一遍用ASE的match_cells或者手动把初胞晶格平移到超胞包含的位置上确保坐标一致。还有一种常见情况是结构弛豫后掺杂原子附近畸变太大初胞的刚性晶格结构和超胞实际结构差得太多这种“物理错位”不会被工具修正只能接受它是体系本身的真实特征。6.2 能量对齐是个隐形炸弹VASP输出的本征能量没有绝对零点展开结果是否准确很大程度上取决于你把超胞能带和初胞能带放在同一个能量尺度上比较。如果超胞计算用的是掺杂体系的势场杂质会让整体静电势漂移能带整体平移零点几个电子伏特是常见的事。我的做法是做一个无掺杂超胞的对照计算然后把掺杂超胞的深能级态比如半芯态能量对齐到无掺杂超胞的对应能级上。PyProcar里有费米能参数可以微调BandUP也有能量对齐的输入项。这一步偷懒展开出来的能带位置可能整体偏了后续分析全错。6.3 LREAL和自旋极化的连锁反应前面说了展开要求LREAL.FALSE.但大超胞用LREAL.FALSE.计算会明显变慢。有些同学贪快在能带计算时开了LREAL.TRUE.结果展开权重变成一堆噪声完全没法用。你要是算很大的体系那就提前规划好预算别在这一步省时间。自旋极化也要小心。掺杂体系开ISPIN2之后能带会劈裂成自旋上、自旋下两套。展开时如果你忘了指定自旋通道工具默认只展开其中一个另一半被丢掉看上去能带就残缺了。我的习惯是静态和能带阶段都显式设置MAGMOM并且在展开脚本里检查is_spin的通道数。6.4 K路径坐标定义错导致的“断带”还有一次我把KPOINTS里的高对称坐标写成了超胞坐标展开出来的能带图断成好几截完全对不上初胞的对称性。原因是超胞BZ的高对称点和初胞BZ的高对称点坐标不同你沿超胞路径跑了一遍再按初胞路径去展开映射自然错乱。解决方式很直接KPOINTS只用初胞高对称路径坐标并且在展开工具里确认一遍输入的初胞POSCAR是没掺过杂的初胞结构。这样工具才知道你要还原到哪个布里渊区。6.5 我最后固定下来的一套流程踩完这些坑以后我现在每做一个掺杂超胞展开都固定走这套流程初胞结构先做一遍无掺杂的能带计算作为对照组。按目标掺杂浓度选择超胞倍数用ASE建超胞、替换原子。结构弛豫ISIF2、LREAL.FALSE.收尾时检查原子位移。静态自洽EDIFF1E-6、ISMEAR-5确认WAVECAR和CHGCAR齐全。沿初胞高对称路径做能带计算ICHARG11、LORBIT11、LWAVE.TRUE.。用PyProcar快速展开再用BandUP复核关键体系。能量对齐用无掺杂深能级做基准画图时权重阈值设在0.5。这套流程我大概是从去年开始固定的之后大部分替位掺杂体系的能带分析都顺了很多。尤其是石墨烯、过渡金属硫族化物这类二维材料展开之后看杂质能级和宿主带的杂化关系非常直观。你如果刚把这些概念串起来建议先挑一个简单体系完整跑一遍把每一阶段输出文件都翻一遍之后换成复杂的掺杂体系就稳了。