QE：状态方程模型与单位对照

下载五份原始SCF日志、拟合输入与实际输出、CSV 和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-state-equation
module load ase/3.29.0 qe/7.6-8173m
bash read.sh
bash plot.sh

使用 QE 7.6 的 ev.x；CPU 构建为 Intel oneAPI 2024.2 / Intel MPI 2021.13 与 MKL。数据来自此前实际算出的五组 Si₂ 总能和原胞体积，本篇只做状态方程拟合，不重新跑 SCF。后处理使用 NumPy 和 Matplotlib。

1. 先把模型差异限定在同一份数据#

所有模型都使用同样的五个 Si₂ 原胞体积和总能，只改变 ev 输入第三行：

mkdir birch1-new birch2-new
cp inputs/volume-energy.dat birch1-new/
cp inputs/volume-energy.dat birch2-new/
(cd birch1-new && ev.x < ../inputs/birch1.in > native.out)
(cd birch2-new && ev.x < ../inputs/birch2.in > native.out)
head -n 5 birch1-new/fit.out
head -n 5 birch2-new/fit.out

输入编号 1 和 2 分别对应程序提示中的 birch1 与 birch2。该版本不同输出位置的阶数文字标签并不完全相同，复现时应同时记录输入编号、公式约定和参数，不要只按一个阶数字样判断模型。

实际对照：

本例三种模型在采样区间附近相近；这不是对任意材料、任意体积区间都成立的保证。

2. 高阶模型残差很小也可能不可靠#

本例只有五个点，带更多参数的模型可以把残差压到接近零，但 B′、B″ 等参数可能仍不稳定。不要用“拟合线穿过所有点”替代更多体积采样、电子参数收敛和不确定度评估。

本包同时检查了压力是能量对体积的负导数、体模量对应平衡点能量曲率，帮助识别只看残差时容易漏掉的单位或模型问题。

3. FCC也可以输入常规胞晶格常数#

fcc 模式第一列读取 FCC 常规立方胞的 a，程序按 V=a³/4 转为原胞体积。不是把第一列原胞矢量长度直接立方，也不是把 V 当作 a。

本包从同一组原始体积换算出 fccAng-energy.dat，晶格常数单位为 Å；完整输入是：

mkdir fcc-unit-new
cd fcc-unit-new
cp ../inputs/fccAng-energy.dat .
cat > ev.in <<'EOF'
Ang
fcc
4
fccAng-energy.dat
fit.out
EOF
ev.x < ev.in > native.out
head -n 7 fit.out
cd ..

只改变第一列的几何输入约定，能量列仍然是 Ry。

此条件恢复约 41.3081 ų 的原胞体积和 88.1511 GPa 的体模量，与直接体积输入一致。用于换算的原胞和常规胞关系可先按晶胞转换工具核对。

4. 当前Keane输出的使用范围#

本包还实际运行了编号 3 的原生 Keane 条件，原始文件保留在 reference-results/keane/。它的五点能量拟合残差很小，但该版报告的 B₀ 和压力与其能量曲线的曲率、导数不一致。

这组结果中，原输出 B₀ 约 125.03 GPa，同一能量曲线的平衡曲率却约 88.10 GPa；压力最大差约 1.73 GPa。因此本篇不把它的参数作为可直接报告的物性，默认快捷计算采用前面三种已对照的模型。

cat analysis-read/model-boundary.json
head -n 8 reference-results/keane/fit.out

这是当前 QE 7.6 程序输出的具体范围，不表示所有 Keane 理论公式都不可用。换版本时应重新检查对应实现。当前源码可见 QE 7.6 ev。

5. 读取新拟合#

bash run.sh
bash read.sh work-new analysis-new
bash plot.sh analysis-new
cat analysis-new/fit-summary.csv

fit-summary.csv 同时保留原生 B₀ 和由能量曲率重建的值；所有曲线从指定结果读取。资料包保留原始输出，便于核对具体列与单位。