CP2K 表面功函数：Al薄膜、真空势与k点

功函数是把电子从材料的费米能级移到真空所需的能量。这里用理想对称Al(001)薄膜演示：准备结构和输入，提交单点计算，读取两侧真空平台，计算功函数，最后比较真空厚度、k点与薄膜层数。结构固定、未经弛豫；表中数值用于理解流程，不能直接作为已收敛的Al材料常数。

1. 下载后先看已有结果#

下载完整输入、全分辨率网格与绘图脚本，解压到新目录：

python3 -m zipfile -e cp2k-Al-surface-workfunction-wave299-*.zip .
cd cp2k-Al-surface-workfunction-wave299
bash preview.sh
bash run.sh

preview/重画参考图，results/重新读取12组完整网格并输出CSV和JSON，plots/保存结果图。两条命令只处理已有数据，不重新提交DFT。再次运行请换新解压目录，避免覆盖。

资料包包括全部12组input.inp、slab.xyz，CP2K基组/赝势数据及许可说明、完整自有NPZ网格和最小读取程序。不需要另找第三方画图插件。后处理用pymatgen/2026.9.24中的Python3.11、NumPy2.3.5、Matplotlib3.11.2与ASE3.29。

2. 建立自己的计算目录#

以下选3层、18 Å真空间隔、12×12×1网格的实际输入。每层2个Al，共6个原子。vacuum=18 Å指周期重复薄膜的最外原子面之间的间隔；薄膜自身厚4.05 Å，因此晶胞z长22.05 Å。

