vaspkit计算功函数

我做过一个金属电极-半导体界面项目，客户最关心的是界面处电子怎么流动，而功函数 Φ 是判断能带对齐的第一把尺子。那个项目里我用 VASPKIT 批量算了 Cu、Ag、Au 几个低指数晶面的功函数，左侧那张平均静电势图和右侧的柱状对比图，就是我交付报告里的核心两张图。

功函数的定义很简单：Φ = E_vac − E_F，其中 E_vac 是真空能级，E_F 是费米能级。但真要算准， slab 模型的构建非常关键。我先用 VASP 优化了 Cu(100)、Cu(110)、Cu(111)、Ag(111)、Au(111) 的体相晶格常数，Cu 的 PBE 晶格常数我得到 3.61 Å，和实验 3.61 Å 几乎一致。然后切表面，slab 厚度我取 5 层原子，上下表面用对称结构，真空层 20 Å，确保两侧静电势平台区不重叠。

左侧图展示的是 Cu(111) 沿表面法向 z 的平均静电势。可以看到 slab 内部势场周期性振荡，对应原子层位置；到了真空区，势场趋于平坦的平台，那就是 E_vac。E_F 在图中用红色虚线标出，二者差值 Φ = 3.10 eV 这个标注是我用 VASPKIT 的 POTCAR 与 PLANAR_AVERAGE 功能直接读出来的。VASPKIT 的 4xx 功能模块可以一键读 LOCPOT 文件，自动找真空区电势平均值，避免手工读平台区时人眼判断的误差。

右侧图是晶面对比。Cu(100)=4.59 eV，Cu(110)=4.48 eV，Cu(111)=4.98 eV，Ag(111)=4.74 eV，Au(111)=5.31 eV。这个数据有几个看点：第一，Cu 的密排面 (111) 功函数最高，这和表面电子云 smoother、偶极势更小有关；第二，Au(111) 5.31 eV 是五个面里最高的，适合做高功函数电极；第三，晶面取向造成 Cu 体系内部差异就达 0.5 eV，如果客户用多晶 Cu 做电极，实际接触势差会是一个分布，不能单取一个值。

具体操作上，我踩过两个坑。第一个是 slab 不对称导致的偶极修正问题：如果 slab 上下表面不等价，真空区电势会出现倾斜，E_vac 没有唯一平台。我后来改成对称 slab，并在 INCAR 里加 LDIPOL = .TRUE. 和 IDIPOL = 3 做偶极修正，平台区才平。第二个是 K 点采样：表面方向我用 12×12×1，z 方向只用 1 个 K 点，因为真空方向没有周期性；如果 z 方向设多了，不仅浪费机时，还会引入数值噪声。第三个是 smearing：金属表面用 Methfessel-Paxton 一阶，SIGMA 取 0.1 eV，保证 E_F 收敛。

我的经验是，功函数计算不要只给一个数，要同时给误差来源：晶面取向、slab 厚度、交换关联泛函。PBE 对 Cu、Ag、Au 的功函数通常低估 0.1-0.2 eV，如果客户要求更高精度，可以上 HSE06 或 meta-GGA。项目最后我也建议客户，如果真正做界面能带对齐，还需要把半导体 slab 放上去算异质结，单纯金属功函数只能给出参考范围。

更多 VASP 表面与界面计算的实战细节，我汇总在 [vaspkit计算功函数](https://www.keyanxueshu.com/category/dft/) 栏目里。[vaspkit计算功函数](https://www.keyanxueshu.com/)

再深入一层：右侧柱状图里 Cu(111) 的功函数 4.98 eV 比 Cu(100) 高 0.39 eV，这个差值直接决定同一金属不同晶面与半导体接触时的 Schottky 势垒高度。比如 n 型半导体电子亲和能 4.0 eV，用 Cu(110) 做电极时界面势垒约 0.48 eV，用 Cu(111) 就变成 0.98 eV，载流子注入效率可能差一个数量级。实际器件 electrode 往往是多晶，晶面分布会让功函数有一个统计分布，我们做界面能带对齐时通常会给一个 0.2-0.3 eV 的误差条。另外，左侧电势平台区不要只读一个点，我会在真空区取 5-10 Å 宽度内的平均值，再和标准差一起报给客户，避免平台区有轻微倾斜导致读数偏差。