Quantum Espresso中PDOS计算与轨道分析全流程详解

做第一性原理计算，能带图和态密度图基本是标配。上一期我们聊了能带计算，这一期把话题转到Quantum Espresso里的PDOS计算与轨道分析。简单说，总态密度只能告诉你某个能量区间里有多少电子态，而PDOS会把这一大坨态拆开，明确告诉你它们分别来自哪个原子、哪种轨道（s、p、d还是f）。这个信息在做催化、半导体、电池材料分析时几乎是刚需，比如判断价带顶是金属d轨道还是氧p轨道主导，吸附前后哪条轨道发生了杂化，都要靠它说话。

适合看这篇内容的朋友分两种：一种是用QE跑过基本SCF，但代码只停留在“跑完能算个总能量”的阶段；另一种是已经在其他软件里画过态密度，想换到QE里把投影流程彻底搞清楚。文章里的命令和输入文件都是可以直接抄作业的，我另外会把我踩过的坑一并写出来，至少在“为什么这么设参数”这件事上，尽量让你看完之后就通透了。

1. 为什么要算PDOS：从总态密度到轨道投影

1.1 总DOS能告诉我们的和不能告诉我们的

态密度（DOS，Density of States）的物理含义很直白：在某个能量附近，单位能量区间内有多少个可占据的电子态。它直接反映了体系的电子结构特征，比如半导体带隙、金属的费米面平台、范霍夫奇点，这些都能从总DOS上看出来。

但总DOS有一个先天短板：它把整个体系所有原子、所有轨道的信息全揉在一起了。你看着总DOS里一个很高的峰，根本不知道它是第几个原子贡献的，更别说是哪种轨道了。以前有朋友拿总DOS的峰位去硬套文献里的轨道归属，套了半天对不上，后来用PDOS一算才发现，原来那个峰主要是表面原子贡献的，体相原子根本没有参与。这就是总DOS在“归属问题”上的局限性。

PDOS，全称Projected Density of States，也就是投影态密度，正好解决这个归属问题。它的处理方式是对每个本征态按原子位置和角动量做投影，然后看每个投影分量的态密度分布。这样你就能得到“第1个Mo原子的d轨道态密度”“第2个S原子的p轨道态密度”，每一层都拆得明明白白。

1.2 PDOS的物理图像：把“乐手”单独拎出来

用一个生活化的类比来理解PDOS：总DOS就像一场演出的现场录音，听起来很热闹，但你分不清吉他和贝斯各自弹了些什么；PDOS则相当于给每个乐手单独开了一个音轨，想听哪个就调哪个。你要判断这首歌的灵魂是吉他的扫弦还是键盘的铺底，直接切对应音轨就行。

在代码层面，PDOS的计算基于一个简单的投影权重：

[ P_{\alpha,l,m}(E) = \sum_{n} \int_{BZ} d\mathbf{k} \left|\langle \phi_{\alpha,l,m} | \psi_{n\mathbf{k}} \rangle\right|^2 \delta(E - \varepsilon_{n\mathbf{k}}) ]

通俗解释一下：(|\psi_{n\mathbf{k}}\rangle)是第(n)条带在(\mathbf{k})点的本征态，(|\phi_{\alpha,l,m}\rangle)是以第(\alpha)个原子为中心的投影基函数（带角动量(l)和磁量子数(m)）。括号里的投影权重越大，说明这个本征态越像该原子的那条轨道。把所有本征态的贡献按能量排开，就得到了对应轨道的PDOS。

实际操作中，我们不会逐个态去手算这个权重，Quantum Espresso里的 projwfc.x 就是专门干这个的。它读入nscf步骤产生的波函数，完成投影，再按能量统计，输出一系列文本文件。画图的时候把这些文件拉出来，该填色填色，该叠加叠加，一张可以放进论文的PDOS图就出来了。

1.3 工具链梳理：pw.x 与 projwfc.x 的分工

在Quantum Espresso里做PDOS，标准流程是三步：先用 pw.x 做SCF自洽得到电荷密度，然后用 pw.x 做一次nscf计算，在一个更密的k网格上展开波函数，最后用 projwfc.x 读取nscf的波函数做轨道投影。如果想把同一元素不同原子的贡献合并，还需要用到 projwfc.x 家族的后处理工具 sumpdos.x 。

