使用DFT计算功函数

本教程采用AMS2023完成，采用Python脚本处理计算生成的band.rkf文件得到。AMS2024版以后，将集成到图形界面中，通过图形窗口操作即可可以直接查看结果。

对于同种材料，BAND与平面波方法相比，区别不大，但对于异质结，二者区别很大。

平面波方法无法真正处理二维周期性体系，因此不得不通过在Z方向添加足够大的真空，来模仿二维周期性边界。但由于静电势是非常长程的作用，乃至50埃以上仍然有所影响，因此平面波方法不得不进行偶极较正。但这并不能真正解决问题，得到的静电势分布，离表面很近的区域静电势就处于水平状态，这意味着电子在这些区域（例如下图的例子，在z轴±13埃以外，即仅仅距离材料表面几埃以外）上下移动，既不消耗能量也不获得能量，而这应该是违背常识的。例如下图：

BAND得到的静电势分布，更接近真实。

计算功函数，模型一般为二维周期性结构，即在XY方向无限延申、重复，而Z方向以外，则为半无限大真空。如果对周期性边界条件理解不充分，建议参考：周期边界条件。

本例中演示一个界面体系：Al的(100)面与LiF的(100)面形成界面，二者各自存在一定的厚度，然后计算两个面各自的功函数。

形成界面，对宏观视角来讲，并没有太复杂的图景，但是在微观层面，尤其是在计算层面则稍微复杂、困难，因为两种物质有各自不同的周期性，而第一性原理计算、分子动力学计算等原子层面、电子层面的计算模拟，都要求体系只有一种周期，不允许这个体系存在两种周期，两种周期是无法进行计算的。

首先界面的建模是基于两个二维晶格，晶格常数分别为A1、B1与A2、B2，夹角为α1、α2。

如何把两种周期变为一种周期呢？思想非常简单：各自晶格常数，分别取最小公倍数形成超胞，得到各自的超胞晶格矢量，如果两种物质各自的超胞晶格矢量相同，那就可以合并为相同的晶格矢量，从而能够进行计算。在实际操作中，两种物质的晶格矢量毕竟不是整数，所以很难、几乎不可能真正得到完全相同的超胞晶格矢量，因此在实际计算中，只要二者之间的差异不太大，从而带来的误差不太大，也就可以忍受了。这种情况下，往往以其中一个超胞晶格矢量为准，另一个可以微调（稍微等比放大/缩小键长），最终的键长放大/缩小量小于0.1 Å，问题就不太大，当然越小越好。

上面说的是最简单的情况，也就是两种物质的晶格夹角 α 相同，因此较为简单。而不同夹角，也需要类似调整，最终的效果就是二者的 A、B 与 α 相同或非常接近。

本例中的模型，采用已经创建好的模型（点击打开链接，复制全部文字，然后在AMSinput中Ctrl V粘贴即可）。

下面采用了DZP基组，对于重一些的元素，用户可以自行设置较高基组，积分精度设置的较高，精度会更高一些：

设置一定的电子温度，促进收敛。一般电子温度越低越好，对精度的损坏越小，不过电子温度越高，促进收敛的效果越明显，因此一般尽量用小的电子温度，只要能收敛就可以，不要再提高电子温度了。

保存并运行作业。

SCM → View可以显示静电势（对xy平面求均值）沿着z方向的分布情况，如此可以得到静电势在Slab两侧的最高值，然后结合费米能级的数值，即可计算得到功函数。

首先我们需要修改格点的情况，因为本例需要离Slab很远区域，例如100埃范围的数据，而格点并不需要很细，因此我们将默认格点的密度降低为原先的1/4（即默认的0.0529177210903乘以5，得到0.2645886054515，当然不降低格点密度也可以，但是就会特别特别慢，而且计算静电势的分布的XY平面的均值，实际并不要求很细的格点），因此View窗口 → Field → Grid → Details，修改数据如下：

点击Use即生效。

菜单栏Add → Graph along line，得到一个窗口，该窗口可以生成各种函数沿着某个指定空间线段内的分布情况。首先我们设定线段的情况：

含义如下：

• 80，默认显示函数的范围一般比较小，比如原子外围几埃的范围，我们需要显示很远范围的静电势分布，因此设置了80埃，即Slab往上80埃、往下80埃，加上Slab本身的厚度，应该约170多埃的范围。注意不要超过我们前面设定的100埃上限。
• 1000，这条静电势分布曲线，我们希望用1000个点，这样点数较多，曲线会足够光滑。
• line along z-axes，曲线沿着Z轴画。因为这是一个二维周期体系，Z轴正好垂直于Slab。
• Average，对静电势在XY平面求平均，从而每一个坐标z处都对应一个静电势的均值。

设置完毕，上方的Select Field选择Potential → Coulomb Potential，于是显示静电势的分布情况，默认只显示了约20埃的范围：

点击Do it 按钮，将在现在的基础上，左右各延长80埃：

通过鼠标滚轮放大缩小，左键拖动，可以找到左侧的最高点和右侧的最高点，双击纵坐标，（可以将默认的单位Ha修改为eV）分别为：0.4995eV、-0.2386eV。

SCM → logfile，尾部可以看到：

<Nov06-2024> <22:02:40>  FERMI ENERGY:          -0.1425 A.U.
<Nov06-2024> <22:02:40>                         -3.8776 E.V

二者分别减去费米能级：

• 0.4995eV -( -3.8776) eV = 4.3771 eV
• -0.2386eV -( -3.8776) eV = 3.6390 eV

因为曲线是沿着z轴，因此左侧对应Al的区域（Z轴的下方），右侧对应LiF（Z轴的上方部分），即知Al(100)表面与真空接触一侧功函数为4.38 eV，LiF(100)表面与真空接触一侧功函数为3.64 eV。

而其他文献中计算（前两行数据，分别使用VASP、Quantum ESPRESSO与本教程对比（第三行数据）：

两个其他文献的计算结果来自文献：

• Prada, U Martinez, G. Pacchioni. Work function changes induced by deposition of ultrathin dielectric films on metals: A theoretical analysis. Phys. Rev. B. 78, 235423. DOI: 10.1103/PhysRevB.78.235423
• Kondo, T. Matsushista. Vacuum-Level Shift at Al/LiF/Alq3 Interfaces: A First-Principles Study. ACS Omega 2019, 4, 8, 13426–13434. DOI: 10.1021/acsomega.9b01667