Critic2 声子态密度与振动热力学
从真实力常数计算q网格、态密度、零点能和F/S/Cv，分清每原胞、每化学式单元与Gamma点近似。
阅读本篇目录
声子色散用于看路径上的频率；振动热力学需要对整个布里渊区采样。本节从同一组实际Si16 NEP拟合力常数开始，计算移位q网格、声子态密度和谐近似热力学表。
1. 准备与运行#
下载声子与振动热力学实例包，在已分配的CPU资源中解压：
unzip critic2-vibrations-wave145-*.zip -d critic2-vibrations
cd critic2-vibrations
bash run.sh list
本例使用 critic2/1.4（内部1.4.155，GCC/GFortran13.3、OpenMP、LAPACK3.12.1/BLAS）；对照与绘图使用 phono3py/4.8.2 的独立Python3.11.16环境，其中Phonopy为4.8.1。全部使用1个CPU线程，无需GPU。
先读力常数与声子色散，确认结构、力常数的原子序和eV/Ų单位。下面使用已有实际NEP/symfc结果，不重新进行DFT或求力。
bash run.sh 9
head -n 14 work/thermo6/thermo.dat
完整6³网格输入：
crystal ../../data/unit.vasp
vibrations calculator vasp
vibrations load_fc2 ../../data/FORCE_CONSTANTS 2 2 2
vibrations thermo mesh 6 6 6 shift .5 .5 .5 tmin 0 tmax 1000 tstep 50 cutoff 1 noxdebye file thermo.dat
MESH 6 6 6是216个q点的均匀网格；SHIFT .5 .5 .5指各方向 半个网格步长，坐标为 (i+0.5)/6，不是把每个q点在倒空间直接平移0.5。内部可以利用对称性减少实际对角化次数，但热力学仍按完整网格权重求和。
温度从0到1000K，步长50K。CUTOFF 1的单位为cm⁻¹，低于此频率的模式不纳入热力学求和；本例移位网格所有1296=216×6个模式均纳入。NOXDEBYE关闭额外Debye–Einstein拟合，使用直接谐振求和。
2. 读取零点能和F、S、Cv#
表头明确两种归一方式。第2–4列按 2原子的原胞，第5–7列按 Si化学式单元，即1个Si原子。两组差2倍；不是把整个16原子超胞当成输出归一单位。
6³移位网格实际结果：
| 温度 / K | Fvib / kJ·mol⁻¹原胞 | S / J·K⁻¹·mol⁻¹原胞 | Cv / J·K⁻¹·mol⁻¹原胞 |
|---|---|---|---|
| 0 | 11.4723648 | 0 | 0 |
| 300 | 5.7623454 | 41.8464534 | 40.1204460 |
| 1000 | −46.0901354 | 97.0450807 | 48.8323401 |
0K的Fvib等于零点能ZPE；负的高温Fvib包含熵贡献，并不代表结构不稳定。该表只有固定体积下的 振动贡献，没有电子总能、pV项和热膨胀。这里的Cv也不是恒压热容Cp。
grep 'ZPE' work/thermo6/thermo.dat
grep '300.000' work/thermo6/thermo.dat
3. 用4³网格检查采样影响#
bash run.sh 7
grep '300.000' work/thermo4/thermo.dat
实际300K结果：
| q网格 | Fvib / kJ·mol⁻¹原胞 | S / J·K⁻¹·mol⁻¹原胞 | Cv / J·K⁻¹·mol⁻¹原胞 |
|---|---|---|---|
| 4³移位 | 5.7648466 | 41.8381159 | 40.1204460 |
| 6³移位 | 5.7623454 | 41.8464534 | 40.1204460 |
这是网格敏感性的实操示例；正式计算需按目标温区和精度继续检查，不能仅凭两个网格就宣布热力学全部收敛。
4. 计算声子态密度#
bash run.sh 10
bash run.sh 11
head -n 6 work/dos-gaussian/dos.dat
head -n 6 work/dos-tetrahedron/dos.dat
高斯展宽完整输入：
crystal ../../data/unit.vasp
vibrations calculator vasp
vibrations load_fc2 ../../data/FORCE_CONSTANTS 2 2 2
vibrations thermo mesh 6 6 6 shift .5 .5 .5 temperatures 0 300 1000 noxdebye dosfile dos.dat sigma 8 npts 1001 file thermo.dat
四面体输入：
crystal ../../data/unit.vasp
vibrations calculator vasp
vibrations load_fc2 ../../data/FORCE_CONSTANTS 2 2 2
vibrations thermo mesh 6 6 6 shift .5 .5 .5 temperatures 0 300 1000 noxdebye dosfile dos.dat tetrahedron npts 1001 file thermo.dat
二者同用6³移位网格。SIGMA 8表示 8cm⁻¹高斯宽度；NPTS 1001表示DOS输出频率网格点数，不是q网格点数。高斯曲线按输出区间积分，并在约±5sigma处截尾；四面体利用q网格连接关系，不是给高斯换一个宽度。
实际DOS积分分别为 5.99999797 和 6.00000001 模/原胞，接近 3×2=6。本包保留原生数值，不把曲线再强制归一。两种绘图方式的原生F、S、Cv表相同，因为直接热力学求和使用同一组实际模式。
5. 不要用Γ点替代全布里渊区#
bash run.sh 13
grep '300.000' work/clear-reload/gamma-thermo.dat
本例Γ点单独求和在300K得到Fvib=7.9768829kJ/mol原胞、Cv=16.1365831J/K/mol原胞，而6³网格分别为5.7623454和40.1204460。Γ点缺少全布里渊区的声学支采样；即使程序正常输出，也不能当作完整晶体热力学。
完整重算与画图用 bash run-all.sh，只画随附数据用 bash analyze.sh。理论定义可参看Phonopy官方公式。
只分析附带结果时，检查表在 analysis-reference/CHECKS.json，图片在 figures-reference/；完整重算分别写入 analysis/ 和 figures/。
实际求力到力常数#
固定Si16真实QE单点求力、原生有限/随机配对FC2与Phonopy/Symfc对照。QE导出的FC2按Ry/bohr²理解；finite reader扣除单帧整体力漂移，但不自动扣未位移参考。随机正负顺序1/3、2/4，少量数据显式NORIDGE或RIDGE；自动CV需要更多配对数据。