sisl 非正交矩阵与电子态分析

使用集群 sisl/0.16.4，环境为 Python 3.11、NumPy、SciPy、netCDF4 和 Matplotlib，不需要加载 SIESTA 的 MPI 编译环境。本文是明确给定参数的教学模型，在集群实际运行并与独立公式比较；不是新的材料 DFT 或 TBtrans 自洽计算。

下载本组完整脚本、原生矩阵、CSV和配图。解压后进入 sisl-core-model-examples：

module load sisl/0.16.4
bash results.sh
bash plot.sh

参考数据在 reference-results、CSV在 reference-data。要运行全部短模型，用 bash run.sh，新结果写入 work/results，已完成小项自动跳过；图由参考数据重画到 figures。脚本限制数学库为1线程。

建一个带重叠的链模型#

下面是完整的最小输入，保存后可直接运行：

cat > chain.py <<'EOF'
import sisl
import numpy as np
g = sisl.Geometry([[0,0,0]], sisl.Atom(6, R=2.01),
    lattice=sisl.Lattice([2,10,10], nsc=[3,1,1]))
H = sisl.Hamiltonian(g, orthogonal=False)
H.construct([[0.01,2.01], [(0.4,1.0),(-1.0,0.12)]])
k = [0.24,0,0]
state = H.eigenstate(k, dtype=np.complex128)
print('energy / eV:', state.eig)
print('S norm:', state.norm2())
print('Euclidean norm:', np.sum(abs(state.state)**2, axis=1))
print('Cartesian dE/dq / eV Angstrom:', state.derivative())
H.write('chain.nc')
EOF
python chain.py

晶格长度2 Å；nsc=[3,1,1]让实空间矩阵包含左右相邻胞。每原子一个轨道，位能0.4 eV、最近邻跃迁−1 eV；重叠矩阵的对角元1、最近邻0.12。这里的C只是模型位点标签，不能把这组参数当成真实碳链势函数。

非正交归一化为什么不是系数平方和1#

分数k对应相位θ=2πk：H(k)=0.4−2cosθ、S(k)=1+0.24cosθ，因此本例能量为 E=H/S。程序求解的是 Hc=ESc，正确归一化是 c†Sc=1。

实际17个k点中，k=0.48得到 E=3.1293515703 eV，S=0.7618924717，c†Sc≈1，但c†c=1.3125211722。不要为了“归一化”再把这个正确态强制缩放成欧氏范数1。

state.derivative()对笛卡尔倒空间变量求导，结果为eV·Å；不是对无量纲分数k的导数。速度还需要除以ℏ并换单位，不能把这个数组直接标成m/s。

python - <<'EOF'
import sisl
H = sisl.get_sile('chain.nc').read_hamiltonian()
shifted = H.copy()
shifted.shift(0.3)
print(shifted.eigh([0.24,0,0]) - H.eigh([0.24,0,0]))
EOF

本例每条能量都提高0.3 eV。非正交情况下对应 H→H+0.3S，不是只给H的实空间对角元加0.3。

从二聚体态计算IPR和密度矩阵#

本包还有一个正交二聚体，H的对角元为+0.2、−0.2 eV，非对角元−1 eV。实际能量为±1.0198039027 eV；低能态概率为0.4019419324、0.5980580676，state.ipr()=0.5192307692。

IPR使用系数定义 Σ|c|⁴/(Σ|c|²)²。这个例子特意使用正交基；对于非正交原子轨道基，不应未经说明就把系数IPR当成唯一的实空间局域性指标。

python - <<'EOF'
import numpy as np
x = np.load('reference-results/state-and-density/finite-dimer-matrices.npz')
print('H:', x['H'])
print('density:', x['density'])
print('electrons:', np.trace(x['density']))
print('band energy / eV:', np.trace(x['energy_density']))
EOF

本例让低能态占据2电子、高能态0电子。密度矩阵迹为2；能量密度矩阵迹为−2.0396078054 eV。它是指定占据的模型带能，不是含双计数修正的DFT总能。完整 DensityMatrix、EnergyDensityMatrix赋值命令在 run_core.py 的 state-and-density 部分。

自旋块与普通轨道索引#

spin='non-collinear'与spin='spin-orbit'决定内部矩阵存储格式。资料包对两者分别构造同一个 H=B·σ，B=(0.2,0.3,0.4) eV；两本征值为±0.5385164807 eV，自旋期望值为±B/|B|。

python - <<'EOF'
import sisl
for name in ['non-collinear','spin-orbit']:
    H = sisl.get_sile('reference-results/spin-blocks/'+name+'.nc').read_hamiltonian()
    state = H.eigenstate()
    print(name, state.eig, state.spin_moment())
EOF

本例没有计算真实材料的SOC常数或磁各向异性；仅验证自旋矩阵装配、复数非对角元和态分析。读取自己的DFT输出时，应保留其真实自旋类型，不能用重命名把非磁矩阵变成SOC计算结果。

稀疏矩阵扩胞后要看能带折叠#

本包 sparse-supercell 以正交链调用 H.tile(2,0)，轨道数变2、胞长变4 Å；折叠后两条带分别对应原链k/2和k/2+1/2。untile(2,0)可还原这个规则超胞。

H.sub([0])会同时删去与其他被移除轨道的耦合，在本例留下只有0.4 eV的孤立位点；它不等于对被删区域做严格的自能消去。更完整的开边界处理见递归自能。