（一）计算所需程序

1. Quantum ESPRESSO 5.x或6.x，自行下载并编译

2. plotband (https://gitee.com/yhli/misc/tree/master/qe/plotband)  用于处理能带数据

3. analydos (https://gitee.com/yhli/misc/tree/master/qe/analydos)  用于分析投影信息

若要结合QE 6.x使用，编译analydos时需在Makefile的FPPFLAGS中加入-DQE6参数。若用gfortran编译，FPPFLAGS中还需加入-cpp参数，用ifort编译则无需此参数。

(二) 计算能带和分波态密度

1. 用pw.x进行自洽计算(calculation="scf")
   

2. 用pw.x进行能带计算(calculation="bands")

3. 用bands.x提取能带数据

4. 用projwfc.x计算分波态密度

以画黑磷pz轨道的fatband为例，每一步对应输入文件如下：

scf.in:

&CONTROL
    prefix           = '0'
    calculation      = 'scf'
    restart_mode     = 'from_scratch'
    pseudo_dir       = './'
    outdir           = './'
    wfcdir           = './'
    verbosity        = 'high'
    wf_collect       = .TRUE.
/
&SYSTEM
    ibrav            = 0
    celldm(1)        = 1.889726125
    ntyp             = 1
    nat              = 4
    ecutwfc          = 60
/
&ELECTRONS
    electron_maxstep = 250
    conv_thr         = 1.0d-10
    mixing_mode      = 'plain'
    mixing_beta      = 0.7
    mixing_ndim      = 8
    diagonalization  = 'david'
    diago_david_ndim = 4
    diago_full_acc   = .TRUE.
/
ATOMIC_SPECIES
   P  30.97   P.pbe-mt_fhi.UPF
CELL_PARAMETERS {alat}
   3.300150259   0.000000000   0.000000000
   0.000000000   4.617311378   0.000000000
   0.000000000   0.000000000  20.000000000
ATOMIC_POSITIONS {crystal}
   P   0.799366210   0.081375072   0.572085817
   P   0.798588284   0.259277753   0.466862505
   P   0.298663408   0.581549059   0.466844226
   P   0.299372098   0.759058116   0.572047452
K_POINTS {automatic}
 10 8 1 0 0 0

bands.in:

&CONTROL
    prefix           = '0'
    calculation      = 'bands'
    restart_mode     = 'from_scratch'
    pseudo_dir       = './'
    outdir           = './'
    wfcdir           = './'
    verbosity        = 'high'
    wf_collect       = .TRUE.
/
&SYSTEM
    ibrav            = 0
    celldm(1)        = 1.889726125
    ntyp             = 1
    nat              = 4
    ecutwfc          = 60
    nbnd             = 60
/
&ELECTRONS
    electron_maxstep = 250
    conv_thr         = 1.0d-10
    mixing_mode      = 'plain'
    mixing_beta      = 0.7
    mixing_ndim      = 8
    diagonalization  = 'david'
    diago_full_acc   = .TRUE.
    startingwfc      = 'random'
/
ATOMIC_SPECIES
   P  30.97   P.pbe-mt_fhi.UPF
CELL_PARAMETERS {alat}
   3.300150259   0.000000000   0.000000000
   0.000000000   4.617311378   0.000000000
   0.000000000   0.000000000  20.000000000
ATOMIC_POSITIONS {crystal}
   P   0.799366210   0.081375072   0.572085817
   P   0.798588284   0.259277753   0.466862505
   P   0.298663408   0.581549059   0.466844226
   P   0.299372098   0.759058116   0.572047452
K_POINTS {crystal_b}
   5
   0.0   0.0   0.0  40   ! Gamma
   0.5   0.5   0.0  40   ! M
   0.5   0.0   0.0  40   ! X
   0.0   0.0   0.0  40   ! Gamma
   0.0   0.5   0.0  40   ! Y

bands.x.in:

&BANDS
    prefix     = '0'
    outdir     = './'
    filband    = 'bndstr.dat'
    no_overlap = .TRUE.
/

pdos.in:

&PROJWFC
    prefix  = '0'
    outdir  = './'
    ngauss  = 0
    degauss = 0.005
    Emin    = -30
    Emax    =  30
    DeltaE  = 0.05
    filpdos = '0'
    filproj = 'proj.dat'
/

(三) 后期处理

1. 用plotband.x处理能带数据，得到bndstr.pd（也可以取其它名字）

2. 打开proj.dat（如果用QE 6.x计算，文件名为proj.dat.projwfc_up），删掉文件头，保留如下内容（第一个数字为计算投影所用的原子轨道数目，第二个为K点数目，第三个为能带数目）

3. 新建一个pz.inp，内容如下：

'proj.dat'
1 161
1 60
1
P 2 1 1

其中第一行是投影文件名称。第二行是要分析的K点编号最小值和最大值，第三行是要分析的能带编号最小值和最大值，可以取所有能带和所有K点。第四行表示要分析一种原子轨道的贡献，第五行表示该原子轨道为P原子的pz轨道，具体取值根据赝势文件确定：

    Valence configuration:
    nl pn  l   occ       Rcut    Rcut US       E pseu
    3s  3  0  2.00      0.000      0.000     0.000000
    3p  3  1  3.00      0.000      0.000     0.000000
    3d  3  2  0.00      0.000      0.000     0.000000
    4f  4  3  0.00      0.000      0.000     0.000000

由赝势文件中给出的信息可知，2对应3p轨道，第一个1表示角量子数l=1，第二个1表示与l=1对应的第一个p轨道。根据projwfc.x的帮助文档，这个轨道就是pz。如果要计算px+py轨道对应的贡献，输入文件应该这样写：

'proj.dat'
1 161
1 60
2
P 2 1 2

P 2 1 3

然后运行analypdos.x: analypdos.x pz.inp pz.out，输出文件pz.out格式如下：

4. 运行analydos_pp.py，将bndstr.pd和pz.out拼成便于画图的格式：analydos_pp.py pz.out bndstr.pd pz.xyz，输出文件pz.xyz格式如下：

第一列是K路径长度，第二列是能量，第三列是pz轨道对能带的贡献。取第一列和第二列作图，即可得到能带，再以第三列为权重，即可得到fatband。fatband可以用Origin画，可以用matlab或matplotlib的scatter函数画。下面是用gnuplot画fatband的脚本：

reset

# terminal settings
width_inch = 4
height_inch = 3

font_file = "font 'Arial, 10' fontscale 1.0"

# shared settings for all figures
width_line = 2
width_border = 2
width_arrow = 2

# band structure settings
x1 = 0.000000
x2 = 0.619194
x3 = 0.979243
x4 = 1.482994
x5 = 1.843043

y1 = -15
y2 = 5
dy = 5
nmytics = 5

symbol_G = '{/Symbol G}'

# setup terminal
set term pdfcairo size width_inch, height_inch @font_file
set output 'bp.pdf'

# plot band struture
unset key

set border lw width_border

set xrange [x1:x5]
set yrange [y1:y2]

set xtics nomirror  (symbol_G x1, 'M' x2, 'X' x3, symbol_G x4, 'Y' x5)
set ytics dy nomirror out
set mytics nmytics
unset y2tics

set ylabel 'Energy (eV)' offset char 1

set arrow from x2, y1 to x2, y2 nohead lw width_arrow
set arrow from x3, y1 to x3, y2 nohead lw width_arrow
set arrow from x4, y1 to x4, y2 nohead lw width_arrow

plot 'pz.xyz' u 1:2:3 w lp lw width_line lc 'blue' ps var pt 7

# flush output
set output

最终效果如下：

将不同轨道对应的fatband画在一张图中，便可得到在文章中经常见到的五颜六色的fatband图。但这种图有一种致命缺陷，就是如果两种原子轨道对同一条能带贡献相近，后来画的轨道就会覆盖前一个轨道，导致错误的结论。