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. 先区分系数和网格#

LCAO写作ψ=Σμ Cμ φμ。WFC_NAO_K*.txt给出原子轨道基上的复数系数，不是电子密度图。原子轨道彼此不正交，归一化必须检查C†SC=1；不能只算Σ|Cμ|²。

本例Si2采用每个Si的2s2p1d数值轨道：2个s径向函数、2组p和1组d，共2+6+5=13个轨道；双原子共26个。14条带在8个完整k点上输出，同一SCF还导出H(R)/S(R)。

2. 第一步：自洽并保存匹配文件#

INPUT_PARAMETERS
suffix tutorial
calculation scf
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

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

这里使用out_wfc_lcao 1输出文本系数，out_mat_hs2 1输出实空间稀疏矩阵。symmetry -1保留完整8个k点，便于之后的get_wf和get_pchg读取同一套波函数。若前一步只保留3个不可约k点，后一步要求8个完整k点，不能只复制现有3个文件然后继续。

3. 第二步：从匹配SCF提取波函数#

附件run.sh会依次完成以下操作。在已经分配到计算资源、且已加载相应ABACUS模块的环境中，单独运行该阶段的核心命令为：

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

后处理INPUT与自洽阶段保持同一STRU、KPT、轨道及14条带，主要设置为：

calculation get_wf
init_wfc file
symmetry -1
out_wfc_r 1
out_wfc_norm 1 13*0
out_wfc_re_im 1 13*0

当前3.10.1 LCAO路径的选带使用out_wfc_norm等对应参数；不要把PW的bands_to_print直接用于这里并假定已选中。本例只用Γ点第1个非简并态说明实空间振幅，文件为BAND1_k_1_s_1_ENV.cube、..._REAL.cube、..._IMAG.cube。

ENV是振幅而非密度；图中REAL/IMAG组合的积分∫|ψ|²dV约为0.999999827。与独立get_pchg的Γ点密度逐点最大差约2.205e-05 Bohr⁻³，来自这一路径的网格投影差异，不能当作两个输出逐点完全相同。

本教程没有把非Γ点get_wf的REAL/IMAG输出解释成已核对的完整复Bloch波函数；查看非Γ态的空间概率分布，请使用下一篇的get_pchg。Γ点整体正负号或复数整体相位不影响概率。

4. H/S与系数必须在同一个坐标和基底里#

WFC_NAO_K*.txt头部k点采用以2π/lat0为单位的笛卡尔坐标；H(R)/S(R)的傅里叶变换需要倒格矢分数坐标。本例原胞矩阵A的各行是无量纲晶格矢量，所以k_fractional = k_cartesian @ A.T。例如(-0.5,0.5,0.5)转为(0.5,0,0)。直接拿头部三数当分数坐标，会得到明显错误的S归一化。

两套网格数据中，每套8个k点×14带的最大C†SC−I误差约4.971e-09，广义本征方程HC=SCε残差约7.095e-09 Ry。Γ点第1带的普通C†C约0.455534，不应为了凑成1而自行缩放系数。

python - <<'PY'
import csv
with open('reference/analysis/lcao-metric.csv') as f:
    for r in csv.DictReader(f):
        if r['source']=='lcao-scf':
            print(r['kpoint'],r['metric_error'],r['residual_Ry'])
PY

系数中能量以Ry给出，本例采用该版本常数1 Ry=13.605698 eV；矩阵的轨道顺序、赝势、轨道文件、k点和结构必须与该次SCF一致。跨体系、跨基组拼接系数没有物理意义。

实现参考：LCAO系数读取、LCAO波函数提取。矩阵本身的处理见H/S矩阵重建能带。