DFT计算表面能

项目过程：Cu 低指数晶面表面能的系统评估

这个项目的客户做铜基催化载体，需要知道 Cu(111)、(100)、(110) 三个低指数面的表面能，用来判断哪个晶面在反应条件下更稳定、更适合作为活性位点载体。他的核心诉求不是算一个绝对精确的数值，而是要一个相对可靠的排序和厚度收敛依据。

我采用的是 VASP + PBE 的平板模型方案。先对 Cu 块体做优化，得到晶格常数 a≈3.61 Å，再切出 2×2 的 (111)、(100)、(110) 平板。关键参数我设置为：平面波截断能 520 eV，K 点取 11×11×1，真空层厚度 15 Å，平板层数从 3 层扫到 12 层，每次上下表面各留一层固定来模拟体相约束。能量收敛判据是 10⁻⁵ eV，离子弛豫力收敛到 0.02 eV·Å⁻¹。

表面能 γ 按 γ = (E_slab − n·E_bulk)/(2A) 计算，其中 A 是单胞表面积，分母 2 是因为平板有上下两个表面。我还做了偶极校正，避免不对称弛豫带来的人工电场。

案例分析：从收敛曲线看哪个晶面最稳定

第一张案例图给出的是 γ 随平板层数的变化。Cu(111) 从 3 层的约 1.69 J·m⁻² 下降到 9 层的 1.28 J·m⁻²，之后基本不再变化；Cu(100) 从 1.82 收敛到 1.40 J·m⁻²。两条曲线都在 9 层附近进入平台区，说明 9 层是这套参数下的收敛厚度。

第二张案例图把 DFT-PBE 结果和实验值做了柱状对比。PBE 算出来的 (111)、(100)、(110) 表面能分别是约 1.28、1.40、1.47 J·m⁻²，实验值则分别是约 1.85、1.88、1.95 J·m⁻²。PBE 整体低估了 30% 左右，这是典型的 GGA 自洽误差；但各晶面之间的相对顺序和各向异性比 γ₁₁₀/γ₁₁₁=1.15 与实验非常一致。

总结教训：PBE 绝对值偏低，但相对排序可信

第一，平板厚度必须做收敛测试。3 层平板会把表面能高估 30% 以上，9 层基本可靠，再往上加层数性价比很低。

第二，真空层 15 Å 对金属表面足够。我试过 10 Å，由于镜像相互作用，γ 会偏高约 0.05 J·m⁻²，15 Å 是一个稳妥的折中。

第三，绝对表面能不能和实验直接比较。PBE 低估了约 0.5 J·m⁻²，但在催化剂筛选里我们更关心 γ₁₁₀/γ₁₁₁ 这种相对量，这个数据可以直接用。

项目经验：把表面能变成晶面选择依据

如果客户的目的是预测纳米颗粒暴露面，可以把算出来的 γ 接到 Wulff 构造里。Cu(111) 表面能最低，所以在热力学平衡下它会占最大暴露面积；(110) 能量最高，往往出现在台阶或缺陷处，反而可能是活性位点。这个项目让我再次确认：DFT 表面能的强项不是给出实验值的精确复现，而是给出能量排序和物理趋势。

这篇文章把我在这个项目里用到的计算流程、参数取舍和图上读数做了完整复盘。如果你也在做类似方向，可以参考我的思路：DFT计算表面能。