sisl 声子模型：质量约化与单位

使用集群 sisl/0.16.4，环境为 Python 3.11、NumPy、SciPy、netCDF4 和 Matplotlib，不需要加载 SIESTA 的 MPI 编译环境。本文是明确给定参数的教学模型，在集群实际运行并与独立公式比较；不是新的材料 DFT 或 TBtrans 自洽计算。

下载本组完整脚本、原生矩阵、CSV和配图。解压后进入 sisl-core-model-examples：

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

参考数据在 reference-results、CSV在 reference-data。要运行全部短模型，用 bash run.sh，新结果写入 work/results，已完成小项自动跳过；图由参考数据重画到 figures。脚本限制数学库为1线程。

这个例子计算什么#

一个周期位点、质量28.085 u、相邻间距2 Å，纵向弹簧常数5 eV/Ų。只连接x方向，y/z是未加约束的自由零分量。这是用于理解矩阵接口的一维谐振模型，不是硅晶体DFT声子谱，也不能用另外两个零分量评价真实材料动力学稳定性。

为什么不能把力常数直接写进去#

sisl的 DynamicalMatrix 保存质量约化的矩阵，eigh()返回声子能量的平方；eigenmode().hw再取带符号平方根，单位为eV。把eV/Ų力常数直接当eV²赋值，会让整个频率尺度错误。

本例用下列转换把K/M变为eV²，完整常数写在脚本中：

cat > phonon.py <<'EOF'
import sisl
import numpy as np
hbar = 6.582119569e-16  # eV s
unit = 1.602176634e-19 / (1e-20 * 1.66053906660e-27) * hbar**2
K, mass = 5.0, 28.085
C = K / mass * unit
atom = sisl.Atom(14, R=[2.01,2.01,2.01])
g = sisl.Geometry([[0,0,0]], atom,
    lattice=sisl.Lattice([2,10,10], nsc=[3,1,1]))
D = sisl.DynamicalMatrix(g)
D[0,0] = 2*C
D[0,0,(1,0,0)] = -C
D[0,0,(-1,0,0)] = -C
for k in [0.0,0.25,0.5]:
    hw = D.eigenmode([k,0,0]).hw
    print(k, hw, 'THz:', hw / (2*np.pi*hbar) / 1e12)
EOF
python phonon.py

一个位点设置3个分量作为x/y/z自由度。这里下标0为x分量，左右胞耦合各−C，对角元2C；因此Γ点总和为0，满足这一简单模型的平移和规则。

实际输出#

实际计算21个k点，Γ点纵向能量为0，区边界k=0.5得到 0.05455997959 eV = 54.55997959 meV = 13.19254437 THz。

独立解析式为 ħω(k)=sqrt(2C[1−cos(2πk)])，与原生输出一致。频率ν=ω/(2π)，所以从ħω换为THz要除以2πħ；直接除以ħ得到角频率，二者相差2π。

head -n 8 reference-data/Sisl-phonon-chain.csv
bash plot.sh

CSV同时保存eV、THz、解析结果、质量和弹簧常数，便于改参数后对照。增加质量而保持弹簧常数，频率按1/√M减小；这里不是通过重跑DFT得到同位素力常数。

迁移到真实声子数据时#

先确认输入是未约化力常数还是质量约化矩阵、质量单位、坐标顺序、超胞映射与声学和规则，再传给sisl。若已有Phonopy/DFPT结果，优先使用对应的正式转换工具；不能仅因为矩阵尺寸相同就直接赋值。DynamicalMatrix.eigh()的eV²与eigenmode().hw的eV也不能混在一列画图。

真实晶体计算还应检查力常数与超胞/采样收敛。本文用于理解接口和单位，不替代材料声子计算的前处理。