替位掺杂加超胞再unfold几乎是做材料计算的人都会碰到的组合。先建一个够大的超胞在原胞的某个格位上放一个替位杂质跑完VASP之后再把能带从超胞布里渊区“拉”回原胞布里渊区好判断杂质对带隙、载流子行为、能级位置的影响。这套流程本身不复杂但坑不少尤其是双替位体系两个杂质之间的距离、格位选择、磁矩初值都会直接影响计算结果。我一开始接触这个课题时以为只要把超胞建好、原子一替换、VASP一跑、能带一画就完事。实际走一遍才发现真正的重头戏在最后一步unfold不处理超胞里的能带折叠得密密麻麻根本分不清哪些是宿主能带、哪些是杂质态处理不好图像又会出现一堆权重极低的伪带白算一场。这篇文章把我的完整流程、参数选择思路、常见坑点整理出来希望能让后来人少走点弯路。1. 替位掺杂计算的整体思路为什么绕不开超胞和unfold1.1 替位掺杂为什么必须用超胞替位掺杂的本质是把晶体原有格位上的一个或者几个原子换成杂质原子。这里的“换”不是像相变那样重构整个原子排布而是在周期性晶格上引入一个局部扰动。由于VASP这类第一性原理软件建立在周期性边界条件上模拟对象必须是有限大小的晶胞在三维空间无限重复想在一个所有格位都被占据的完美原胞里“单独”放一个杂质原子物理上不存在每一个原胞里都会被塞进一个杂质相当于100%的高浓度掺杂根本不是稀掺杂场景。所以超胞是绕不开的。把原胞在某个方向重复几次得到一个包含大量原子的周期盒子再在这个盒子的某个格位上替换原子杂质仍然在三维周期中出现但两个最近邻杂质镜像之间隔着一段晶格距离这样才有可能近似真实的稀掺杂环境。常见的做法是2×2×2、3×3×3或者更大具体看宿主晶格和想要的掺杂浓度。比如金刚石结构一个原胞8个原子2×2×2超胞64个原子替换1个掺杂浓度为1/64约1.56%替换2个浓度约3.125%。这个数字在实验里已经算很高的掺杂量了如果要模拟千分之一级别的掺杂超胞要建到非常大机时往往不允许。实际研究里更关注的是“杂质浓度对结果的影响趋势”而不是严格复现某个实验浓度。1.2 双替位体系的特殊考量标题里特意写了“双替位”说明涉及两个杂质原子。双替位分为两种常见情况两个相同元素占据两个等价格位或者两个不同元素分别占据不同子晶格格位。前者常见于同族元素掺杂比如把两个最近的阳离子同时替换成另一种金属后者常见于联合掺杂比如同时替换一个阳离子和一个阴离子来平衡电荷。双替位比单替位麻烦的不是VASP计算本身而是建模阶段的“距离控制”。两个杂质如果放得太近会形成强的相互作用杂质能级展宽成杂质带和单杂质的结果对不上放得太远又可能超出超胞容纳能力。我个人的经验是杂质镜像间距和杂质-杂质间距都最好大于10 Å严谨一点做到12-15 Å。这个距离下局域杂质态之间的波函数交叠已经很小结果基本可以代表稀掺杂行为。另外双替位还牵涉自旋排列问题。如果两个杂质都有未配对电子需要设置合理的初始磁矩并且在静态计算后检查OUTCAR里的总磁矩是否收敛到预期的物理状态。铁磁耦合和反铁磁耦合可能给出完全不同的能带结果这一点比单替位更容易被忽视。1.3 unfold能带反折叠解决什么问题超胞建好后直接算能带输出的能带图会让人一头雾水。原因在于超胞的倒空间比原胞小原胞布里渊区的高对称点会折叠到超胞布里渊区的不同位置原本一条清晰的色散关系被“折叠”成许多条能带。杂质含量越低、超胞越大折叠越严重。从折叠后的能带里很难回答一个最核心的问题掺杂之后原来的带边位置怎么移动杂质能级在带隙里什么位置有效质量有没有明显变化unfold做的事情就是把超胞的波函数投影回原胞的基组上给每一条超胞能带赋一个“权重”。权重高说明这条带的电子态本质上是原胞能带的延伸权重低说明它是折叠伪态或杂质局域态。于是我们可以把权重画成颜色或者气泡大小还原出掺杂体系的“有效原胞能带”再和未掺杂的原胞能带对比就能直观看出杂质的影响。这套方法的物理基础并不复杂超胞波函数可以表示为原胞波函数在超胞倒格矢下的折叠展开unfold就是把这个展开过程中的系数提取出来。实际做的时候不需要自己写公式VASPKIT和BandUp都实现了完整流程但理解了这个思想才不会在结果异常时一头雾水。2. 建模实操替位掺杂超胞的构建与检查2.1 从基础结构到超胞工具与代码示例我习惯用ASE和Pymatgen做建模两者都可以直接从Materials Project或实验CIF文件读结构也可以手动构建。不管用哪个都要注意第一步是拿到“原胞”或至少是合理的单胞因为后续unfold需要把超胞和原胞的晶格关系对应起来。以ASE为例从晶格常数构建一个闪锌矿结构的GaAs原胞然后扩成2×2×2超胞再把0号原子替换为Sifrom ase.build import bulk atoms bulk(GaAs, zincblende, a5.65) supercell atoms * (2, 2, 2) # 替换0号原子Ga位被Si替位 supercell[0].symbol Si supercell.write(POSCAR, directTrue)用Pymatgen的方式也差不多from pymatgen.core import Structure struct Structure.from_file(POSCAR) sup struct * (2, 2, 2) sup.replace(0, Si) sup.to(filenamePOSCAR_supercell)这里有个容易被忽略的点替换原子时要确认0号位点确实是你想替换的子晶格。如果不确定可以先用supercell.symbols打印元素顺序或者通过Pymatgen的indices筛选特定元素的所有格位。2.2 替位格位选择与掺杂浓度设计格位选择不能说“看着像就换”。得先明确宿主材料中哪个元素在哪个Wyckoff位置杂质元素进入后更倾向于占据哪种配位环境。比如掺杂钙钛矿A位和B位配位数完全不同替位难度和电子结构影响也不同。一个很实用的办法是算一下掺杂形成能比较不同格位方案的能量高低能量最低的往往是优先占据的格位。掺杂浓度也不是随便定的。我通常先做一组不同超胞尺寸的计算比如2×2×2、3×3×3、4×4×4把杂质能级位置或带隙随浓度变化的趋势画出来找到基本稳定的超胞尺寸。这样虽然多花一点机时但能避免结果被“伪浓度效应”污染。这里也顺带说一下周期镜像很多人只盯着超胞里的两个杂质忘了它们通过周期性边界条件会形成无限排列。如果超胞不够大哪怕两个杂质放在超胞两端也会通过镜像产生短距离相互作用。所以我在确定超胞尺寸时会先在ASE里直接测量最近镜像距离确认大于10 Å之后才进入下一步。2.3 双替位体系的间距控制与位点组合双替位的建模细节值得单独说。第一步明确两个杂质分别占据什么格位第二步在满足周期性边界条件的前提下尽量让两个杂质在超胞内的距离最大化第三步用周期性条件检查杂质与自己在相邻镜像中的距离以及两个杂质之间的真实最短距离。实际操作里我常常在超胞里列出所有候选格位计算任意两个格位之间的周期最小距离然后挑出距离最远的格位组合。这一步用ASE的get_distance加上周期图像判断或者直接手动构造多个候选结构比较。别偷懒随便选两个格位可能会导致两个杂质在周期上几乎挨在一起算出来的“双替位”其实等价于一个杂质二聚体。如果做的是异种元素双替位还要额外注意电荷补偿。比如在氧化物里用一个低价阳离子替换高价阳离子会造成净空穴掺杂同时用一个高价阴离子替换低价阴离子两者可能形成电荷补偿能带里浅能级的行为会很不一样。建模阶段就要想清楚电荷平衡和载流子类型否则跑完才发现电子数不对。2.4 建模后的结构检查清单结构文件写出来之后我不会立刻丢进VASP。先花两分钟做一轮静态检查能省后面可能几天的debug时间元素种类和数量是否符合预期特别是替换原子后有没有误伤其他原子掺杂位点周围是否有过近的原子间距小于共价半径之和往往意味着结构不合理超胞内总电子数是多少这关系到NELECT是否需要手动调整也影响自洽收敛对称性是否被破坏是否降为P1K点必须用Gamma-centered网格是否需要额外添加真空层二维体系、表面体系如果需要确保真空层足够厚且后续unfold按二维布里渊区处理。我见过不少人在第三步栽跟头用Pymatgen的replace后结构里某两个原子距离短到1.2 ÅVASP跑起来以后不是氧原子飞了就是自洽震荡。建模阶段多看一眼后面轻松很多。3. VASP计算参数设置与计算流程3.1 三步走优化、静态、能带的衔接逻辑替位掺杂的计算流程我基本固定为三步结构优化、静态自洽、能带非自洽。结构优化解决“杂质原子进去后周围原子怎么松弛”的问题静态自洽得到准确的电荷密度CHGCAR并作为能带计算的基础能带非自洽不重新生成电荷密度而是沿着给定K路径求解本征值速度快且结果与静态计算自洽。优化这一步我会用ISIF3配合IBRION2让晶胞体积形状和原子位置一起松弛。掺杂体系里杂质会引起局部应力只放开原子位置不放开晶胞容易得到一种“人为的弹性应力状态”不真实。不过如果宿主晶格非常大或者确定要用实验晶格常数ISIF2只优化原子位置也可以但要写清楚说明理由。静态自洽要把K点网格加密比优化时至少每个方向多一倍。自洽结束后检查OUTCAR里的总能量变化和磁矩是否稳定再进入能带计算。能带计算注意设置ICHARG11从CHGCAR读取电荷密度同时设NSW0、IBRION-1避免任何多余的结构修改。3.2 INCAR关键参数自旋、U与smearing如何处理掺杂体系最容易被参数坑到的三个点自旋极化、强关联修正、smearing方式。我逐个说。含过渡金属或稀土元素的掺杂几乎都要开ISPIN2。比如CeO2里掺Sn、铁基材料掺Co杂质附近电子局域性强不限制自旋方向会得到错误的基态。初始MAGMOM怎么给对于未知体系我通常给过渡金属离子一个合理的局域磁矩初值比如3-5 μB然后看收敛后的总磁矩是否落在合理范围。如果计算稳定后某个过渡金属磁矩变成接近0可能是初始设置导致陷入非磁性亚稳态需要换初值重新跑。DFTU的取舍宿主体系本身是强关联、带隙靠U才能算对的时候掺杂计算必须带着U跑。U值不能自己乱编参考同族材料已发表的计算数据。这里强调一点U值不同杂质能级在带隙里的相对位置可能有明显差异做对比研究时一定要保持相同的U值设定否则趋势比较没有意义。smearing的选择半导体或绝缘体用ISMEAR0Gaussian smearing就够如果只关心总能和带隙ISMEAR-5tetrahedron方法更准但-5不适用于自洽后的能带计算。金属性体系用ISMEAR1配合合适的SIGMA。掺杂体系最容易出现的情况是本身是半导体掺入杂质后出现部分占据的杂质带自洽过程中出现类似金属的态密度尾部导致smearing振荡。这时候把SIGMA稍微调大到0.1-0.2 eV同时提高NELM能明显提升收敛稳定性。一个实用的INCAR模板大致如下细节根据体系微调Sys doped_supercell ISTART 0 ICHARG 2 ENCUT 520 EDIFF 1E-5 EDIFFG -0.02 IBRION 2 ISIF 3 NSW 200 NELM 150 ISMEAR 1 SIGMA 0.10 ISPIN 2 MAGMOM 8*3 2*1 6*0 LORBIT 11 NCORE 16最后一行NCORE根据并行核数和节点配置设一般取每个节点的物理核数附近。含过渡金属时ENCUT建议取POTCAR里最大ENMAX的1.3倍以上别省那点机时。3.3 KPOINTS与K路径的对应关系超胞的K点网格要按超胞倒空间来生成。优化和静态计算用Gamma-centered网格密度根据超胞大小调整一般2×2×2超胞先试4×4×4收敛测试后加密。到了能带计算阶段关键就变了这条K路径不是随便选的而是由unfold需求驱动的。你的目标是最后在原胞布里渊区里画出分辨率的能带所以非自洽计算的K点必须沿着原胞高对称K路径映射到超胞倒空间后的路径来取。映射关系由超胞矩阵和原胞矩阵的变换决定工具会自动处理但你要确保输入的原胞POSCAR和扩胞时用的那个原胞完全一致包括基矢方向否则映射关系错乱后面unfold结果全是错的。我用的办法是先用VASPKIT的K点路径生成功能给出原胞的POSCAR生成标准高对称路径再让VASPKIT根据超胞关系把这条路径换算成超胞的K坐标写入能带计算的KPOINTS文件。3.4 双替位体系的计算注意点双替位体系里两个杂质位置不同容易形成对称性很低的构型自洽计算比单替位更容易不收敛。建议优化阶段先用粗K点和中等ENCUT快速跑到一个粗糙的收敛态再用精参数从粗糙结果继续跑。VASP的ISTART和ICHARG可以配合好从WAVECAR或CHGCAR续算不需要每次都从头开始。另外如果两个杂质是相同元素还要考虑它们在超胞里是等价位点还是不等价位点。等价位点会让计算保持对称性能带里杂质态色散明显一些不等价位点则可能产生两套杂质能级。当你看到能带里出现两个靠得很近的杂质带先别急着下物理结论检查一下是不是两个位点的能量差异导致的而不是掺杂浓度效应。4. unfold实操全流程从超胞能带到原胞能带4.1 工具选择与输入数据准备unfold的工具有好几个我用得最多的是VASPKIT偶尔也用BandUp。VASPKIT的优势是集成在常见工作流里交互式菜单方便BandUp是独立工具逻辑更透明支持自旋极化体系输出权重数据。PyProcar也有unfold能力适合喜欢Python画图的人。原则就一条哪个工具能让你最快看到正确图像就用哪个。无论用哪个工具都需要准备好以下输入超胞能带计算的PROCAR。这个文件包含每条能带在每个K点的波函数展开系数是unfold的数据来源非自洽计算的KPOINTS文件里面的K点要符合映射关系超胞POSCAR和原始胞POSCAR。原始胞可以是未被掺杂的原胞但必须保证它就是扩胞前那个胞如果计算了自旋极化确认PROCAR里两个自旋通道的分量都完整。准备数据时最常遇到的坑是K点路径不是从原胞路径映射来的而是自己在超胞布里渊区里随意选的。这样unfold工具根本没有办法把超胞倒空间坐标映射回原胞坐标结果自然是一堆乱码。先把路径映射关系跑通再谈结果。4.2 用VASPKIT完成unfold并得到权重能带以VASPKIT为例实际过程大致如下。我建议在能带计算输出目录里建一个子目录把POSCAR、KPOINTS、PROCAR、OUTCAR、WAVECAR如果有拷贝进去然后运行VASPKITmkdir unfold cd unfold cp ../POSCAR ../KPOINTS ../PROCAR ../OUTCAR . vaspkit交互界面里进入能带反折叠功能菜单具体编号以你使用的VASPKIT版本为准一般与“Band structure unfolding”相关。按提示输入原始胞POSCAR路径、超胞POSCAR路径、能量范围、是否需要自旋分解等。VASPKIT会读取PROCAR完成投影输出包含权重信息的能带数据文件。如果你不小心把原始胞POSCAR给成了超胞自身或者给成了conventional cell而不是primitive cell映射关系就会错。最直观的表现是unfold之后的权重图上本来应该高权重的能带也变成了散碎的亮斑。这时候第一件事不是怀疑物理而是检查两个POSCAR是否存在准确的倍数关系。4.3 能带图绘制与结果判读unfold输出的是一个K点序列每个K点有若干条能量本征值和对应的权重。画图推荐散点图横坐标是原胞K路径的累积距离纵坐标是能量点的大小或颜色映射权重。权重低于某个阈值比如0.1的点可以淡化它们主要是折叠伪态不值得关注。以Python的matplotlib为例画权重能带的一小段示意逻辑import numpy as np import matplotlib.pyplot as plt data np.loadtxt(unfold_band.dat) k_dist data[:, 0] energy data[:, 1] weight data[:, 2] plt.figure(figsize(6, 8)) plt.scatter(k_dist, energy, cweight, s1, cmapviridis, vmin0, vmax1) plt.ylim(-5, 5) plt.xlabel(K path) plt.ylabel(E - E_F (eV)) plt.colorbar(labelWeight) plt.tight_layout() plt.savefig(unfold_band.png, dpi300)判读结果时先看未掺杂原胞能带的特征能不能在权重图里复现。比如原本有个带隙权重图上在原来带隙位置依然有空白区说明掺杂没有破坏基本半导体特征带隙内出现了高权重色散或低权重平带分别对应类宿主杂质带和局域杂质态。如果发现“带隙里出现一条完全平整、权重又很低的带”大概率是数值问题参考第5节排查。另外建议把未掺杂原胞的能带也画在同一张图里作为参照能非常直观地看到掺杂引起的带边相对移动。很多论文里那些“掺杂后带隙变小0.3 eV”的结论就是这么一层一层对比出来的。5. 常见问题与排查技巧实录5.1 自洽不收敛先查电荷密度振荡掺杂体系自洽不收敛几乎都发生在带隙内出现半占据杂质态的情况。电荷密度在杂质位点附近来回振荡电子数在两个相近的电子态之间反复跳变。我排查的顺序是先看OUTCAR里最近几步的能量差确认是单调变化还是振荡。如果是振荡把AMIX从默认的0.2降到0.1左右BMIX也同步降低或者改用线性混合。再把SIGMA从0.05往上调到0.1甚至0.15软化费米面附近的占据函数。还有一招是把NELM调到150以上给足迭代次数。如果自洽就是稳定不下来而且你加了自旋很可能是初始MAGMOM给错导致自旋通道之间来回切换。试着用不同的磁矩初值跑几次看最终能否收敛到同一个态。别小看这一步我遇到过因为MAGMOM里把两个杂质的自旋方向设成同一方向而实际基态是反铁磁排列结果算了很久都收敛不到物理正确的解。5.2 unfold结果出现大量伪带unfold图一出来满屏都是权重不到0.05的亮点这种“脏图”我见过太多次。最根本的原因通常是K点路径没有正确映射到超胞倒空间。检查KPOINTS文件里K点是否写成了超胞倒空间的真实坐标如果K点坐标看起来像是原胞的高对称点直接沿用那就是映射关系错了。另一个常见原因是输入的原胞POSCAR不是真正的prime cell而是conventional cell。比如面心立方晶体的conventional cell包含4个原子而primitive cell只包含1个原子。你用conventional cell扩成2×2×2超胞然后unfold时又把conventional cell当作原胞输入映射矩阵的根号关系就会不对。我踩过这个坑之后每次都会先确认POSCAR的原子数是否和该结构的最小原胞一致。5.3 杂质能级或带边位置不合理的处理unfold后杂质能级的位置和预期差很远比如应该出现浅施主能级结果算出来是深能级。这类问题大概率不是unfold工具造成的而是前期的掺杂浓度设计太密杂质波函数之间相互作用导致杂质能级展宽或移动。处理方法增大超胞、降低掺杂浓度重复计算。如果增大的超胞太大机时吃不消可以先用粗K点和低ENCUT跑一个趋势性结果确认杂质能级位置开始收敛再用精参数跑最终构型。另外半导体体系的DFT带隙本身偏低杂质能级相对带边的位置会受带隙误差影响必要时用杂化泛函或GW做交叉验证但那样机时开销非常大一般放在趋势性研究之后再做。5.4 资源消耗与并行策略问题超胞越大K点越多自旋又翻倍计算量很容易到几万核时。我常用的策略是分阶段控制成本优化用低密度K点和中等ENCUT静态计算时K点加密、ENCUT提到最终值能带计算用非自洽跑它本身不生成电荷密度K点数量虽然多但每个K点只需迭代几次整体不会太夸张。并行参数上NPAR或NCORE不要设得太小否则频繁的多节点通信会拖垮计算效率。我通常把NCORE设成单节点物理核数节点内共享内存效率最高。内存不足时可以通过降低KP点的并行分组来缓解但一般情况下掺杂体系的内存压力主要来自大基组和大量K点适当控制ENCUT比盲目加节点更有效。跑完后第一时间检查OUTCAR里的计算时间和并行效率如果加速比很低及时调整参数再继续而不是干等。最后再分享一点实际经验这套流程我前前后后跑了几十次最大的体会是unfold最耗时间的反而不是计算而是前期的K路径映射和模型准备。超胞建好、原子替换好不代表万事大吉KPOINTS不对后面算得再准也没法画图。建议第一次接触这个流程的朋友先拿一个原子数很少、能带关系简单的体系整体走通一遍比如8原子原胞扩成64原子的超胞替换一个原子把unfold跑出图来确认图像合理再上自己真正关心的复杂体系。这个预演成本很低但能帮你把工具链里的每一个环节都验证到位真到了几万元素的大超胞计算心态会稳很多。