sisl终端处理结构与密度网格

不必每次都写Python脚本，sisl也提供处理结构和网格的终端命令。本篇读取鸢算科技集群实际SIESTA5.4.2的Si2输出，用sisl0.16.4运行；环境为Python3.11及NumPy/SciPy/Matplotlib，不需要加载SIESTA MPI编译环境，也不重新计算DFT。

1. 下载和运行完整示例#

下载实际Si2文件、命令与参考结果，保存为sisl-cli.zip：

unzip sisl-cli.zip -d sisl-cli
cd sisl-cli
ls inputs reference-results reference-figures
module load sisl/0.16.4
sgeom --help
sgrid --help

reference-results/、reference-analysis/、reference-figures/是集群实际结果。要在新目录重新运行全部示例，可在可用CPU分区提交：

sinfo
sbatch job.slurm

模板按集群CPU分区56物理核整节点规则申请，工具顺序以1线程运行。若已在分配到的CPU作业中，直接bash run.sh。新结果写入results/、analysis/和figures/，不覆盖参考。

2. sgeom：转换、超胞与原子顺序#

module load sisl/0.16.4
mkdir work
sgeom inputs/si.XV --format .14f --out work/si.xyz
sgeom inputs/si.XV --format .14f --tile 2 0 --tile 2 1 --tile 2 2 --out work/si-222.xyz
sgeom inputs/si.XV --format .14f --repeat 2 0 --repeat 2 1 --repeat 2 2 --out work/si-repeat.xyz

XV中有实际晶胞、原子种类和坐标。sgeom识别其格式，写出的XYZ保留晶胞信息。本例2×2×2从2原子变成16原子，晶胞体积变为8倍。tile和repeat得到相同的坐标集合，但排列顺序不同；如果后续文件按原子序号关联数据，不能只看结构图相同就直接替换。

3. 平移、旋转与逆变换#

sgeom inputs/si.XV --format .14f --translate 1,2,3 --out work/si-shift.xyz
sgeom inputs/si.XV --format .14f --rotate-z 90 --out work/si-rotated.xyz
sgeom work/si-222.xyz --format .14f --untile 2 0 --untile 2 1 --untile 2 2 --out work/si-untiled.xyz

平移量为笛卡尔坐标Å；旋转默认角度制，晶胞矢量也随之旋转。untile针对由相应tile生成的超胞还原，本例恢复原始2原子坐标。它不等于对任意缺陷超胞自动识别原胞。

需要Python时，同一文件也能直接读取：

python - <<'PY'
import sisl
g = sisl.get_sile('work/si-222.xyz').read_geometry()
print('atoms:', g.na)
print('cell / Angstrom:', g.cell)
print('volume / Angstrom^3:', g.volume)
PY

4. sgrid：缩放、相减和复制#

sgrid inputs/si.RHO --geometry inputs/si.XV --out work/rho.cube
sgrid inputs/si.RHO --scale 2 --out work/double.RHO
sgrid inputs/si.RHO --diff inputs/si.RHO --out work/zero.RHO
sgrid inputs/si.RHO --tile 2 0 --out work/tiled.RHO

本例RHO网格为24³，sisl读入后的密度单位为e/ų，原始积分为8.000000162个价电子。乘2后约16；自减为零；沿方向0复制两份后晶胞与网格扩展，总积分约16。

python - <<'PY'
import sisl
for name in ['inputs/si.RHO','work/double.RHO','work/zero.RHO','work/tiled.RHO']:
    g = sisl.get_sile(name).read_grid()
    print(name, g.shape, g.grid.sum()*g.dvolume)
PY

这里的CUBE用于保存一般标量场，sisl写出后数值仍为e/ų，不会仅因为扩展名变成.cube就自动乘上密度的Bohr单位转换。交给要求e/Bohr³的外部电荷分析程序前，应另做单位转换；不要直接把此文件与原子单位密度CUBE混用。

5. average、sum为什么不能只看文件名#

sgrid inputs/si.RHO --average 0 --out work/average.RHO
sgrid inputs/si.RHO --sum 0 --out work/sum.RHO

安装版sisl0.16.4的这两种操作都会把方向0的网格数变成1，对应晶胞方向也缩为原来的1/24。因此：

如果要把平均场理解为原晶胞上的投影，需使用原体积解释它，而不能直接套用缩小后的文件体积。本例average的均值乘原晶胞体积约为8.000000954。图中的曲线从原网格沿另外两个分数坐标方向求平均，不把结果当成重新计算的三维密度。

6. 平滑与插值后必须再核对积分#

sgrid inputs/si.RHO --smooth 0.2 gaussian wrap --out work/smoothed.RHO
sgrid inputs/si.RHO --interp 48 48 48 1 wrap --out work/interpolated.RHO

wrap按周期边界处理。这个小例子中，周期Gaussian平滑后积分为8.000000869；24³线性插到48³后为7.882768320，相对原密度减少约1.47%。插值不是自动守恒，也不能替代DFT重新计算更细的网格。

python - <<'PY'
import sisl
for name in ['inputs/si.RHO','work/smoothed.RHO','work/interpolated.RHO']:
    g = sisl.get_sile(name).read_grid()
    print(name, g.grid.sum()*g.dvolume)
PY

需要定量电荷分区或高精度积分时，优先从原计算程序增加网格精度并做收敛比较；不要为了接近预期电子数而悄悄归一化插值结果。继续看从实际EIG/KP提取DOS。

命令入口也可参考sisl官方CLI说明，具体参数与行为以上述安装版0.16.4的帮助和实际文件为准。