pymatgen：E(V)到Debye热力学

已有一组能量—体积E(V)数据，想先估计温度变化下的体积和振动自由能，可以尝试各向同性准谐Debye近似。本例将鸢算科技集群实际计算的5个Si静态DFT点交给pymatgen，不重跑DFT。

这是一条使用E(V)、原子质量和假定泊松比的Debye近似，不是对每个体积计算声子后的完整准谐声子工作流，也不包含显式非谐MD统计。

下载实际E(V)、结构、完整脚本、热力学CSV与绘图脚本。

unzip pymatgen-debye-wave708-*.zip
cd pymatgen-debye
bash results.sh
bash plot.sh

1. 先统一原子数和单位#

输入为8原子Si常规胞，原E(V)表按每原子输出：

volume_A3_atom,E0_eV_atom
1.830124999999999957e+01,-5.349305151250000279e+00
1.954906250000000156e+01,-5.413179908750000102e+00
2.079687500000000000e+01,-5.423328306250000175e+00
2.204468749999999488e+01,-5.394624311250000304e+00
2.329250000000000398e+01,-5.338050218750000298e+00

第一列ų/atom，第二列eV/atom。传入8原子Structure时，两列都乘8，得到ų/cell和eV/cell；不能把每原子能量与整胞体积混在一起。

cat inputs/Si-EV-per-atom.csv
cat inputs/POSCAR-Si8

这些5个点来自同方法固定体积Si计算，覆盖相对参考体积88%、94%、100%、106%、112%。先前E(V)教学拟合见状态方程实操。

2. 设置Debye模型#

环境为pymatgen2026.9.24、Python3.11、NumPy/SciPy和Matplotlib，不需要重新编译VASP。使用当前类名QuasiHarmonicDebyeApprox：

import numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.quasiharmonic import QuasiHarmonicDebyeApprox

s = Structure.from_file('inputs/POSCAR-Si8')
ev = np.loadtxt('inputs/Si-EV-per-atom.csv', delimiter=',', skiprows=1)
n = len(s)
q = QuasiHarmonicDebyeApprox(
    energies=(ev[:, 1] * n).tolist(),
    volumes=(ev[:, 0] * n).tolist(),
    structure=s,
    t_min=100, t_max=500, t_step=100,
    eos='birch_murnaghan', pressure=0.0,
    poisson=0.25)
for T, V, G in zip(q.temperatures, q.optimum_volumes, q.gibbs_free_energy):
    print(T, V/n, G/n)

温度单位K，压力GPa，poisson=0.25是本例假定的泊松比，不是从这5个能量点自动拟合出的结果。有可靠弹性数据时，应使用与目标状态相符的值并检查敏感性。

3. 实际执行并查看输出#

在新的解压目录运行：

module load pymatgen/2026.9.24
bash run.sh
cat analysis/Si-Debye-temperature-volume-free-energy.csv

主要结果为：

这里振动自由能按8原子整胞计；换成每原子要除以8。它包含模型的零点振动贡献，所以低温下不必为零。

在零压下，通过E(V)+F_vib(V,T)寻找最小值。实际拟合体模量约为88.5651 GPa。结果是这个有限体积区间及Debye假设的预测，不作为已收敛Si热物性常数。

4. 自由能曲线如何画#

bash plot.sh
ls figures

三条曲线使用相同的静态能量参照，没有逐条平移最低点。右图放大最小值附近，便于观察温度引起的体积变化。

资料包保留直接连续最小化与原生拟合最小值的比较；本例5个温度的G差低于3.4E-7 eV/cell，体积差低于0.0001ų/cell。这些只说明本组数值处理一致，不消除有限E(V)数据和Debye物理近似的误差。

5. F、U、S分别是什么#

T = 300
V = q.optimum_volumes[2]
F = q.vibrational_free_energy(T, V)
U = q.vibrational_internal_energy(T, V)
S = (U - F) / T
print(F, U, S)

F和U为eV/cell，S为eV/(K·cell)。使用当前温度对应的体积，不要误取其他温度的索引。总G还含静态E(V)及有压时的PV项，不能把振动F单独当作整个晶体总自由能。

想得到更可靠的热膨胀、相间自由能差或精细相变温度，应增加体积点和收敛检查，结合真实声子或其他适用方法。Debye模型不提供声子色散、特定软模或完整各向异性响应；相关操作见声子与热力学。

接口与模型范围见pymatgen准谐模块。