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作业。
先创建三个区域#
sisl形状对象可直接判断笛卡尔点是否落入区域，不必先生成原子模型。以下脚本可直接运行：
cat > select.py <<'EOF'
import numpy as np
import sisl
xyz = np.array([[0,0,0],[0.9,0,0],[0,1.8,0],[0,0,2.8],[1.5,0,0]])
sphere = sisl.Sphere(1.0,center=[0,0,0])
box = sisl.Cuboid([2,4,6],center=[0,0,0])
ellipsoid = sisl.Ellipsoid([1,2,3],center=[0,0,0])
for name,shape in [('sphere',sphere),('box',box),('ellipsoid',ellipsoid)]:
mask = shape.within(xyz)
print(name,mask,xyz[mask])
EOF
python select.py
球的1是半径；盒子的 [2,4,6] 是三个完整边长；椭球的 [1,2,3] 是半轴长度。把盒边长也当成半轴，选中的范围会大一倍。
本例所有点和尺寸都按Å定义。形状对象本身不会从一组无标签数字推断你使用的是Å、Bohr还是分数坐标，调用前要统一单位。
本次实际采样与独立几何判断#
完整实例使用15×17×19=4845个自有笛卡尔点。分别与下列独立条件比较，所有采样点的选择结果一致：
| 形状 | 独立条件 |
|---|---|
| 半径1的球 | x²+y²+z²≤1 |
| 2×4×6盒 | |x|≤1、|y|≤2、|z|≤3 |
| 1/2/3半轴椭球 | x²+(y/2)²+(z/3)²≤1 |
这是几何掩码示例，不是原子密度云或电子密度结果。图上点数由离散采样决定，不能直接拿入选点数当成连续体积。
并集、交集与差集#
python - <<'EOF'
import sisl
import numpy as np
xyz = np.array([[0,0,0],[0.9,0,0],[0,1.8,0],[0,0,2.8],[1.5,0,0]])
s = sisl.Sphere(1.,center=[0,0,0])
b = sisl.Cuboid([2,4,6],center=[0,0,0])
e = sisl.Ellipsoid([1,2,3],center=[0,0,0])
print('union:',(s|b).within(xyz))
print('intersection:',(s&e).within(xyz))
print('box minus sphere:',(b-s).within(xyz))
EOF
实际完整采样同时比较了 |、&、- 与NumPy布尔掩码的并、交、差，结果一致。当前球本身落在给定盒/椭球中，所以部分组合会等同于较大或较小区域；换成不同中心后应重新看几何关系。
用到自己的原子结构#
将真实结构中的笛卡尔坐标数组传给 within，得到布尔掩码，再据此取原子。例如对通过sisl读取的几何对象g，可用 g.sub(np.flatnonzero(shape.within(g.xyz)))生成子结构。
周期体系要先决定是否展开超胞、把中心放在哪个周期像。这里的简单区域选择不会自动替你搜索所有周期像；直接在一个原胞中剪球，可能只得到球的一部分。
原子被删除后，原轨道编号、耦合矩阵、原子分组与化学计量也会变化。截取原子结构不等于对电子矩阵完成严格的开边界消去，后续计算需要重新准备相匹配的输入。
head -n 10 reference-data/Sisl-shape-selections.csv
每个区域的实际入选坐标在 reference-results/shape-selections，可以直接检查选了哪些点，再把同一操作应用到自己的结构。