ABACUS H/S矩阵重建能带实操
读取原生CSR、恢复H(k)和S(k)、解广义本征值，并与真实ABACUS能级逐项比较。
阅读本篇目录
已经输出H(R)和S(R)之后，怎样把它们用起来？本篇使用集群真实Si双原子LCAO计算的原生矩阵，读取CSR、做傅里叶变换，再与原程序的42个能级比较。它是矩阵后处理实例，不是新增一次SCF，也不把矩阵导出当作已经完成机器学习训练或拓扑计算。
1. 先运行附件#
下载原生矩阵、输入、能级和分析脚本，解压进入附件根目录：
bash run.sh
只用一个CPU核心，无需GPU或新的DFT作业。脚本加载ase/3.29.0，利用ASE、NumPy、SciPy和Matplotlib处理数据；原始DFT采用ABACUS3.10.1、GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8和OpenBLAS0.3.30。work已存在时会停止，重复体验请重新解压。
head -8 work/eigenvalues.csv
native_eV是istate.info原生能级，reconstructed_eV是从H/S矩阵重新求出的能级。三组k点各14个能级，最大差约1.14×10⁻⁷ eV，误差包含矩阵文本输出精度。
左图是在本次固定H(R)、S(R)模型上沿Γ—L—X—Γ采样，横轴按各段采样点排列；不是重新做SCF或证明k网格已经收敛。右图只比较原始计算实际输出的三组k点、42个能级。
2. 自己计算时先输出矩阵#
在自己的LCAO输入中设置out_mat_hs2 1；本例还设置out_dm1 1，供下一篇检查电子数。先完成SCF，再读取OUT.tutorial/中的矩阵与istate.info。需要重算原例时，可下载原始完整ABACUS算例包，解压进入abacus-examples，运行：
python new-case.py density-matrices
cd work/density-matrices
sbatch run.sh
提交脚本中的分区和账户按自己的核时/套餐权限使用yskit调整。下面保留原例完整输入；ks_solver elpa负责局域轨道本征问题，不能把文件中出现device gpu理解成所有步骤都在GPU上执行。LCAO本例不靠out_wfc_pw输出PW系数，该项不是本篇矩阵重建的前提。
完整INPUT、STRU、KPT
INPUT#
INPUT_PARAMETERS
suffix tutorial
calculation scf
basis_type lcao
ecutwfc 40
scf_thr 1e-8
scf_nmax 150
ks_solver elpa
device gpu
precision double
cal_force 1
cal_stress 1
pseudo_dir ../../PP_ORB
orbital_dir ../../PP_ORB
out_chg 1
out_pot 1
out_wfc_pw 1
out_dm1 1
out_mat_hs2 1
out_mat_r 1
STRU#
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.0.upf
NUMERICAL_ORBITAL
Si_gga_8au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
10.26
LATTICE_VECTORS
0 0.5 0.5
0.5 0 0.5
0.5 0.5 0
ATOMIC_POSITIONS
Direct
Si
0.0
2
0 0 0 1 1 1
0.25 0.25 0.25 1 1 1
KPT#
K_POINTS
0
Gamma
2 2 2 0 0 0
3. 矩阵文件怎样读#
本例有26个局域轨道，即两个Si原子各13个轨道。原生文件为：
| 文件 | 本例R块数 | 内容 |
|---|---|---|
| data-HR-sparse_SPIN0.csr | 177 | H(R)，能量单位Ry |
| data-SR-sparse_SPIN0.csr | 177 | S(R)，无量纲 |
| data-DMR-sparse_SPIN0.csr | 183 | 实空间密度矩阵，下一篇使用 |
前三行是步数、矩阵维数、R块数。每块先写Rx Ry Rz 非零数；非零数大于零时，随后三行分别是数值、列索引、行指针。列索引从0开始，行指针长度为27，最后一个值必须等于非零数；零块没有这三行。不能把所有行直接交给numpy.loadtxt当成一个普通二维表。
傅里叶变换为H(k)=Σ H(R)exp(−2πi k·R)，S同样处理；这里k使用倒格矢分数坐标，R使用晶格整数坐标，所以有2π。不要将笛卡尔k坐标原封不动放进这个式子。格式与符号依据ABACUS官方矩阵说明，本例文件名按3.10.1保留；新版3.11的输出名和参数需另按版本确认。
4. 为什么必须把S一起带上#
NAO通常不是正交基，要求解H C = S C ε，而不是直接对H求普通本征值。附件核心操作如下：
from scipy.linalg import eigh
energy_Ry, coefficients = eigh(Hk, Sk)
energy_eV = energy_Ry * 13.605698
单位换算使用ABACUS3.10.1自身的常数，便于和原生输出逐项核对。脚本同时检查H/S的Hermitian关系、S正定、C†SC=I和本征方程残差；对H直接求普通本征值的差异也记录在结果里，不能用作该NAO基下的能级。
在附件根目录可查看结果：
python - <<'PY'
import json
r = json.load(open('work/RESULTS.json'))
print('orbitals:', r['orbitals'])
print('max band error/eV:', r['max_band_error_eV'])
print('wrong ordinary-eigenvalue errors/eV:',
r['wrong_ordinary_eigenvalues_max_errors_eV'])
PY
附件还做了两个数学对照：同时重排H和S的轨道顺序，能级不变；把H改成H+cS，所有能级平移c。它们说明矩阵读写与基底处理的含义，不是额外的材料计算。接着看密度矩阵、电子数与布居。