Phonopy 原生QHA热膨胀实操

下载位移准备、真实力数据与QHA资料包。在已分配CPU资源内，用普通账号解压到新目录运行：

unzip phonopy-preparation-qha-wave242-*.zip -d phonopy-qha-practice
cd phonopy-qha-practice
bash run.sh
ls preparation thermal figures

脚本使用1个物理CPU，无需GPU。先实际生成位移和晶格采样，再从已提供的真实数值力与能量计算热膨胀，最后生成四张高清图；不会提交新的VASP计算。输出目录已存在时请另解压一份。

Phonopy为4.8.3，Python3.11.16、C/OpenMP由GCC13.3编译、默认phonors0.6.0后端，symfc为1.7.3。源数据来自本集群已经完成的VASP6.6.1 GPU计算：NVHPC25.3、CUDA12.8、HPC-X2.22.1、MKL2023.2。电子计算与本包的CPU后处理是两个阶段。

1. 先明确要准备哪些计算#

QHA把不同晶格下的静态能量与声子振动自由能相加，在每个温度寻找自由能最低的位置。对本例立方Si，只有一个独立晶格参数a；用各向异性QHA程序演示的是立方一维特例，不是一个多轴独立变化的材料实例。

完整流程：

1. 确认平衡参考结构和电子精度。
2. 选择五个体积，每个体积计算一次静态总能量。
3. 每个体积生成位移结构并计算原子力，本例共10个受力计算。
4. 收集位移、力和对应静态能量，统一原胞归一化。
5. 生成QHA数值数据集，设置温度和声子网格。
6. 提取a(T)、体积和热膨胀系数，比较网格和拟合设置。

全部15份电子计算是本集群已有真实Si结果。本篇新运行的是收集与后处理；包内保留原始INCAR/POSCAR/KPOINTS和数值力/能量。较早的Vinet体积路径QHA是拟合F(V)；这里展示直接拟合F(a)的原生程序，二者拟合形式不同。

2. 静态能量与位移受力分别设参数#

源计算均ENCUT=550 eV、EDIFF=1E-6、PBE PAW64。静态原胞2Si使用9×9×9电子k网格：

Gamma mesh
0
Gamma
9 9 9
0 0 0

其INCAR的关键部分：

ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
IBRION = -1
NSW = 0
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
ISYM = 2
LWAVE = .FALSE.
LCHARG = .FALSE.

16Si位移超胞使用4×4×4电子k网格、ISYM=0，完整输入见位移与力收集和附件inputs/DFT-inputs/。不同大小晶胞不能机械套用同一个k网格数；本例精度是操作参考，自己的材料仍需检查。

五个参考体积比例是0.96、0.98、1.00、1.02、1.04。静态能量读VASP6的energy(sigma→0)，按每个2Si原胞计；资料包STATIC-ORACLE.json记录实际数值与来源哈希。

3. 如何从自己的真实输出生成数据集#

先按每个晶格点收集力，保留成phonopy_params.yaml。示例目录按01至05编号，使静态输出与声子数据按同一顺序排列：

reference.yaml
static/01.xml ... static/05.xml
phonons/01/phonopy_params.yaml
phonons/02/phonopy_params.yaml
...
phonons/05/phonopy_params.yaml

reference.yaml与声子数据需使用符合晶轴约定的常规胞；原胞和超胞对应关系见晶格采样与晶胞对应。有了自己的真实结果后，命令形式为：

phonopy-anisotropic-qha-dataset reference.yaml \
  --static static/{01..05}.xml \
  --phonon phonons/{01..05}/phonopy_params.yaml \
  --no-electronic -o aniso_qha_dataset.hdf5

Bash展开后的两组列表一一对应。当前版本会检查静态晶胞和声子晶胞是否相配；目录排序错误不能靠后面拟合来修复。若只提供静态能量、不提供声子数据，生成的数据集也不会凭空包含声子自由能。

本集群已用本例真实完整输出执行过数据集生成。下载包提供其数值HDF5及带力YAML，不分发完整原始电子计算输出；用户可直接进入下一步复现热膨胀，也可用自己的原始输出重建数据集。

4. 先用一组声子网格运行#

在初始解压目录创建新的手动练习目录：

mkdir manual-qha
cp inputs/aniso_qha_dataset.hdf5 manual-qha/
cd manual-qha
phonopy-anisotropic-qha aniso_qha_dataset.hdf5 \
  --mesh 24 --tmax 500 --dt 100 --polynomial-degree 2 \
  --no-electronic --smooth-lattice none
ls *temperature.dat axial_thermal_expansion.dat

这里的--mesh 24是网格密度长度参数，本例实际生成8×8×8声子q网格；不是24×24×24网格，也不是VASP的电子k网格。第二组密度长度48实际是15×15×15声子网格。

--polynomial-degree 2用二次自由能面拟合五个点，--no-electronic不加电子温度自由能，--smooth-lattice none不做晶格温度曲线平滑。原生热统计排除Γ点三个声学模；应先检查实际声子，再选择低频处理。

当前程序对本例单原子位移的type-1数据，先传统有限差分求力常数，再用symfc投影对称化。--fc-calculator选择拟合器主要作用于type-2密集/随机位移数据，不能把命令里显示symfc就理解成所有力常数都由随机位移拟合得到。这也不同于原始力常数C/Rust比较中不加对称化的处理。

5. 看实际结果#

head -10 lattice_parameters-temperature.dat
head -10 volume-temperature.dat
head -10 axial_thermal_expansion.dat

输入温度0–500 K，每100 K一个点；当前中央差分流程不返回最高温度点，因此最终表格是0–400 K。不是500 K计算被遗漏。示例0 K热膨胀系数行按程序边界约定写0，温度导数应结合有限差分和温度步长理解。

两组a的最大差约9.98×10⁻⁶ Å。本例有低温收缩和更高温度的膨胀趋势；它用于说明提取流程，不作为收敛材料参数或实验定量比较。

曲线对每个温度减去该温度采样点的最低值，纵轴为每2Si原胞的meV。拟合最低点可以落在两个采样点之间，所以曲线局部低于零并非负的绝对自由能错误。

a,b,c是常规晶格长度；本例volume-temperature.dat为2Si原胞体积，满足V=a³/4。热膨胀表列出各轴系数和轴向系数之和，不应再把它们乘上原子数。

6. 自己算材料时优先改什么#

• 体积/晶格采样范围应包含各温度的极小值；极小值落在范围外时，先补真实采样点。
• 检查电子k网格、超胞大小、位移幅度与声子q网格，不能只比较两条画得接近的曲线。
• 温度步长会影响热膨胀导数，本例100 K是简明操作演示，可在自己的数据上加密。
• 多独立晶格参数需要多维采样点；不要把立方Si的一维五点数据套到六方或正交材料。
• 对复杂结构，需要在每个固定晶格点处理内部原子位置，然后再做位移求力；不要让不同体积的声子数据来自不同相或不同原子顺序。

原始数值、两组计算、独立拟合及绘图脚本均在资料包中。更多配置含义见Phonopy官方设置说明；本篇以实际4.8.3命令为准。