3、第一性原理计算：平面波方法、赝势与投影缀加波（PAW）方法、k点采样、收敛性测试与参数设置

第一性原理计算（First-Principles Calculations）的核心目标是在不引入经验参数的前提下，通过求解量子力学方程来预测材料的电子结构、总能量、力学性质与反应路径。在众多实现路径中，基于Kohn-Sham密度泛函理论（DFT）的平面波方法已成为凝聚态物理与材料科学领域最主流的工具。本章将系统阐述平面波基组的数学基础、赝势与PAW方法的物理图像、k点采样的布里渊区积分策略，以及收敛性测试的标准流程与参数设置规范，帮助学习者建立从理论到实操的完整知识链条。

3.1 平面波基组与截断能

核心原理：从Kohn-Sham方程到平面波展开

在Kohn-Sham DFT框架下，单电子波函数满足：

[-ħ²/(2m)∇² + V_eff(r)] ψ_i(r) = ε_i ψ_i(r)

对于周期性晶体，Bloch定理允许将波函数写为平面波与周期函数的乘积：

ψ_{nk}(r) = exp(ik·r) · u_{nk}(r) = Σ_G c_{n,k+G} exp(i(k+G)·r)

其中G为倒格矢，k为布里渊区内的波矢。平面波基组的优势在于：完备性好、无基组叠加误差（BSSE）、力与应力计算简单、收敛性系统可控。然而，由于波函数在原子核附近振荡剧烈，直接使用平面波展开需要极高的截断能（Ecut），计算成本难以承受。因此，平面波方法必须与赝势或PAW方法配合使用。

截断能Ecut的物理意义与选择策略

截断能定义了平面波展开的最大动能：

Ecut = ħ²|k+G|²/(2m)

Ecut越大，基组越完备，但计算量按O(N³)增长。实际选择需平衡精度与效率：

• 经验规则：对于大多数体系，Ecut取势函数推荐值的1.2~1.5倍即可获得良好收敛；
• 系统测试：以总能变化<1 meV/atom为标准，逐步增加Ecut（如从400 eV到600 eV，步长50 eV）；
• 性质依赖：结构优化和弹性常数计算需要更严格的Ecut（力收敛），而能带结构对Ecut的敏感度相对较低。

案例拆解：硅晶体的Ecut收敛测试

以金刚石结构的Si晶体为例（晶格常数a=5.43 Å），使用PAW势Si（标准版本，ENMAX=245 eV）。测试流程如下：

1. 初始设置：使用Gamma-centered 8×8×8 k点网格，ISMear=0（高斯展宽），SIGMA=0.1 eV；
2. Ecut扫描：分别计算Ecut=200, 250, 300, 350, 400, 450, 500 eV下的总能；
3. 数据分析：记录总能变化量，当Ecut从400 eV增至450 eV时，总能变化为0.8 meV/atom，满足<1 meV/atom的收敛标准；
4. 最终选择：取Ecut=450 eV（约1.8×ENMAX），确保结构优化时力收敛精度达到0.01 eV/Å。

该测试耗时约30分钟（8核并行），但避免了因Ecut不足导致的结构参数偏差（如晶格常数误差>0.5%）和弹性常数计算失败。

3.2 赝势与投影缀加波（PAW）方法

核心原理：从全电子到有效势

赝势的基本思想是：原子核附近的芯电子不参与化学成键，其波函数可以用平滑的赝波函数替代，同时保持价电子在芯区外的散射性质不变。根据构造方式，赝势分为：

PAW方法通过线性变换将全电子波函数与赝波函数联系起来：

|ψ_nk⟩ = |ψ̃_nk⟩ + Σ_i (|φ_i⟩ - |φ̃_i⟩) ⟨p̃_i|ψ̃_nk⟩

其中|φ_i⟩为全电子部分波，|φ̃_i⟩为赝部分波，|p̃_i⟩为投影函数。PAW方法在保持平面波高效性的同时，恢复了全电子计算的精度，尤其适用于过渡金属、稀土元素和磁性体系。

判断方法：如何选择赝势/PAW势

• 检查势文件版本：VASP的POTCAR文件应注明生成日期和交换关联泛函（如PBE、PBEsol、SCAN）；
• 核对电子配置：确保势文件包含正确的价电子数（如Fe的PAW势应包含3p、3d、4s共16个电子）；
• 对比测试：对同一体系分别使用不同势文件计算晶格常数和体模量，与实验值对比；
• 注意软硬程度：软势（低ENMAX）适合大体系初步筛选，硬势（高ENMAX）适合高精度计算。

实战场景：过渡金属氧化物的PAW势选择

以α-Fe₂O₃（赤铁矿）为例，其反铁磁基态和强d电子关联对势函数敏感。推荐流程：

