pymatgen：实际波函数的SOC spillage

SOC打开前后，波函数的占据子空间可能发生变化。pymatgen的SOCSpillage可以读取两套相匹配的WAVECAR，计算其原生定义下的子空间spillage。

本例在鸢算科技集群实际完成同一个Bi₂结构的两套VASP6.6.1单点：标量双自旋计算用P100，SOC计算用V10016G；均为单MPI、单科学线程，NVHPC25.3 GPU构建，EDIFF=1E-6。没有进行新的结构优化。

下载完整输入、提交/分析脚本、实际数值与绘图资料。POTCAR和WAVECAR不随包发布；可先读取已有数值，再在已获VASP授权的环境中按输入重算。

unzip pymatgen-SOC-spillage-wave720-*.zip
cd pymatgen-SOC-spillage
module load pymatgen/2026.9.24
bash results.sh
bash run.sh
bash plot.sh

run.sh只从实际主夹角数据复算各k点spillage，写入analysis-reference，不启动DFT。plot.sh从参考CSV生成figures。

1. 两套计算必须真正匹配#

本例使用相同POSCAR、Bi的PBE64赝势、ENCUT=300eV、24带和四个显式k点；仅改变是否包含SOC与相应波函数表示。

inputs/scalar/POSCAR与SOC目录中的文件相同：

Bi buckled honeycomb layer             
   1.00000000000000     
     4.3888421599896201    0.0000000000000000    0.0000000000000000
     2.1944210799948101    3.8008488037511787    0.0000000000000000
     0.0000000000000000    0.0000000000000000   15.0000000000000000
   Bi
     2
Selective dynamics
Direct
  0.1666666666666643  0.1666666666666643  0.0579618666614847   F   F   T
  0.8333333333333357  0.8333333333333357  0.9420381333385153   F   F   T
 
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

结构为固定Bi₂层状模型，真空方向15Å。是否用于真实材料研究，还需另外检验几何、真空、截断能和k网格；本例只展示四点后处理流程。

共同KPOINTS：

Matched explicit shifted mesh
4
Reciprocal
0.125 0.125 0 0.25
0.625 0.125 0 0.25
0.125 0.625 0 0.25
0.625 0.625 0 0.25

显式列出并关闭对称性，是为了两套WAVECAR保留相同k点及顺序。不能把两个程序产生的不同不可约k点列表按行号直接比较。

2. 两份完整INCAR#

标量双自旋：inputs/scalar/INCAR#

SYSTEM = Bi2 matched occupied-subspace comparison
ENCUT = 300
PREC = Accurate
EDIFF = 1E-6
NELM = 120
ALGO = Normal
ISYM = -1
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.02
IBRION = -1
NSW = 0
NBANDS = 24
LASPH = .TRUE.
ADDGRID = .TRUE.
LWAVE = .TRUE.
LCHARG = .FALSE.
NCORE = 1
KPAR = 1
ISPIN = 2
MAGMOM = 2*0

SOC旋量：inputs/SOC/INCAR#

SYSTEM = Bi2 matched occupied-subspace comparison
ENCUT = 300
PREC = Accurate
EDIFF = 1E-6
NELM = 120
ALGO = Normal
ISYM = -1
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.02
IBRION = -1
NSW = 0
NBANDS = 24
LASPH = .TRUE.
ADDGRID = .TRUE.
LWAVE = .TRUE.
LCHARG = .FALSE.
NCORE = 1
KPAR = 1
LSORBIT = .TRUE.
LNONCOLLINEAR = .TRUE.
SAXIS = 0 0 1
MAGMOM = 6*0

标量计算明确设ISPIN=2，保留自旋上下两个通道；本例初始磁矩为0，最终每个k点上下自旋各有5条明确占据带。SOC旋量则选10条占据带，对应同样10个价电子。

MAGMOM长度和波函数表示要对应：标量每原子一个数，SOC非共线每原子三个分量。这里的单点不移动原子，POSCAR选择性约束不会变成实际几何优化。

3. 准备赝势并用Slurm提交#

bash prepare-potcar.sh
bash submit.sh
squeue -u "$USER"

脚本分别进入inputs/scalar和inputs/SOC提交job.slurm。模板对应P100-HS的vasp/6.6.1-p100、V10016G-HS的vasp/6.6.1-v100，请按自己账户可用资源选择适用模块、分区和GPU类型。

提交脚本使用一个MPI进程和mpirun --bind-to none，仍由Slurm分配CPU/GPU；计算输出包括各自的OUTCAR、vasprun.xml和WAVECAR。

两套都完成后，从包根目录执行：

bash analyze-new.sh
cat analysis-native/Bi2-SOC-spillage-kpoints.csv

实际前置计算分别在19步和22步电子迭代达到EDIFF，Slurm用时约39秒和80秒；不同负载与体系下耗时会变化。

4. 原生接口如何调用#

python - <<'PY'
from pymatgen.analysis.topological.spillage import SOCSpillage
comparison = SOCSpillage(
    wf_noso='inputs/scalar/WAVECAR',
    wf_so='inputs/SOC/WAVECAR')
print('采样点最大spillage:', comparison.overlap_so_spinpol())
PY

该版本原生接口主要返回采样点中的最大值；完整脚本另外逐k点求独立主夹角结果，并核对原生最大值。

本例严格检查晶格、ENCUT、k点和逐点G向量及其顺序一致，没有用截短或取最小长度的方式拼凑不同基底。不同k点各自可以有不同平面波数量，但同一个k点的两套计算必须对应。

5. 读实际数值#

cat reference-data/Bi2-SOC-spillage-kpoints.csv
cat reference-data/Bi2-sampled-direct-gaps.csv

原生最大值为0.3512587923324162，独立QR正交化和主夹角算法得到0.3512587923324011，差异约1.5×10⁻¹⁴。

四点的标量直接带隙约1.64–3.22eV、SOC直接带隙约1.02–2.80eV，仅用于检查这里选择的占据子空间清楚分离；不能当作充分采样后的材料基本带隙。

6. 主夹角数据与含义#

把两套占据波函数分别正交化，取重叠矩阵的奇异值sᵢ，本例逐点计算：

γ(k)=Σᵢ(1−sᵢ²)。

head -n 12 reference-data/Bi2-SOC-spillage-principal-values.csv

N=10是本例占据子空间维数，不是把所有24条带都放进去。图中每个小格是一个1−sᵢ²贡献；同一k点求和对应左侧柱状图。

原生实现比较的是正交化后的赝平面波系数子空间，这里没有另加入PAW增广重叠修正。spillage反映所选子空间变化，不是Z₂不变量，也不能仅凭一个大于或小于某值的结果判定拓扑相。

7. 换自己的体系#

先保持两套结构、赝势、截断能和k网格一致，分别收敛，确认占据数和子空间选择适用。随包分析脚本按本例10电子、5+5/10占据维数检查；换材料时应修改相应维数与占据判断，不能强行取前10条带。

金属、分数占据或不同k点占据维数变化的情形，需要另外定义适用的投影空间与方法。进一步做拓扑判断时，应结合适用的对称性、Berry/Wilson loop等分析，并进行必要的网格和参数收敛。

接口见pymatgen SOC spillage。