CP-VASP吸附前后的电子巨正则能差

下载输入文件、原生数值摘录和提取绘图脚本（ZIP）。

unzip cpvasp-reference-adsorption-wave253-*.zip
cd cpvasp-wave253
bash run.sh

bash run.sh 只读取附件中已完成的集群计算数值摘录并画图，不启动VASP、不提交新作业。输出在 analysis/ 和 figures/。如需自己计算，下面按步骤准备输入、提交任务，再提取自己的输出。附件中的 .selected.txt 是从真实原生结果逐行选取的内容，不是可供接续计算的完整OUTCAR。

1 要比较什么#

反应写成 Cu₄＋½H₂→Cu₄H。每个目标电子化学势下分别比较裸电极、top吸氢和fcc吸氢，再减去半个气相H₂的能量。裸电极与吸氢体系必须使用相同TARGETMU、表面晶胞、截断能和溶剂参数。

这里Cu₄和Cu₄H保持固定几何，H只在上表面；H₂单独优化。结果是这个模型下的吸附电子巨正则能差，没有加入零点能、熵、标准态或SHE/RHE参照，不能称作完整吸附自由能或自由能台阶图。

2 准备吸附模型并提交#

先按中性参照教程复制 inputs 为 work，计算中性电极并运行 python prepare_targets.py。吸附目录提供Cu₄H的完整输入，POSCAR为Cu、H顺序，需准备相同顺序的合法POTCAR。top和fcc的完整结构与收敛参数也见氢吸附位点教程。

cat work/cpvasp-H-top-mu+0.0/POSCAR
cat work/cpvasp-H-fcc-mu+0.0/INCAR
for name in work/cpvasp-H-*; do
  (cd "$name" && sbatch submit.sh)
done

与中性Cu₄相比，中性Cu₄H多一个H价电子，故其 N₀=45，输入初始 NELECT=45。真实最终电子数由目标电势决定。此处采用 INIT_ECHANGE=0.0005、ECHANGE_MAX=0.003、CAP_MAX=0.05 的较小更新幅度，最终误差仍按0.001 eV检查。

3 单独计算气相H₂参照#

H₂不加连续介质溶剂，也不开恒电势控制，2个电子、非自旋极化单重态。15 Å立方盒、Γ点；ISIF=2保持盒子固定，IBRION=2优化原子位置，受力标准0.02 eV/Å。

INCAR

SYSTEM = H2 reference for CP CuH adsorption
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
NELM = 160
ALGO = Normal
ISMEAR = 0
SIGMA = 0.1
ISPIN = 1
ISYM = 0
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
NCORE = 1
NBANDS = 8
IBRION = 2
POTIM = 0.2
NSW = 40
ISIF = 2
NELECT = 2
LCEP = .FALSE.
LSOL = .FALSE.

POSCAR

H2 isolated singlet
1
15 0 0
0 15 0
0 0 15
H
2
Cartesian
7.5 7.5 7.13
7.5 7.5 7.87

KPOINTS

Gamma molecule
0
Gamma
1 1 1
0 0 0

按H元素准备同一PAW_PBE H赝势，然后提交：

cd work/H2
sbatch submit.sh
grep "reached required accuracy" OUTCAR
grep "energy(sigma->0)" OUTCAR | tail -1

本次实际优化后H—H键长约 0.7501 Å，最大受力约 0.00012 eV/Å，E(H₂)=−6.77194887 eV。分子键长和气相能量与当前泛函、赝势和盒子设置对应，不把它作为实验基准值。

4 为什么不能只减两个总能#

当前补丁的输出采用：

n_excess = N - N0
GCE = E(sigma->0) - TARGETMU * n_excess
Δn_excess = (N_Cu4H - 45) - (N_Cu4 - 44)
ΔE = E_Cu4H - E_Cu4 - 0.5 E_H2
ΔΩ = ΔE - TARGETMU * Δn_excess

用过量电子数之差，不是直接用两个总电子数之差。 Cu₄H本来就比Cu₄多一个中性价电子；若错把这个H价电子也当作从电子库额外取来的电子，会产生很大的参照误差。也可用同一TARGETMU下 GCE(Cu₄H)−GCE(Cu₄)−½E(H₂) 交叉检查。不同程序或补丁版本若输出不同定义的巨正则量，应先核对定义。

以TARGETMU=−4.703 eV的top位点为例：

ΔE                    = 0.803935865 eV
Δn_excess             = 0.016003 electron
-TARGETMU × Δn_excess = 0.075262109 eV
ΔΩ                    = 0.879197974 eV/H

原生CPM行的打印精度有限，所以直接用原生GCE相减时最后几位可有微eV级差别；附件同时保存两种计算，不能把打印位数差当作新的物理效应。

5 提取、检查和画自己的数据#

回到附件根目录，11个计算目录均应保留完整OUTCAR和 stdout.log：

cd ../..
module load ase/3.29.0
python analyze.py --results work --output own-analysis
python plot.py --analysis own-analysis --structures work --output own-figures
cat own-analysis/adsorption.csv

真实示例结果为：

本例固定几何下top的ΔΩ低于fcc。两个结果均为正，仅表示按上述电子能量和固定模型定义的反应能差为正；不能据此直接判断实际电催化反应的速率、平衡覆盖度或实验过电位。若研究吸附自由能，应在相同电势条件下优化结构，继续补零点能、振动和熵等修正；若画自由能台阶图，还需准备全部反应中间体与一致参照。

裸Cu₄电容使用两面归一，但这里H只在上表面，ΔΩ按每个吸附H报告，不再除以两个表面。更大表面、覆盖度、显式水、偶极影响和不同溶剂参数需要单独检查。

参数、短MD与原子位置优化#

继续参考14个补丁参数与电子数方案、短程MD与溶剂切换和原子优化及溶剂版本。这些实例区分每步收敛与动态电子数波动、离子和电子温控、晶胞固定下的原子优化，并给出6.4.2经典溶剂E-fermi+FERMI_SHIFT参照与实际终端提取方式。当前动态4/5方案不直接使用原生GCE列做反应能比较；参数及版本适用范围见对应教程。