1. 使用Fe_pv（3p半芯态）和O标准势，设置ISPIN=2，初始磁矩Fe=4.2 μB；
2. 测试Fe（仅3d、4s）与Fe_pv的差异，发现Fe势高估晶格常数约1.2%；
3. 最终采用Fe_pv + U=4.0 eV（J=0.9 eV）的GGA+U方案，获得与实验吻合的带隙（2.1 eV）和磁矩（4.1 μB）。

3.3 k点采样与布里渊区积分

核心原理：离散化布里渊区

周期性体系中的物理量（总能、电荷密度、态密度）需要对布里渊区（BZ）内的k点进行积分。由于波函数和本征值在k空间平滑变化，可用有限k点网格近似积分：

E_total = Σ_n ∫_BZ ε_nk f(ε_nk) dk ≈ Σ_n Σ_k w_k ε_nk f(ε_nk)

常用的采样方案包括：

• Monkhorst-Pack网格：均匀分布，适用于大多数半导体和绝缘体；
• Gamma-centered网格：包含Γ点，适用于金属体系（避免费米面附近奇异性）和六方晶系；
• 四面体方法（Blöchl修正）：对绝缘体态密度计算更精确。

k点密度与体系尺寸的关系

k点网格密度应随体系尺寸增大而减小。经验公式：

N_k ∝ 1/L（L为超胞边长）

例如，对于原胞（L≈5 Å）需要8×8×8网格，而2×2×2超胞（L≈10 Å）仅需4×4×4网格即可达到相同精度。对于表面体系，仅在平行表面方向需要密集k点，垂直方向只需1个k点（或少量k点）。

案例拆解：铜金属的k点收敛测试

Cu是典型金属，费米面穿越d带，需要较密k点。测试方案：

结果表明，10×10×10网格已满足<0.1 meV/atom的收敛标准。但考虑到后续需要计算弹性常数（对k点敏感），最终选择12×12×12网格。同时设置ISMEAR=1（Methfessel-Paxton）和SIGMA=0.2 eV，确保费米面附近占据数平滑。

3.4 收敛性测试与参数设置规范

标准测试流程

1. 结构预处理：使用较低精度（Ecut=1.2×ENMAX，中等k点）进行初步结构优化，获得合理初始构型；
2. Ecut收敛测试：固定k点（取较高密度），扫描Ecut；
3. k点收敛测试：固定Ecut（取收敛值），扫描k点密度；
4. 联合验证：使用收敛的Ecut和k点，计算目标性质（如晶格常数、带隙），与实验或文献对比；
5. 其他参数检查：包括smearing宽度（SIGMA）、力收敛阈值（EDIFFG）、电子步收敛（EDIFF）、真空层厚度（表面体系）等。

参数设置速查表

常见误区与复盘

实战案例：钙钛矿SrTiO₃的完整参数设置

SrTiO₃是典型的立方钙钛矿绝缘体（带隙约3.2 eV），计算流程如下：

1. 结构优化：使用PBE泛函，Ecut=500 eV（1.4×ENMAX），Gamma-centered 6×6×6 k点，ISMEAR=0，SIGMA=0.05 eV，EDIFFG=-0.02 eV/Å；
2. 收敛性验证：将Ecut增至600 eV，总能变化0.5 meV/atom；k点增至8×8×8，总能变化0.3 meV/atom，确认收敛；
3. 性质计算：在收敛参数下计算能带结构（沿高对称路径Γ-X-M-Γ-R），态密度（使用四面体方法ISMEAR=-5），以及弹性常数（施加应变±1%，计算应力-应变关系）；
4. 结果对比：优化晶格常数a=3.905 Å（实验值3.905 Å），带隙2.1 eV（PBE低估，符合预期），体模量B=183 GPa（实验值183 GPa）。

该案例表明，系统化的收敛性测试虽然增加了前期计算时间（约2小时），但确保了后续所有性质计算的可靠性，避免了因参数不足导致的返工。

3.5 总结与工程实践建议

平面波DFT计算的参数设置本质上是精度-效率权衡的决策过程。核心原则包括：

• 先收敛后生产：任何生产性计算前必须完成Ecut和k点收敛测试，并记录测试数据；
• 体系定制化：金属、半导体、绝缘体、磁性体系、表面体系各有不同的参数敏感点；
• 势函数一致性：同一计算系列中，所有结构必须使用相同的势文件、泛函和截断能；
• 文档化记录：保存每次计算的INCAR、KPOINTS、POTCAR版本和收敛测试数据，便于复现和排错。

在实际项目中，建议建立标准化的输入文件模板，并编写脚本自动执行收敛性测试。例如，使用Python脚本循环修改Ecut和k点密度，提取总能数据并生成收敛曲线。这不仅提高了效率，也减少了人为错误。最终，通过严格的收敛性测试和合理的参数设置，第一性原理计算能够为材料设计与性能优化提供可靠的定量预测。