我见到不少第一次上手的人会问：为什么不能直接拿SCF的结果算PDOS？为什么非得再跑一次nscf？这个问题很关键，下一节专门解释。

另外说一下工具对比。用过VASP的人可能知道，VASP里看PDOS需要LORBIT和ISMEAR配合，然后从 DOSCAR 里自己解析轨道分量。QE的做法其实更灵活， projwfc.x 直接给出原子-角动量-磁量子数逐层分解的文件，还能配合 filproj 输出投影矩阵供后续Wannier或fatband分析使用。缺点是有一定的学习成本，文件多，命名规则需要适应。这篇文章后面会把这些文件规则捋清楚。

2. 计算思路拆解：SCF之后为什么还要多算一次nscf

2.1 SCF自洽只是入场券

很多教程会把PDOS的输入文件直接丢给你，但不说为什么。我先把这个关键问题讲透：SCF自洽计算的核心目标是让电荷密度自洽收敛，它只需要保证总能量和电荷密度达到精度要求就行，所以k点网格通常不会开得特别密。对于大部分体系，一个中等密度的k网格就能让总能量收敛得很好。

但态密度是k空间分布的统计结果，它对k点采样的密度非常敏感。打个比方，SCF的k网格像是在地图上每隔10公里设一个观测站，总能量这种“大尺度平均量”用这些站点就能估出来；而DOS需要画出能量区间的细节，10公里一个站点太稀了，山谷和山峰可能直接被漏掉。所以必须在nscf步骤里把k网格加密，通常比SCF加密一倍甚至更多，这样态密度曲线才会平滑、可靠。

nscf和scf的区别在于：nscf不再自洽迭代，它利用SCF已经收敛的电荷密度，固定住这个势场，只在一个新的k点网格上重新求解Kohn-Sham方程，得到本征值与本征函数。由于不需要反复迭代，即使k点数增加很多，计算代价也可控。这步就是为PDOS和能带计算提供“高质量波函数库”的。

2.2 nscf的输入文件与关键参数

直接给一个我在MoS2单层上常用的nscf输入文件，细节我逐条解释。MoS2是典型的二维半导体，价带顶主要包含Mo 4d和S 3p成分，导带底几乎被Mo 4d主导，用来演示轨道分析非常经典。整个体系3个原子，计算量不大，适合拿来跑通全流程。

&CONTROL
   calculation = 'nscf'
   prefix      = 'mos2'
   outdir      = './out'
   pseudo_dir  = './pseudo'
/
&SYSTEM
   ibrav       = 0
   nat         = 3
   ntyp        = 2
   ecutwfc     = 60
   ecutrho     = 480
   occupations = 'tetrahedra'
   nbnd        = 40
/
&ELECTRONS
   conv_thr    = 1.0d-8
/
ATOMIC_SPECIES
Mo 95.94 Mo.pbesol-spn-rrkjus_psl.1.0.0.UPF
S  32.06  S.pbesol-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (angstrom)
  3.180000   0.000000   0.000000
 -1.590000   2.754000   0.000000
  0.000000   0.000000  15.000000
ATOMIC_POSITIONS (crystal)
Mo  0.000000  0.000000  0.000000
S   0.333333  0.666667  0.105000
S   0.666667  0.333333 -0.105000
K_POINTS (automatic)
24 24 1 0 0 0

几个重点说一下。

第一， calculation = 'nscf' ，这是最核心的一行，告诉pw.x不要做自洽迭代，直接用已有的电荷密度去解本征问题。运行这个输入文件之前，必须先在同一 outdir 下跑完SCF，否则它会报错找不到电荷密度文件。

第二， occupations = 'tetrahedra' 。这一步是让态密度更平滑的关键。QE使用四面体方法对布里渊区积分做处理，能有效避免直接把离散k点上的态用高斯函数涂抹带来的误差。如果体系是金属，SCF阶段通常用 smearing ，nscf阶段也推荐用 tetrahedra 或保持 smearing ，具体看你想在DOS里保留多大的展宽控制权。但要提醒一句： tetrahedra 要求k点网格是完整的、规则分布的，自动生成的Monkhorst-Pack网格没问题。

