GPAW有限位移声子

下载完整Python输入、原始结果、CSV和绘图脚本。以下短命令在解压目录、已分配的CPU资源中执行，输出写入 work/。本模块为CPU构建，加载GPU资源不会自动启用GPU计算。

1 准备并运行13个力条件#

bash run.sh phonons

该入口在已分配CPU资源中执行；串行会依次计算全部位移，缓存位于 work/response-phonons-three/forces/。核心代码：

from ase.build import bulk
from ase.phonons import Phonons
from gpaw import GPAW, PW, FermiDirac
atoms = bulk('Si', 'diamond', a=5.43)
calc = GPAW(mode=PW(300), xc='PBE', kpts=(1,1,1), nbands=120,
    occupations=FermiDirac(0.01),
    convergence={'energy':1e-7,'density':1e-6,'eigenstates':1e-9},
    symmetry={'point_group':False}, txt='phonon.txt')
p = Phonons(atoms, calc, supercell=(3,3,3), delta=0.01,
            name='forces')
p.run()
p.read(acoustic=True)
path = atoms.cell.bandpath('GXWKGL', npoints=100)
bands = p.get_band_structure(path)

Si原胞2原子，3×3×3超胞54原子。对两个原胞原子的三个方向做正负位移，共12个条件，另有未位移参照，合计13个；不是对54个原子全部做324个独立位移。

2 检查力与缓存#

ls work/response-phonons-three/forces/cache.*.json

本次参考结构最大原子力约1.46×10⁻⁵ eV/Å。每个位移都进行实际DFT力计算，不能用弹簧常数示意值代替这些力。

位移delta=0.01 Å；PW300 eV、超胞Γ点电子采样是此例设置。关闭点群后，任意偏移k网格不一定保留晶体所需对称性，检查理想结构的残余力很有帮助。

3 从力常数到色散#

路径为Γ–X–W–K–Γ–L，纵轴为声子能量meV。最高约66.89 meV；Γ点声学模约10⁻⁹ eV量级的带符号数值零，不解释为显著物理虚频。

此前较小2×2×2超胞会产生明显的数值虚支。加密电子网格减少了残余力，继续扩大到奇数3×3×3超胞后得到上图；因此虚频需要先检查电子精度和有限超胞，而不是直接判断材料不稳定。

4 自己的体系如何继续#

先优化结构，再分别检查位移大小、超胞、电子k网格和PW截断。acoustic=True施加声学和约束，不会自动修复所有不收敛的力常数。极性材料的LO–TO分裂还需要Born有效电荷与介电张量，不能从这个非极性Si例直接得到。

更长程、高阶和热导流程可继续连接Phonopy或phono3py，但应先确保力常数格式、质量和单位匹配。