mkdir my-Al
cp inputs/Al3L-v18-k12/* my-Al/
cp DATA/BASIS_MOLOPT DATA/GTH_POTENTIALS run.slurm my-Al/
cd my-Al
ls input.inp slab.xyz BASIS_MOLOPT GTH_POTENTIALS run.slurm
nano input.inp

完整slab.xyz：

6
Ideal fixed Al(001), unrelaxed symmetric slab
Al 0.000000000000 0.000000000000 9.000000000000
Al 2.863782463806 0.000000000000 9.000000000000
Al 1.431891231903 1.431891231903 11.025000000000
Al 4.295673695708 1.431891231903 11.025000000000
Al 0.000000000000 0.000000000000 13.050000000000
Al 2.863782463806 0.000000000000 13.050000000000

完整input.inp：

&GLOBAL
  PROJECT AL-SURFACE
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME GTH_POTENTIALS
    &QS
      METHOD GPW
      EPS_DEFAULT 1.0E-12
    &END QS
    &MGRID
      CUTOFF 400
      REL_CUTOFF 50
      NGRIDS 5
    &END MGRID
    &POISSON
      PERIODIC XYZ
    &END POISSON
    &KPOINTS
      SCHEME MONKHORST-PACK 12 12 1
      SYMMETRY F
    &END KPOINTS
    &SCF
      EPS_SCF 1.0E-6
      MAX_SCF 200
      SCF_GUESS ATOMIC
      ADDED_MOS 20
      &DIAGONALIZATION
        ALGORITHM STANDARD
      &END DIAGONALIZATION
      &SMEAR ON
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE 300
      &END SMEAR
      &MIXING
        METHOD BROYDEN_MIXING
        ALPHA 0.1
        NBROYDEN 8
      &END MIXING
      &PRINT
        &RESTART OFF
        &END RESTART
      &END PRINT
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
    &PRINT
      &E_DENSITY_CUBE
        STRIDE 1 1 1
        FILENAME =density.cube
      &END E_DENSITY_CUBE
      &V_HARTREE_CUBE
        STRIDE 1 1 1
        FILENAME =potential.cube
      &END V_HARTREE_CUBE
    &END PRINT
  &END DFT
  &SUBSYS
    &CELL
      ABC 5.727564927611 2.863782463806 22.050000000000
      PERIODIC XYZ
    &END CELL
    &TOPOLOGY
      COORD_FILE_NAME slab.xyz
      COORD_FILE_FORMAT XYZ
    &END TOPOLOGY
    &KIND Al
      ELEMENT Al
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q3
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

关键设置：

这块薄膜上下表面对称、中性。单面吸附、不对称薄膜或带电模型需要另外处理偶极、电场与边界条件，不能直接套用当前左右真空窗口。

3. 提交，然后检查真正的输出#

sbatch run.slurm
squeue -u "$USER"

模板实际运行：

module purge
module load cp2k/2026.2-p100
export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1
cp2k.psmp -i input.inp -o output.log

这是CP2K2026.2、GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33构建，当前输入使用1个MPI进程、1线程。模块带GPU支持不等于每个算符都在GPU上执行。

作业结束后：

grep -E 'SCF run converged|Fermi|ENERGY' output.log | tail -n 12
ls -lh density.cube potential.cube

应看到电子收敛以及两份实际CUBE。不能仅凭作业退出或文件存在就认定功函数已收敛。

4. 用脚本提取自己的结果#

回到解压目录，读取新计算：

cd ..
module purge
module load pymatgen/2026.9.24
python read_my_surface.py my-Al
cat my-Al/surface-results/workfunction.json
head my-Al/surface-results/potential-profile.csv

脚本读取my-Al/density.cube、potential.cube和output.log，输出完整z方向势/密度CSV、电子数积分、左右平台与功函数JSON。仅适用于本例这种z为表面法向的正交晶胞和默认真空窗口；换结构后先检查窗口位置。

CP2K2026.2原生CUBE采用固定13字符数值字段，负数有时紧挨前一字段。随包read_cube.py按原始字段完整读取，并用ASE交叉检查；不会修改原始文件，也不截掉密度中极小的负数尾部。

计算顺序：

1. 对每个z网格点，在整个xy平面平均势，得到V(z)。
2. 将CUBE势的Hartree单位乘27.211386245988换成eV。
3. 读取output.log中的费米能，同样统一为eV。
4. 取左侧z=1–3 Å与镜像右侧真空窗口的平均电子势能；图中阴影对应实际选区。
5. 分别计算W_left = Vvac_left - EF和W_right = Vvac_right - EF。

势的符号不能凭文件名猜。 2026.2的V_HARTREE_CUBE写出的量与通常静电势φ符号相反；对电子而言的势能−eφ正对应此输出乘Hartree→eV的系数。因此这里不再额外取负。它也不是把含交换关联项的任意总势直接当作真空参照。

3层18 Å、12×12×1本例得到左右平均功函数约4.414262 eV。左右表面对称，所以两个值接近；左右相同本身并不能证明面内k点或厚度已收敛。Al使用q3赝势，6个原子的价电子积分应接近18，本例为18.000089；不要拿包含冻结芯电子的总核电荷78作积分目标。

5. 真空、k点与层数要分开比较#

全部条件均为同一种PBE、DZVP-MOLOPT-SR-GTH/GTH-PBE-q3、400/50 Ry、300 K电子展宽。下面直接列出12次实际结果；k列表示k×k×1。

在3层、2×2×1固定时，真空18→24 Å变化约0.000243 eV，说明这项真空对照变化很小。但2→4面内网格却变化约0.592609 eV：真空足够不代表k点足够。

采用0.05 eV作为这组示例的相邻网格变化参照：

所以主示例3层最新一对变化已变小，但5层尚需继续检查。不能据此声称厚度、基组、展宽或材料功函数全部收敛。实际研究还要优化表面几何、增加层数、检查更密面内网格、基组/截断与占据展宽，并固定同一能量参照。

本组两侧平台范围均小于0.05 eV、左右差小于0.02 eV，价电子积分误差小于0.005；这些检查能发现读取/选区问题，却不能替代上述物理收敛。

6. 最后应得到哪些文件#

my-Al/保存自己的原生计算；my-Al/surface-results/workfunction.json给出左右真空、费米能和功函数；CSV保留完整平面平均势与密度。资料包的results/和三张汇总图则用于逐项对照本页12组参考条件。两者分开保存，避免把附件参考值当成自己新任务的结果。

需要继续读实空间数据，可看ELF、PDOS与电荷和局域占据轨道。势定义可对照2026.2势CUBE手册与官方功函数练习。