sisl 轨道、电子密度与终端作图

使用集群 sisl/0.16.4 的 Python 3.11 环境、NumPy/SciPy、netCDF4 和 Matplotlib。本组是实际执行的解析轨道/网格与几何操作，不是新的SIESTA或外部ATOM自洽计算。

下载完整脚本、原生数据、网格及参考图，解压后进入 sisl-orbitals-grid-examples。

module load sisl/0.16.4
bash results.sh
bash plot.sh

新目录中执行 bash run.sh 运行全部本组短模型，新结果写入 work，参考文件保持在 reference-results。脚本以单线程运行，不提交新的DFT/MD作业。

生成一个解析1s轨道#

sisl的 HydrogenicOrbital由量子数和有效核电荷给出解析函数。下面的H1s是一电子氢模型，不需要赝势或SCF：

cat > radial.py <<'EOF'
import sisl
import numpy as np
from scipy.integrate import trapezoid
a0 = float(sisl.units('Bohr','Ang'))
r_bohr = np.linspace(0,150,15001)
orb = sisl.HydrogenicOrbital(1,0,0,1.0,R=120.0)
R_bohr = orb.radial(r_bohr*a0) * a0**1.5
u = r_bohr * R_bohr
print('integral |u|^2 dr:',trapezoid(abs(u)**2,r_bohr))
np.savetxt('H1s-u.dat',np.c_[r_bohr,u])
EOF
python radial.py

radial输入长度为Å；为与Bohr单位的独立公式比较，振幅还需乘a₀^(3/2)。归一化的是 ∫|R|²r²dr 或 ∫|u|²dr，不是直接积分R²。

本例实际计算1s、2p、3d、4f四个氢样轨道，并与独立Laguerre多项式公式比较。四者的径向积分均接近1。各个轨道的整体正负相位没有物理区别，比较振幅时需要考虑这一点。

这些轨道分别是不同激发量子数的数学基函数，不表示一个氢原子同时在四个轨道各放一个电子。

用stoolbox绘制径向文件#

资料包提供 AEWFNR0到AEWFNR3样式文件，内容就是上述实际解析径向函数u(r)，文件头明确注明来源。它们不是外部ATOM计算输出，也不是本批生成的赝势；这里只用它们演示sisl已有的读取/作图接口。

cd reference-results/hydrogen-analytic
stoolbox atom-plot --plot wavefunction -l spdf --save hydrogen.png INP
cd ../..

实际原生命令已执行并生成PNG。当前文件只有解析全电子氢样轨道，没有伪波函数，程序提示缺少PSWFNR文件是预期的；不能据此评价赝势质量。INP用于提供元素和通道元信息，没有被提交给ATOM程序执行。

若你有自己通过合法ATOM安装算出的真实输出，命令相同，改到对应输出目录执行即可。默认还可画charge/log/potential，但需要相应真实输出文件；不能拿本包的解析数据补成不存在的赝势诊断结果。

从密度矩阵得到实空间密度#

下面给一个H1s基函数占据1电子，在12 Å立方盒中心计算：

cat > density.py <<'EOF'
import sisl
import numpy as np
orb = sisl.HydrogenicOrbital(1,0,0,1.0,R=4.5)
g = sisl.Geometry([[6,6,6]], sisl.Atom(1,[orb]),
    lattice=sisl.Lattice([12,12,12],nsc=[1,1,1]))
DM = sisl.DensityMatrix(g)
DM[0,0] = 1.0
grid = sisl.Grid(0.12,geometry=g)
DM.density(grid)
print('grid:',grid.shape)
print('electrons:',np.sum(grid.grid)*grid.dvolume)
grid.write('hydrogen-density.grid.nc')
EOF
python density.py

实际网格为 100×100×100，积分得到 1.0004397925电子。这是有限网格积分的结果，与1电子相差约0.044%；不能把这点数值误差解释成体系额外带了电。

密度来自指定的解析基函数和密度矩阵，并没有运行原子DFT。矩阵占据从1改为2会改变数学密度，不能据此认为物理氢负离子的电子结构已经自洽求解。

截断和单位要一并记录#

本例轨道截止半径4.5 Å。中心线与解析式比较时，截止之外也要设为0；与无限延伸的氢1s函数比较，尾部自然会有差别。忽略这个约定，会把正确的局域基截断当成密度程序错误。

本组密度单位为电子/ų，长度为Å，积分乘体元ų后得到电子数。这些单位来自本例的构造。当前通用NC写出没有嵌入足够明确的单位标记，sisl读回会提示无法确定单位；本包用 GRID-UNITS.json记录定义，不能把任意同扩展名外来文件都按同一单位处理。

cat GRID-UNITS.json
head -n 8 reference-data/Sisl-Hydrogen-density-line.csv

读取真实SIESTA密度时还需区分原始文件中的Bohr、Å与电子密度变换，见真实结构与密度网格操作。更多矩阵占据与归一化见电子态分析。