1. 项目概述与核心思路
1.1 为什么需要替位掺杂和能带折叠
做第一性原理计算的人，尤其是用VASP做半导体材料研究的，几乎都会撞上同一堵墙：超胞。纯原胞计算固然简单高效，但掺杂问题绕不开超胞。以替位掺杂为例，你要在晶格中替换掉一个原子，如果只用原胞，掺杂浓度动辄25%甚至50%，这跟实际实验条件下的掺杂浓度差了十万八千里。实验里掺杂浓度往往在百分之几甚至千分之几的量级，所以就必须把原胞扩展成2x2x2、3x3x3甚至更大的超胞，让被替换的那个原子所占的比例降到可接受的范围。
但超胞一建立，另一个麻烦立刻出现：能带折叠。晶体周期性变了，第一布里渊区缩小，原本在原胞布里渊区里的能带会被折回小布里渊区里，能带图变得一团乱麻，直接看超胞的能带根本无法判断杂质能级的真实位置、带边变化趋势这些关键信息。很多人一开始不知道这个问题，拿着折叠后的能带图硬分析，结论自然是错的。就这一条，足以让一个月的计算白做。
UNFOLD技术解决的就是这件事。它的核心逻辑是把超胞能带“展开”回原胞布里渊区，利用超胞中每个原子的相位信息还原出它对应的原胞k点上的谱权重，最终得到一张能反映超胞计算真实物理信息的有效能带图。掺杂后杂质能级有没有进入带隙、是否形成局域态、带边是否出现明显的扰动，全都要靠unfold之后的图来看。
1.2 这套流程适合谁
这套“自用vasp单替代掺杂超胞unfold计算流程”不是论文里的套话，是我自己在实际项目中反复用过的完整流程。适合要做以下事情的人参考：
- 研究掺杂对半导体带隙、带边位置的影响，比如找浅能级掺杂、深能级缺陷等。
- 做合金或者固溶体体系，需要还原有效能带结构。
- 分析超胞计算中因为能带折叠产生的“假带隙”或者“伪能带”，想得到干净的、与原胞可对照的band structure。
整个过程大致分四步：准备结构文件 → 自洽计算 → 非自洽能带计算 → 能带展开处理。接下来我按实际操作顺序拆开讲，每一步会说明背后的原理、命令行参数怎么选、容易踩哪些坑。
2. 建模细节与超胞构建实操
2.1 从原胞出发构建替位掺杂超胞
构建超胞最稳妥的方式，是先在原胞上做文章。很多人一上来就手动台账式地复制原子，不仅效率低，而且对称性容易搞错。我习惯用VESTA或者ASE来完成这一步。
以ASE为例，先读取原胞结构，然后做一个规则超胞：
from ase.io import read from ase.build import supercell atoms = read('POSCAR_prim') # 构建2x2x2超胞 atoms_sc = supercell(atoms, [2, 2, 2]) atoms_sc.write('POSCAR_sc', format='vasp')超胞尺寸的选择要讲策略。不是越大越好，而是要看掺杂浓度和计算成本之间的平衡。一般来说：
- 2x2x2的超胞对应掺杂浓度约12.5%（每8个原子替换1个），适合做定性判断。
- 3x3x3的超胞掺杂浓度约3.7%，更接近实验场景，但原子数会增加到几十甚至上百个，计算成本显著上升。
- 对于某些特殊体系，比如层状材料，可能还需要考虑单层内的扩展方式，不能简单做三维超胞。
我的建议是，先做一个中等级别的超胞做快速测试，确认模型没有问题、物理趋势符合预期后，再上大超胞做精细计算。一上来就用3x3x3甚至4x4x4，一旦发现初始结构就有问题，返工成本非常吓人。
2.2 替位原子的选择与结构优化注意事项
替位掺杂的核心动作是“替换哪个原子、用什么原子替换”。这一步看起来简单，实际操作中有一个非常需要留意的地方：缺陷位置的选择。
在一个超胞里，同一个元素可能有多个不等价位置。比如在化合物半导体中，A位和B位的化学环境完全不同，替换A位和替换B位得到的掺杂性质可能天差地别。另一方面，如果超胞足够大，同种原子还有可能因为对称性降低而产生微小的环境差异。不考虑这些直接替换，结果很可能不具备代表性。
替换之后的结构需要做结构优化，但优化之前先做一个简单的检查：用VESTA打开新的POSCAR，确认替换原子周围没有重叠原子，键长是否合理。另外，掺杂后体系带上了额外的电子或失去电子，计算里要对应地调整NELECT或者用背景电荷中和，这一点做缺陷计算的老手都有体会，但新手特别容易忽略。
经验之谈：替位掺杂后的初始结构，最好先做ISIF=2的离子位置优化（只优化原子坐标，不变晶格），避免因为原子受力过大导致优化崩溃。对于晶格常数变化较大的体系，再考虑ISIF=3（同时优化原子坐标和晶格参数），也很有必要。
2.3 POSCAR准备中容易被忽略的原子顺序问题
POSCAR里原子顺序不是小事，直接决定VASP怎么读你的结构。尤其是掺杂后原子种类发生了变化，POSCAR第二行的元素顺序必须和坐标块顺序严格一致，否则VASP会按错误的顺序分配元素类型，计算结果完全不可用。
此外，推荐把被替换的杂质原子放在坐标块的最后，并顺手在文档里记录一下它是第几个原子。这听起来像废话，但等后面要分析局域态、做投影能带的时候，你会发现一个准确的原子序号能救命。VASPKIT、pymatgen这些工具在做unfold时往往需要指定特定原子或特定元素，坐标顺序混乱会导致后续处理脚本全部报错。
3. 自洽计算与非自洽计算的关键参数
3.1 INCAR参数配置：从自洽到非自洽
自洽计算没有太多花活，收敛标准往严了设就行。我常用的参数组合大概是这样的：
自洽INCAR：
SYSTEM = doped_supercell ENCUT = 520 EDIFF = 1E-6 EDIFFG = -0.01 ISMEAR = -5 LORBIT = 11 PREC = Accurate LREAL = Auto ISPIN = 2 NELECT = 特殊体系需要调整 # 掺杂后磁性原子或局域态需要考虑ISPIN=2其中ISMEAR=-5（四面体方法）适合半导体和绝缘体，金属体系要换成ISMEAR=1并配合小的smearing宽度。掺杂体系里如果引入了局域磁矩，一定要开ISPIN=2，让电子自旋自由优化。很多过渡金属掺杂体系，自旋极化打开和关掉，能带结构差异非常大。
结构优化完成之后，非自洽计算直接在优化好的结构基础上做。INCAR变化点主要是：
ICHARG = 11 # 读取自洽计算的CHGCAR，不做电荷自洽迭代 ISMEAR = 0 SIGMA = 0.05这里有一个关键点：非自洽能带计算时，K点路径必须用原胞的能带路径，而不能用超胞的。很多人在这里犯迷糊，认为计算的是超胞，K点也要按超胞来走，结果做出来的band structure根本没法用。
为什么？因为展开unfold的目的就是把超胞能带回映射回原胞布里渊区。如果你选取的超胞的高对称路径本身就不是原胞路径的子集，映射关系就会非常混乱。实际操作中，应该在原胞的高对称k点路径上取点，然后根据超胞的倒格矢关系把原胞k点换算成超胞可接收的k点坐标。大多数时候，直接用原胞路径的分数坐标作为超胞计算的k点路径就能work，但前提是超胞的倒空间能覆盖这些点。
3.2 KPOINTS文件准备：自动生成还是手写
对于自洽计算，K点网格使用自动生成即可，密度要高一些。常见做法是用KSPACING控制自动密度，或者手动指定Monkhorst-Pack网格：
Automatic mesh 0 Gamma 8 8 8 0 0 0超胞变大后，倒空间的k点密度需求其实变小了，所以2x2x2的超胞用4x4x4的k点网格通常就够了，甚至3x3x3超胞用3x3x3都行。关键是测试收敛，你可以在固定其他参数的前提下测试k点从3到6的变化，看总能和带边的变化是否小于0.01 eV。
非自洽计算的KPOINTS要沿着能带路径走，通常用line-mode：
Band path 20 Line-mode rec 0.000 0.000 0.000 Gamma 0.500 0.000 0.000 M ...注意这里我写的是高对称点的坐标，实际路径要根据材料体系来定。常见的做法是用VASPKIT或者pymatgen的bandstructure模块生成标准路径。另外，路径上每个高对称段之间的取样点数（示例中的20）决定了能带曲线的平滑程度，通常20~30个点就够画出一张不错的图。如果你想做更精细的分析，可以增加到40，但计算时间也会翻倍。
这里有一个我踩过很多次的坑：超胞计算的能带路径，如果直接套用原胞路径的分数坐标，会出现某些k点不在超胞倒格矢整数倍网格上的问题。VASP允许任意k点坐标，但后面的unfold程序会要求每个k点有对应的投影信息，k点坐标不在预期网格上时，展开公式里的权重计算容易出错。稳妥的做法是先用原胞路径坐标计算一遍，确认能带展开后的高对称点位置与原始材料一致，再做精细调整。
4. 能带展开的核心工具与计算原理
4.1 主流unfold工具怎么选
我常用的是VASPKIT自带的能带展开功能，以及开源的BandUP程序。两个工具各有优劣。
VASPKIT的unfold模块操作简单，功能集成在程序里，不需要额外写脚本，适合快速出图。BandUP则是专门做能带展开的专业工具，支持投影权重计算，结果更精细，但需要多一步准备输入文件。两者都需要你提供超胞计算得到的WAVECAR或PROCAR文件，区别在于处理方式。
我个人的习惯是，先在VASPKIT里快速做一次unfold，看看整体趋势是否正确，如果只是发文章需要配图、或者做趋势判断，这已经足够。若涉及到复杂的轨道投影分析、边带色散关系的研究，再上BandUP做完整展开。
4.2 展开原理简述
展开原理可以用一个简单的类比来解释。你建超胞，相当于把原胞的倒空间折叠了N倍，原本在第一布里渊区内的一条能带，会被折叠成N条不同的分支。超胞能带图里那些交叉、简并，很多不是真实的物理行为，而是折叠的几何效应。
Unfold要做的事情，就是给每条超胞能带赋予一个“原胞对应k点上的概率权重”。这个权重来源于超胞波函数按原胞波函数展开的投影系数。权重越大，说明这条能带在展开后的原胞图上越重要；权重很小，就是折叠带来的伪影。最终画出来的有效能带，颜色深浅表示谱函数权重，能直观显示杂质能级、能带色散和带边位置。
用BandUP做展开时，输入文件需要提供原胞和超胞的晶格矩阵、原子坐标映射关系，以及超胞k点路径上各k点对应的原胞k点坐标。程序会自动计算匹配关系，输出包含权重的能带数据。
4.3 从WAVECAR到有效能带：完整操作序列
这里给出BandUP的全流程，直接照着操作即可：
准备三个文件：
- 原胞的POSCAR（primitive cell）
- 超胞的POSCAR（supercell）
- 超胞计算得到的WAVECAR和OUTCAR，以及能带计算的KPOINTS
第一步，用BandUP中的bandup命令计算展开权重。需要输入超胞结构、原胞结构以及计算得到的能带信息。对于非自洽计算，这些信息在OUTCAR和KPOINTS里都能找到。
第二步，用plot_band工具绘制展开能带。这个工具读入权重文件，输出权重加权的能带图。输出的数据可以直接用gnuplot或者Python的matplotlib画图，我更喜欢后者，控制力更强。
第三步，把展开后的能带图跟原胞未掺杂的能带图放在一起对比。这是整个流程中最有信息量的一步，能带展开的价值在这里展现无遗。
如果你用VASPKIT，会更简单些。在能带计算完成后，直接运行VASPKIT的unfold相关选项，程序会提示输入必要参数，包括原胞和超胞的对应关系，输出一个可以直接出图的数据文件。整个过程不涉及复杂的编译和脚本，对新手更友善。
5. 实操过程中的关键坑与排查方法
5.1 展开后能带出现大量伪能带
这是unfold流程中最常见的现象。通常情况下，伪能带的出现源于K点路径选取不当或权重阈值设置太低。BandUP画图时可以设置一个最小权重阈值，把权重小于某个值的点直接丢弃，这样图面会干净很多。我一般取0.05到0.1之间，太小的值会让图片布满杂点，太大则可能丢掉真实物理信息。
还有一类特殊情况：如果掺杂浓度太高（比如2x2x2超胞只有8个原子就替换1个），杂质原子间的相互作用会比较强，展开后的能带会呈现明显的色散，而不是一个平缓的杂质态。这时用unfold也救不了太多，物理本质上已经不再是“稀掺杂”极限了。想看到类似实验中的孤立杂质能级，最好用3x3x3以上的超胞。
5.2 展开前后的能带位置对不上
如果你比较展开后的能带和原胞能带，发现带隙大小差异很大，先不要急着怀疑unfold程序出了bug。最可能的原因是自洽计算和能带计算没有使用同一个结构。这就又回到前面强调的：非自洽计算必须用结构优化后的CONTCAR，而不是最开始构建超胞时的POSCAR。这两者差异大的时候，能带整体偏移可以达到零点几个电子伏特。
另外提醒一点：OUTCAR里输出的费米能级位置也要注意。展开后能带的零点通常取到费米能级或VBM，需要根据体系性质来选择。对半导体，推荐用VBM做参考零点；对金属，直接用费米能级更合理。两种选择对能带图的阅读体验影响很大，文章里一定要保持一致。
5.3 PROCAR与WAVECAR文件缺失或截断
能带展开依赖的是投影信息或波函数数据。如果计算过程中磁盘空间不足导致WAVECAR没有完整写出，后面的展开会直接失败或者产生错误结果。检查磁盘空间是经验之谈，而这里有一个隐藏注意点：VASP在非自洽能带计算时，WAVECAR会很大，尤其超胞原子数多、k点路径又长的情况下，几十GB都是可能的。建议在计算前就规划好存储位置。
另外，VASPKIT的unfold模块在读取数据时，要求能带计算时LORBIT=11这一项已开。这个参数会让PROCAR里保存每个能带在原子和轨道上的投影分量，是展开程序计算权重的数据基础。如果忘了设置，我遇到过程序直接报错或者输出空白数据的情况，排查半天才发现是LORBIT的问题。
5.4 一个与原子坐标序相关的典型报错
BandUP在计算映射时对原子序号非常敏感。超胞中原子排列顺序与原胞的周期性格点不匹配时，程序会提示找不到合适映射。解决办法是先对原胞和超胞做对称性分析，确认超胞确实是由该原胞扩展而来。有时候构建超胞时做了原子重排或者用了非标准设定，表面上看没问题，实际上行不行得通要验证。
我在一个层状材料体系上踩过一次这样的坑。当时超胞是在原胞基础上旋转了一个非标准角度构建的，VESTA里看起来正常，但BandUP一直报错。后来改用ASE的make_supercell方法重建结构，问题直接消失。所以建议你构建超胞时尽量用程序化的标准方法，避免在建模时引入非预期的对称性破坏。
6. 结果分析与可视化：从数据到结论
6.1 用Python快速绘制展开能带
展开之后的数据通常存成每行包含“能量、权重、k点路径距离”的文本文件。用Python画图非常方便：
import numpy as np import matplotlib.pyplot as plt data = np.loadtxt('unfolded_band.dat') kdist = data[:, 0] energy = data[:, 1] weight = data[:, 2] fig, ax = plt.subplots(figsize=(5, 6)) sc = ax.scatter(kdist, energy, c=weight, cmap='hot_r', s=1, vmin=0.05, vmax=1.0) ax.set_ylim([-3, 3]) ax.set_ylabel('Energy (eV)') plt.colorbar(sc, label='Spectral weight') plt.tight_layout() plt.savefig('unfolded_band.png', dpi=300)画图时颜色映射建议用hot_r这类深色背景的方案，谱权重高的点呈现亮色，低权重的伪能带自然藏进背景里，图面层次非常分明。
6.2 怎么从展开能带里读取关键信息
拿到展开能带图之后，要看三个核心信息：
第一，带隙大小和带边位置。展开后的价带顶和导带底位置，直接回答了掺杂是否显著改变带隙这个问题。与未掺杂原胞的带隙对比，差值可以量化掺杂效应。
第二，杂质能级的出现。如果掺杂在带隙中引入了能级，展开图上会看到一条权重很高但色散很小的平带。这条带的能量位置对应杂质能级的深度，平带色散小说明杂质态高度局域化，符合替位掺杂浅能级的物理图像。
第三，能带色散变化。展开后的有效能带如果相比原胞能带出现明显的“模糊”或者带边增宽，说明掺杂原胞间的相互作用破坏了原有的晶体周期性。这种信息在折叠能带里是完全看不出来的。
6.3 投影能带与分波态密度的呼应
习惯上我还会把展开能带与分波态密度(PDoS)画在一起。展开能带中识别出的杂质能带，在PDoS中会对应一个明显的杂质原子轨道峰。两者互相印证，结论会更扎实。VASPKIT输出PDoS的方式很简单，LORBIT=11已经记录了各原子的分波信息，直接用vasprun.xml导出即可。
例如你掺杂的是氮元素替代氧化物中的氧，PDoS里氮的2p轨道在带隙中出现的尖峰位置，应当与展开能带里杂质平带的能量完全对应。要是对不上，那就要回查是不是原子坐标或者原子种类弄错了。
7. 个人实操体会与后续扩展思路
做完整个流程之后，我最深的体会是：掺杂超胞计算，真正的门槛不在VASP命令行参数，而在于对“周期性近似”这个大前提的把握。超胞是人为构造的，周期性也是人为重复的，因此引入的任何缺陷本质上都是周期性缺陷阵列，而不是真正的孤立缺陷。unfold能消解能带折叠这个几何问题，但无法消除掺杂原子之间的残余相互作用。理解这层物理，才不会对计算结果做过度解读。
另外，作为长期跑VASP的实践者，我强烈建议每一步计算都留下完整的作业记录。哪个INCAR配哪一版结构、KPOINTS路径怎么定义的、展开用的哪个工具配置，都记录下来。这个项目隔两周再回来看，你大概率会忘掉一多半细节。别问我怎么知道的。
后续扩展方面，如果掺杂体系涉及磁性，可以进一步做自旋极化的展开能带，区分自旋向上和自旋向下的杂质态分布。这样就能判断杂质能级是不是自旋劈裂，对磁性掺杂和稀磁半导体研究尤其重要。另外还能把展开后的能带和光电性质计算（比如介电函数、吸收谱）联合分析，把电子结构和宏观物性串成一条完整的证据链。
如果你刚接触这套流程，建议先用一个简单体系完整跑通一遍，比如硅的替位掺杂，把流程中每个环节的输入输出都摸清楚，再切换到自己的目标材料上。这样即使后续遇到问题，排查范围也能快速缩小，不至于在工具链和物理分析两头同时抓瞎。
