sisl：提取DOS和轨道投影态密度

DOS需要对整个布里渊区加权求和，不能用高对称能带路径代替。以下衔接模型能带教程，用24×24×1完整k网格和0.08 eV高斯sigma。

1. 生成完整数据#

module load sisl/0.16.4
python calculate.py
python analyze.py
head -6 analysis/dos.csv

analysis/dos.csv依次为能量、总DOS、A轨道PDOS、B轨道PDOS。单自旋块每个原子一个轨道，本模型每胞共2态，没有另乘自旋简并。

2. 原生DOS与PDOS操作#

下面代码接在前文H的定义之后：

import numpy as np
from sisl.physics import MonkhorstPack,get_distribution
E=np.linspace(-9,9,2001)
bz=MonkhorstPack(H,[24,24,1],trs=False)
distribution=get_distribution("gaussian",smearing=0.08)
dos=np.zeros_like(E)
pdos=np.zeros((2,len(E)))
for k,w in zip(bz.k,bz.weight):
 state=H.eigenstate(k)
 dos+=w*state.DOS(E,distribution=distribution)
 p=state.PDOS(E,distribution=distribution)
 assert p.shape==(1,2,len(E))
 pdos+=w*p[0]

本版PDOS多一个显式自旋轴，当前形状为(1,2,NE)。先检查再取第0个自旋块，避免把自旋轴当成轨道。trs=False保留本例完整576个k点，权重总和为1。

3. 核对态数与投影#

python - <<'PY'
import numpy as np
a=np.loadtxt("analysis/dos.csv",delimiter=",",skiprows=1)
print("DOS总态数:",np.trapezoid(a[:,1],a[:,0]))
print("A轨道:",np.trapezoid(a[:,2],a[:,0]))
print("B轨道:",np.trapezoid(a[:,3],a[:,0]))
print("投影求和最大误差:",abs(a[:,2]+a[:,3]-a[:,1]).max())
PY

本次实际积分总态数2.0，A/B各1.0，投影求和与总DOS一致。另用解析能带直接做高斯求和，和原生DOS的最大差约9.16×10^-16。这里核对的是同一模型和同一网格的算法一致性，不能据此宣称真实材料DOS已收敛。

4. 看图并调整参数#

cat analysis/SUMMARY.json
ls figures/graphene-bands-dos.png

DOS图位于能带教程右侧。当前网格与展宽下仍有明显离散网格起伏；研究计算需要增加网格、改变展宽作对照，并保持能量参考与态数归一化一致。不要把人为展宽或有限网格峰解释成材料的新能级。

下载可复跑脚本，新目录解压后执行bash run.sh。使用独立sisl/0.16.4模块，Python3.11、NumPy/SciPy/netCDF4与Matplotlib，不需要加载SIESTA MPI工具链。results/是原生输出，analysis/为CSV，figures/为三张图。本组是明确构造的紧束缚模型，没有运行SIESTA自洽计算、DFT或TBtrans输运。