sisl：读取实际SIESTA电子结构输出

前面的sisl模型教程使用显式构造的矩阵。本篇改为读取本次SIESTA Si2实算产生的真实文件，解释文件精度和单位，而不把模型文件冒充DFT结果。

1. 先得到Si2输出#

下载文末Si2资料包、在新目录解压：

sbatch job.slurm
squeue -u "$USER"
ls cases/dzp150/si.{XV,EIG,KP,HSX,RHO}

在作业完成后再读取。本例程序SIESTA5.4.2，后处理sisl0.16.4，DZP每个Si13个轨道，Si2共26个轨道，原生k点表有44条记录。

2. 读坐标、本征值和密度#

module load sisl/0.16.4
python - <<'PY'
import sisl
d="cases/dzp150/"
g=sisl.get_sile(d+"si.XV").read_geometry()
e=sisl.get_sile(d+"si.EIG").read_data()
rho=sisl.get_sile(d+"si.RHO").read_grid()
print("atoms",g.na,"cell Angstrom",g.cell)
print("EIG spin,k,band",e.shape)
print("valence electrons",rho.grid.sum()*rho.dvolume)
PY

实际EIG形状为(1,44,26)；RHO积分约 8.00000016，对应两个Si的8个价电子，不包含赝势去掉的芯电子。密度截面与实际结构见入门教程配图。

3. 正确读取当前HSX#

当前Intel Fortran写HSX的逻辑真值为−1，sisl0.16.4默认精度标志读取会遇到该值。本包只对当前读取对象增加一个兼容子类，原文件和公共软件都保持不变：

"""Process-local reader for Intel Fortran's -1 logical true in own HSX output."""
import numpy as np
from sisl.io.siesta.binaries import hsxSileSiesta,_siesta
class IntelLogicalHSX(hsxSileSiesta):
    @property
    def _stored_dtype(self):
        value=int(_siesta.read_hsx_is_dp(self.file))
        assert value in [-1,0,1],('Unexpected HSX precision logical',value)
        return np.float64 if value!=0 else np.float32

python - <<'PY'
import runpy
reader=runpy.run_path("hsx-reader.py")["IntelLogicalHSX"]
H=reader("cases/dzp150/si.HSX").read_hamiltonian()
print(H.no,H.eigh([0,0,0]))
PY

这只处理精度字段的−1/0/1逻辑约定，不改变Hamiltonian数值，也不放宽电子收敛要求。

4. 比较HSX与EIG时注意k点文本精度#

python - <<'PY'
import sisl,numpy as np,runpy
d="cases/dzp150/"
H=runpy.run_path("hsx-reader.py")["IntelLogicalHSX"](d+"si.HSX").read_hamiltonian()
bz=sisl.get_sile(d+"si.KP").read_brillouinzone(H.geometry)
E=sisl.get_sile(d+"si.EIG").read_data()[0]
k=np.rint(bz.k*4)/4
reconstructed=np.array([H.eigh(q) for q in k])
print("max difference eV",abs(reconstructed-E).max())
PY

KP文本坐标只有有限有效位数，直接读回的高能带可出现约10^-4 eV差异。这里已知原输入是4×4×4网格，因此按该输入重建同一分数k点；DZP150误差降至约3.99×10^-7 eV，原KP文件完整保留。未知网格或非均匀路径不能直接照搬这个四舍五入操作。

完整脚本还独立解析EIG文本并逐项核对，检查RHO电子数、最终晶胞和受力。运行后看analysis/cases.csv、analysis/SUMMARY.json与三张图。

下载Si2实算资料包，在新目录解压，执行sbatch job.slurm申请CPU资源。分区按56核整节点分配；本例小任务每条件2个MPI进程，最多四条件并发。参考数据同样采用2个MPI进程／条件。软件为SIESTA5.4.2 CPU MPI（oneAPI2024.2、IntelMPI2021.13、MKL LP64顺序库、AVX512/O2），后处理sisl0.16.4。包内无软件程序或赝势，准备步骤读取公共PseudoDojo0.4 PBE标量相对论Si PSML，并记录哈希。