TDEP三声子散射曲面实操
从真实二阶力常数生成X点与非高对称q点的三声子散射曲面，读取原生POV顶点并画图。
阅读本篇目录
本例使用鸢算科技集群实际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 曲面在说明什么#
固定初始声子 (q₁,b₁)后，程序寻找第二个波矢q₂，使指定三条声子分支满足能量条件。当前25.03代码使用 q₃=q₁+q₂，按倒格矢折回布里渊区；两个分支的频率条件分别为：
| 通道 | 当前程序的零面条件 |
|---|---|
| plus | ν₁(q₁)+ν₂(q₂)−ν₃(q₁+q₂)=0 |
| minus | ν₁(q₁)−ν₂(q₂)−ν₃(q₁+q₂)=0 |
minus分支涉及相反方向的参与波矢约定，不要只看“q₁+q₂”就把它当成吸收过程。这里的图展示几何允许区域；有曲面不等于有很强散射，也不直接等于某个寿命。
2 准备输入文件#
几何曲面需要 infile.ucposcar和 infile.forceconstant；加上三阶强度时需要 infile.forceconstant_thirdorder。本例原胞为2原子，完整输入如下：
Si
1.0000000000000000
0.0000000000000000 2.7149999999999999 2.7149999999999999
2.7149999999999999 0.0000000000000000 2.7149999999999999
2.7149999999999999 2.7149999999999999 0.0000000000000000
Si
2
Direct
0.0000000000000000 0.0000000000000000 -0.0000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
力常数较长，使用包内完整文件，不要只复制前几行。查看其头部可确认原子数和实际截断：
head -n 10 inputs/infile.forceconstant
head -n 10 inputs/infile.forceconstant_thirdorder
二阶实际截止到4.502318 Å这一壳层；命令请求5 Å不表示模型里一定出现恰好5 Å的原子对。原子编号、原胞与力常数必须匹配，不能把64原子超胞POSCAR直接替代这里的2原子原胞。
3 先列出X点第4支的实际通道#
mkdir surface-X
cp inputs/infile.* surface-X/
cd surface-X
module load tdep/25.03-cgal-8173m
phasespace_surface --highsymmetrypoint X --nband 4 -qg 4 4 4 --verbose > native.log
tail -n 22 native.log
日志的 plus/minus下面，非空通道一行给出 b₁ b₂ b₃ 三角形数量 幅度汇总；不加 --intensities时最后一项为0，不能因此推断真实散射矩阵元一定为0。
本例X点、密度参数4时，plus的4/1/5通道有240个三角形，minus的4/2/2有1408个。先确认实际存在的通道，再选择画图。
4 选择两个通道输出POV#
用新目录保留每次结果，避免同名文件覆盖。在包含同一组输入的目录中：
phasespace_surface --highsymmetrypoint X -qg 12 12 12 --povray --modespec 4 1 5 0 1 70 4 2 2 1 3 70 > native.log
ls outfile.phasespacesurface.*
head -n 15 outfile.phasespacesurface.ini
--modespec每组6个整数：初始支 第二支 第三支 plus/minus 颜色 不透明度。plus为0，minus为1；颜色1/2/3，不透明度0至100。这里选4/1/5的plus与4/2/2的minus。
密度参数12时，本例两通道分别为576和1904个三角形。图直接读取原生POV的全部顶点，单位转为cycles/Å；没有用手工球面或平滑网格替换它们。两面板保持各自真实坐标范围，不根据图上大小比较散射强弱。
5 也可以指定非高对称q点#
phasespace_surface --qpoint .2 .1 .3 -qg 12 12 12 --povray --modespec 4 1 5 0 1 70 4 1 3 1 3 70 > native.log
ls outfile.phasespacesurface.pov
--qpoint是原胞倒格基矢下的分数坐标，不是Å⁻¹笛卡尔坐标。本例对应Cartesian cycles/Å约 (0.0368324,0.0736648,0)。三角形顶点原生写为Bohr⁻¹，提取时除以0.529177210903转换为Å⁻¹；不要额外乘2π后仍标成cycles/Å。
在这个q点直接沿用X点的minus422时，本例密度4/6没有该通道曲面，POV只含另一张plus曲面；改为实际存在的minus413后可生成完整两面。包中保留两种输入和输出。没有指定通道不等于程序报错，不能靠给空数组补点来“补全”曲面。
6 不安装POV-Ray也能看图#
原生输出为 outfile.phasespacesurface.pov/.ini。包内使用Python读取其真实三角形并画图，集群终端无需图形界面：
cd /path/to/tdep-scattering-surfaces-wave287
bash extract.sh
module load pymatgen/2026.9.24
python plot.py results/data my-results
ls my-results/figures
将目录改成实际解压位置，通过文件传输下载PNG查看即可。需要更进一步理解网格与强度，请继续下一篇：散射曲面的网格与强度。程序用途也可对照官方程序目录。