ASE原生有限位移声子与晶体热力学

下载全部输入、可运行脚本、实际数据和四张高清图。在已分配CPU资源内，使用新解压目录执行 bash run.sh；ASE3.29.0、单进程单线程，生成 work/。本包不调用GPU或外部DFT程序。

1 先看完整流程#

建立参考晶胞 → 选择计算器 → 超胞有限位移求力 → 读取力常数并施加声学求和规则 → 声子带和DOS → 晶体热力学。 本例使用ASE原生Phonons接口，不调用Phonopy。

module load ase/3.29.0
bash run.sh
cat work/phonon-output/RESULTS.json
head work/phonon-output/dos.csv

Cu用内置EMT经验势，晶格常数3.6Å，原胞1原子，3×3×3超胞共27原子。位移幅度0.01Å。结果用于学习流程，EMT数值不当作DFT或实验值。

2 准备结构并计算位移力#

完整程序是 phonons_thermo.py，关键部分：

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.phonons import Phonons

atoms = bulk('Cu', 'fcc', a=3.6)
phonons = Phonons(atoms, EMT(), supercell=(3,3,3),
                  delta=0.01, name='phonon')
phonons.run()
phonons.read(acoustic=True)

delta单位Å，acoustic=True施加平移不变性约束。缓存保存各位移实际计算的力，后续可读取同一组缓存；结构、势函数或位移设置改变后，应使用新目录重新计算，不能沿用不匹配的缓存。

本例参考原胞只有一个原子，三个方向的正负位移加参考共七组力调用，合起来是一个声子工作流，不是七种材料。

3 计算声子带和DOS#

path = atoms.cell.bandpath('GXWKGLUWLK,UX', npoints=100)
band_energies = phonons.band_structure(path.kpts, verbose=False)
raw = phonons.get_dos(kpts=(8,8,8), verbose=False)
dos = raw.sample_grid(npts=800, width=0.0004)
energy = dos.get_energies()
weight = dos.get_weights()

ASE这里的声子能量单位为eV；图中乘1000显示meV。DOS展宽0.0004eV即0.4meV，不能把它当作电子温度。路径中的逗号表示不连续线段，绘图应分段连接。

Γ点三个声学模的能量绝对值低于5×10⁻¹⁰eV；沿路径最高声子能量约 33.0817meV。DOS积分为 2.9999708态/原胞，与一个原子三个振动自由度相符。Γ点接近零并不能独自证明整条声子谱稳定，仍需检查整个q空间和数值设置。

4 从DOS得到晶体自由能#

from ase.thermochemistry import CrystalThermo
atoms.calc = EMT()
thermo = CrystalThermo(phonon_DOS=weight, phonon_energies=energy,
                      potentialenergy=atoms.get_potential_energy(),
                      formula_units=1)
print(thermo.get_helmholtz_energy(300, verbose=False))

这里DOS按一个Cu原胞归一，formula_units=1。零点能为 0.03247191eV/atom；温度相关U、S、F写入 crystal-thermo.csv，其中 F=U−TS。不要把总DOS、投影DOS和每原子归一混用，也不要在热力学公式中把meV当eV。

5 换材料时检查什么#

修改晶胞、计算器、超胞和位移；检查更大的超胞、更多q点及不同位移幅度。Si不适合直接使用这套Cu EMT参数，不能只把元素名改成Si；应选择对应DFT、机器学习势或其他可靠势函数。

本例没有完成超胞与q网格收敛，也没有给出非谐寿命和热导率；谐性声子DOS与晶体自由能并不包含这些量。气相分子的平移/转动自由能另见ASE分子热化学。接口参考：ASE Phonons。