第三， nbnd = 40 。这个参数经常有人忽略。QE默认的能带数只保证覆盖占据态再加一点空带，但PDOS要看的能量范围可能远高于费米能级，尤其是想看导带较高区域的轨道成分时，默认的 nbnd 根本不够。一个粗略的经验是：先看SCF输出的最高占据带序号，然后在此基础上增加10到20条空带。我这里设40，是给MoS2留了比较充足的导带空间。

第四，k点网格。SCF阶段我用的是 12 12 1 ，这里改成 24 24 1 ，相当于在二维布里渊区每个方向加密了一倍。如果发现PDOS曲线还是带毛刺，可以继续加到 36 36 1 或者更密。二维材料只需要一个k点沿z方向，因为z方向被真空层隔开，色散很小。

2.3 赝势、prefix与outdir的一致性

这一小节看起来是“常识”，但我在各种教学场合见过太多人栽在这里，所以值得单独强调：nscf和scf的输入文件里， prefix 、 outdir 、赝势文件、 ecutwfc 、 ecutrho 、晶格参数、原子坐标，必须完全一致。

prefix 如果写错了，pw.x会去另一个前缀的文件目录里找电荷密度，结果必然是找不到。 outdir 如果不同，比如一个写 ./out 一个写 ./tmp ，也一样抓瞎。赝势不一致更隐蔽，SCF用A赝势，nscf换成了B赝势，程序通常不会立刻报错，但算出来的波函数和电荷密度根本不是同一个体系，出来的PDOS就是张冠李戴。

还有一个容易踩的坑是：如果体系有磁性（比如自旋极化计算），nscf里 nspin 的设置要和SCF保持一致。SCF里开了 nspin = 2 ，nscf里忘记写，程序会按非磁处理，报错或者结果完全不对。

这里多说一句，虽然我给的示例是在本机用串行或单节点mpi运行，但实际生产环境通常需要并行。命令格式一般是 mpirun -np 16 pw.x -in nscf.in > nscf.out 。并行跑的时候要注意QE对不同k点的并行划分，只要输入文件没问题，这部分很少出幺蛾子。

3. projwfc.x实操：从命令到输出文件

3.1 一个可以直接用的projwfc输入文件

SCF和nscf都跑完之后，进入今天的主角： projwfc.x 。它的输入文件比pw.x简洁得多，最核心的内容是告诉它波函数放在哪、要投影哪些能量范围、用什么展宽方式。下面是我在MoS2例子上用的文件。

&PROJWFC
   prefix   = 'mos2'
   outdir   = './out'
   filpdos  = 'mos2_pdos'
   Emin     = -20.0
   Emax     = 10.0
   DeltaE   = 0.05
   ngauss   = 0
   degauss  = 0.02
/

运行命令一行就够：

mpirun -np 4 projwfc.x -in projwfc.in > projwfc.out

参数逐个解释。

prefix 和 outdir 不用多说，和nscf保持一致就好。 filpdos 是输出文件的前缀，我习惯写成 mos2_pdos ，这样后面生成的所有文件都带这个前缀，方便管理。

Emin 和 Emax 是能量窗口，单位是电子伏特。这里我设成-20到10 eV，意思是只计算费米能级以下20 eV到以上10 eV这个范围内的投影态密度。为什么要设这个窗口？因为PDOS文件是逐能量点输出的，范围越大，计算量和文件体积越大。如果你只关心价带顶和导带底附近的轨道成分，窗口可以再收窄，比如-5到5 eV。

DeltaE = 0.05 是能量步长，单位也是eV。这个值决定最后曲线的分辨率，0.05 eV已经比较细腻，画出来的图足够写论文。如果你只是快速看一眼趋势，用0.1 eV也行，文件会小一半。

degauss = 0.02 是高斯展宽宽度，注意单位是Ry。0.02 Ry大约等于0.27 eV，这个展宽能把离散能级抹成比较光滑的曲线。如果是金属体系或者k点不够密，可以适当加大到0.03甚至0.05 Ry，代价是峰会被压扁，细节丢失。 ngauss = 0 表示使用普通高斯展宽函数。对大多数体系，这个配置比较稳。

