ABACUS Onsager提取与低频外推实操

下载18个条件的输入、原生谱线、电流关联与速度矩阵，解压进入附件根目录。先读取已有结果：

bash analyze.sh

表和数组在reference/analysis，高清图在reference/figures。整套新计算可提交：

sbatch submit.slurm

附件在新的work中运行，已有work时请新解压。当前复现模板分配1张P100，但本功能实际使用该主机的CPU执行路径：它不会因申请了GPU而获得GPU加速。此模板用于复现当前案例；单独的生产电导任务应选择合适CPU资源及匹配构建，勿把device改成gpu就视为已验证。本例没有8173M节点性能数据。

使用ABACUS3.10.1模块abacus/3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2环境，计算指定CPU双精度、1个MPI进程和1个线程；ASE3.29.0/NumPy负责后处理。P100模块的CUDA环境为12.8，这批KG路径未用GPU。

1 先辨认原生输出的DC算法#

当前3.10.1先得到正频率中点谱线，再用前两个点线性外推ω=0：σDC≈1.5σ(0.1eV)−0.5σ(0.3eV)。它不是自动完成了k网格、尺寸、构型和散射模型收敛。

在本例kBT=0.01Ry条件，有限频率σ为正，但原生首两点外推DC约−28502.1S/m。这个负值表明本次外推不可作物理直流值，不能解释成材料有负电导，也不能把数据截成0后宣称结果合理。

module load ase/3.29.0
python - <<'PY'
import numpy as np
s=np.loadtxt('reference/cases/al-kg-Te0.01/Onsager.txt')
print('first two energies/eV',s[:2,0])
print('positive finite-frequency sigma',s[:2,1])
print('two-point extrapolate/S per m',1.5*s[0,1]-0.5*s[1,1])
PY

2 从Onsager列提取电子热导#

在文件给出的列定义下，κ=[L22/e²−(L12/e)²/σ]/T，T为费米–狄拉克电子温度。注意文件的L12/e和L22/e²已经按相应单位输出，不要再凭印象乘一个电子电荷。

module load ase/3.29.0
python - <<'PY'
import numpy as np
s=np.loadtxt('reference/cases/al-kg-k4/Onsager.txt')
T=0.1*13.605698*11604.518026
kappa=(s[:,4]-s[:,3]**2/s[:,1])/T
print('electronic T/K',T)
print('native kappa first3',s[:3,2])
print('from printed Onsager columns',kappa[:3])
PY

原生Onsager表一般只有6位有效数字，两个很大的热输运项相减时，末几位会被放大。因此用打印列重建κ，不能强求逐位一致。附件按各列输出精度传播误差界；较低电子温度例中约−1.5×10⁻⁸W/m/K的极小负κ也保留，按接近零时的数值精度解释，不擅自修图。

程序还打印Lorenz number，单位为kB²/e²的倍数；它不是直接以W·Ω·K⁻²输出。当前未收敛DC及人为展宽条件下，不能用一个Lorenz数就认定已验证Wiedemann–Franz关系。

3 用实际速度矩阵和关联函数交叉检查#

OUT.tutorial/vmatrixN.dat保存第N个原生k点所有i<j带对的三方向速度矩阵元模方之和。它是本版本带记录长度的二进制文件，并不是波函数文件。下面读取主要参考的第一个k点，核对32带应有496个带对：

python - <<'PY'
import struct
from pathlib import Path
import numpy as np
p=Path('reference/cases/al-kg-k4/OUT.tutorial/vmatrix1.dat')
b=p.read_bytes()
nbytes=struct.unpack('<i',b[:4])[0]
assert nbytes==struct.unpack('<i',b[-4:])[0] and len(b)==nbytes+8
v=np.frombuffer(b[4:-4],dtype='<f8')
print('band pairs',len(v),'expected',32*31//2)
print('finite nonnegative norms',np.isfinite(v).all() and (v>=0).all())
PY

这段对应本批x86/3.10.1实际记录格式，换版本或平台应重新确认。附件analyze.py用这些矩阵、真实本征值、带k权重的占据和费米能重建若干时间点的c11/c12/c22；再对原生打印关联函数施加解析Gaussian/Lorentzian核，独立作正弦变换，恢复频域Onsager系数。

module load ase/3.29.0
python - <<'PY'
import numpy as np
with np.load('reference/analysis/spectra.npz') as a:
    native=a['al-kg-k4__native'][:,[1,3,4]]
    rebuilt=a['al-kg-k4__independent_Onsager']
    error=np.max(abs(rebuilt-native),axis=0)/np.max(abs(native),axis=0)
    print('relative peak errors of sigma,L12/e,L22/e2',error)
PY

这些交叉检查用于确认文件、权重、单位和提取方法一致，并不替代材料模型与参数收敛。实际生产仍需合适的离子构型、电子温度、k点/带数/展宽及低频极限判断。

原生算法和单位常数见3.10.1 EleCond。