ABACUS 带电晶胞与门电极实操

下载完整输入、原生密度/电势/日志与绘图脚本，解压进入附件根目录。先读已有数据：

bash analyze.sh

它只用一个CPU核心分析参考数据，生成CSV和四张图，不启动DFT。需要完整复算时：

sbatch submit.slurm

脚本申请一张P100，顺序处理18个目录，其中3个是读取匹配SCF密度的静电势导出控制，1个是明确标识的门区电子聚集反例。研究自己的体系只需选用对应模型输入，不必重算所有对照；提交资源按自己的账户通过yskit核对。参考环境为ABACUS3.10.1、GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8、OpenBLAS0.3.30，每个计算1个MPI进程/1个CPU线程、GPU双精度。不要在登录节点执行DFT。

本批是固定Pt12薄层、Γ点、60 Ry和scf_thr=1e-6的操作示例；电子展宽0.015 Ry用于金属占据平滑，不能解释成溶液或离子温度。真空厚度、横向尺寸、k点和数值精度未作材料性质收敛，不输出精确功函数、恒电势反应能或溶剂化自由能。

1. 先数电子，再定义补偿方式#

本例Pt赝势每个原子含18个价电子，12个Pt的中性电子数是216。nelec 217表示增加一个电子，晶胞电子电荷变化为−1e；nelec 215则少一个电子、对应+1e。不要把电子数变化的正负号直接当成晶胞电荷的正负号。

周期带电晶胞和带补偿平面的栅极模型，边界条件及电势参考不同。它们都使用固定nelec，不等于恒电势算法；不能直接用两个模型的绝对总能量差解释电化学反应能或实验电极电位。

2. 温和充电的完整栅极输入#

INPUT_PARAMETERS
suffix tutorial
calculation scf
dft_functional pbe
basis_type pw
device gpu
precision double
pseudo_dir ../../PP_ORB
orbital_dir ../../PP_ORB
symmetry -1
ecutwfc 60
scf_thr 1e-6
scf_nmax 200
ks_solver cg
nbands 144
gamma_only 0
smearing_method gauss
smearing_sigma 0.015
mixing_type broyden
mixing_beta 0.1
cal_force 1
cal_stress 0
out_chg 1
out_pot 2
efield_dir 2
efield_pos_max 0.8
efield_pos_dec 0.1
nelec 216.1
gate_flag 1
zgate 0.7
block 1
block_down 0.65
block_up 0.75
block_height 1.0

ATOMIC_SPECIES
Pt 195.078 Pt_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.889716

LATTICE_VECTORS
5.54844 0 0
-2.77422 4.80509 0
0 0 20

ATOMIC_POSITIONS
Cartesian

Pt
0.0
12
0.0000000000 0.0000000000 5.0000000000 1 1 1
2.7742217000 0.0000000000 5.0000000000 1 1 1
-1.3871108000 2.4025464000 5.0000000000 1 1 1
1.3871110000 2.4025461000 5.0000000000 1 1 1
2.7742215000 1.6016975000 7.2651424000 1 1 1
0.0000000000 1.6016975000 7.2651424000 1 1 1
1.3871107000 4.0042440000 7.2651424000 1 1 1
-1.3871108000 4.0042440000 7.2651424000 1 1 1
0.0000000000 3.2033951000 9.5302848000 1 1 1
2.7742217000 3.2033951000 9.5302848000 1 1 1
1.3871108000 0.8008486000 9.5302848000 1 1 1
4.1613325000 0.8008486000 9.5302848000 1 1 1

K_POINTS
0
Gamma
1 1 1 0 0 0

本例nelec 216.1增加0.1个电子，gate_flag 1加入补偿平面，方向由efield_dir 2确定。zgate 0.7是单胞分数位置，不是0.7 Å。block_down/up 0.65/0.75定义门区势垒，block_height 1.0单位Ry，约13.606 eV；这个数值势垒不是某种真实绝缘体的带隙。

pt-uniform-mild使用相同216.1电子数、相同几何和k点，但不启用栅极；另有215和217电子的均匀周期带电对照。研究时需要同时核对电荷、补偿方式、真空、门区和电子密度，不能只改一个电子数。

3. 电子实际去了哪里#

本例保留一个清楚标识的反例pt-gate-plus：增加1个电子、势垒只有0.01 Ry。虽然它达到电子收敛，门区却聚集了大量电子，不能解释成金属表面获得了目标的1e电荷。温和充电/较高势垒组合是另一个模型，两者同时改变了充电量和势垒，不能当作严格的单参数对照。

门区z/L=0.65–0.75的实际积分：

中性模型也有电子尾部，所以不能把门区积分全部称为额外泄漏电子；应结合中性参照和完整密度差剖面。本例温和模型在更宽空区z/L=0.6–0.85相对中性的电子数变化约0.0016e，也一并检查，而不是只检查势垒中心。

4. 用ASE自己积分#

ASE读取cube后原子坐标/体积以Å表示，数据仍保持原生Bohr⁻³。下面在附件根目录直接读取参考文件：

module load ase/3.29.0
python - <<'PY'
from pathlib import Path
import numpy as np
from ase.io.cube import read_cube
from ase.units import Bohr
for name in ['pt-vacuum','pt-gate-plus','pt-gate-mild']:
    p=Path('reference/cases')/name/'OUT.tutorial/SPIN1_CHG.cube'
    with p.open() as f:
        a=read_cube(f)
    rho=a['data']
    zfrac=np.arange(rho.shape[2])/rho.shape[2]
    mask=(zfrac>=0.65)&(zfrac<=0.75)
    voxel=a['atoms'].get_volume()/rho.size/Bohr**3
    print(name,'total',rho.sum()*voxel,'gate',rho[:,:,mask].sum()*voxel)
PY

本例全部cube的总电子数与目标相差小于0.003e，同时用原生带权占据数独立核对。小量充电尤其需要注意网格和文本输出精度，不要为了让图好看强行缩放密度差到目标电子数。

python - <<'PY'
import csv
with open('reference/analysis/conditions.csv') as f:
    for r in csv.DictReader(f):
        if r['name'] in ['pt-uniform-mild','pt-gate-mild','pt-gate-plus']:
            print(r['name'],r['electron_target'],r['cube_electrons'],
                  r['gate_region_electrons'],r['role'])
PY

5. 自己建模时的顺序#

先完成中性薄层并确认真空；再选择充电量和补偿方式；检查门平面及势垒没有穿过原子；计算后核对电子数和电子分布；最后才讨论表面响应。若出现门区聚集，先重新审视模型和势垒/真空设置，而不是把结果直接用于电容、吸附能或电极电位。

参数定义见官方补偿平面说明。