运行结束后， projwfc.out 尾部会有一段统计信息，告诉你总共有多少个能级被投影、每个原子每种角动量的投影总和是多少。其中有个指标值得看一眼：每个原子上各角动量投影之和，应该大致等于该原子的价电子数。如果偏差很大，说明投影半径、赝势或者k点可能有问题。

3.2 输出文件里的每一列是谁

跑完 projwfc.x 之后，目录下会冒出一堆文件，命名规则大致是这样的：

mos2_pdos_tot.dat
mos2_pdos_atm#1(Mo)_wfc#1(s).dat
mos2_pdos_atm#1(Mo)_wfc#2(p).dat
mos2_pdos_atm#1(Mo)_wfc#3(d).dat
mos2_pdos_atm#1(Mo)_wfc#4(f).dat
mos2_pdos_atm#2(S)_wfc#1(s).dat
mos2_pdos_atm#2(S)_wfc#2(p).dat
mos2_pdos_atm#2(S)_wfc#3(d).dat
...

用户名里面的 wfc#1 到 wfc#4 ，分别对应s、p、d、f这条角动量通道。每个原子都会生成这么一组文件，和原子在输入文件里的出现顺序一一对应。文件名里括号中的元素符号是代码帮你标注的，非常贴心。

打开一个文件看，比如Mo的d轨道投影文件 mos2_pdos_atm#1(Mo)_wfc#3(d).dat ，头部注释会写明每一列的含义，大致长这样：

# E (eV)    l=2 m=-2    l=2 m=-1    l=2 m=0    l=2 m=1    l=2 m=2
-20.0000   0.1234E-03  0.5678E-03  ...

第一列是能量，后面的列是同一个角动量通道下不同磁量子数m分量的投影态密度。也就是说，对每个原子，d轨道被进一步拆成了(d_{xy})、(d_{yz})、(d_{z^2})、(d_{xz})、(d_{x^2-y^2})这几个分量，具体对应关系在注释行里写得很清楚。

这里有个常见疑惑：如果我想画“Mo的d轨道总PDOS”，是不是直接用这个文件就行？不行。这个文件里的5个m分量是分开的列，你需要先把它们加起来。好在很多绘图脚本都支持直接按列求和，后面我给出一个现成的Python脚本。

mos2_pdos_tot.dat 是总态密度文件，第一列能量，第二列总DOS（如果是自旋极化，还会有第三列，对应另一个自旋方向）。注意：各原子PDOS加起来的和，并不严格等于总DOS，因为projwfc的投影只计算了原子球（Muffin-tin球内）的部分，间隙区域的电子态没有分配到任何原子头上。这个差异是正常的，不要拿它当错误去排查半天。

3.3 用sumpdos.x把同元素轨道累加

如果你的体系里某种元素有多个不等价原子，比如表面和体相的原子，直接看单个原子上某条轨道的PDOS可能看不出整体规律。这时候可以用 sumpdos.x 把多个原子的同一角动量通道合并起来，输出成一个文件。

下面是一个合并Mo所有s轨道和所有d轨道的输入示例（不同版本输入格式略有差异，以你安装目录下 Doc/INPUT_SUMPDOS.html 里的模板为准）：

&DOS
   Emin   = -20.0
   Emax   = 10.0
   DeltaE = 0.05
   ngauss = 0
   degauss = 0.02
/
2
1 1
3 3

解释一下：分段里设置的能量范围和展宽参数需要与projwfc一致，然后是合并数量2，下面两行分别表示“原子1的wfc#1（s通道）”和“原子3的wfc#3（d通道）”。在MoS2例子里，如果Mo是第1个原子，这个配置就是合并Mo的s轨道和d轨道，输出文件的格式和单独原子的一样，画图直接用。

实际用的时候，最常见的需求其实是把同一个元素所有原子的同一条轨道全部合并，比如把两个S原子的p轨道合并成“S-p”。那就把两个S原子的p通道都列进去，数量改成2，画出来就是整个S元素的p轨道贡献。

