TBtrans能量网格与电流核对

算完一条I–V曲线后，应先检查电流如何从透射积分得到。本篇在相同器件哈密顿量上比较两种TBtrans实轴网格，再用独立的费米函数积分核对sisl输出。

下载完整输入、实际NC/HSX数据和绘图脚本。解压进入资料包目录，先画鸢算科技集群的实际参考结果：

bash view.sh

figures下得到结构与电极k点、七个偏压的透射/I–V，以及两种能量网格的电流比较三张图；analysis/cases.csv和各偏压透射CSV保存数值。参考数据位于 reference/cases，无需先重新计算。

环境为SIESTA/TranSIESTA/TBtrans5.4.2，Intel oneAPI2024.2、Intel MPI2021.13、MKL顺序LP64；后处理使用sisl0.16.4、NumPy、SciPy和Matplotlib。SIESTA/TBtrans在本例采用CPU路径。新计算先确认CPU队列可用，再提交：

sinfo -p cpu-8173m
sbatch job.slurm
squeue -u "$USER"

模板按CPU队列整节点规则申请56核，工作流最多同时运行4个条件、每个条件2个MPI进程、每进程1线程，实际最多使用8个物理核心。模块名称说明构建目标，本文数据不用于特定CPU节点的性能比较。请在分配的计算资源内运行完整DFT流程。

1 分清k点和能量采样#

电极的32/64k控制纵向布里渊区采样；TBtrans的能量网格用于透射和电流积分。它们是两种不同的数值参数，不能只提高一个就认为另一个已经收敛。

本例实轴积分窗口为−1.5至+1.5 eV，比较0.010 eV与0.005 eV两种中点网格，实际分别有300和600个能量点。只改变下面TBtrans段的delta：

%block TBT.Contour.line
part line
from -1.5 eV to 1.5 eV
delta 0.005 eV
method mid-rule
%endblock TBT.Contour.line

TranSIESTA自洽过程中使用的非平衡积分轮廓保持原设置，这次没有随之细化。完整工作流及FDF见七偏压I–V教程。

2 保留两份NC文件#

TBtrans发现同名 device.TBT.nc时会停止，不能用同一个文件名连续覆盖。脚本完成粗网格后，将其移动为 device-coarse.TBT.nc，然后运行细网格，新的结果保留为 device.TBT.nc。

在已经分配的计算资源中，程序阶段依次读取 input.fdf和 input-fine.fdf；run.py负责MPI与绑核。两份NC共用同目录的 device.TS.HSX，重新计算TBtrans能量网格不必重复求解这一偏压的DFT。

ls reference/cases/bias-plus-0p20/*TBT.nc

3 独立读出网格、化学势和电流#

module load sisl/0.16.4
cd reference/cases/bias-plus-0p20
python - <<'PY'
import sisl
import numpy as np
from scipy.special import expit
from scipy.constants import elementary_charge,Planck
currents=[]
for filename in ["device-coarse.TBT.nc","device.TBT.nc"]:
    f=sisl.get_sile(filename)
    E=f.E
    T=f.transmission("high","low")
    dE=float(np.diff(E)[0])
    muL,muR=f.mu("high"),f.mu("low")
    kTL,kTR=f.kT("high"),f.kT("low")
    manual=np.sum(T*dE*(expit(-(E-muL)/kTL)-expit(-(E-muR)/kTR)))*elementary_charge**2/Planck
    native=float(f.current("high","low"))
    reverse=float(f.current("low","high"))
    print(filename,"points=",len(E),"dE/eV=",dE)
    print("Single spin block: API/manual/reverse current (A):",native,manual,reverse)
    print("Nonmagnetic total current (microA):",2*native*1e6)
    currents.append(native)
print("Relative mesh change:",abs(currents[1]-currents[0])/max(abs(currents[1]),1e-12))
PY

脚本分别打印点数、相邻能量间隔、sisl电流、独立积分值和反向电流。E、mu与 kT在这里均使用eV；kT是热能，并非以K为单位的温度。上述API和独立积分均为单自旋块电流。积分乘 e²/h得到安培；非磁性总电流再将两个相同自旋块相加。

独立积分采用文件中实际的中点均匀网格。换成非均匀或自适应网格时，应使用相应积分权重，不能直接沿用统一dE。

4 观察本组网格敏感性#

本组七个偏压点中，两种网格的相对电流差最大约为 1.26×10⁻⁸。+0.20 V时，非磁性总电流分别为30.97762194与30.97762227 μA。零偏压两者均为0，表格和图保留原始结果。

这说明本例实轴后处理积分对这两种步长不敏感。窄共振、低温或更大偏压的器件需要重新选择窗口和能量分辨率，不能将这里的微小差值当作任意体系的误差保证。

5 电流异常时按顺序检查#

先从NC读取两个电极的实际化学势，确认差值与目标偏压一致；再确认该偏压使用自己的自洽HSX。然后检查窗口是否覆盖两端费米函数有差异的区域，比较细化网格后的电流，并核对电流方向和自旋约定。

同一个两电极结果中，high→low和low→high电流应大小相近、符号相反。本例同时检查正反透射、本征通道总和及独立费米积分。若这些量不一致，应回到文件、单位与读取方式核对；仅凭曲线平滑不足以判断计算可信。

参考：TBtrans官方手册、sisl的TBtrans文件接口。