📚 课程目录第4章声子色散与态密度计算

4. 声子色散与态密度计算

声子色散曲线和声子态密度（Phonon Density of States, PDOS）是晶格动力学计算中最常交付的两类结果。色散曲线回答“哪些波矢下存在哪些声子模式、频率如何随 q 变化”，态密度回答“在某一频率区间内有多少声子模式可供占据”。两者互为补充：色散是 q 空间中的带状结构，态密度是频率轴上的模式计数。本节围绕高对称路径选取、色散曲线绘制、态密度计算以及频率归一化与单位换算四个关键环节展开，并给出可复现的操作与判读方法。

一、高对称路径的选取

色散曲线只能沿布里渊区中的一维路径绘制，因此路径必须经过能代表整个布里渊区特征的高对称点。高对称点是那些在晶体点群操作下保持不变或仅发生等价变换的波矢点，常见标记为 Γ（q=0）、X、M、K、L、W、A、H 等，具体符号取决于布拉维格子类型。例如面心立方（FCC）常用路径为 Γ–X–W–K–Γ–L–U–W–L–K，体心立方（BCC）常用 Γ–H–N–Γ–P，六角格子常用 Γ–M–K–Γ–A–L–H–A。

路径选取不是随意的，应遵循三条原则：第一，路径必须闭合或分段覆盖布里渊区的独立区域，避免遗漏关键极值点；第二，路径应尽量经过已知的对称性特殊点，使色散曲线在端点处出现可解释的简并或分裂；第三，路径分段数（每段 k 点密度）要足够，通常每段取 30～60 个 q 点，总点数在 200～500 之间，既能保证曲线平滑，又不至于计算量过大。以六角晶格为例，若每段取 40 个点，Γ–M、M–K、K–Γ、Γ–A、A–L、L–H、H–A 共 7 段，总点数约 280 个，属于合理范围。

路径选取完成后，建议先检查 q 点是否严格落在高对称点上，尤其是 Γ 点必须精确为 (0,0,0)，否则声学支在 Γ 点附近会出现非物理的微小频率偏移。

二、色散曲线的绘制与判读

色散曲线的横轴是沿路径累积的波矢距离，纵轴是声子频率。绘制时通常将各高对称点作为横轴刻度，段与段之间用竖线分隔。判读色散曲线时重点关注三类特征：

声学支在 Γ 点归零：三维晶体有 3 支声学支，在 Γ 点频率应严格为 0（数值计算中通常小于 1 cm⁻¹ 或 0.1 THz）。若明显不为零，说明力常数或非解析项处理有问题。
光学支与声学支的间隙：若存在带隙，说明该频率区间没有声子模式，对应态密度为零；若带隙消失，则声学支与光学支交叉或简并。
虚频（负频）的出现：若某支在特定 q 点出现负频率，说明该结构在对应波矢下动力学不稳定，通常对应软模相变或结构弛豫不充分。

绘制时纵轴单位常用 THz 或 cm⁻¹，两者换算关系为 1 THz ≈ 33.356 cm⁻¹。若计算输出为 eV，则 1 eV ≈ 241.799 THz ≈ 8065.54 cm⁻¹。绘图前应统一单位，避免同一张图中混用。

三、态密度计算

态密度 g(ω) 的定义是单位频率区间内的声子模式数，满足 ∫g(ω)dω = 3N，其中 N 为原胞中的原子数。实际计算中，态密度通过对布里渊区进行均匀 q 点采样后，将各 q 点的本征频率用高斯或洛伦兹展宽函数累加得到：

g(ω) = (1/N_q) Σ_{q,ν} δ(ω - ω_{q,ν})
     ≈ (1/N_q) Σ_{q,ν} (1/(σ√(2π))) exp[-(ω-ω_{q,ν})²/(2σ²)]

