pymatgen弹性张量与反应能

下载输入、实测结果和完整脚本。在解压目录运行 bash run.sh，结果写入 work/；reference/保留原始结构与CSV，脚本不会自动提交新的电子结构计算。

1 从实际弹性常数构造张量#

module load pymatgen/2026.9.24
python - <<'ELASTIC'
import json
import numpy as np
from pymatgen.analysis.elasticity.elastic import ElasticTensor
r = json.load(open('inputs/elastic.json'))['relaxed']
c = np.zeros((6,6))
c[:3,:3] = r['C12']
np.fill_diagonal(c[:3,:3], r['C11'])
np.fill_diagonal(c[3:,3:], r['C44'])
t = ElasticTensor.from_voigt(c)
print('B/G，GPa:', t.k_vrh, t.g_vrh)
print('Young，GPa:', t.y_mod / 1e9)
print('Poisson:', t.homogeneous_poisson)
ELASTIC

本例Si数据来自集群既有VASP应力–应变计算。输入采用GPa；此接口的 k_vrh、g_vrh为GPa，y_mod返回Pa，因此显示GPa时除以10⁹。

2 固定离子与松弛离子要分开#

两组同为Si，但应变后是否允许内部原子位置松弛不同。不能混用两组C11/C12/C44拼出一张张量。上面按立方对称性构造矩阵的代码也不能直接用于任意低对称材料，后者需要实际完整张量。

Voigt次序为xx、yy、zz、yz、xz、xy；涉及剪切应变时应确认工程剪切与张量剪切的因子约定。先核对原始应力单位、正负号、平衡残余应力和应变区间，再解释拟合值。

3 先配平化学式#

python - <<'BALANCE'
from pymatgen.core import Composition
from pymatgen.analysis.reaction_calculator import Reaction
r = Reaction([Composition('H2'),Composition('O2')],
             [Composition('H2O')])
print(r)
BALANCE

输出 H2 + 0.5 O2 -> H2O。这一步只配平，没有给出水生成能；需要同一套可比能量后才能计算反应能。

4 计算真实Cu/Au反应能#

资料包的 run_workflows.py用同方法Cu、Au和CuAu总能建立 ComputedEntry，再调用：

from pymatgen.analysis.reaction_calculator import ComputedReaction
rx = ComputedReaction([cu_entry, au_entry], [cuau_entry])
print(rx)
print(rx.calculated_reaction_energy)

实际输出对应 Cu + Au -> CuAu，反应能−0.1129197775 eV，每个所打印的CuAu化学式含2个原子；除以2后为−0.05645988875 eV/atom，与凸包教程的CuAu形成能一致。

用户自己的输入可能含不同大小超胞；对象会按配平反应处理。请始终将打印的反应式、总能和归一化方式一起保存，避免把“每次反应”“每化学式”和“每原子”混为一谈。