TDEP散射曲面的网格与强度

本例使用鸢算科技集群实际Si64 DFT短轨迹拟合的300 K模型A：二阶请求截断5 Å、三阶4.2 Å。它与TDEP力常数教程来自同一组数据；这里直接读取已有力常数，不重新运行DFT或拟合。

可选模块为 tdep/25.03-cgal-8173m，基于25.03固定提交d38f435，ifx/icx/icpx 2024.2.1、Intel MPI 2021.13、顺序MKL、FFTW 3.3.11与HDF5 1.14.6；增加CGAL 5.6.1、Boost 1.86、GMP 6.3、MPFR 4.2.2。它单独增加 phasespace_surface，原14个程序与默认模块保持原入口。使用兼容AVX512的CPU；前置DFT使用GPU不代表本程序是GPU程序。

下载真实力常数、结构、数值、3张图与操作脚本，保存为 tdep-surfaces.zip：

python3 -m zipfile -e tdep-surfaces.zip .
cd tdep-scattering-surfaces-wave287
bash preview.sh
ls preview/figures

想自己从已有力常数运行原生程序，依次执行：

bash prepare.sh
bash calculate.sh
bash extract.sh

原生结果在 work/surfaces/，提取后的三角形和数值表在 results/data/。包中包含11个很小的单进程CPU条件，其中两个保留无指定通道曲面的例子，帮助理解如何区分“没有该通道”和文件读取失败。正式生产请申请相应Slurm CPU资源；本教程不提供8173M节点性能数据。

1 qgrid不是简单的n³总点数#

phasespace_surface调用CGAL构造不可约楔形区四面体网格，再利用对称性扩展。这里的 -qg 4 4 4等是密度控制参数，不能直接报告为“全区只有64点”。查看真实网格应读日志：

grep -E 'points|tetra|reduced from' work/surfaces/X-q12-geometry/native.log

同一个选定通道在两个密度参数下出现相同三角形数量，也不代表参数未生效；其他区域的网格可以已经变化。本例X点plus415在4/6时均240个，8时352个，12时576个。

2 用原始力常数把频率代回曲面#

原生曲面是四面体内频率差的线性插值零面。对其顶点重新用真实二阶力常数计算声子，能量差通常不会严格为0；这反映插值误差，不能把插值式中强制为0当成独立精度检查。

下载包提供单独的实空间力常数矩阵读取器：先与原生80个q点、6支频率对照，最大差约4.36×10⁻⁷ THz，再计算曲面顶点的真实频率差。它使用原生自然同位素加权Si质量28.0854987064 u，不把质量与频率单位混用。

module load pymatgen/2026.9.24
python - <<'PY'
import json
a=json.load(open('results/data/surface-summary.json'))
for r in a:
    if r['case'].startswith('X-'):
        print(r['case'], r['kind'], r['triangles'], r['RMS_residual_THz'])
PY

残差不是每次加密都单调减小，图形变细也不等于物理量收敛。本例只用于理解曲面生成和检查流程；不能据此认为网格12已经足够准确。超胞、力常数截断和原始短采样误差也还存在。

左下按三角形顶点出现次数统计，重复顶点未去重；它不是体积加权的态密度。图和数值保留实际粗网格差异。

3 实际运行三阶幅度#

先确认同一原胞的三阶文件已准备好，进入相应目录：

ls infile.ucposcar infile.forceconstant infile.forceconstant_thirdorder
phasespace_surface --highsymmetrypoint X --nband 4 -qg 4 4 4 --intensities --verbose > native.log
tail -n 22 native.log

下载包对应 work/surfaces/X-q4-intensities/。本例实际得到非零的原生幅度汇总；例如minus413约7.6879×10⁻⁹，plus426约4.4664×10⁻⁹，具体值以随包JSON为准。

cat results/data/native-amplitude-norms.json

**这个版本打印的 psisq_norm是顶点 |V|²数组的欧氏L2范数。**它不是曲面积分，更不是已乘Bose占据、频率δ函数和积分权重的散射率或寿命。顶点数量改变会影响范数，不能拿不同网格的这一列直接比较材料散射强弱。

4 强度开关改变哪些输出含义#

在同样输入、q点和网格下，本例打开 --intensities前后的曲面顶点完全相同；它读取三阶力常数，计算附加幅度信息。

POV中跟在顶点后面的向量则有不同含义：不加开关时是归一化插值梯度；加开关后，此分支改写为原生 v₁±v₂−v₃。不要把两者当成相同物理量逐值比较，也不能假定后者仍是单位法向。

25.03的逐点强度数组存在plus/minus存储混用，因此这里读取原生标量范数，不把该数组当成可靠的逐三角形着色接口。需要定量散射与寿命时，继续使用热导程序中单位明确的模式结果，并完成相应网格与模型检查。

5 换成自己的体系#

先在相同原胞约定下完成二/三阶拟合及独立力误差检查，再列出目标q点实际通道，选择模式输出曲面，检查顶点能量差和网格敏感性，最后才讨论三阶幅度与输运。出现虚频、不同质量、原子映射或长程极性项时，要先处理相应模型条件，不能直接套用Si的图和阈值。

module load tdep/25.03-cgal-8173m
phasespace_surface --help
module list

本模块增加了基础模块原先缺少的第15个原生程序；15个程序均已有操作例，但不代表所有参数组合或每一种材料都适用。