3.4 画一张能放进论文的PDOS图

PDOS画图我习惯直接用Python + matplotlib，读文本文件、按列求和、填色，一套组合拳下来非常顺手。下面这段脚本可以直接改着用：

import numpy as np
import matplotlib.pyplot as plt

def read_pdos(filename, ef=0.0):
    data = np.loadtxt(filename, comments='#')
    energy = data[:, 0] - ef
    dos = np.sum(data[:, 1:], axis=1)
    return energy, dos

ef = 0.0  # projwfc默认以费米能级为零点，这里直接用

e_mo_d, dos_mo_d = read_pdos('mos2_pdos_atm#1(Mo)_wfc#3(d).dat', ef)
e_s_p, dos_s_p = read_pdos('mos2_pdos_atm#2(S)_wfc#2(p).dat', ef)

plt.figure(figsize=(6, 4))
plt.plot(e_mo_d, dos_mo_d, label='Mo-d', color='crimson')
plt.fill_between(e_mo_d, dos_mo_d, color='crimson', alpha=0.3)
plt.plot(e_s_p, dos_s_p, label='S-p', color='steelblue')
plt.fill_between(e_s_p, dos_s_p, color='steelblue', alpha=0.3)
plt.axvline(0, color='black', ls='--', lw=0.8)
plt.xlabel('E - E_F (eV)')
plt.ylabel('PDOS (states/eV)')
plt.legend()
plt.xlim(-10, 5)
plt.tight_layout()
plt.savefig('pdos_mos2.png', dpi=300)

注意脚本里的 np.sum(data[:, 1:], axis=1) ，它把同一原子同一角动量下所有m分量求和了，得到的就是该角动量的总PDOS。如果你只想画某个特定的(d_{z^2})分量，就单独取那一列，不需要求和。

如果系统是自旋极化，文件里会有自旋上、自旋下两套数据，画图时通常一个画在上半部分、一个画在下半部分，或者用两条曲线区分。具体要看你的体系是否需要。

只画PDOS曲线还不够，我强烈建议把费米能级对齐检查一下。projwfc的输出默认已经把能量零点放在费米能级上，所以横轴在0附近就是费米能级的参考位置。但如果用了 sumpdos.x 合并，它不会自动帮你搬移零点，需要在脚本里手动减掉费米能量。费米能级本身可以从nscf输出文件里找，搜 Fermi 关键字即可。

4. 轨道分析：从数值到结论

4.1 看带边：MoS2的Mo-d与S-p杂化

拿到PDOS图之后，最关键的动作是“读图”。以MoS2为例，价带顶和导带底的轨道成分，直接决定了它的载流子输运特性和催化活性。

如果你把Mo-d和S-p两条PDOS曲线画在同一张图上，会看到几个非常典型的特征：价带顶附近Mo-d和S-p都有明显峰，而且在同一个能量区间内重叠；导带底几乎只有Mo-d的峰，S-p的贡献非常弱。这说明什么？说明MoS2的价带顶是Mo 4d和S 3p杂化形成的成键态，而导带底是以Mo 4d为主的反键态。

这个结论对应到物理图像上就是：光激发一个电子从价带到导带，本质上是在把电子从成键态转移到以金属d轨道为主的反键态，整个过程中硫的p轨道成分会发生变化。这也是为什么MoS2这类过渡金属硫化物在催化析氢反应中表现出色，因为S-p和Mo-d的杂化直接影响了表面活性位点的电子结构。

读PDOS时有一个小技巧：不要只盯着峰所在的能量位置，还要看看带宽。如果某条轨道在很宽的能量范围内都有分布，说明它和周围原子的轨道杂化得很厉害，电子比较“离域”。如果PDOS在某个能量区间出现又尖又窄的峰，说明这条轨道非常局域，类似原子能级，通常对应孤对电子或者弱相互作用的芯态区域。

4.2 PDOS与能带联动的分析习惯

PDOS单独看能说明“哪些轨道在哪些能量上有贡献”，但要说清楚一条能带是由哪种轨道组成的，最好把PDOS和能带图放在一起看。

