pymatgen：从分子与碎片能量计算断键能

要计算一根键的断裂能，先明确断裂后得到什么碎片，以及碎片的电荷和自旋。本例计算气相中性均裂：

H₂O → OH（双重态）＋ H（双重态）。

在鸢算科技集群用xTB6.7.1、GFN2实际优化水和OH，计算H原子单点，并对优化后的水和OH计算Hessian。再用pymatgen2026.9.24的原生BondDissociationEnergies匹配断键与碎片，独立整理正号的电子解离能Dₑ和含零点能修正的D₀。

这是展示操作流程的GFN2近似结果。自由基能量和断键能对方法、自旋态与几何优化敏感；本例数值不作为高精度量化结果或实验水解离能的替代。

下载完整结构、能量与频率CSV、分析、计算及绘图脚本。

unzip pymatgen-bond-energy-wave714-*.zip
cd pymatgen-bond-energy
bash results.sh
bash plot.sh

1. 准备个人目录的后处理依赖#

module load pymatgen/2026.9.24
bash prepare.sh
bash run.sh
cat analysis/Water-OH-electronic-cleavage-energy.csv
cat analysis/zero-point-result.json

prepare.sh把固定版本、固定哈希的社区OpenBabel预编译包放到当前目录的toolchain/python，运行库为OpenBabel3.1.0；公共Python环境不变。Python3.11、NumPy、SciPy、pymatgen及Matplotlib使用集群现有环境。

run.sh读取包内已经实际计算的能量和零点能，完成原生碎片分析及单位换算，不再运行量化计算。首次执行在analysis生成结果；保留已有目录后再分析时，请在另一份解压目录操作。

2. 认识三个输入结构与电子态#

水：电荷0、单重态#

inputs/water.xyz，坐标单位Å：

3
H2 O1
O 0.006958 0.008952 -0.000000
H 0.965819 -0.013851 0.000000
H -0.255564 0.931526 0.000000

OH：电荷0、双重态#

inputs/hydroxyl.xyz，坐标单位Å：

2
H1 O1
O 0.003326 -0.000000 -0.000000
H 0.966674 0.000000 -0.000000

H：电荷0、双重态#

inputs/hydrogen.xyz，坐标单位Å：

1
H1
H 0.000000 0.000000 0.000000

水和OH的几何分别优化，不能只从水的坐标删除一个H，就把未松弛的OH当作绝热解离产物。本例H原子只有单点，不做没有内部自由度的几何优化。

xTB参数--uhf表示未成对电子数：水用0，OH和H用1；它不是自旋多重度。转换到pymatgen分子对象时，则分别设置多重度1、2、2。

3. 如果要重新计算这三个体系#

包内保留已优化XYZ作为起始结构。重新计算时先确认CPU分区可用；用Slurm提交：

sbatch calculate.slurm
squeue -u "$USER"

作业依次运行calculate.sh、collect-new.py和analyze-new.sh，新量化输出在work-native，新分子条目在inputs-new，新分析在analysis-new。原参考数据不覆盖。脚本固定单线程、顺序执行；若CPU分区未开放，等待可用CPU资源后提交。

水的电子计算命令为：

module load xtb/6.7.1
export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
xtb input.xyz --gfn 2 --acc 0.1 --chrg 0 --uhf 0 --opt tight --json > native.out 2>&1

OH使用--uhf 1 --opt tight；H使用--uhf 1，去掉--opt tight。这些命令在三个独立工作目录中运行，不能共用同名临时文件。

检查电子计算，再把优化后的坐标用于频率计算：

rg 'TOTAL ENERGY|GEOMETRY OPTIMIZATION CONVERGED|normal termination' work-native/*/native.out
cat inputs-new/molecule-entries.json

水的Hessian使用：

xtb input.xyz --gfn 2 --acc 0.1 --chrg 0 --uhf 0 --hess --json > native.out 2>&1

这里的input.xyz应为水优化得到的xtbopt.xyz；OH相应改为--uhf 1并使用OH优化结构。calculate.sh已经按这一顺序复制文件，日常直接提交脚本即可。

4. 读取同一种能量#

cat reference-data/Water-OH-H-GFN2-energies.csv

这三项均读取xTB输出的TOTAL ENERGY，包含同一GFN2模型的相应能量项。不要把JSON中名为electronic energy的另一项混进此表，也不要混用不同方法、溶剂模型或不同电荷参考的能量。

5. 原生断键接口与正负号#

bash python.sh - <<'PY'
import json
from pymatgen.analysis.bond_dissociation import BondDissociationEnergies
with open('inputs/molecule-entries.json') as f:
    entries = json.load(f)
analysis = BondDissociationEnergies(
    entries['water'], [entries['hydroxyl'], entries['hydrogen']]
)
for row in analysis.bond_dissociation_energies:
    print(row)
PY

每个条目包含formula_pretty、initial_molecule、final_molecule、final_energy和smiles。本例明确约定final_energy为Hartree；接口不会替你把输入单位自动换成eV。

本版本原生结果的首个能量数值为母分子能量减碎片能量之和，本例−0.248686681978 Hartree。若要展示正的断键代价，应按下式计算：

Dₑ = E(OH) + E(H) − E(H₂O) = 6.767109 eV = 652.926794 kJ/mol。

run.py同时检查原生符号与独立差值，避免只给输出贴上“断键能”就忽略方向。方法更换后，各碎片仍须单独优化，选好相应电子态。

本例只提供中性OH和H。原生接口还会寻找其他电荷分离组合，因此会提示缺少离子碎片；提示保存在analysis/native-notes.txt。这里没有计算OH⁻＋H⁺或OH⁺＋H⁻，不能据此判断异裂通道。

6. 加入零点振动能#

cat reference-data/Water-OH-positive-vibrational-modes.csv
cat reference-data/Water-OH-zero-point-cleavage-energy.csv

水的三个内部频率为1539.26、3643.59、3651.90 cm⁻¹，OH的一个内部频率为3575.88 cm⁻¹。水的6个、OH的5个刚体平移/转动模不计入内部振动；本例没有显著的虚频。真实新体系若有内部虚频，先检查几何是否为适用的局部极小值。

D₀ = Dₑ + ZPE(OH) + ZPE(H) − ZPE(H₂O) = 6.441101 eV = 621.471729 kJ/mol。

H原子没有分子内部振动，零点能取0。图中两条含零点能的绿色能级共用同一电子能量零点，D₀是两者的间距；不能把D₀与电子母分子能量之间的高度混为一谈。

D₀仍不是298K断键焓或反应自由能。本例未加入平移、转动、温度、压力、熵或溶剂修正。需要气相分子自由能时可参考分子RRHO与准RRHO，并保持每个物种的标准态一致。

7. 换成自己的断键问题#

先写清守恒的反应式、每个物种电荷与多重度；以相同方法分别优化母分子和碎片，检查收敛和振动，再收集同一种总能量。检查成键图是否对应准备断开的键，最后计算能量差和所需热力学校正。

环体系断一根键可能只得到一个开环分子；多重断键、带电断裂和不同自旋通道需要相应产物及能量，不能用本例三个条目代替。生成几何碎片见分子碎片与开环；碎片图本身不提供量化能量。

接口参考pymatgen断键分析；量化参数见xTB官方手册。