SIESTA专题

SIESTA碳链输运：电极、有限偏压自洽与电流提取

用8原子电极和32原子器件，逐步完成零偏压、正负偏压、自洽输运、透射与电流曲线。

更新于 2026-11-05约 13 分钟阅读按任务与环境选择适用配置

阅读本篇目录

本例搭建一条均匀碳链：先算周期电极，再算开放边界器件，最后提取透射和电流。每个有限偏压都重新进行TranSIESTA自洽计算；计算出的哈密顿量交给TBtrans，不把零偏压透射直接当作所有偏压的结果。

1. 下载后先提交整套工作流#

将下载包保存为carbon-chain.zip：

unzip carbon-chain.zip -d carbon-chain
cd carbon-chain
sinfo
sbatch job.slurm
squeue -u "$USER"

模板使用CPU分区，提交前确认分区可用。整套脚本最多并发4组、每组2个MPI进程和1线程。环境为SIESTA/TranSIESTA/TBtrans5.4.2、Intel oneAPI2024.2、Intel MPI2021.13和MKL顺序LP64；sisl0.16.4负责读取NC文件、提取数据和绘图。运行成功后查看：

ls cases
cat analysis/cases.csv
ls figures

下面的分步命令用于理解文件传递，应在已申请的计算资源内执行。自动脚本已经完成同样的步骤，不需要在同一目录重复运行。

2. 准备结构和周期电极#

module purge
module load siesta/5.4.2-8173m-mpi sisl/0.16.4
export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1
python prepare.py
cat cases/lead-k64/input.fdf

电极包含8个C原子，C–C间距1.27Å，沿x方向周期长度10.16Å，另外两个方向为12Å。器件包含32个C，长度40.64Å，两端各8个原子对应电极区域。电极主层必须足够长以容纳基组轨道的耦合范围，不能只为缩小算例任意改成两个原子。

电极和器件结构，以及电极k点对照

展开：64k电极的完整input.fdf

SystemLabel electrode
SystemName Own uniform carbon chain NEGF tutorial
NumberOfAtoms 8
NumberOfSpecies 1
LatticeConstant 1.0 Ang
%block ChemicalSpeciesLabel
1 6 C
%endblock ChemicalSpeciesLabel
%block LatticeVectors
10.1600000000 0.0 0.0
0.0 12.0 0.0
0.0 0.0 12.0
%endblock LatticeVectors
AtomicCoordinatesFormat Ang
%block AtomicCoordinatesAndAtomicSpecies
0.0000000000 6.0000000000 6.0000000000 1
1.2700000000 6.0000000000 6.0000000000 1
2.5400000000 6.0000000000 6.0000000000 1
3.8100000000 6.0000000000 6.0000000000 1
5.0800000000 6.0000000000 6.0000000000 1
6.3500000000 6.0000000000 6.0000000000 1
7.6200000000 6.0000000000 6.0000000000 1
8.8900000000 6.0000000000 6.0000000000 1
%endblock AtomicCoordinatesAndAtomicSpecies
PAO.BasisSize DZP
PAO.EnergyShift 200 meV
MeshCutoff 150 Ry
XC.functional GGA
XC.authors PBE
ElectronicTemperature 300 K
MaxSCFIterations 500
MinSCFIterations 3
SCF.DM.Converge true
SCF.DM.Tolerance 1.e-7
SCF.FreeE.Converge true
SCF.FreeE.Tolerance 1.e-6 eV
DM.MixingWeight 0.05
DM.NumberPulay 6
DM.UseSaveDM false
WriteDM true
Save.HS true
TS.HS.Save true
TS.DE.Save true
%block kgrid_Monkhorst_Pack
64 0 0 0.0
0 1 0 0.0
0 0 1 0.0
%endblock kgrid_Monkhorst_Pack

cd cases/lead-k64
mpirun -np 2 siesta --electrode < input.fdf > siesta.out
ls electrode.HSX
grep 'siesta: E_KS(eV)' siesta.out | tail -1
cd ../..

这里PBE/DZP、网格150Ry、PAO能量位移200meV、电子温度300K。输入密度矩阵容差为1.e-7，能量收敛为1.e-6 eV，用于稳定后面的有限偏压自洽。C.psml由脚本从集群已部署的PseudoDojo0.4 PBE SR目录复制，资料包不分发赝势。

另一个lead-k32目录只把沿链k点数改为32。两次电极总能量分别约−1298.360152和−1298.362000eV；后面的器件统一使用64k电极。两个网格只是操作对照，不能据此宣称全部参数收敛。