具体操作上，先把能带计算跑出来（能带计算本身需要单独一次nscf，沿高对称路径设k点），然后把能带图的能量轴和PDOS图的能量轴对齐，通常以费米能级为0。看某一条能带在VBM处色散很平缓，对应到PDOS里是一个局域峰，说明这条带主要由某条局域轨道构成；如果某条带色散很大、带宽几个eV，对应PDOS里是一条很宽的分布，说明这条带是杂化离域带。

这里有一个通用的分析框架：把PDOS的最高占据峰和能带图的价带顶对齐，把最低空峰和导带底对齐，然后逐个能量区间比较轨道成分。比如你看到一个能量区间里同时有Fe-d和O-p的峰，而且位置基本重合，就可以初步判断Fe 3d和O 2p在发生杂化。如果再结合晶体场理论，还能进一步判断是σ型杂化还是π型杂化，这部分内容如果展开会写成长文，这里先点到为止。

4.3 投影权重的定量化使用

PDOS不只是用来“看图说话”，投影权重本身也可以做定量的分析。一个比较常用的量化方式是积分特定能量区间的PDOS，看看某种轨道在这个窗口里的相对占比。举个例子，你想知道MoS2价带项下面2 eV范围内，Mo-d和S-p各占多少贡献，可以这样算：

mask = (e_mo_d >= -2.0) & (e_mo_d <= 0.0)
area_mo_d = np.trapz(dos_mo_d[mask], e_mo_d[mask])
area_s_p = np.trapz(dos_s_p[mask], e_s_p[mask])
print('Mo-d contribution:', area_mo_d / (area_mo_d + area_s_p))

这个数值可以直接写进文章里，作为“Mo-d与S-p杂化强度”的一个直观指标。不过要提醒一句，QE的projwfc投影是在原子球内做的，球外电子态没有分配，所以这个占比是半定量的，适合同一个体系内部的对比。比如比较不同吸附构型、不同掺杂位置之间的轨道变化，同一套参数下比值是有参考意义的。

另外， projwfc.x 还支持输出 filproj ，把每个本征态的投影权重保存下来。这个文件是后面Wannier90计算或者fatband投影能带图的重要输入。如果你想画那种“能带上某一段加粗表示d轨道贡献大”的fatband图，数据基础就在这里，这一步很多教程都是另外开一篇讲的，这里提一下你至少知道有这么个出口。

5. 常见问题与排查技巧实录

5.1 报错速查表

我把实际使用中遇到最多的报错和问题整理成了一张表，建议收藏。

这几类问题里，第二类比较容易被忽视。 projwfc.x 不会因为你把Emin设得太靠近费米能级而报错，它只会默默输出一堆0，或者整条曲线都是平的。所以如果刚跑完发现PDOS全是零，第一时间去检查能量窗口。

5.2 PDOS太锯齿怎么办

锯齿是所有算态密度的人都会遇到的事。原因不外乎两个：k点采样不够密，或者展宽太小。

先说k点。态密度本质上是对布里渊区大量k点本征值的统计，k点太少，本征值就只是离散的一串点，展宽成曲线后自然处处是尖刺。解决办法很简单：把nscf的k网格加密。从 12 12 1 改成 24 24 1 如果还不够，继续往上加。二维材料因为z方向不需要密k点，所以加x-y方向的点基本不心疼时间。

再说展宽。 degauss 设置得太小，等同于把每个离散能级用很细的针状高斯去涂抹，结果当然不光滑。建议先从0.02 Ry起步，如果锯齿还是明显，可以尝试0.03或0.05 Ry。但要注意平衡：展宽越大，峰会被压得越矮，本该分开的两个峰可能被抹成一个。我一般的策略是：先固定k点在比较密的程度，然后调节degauss，直到曲线看起来光滑且保留关键峰型。

另一个不太常提的细节是 DeltaE 。如果能量步长设得太大，比如0.2 eV，即使DOS本身很光滑，采样点太少也会造成视觉上的“折线感”。想要精致一点的图，直接降到0.02到0.05之间。

5.3 费米能级对齐

