ABACUS 密度矩阵、电子数与布居实操

密度矩阵不能只看对角线再求和：非正交局域基还需要重叠矩阵。本篇接着H/S能带重建，从同一份真实Si2计算的数据得到电子数、轨道与原子Mulliken布居以及带能，再和原生能级占据、电荷密度网格核对。

1. 运行与查看#

下载矩阵、原生数据和脚本，解压进入根目录：

bash run.sh

一个CPU核心即可；脚本读取已有DFT结果，不运行新的SCF。原生ABACUS3.10.1的完整输入与编译环境见上一篇，分析环境为ase/3.29.0。

head -10 work/orbital-populations.csv

中图画的是Σ D(R)S(−R)，不是单独某一块D(R)。右图前13个轨道属于第一个Si原子，后13个属于第二个；颜色仅区分原子。NAO基不必满足S的每个对角元都是1，也不能按“每个轨道正好一个电子”解释布居。

2. 输出和轨道顺序#

在LCAO多k点SCF中，本例设置：

out_dm1 1
out_mat_hs2 1
out_chg 1

原生data-DMR-sparse_SPIN0.csr含183个R块，H/S各含177个。R块数量不同不一定是数据丢失；组合时按实际R坐标匹配，不能把第i块D和第i块S机械相乘。文件的轨道顺序和基函数必须相同，本例每个Si为2s、2p、1d组径向轨道，共2+2×3+1×5=13个。

本篇只针对实际验证过的nspin=1实数CSR。自旋极化、SOC自旋量、复数格式和新版接口不能套用一个额外乘2就直接处理。当前非自旋极化D及原生占据已经包含相应占据因子，不再重复乘自旋简并度。

3. 实空间电子数怎么得到#

本例按照相反的R相乘：

N = Σ_R Tr[D(R) S(−R)]

对比结果如下：

两个Si赝势各提供四个价电子，因此这里核对的是8个价电子，不是包含芯电子的28个全电子。istate.info本例的占据已经包含k权重；不能再乘一次k权重或自旋因子。cube差别来自文本输出和积分精度，不能把最后几位当成电荷转移。

矩阵电子数与cube积分的区别，也可结合电荷、ELF与势的网格分析理解。cube的电荷密度单位为Bohr⁻³，ASE读取的晶胞是Å，积分必须保持体积单位一致；附件用ASE读取并换回Bohr³，避免直接相乘导致单位错误。

4. 原子Mulliken布居#

将Σ D(R)S(−R)的对角元素按原子所属轨道求和，得到两个Si的布居约为：

原子布居总和等于价电子总数。两个对称Si之间约10⁻⁴电子量级的数值差异，不应被解释成可靠的电荷转移。Mulliken依赖局域轨道基，不等于形式氧化态，也不等于Bader划分；要比较不同材料或基组，应保持方法和参考一致。

5. 带能不等于DFT总能量#

同样可以计算：

E_band = Σ_R Tr[D(R) H(−R)]

本例矩阵结果为15.4081919722 eV，原生E_band为15.4081920055 eV；按istate.info的本征值乘加权占据求和也得到后者。DFT总能量还包含Hartree、交换关联等修正，不能拿这项带能直接计算形成能、吸附能或相稳定性。

在附件根目录提取结果：

python - <<'PY'
import json
r = json.load(open('work/RESULTS.json'))
for key in ['matrix_electrons', 'occupation_electrons', 'cube_electrons',
            'atom_Mulliken_electrons', 'band_trace_eV', 'native_band_eV',
            'total_energy_eV']:
    print(key, r[key])
PY

完整结果在work/RESULTS.json，轨道布居在work/orbital-populations.csv。先用这种电子数和能级交叉检查确认接口正确，再把矩阵送入自己的后处理程序；这一步不会自动完成DeePH、DeePTB、TB2J或其他模型训练。