VASP带电空位计算实操

本例回答一个具体问题：删除一个Si后，怎样分别计算q=−1、0、+1、+2的空位，并画出随费米能变化的形成能？操作顺序是准备完整晶体与空位模型→分别计算各电荷态→提取能量和原子位点电势→做eFNV修正→画形成能和转变能级。

1 先用实际结果跑通分析#

下载带电硅空位输入与实际结果，保存为charged-defects.zip：

python3 -m zipfile -e charged-defects.zip .
cd vasp-charged-defects
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

这组命令读取本集群已完成的15个计算：16、54、128位点的完整超胞及各四个空位电荷态。它运行原生pydefect电荷修正并重画4张图，不会重新提交DFT。数值在results/CORRECTIONS.json，尺寸比较在results/SIZE-COMPARISONS.json。

分析使用pymatgen环境，setup.sh把附带的pydefect 0.10.1、vise 0.9.6和PyYAML 6.0.3安装到示例自己的.deps，不修改共享环境。包内保留这些开源依赖的许可证。

2 准备结构：先明确优化哪些自由度#

128位点模型来自2原子Si原胞的4×4×4扩胞，删去原点Si后有127原子。完整晶体采用统一晶格常数5.468925912596983 Å。空位优化采用ISIF=2，只优化原子位置，不改变晶胞形状和体积；完整参考与所有缺陷电荷态使用同一晶胞。

下面生成与附件中128位点起始输入一致的POSCAR，便于看清模型构造：

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
from itertools import product
import hashlib
out = Path('generated-models'); out.mkdir(exist_ok=True)
a = 5.4689259125969834
lattice = [[0,2*a,2*a],[2*a,0,2*a],[2*a,2*a,0]]
positions = []
for t in product(range(4), repeat=3):
    for v in [(0,0,0), (.25,.25,.25)]:
        positions.append([(t[k]+v[k])/4 for k in range(3)])
for name, xyz, reference in [
    ('POSCAR-bulk128', positions, 'bulk128'),
    ('POSCAR-vacancy127', positions[1:], 'vacancy127-q0')]:
    text = 'Si primitive4x4x4\n1.0\n'
    text += '\n'.join(' '.join(f'{x:.15f}' for x in row)
                      for row in lattice)
    text += '\nSi\n' + str(len(xyz)) + '\nDirect\n'
    text += '\n'.join(' '.join(f'{x:.15f}' for x in row)
                      for row in xyz) + '\n'
    target = out/name; target.write_text(text)
    assert target.read_bytes() == Path(
        f'inputs/size128/{reference}/POSCAR').read_bytes()
    print(name, len(xyz), hashlib.sha256(target.read_bytes()).hexdigest())
PYCODE

附件保留原始输入和提取的优化后结构。图中把优化后的原子周期折回空位附近，叉号表示被删除的位点。

本例各电荷态采用一套共同初始构型和MAGMOM=127*0.5。这适合展示流程，尚未穷举不同初始磁态、局部对称破缺或Jahn–Teller畸变；自己的定量研究需要比较多个合理起点，不能仅因力收敛就认定找到全局最低能构型。

3 设置电荷与INCAR#

本例Si PAW每原子有4个价电子，因此127原子的中性电子数为508。正电荷表示移走电子，NELECT=508−q：

q=+1的真实输入如下：

SYSTEM = vacancy127-q+1
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
IBRION = 2
NSW = 160
ISIF = 2
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .TRUE.
NELM = 150
ALGO = Normal
ISYM = 0
LVHAR = .TRUE.
ISPIN = 2
MAGMOM = 127*0.5
NELECT = 507

EDIFF=1E-6控制电子收敛，EDIFFG=-0.02控制结构优化的力阈值。ISPIN=2进行共线自旋极化；ISYM=0关闭空间对称操作，但它本身不会自动给初始结构添加扰动。LVHAR=.TRUE.保存静电势相关输出，LCHARG=.TRUE.保存电荷密度，方便后续核对。高斯展宽下总磁矩不一定为整数，应同时检查占据、磁态与展宽敏感性。

