Phonopy 后端、HDF5与低频精度

下载输入、力常数、转换脚本与绘图资料包。在已分配的 CPU 资源内解压，一份新目录运行一次：

unzip phonopy-special-formats-wave241-*.zip -d phonopy-practice
cd phonopy-practice
bash run.sh
ls analysis figures

run.sh 每次使用 1 个 CPU，不需要 GPU。它完成结构生成、原生格式转换、读回比较和既有 Si 力常数后处理，最后得到 CSV 和三张高清图；不会启动 VASP 或其他电子结构程序。输出目录已存在时，另解压到新目录继续练习。

环境为 phonopy/4.8.3、Python 3.11.16、NumPy 2.3.5、SciPy 1.16.3、spglib 2.7.0；C/OpenMP 部分使用 GCC 13.3，默认 Rust 后端为 phonors 0.6.0。CRYSTAL 外部几何独立读回使用集群的 ASE 3.29.0 环境。后面的短命令在 module load phonopy/4.8.3 后使用。

1. 本例从哪一步开始#

本例使用本集群已经完成的Si16 VASP DFPT数值力常数，不重新求DFT。超胞16个原子、原胞2个原子，源计算为VASP6.5.1 CPU版、IBRION=8、2×2×2电子k网格。完整DFPT工作流见qvasp与Phonopy声子。这是小体系操作参考，不把一次小网格计算称作材料收敛数据。

从 inputs/Si16/phonopy_disp.yaml 和 FORCE_CONSTANTS 读入、生成完整数据文件后，分别用默认Rust和C后端计算相同的6个q点、4×4×4声子网格；保持原始力常数，不对称化、不额外施加ASR、不裁掉负频。

2. 准备同一份数据#

完整执行 bash run.sh 后，也可从初始解压目录单独准备新的手动目录：

mkdir manual-backend
cp work/backends/complete.yaml manual-backend/phonopy.yaml
cd manual-backend

complete.yaml 包含结构与力常数。本例先加载数值文件再保存，所用Python方法：

import phonopy
p = phonopy.load('inputs/Si16/phonopy_disp.yaml',
                 force_constants_filename='inputs/Si16/FORCE_CONSTANTS',
                 is_nac=False, symmetrize_fc=False)
p.save(filename='complete.yaml', settings={'force_constants': True})

这段路径以原始解压目录为起点。在自己的体系中，必须让结构、超胞矩阵、原子顺序与力常数完全匹配；只复制一个FORCE_CONSTANTS而换了结构，程序能运行也未必有物理意义。

3. 默认Rust与C怎么选择#

在新练习目录先运行默认Rust：

phonopy phonopy.yaml --no-fc-symmetry --nonac \
  --qpoints '0 0 0 .5 0 .5 .1 .2 .3 .5 .5 .5 .5 .25 .75 .25 .25 .25' \
  --qpoints-format hdf5 --writedm --eigenvectors \
  --writefc --writefc-format hdf5
mv qpoints.hdf5 qpoints-Rust.hdf5

再指定C后端：

phonopy phonopy.yaml --legacy-backend --no-fc-symmetry --nonac \
  --qpoints '0 0 0 .5 0 .5 .1 .2 .3 .5 .5 .5 .5 .25 .75 .25 .25 .25' \
  --qpoints-format hdf5 --writedm --eigenvectors
mv qpoints.hdf5 qpoints-C.hdf5

4.8.3默认已经使用Rust。--rust在这个版本是弃用的空操作参数，选择C用--legacy-backend。这与VASP用CPU还是GPU没有关系，也不改变源力常数来自哪次电子结构计算。

4. 直接提取HDF5频率#

python - <<'PYCODE'
import h5py
import numpy as np
with h5py.File('qpoints-Rust.hdf5') as h:
    print('datasets:', list(h.keys()))
    q = h['qpoint'][:]
    f = h['frequency'][:]
with h5py.File('qpoints-C.hdf5') as h:
    fc = h['frequency'][:]
print('qpoints:', q.shape)
print('frequencies:', f.shape)
print('max difference / THz:', np.max(np.abs(f-fc)))
np.savetxt('frequencies.csv', np.c_[q, f], delimiter=',',
           header='qx,qy,qz,f1,f2,f3,f4,f5,f6', comments='')
PYCODE

本例q点数组为(6,3)，频率为(6,6)。实测C/Rust频率最大差约3.24×10⁻¹⁰ THz；完整64点网格的差也约3.24×10⁻¹⁰ THz。不同输出格式的小数位数会让末位略有差别，不能据此宣称某个后端提高了材料计算精度。

左图是六个独立q点的六支频率，用线辅助辨认分支，不是高对称能带路径。右图逐个网格q点比较两后端，原始CSV和绘图脚本在资料包中。

5. 为什么很小的Γ点负频要看全精度矩阵#

本例Γ点三支声学模约为−1.0×10⁻⁴ THz，即约−0.00334 cm⁻¹。保留它们是为了展示原始数值，不代表判定硅存在真实动力学不稳定性。物理判断还要检查超胞、电子精度、声学和规则、网格与适当的对称化处理。

原生YAML动力学矩阵显示到10位小数时，极小的本征值可能在打印中被舍入掉；用这张文本矩阵重新对角化会得到零附近的数，与原生全精度频率相差约10⁻⁴ THz。这次差异来自矩阵输出精度，不能直接认为默认后端自动删除了负频。

使用HDF5的dynamical_matrix保留全精度，再独立计算本征频率即可比较。资料包实际操作：

with h5py.File('qpoints-Rust.hdf5') as h:
    dm = h['dynamical_matrix'][:]
    native = h['frequency'][:]
eigenvalue = np.linalg.eigvalsh(dm)
# 当前安装版本、当前VASP单位体系的换算系数；其他接口不能照抄。
factor = 15.633302300230191
independent = np.sign(eigenvalue) * np.sqrt(np.abs(eigenvalue)) * factor
print(np.max(np.abs(independent-native)))

换算系数应从当前Phonopy对象的p.unit_conversion_factor读取；原始数据文件中的旧版本显示值可能少几位小数。不同计算软件的力与长度单位不相同，不能一律套用VASP系数。

6. 什么时候改参数，什么时候保留原始结果#

• 仅比较后端：结构、质量、力常数、NAC、网格、ASR等全部保持一致。
• 关注很低的频率：保留HDF5矩阵与原始频率，别从低精度文本反推极小本征值。
• 做热性质：先检查低频与网格收敛，再按具体物理需求处理声学模；不要为让图看起来平滑直接裁负值。
• 求力常数：symfc/ALM是另一层求解方法，不能用切换C/Rust代替更好的位移—力数据。

本例所有数据、转换和绘图均为CPU轻量后处理；没有GPU加速比或CPU节点性能结论。官方接口单位和设置标签可进一步查询。