Phonopy 热位移矩阵与统计实操
提取笛卡尔和CIF热位移矩阵、方向投影、声子三阶矩，比较量子与经典谐统计。
阅读本篇目录
下载输入、数值参考与绘图脚本。在已分配的CPU资源中新建目录：
unzip phonopy-masses-modes-tensors-wave245.zip -d phonopy-modes
cd phonopy-modes
bash run.sh
run.sh加载phonopy/4.8.3，依次执行run.py、analyze.py、plot.py，生成source/、analysis/和figures/；同名结果目录存在时不覆盖，重做请换新的解压目录。使用Python3.11.16、NumPy2.3.5，C/OpenMP部分由GCC13.3构建，当前默认Rust数值后端。每个程序使用1个CPU线程，不需要GPU。
本批复用鸢算科技集群实际Si16的VASP6.6.1 DFPT力常数，2Si原胞、2×2×2超胞，ENCUT=450 eV、EDIFF=1E-6、IBRION=8，电子k网格2×2×2。与其他教程中的550 eV有限位移数据分开。本包只含安全输入、真实数值力常数和后处理脚本，运行时不启动VASP。保留原始力常数，不施加声学和规则或额外对称化，不把微小负频改成正频。下面的结果用于操作参考，正式材料参数另做收敛。
1 输出热位移矩阵#
cd source/thermal-tensor
phonopy ../reference.yaml --no-fc-symmetry \
--mesh '8 8 8' --nomeshsym --gc --tdm \
--fmin 0.1 --tmin 0 --tmax 600 --tstep 300
cd ../..
cat analysis/thermal-tensors.csv
--tdm生成thermal_displacement_matrices.yaml，在0、300、600 K计算每个原子的热均方位移矩阵，单位Ų。--fmin 0.1排除低于0.1 THz的模式，包括本例Γ点三个平移模式；它是统计截断，不是给负频改号。--nomeshsym使本征矢数据可用于完整热位移求和，--gc采用Γ中心网格。
原生YAML以xx、yy、zz、yz、xz、xy排列六个独立分量，保留5位小数；analysis/thermal-tensors.csv和NPZ另存高精度API数值。两个等价Si的笛卡尔矩阵相同，各方向均方位移为：
| 温度/K | Uxx=Uyy=Uzz/Ų | 单方向RMS/Å |
|---|---|---|
| 0 | 0.00265501 | 0.05153 |
| 300 | 0.00843072 | 0.09182 |
| 600 | 0.01618281 | 0.12721 |
RMS是均方位移的平方根，不能把Ų直接写成Å；0 K非零来自零点运动。
2 方向投影与CIF导出#
cat source/thermal-projection/thermal_displacements.yaml
cat source/thermal-cif/tdispmat.cif
投影调用在相同网格、温度和下限基础上使用--td --pd '1 1 0'；CIF调用使用--tdm-cif 300，输出300 K的矩阵及tdispmat.cif。
方向投影应与dᵀUd相符，d是归一化方向。本例Si的笛卡尔矩阵各向同性，因此沿(1,1,0)的均方位移与Uxx相同。
**CIF矩阵和笛卡尔矩阵不能逐分量混用。**本例原胞基矢夹角60°，300 K时笛卡尔非对角分量近零，CIF非对角分量却约−0.00281024 Ų，对角为0.00843072 Ų。这是非正交晶胞的坐标表示效果，不表示硅在笛卡尔空间出现了剪切型热运动或各向异性。负的非对角项也不意味着方差为负，应看完整矩阵及其物理坐标定义。
3 提取声子频率分布的三阶矩#
cat source/moment3/native.out
本例调用--moment --moment-order 3 --fmin 0.1，同样使用8×8×8网格。原生输出：
Order| Total | Projected to atoms
3 |1511.37413 | 1511.37413 1511.37413
这里是满足频率范围的声子态的归一化频率三阶矩，单位THz³；不是三阶力常数，不包含三声子散射，也不能当作非谐热导率。两个Si等价，因此两个原子投影相同。需要其他阶数时改--moment-order，比较时保持频率范围一致。
4 比较量子与经典谐统计#
head -n 8 analysis/classical-quantum.csv
cat source/thermal-quantum/thermal_properties.yaml
cat source/thermal-classical/thermal_properties.yaml
两组调用使用相同力常数、8×8×8网格和0.1 THz下限，温度100–2000 K、步长100 K。第二组增加--classical。300 K时量子Cv约39.73750 J·mol⁻¹·K⁻¹；经典值49.8381172 J·mol⁻¹·K⁻¹，不随温度改变，量子曲线在高温趋近经典曲线。
这里“每摩尔”是每摩尔2Si原胞，若按每摩尔Si原子应除以2。由于有限8³网格上三个Γ平移模被排除，经典保留模式数为6−3/512，结果略低于6R；不能为了得到整齐的6R把这三个模式硬加回统计。更密网格会改变这项有限网格差别。
--classical只改变统计表达式，不产生MD轨迹；两条曲线都采用固定谐力常数，不包含热膨胀或完整非谐效应。若需要热膨胀，看原生QHA热膨胀实操；若需要结构模式、质量修改和动画，看质量替换与模式结构。设置含义参考Phonopy官方说明，使用时核对本机版本。
继续：可选拟合后端与温度相关流程#
ALM力常数拟合实操提供真实位移受力、type-2数据和同条件symfc对照；多项式势训练与声子演示48/8/8划分、模型保存、独立误差及原生求力；MLP与SSCHA实操给出两组实际300K短流程、误差和数据对应。可选依赖在个人环境创建，NEP标签和DFT标签来源清楚区分，热采样范围与收敛要求见正文。
方向参数的坐标定义#
--pd和Python接口direction使用原胞实空间基矢下的约化坐标。与笛卡尔热位移矩阵做dᵀUd比较时，先计算cart_direction = reduced_direction @ phonon.primitive.cell，再归一化该笛卡尔向量。本例原胞基矢不正交，约化方向(1,1,0)对应笛卡尔方向比例(1,1,2)，不能直接把两个三元组当作相同坐标。Si的笛卡尔矩阵各向同性，因此本例不同方向的投影恰好相同；换成各向异性晶体时必须先做这一坐标变换。