128位点和127原子模型都用以下KPOINTS：

Gamma mesh
0
Gamma
2 2 2
0 0 0

16、54位点模型分别采用4×4×4与3×3×3网格；同一尺寸的完整晶体和缺陷必须使用一致的数值设置。不同尺寸的倒空间密度近似对应，仍有残余k点误差，详见尺寸与转变能级。

4 自己提交计算时怎么做#

先复制附件中某个实际输入到新的计算目录，再按元素顺序准备本人有权使用的Si PAW PBE 64版POTCAR。下面先准备q=+1输入并检查Slurm脚本语法：

mkdir -p my-vacancy-qplus1
cp inputs/size128/vacancy127-q+1/{INCAR,POSCAR,KPOINTS} my-vacancy-qplus1/
cp run-p100.slurm my-vacancy-qplus1/run.slurm
bash -n my-vacancy-qplus1/run.slurm
cat my-vacancy-qplus1/INCAR

准备POTCAR后，在登录节点切入该目录提交：

cd my-vacancy-qplus1
sbatch run.slurm
squeue -u "$USER"

完整超胞和其余三个电荷态分别建目录提交，不能把不同电荷态接在同一个旧OUTCAR目录中。附件同时提供P100与V100 16G脚本；更换平台时，队列、GPU资源名和VASP模块要一起对应。两种脚本均为1 MPI进程、1 CPU线程、1 GPU。本例并非CPU/GPU性能对比。

VASP为6.6.1 GPU版，NVHPC 25.3、CUDA 12.8、HPC-X 2.22.1、MKL 2023.2；所有模型统一PBE、PAW PBE 64、ENCUT=550 eV。POTCAR、WAVECAR、CHGCAR和完整VASP原生输出不随公开包分发。

5 从自己的输出提取数据#

extract.py实际读取vasprun.xml、OUTCAR与完整晶体的EIGENVAL；它检查电子收敛、最大力、原子数和晶胞，再提取最终能量、价带顶与逐原子位点电势。假设五个计算目录放在my-calculations/中：

source env.sh
"$YS_PY" extract.py   --bulk my-calculations/bulk128   --defects my-calculations/vacancy127-q0             my-calculations/vacancy127-q+1             my-calculations/vacancy127-q-1             my-calculations/vacancy127-q+2   --epsilon dielectric.json --output my-size128.json

这个读取脚本已对本例15份实际输出逐项运行。它目前专用于非磁性完整Si晶体+单个Si空位，不适用于任意多组分、磁性宿主、带隙闭合或一般复杂缺陷；更换材料要先调整价电子、缺陷构造、原子对应、介电张量和带边判定。

dielectric.json来自相同PBE、PAW、晶格与550 eV设置的Si原胞LEPSILON计算，6×6×6网格，电子介电张量对角元约16.41154。Si为非极性材料；这是本例计算参照，不是实验介电常数，也不代表介电参数已经充分收敛。

6 电荷修正与形成能#

本例使用正单位试探电荷的电势约定，等于VASP原子核位点静电势的负值。脚本把实际逐位点电势送入pydefect原生eFNV接口，并用独立的介电度量Ewald实现核对所有位点模型电势及点电荷项，误差不超过1e-6 V或eV。原生位点电势本身的打印精度约1e-4 V，不能把模型计算的更多小数当成DFT输入精度。

形成能取：

Ef(q, EF) = Edef(q) - (127/128)*Ebulk
            + q*(EVBM + EF) + Ecorr(q)

这里EF从完整晶体价带顶起算；Si化学势为同尺寸完整晶体每原子能量。统一使用最终e_0_energy，不混入展宽自由能。eFNV的Ecorr已经包含本例的模型电荷与远区对齐项，不要另外再加一次同样的电势对齐。

128位点结果如下：

修正不一定总为正，尤其对齐项可能改变总符号。继续阅读尺寸与转变能级，再判断哪些交点在带隙内且属于最低能包络。

原理与接口可参照VASP NELECT和pydefect项目。本例为三维周期体相缺陷，不把这一修正流程直接套到带电薄层或含真空的二维模型。