企业官网建设流程全解析

替位掺杂加超胞再unfold，几乎是做材料计算的人都会碰到的组合。先建一个够大的超胞，在原胞的某个格位上放一个替位杂质，跑完VASP之后，再把能带从超胞布里渊区“拉”回原胞布里渊区，好判断杂质对带隙、载流子行为、能级位置的影响。这套流程本身不复杂，但坑不少，尤其是双替位体系，两个杂质之间的距离、格位选择、磁矩初值都会直接影响计算结果。

我一开始接触这个课题时，以为只要把超胞建好、原子一替换、VASP一跑、能带一画就完事。实际走一遍才发现，真正的重头戏在最后一步unfold：不处理，超胞里的能带折叠得密密麻麻，根本分不清哪些是宿主能带、哪些是杂质态；处理不好，图像又会出现一堆权重极低的伪带，白算一场。这篇文章把我的完整流程、参数选择思路、常见坑点整理出来，希望能让后来人少走点弯路。

1. 替位掺杂计算的整体思路：为什么绕不开超胞和unfold

1.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号原子替换为Si：

from ase.build import bulk atoms = bulk('GaAs', 'zincblende', a=5.65) supercell = atoms * (2, 2, 2) # 替换0号原子：Ga位被Si替位 supercell[0].symbol = 'Si' supercell.write('POSCAR', direct=True)

用Pymatgen的方式也差不多：

from pymatgen.core import Structure struct = Structure.from_file('POSCAR') sup = struct * (2, 2, 2) sup.replace(0, "Si") sup.to(filename='POSCAR_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是否需要手动调整，也影响自洽收敛；
• 对称性是否被破坏，是否降为P1，K点必须用Gamma-centered网格；
• 是否需要额外添加真空层（二维体系、表面体系），如果需要，确保真空层足够厚且后续unfold按二维布里渊区处理。

我见过不少人在第三步栽跟头：用Pymatgen的replace后，结构里某两个原子距离短到1.2 Å，VASP跑起来以后不是氧原子飞了就是自洽震荡。建模阶段多看一眼，后面轻松很多。

3. VASP计算参数设置与计算流程

3.1 三步走：优化、静态、能带的衔接逻辑

替位掺杂的计算流程我基本固定为三步：结构优化、静态自洽、能带非自洽。结构优化解决“杂质原子进去后周围原子怎么松弛”的问题；静态自洽得到准确的电荷密度（CHGCAR），并作为能带计算的基础；能带非自洽不重新生成电荷密度，而是沿着给定K路径求解本征值，速度快且结果与静态计算自洽。

优化这一步，我会用ISIF=3配合IBRION=2，让晶胞体积形状和原子位置一起松弛。掺杂体系里杂质会引起局部应力，只放开原子位置不放开晶胞，容易得到一种“人为的弹性应力状态”，不真实。不过，如果宿主晶格非常大，或者确定要用实验晶格常数，ISIF=2只优化原子位置也可以，但要写清楚说明理由。

静态自洽要把K点网格加密，比优化时至少每个方向多一倍。自洽结束后检查OUTCAR里的总能量变化和磁矩是否稳定，再进入能带计算。能带计算注意设置ICHARG=11，从CHGCAR读取电荷密度，同时设NSW=0、IBRION=-1，避免任何多余的结构修改。

3.2 INCAR关键参数：自旋、+U与smearing如何处理

掺杂体系最容易被参数坑到的三个点：自旋极化、强关联修正、smearing方式。我逐个说。

含过渡金属或稀土元素的掺杂，几乎都要开ISPIN=2。比如CeO2里掺Sn、铁基材料掺Co，杂质附近电子局域性强，不限制自旋方向会得到错误的基态。初始MAGMOM怎么给？对于未知体系，我通常给过渡金属离子一个合理的局域磁矩初值，比如3-5 μB，然后看收敛后的总磁矩是否落在合理范围。如果计算稳定后某个过渡金属磁矩变成接近0，可能是初始设置导致陷入非磁性亚稳态，需要换初值重新跑。

DFT+U的取舍：宿主体系本身是强关联、带隙靠+U才能算对的时候，掺杂计算必须带着U跑。U值不能自己乱编，参考同族材料已发表的计算数据。这里强调一点，U值不同，杂质能级在带隙里的相对位置可能有明显差异，做对比研究时一定要保持相同的U值设定，否则趋势比较没有意义。

smearing的选择：半导体或绝缘体用ISMEAR=0（Gaussian smearing）就够，如果只关心总能和带隙，ISMEAR=-5（tetrahedron方法）更准，但-5不适用于自洽后的能带计算。金属性体系用ISMEAR=1配合合适的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（如果有）拷贝进去，然后运行VASPKIT：

mkdir 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, c=weight, s=1, cmap='viridis', vmin=0, vmax=1) plt.ylim(-5, 5) plt.xlabel('K path') plt.ylabel('E - E_F (eV)') plt.colorbar(label='Weight') plt.tight_layout() plt.savefig('unfold_band.png', dpi=300)

判读结果时，先看未掺杂原胞能带的特征能不能在权重图里复现。比如原本有个带隙，权重图上在原来带隙位置依然有空白区，说明掺杂没有破坏基本半导体特征；带隙内出现了高权重色散或低权重平带，分别对应类宿主杂质带和局域杂质态。如果发现“带隙里出现一条完全平整、权重又很低的带”，大概率是数值问题，参考第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跑出图来确认图像合理，再上自己真正关心的复杂体系。这个预演成本很低，但能帮你把工具链里的每一个环节都验证到位，真到了几万元素的大超胞计算，心态会稳很多。