pymatgen：表面能与功函数实操

表面能来自slab相对体相的能量差；功函数来自真空势与费米能之差。本例实际复用鸢算科技集群两个Cu(111)计算目录的结果，分别调用pymatgen表面分析接口。

两种性质使用各自对应的计算数据，不能把一个slab的LOCPOT和另一个目录的费米能混在一起。

下载实际数据、结构、完整后处理和绘图脚本。

unzip pymatgen-surface-local-tools-wave721-*.zip
cd pymatgen-surface-local-tools
module load pymatgen/2026.9.24
bash results.sh
bash plot.sh

环境为pymatgen2026.9.24、Python3.11、NumPy/SciPy和Matplotlib。run.sh重新后处理到新建的analysis，不提交DFT；reference-data和原图保留，plot.sh画图到figures。再次分析请保留已有目录，在另一份解压目录运行。

1. 准备表面能输入#

cat inputs/source-reference.json
cat inputs/POSCAR-cu-bulk
cat inputs/POSCAR-cu-surface-symmetric

source-reference.json保留实际体相和slab晶胞总能量；POSCAR保存对应几何。使用的是5原子对称Cu(111) slab及4原子Cu体相。

python - <<'PY'
import json
from pymatgen.core import Structure
from pymatgen.entries.computed_entries import ComputedStructureEntry
from pymatgen.analysis.surface_analysis import SlabEntry
with open('inputs/source-reference.json') as f:
    e = json.load(f)['energies_eV']
bulk = ComputedStructureEntry(
    Structure.from_file('inputs/POSCAR-cu-bulk'), e['cu-bulk'])
slab = SlabEntry(Structure.from_file('inputs/POSCAR-cu-surface-symmetric'),
    e['cu-surface-symmetric'], miller_index=(1,1,1))
gamma = float(slab.surface_energy(bulk))
print('表面能 eV/A^2:', gamma)
print('表面能 J/m^2:', gamma*16.021766208)
PY

本例得到约1.327645 J/m²。这里是两个相同表面的对称slab：γ=(E_slab−N_slab E_bulk/N_bulk)/(2A)，A是单个表面的面积。遗漏分母2会得到两倍结果。

这一数值只对应给定slab厚度、真空和原方法，不代表已经完成表面厚度/k网格收敛。不同终止面、非化学计量或吸附体系，还要加入适用的化学势与参考条目。

2. 从同一个计算读取LOCPOT和EF#

若处理自己的slab计算，在包含LOCPOT、OUTCAR、CONTCAR的目录运行：

python - <<'PY'
from pymatgen.io.vasp.outputs import Locpot, Outcar
from pymatgen.core import Structure
from pymatgen.analysis.surface_analysis import WorkFunctionAnalyzer
locpot = Locpot.from_file('LOCPOT')
values = locpot.get_average_along_axis(2)
ef = Outcar('OUTCAR').efermi
structure = Structure.from_file('CONTCAR')
analysis = WorkFunctionAnalyzer(structure, values, ef)
print('EF:', ef)
print('真空势最大值:', analysis.vacuum_locpot)
print('功函数:', analysis.work_function)
print('原生平坦性提示:', analysis.is_converged())
PY

先确认真空方向确实是晶胞c方向。用Vasprun(parse_dos=False)时，efermi可能没有被读取；本例明确从同目录OUTCAR读取，而不把None或其他目录EF代进去。

下载包保存的是该实际LOCPOT的平面平均数据，因而学习后处理无需重新读整个三维网格：

python - <<'PY'
import json, numpy as np
from pymatgen.core import Structure
from pymatgen.analysis.surface_analysis import WorkFunctionAnalyzer
with open('inputs/source-reference.json') as f:
    ef = json.load(f)['workfunction_EF_eV']
d = np.genfromtxt('inputs/Cu111-planar-potential.csv',delimiter=',',names=True)
s = Structure.from_file('inputs/POSCAR-cu-workfunction')
a = WorkFunctionAnalyzer(s, d['potential_eV'], ef)
print(a.work_function)
PY

3. 本例实际功函数#

cat reference-data/Cu111-work-function.csv
head -n 8 reference-data/Cu111-planar-potential.csv

即Φ=V_vac−EF，约4.884 eV。本例另外在离开最外层原子2Å之外取真空窗口，检查平台波动和平台平均值；最大值与该平台相符。

原生算法取一维势的最大值。若最大值出现在slab内部、真空没有平台，或两侧因偶极产生不同平台，就不能直接把这个单数值当作可靠的双面功函数。应检查电势、真空厚度、偶极修正及两侧参考。

CSV中z坐标按周期网格i·c/N生成。库内部某些显示坐标采用含端点的线性数组，显示位置可能有一个网格量级差别，能量相减本身不受这一位置标签影响。

4. 换自己的体系#

准备同方法的体相/slab能量及对应原子数，明确两个表面是否等价；功函数则保持LOCPOT、结构和EF来自同一计算。确认真空方向和平台后，再比较不同表面或吸附前后的变化。

已有实际输入与计算步骤见VASP功函数与电荷，形貌构造见Wulff实操。接口见pymatgen表面分析。