3. 先完成零偏压器件#

cd cases/bias0
mpirun -np 2 siesta < input.fdf > siesta.out
ls device.DM device.TSDE device.TS.HSX
mpirun -np 2 tbtrans < input.fdf > tbtrans.out
ls device.TBT.nc
cd ../..

电极使用普通DFT输出electrode.HSX；TBtrans使用器件自洽后的device.TS.HSX。两类文件用途不同。器件采用SolutionMethod transiesta，不是普通周期DFT。

展开：零偏压器件的完整input.fdf（含32个原子坐标和积分轮廓）

SystemLabel device
SystemName Own uniform carbon chain NEGF tutorial
NumberOfAtoms 32
NumberOfSpecies 1
LatticeConstant 1.0 Ang
%block ChemicalSpeciesLabel
1 6 C
%endblock ChemicalSpeciesLabel
%block LatticeVectors
40.6400000000 0.0 0.0
0.0 12.0 0.0
0.0 0.0 12.0
%endblock LatticeVectors
AtomicCoordinatesFormat Ang
%block AtomicCoordinatesAndAtomicSpecies
0.0000000000 6.0000000000 6.0000000000 1
1.2700000000 6.0000000000 6.0000000000 1
2.5400000000 6.0000000000 6.0000000000 1
3.8100000000 6.0000000000 6.0000000000 1
5.0800000000 6.0000000000 6.0000000000 1
6.3500000000 6.0000000000 6.0000000000 1
7.6200000000 6.0000000000 6.0000000000 1
8.8900000000 6.0000000000 6.0000000000 1
10.1600000000 6.0000000000 6.0000000000 1
11.4300000000 6.0000000000 6.0000000000 1
12.7000000000 6.0000000000 6.0000000000 1
13.9700000000 6.0000000000 6.0000000000 1
15.2400000000 6.0000000000 6.0000000000 1
16.5100000000 6.0000000000 6.0000000000 1
17.7800000000 6.0000000000 6.0000000000 1
19.0500000000 6.0000000000 6.0000000000 1
20.3200000000 6.0000000000 6.0000000000 1
21.5900000000 6.0000000000 6.0000000000 1
22.8600000000 6.0000000000 6.0000000000 1
24.1300000000 6.0000000000 6.0000000000 1
25.4000000000 6.0000000000 6.0000000000 1
26.6700000000 6.0000000000 6.0000000000 1
27.9400000000 6.0000000000 6.0000000000 1
29.2100000000 6.0000000000 6.0000000000 1
30.4800000000 6.0000000000 6.0000000000 1
31.7500000000 6.0000000000 6.0000000000 1
33.0200000000 6.0000000000 6.0000000000 1
34.2900000000 6.0000000000 6.0000000000 1
35.5600000000 6.0000000000 6.0000000000 1
36.8300000000 6.0000000000 6.0000000000 1
38.1000000000 6.0000000000 6.0000000000 1
39.3700000000 6.0000000000 6.0000000000 1
%endblock AtomicCoordinatesAndAtomicSpecies
PAO.BasisSize DZP
PAO.EnergyShift 200 meV
MeshCutoff 150 Ry
XC.functional GGA
XC.authors PBE
ElectronicTemperature 300 K
MaxSCFIterations 500
MinSCFIterations 3
SCF.DM.Converge true
SCF.DM.Tolerance 1.e-7
SCF.FreeE.Converge true
SCF.FreeE.Tolerance 1.e-6 eV
DM.MixingWeight 0.05
DM.NumberPulay 6
DM.UseSaveDM false
WriteDM true
Save.HS true
TS.HS.Save true
TS.DE.Save true
SolutionMethod transiesta
TS.Voltage 0.0 eV
%block TS.ChemPots
high
low
%endblock TS.ChemPots
%block TS.ChemPot.high
mu V
contour.eq
begin
C-high
T-high
end
%endblock TS.ChemPot.high
%block TS.ChemPot.low
mu 0.0 eV
contour.eq
begin
C-low
T-low
end
%endblock TS.ChemPot.low
%block TS.Contour.C-high
part circle
from -40.0 eV + V to -10 kT + V
points 34
method g-legendre
%endblock TS.Contour.C-high
%block TS.Contour.T-high
part tail
from prev to inf
points 12
method g-fermi
%endblock TS.Contour.T-high
%block TS.Contour.C-low
part circle
from -40.0 eV to -10 kT
points 34
method g-legendre
%endblock TS.Contour.C-low
%block TS.Contour.T-low
part tail
from prev to inf
points 12
method g-fermi
%endblock TS.Contour.T-low
%block TS.Contours.nEq
neq
%endblock TS.Contours.nEq
%block TS.Contour.nEq.neq
part line
from -|V| - 5 kT to |V| + 5 kT
delta 0.01 eV
method mid-rule
%endblock TS.Contour.nEq.neq
%block TS.Elecs
high
low
%endblock TS.Elecs
%block TS.Elec.high
HS ../lead-k64/electrode.HSX
chemical-potential high
semi-inf-direction -a1
electrode-position 1
%endblock TS.Elec.high
%block TS.Elec.low
HS ../lead-k64/electrode.HSX
chemical-potential low
semi-inf-direction +a1
electrode-position end -1
%endblock TS.Elec.low
TBT.HS device.TS.HSX
TBT.T.All true
TBT.DOS.A true
TBT.DOS.Gf true
TBT.T.Eig 12
%block TBT.Contours
line
%endblock TBT.Contours
%block TBT.Contour.line
part line
from -1.5 eV to 1.5 eV
delta 0.01 eV
method mid-rule
%endblock TBT.Contour.line

