替位掺杂超胞能带展开实操：从折叠到有效能带

1. 为什么要做替位掺杂超胞的unfold1.1 超胞能带折叠是超胞的“副作用”我最早接触掺杂计算时也有过很迷惑的阶段单胞能带明明很干净为什么一把体系扩成 2x2x2 或者 3x3x3 的超胞能带图就完全变了个样后来才意识到问题不是算错了而是能带折叠band folding在作怪。单胞变成超胞实空间体积变大倒空间布里渊区就对应缩小。原来在单胞里两个不同 k 点上的能带落到超胞布里渊区后可能被压到同一个 k 点上原本一条完整的色散关系被“折”成了好几段。超胞越大折叠越厉害。函数关系上一个超胞计算出的本征态实际上是多条原始能带的叠加只是我们拿到手的是一堆“混在一起”的本征值看不出哪个属于原来的哪条带。这时候如果直接拿超胞能带去和实验数据比或者去讨论带隙、有效质量、杂质能级位置基本是白搭。必须做能带展开unfold把超胞能带按单胞的晶格周期重新投影回去恢复出“有效能带结构”。1.2 替位掺杂到底改了谁的能带替位掺杂就是把超胞里某个原子换成另一种元素。比如硅超胞里换上一个磷原子或者氧化物里换一个过渡金属。这种做法的核心意图是在一个尽可能接近真实晶体的周期体系中看杂质原子带来的局域能级、电荷转移、晶格弛豫以及对原有能带边的影响。但问题在于掺杂超胞本身已经破坏了原胞的平移对称性。严格来说超胞中不再存在“原胞能带”只有超胞能带。可我们又想回答一个物理问题这个杂质到底在带隙里引入了什么态它跟宿主晶格的主带之间是什么关系要回答这个问题就需要把超胞能带的波函数投影到单胞的布洛赫态上看它像原来的哪条带。这正是 unfold 算法做的事情。所以替位掺杂体系和 unfold 几乎是绑定的。单胞能带只是参考超胞能带只是中间量unfold 之后的带才是能讨论物理的那条“有效能带”。1.3 什么场景必须做unfold如果你的掺杂浓度很低比如 100 多个原子的超胞只换一个原子直接看超胞能带会出现海量折叠带杂质态细节被完全淹没这时候必须 unfold。下面这些场景我实际都遇到过都属于“不做 unfold 就无法正确解读”的情况掺杂后隙态in-gap state的判据想知道某个杂质能级是深能级还是浅能级在超胞能带里几乎看不清楚。载流子有效质量超胞能带的色散被折叠打断局部斜率失真展开后才能得到接近单胞的有效质量。能带边对齐比较掺杂前后带边位置变化同样需要展开到同一套晶格坐标系下。合金或随机占位体系虽然标题是替位掺杂但原理也适用于无序合金的虚拟晶胞。一句话总结只要你想把“超胞的结果”映射回“单胞的语言”unfold 就是必经之路。2. 环境准备与核心工具选型2.1 软件栈VASP VASPKIT BandUP这套流程的底层计算还是 VASP稳。我们用到的辅助工具主要有VASPKIT做结构文件转换、KPOINTS 生成、基础后处理都非常方便。BandUP做 unfold 的老牌开源工具物理上很靠谱适合进阶用户。VESTA可视化掺杂结构确认原子位置没有问题。ASE 或 phonopy生成超胞很方便尤其是需要把单胞扩成特定尺寸超胞时。我个人的分工是建模和常规数据准备用 VASPKITunfold 后处理优先用 BandUP因为它的谱权重计算更透明。如果你只想快速出图VASPKIT 里也集成了 unfold 模块虽然运算细节不如 BandUP 可调但胜在省事。需要提醒一句VASPKIT 和 BandUP 的版本差异会影响菜单位置和命令行参数下面提到的交互菜单数字只是我机器上的情况换版本时先看一遍自带的帮助文档。2.2 单胞、超胞和坐标系的统一做 unfold 的前提是有一套干净的“参考单胞”。这里说的单胞不是随意选的结构基元而是你要把超胞能带展开回去的那个最小晶胞。换句话说参考单胞决定了展开后能带路径的布里渊区。对于替位掺杂参考单胞通常取宿主晶体的原胞或者惯用单胞。另一个关键是超胞生成方式要记录清楚如果你用phonopy -d --dim2 2 2 -c POSCAR生成超胞phonopy 会在输出里保留超胞矩阵信息。如果用 ASE 的make_supercell建议顺手打印一下cell_to_supercell_matrix后面 unfold 时要用。如果用 VESTA 手动扩胞也可以但一定要保证原始晶胞和超胞的坐标轴定义一致。很多 unfold 报错都出在坐标系不一致上比如参考单胞是惯用晶胞超胞却是在原胞基础上扩的或者两种文件的原子排序不同导致投影时波函数相位对不上。所以在开始计算之前就把单胞、超胞两类 POSCAR 放到同一个文件夹里保证原子类型顺序和坐标系完全对应能省掉后面一半的排错时间。2.3 输入文件的组织方式我的目录习惯是按步骤分目录避免后面搞混project/ ├── 1_primitive_cell/ │ └── POSCAR_prim ├── 2_supercell_build/ │ └── POSCAR_sup ├── 3_optimization/ │ ├── POSCAR │ ├── POTCAR │ ├── INCAR │ └── KPOINTS ├── 4_scf/ ├── 5_band/ │ ├── KPOINTS_band │ └── ... └── 6_unfold/这样每次计算都是独立目录算坏了不用从头来。尤其是静态自洽和能带计算最好不要在同一个目录里反复覆盖 WAVECAR 和 CHGCAR。算完之后把这个文件夹留好写论文补数据时能少走很多弯。3. 从建模到结构优化单替位掺杂超胞的常规流程3.1 掺杂位点选择和POSCAR构建我以硅的超胞为例假设单胞是金刚石结构的 2 原子原胞做一个 2x2x2 超胞共 64 个原子。单替位掺杂就是把这 64 个原子中的某一个硅原子替换成杂质原子比如磷。具体操作步骤是用 VASPKIT 或者 ASE 生成 64 原子的超胞得到POSCAR_sup。用 VESTA 打开超胞选中一个想要替换的位点。这一步可以看成“把某个原子换成另一种元素”但不要直接在 VESTA 里改容易在导出时破坏原子索引顺序。回到文本编辑器手动改 POSCAR。把行末的元素标签从Si改成Si P对应的原子数量从64改成63 1然后在坐标列表里把目标行的原子类型对应上。比如原始 POSCAR 可能是Si 1.0 5.4307 0.0000 0.0000 0.0000 5.4307 0.0000 0.0000 0.0000 5.4307 Si 64 Direct ...替换后变成Si P 1.0 5.4307 0.0000 0.0000 0.0000 5.4307 0.0000 0.0000 0.0000 5.4307 Si P 63 1 Direct ...注意VASP 的 POSCAR 中每种元素的坐标顺序必须和元素名称一一对应。也就是说前面 63 个坐标是 Si最后一个坐标必须正是你要替换成 P 的那个原子。改完以后用 VESTA 重新打开一次确认掺杂原子周围的键长、配位环境是否合理。选择掺杂位点时还要关注周期性镜像问题。超胞里原子是周期重复的如果两个掺杂位点在相邻镜像中距离太近实际上会变成“双掺杂”效应单替位掺杂的化学配置就失效了。一句话经验选位点时尽量别挑超胞边界上的原子。3.2 INCAR/POTCAR/KPOINTS参数建议POTCAR 的拼接顺序必须与 POSCAR 元素顺序一致。比如上面是Si P那就先 cat 硅的 POTCAR再 cat 磷的 POTCARcat POTCAR.Si POTCAR.P POTCAR然后检查一下grep TITEL POTCAR确保顺序正确。很多人第一次做掺杂就挂在这一步VASP 直接报错或者算出来的电荷密度混乱。INCAR 我给的典型参数是SYSTEM doped Si supercell ISTART 0 ICHARG 2 PREC Accurate ENCUT 520 EDIFF 1E-6 EDIFFG -0.02 IBRION 2 ISIF 3 NSW 100 ISMEAR 0 SIGMA 0.05 LREAL Auto LORBIT 11 LWAVE .TRUE. LCHARG .TRUE. NPAR 4对半导体体系ISMEAR 0 比 ISMEAR -5 更适合结构优化因为可以更平滑地处理带隙附近的占据态。ISIF 3 表示同时优化晶胞参数和原子位置对替位掺杂引起的晶格膨胀或收缩比较重要。如果体系原子数多NSW 可以放宽到 200但一般 100 步内能收敛。KPOINTS 我建议采用 VASP 6 支持的自动 k 间距格式更省心K-Spacing 0.03 Gamma这个 0.03 单位是 1/Å数值越小 k 点越密。对 2x2x2 的 64 原子超胞实际生成的 k 网格大概在 4x4x4 到 6x6x6 之间已经足够优化了。如果你的 VASP 版本旧直接写显式 Monkhorst-Pack 也行Monkhorst-Pack 5 5 5 0 0 0结构优化计算量大我用 32 核跑 64 原子体系大概几个小时就能收敛这取决于体系大小和赝势精度。3.3 结构优化完成后的检查清单优化完不建议直接进入能带计算先检查下面几项别急着往下走查收敛标志grep reached required accuracy OUTCAR没有这个标志说明离子步没跑完。查受力grep RMS OUTCAR一般优化结束后力的 RMS 在 0.01 eV/Å 以下。查压力grep external pressure OUTCAR压力过大说明晶胞参数没有完全释放ISIF 设置有问题。查最终构型cp CONTCAR POSCAR作为下一步的输入千万别用优化前的 POSCAR。另外把优化前后掺杂原子附近几个键长对比记录一下。比如 P 的共价半径和 Si 不同替位后周围 Si 原子会有微小的结构弛豫这个弛豫量直接影响杂质能级的位置论文里通常会提到。4. 能带计算与原胞K路径设计4.1 静态自洽和能带计算的差异结构优化完成后不要直接拿优化后的结构算能带。我自己的习惯是先做一次静态自洽计算把电荷密度CHGCAR和波函数WAVECAR稳定下来再用这个电荷密度做能带计算。静态计算的 INCAR 和优化时比改动不大ISTART 1 ICHARG 2 NSW 0 IBRION -1 ISIF 2 ISMEAR 0 SIGMA 0.05 LORBIT 11 LWAVE .TRUE. LCHARG .TRUE.这里的核心逻辑是能带计算要在自洽电荷密度基础上一步对角化得到本征值不要边迭代边出能带。如果直接从优化那一步继续算能带SCF 收敛质量参差不齐出来的能带会有噪声unfold 后的权重也会很乱。4.2 KPOINTS高对称路径的换算这是整个流程里最容易踩坑的部分。对于超胞能带计算KPOINTS 不能直接沿用单胞的路径。原因很简单超胞和单胞的倒空间基矢不同同样的高对称点在两套坐标里的分数坐标不一样。如果只是按单胞的习惯写G 0.0 0.0 0.0、X 0.5 0.0 0.0这类分数坐标给超胞VASP 实际上会把它当成超胞倒空间里的不同点算出来的路径根本不是单胞里你想要的那条直线。正确的做法是先把单胞的高对称路径转换成超胞倒空间中的坐标。操作上我会借助 unfold 工具的前处理模块自动生成KPOINTS也可以手动换算。比如单胞的倒格子基矢b_prim和超胞倒格子基矢b_sup之间的关系是b_sup b_prim · M^{-1}严格来说还要看晶胞矩阵正倒空间的定义但关键是超胞分数坐标 k_sup 与原始分数坐标 k_prim 之间有一个线性变换如果你不想手动算矩阵最简单的方案是把原始单胞能带路径的每个 k 点坐标写成笛卡尔坐标直接塞给超胞计算。因为能带路径在绝对倒空间里是连续的直线这种写法不依赖具体分数坐标系。格式大致是Line-mode Cart 0.0000 0.0000 0.0000 G 0.5000 0.0000 0.0000 X不过更省事的是用 BandUP 的前处理工具它能根据你提供的单胞和超胞直接生成一个可用的 band KPOINTS 文件。4.3 能带计算的INCAR调整能带计算本身不需要再优化结构INCAR 里打开波长和波函数相关的选项为 unfold 做准备SYSTEM band calculation ISTART 1 ICHARG 11 NSW 0 IBRION -1 ISIF 2 ISMEAR 0 SIGMA 0.05 LORBIT 11 LREAL Auto LWAVE .TRUE. LCHARG .FALSE. NELMIN 5LORBIT 11是输出投影权重用的后面即便不做 unfold也可以用于分析杂质态的轨道成分。LWAVE .TRUE.是因为部分 unfold 工具需要读取波函数信息。LCHARG .FALSE.不输出电荷密度省点磁盘空间。这里有个细节ICHARG 11表示读取 CHGCAR 后不再自洽直接算本征值和能带。但要注意 VASP 仍然需要一个初始电荷密度所以静态计算目录里的 CHGCAR 一定要拷贝过来。NBANDS 这里我建议检查一下。默认的NBANDS可能只包含占据带和少量空带但 unfold 后想看清导带结构空带要留够。我一般会在OUTCAR里看实际自洽时的 NELECT 和 NIONS然后指定NBANDS 128具体数值取决于体系电子数和要不要算高能导带。宁可多算一些后面画图时裁掉不用也不要因为空带不够重跑。5. unfold能带展开实操细节5.1 用BandUP展开超胞能带BandUP 是我个人用得最多也最放心的 unfold 工具。它的流程分两步前处理生成 KPOINTS 和变换矩阵后处理算谱权重。前提条件是准备好两个文件POSCAR_prim参考单胞的晶胞和原子位置。POSCAR_sup掺杂超胞的晶胞和原子位置。前处理时BandUP 会比对两个结构的原子列表算出从单胞到超胞的变换矩阵确定每个超胞 k 点能映射到单胞里的哪个 k 点并按你的高对称路径在超胞倒空间中生成合适的 KPOINTS。以我的版本为例交互式命令大概是这样bandup_pre然后按提示输入 primitive cell 和 supercell 的文件名再输入要计算的路径点比如 G-X-W-L。如果没有交互模式参考它的 README 用命令行参数传入。前处理结束后目录里会出现一个新的 KPOINTS 文件和变换矩阵文件。把这个 KPOINTS 复制到能带计算目录跑完 VASP 之后再执行后处理bandup_post它会读取 VASP 的输出波函数信息和本征值输出一堆带权重的能带数据文件。这些文件的格式一般是每一行对应一个 k 点一个能量值以及一个谱权重。不同版本的 BandUP 命令名可能略有不同有的叫BandUP有的分pre_process和post_process。建议运行前先ls看看安装目录里的可执行文件或者bandup --help。5.2 用VASPKIT的unfold流程如果你的 VASP 实在版本低、BandUP 又装不上VASPKIT 的 unfold 模块是另一个选择。我在自己机器上也用过方便是方便但可控性差一点。操作思路如下准备好参考单胞的POSCAR_prim和掺杂后的POSCAR_sup。进入能带计算目录运行vaspkit。在电子结构相关菜单里选择 “Unfold band structure” 选项输入单胞和超胞文件。VASPKIT 会生成一个针对超胞路径的 KPOINTS 文件替换掉原来的 KPOINTS。跑 VASP。再次运行 VASPKIT选择 unfold 后处理得到可绘制权重能带的.dat文件。因为 VASPKIT 版本迭代很快菜单位置会有变化。遇到找不到的情况用vaspkit -h看任务列表里的 task name搜索 “unfold” 关键字最直接。VASPKIT 的优势是跟 VASP 结合紧密生成的 KPOINTS 一般不会出错。缺点是如果你想调 unfold 算法里的细节比如投影基底的选择或布里渊区划分就不如 BandUP 灵活。5.3 权重图怎么读颜色与线宽unfold 下来最重要的结果就是谱权重有的工具叫 spectral weight有的叫W(k, E)。这个值的物理含义是超胞本征态中包含的参考单胞波函数成分有多大。画图的时候一般把 k 点放在横轴能量放在纵轴每个点的颜色或大小对应该点权重。常见两种画法散点图权重越高点越大。热力图权重越高颜色越亮。我实际读图的经验是权重接近 1 的带就是“宿主能带”可以直接和单胞能带对应权重很低但仍出现在本征值列表里的带大多是折叠产生的假带不参与物理讨论掺杂引入的杂质态通常在带隙里显示为一条权重不高的窄带或平带说明它来自超胞里被破缺的局域原子而不是单胞的 Bloch 态。有一个实操误区看到杂质能级权重低就觉得可以忽略。实际上浅能级杂质态在带边附近会和宿主能带杂化权重有分布不能只看峰值。要找杂质态的位置配合LORBIT 11的投影数据看掺杂原子的 p 轨或 d 轨成分比单独看权重更稳。6. 常见问题与排错经验6.1 能带出现“断带”或杂乱点我碰到最多的问题就是 unfold 完的能带图里出现一堆不连续的短线。排查的思路一般是这样的先确认 KPOINTS 是不是 unfold 工具生成的。手动写的超胞 KPOINTS 经常因为坐标换算错误导致 k 点路径不对unfold 出来就是断的。再确认能带计算的CHGCAR是否来自静态计算。如果直接用了优化过程的中间电荷密度本征值会有噪声权重分布也跟着乱。最后看是不是空带不够。有些高能导带在超胞里发生了反交叉NBANDS 太少会截断导致只有部分能带可展开。如果只想看带边附近的谷位置断带问题不大但如果你想画完整带结构必须保证 k 点路径顺滑、NBANDS 足够。6.2 权重与投影不一致有时候 unfold 出来的带看起来色散没问题但杂质原子的投影权重却对不上。比如明明替位的是 P 原子LORBIT投影却显示 P 附近没有态。这种问题多半是原子索引错位。替位掺杂后P 原子在 POSCAR 中排在最后一位但超胞里其它原子的顺序和参考单胞不完全一致导致 unfold 工具在映射原子时匹配错位。解决办法是构造超胞时统一用同一套基向量不要混用惯用晶胞和原胞在跑 unfold 前先对比参考单胞和超胞的前几个原子坐标确保原子物种顺序一致。还有一个小坑如果用了非共线的磁性体系或者自旋极化体系需要分别对 spin-up 和 spin-down 做 unfold。不要把两个自旋通道的能量混在一起画否则看起来权重永远对不上。6.3 掺杂后费米能级漂移掺杂体系的费米能级往往不在带隙中间施主掺杂会把费米能级推向导带底附近受主掺杂则推向价带顶。这在 unfold 图里很直观但要注意一个问题VASP 的 Kohn-Sham 能级本身没有绝对意义掺杂前后的费米能级不能直接比较最好通过一个深层参考态对齐。我在实践中常用的做法是在超胞里选一个远离掺杂位点的原子记录它的最低占据芯能级比如 s 半芯态。把掺杂前后的芯能级对齐到同一数值。再比较带边或者杂质能级的相对变化。如果不做这一步很容易得出“掺杂使能带整体移动了 0.3 eV”的错误结论。对于替位掺杂这种误差完全可以避免。还有一个很实在的建议纯 DFT 会低估带隙杂质能级的位置可能不准。如果你要讨论深能级准确位置最好在 unfold 流程之后再用杂化泛函 HSE 重算一次筛选关键掺杂态。HSE 的计算贵但至少可以用 unfold 算出来的候选能级角度先做一轮筛选挑重点体系算耗时会可控很多。我自己跑这套流程最深的体会是unfold 不是一个可有可无的后处理而是掺杂超胞能带分析的“标配”。刚开始图快跳过前处理直接用单胞 KPOINTS 算折腾两天发现能带全乱后来老老实实把单胞、超胞、变换矩阵准备好再交给 BandUP 展开半小时就能出干净的图。如果你也正准备做替位掺杂的能带分析建议按这个顺序走先想清楚参考单胞再严格生成超胞优化后务必静态自洽最后算能带前就让 unfold 工具生成 KPOINTS。这样一轮跑下来比反复试错省下的时间要多得多。