pymatgen：声子热容、熵和自由能实操
从Si真实力常数模型的DOS积分热力学函数，区分每摩尔原胞与每摩尔化学式。
阅读本篇目录
先按Phonopy数据互通准备total_dos.dat。这里计算的是固定谐力常数的振动贡献；没有加入静态电子总能、热膨胀、非谐效应或实验定压条件。晶体稳定性也不能仅从几个正频率点判断。
unzip pymatgen-phonon-wave665-*.zip
cd pymatgen-phonon
bash results.sh
bash plot.sh
reference-data/保存实际计算的数值，figures/reference.png为重绘图；读取和绘图不会启动VASP。若要运行本篇全部后处理，普通账号提交sbatch job.slurm。脚本申请cpu-8173m独占整节点1任务×56核，小例内部单线程，计费仍按整节点；不要在管理节点启动计算。
1. 低频端点先说明#
本例Gaussian展宽使0 THz处仍有0.01385216 states/(THz·原胞)的小密度。谐声子的熵和自由能权重在零频处有对数奇异；不要无说明地把零频代入对数，也不要把展宽尾部误认成扩散模。
本篇明确采用正频率0.01–20 THz作温度积分，原始0点和完整DOS仍保留。正频区积分5.99899479模/原胞，未重新缩放为6。这是有限网格和展宽下的教学处理；低温精确研究应改善低频采样并检查截断与展宽。
from pymatgen.io.phonopy import get_ph_dos
from pymatgen.phonon.dos import PhononDos
raw = get_ph_dos('total_dos.dat')
positive = raw.frequencies > 0
dos = PhononDos(raw.frequencies[positive], raw.densities[positive])
T = 300
print(dos.cv(T))
print(dos.entropy(T))
print(dos.internal_energy(T))
print(dos.helmholtz_free_energy(T))
2. 默认单位是每摩尔原胞#
| T / K | Cv / J·K⁻¹·mol-cell⁻¹ | U / kJ·mol-cell⁻¹ | F / kJ·mol-cell⁻¹ | S / J·K⁻¹·mol-cell⁻¹ |
|---|---|---|---|---|
| 50 | 6.079744 | 11.771268 | 11.641631 | 2.592734 |
| 100 | 15.813663 | 12.323336 | 11.337251 | 9.860852 |
| 300 | 39.876038 | 18.398005 | 6.169937 | 40.760229 |
| 600 | 46.975492 | 31.732229 | −11.044069 | 71.293830 |
| 1000 | 48.797110 | 50.973921 | −44.854194 | 95.828114 |
U与F包含该谐模型的零点振动贡献，所以低温U并不趋于零。上述函数返回J而不是kJ，表里U和F除以1000展示；F=U−TS保持同一归一化才成立。
3. 两原子原胞和一个Si化学式#
from pymatgen.io.phonopy import get_ph_bs_symm_line
structure = get_ph_bs_symm_line('band.yaml').structure
print('per mole-cell:', dos.cv(300))
print('per mole Si:', dos.cv(300, structure=structure))
本例原胞含2个Si化学式，提供结构后自动除以2；300 K得到19.93801920 J/(K·mol-Si)。若把39.876直接和以mol-Si报告的实验热容比较，会多出一倍。对于化合物，按其约化化学式数处理，不能总是除以原子数。
4. 高温极限用于核对归一化#
本例额外在10000 K检查形式上的高温极限，Cv=49.86739339 J/(K·mol-cell)，接近两原子6R。这个温度只是谐振子积分的数值检查，不是Si仍为稳定晶体的物理预测。
head reference-data/Si-pymatgen-phonon-thermal.csv
自己的材料要同时核对每胞声子模数、频率单位、负频与零模、归一化、采样和适用温度。这里没有用DOS拟合生成新的DFT数据。