TS.Elec.high/low从同一个电极读取左右边界，分别沿-a1和+a1延伸。器件坐标必须与接触区域的电极排列匹配。相对路径依赖整个cases目录结构，下载后不要只移动单个器件文件夹。

4. 复制零偏压状态，再算正负偏压#

cp cases/bias0/device.DM cases/bias-plus/
cp cases/bias0/device.TSDE cases/bias-plus/
cp cases/bias0/device.DM cases/bias-minus/
cp cases/bias0/device.TSDE cases/bias-minus/
grep -E 'TS.Voltage|DM.UseSaveDM' cases/bias-plus/input.fdf

bias-plus使用TS.Voltage 0.1 eV，bias-minus使用-0.1 eV；对于单电子，这对应电压差±0.1V。两个非零偏压目录都设DM.UseSaveDM true，从零偏压DM和TSDE开始。不能从没有零偏压状态的空目录直接开始有限偏压。

cd cases/bias-plus
mpirun -np 2 siesta < input.fdf > siesta.out
mpirun -np 2 tbtrans < input.fdf > tbtrans.out
cd ../bias-minus
mpirun -np 2 siesta < input.fdf > siesta.out
mpirun -np 2 tbtrans < input.fdf > tbtrans.out
cd ../..

TS.Contour.nEq.neq定义非平衡积分，TBT.Contour.line定义−1.5至1.5eV、步长0.01eV的后处理能量网格。改变偏压必须重新自洽，不是只修改画图横坐标。检查末尾电子迭代、密度矩阵变化和程序正常结束，再提取电流；偏压下打印的能量迭代量不当作平衡自由能比较。

5. 提取透射、本征通道和电流#

python analyze.py
cat analysis/cases.csv
ls analysis/*-transmission.csv
ls figures

程序用sisl读取TBtrans的NC数据，提取左右透射、通道本征值、谱密度和电流，并独立进行费米分布积分。最小读取形式：

import sisl
f = sisl.get_sile('cases/bias-plus/device.TBT.nc')
T = f.transmission(0, 1)
I = f.current(0, 1)
print(T.shape, I)

本例为非磁计算，NC/sisl读出的电流按单自旋块计。物理总电流应将两个相同自旋块相加，即乘2；自旋极化计算应分别提取并求和，不能照搬乘2。

偏压/V	费米能附近透射/单自旋块	单自旋块电流/μA	非磁总电流/μA
−0.1	1.999061	−7.744445	−15.488890
0	1.999061	0	0
+0.1	1.999055	+7.744444	+15.488888

各偏压自洽后的透射和非磁总电流曲线

接近2的单自旋透射来自该模型的导电通道数，不能直接理解为自旋简并已经乘过2。图中的电流已经相加两个自旋块；保存CSV时要同时记录单位和自旋口径，避免后续重复乘2。

6. 换成自己的体系时先检查什么#

先分别优化或选择合理的电极与器件结构，保证接触区域原子和基组匹配；再检查电极k点、基组、实空间网格、屏蔽区长度与积分轮廓。新偏压优先从已收敛的邻近偏压状态开始，逐个检查收敛，再汇总I–V。

本例原子位置固定，只有三个偏压和两个电极k点，目的是演示可操作的自洽输运流程，不给出真实碳链稳定性、接触电阻或收敛材料I–V结论。更多背景见SIESTA官方碳链教程与5.4输入参考。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。