Phonopy QHA状态方程与压力
用实际五体积Si数据运行原生phonopy-qha，比较三种EOS和外压，解释文件、单位与温度范围。
阅读本篇目录
下载原生QHA、随机位移与结构优化实例包，在普通账号已分配的CPU资源中运行：
unzip phonopy-qha-random-relax-wave247.zip -d phonopy-practical
cd phonopy-practical
bash setup.sh
bash run.sh
setup.sh创建本目录个人.env，需要联网下载固定版本的conda-forge Linux x86_64包；公共模块不变。环境为Phonopy4.8.3、Python3.11.18、NumPy2.4.6、pypolymlp0.22.1和symfc1.7.3，实际预编译C/C++运行库libgcc/libstdcxx16.2.0。绘图另用集群phonopy/4.8.3模块。所有计算使用1个CPU线程，不需要GPU，不启动DFT。
一键流程输出source/、analysis/和figures/；已有结果目录不覆盖，重新练习请换新的解压目录。以下手动命令都放在独立新目录，便于自己修改参数。
1 先认识两类输入#
cat inputs/e-v.dat
head -n 24 inputs/thermal-v100.yaml
e-v.dat每行是原胞体积ų和电子静态能eV；五个thermal-v*.yaml分别是相同体积的振动自由能、熵和热容。本例按体积从小到大排列：
39.2570229468125 -10.82886164
40.0748775915378 -10.84363673
40.8927322362631 -10.84830721
41.7105868809883 -10.84388942
42.5284415257136 -10.83130901
数据复用鸢算科技此前实际Si计算：VASP6.6.1、ENCUT550 eV、EDIFF1E-6，2Si静态电子k网格9³，16Si位移受力电子网格4³；声子热性质为16³ q网格，排除Γ点三个声学平移。五个体积为参考体积96%、98%、100%、102%、104%。
本包没有把新模型力代替这些DFT力；thermal-v*.yaml每摩尔指每摩尔2Si原胞。电子静态能和振动项必须采用相同原胞归一，不能将16Si超胞总能直接与2Si振动自由能相加。
2 用原生QHA命令拟合#
mkdir manual-qha
cd manual-qha
../.env/bin/phonopy-qha ../inputs/e-v.dat \
../inputs/thermal-v{96,98,100,102,104}.yaml \
--eos vinet --tmax 500
cd ..
花括号在Bash中展开为五个顺序确定的文件，不用无序的通配符。命令在每个温度拟合F(V,T)，寻找体积最小值，并写出体积、体膨胀、体模量及相关热性质文件。
cat source/qha-vinet/volume-temperature.dat
cat source/qha-vinet/bulk_modulus-temperature.dat
cat source/qha-vinet/thermal_expansion.dat
3 三种EOS与1 GPa压力的实际结果#
一键流程还运行--eos birch_murnaghan、--eos murnaghan，以及Vinet加--pressure 1。压力单位GPa，加入PV项后最小化相应势函数，不需要把压力数值写进e-v文件。
| 300 K条件 | 平衡体积/ų·(2Si)⁻¹ |
|---|---|
| Vinet，0 GPa | 41.07859106 |
| Birch–Murnaghan，0 GPa | 41.07864957 |
| Murnaghan，0 GPa | 41.07877825 |
| Vinet，1 GPa | 40.61113098 |
三条零压曲线很接近，是这组输入和体积区间的结果，不表示EOS选择总是无关紧要，也不替代电子、超胞和声子网格收敛。所有本例最小值都在五个原始体积范围内；如果最小值跑到采样边界外，应补算体积点，而不是外推后直接作为结果。
4 读文件时注意什么#
volume-temperature.dat是ų/2Si原胞；bulk_modulus-temperature.dat是GPa；thermal_expansion.dat是体膨胀系数K⁻¹。立方、各向同性条件下线膨胀系数约为体膨胀的1/3，不能直接把体膨胀曲线当作晶格常数的相对变化率。
这里输出温度0–500 K、步长50 K。日志中用于差分的额外550 K拟合点不意味着输出范围已扩展；看实际.dat文件。本例0 K膨胀系数0是程序边界赋值，不是额外计算出的收敛材料结论。低温体积小幅下降同样只对应当前有限设置。
需要重新准备位移和受力时看硅有限位移；从常规晶胞和三次晶格拟合出发时看原生QHA热膨胀。命令说明见Phonopy官方QHA文档。