ABACUS 选带部分电荷与网格检查

下载完整输入、原生结果与绘图脚本。解压进入附件根目录，先运行：

bash analyze.sh

这一步读取参考结果，用一个CPU核心提取CSV并生成四张高清图，不启动DFT。完整计算包含6个步骤，提交：

sbatch submit.slurm

CPU模板预留1个MPI进程×56核，符合整节点规则；此微型流程逐个条件实际仅用1线程。参考实算使用ABACUS3.10.1的P100模块、CPU后端，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2及OpenBLAS0.3.30；不是GPU加速或8173M节点性能数据。附件的submit-p100-cpu.slurm对应参考模块和资源环境，申请一张P100资源后使用CPU计算；按自己的资源权限选择，不能把CPU模块和GPU节点模块混用。不要在登录节点执行DFT。

1. 哪些电子态分布在哪里#

总电荷回答所有占据态的整体分布，部分电荷可单独查看某条带、某个k点。例子仍用Si2，只选择第1条带，完整8个k点分别输出；它是理解文件和流程的小体系，不是特定缺陷态或界面电荷转移研究。

2. 先做同一结构的LCAO自洽#

先用局域轨道波函数教程中的完整lcao-scf输入，保留WFC_NAO_K1.txt至WFC_NAO_K8.txt。本例要求相同结构、赝势、轨道、k点、带数，前后都用symmetry -1，避免不可约k点和完整k点混用。

3. 再做get_pchg#

INPUT_PARAMETERS
suffix tutorial
calculation get_pchg
basis_type lcao
ecutwfc 80
scf_thr 1e-8
scf_nmax 150
ks_solver elpa
device cpu
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
out_wfc_lcao 1
nbands 14
symmetry -1
init_wfc file
out_wfc_r 1
if_separate_k 1
out_pchg 1 13*0

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

K_POINTS
0
Gamma
2 2 2 0 0 0

关键参数为：

• calculation get_pchg：读取已保存波函数，提取部分电荷；不重新求一套不同的SCF。
• init_wfc file：从匹配文件初始化。
• out_pchg 1 13*0：14条带的0/1掩码，只选择第1条。它不是“输出第1到第13条带”。
• if_separate_k 1：每个k点独立输出。本例应恰好有8个BAND1_K*_SPIN1_CHG.cube。
• symmetry -1：保留完整k点并避免部分电荷对称性/权重歧义。

在已分配资源并加载模块的环境中，单独运行该阶段：

cd work/cases/lcao-get-pchg
mkdir -p OUT.tutorial
cp ../lcao-scf/OUT.tutorial/WFC_NAO_K*.txt OUT.tutorial/
mpirun -np 1 abacus > partial-charge.log

这里应使用out_pchg。本例实际测试中，仅设置bands_to_print时LCAO仍按默认费米能附近范围输出；所以需要检查最终文件名和数量，不能只看输入里是否写了选带字符串。

4. 分别输出的cube如何读#

这一3.10.1路径在if_separate_k 1时使用单态归一化，不包含k点权重或自旋占据因子，每张cube的积分应接近1，不能直接把八张相加称为总电子密度。使用ASE读取cube时，原子坐标会转成Å，数据仍是原生Bohr⁻³，因此积分需要用单胞体积（ų）除以Bohr³。

module load ase/3.29.0
python - <<'PY'
from pathlib import Path
from ase.io.cube import read_cube
from ase.units import Bohr
root=Path('reference/cases/lcao-get-pchg/OUT.tutorial')
for p in sorted(root.glob('BAND1_K*_SPIN1_CHG.cube')):
    with p.open() as f:
        a=read_cube(f)
    integral=a['data'].mean()*a['atoms'].get_volume()/Bohr**3
    print(p.name, integral)
PY

5. 为什么还要检查网格#

LCAO的ecutwfc控制网格精度，数值轨道本身由轨道文件定义。相同结构分别用40和80 Ry完成自洽，并用各自对应系数提取部分电荷，实际积分如下：

40 Ry的第4点约0.9928，不能仅因文件生成就当作归一化精确。80 Ry本例全部点接近1，最大偏差约6.6×10⁻⁸。这是本例的实际网格对照，不代表所有材料都只需80 Ry；换赝势、轨道或体系后仍要检查。

对照目录是lcao-grid40-scf及lcao-grid40-pchg。每套后处理读取自己的SCF系数，不把40 Ry和80 Ry结果混在一个OUT目录里。为绘图统一色标时，保持数据单位、原胞坐标和所选态相同；单条简并带的图像也不唯一，需要按研究目的合并完整子空间。

原理与实现参考：官方波函数教程、3.10.1 LCAO部分电荷实现。