📚 课程目录第4章计算前的结构准备与收敛性工程

4. 计算前的结构准备与收敛性工程

介电谱计算对结构、基组和布里渊区采样极其敏感。同一材料，若截断能偏低、k点过稀或展宽参数选择不当，虚部峰位可能偏移0.2 eV以上，峰高甚至相差数倍。本章围绕晶体结构建模、原胞与超胞选择、截断能、k点网格、展宽参数五个关键输入，说明它们如何影响介电谱收敛性，并给出一套可自动化的收敛测试方案。

一、结构建模：从晶体对称性到计算原胞

介电张量由晶体对称性决定。立方晶系只需一个独立分量ε_xx；六方晶系有ε_xx=ε_yy与ε_zz两个独立分量；单斜与三斜晶系则需完整张量。因此建模第一步不是"把原子放进去"，而是确认空间群与点群，据此决定后续需要计算哪些分量、是否需要非对角元。

结构文件通常来自CIF或实验文献。导入后必须检查：晶格常数与角度是否与文献一致、原子占位是否完整、是否有未识别的对称操作。以金红石TiO₂（P4₂/mnm，四方）为例，原胞含2个Ti和4个O，晶格常数a=b=4.594 Å、c=2.959 Å。若误用惯用胞（含4个Ti、8个O）计算，虽然能量收敛值相同，但k点路径与介电张量分量会因坐标系选择不同而需要重新对应。

二、原胞与超胞：不是越大越好

介电谱计算原则上使用原胞。原因有三：原胞原子数最少，相同截断能下平面波基组规模最小；原胞的倒格矢最短，k点网格密度更容易控制；介电张量的对称性分析在原胞下最直接。

超胞仅在两种情况下使用：一是掺杂或缺陷体系，需要构造足够大的超胞以避免周期性镜像相互作用；二是需要引入真空层模拟二维材料或表面。此时必须注意，超胞会折叠布里渊区，原本在原胞中不同k点的跃迁会被折叠到同一k点，介电谱的峰位识别会变得困难。若必须用超胞，建议同时计算原胞作为参照，确认主要跃迁峰位置一致后再分析超胞结果。

三、截断能与k点网格：两个必须联调的参数

截断能E_cut决定平面波基组大小，直接影响价带和导带波函数的描述精度。E_cut不足时，导带高能态首先失真，介电谱虚部在高能区（通常>10 eV）出现虚假振荡或峰位偏移。k点网格则决定布里渊区积分精度，网格过稀时，介电谱虚部的低能峰（带边附近）峰高被低估、峰形不对称。

两者不能独立收敛。E_cut提高后，波函数更精确，k点积分误差的相对重要性会变化；反之，k点加密后，E_cut不足带来的误差可能被部分平均掉，造成"假收敛"。因此必须做二维收敛测试：固定k点，扫描E_cut；再固定E_cut，扫描k点；最后在两者都较高的水平上交叉验证。

四、展宽参数：物理展宽与数值展宽的平衡

介电谱虚部ε₂(ω)在独立粒子近似下是一系列δ函数的叠加。实际计算中必须引入展宽η将其平滑为连续曲线。η的物理来源包括电子-声子散射、缺陷散射和仪器分辨率，但在第一性原理计算中，η更多是一个数值参数。

η过小（如0.01 eV），谱线保留大量尖锐毛刺，峰位难以自动识别；η过大（如0.5 eV），相邻峰合并，带边信息丢失。经验做法是取0.05～0.15 eV，并做η扫描：若η从0.1 eV降到0.05 eV时峰位漂移小于0.05 eV、峰高变化小于5%，则认为谱线对η不敏感，可固定η=0.1 eV用于后续生产计算。

五、自动化测试方案

手工逐组参数测试效率低且易漏。推荐用脚本生成参数矩阵并批量提交。以下为Python伪代码示例，展示如何组织E_cut与k点的二维扫描：

import itertools, subprocess, os

ecut_list = [400, 500, 600, 700, 800]      # eV
kmesh_list = [(2,2,2), (4,4,4), (6,6,6), (8,8,8)]
eta_list = [0.05, 0.10, 0.20]              # eV

for ecut, kmesh, eta in itertools.product(ecut_list, kmesh_list, eta_list):
    tag = f"E{ecut}_K{kmesh[0]}_eta{eta}"
    workdir = f"./conv/{tag}"
    os.makedirs(workdir, exist_ok=True)
    # 生成输入文件：结构 + 截断能 + k点 + 展宽
    with open(f"{workdir}/input.in", "w") as f:
        f.write(f"ecut = {ecut}\n")
        f.write(f"kmesh = {kmesh[0]} {kmesh[1]} {kmesh[2]}\n")
        f.write(f"eta = {eta}\n")
    # 提交计算（此处以占位命令示意）
    subprocess.run(["dft_dielectric", "input.in"], cwd=workdir)

计算完成后，用脚本提取每组参数下的ε₂(ω)曲线，自动识别主峰位置与峰高。判定收敛的量化标准建议为：相邻两档E_cut（如600与700 eV）的峰位差小于0.05 eV，相邻两档k点（如6×6×6与8×8×8）的峰高差小于5%。若同时满足，则取较低档参数作为生产参数，兼顾精度与效率。

实用建议：对于含过渡金属氧化物的体系（如TiO₂、ZnO），E_cut通常需要600 eV以上才能收敛；对于主族半导体（如Si、GaAs），400～500 eV即可。k点网格方面，金属体系需要比半导体更密的网格（如12×12×12以上），因为费米面附近的积分对k点更敏感。

常见陷阱：不要用"总能收敛"代替"介电谱收敛"。总能对E_cut和k点的敏感度远低于介电谱虚部的高能区。若只检查总能，可能在E_cut=400 eV、k=4×4×4时就误判收敛，导致介电谱在8 eV以上出现明显误差。

完成上述收敛测试后，应保存一份参数记录表，至少包含：材料名称、空间群、原胞原子数、E_cut、k点网格、展宽η、主峰位置、主峰高度、收敛判定结论。这份记录既是后续计算的起点，也是结果可复现性的依据。