做材料计算的人应该都吃过费米能级对不齐的亏。不同代码、不同版本、不同软件包里，能带图和DOS图横轴的零点定义可能不一样。QE的projwfc.x输出时默认会把能量零点设在费米能级上，但你用 sumpdos.x 合并之后，文件里的能量可能就变成绝对能量了。

稳妥做法是：每次画图前，都从nscf的输出文件里把费米能级抓出来。命令可以这样：

grep "Fermi" nscf.out

然后把你导入脚本的能量列减去这个值。这样即使你换了QE版本，输出文件格式有了变化，画图逻辑也不会崩。

对于半导体和绝缘体，费米能级严格来说应该在带隙中间，scf输出的费米能级在半导体体系里会有一点点不确定性，这会导致你画的图零点略有偏移。这种情况下，建议用价带顶作为参考零点，而不是程序输出的费米能级。具体做法是：先看总DOS或能带图确定VBM位置，再统一平移。

5.4 关于投影方法的几个提醒

用projwfc的PDOS，有一些默认前提需要心里有数。

projwfc的投影是在以原子为中心的球形区域里做的，球的半径由赝势或代码内部设定，而不是随便给一个值。这意味着，如果你关注的轨道在原子间有很大延展，比如某些共价体系的成键态，投影权重会有一部分落在球外没被算上，导致PDOS夸大了局域特征。这是投影方法的物理限制，不是QE的bug。

另外，PDOS曲线非常依赖赝势。不同赝势对“价电子”的定义不同，比如有的Mo赝势把4p半芯态算进价电子，有的不算。这会导致你在PDOS的低能区域看到完全不同的峰。所以做横向对比的时候，比如对比文献里的PDOS，一定确认对方用的赝势和能量范围。

最后，projwfc的投影权重本身没有唯一的标准定义，它是把平面波本征函数往原子轨道基函数上投影的结果。不同软件采用不同的投影基组，画出来的PDOS曲线形状可能有细节差异。所以PDOS更适合用来看趋势、看相对变化，而不是抠某个具体数值的绝对值。

5.5 我踩过的几个坑

第一个坑是跑完projwfc之后发现 pdos_atm#1(Mo)_wfc#3(d).dat 里前几列全为0，只剩高能区有一点数。查了半天发现是nscf里的 nbnd 设得太少，导带空带不够，高能区域的态压根没被算出来。把 nbnd 从24加到40再跑，曲线马上正常。从那以后我养成一个习惯：只要算PDOS，nscf的 nbnd 一定会在默认值基础上再加15到20条空带。

第二个坑是二维材料没加够真空层就跑去算PDOS。当时我只给了10 Å真空，PDOS的形状倒也算得出来，但导带底附近有明显的层间人为色散，多出来的虚假峰误导了我好几天。后来把真空层加到15 Å以上，再重新跑nscf和projwfc，问题消失。做二维体系的朋友，这个参数一定要重视。

第三个坑是使用 sumpdos.x 时没有确认版本差异。我手头有两个QE版本，一个接受的输入格式是 &DOS 分段，另一个则需要固定顺序读取数字，结果用旧模板套新版本，直接报错。后来每次用之前都先去查一下 Doc/INPUT_SUMPDOS.html ，把版本差异避开，就再没出过这种问题。

最后再分享一个实操习惯

前面把PDOS的计算流程和轨道分析讲得差不多了，最后分享一个我用了很久的习惯：每到一个新体系，我不会着急把 degauss 调到很光滑，而是故意先用偏小的展宽和偏密的k点跑一次，看原始数据里有哪些窄峰。这些窄峰往往是某个原子局域轨道的指纹，比如过渡金属的d-d跃迁峰、分子轨道的孤对电子峰。等我知道哪些峰真实存在之后，再调大展宽画正式图片。这样既不会因为展宽太大丢掉真实细节，也不会因为数据太锯齿而误判峰位。

在Quantum Espresso里做PDOS计算本身并不难，难点在于理解每一步设置背后的原因，以及会读输出文件。掌握这套流程后，你会发现从总DOS到轨道分解，其实就是多跑一次nscf、多运行一回projwfc的事。后面如果再碰到复杂体系，比如异质结、掺杂表面、磁性材料，这套流程基本不用改，换一下输入文件和赝势就能直接套用。