QE：从声子DOS计算振动自由能

下载已有计算数据、完整输入、实际后处理输出、CSV和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-phonon-band-neb-tools
module load ase/3.29.0 qe/7.6-8173m
bash read.sh
bash plot.sh

本例使用鸢算科技集群的 QE 7.6 CPU 构建，Intel oneAPI 2024.2、Intel MPI 2021.13 与 MKL 环境。这里直接读取已实际完成的电子结构、声子或 NEB 数据，只运行后处理，不重复 DFT。

1. 先明确自由能包含什么#

fqha.x 从声子态密度计算固定体积下的谐振动 Helmholtz 自由能，包含零点能。本例只提供振动部分，没有加电子总能量、电子热激发或 pV，也没有扫描多个体积，因此不是完整准谐近似热膨胀流程或 Gibbs 自由能。

源数据是两原子 GaN 原胞实际计算的力常数，电子 k 网格4×4×4、声子 q 网格2×2×2；这里采用晶体声学和规则，在8×8×8网格插值 DOS，并关闭长程LO–TO项。网格较粗，适合说明操作，不作为收敛的 GaN 热力学数据。

2. 从力常数生成DOS#

完整 inputs/dos.in：

&input
 asr='crystal',flfrc='gan.fc',dos=.true.,
 nk1=8,nk2=8,nk3=8,deltaE=2.0,
 fldos='gan-dos.dat',loto_disable=.true.
/

在自己的 CPU 资源中执行：

mkdir thermo-work
cp source-data/gan.fc thermo-work/
cd thermo-work
matdyn.x < ../inputs/dos.in > dos.out
head gan-dos.dat
cd ..

前两列是频率/cm⁻¹和声子 DOS/每cm⁻¹的态数；附加列是投影信息。fqha.x 需要没有注释头的两列数值文件，转换时保留所有原始频率和权重：

python3 - <<'PYCODE'
import numpy as np
a = np.loadtxt('thermo-work/gan-dos.dat')
np.savetxt('thermo-work/gan-dos-numeric.dat', a[:, :2], fmt='%.12e')
PYCODE

不要将电子 DOS 的 eV 横轴文件直接代入声子热力学公式，也不要用相对强度谱替代按振动模态归一化的 DOS。

3. 输入温度并计算#

inputs/fqha.in：

gan-dos-numeric.dat
GaN-vibrational-free-energy.dat
0
50
100
200
300
500
800
1000
-1

前两行是输入DOS和输出文件；随后每行一个温度。末尾 −1 是结束输入的控制值，不是要计算负温度。

cd thermo-work
fqha.x < ../inputs/fqha.in > fqha.out
cat fqha.out
cat GaN-vibrational-free-energy.dat
cd ..

4. 实际结果与DOS积分#

DOS积分为 5.989705，两原子原胞应对应6个振动模态。本例保留离散积分误差，没有强制归一化到6。最小频率约 −1.14×10⁻⁵ cm⁻¹，对应DOS为0；此极小数值在本版本程序中按0处理，不表示处理了有实际权重的虚频。若体系存在真实虚频，不能靠删点或取绝对值获得稳定相自由能。

T=0 的值就是本次DOS积分得到的零点能。F变负不等于体系失稳；这里比较的是振动能与熵项组合后的量。

5. 导出与正式计算#

cat analysis-read/vibrational-thermodynamics.csv
cat analysis-read/SUMMARY.json

CSV中的U、S、Cv由同一谐振子分配函数后处理得到；fqha.x 原生输出是温度和F。本版本使用历史物理常数，CSV保持该约定，避免混用转换常数造成细微差别。

完整资料包也能一次重做三个后处理例子，见README；正式热力学研究应先检查结构稳定性、q网格、DOS积分和目标温区。若做准谐近似，需要多个体积下独立的能量与声子数据。