SOC 自旋投影
VASP · SnSe₂/Sr₂N 的 SOC 沿线自旋投影
方法参考能带
SnSe₂/Sr₂N 的 SOC 能带中，费米能附近各态的投影磁化主要沿哪个方向，符号如何随路径改变？把 PROCAR 的三个分量与同一态的能量配对，才能在色散中辨认这些变化。这里沿 SnSe₂/Sr₂N 已计算的 Γ–M–K–Γ 路径读取 PROCAR，保留它的原始投影数值。
Lu 等对门控 MoS₂ 的研究在 Fig. 1A、1C 把相反能谷的面外自旋取向与 SOC 劈裂联系起来，在 Fig. 3C、3D 再用面内上临界场检验超导态的磁场响应。沿不同方向读取费米能附近的自旋分量，是理解这条路线的电子结构起点。下面先把真实 PROCAR 中的投影量读准，再讨论二维自旋纹理需要增加什么数据。
这是一条高对称线上的自旋投影路线。要画二维 k 平面的箭头图，需要额外计算平面网格；不能把下面的线数据摊成一张二维纹理图。SCF 页可用于对照静态输入与输出的读法，能带方法目录说明相应的数据需求。本页直接从这份已经结束的 SOC 能带输出开始。
VASP：PROCAR · LSORBIT · SAXIS · 自旋纹理
下载原始 PROCAR、对应 SCF 输出和提取脚本。包内保留 150×72 组完整的四块投影数据，可重新生成 spin-path.dat；不包含大体积 CHGCAR 或波函数。路径计算继承同一 SOC SCF 的结构、赝势和密度，并关闭 LCHARG。
核对 SOC 密度、路径和自旋基底
[bcgong@localhost snse2_sr2n_spin]$ head -8 POSCAR
11
1.00000000000000
3.9501156207146009 0.0000000001522404 -0.0000000000000000
-1.9750578097205296 3.4209004743098745 0.0000000000000000
-0.0000000000000001 0.0000000000000002 39.4021877938467284
Sn Se N Sr
1 2 1 2
Direct
元素和原子数是 Sn Se N Sr / 1 2 1 2，合计六个原子。部分旧输入中的 SYSTEM 名称没有随着结构更新，判断材料时应以 POSCAR 的元素与坐标为准。
[bcgong@localhost snse2_sr2n_spin]$ grep -E '^[[:space:]]*(ICHARG|LNONCOLLINEAR|LSORBIT|LORBIT|LMAXMIX|ENCUT|EDIFF|ISMEAR|SIGMA)[[:space:]]*=' INCAR
ICHARG = 11 charge: 1-file 2-atom 10+ const
LORBIT = 11 DOSCAR and lm decomposed PROCAR file
LNONCOLLINEAR = T whether perform non-colli-mag calculation
LSORBIT = T whether perform soc calculation
ENCUT = 520 eV
EDIFF = 1E-6 :break condition for electron scf loop
LMAXMIX = 4 d-elements-4 f-elements-6
ISMEAR = 0 :-5-tet -1-fermi 0-gauss
SIGMA = 0.05 :boadening in eV
这份能带计算使用 ICHARG = 11，从相应 SOC SCF 的固定电荷密度求沿线本征态；LORBIT = 11 让 PROCAR 写出原子与轨道投影。LSORBIT 需要非共线可执行程序。读取这种 PROCAR 时，一条 band 下不仅有总权重，还会有三个磁化分量。
固定密度后仍要解这些路径点上的波函数，所以 ICHARG=11 并不跳过本征态求解。这里的 LMAXMIX=4 还涉及 CHGCAR 中保存和读取的 PAW 单中心电荷角动量分量；应与生成密度的 SCF 一起核对。只在路径输入里调大它，不能恢复父 CHGCAR 没有保存的分量。
[bcgong@localhost snse2_sr2n_spin]$ cat KPOINTS
K-Path Generated by VASPKIT.
50
Line-Mode
Reciprocal
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.0000000000 0.0000000000 M
0.5000000000 0.0000000000 0.0000000000 M
0.3333333333 0.3333333333 0.0000000000 K
0.3333333333 0.3333333333 0.0000000000 K
0.0000000000 0.0000000000 0.0000000000 GAMMA
每段 50 个点，三段共 150 个记录，端点在相邻段中重复。这与普通二维均匀 k 网格不同，也没有覆盖整个布里渊区。
[bcgong@localhost snse2_sr2n_spin]$ grep -A4 'transformation matrix from SAXIS to cartesian' OUTCAR
transformation matrix from SAXIS to cartesian coordinates
---------------------------------------------------------
1.0000000 m_x 0.0000000 m_y 0.0000000 m_z
0.0000000 m_x 1.0000000 m_y 0.0000000 m_z
0.0000000 m_x 0.0000000 m_y 1.0000000 m_z
日志记录 VASP 5.4.4 的 complex 版本；父 SOC SCF 与路径计算都采用非共线自旋。父 SCF 的 NELECT=51，最后总磁化约为 (0.0015435, 0.0161994, 0.0000018)；路径输出的末值见下文。这里按实际运行记录解释投影，不能预设这份密度严格满足时间反演对称性。
在当前晶胞、粒子数守恒的自旋子带描述下，51 个电子不能填满时间反演对称绝缘体的 Kramers 成对占据带；奇数填充允许时间反演对称金属，不能单凭它断言时间反演破缺。若讨论分离的低能带子空间，应另说明子空间与实际费米占据的关系。Soluyanov 与 Vanderbilt 的 Sec. II.1、式 (3)给出成对占据态的时间反演关系。
本例自旋基底到 Cartesian 坐标的变换是单位矩阵，三个磁化分量可以依次记作 mx、my、mz。下方 spin_path.py 按本例的默认 SAXIS=(0,0,1) 直接保存这三个分量；使用其他 SAXIS 时，需在读取投影后加入对应的笛卡尔坐标变换。
逐态读取四组投影
[bcgong@localhost snse2_sr2n_spin]$ head -38 PROCAR
PROCAR lm decomposed
# of k-points: 150 # of bands: 72 # of ions: 6
k-point 1 : 0.00000000 0.00000000 0.00000000 weight = 0.00666667
band 1 # energy -36.41469553 # occ. 1.00000000
ion s py pz px dxy dyz dz2 dxz x2-y2 tot
1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
4 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
5 0.089 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.089
6 0.885 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.885
tot 0.974 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.974
1 0.000 -0.000 0.000 -0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000
2 0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000 -0.000 -0.000 0.000
3 0.000 0.000 0.000 -0.000 -0.000 0.000 0.000 -0.000 0.000 0.000
4 0.000 0.000 0.000 -0.000 0.000 0.000 0.000 0.000 0.000 0.000
5 0.016 0.000 0.000 -0.000 -0.000 0.000 0.000 -0.000 0.000 0.016
6 0.159 0.000 0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 0.159
tot 0.176 0.000 0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 0.176
1 0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000 0.000 0.000 0.000
2 0.000 -0.000 0.000 -0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000
3 0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 0.000
4 0.000 -0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
5 0.087 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 0.087
6 0.868 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 0.868
tot 0.956 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 0.956
1 0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 0.000
2 0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000 0.000 0.000 0.000
3 0.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 0.000 -0.000 0.000
4 0.000 -0.000 0.000 -0.000 0.000 0.000 0.000 0.000 0.000 0.000
5 0.007 -0.000 0.000 -0.000 -0.000 0.000 0.000 0.000 0.000 0.007
6 0.064 -0.000 0.000 -0.000 0.000 -0.000 0.000 -0.000 -0.000 0.064
tot 0.071 -0.000 0.000 -0.000 0.000 -0.000 0.000 -0.000 -0.000 0.071
band 2 # energy -36.41344142 # occ. 1.00000000
文件第二行给出 150 个 k 点、72 条带和 6 个原子。一个 k-point 行后面跟 band 编号、能量和占据数；随后有四组投影表。第一组的 tot 是投影权重，后三组的 tot 分别是三个方向的投影磁化。上面的第一个态给出权重 0.974、mx = 0.176、my = 0.956、mz = 0.071。
第一组没有恰好等于 1，是因为这里读的是原子投影子空间中的权重。本例保留后三组磁化原值，颜色表示原子投影子空间中的磁化。完整布洛赫态的归一化自旋期望值需要覆盖全态的自旋矩阵元。
这四组 tot 不能合并成一个“自旋大小”。第一组告诉我们这个态落在所选原子投影空间中的权重，后三组才带方向和正负号；tot 也已经对六个原子求和，层来源要回到逐原子行另作分组。原始权重的变化可以使两处相同取向的态呈现不同颜色深浅。若另存 mz/charge 来比较投影空间内的取向，应保留分子、分母并说明这个比值的定义，不用它替换当前图的原值，更不能把它直接标成全态的 ⟨Sz⟩。
[bcgong@localhost snse2_sr2n_spin]$ tail -3 OSZICAR
DAV: 9 -0.258918997342E+02 -0.18555E-05 -0.18547E-05 26488 0.101E-02
DAV: 10 -0.258918998879E+02 -0.15365E-06 -0.15301E-06 20504 0.363E-03
1 F= -.26828772E+02 E0= -.26827878E+02 d E =-.178699E-02 mag= 0.0017 0.0167 0.0001
[bcgong@localhost snse2_sr2n_spin]$ grep 'aborting loop because EDIFF is reached' OUTCAR
------------------------ aborting loop because EDIFF is reached ----------------------------------------
沿线计算第 10 个电子步满足停止条件。最后仍有约 (0.0017, 0.0167, 0.0001) 的小总磁矩，这份结果并非严格保持零磁矩的验证计算。图中的分裂与颜色可以用于学习如何提取输出；若要把它们解释为某种材料自旋劈裂，还应检查初始磁矩、对称性、SCF 收敛和这种残余磁化的影响。
[bcgong@localhost snse2_sr2n_spin]$ grep 'E-fermi' SCF_OUTCAR OUTCAR
SCF_OUTCAR: E-fermi : -1.4881 XC(G=0): -3.0426 alpha+bet : -2.9160
OUTCAR: E-fermi : -1.5887 XC(G=0): -3.0424 alpha+bet : -2.9160
作图采用生成固定电荷密度的 SCF 费米能 −1.4881 eV。沿线 k 点不是积分网格，不能用沿线输出里重新给出的费米能取代 SCF 的参考。
随例子提供的脚本逐个读取 k-point、band 和四个 tot 行，检查 150 × 72 个块全部齐全。k 的横轴根据 POSCAR 的倒易晶格换算为路径长度，并使用 SCF 的费米能移动能量零点。
[bcgong@localhost snse2_sr2n_spin]$ head -4 spin-path.dat
# ik band kdist_A-1 k1 k2 k3 energy_eV E-Ef_eV occupation charge mx my mz
1 1 0.000000000 0.000000000 0.000000000 0.000000000 -36.414695530 -34.926595530 1.000000000 0.974000000 0.176000000 0.956000000 0.071000000
1 2 0.000000000 0.000000000 0.000000000 0.000000000 -36.413441420 -34.925341420 1.000000000 0.974000000 -0.176000000 -0.956000000 -0.071000000
1 3 0.000000000 0.000000000 0.000000000 0.000000000 -36.277903420 -34.789803420 1.000000000 0.978000000 0.270000000 0.921000000 0.185000000
输出的十三列依次是 k 点编号、带号、累计路径长度、三个分数倒易坐标、原能量、相对 SCF 费米能的能量、占据、投影权重、mx、my、mz。它没有把相邻带号自动当作同一条连续自旋分支；带交叉和近简并处仍需结合波函数连续性判断。
自旋投影的读取与作图
下面的任务从 PROCAR 读取逐态自旋投影，按同一路径和能量参考画图，同时保留自旋基底与投影归一化的说明。
编写 SnSe₂/Sr₂N SOC 路径投影的两步后处理程序，使用 Python 3、NumPy 和 Matplotlib。
第一步提取：输入 PROCAR、POSCAR、SCF_OUTCAR。PROCAR 有 150 点×72 带，每态依次读取 charge、mx、my、mz 四块 tot；按本例默认 SAXIS=(0,0,1) 直接保存磁化分量。由 POSCAR 倒格矢计算累计路径距离（Å⁻¹），从 SCF_OUTCAR 读取 EF=−1.4881 eV，生成十三列 spin-path.dat 与 spin-summary.json。检查每态四块 tot、末态完整和 10800 条记录；保留残余总磁矩 (0.0017,0.0167,0.0001) 的计算条件。
第二步绘图：输入 spin-path.dat、spin-summary.json，并使用同目录 atlas_plot_style.py；三个磁化面板共享 −1…1 色标及 EF±2 eV 窗口，横轴和节点取提取结果。
输出：两步完整源码、依赖、命令、表格与 PNG/PDF。颜色表示未经归一化的原子投影磁化；当前源码按默认 SAXIS 工作，二维纹理需另取平面网格。
后处理源码与运行
完整源码：spin_path.py · plot_spin.py · atlas_plot_style.py。Python 3 依赖：NumPy、Matplotlib。
spin_path.py 的完整源码
from __future__ import print_function
import re, math, json
def cross(a,b): return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]
def dot(a,b): return sum(x*y for x,y in zip(a,b))
p=open('POSCAR').readlines(); scale=float(p[1]); cell=[[float(x)*scale for x in t.split()] for t in p[2:5]]
vol=dot(cell[0],cross(cell[1],cell[2])); rec=[[2*math.pi*x/vol for x in cross(cell[(i+1)%3],cell[(i+2)%3])] for i in range(3)]
ef=float(re.findall(r'E-fermi\s*:\s*([-0-9.]+)',open('SCF_OUTCAR').read())[-1])
rows=[]; dist=0.; prev=None; k=None; b=None; sums=[]
for line in open('PROCAR'):
m=re.match(r'\s*k-point\s+(\d+)\s*:\s*([-0-9.]+)\s+([-0-9.]+)\s+([-0-9.]+)',line)
if m:
k=list(map(float,m.group(2,3,4))); ik=int(m.group(1))
cart=[sum(k[i]*rec[i][j] for i in range(3)) for j in range(3)]
if prev is not None: dist+=math.sqrt(sum((x-y)**2 for x,y in zip(cart,prev)))
prev=cart
m=re.match(r'\s*band\s+(\d+)\s+# energy\s+([-0-9.]+)\s+# occ.\s+([-0-9.]+)',line)
if m:
if b is not None and len(sums)!=4: raise ValueError('Need exactly four tot rows per band')
b=int(m.group(1)); energy=float(m.group(2)); occ=float(m.group(3)); sums=[]
if line.strip().startswith('tot '):
sums.append(float(line.split()[-1]))
if len(sums)==4:
rows.append([ik,b,dist]+k+[energy,energy-ef,occ]+sums)
if len(sums)!=4: raise ValueError('Truncated last band')
if len(rows)!=150*72: raise ValueError('Unexpected k/band block count')
with open('spin-path.dat','w') as f:
f.write('# ik band kdist_A-1 k1 k2 k3 energy_eV E-Ef_eV occupation charge mx my mz\n')
for a in rows: f.write('%d %d '%(a[0],a[1])+' '.join('%.9f'%v for v in a[2:])+'\n')
summary={'nk':150,'nb':72,'nrows':len(rows),'fermi_scf_eV':ef,'axis':'default SAXIS=(0,0,1), Cartesian','quantity':'native PROCAR projected magnetization; no normalization','path_ticks_A-1':[rows[i*72][2] for i in [0,49,99,149]],'max_abs_mz':max(abs(a[-1]) for a in rows)}
json.dump(summary,open('spin-summary.json','w'),indent=2)
print('Parsed %d k-points x %d bands = %d four-block records'%(150,72,len(rows)))
print('SCF E-fermi = %.4f eV'%ef)
print('path ticks / A^-1 = '+str(summary['path_ticks_A-1']))
print('columns charge,mx,my,mz; output spin-path.dat')
plot_spin.py 的完整源码
from atlas_plot_style import install as install_atlas_style
install_atlas_style()
import json
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
x=np.loadtxt('spin-path.dat'); s=json.load(open('spin-summary.json'))
a=x[(x[:,7]>=-2)&(x[:,7]<=2)]
fig,axes=plt.subplots(1,3,figsize=(12,4.5),sharex=True,sharey=True,layout='constrained')
for ax,col,label in zip(axes,[10,11,12],['m_x','m_y','m_z']):
sc=ax.scatter(a[:,2],a[:,7],c=a[:,col],s=5,cmap='coolwarm',norm=Normalize(-1,1),rasterized=True)
ax.axhline(0,color='0.35',ls='--',lw=.7)
for p in s['path_ticks_A-1']: ax.axvline(p,color='0.8',lw=.6,zorder=0)
ax.set_xticks(s['path_ticks_A-1'],['Γ','M','K','Γ'])
ax.set_title(label+' (native PROCAR projection)')
ax.set_xlabel('High-symmetry path')
axes[0].set_ylabel('Energy relative to SCF Fermi level (eV)')
fig.colorbar(sc,ax=axes,label='Projected magnetization',shrink=.75)
fig.suptitle('SnSe2/Sr2N: spin projections along Gamma-M-K-Gamma')
fig.savefig('spin-path.png',dpi=220)
fig.savefig('spin-path.pdf')
解压本页示例包后，在 snse2-sr2n-spin-path 根目录执行：
python3 -m pip install numpy matplotlib
python3 spin_path.py
python3 plot_spin.py
本例保存的提取运行记录如下：
[bcgong@localhost snse2_sr2n_spin]$ python spin_path.py
Parsed 150 k-points x 72 bands = 10800 four-block records
SCF E-fermi = -1.4881 eV
path ticks / A^-1 = [0.0, 0.9183525437502472, 1.4485636269770399, 2.5089857932241535]
columns charge,mx,my,mz; output spin-path.dat
spin_path.py 从原始 PROCAR 检查每态四个 tot 块并生成十三列表；已有 spin-path.dat 与 spin-summary.json 时，直接执行 python3 plot_spin.py。绘图入口使用已提取的表格。
将 spin-path.dat、spin-summary.json、plot_spin.py 和同目录的 atlas_plot_style.py 一起放到本机，执行 python3 plot_spin.py。脚本并排画出 mx、my、mz 三幅着色能带，显示费米能上下 2 eV，三幅图共用 −1 到 1 的颜色标尺，同时输出 PNG 与 PDF。颜色在近简并态之间跳变时，先检查成对态和投影基底，不要把每个带号的突变都解释成独立的物理纹理。
下一步若需更密的 k 平面数据，可到 Wannier90 方法目录查看插值所需的波函数与接口数据，再验证插值能带和自旋矩阵元。拓扑量的数据需求另见 Berry 曲率与 Chern 数方法目录；目录中已有例程使用各自的材料和程序，不能仅凭当前 PROCAR 接续得到这些量。
从沿线投影读到费米面自旋取向
Lu 等的 Fig. 1A、1C（原文 PDF 第 1 页）可以与上图并排读：1A 在六角布里渊区标出 K/K′ 谷及方向相反的面外有效场，蓝、红口袋分别示意自旋向上、向下；1C 画相应能带劈裂，并显示 2H 堆垛相邻层在同一 K 谷的有效场反向。颜色是方向示意，没有连续自旋数值色标。本页则在 Γ–M–K–Γ 路径上用同一个 −1…1 色标分别画 PROCAR 的 mx、my、mz，纵轴统一减去 SCF EF；投影权重没有被归一化掉。对照时应比较谷位置、分量方向与层来源，不能把颜色深浅直接换算成有效场，或将沿线投影当成二维费米面纹理。
颜色是 PROCAR 原子投影空间中的磁化，横轴是 Γ–M–K–Γ。具体看 K 点（第 100 个记录）附近的一对输出态：
| 带号 | E−EF / eV | 投影权重 | mx | my | mz |
|---|---|---|---|---|---|
| 51 | +0.120867010 | 0.644 | −0.021 | −0.016 | −0.641 |
| 52 | +0.124833500 | 0.643 | +0.021 | +0.016 | +0.641 |
两条记录相差 3.96649 meV，投影磁化主要沿 z，符号相反。它们位于 SCF 费米能上方约 0.12 eV，PROCAR 占据数均为零；这两行不能代替费米面上的配对电子。Γ 点的近简并态也可出现相反投影，简并子空间内的基底选择会改变单态颜色。沿线路径、投影权重和父密度的残余磁化共同限定了这张图的读法。
“主要沿 z”还可以用现有列核对。将面内投影幅值与 |mz| 相比，这两态都得到 0.041187，说明在当前投影空间里面内分量约为面外幅值的 4.1%；mz/charge 分别为 −0.995342、+0.996890。这些比值从 PROCAR 已保存的三位小数投影算出，显示六位是为了重现列运算。比值不改变两态都未占据的事实，也没有补回原子投影外的自旋。两个相反颜色来自同一个 K 点的两个态；时间反演比较则要找 −K 所在的 K′ 谷，并核对相应能量、子空间和三分量。本页路径没有 K′，不能把这两行当作已验证的时间反演伙伴。
下面的读取只核对这对记录，不重新拟合能带或绘图：先按 k 点记录号和带号找行，检查两行路径坐标相同，再由原始列计算能量差、mz/charge 和面内/面外幅值比。可交给编程助手的要求是：读取已经生成的十三列 spin-path.dat，默认选第 100 个记录的 51、52 带；文件列数、有限值、记录完整性和坐标不符就停止，原始文件只读。输出必须同时显示 E−EF、占据和投影权重，并标明两个比值不是全态自旋或时间反演检验。
完整读取脚本 inspect_spin_pair.py只使用 Python 标准库。把它放在本页解压目录的 spin-path.dat 旁，运行：
python3 inspect_spin_pair.py spin-path.dat --ik 100 --bands 51 52
本轮从保存的真实 spin-path.dat 读取输出为：
ik band E-EF_eV occupation charge mz/charge inplane/abs(mz)
100 51 +0.120867010 0.000 0.644 -0.995342 0.041187
100 52 +0.124833500 0.000 0.643 +0.996890 0.041187
energy_difference_meV=3.966490
Ratios describe the saved atomic projections; no full-state spin or TR test.
inspect_spin_pair.py 完整源码
#!/usr/bin/env python3
"""Inspect saved projection columns at one path record, without fitting or plotting."""
import argparse
import math
from pathlib import Path
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('data', type=Path)
parser.add_argument('--ik', type=int, default=100)
parser.add_argument('--bands', type=int, nargs=2, default=[51, 52])
args = parser.parse_args()
if args.ik < 1 or min(args.bands) < 1 or args.bands[0] == args.bands[1]:
parser.error("Choose a positive path record and two distinct positive bands")
rows = {}
for line in args.data.read_text().splitlines():
if not line.strip() or line.lstrip().startswith('#'):
continue
values = [float(x) for x in line.split()]
if len(values) != 13 or not all(math.isfinite(x) for x in values):
raise ValueError('Expected thirteen finite spin-path.dat columns')
ik, band = values[:2]
if ik != int(ik) or band != int(band):
raise ValueError('Noninteger record or band index')
if int(ik) == args.ik and int(band) in args.bands:
if int(band) in rows:
raise ValueError('Duplicate requested record')
rows[int(band)] = values
if set(rows) != set(args.bands):
raise ValueError('Requested pair is incomplete')
a, b = [rows[n] for n in args.bands]
if a[2:6] != b[2:6]:
raise ValueError('Pair does not share the same path coordinate')
print('ik band E-EF_eV occupation charge mz/charge inplane/abs(mz)')
for row in (a, b):
q, mx, my, mz = row[9:13]
if q <= 0 or mz == 0:
raise ValueError('Ratio requires positive weight and nonzero mz')
print(f'{int(row[0])} {int(row[1])} {row[7]:+.9f} {row[8]:.3f} '
f'{q:.3f} {mz/q:+.6f} {math.hypot(mx,my)/abs(mz):.6f}')
print(f'energy_difference_meV={(b[6]-a[6])*1000:.6f}')
print('Ratios describe the saved atomic projections; no full-state spin or TR test.')
用 Wannier 插值扩展网格时，还要验证目标能区的色散和自旋矩阵元；仅有能量一致的 hr.dat 不能保证任意轨道基底上的 Pauli 矩阵就是原始 DFT 自旋算符。接口和子空间的准备见 Wannier90。
同文 Fig. 3A、3B（PDF 第 3 页）从电阻曲线与 RN/2 的交点提取面内 Bc₂；3C 比较临界场随温度的变化及不同拟合，3D 再用横轴 T/Tc、纵轴 Bc₂/Bp 比较不同超导态，虚线 Bc₂/Bp=1 是 Pauli 界限。图 1D 的相图使用 90% RN 的 Tc 判据，不能与图 3 的 50% RN 混用。复现这种归一化比较需要原始磁输运曲线、各态一致的 Tc 判据和相应 Bp，逐态提取交点后再画无量纲坐标；本页的 PROCAR 与正常态能带没有这些超导态输入。
正常态 Z₂ 与边界态从占据波函数或经验证的 SOC 哈密顿量出发，不能从这份逐态 PROCAR 的颜色推出，见 Wilson loop 与拓扑判读。
SOC SCF 密度与 EF → 二维 k 网格/费米线 → 三分量自旋与层投影
└─ −k 伙伴、能谷、面内/面外分量的比较
SOC 波函数 → 经验证的自旋子 Wannier 模型 → WCC / Z₂ / 边界谱
EPC 与配对模型 → 超导态及磁场响应