其中 N_q 为 q 点总数，ν 为模式指标，σ 为展宽因子，通常取 0.1～0.3 THz。q 点网格密度直接决定态密度的平滑程度：以 12×12×12 网格为例，若原胞含 2 个原子，则总模式数为 6，N_q = 1728，总采样模式数约 10368，足以得到较平滑的曲线；若网格降到 6×6×6，曲线会出现明显毛刺。实际计算中建议先做收敛测试，比较 8×8×8、12×12×12、16×16×16 三档网格下态密度的差异，当曲线形状基本不变时即可认为收敛。

态密度的判读要点：低频区对应声学支，高频区对应光学支；若某频率区间 g(ω)=0，说明该区间无声子模式；态密度的峰值位置对应色散曲线中平坦的支，即声子群速度接近零的区域。

四、频率归一化与单位换算

不同计算软件输出的频率单位可能不同：VASP 的 phonopy 接口常输出 THz，Quantum ESPRESSO 的 matdyn 输出 cm⁻¹，部分程序输出 eV 或 meV。在绘制色散和态密度之前，必须统一到同一单位。常用换算关系如下表：

单位	换算关系	典型用途
THz	1 THz = 33.356 cm⁻¹ = 4.136 meV	色散曲线纵轴
cm⁻¹	1 cm⁻¹ = 0.02998 THz = 0.12398 meV	红外/拉曼对比
meV	1 meV = 0.2418 THz = 8.066 cm⁻¹	热力学量计算
eV	1 eV = 241.799 THz = 8065.54 cm⁻¹	第一性原理输出

归一化方面，态密度通常按“每原胞每频率区间模式数”归一，即 ∫g(ω)dω = 3N。若软件输出的是“每晶胞”或“每原子”归一，需乘以相应因子。例如原胞含 2 个原子，若输出为每原子态密度，则乘以 2 才得到每原胞态密度。此外，部分程序将态密度按总模式数归一为 1，此时需乘以 3N 才能与色散曲线的模式数对应。

单位换算和归一化错误是声子计算中最常见的低级错误。建议在绘图脚本中显式写出换算因子，并在图注中标注单位与归一化方式，例如“PDOS normalized to 3N modes per primitive cell, frequency in THz”。

五、完整计算流程与结果判读

下面给出一个可复现的示例流程，以六角晶格、原胞含 2 个原子为例。假设已完成力常数计算，得到 FORCE_CONSTANTS 文件。

# 1. 生成高对称路径 q 点
phonopy --dim="4 4 4" --band="0 0 0  0.5 0 0  0.333 0.333 0  0 0 0  0 0 0.5  0.5 0 0.5  0.333 0.333 0.5  0 0 0.5" -c POSCAR

# 2. 计算色散
phonopy --dim="4 4 4" --band="0 0 0  0.5 0 0  0.333 0.333 0  0 0 0  0 0 0.5  0.5 0 0.5  0.333 0.333 0.5  0 0 0.5" --mesh="12 12 12" -c POSCAR

# 3. 提取态密度
phonopy --dim="4 4 4" --mesh="12 12 12" --dos -c POSCAR

计算完成后，检查 band.yaml 中 Γ 点频率是否接近 0，检查 total_dos.dat 的积分是否等于 3N=6。若 Γ 点声学支频率大于 1 cm⁻¹，应检查力常数是否收敛、非解析项修正是否开启。若态密度积分偏离 6 超过 5%，应检查 q 点网格是否足够密、展宽因子是否过大。

最终交付的图通常包含两个面板：左侧为色散曲线，右侧为态密度，两者共享纵轴频率。色散曲线中高对称点用竖线标出，态密度曲线用填充或实线表示。图注中注明计算方法、泛函、q 点网格、展宽因子和单位。

上图展示了六角晶格（2 原子/原胞）的色散曲线与态密度联合判读关系。左侧色散曲线中，低频三条曲线为声学支，在 Γ 点归零；高频三条为光学支，在 Γ 点存在明显分裂。右侧态密度中，低频区对应声学支贡献，高频峰对应光学支的平坦区域，中间带隙处 g(ω)=0。实际计算中，若发现声学支在 Γ 点不归零、态密度积分偏离 3N、或出现虚频，应优先检查力常数收敛性、q 点网格密度和单位换算是否正确。