ASE电位与pH稳定性图实操

想判断材料在什么电位和pH下更倾向保留固相、溶解还是变成其他相，可以先学习Pourbaix电位–pH图。本例在鸢算科技集群实际调用 ASE3.29.0，使用其文献参考形成能数据库计算ZnO的7条竞争反应、4组设置共 90,804个网格点。这是热力学工具实例，不是新做的DFT，也不计算腐蚀速率或钝化膜生长。

1 下载后先看结果，再完整运行#

下载脚本、参考数据、全网格CSV和三张高清图。使用独立Python3.11环境，ASE3.29.0及其NumPy、SciPy、Matplotlib依赖，单进程单线程；不需要GPU。

python3 -m zipfile -e ase-ZnO-Pourbaix-wave312-*.zip .
cd ase-ZnO-Pourbaix-wave312
bash preview.sh

preview/results/保存从参考数值重新提取的CSV，preview/plots/生成说明图。需要从形成能重新计算整张图，再执行：

bash run.sh
cat work/results/VALIDATION.json
ls work/plots
head -n 6 work/results/reaction-coefficients.csv

run.sh依次运行原生Pourbaix网格、反应自由能提取、作图；输出到新的work/，不会覆盖reference/。已有同名工作目录时换一个新解压目录，保留自己的旧结果。

2 第一步：先准备参考形成能#

本例原始参考数据如下，单位是 eV/所写化学式，不是eV/原子：

其中水溶物种来自已安装ASE的solvated('Zn')数据库，Zn、ZnO和ZnO2(aq)补充值采用ASE官方同体系示例。数据库文献出处为Johnson等1992及Pourbaix1966，完整数值以下载包reference/native/DATA.json为准。

读取数据库和补充固相的终端操作：

module purge
module load ase/3.29.0
python - <<'PYCODE'
from ase.phasediagram import solvated
refs = dict(solvated('Zn'))
refs.update({'Zn': 0.0, 'ZnO': -3.323, 'ZnO2(aq)': -2.921})
for name, energy in refs.items():
    print(f'{name:14s} {energy:.12f} eV/formula')
PYCODE

当前原生反应平衡中，水储库实际采用ase.pourbaix.PREDEF_ENERGIES['H2O']=-2.4583 eV；这与旧溶剂数据库中的H2O(aq)数值有约11微eV差别。脚本按真实程序储库计算，不暗中用另一份水能量替换。

3 第二步：设温度、活度与参比电极#

当前ASE参数conc对所有标记(aq)的物种使用同一个活度，包括中性水溶物种，并非只改Zn²⁺浓度。稀理想溶液中可将活度近似理解为c/c°；这不是同时求解真实溶液中各物种浓度与物料守恒的物种分布计算。

四组文件的电位数值轴相同，参比不同。CSV用reference和potential_V两列注明；不要将RHE或AgCl列称为SHE。数据包原始NPZ保留历史键名U_SHE_V作为共用数值网格，解释具体组时须同时读取DATA中的参比。

4 第三步：生成图，并检查一个具体条件#

下面是可直接复制到终端的最小原生示例：

python - <<'PYCODE'
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from ase.phasediagram import solvated
from ase.pourbaix import Pourbaix
refs = dict(solvated('Zn'))
refs.update({'Zn': 0.0, 'ZnO': -3.323, 'ZnO2(aq)': -2.921})
pb = Pourbaix('ZnO', refs, T=298.15, conc=1e-6, reference='SHE')
for pH in [0, 10, 14]:
    energy, reaction = pb.get_pourbaix_energy(U=0.0, pH=pH)
    stable = 'ZnO(s)' if energy <= 0 else reaction.equation()
    print(pH, energy, 'eV/ZnO;', stable)
diagram = pb.diagram(U=np.linspace(-2, 2, 161), pH=np.linspace(0, 14, 141))
fig, ax = plt.subplots(figsize=(12, 7))
diagram.plot(normalize=False, include_water=True, ax=ax)
fig.savefig('my-ZnO-pourbaix.png', dpi=220, bbox_inches='tight')
PYCODE

normalize=False保留eV/式量ZnO。默认归一化会除以ZnO的2个原子，读色条前先分清单位。

本例SHE、活度10⁻⁶的三个实际点：

这里能量 ≤0表示目标ZnO稳定，>0表示存在能量更低的竞争路线。负值不裁成零，也不把它当作扩散势垒或腐蚀速率。

绿区是ZnO(s)，其他色块对应图例物种。两条虚线只是理想H⁺/H₂和O₂/H₂O水窗口，不包含过电位。细边界来自原生反应平面相交，原生端点约保留到0.001；底层CSV保留完整网格值。参考物种有限，因此这些边界不代表已经纳入所有可能的真实化合物。

5 第四步：读反应式，理解浓度和电位为什么改变边界#

python - <<'PYCODE'
from ase.phasediagram import solvated
from ase.pourbaix import Pourbaix
refs = dict(solvated('Zn'))
refs.update({'Zn': 0.0, 'ZnO': -3.323, 'ZnO2(aq)': -2.921})
pb = Pourbaix('ZnO', refs, conc=1e-6)
print('\n'.join(pb.get_equations()))
print('只看生成Zn金属的路线：')
print('\n'.join(pb.get_equations(contains='Zn')))
PYCODE

本例有ZnO向HZnO₂⁻、ZnO₂²⁻、ZnOH⁺、ZnO(aq)、Zn²⁺、Zn和ZnO₂(aq)的7条路线。ASE会补平H₂O、H⁺和电子，不能只比较两种化合物的裸形成能。

例如酸溶解ZnO + 2H+ → Zn++(aq) + H2O，其反应自由能在活度10⁻⁶时为：

ΔG = −1.016672434929 + 0.118318659209 × pH eV/ZnO。

单条路线的ΔG=0位于 pH=8.592664。活度提高到10⁻³，该零点变为 7.092664，相差−1.5。这个零点只是这一条反应的平衡条件；还有ZnOH⁺等路线竞争，不能直接把8.592664画成整张图ZnO的稳定区边界。

另一条还原路线ZnO + 2H+ + 2e- → Zn + H2O给出：

U_SHE = −0.432350 − 0.059159329605 × pH V。

本例参比转换是U_SHE=U_RHE−0.059159329605×pH、U_SHE=U_AgCl+0.222。因此同一还原条件在RHE轴上是−0.432350 V水平线。这里AgCl固定偏移来自当前模型，实际电极的温度、溶液和标定条件需要另行核对。

6 想换成自己的材料，改哪里#

先整理目标材料及竞争固相/水溶物种的同一热力学基准形成自由能，再替换refs和Pourbaix('目标式', refs, ...)。只加一个目标化合物、不提供竞争相，不能得到完整稳定区。DFT总能直接与实验水溶离子形成自由能混用通常缺少一致的参考、溶剂和热修正，需要先建立统一方案。

head -n 8 work/results/grid-values.csv
head -n 10 work/results/reaction-coefficients.csv
ls work/native

完整包保留每组全网格、每条反应常数/电位/pH系数和原生参考输出，可继续用自己的脚本画图。

原生图色条范围为±1 eV，超出部分只发生颜色饱和，完整CSV没有裁值。以上是有限参考集合、理想活度和固定温度下的操作参考；判断真实耐蚀性还要考虑动力学、膜层、溶液组成及更多竞争相。固态成分凸包请看ASE候选合金凸包，结构对齐问题看结构等价比较。