Phonopy 模式Grüneisen参数实操

下载模式Grüneisen与原生绘图资料包，在已分配CPU内用普通账号解压到新目录：

unzip phonopy-gruneisen-cli-wave243-*.zip -d phonopy-mode-practice
cd phonopy-mode-practice
bash run.sh
ls work figures

每次1个CPU，无需GPU。脚本从真实数值位移和力求力常数，再运行12条原生命令，输出六个PDF、三张高清图与CSV。它不会提交新的VASP计算。已存在输出目录时，另解压一份继续练习。

环境是Phonopy4.8.3、Python3.11.16、NumPy2.3.5、GCC13.3 C/OpenMP、phonors0.6.0及symfc1.7.3。源数据来自本集群实际VASP6.6.1 GPU计算，构建环境NVHPC25.3、CUDA12.8、HPC-X2.22.1与MKL2023.2；后处理不改变电子计算来源。

1. 这个参数能看什么#

模式Grüneisen参数描述某支声子频率随体积变化的相对敏感程度，常用写法为γ=−d ln ω/d ln V。正值表示该模式通常随膨胀软化，负值表示相反趋势。它是模式量；热膨胀还需要整体自由能与弹性响应，不能把几个γ直接平均就当作热膨胀系数。

先在参考体积求声子，再在略大和略小体积各求一组声子，最后在相同分数q点比较动力学矩阵。本例使用参考体积的98%、100%、102%，两原子Si原胞、16原子受力超胞，共复用六个既有VASP位移受力结果，没有新增电子计算。

源输入ENCUT=550 eV、EDIFF=1E-6、位移幅度0.02 Å，16原子电子k网格4×4×4。完整INCAR/POSCAR/KPOINTS见位移与力收集及QHA数据资料。本篇输入是数值FORCE_SETS和位移YAML；力常数由传统有限差分求得，并对三个体积统一进行symfc投影对称化。

2. 先确认三个目录#

执行过bash run.sh后：

ls work/v100/POSCAR work/v100/FORCE_CONSTANTS
ls work/v102/POSCAR work/v102/FORCE_CONSTANTS
ls work/v98/POSCAR work/v98/FORCE_CONSTANTS

每个目录有该体积的原胞POSCAR、对应力常数和完整Phonopy数值文件。三个计算必须具有一致的元素顺序、原胞定义、超胞矩阵和计算精度。体积百分比不是晶格长度百分比；对立方材料，晶格长度按体积比的立方根缩放。

3. 用均匀网格计算#

以下在初始解压目录开始，另建一个手动练习目录：

module load phonopy/4.8.3
mkdir manual-gruneisen
cd manual-gruneisen
phonopy-gruneisen ../work/v100 ../work/v102 ../work/v98 \
  --dim '2 2 2' --readfc -c POSCAR \
  --mesh '4 4 4' --nomeshsym --gc
ls gruneisen_mesh.yaml

目录顺序是参考、大体积、小体积，不是按目录名从小到大。--readfc读取各目录的FORCE_CONSTANTS，-c POSCAR指定各体积目录内的结构文件。--gc使用Γ中心网格。

当前4.8.3实际输出名为gruneisen_mesh.yaml，不是部分示例中的通用gruneisen.yaml。

为什么4×4×4只写了36个q点#

这个专用命令当前仍使用时间反演约化，因此有36个代表点。YAML的multiplicity给出每个代表点覆盖的网格点数，总和为64；不要把36行误认为网格缺了点，也不要把另一个工具的weight字段名直接套进来。

python - <<'PYCODE'
import yaml
with open('gruneisen_mesh.yaml') as f:
    d = yaml.safe_load(f)
rows = d['phonon']
print('representatives:', len(rows))
print('total mesh weight:', sum(r['multiplicity'] for r in rows))
r = rows[1]
print('q:', r['q-position'])
print('frequency / THz:', r['band'][0]['frequency'])
print('gamma:', r['band'][0]['gruneisen'])
PYCODE

本例输出代表点36、总权重64；q=(0.25,0,0)的低频支约2.25354 THz、γ约−1.39903。Γ点光学模式约15.07351 THz、γ约0.98117。它们是本组小超胞和体积差下的实际结果，不是通用Si材料参数。

4. 直接用原生工具画图#

phonopy-gruneisenplot gruneisen_mesh.yaml \
  --cutoff .2 -s -o mode-gruneisen.pdf

这里--cutoff .2只用于这张原生网格图的低频显示；它不修改原始力常数或YAML。原始CSV保留全部216个代表点—分支条目，包括Γ附近微小频率的结果。

Γ点声学频率接近零时，相对体积导数会变得很敏感；本例可见约10⁻⁷ THz的原生微小频率，不据此报告真实不稳定性。图只显示频率大于0.2 THz的模式，以便看清有限频率分支。检查实际物理结果时，仍要补足超胞、体积差、电子精度和近Γ点采样的收敛分析。

5. 沿声子路径计算#

回到初始解压目录，在另一个新目录运行：

mkdir manual-gruneisen-band
cd manual-gruneisen-band
phonopy-gruneisen ../work/v100 ../work/v102 ../work/v98 \
  --dim '2 2 2' --readfc -c POSCAR \
  --band '0 0 0 .5 0 .5 .5 .5 .5' --band_points 21
phonopy-gruneisenplot gruneisen_band.yaml \
  -s -o band-gruneisen.pdf

该专用工具使用--band_points下划线形式，主phonopy命令使用--band-points；按各自--help核对。输出是gruneisen_band.yaml。本例路径为Γ→X→L，每段21点；路径表示要适合自己的晶格。

声子支经过简并或交叉时，不要仅把三个体积排序后的第n条频率相减；本例通过动力学矩阵导数与简并子空间进行独立核对。原生路径绘图对Γ附近的处理也可能与直接逐行画原始YAML不同，读数优先保留原始文件并说明低频处理。

6. 最后得到哪些文件#

• work/gruneisen-mesh/gruneisen_mesh.yaml：代表q点、重数、频率与γ。
• work/gruneisen-mesh/mode-gruneisen.csv：所有代表点和六个分支，频率THz、γ无量纲。
• work/gruneisen-band/gruneisen_band.yaml：路径数据。
• 对应目录内的原生PDF与figures/mode-gruneisen.png：可直接查看。

在自己的材料上，先为三个体积计算真正的力，再使用这条流程；不要把结构缩放本身当成新的受力数据。理论与接口参考Phonopy模式Grüneisen说明，文件名和命令以上述实测4.8.3为准。