原子智能

Appearance

第 06 章：现实引力：热力学凸包稳定性与动力学声子谱检验 ​

“在虚拟计算机中，你可以轻松创造出一颗结构令人赞叹的奇幻晶体；
但只要它的自由能位于热力学凸包之上哪怕只有 100 meV/atom，或者在声子谱边缘有一丝微弱的虚频，
自然界的原子就会在飞秒级内自发崩塌，化为一堆毫无用处的粉末。”

1. 引言：算法幻觉与物理现实的残酷碰撞 ​

在近年的 AI 生成晶体热潮中，许多研究论文常常自豪地宣称：“我们的生成模型设计出了 100,000 个全新的晶体材料，它们的 DFT 形成能（Formation Energy ΔHf\Delta H_f）均为负数，因此全部是热力学稳定的！”

在真实的凝聚态物理学家与无机固体化学家眼中，这种论调是极为业余且致命的物理误区！

形成能为负数（ΔHf<0\Delta H_f < 0），仅仅代表该化合物相对于其由单质元素构成的气相或单质基态是有利反应。它丝毫不能保证这个晶体在现实中能够稳定存在！

  • 如果化合物 A2BC3A_2 B C_3 相对于其单质分解是放热的（ΔHf=−1.2 eV/atom\Delta H_f = -1.2 \text{ eV/atom}）；
  • 但它自发分解为已知的稳定二元化合物 2AB+C32AB + C_3 时，能够释放出更多的热量（产生进一步的自由能骤降）；
  • 那么在真实的烧结炉中，A2BC3A_2 B C_3 将永远不会诞生，它会以 100% 的概率分解沉淀为 ABAB 和 C3C_3！

要评判一个新生成的固态材料究竟是“真实的全新基态”还是“注定分解的幻觉假象”，必须经受两道最严苛的物理大审判：

  1. 宏观热力学稳定性（Thermodynamic Stability）：基于多元相图的能量凸包（Convex Hull）分析；
  2. 微观动力学稳定性（Kinetic/Dynamical Stability）：基于全布里渊区声子色散谱（Phonon Dispersion）的无虚频检验。

2. 热力学稳定性与能量凸包（Convex Hull） ​

2.1 吉布斯相平衡与能量凸包几何定义 ​

在一个多元化学反应体系中，自然界在恒温恒压下永远追求体系**总吉布斯自由能（Gibbs Free Energy）**的绝对极小化：

G=H−TS=EDFT+PV−TSG = H - TS = E_{\text{DFT}} + P V - T S

在零温零压的基态量子从头算中，通常以 DFT 总能量 EE 逼近自由能。

设多组分化学体系包含元素组分向量 c=(c1,…,cK)T\mathbf{c} = (c_1, \dots, c_K)^T，满足 ∑k=1Kck=1\sum_{k=1}^K c_k = 1。 假设数据库中已知存在 MM 个热力学稳定或亚稳态相 {p1,p2,…,pM}\{ p_1, p_2, \dots, p_M \}，各相的化学计量组分为 ci\mathbf{c}_i，单位原子形成能为 E(pi)E(p_i)。

在几何上，所有可能存在的相在 (c,E)(\mathbf{c}, E) 空间中构成散点云。能量凸包（Energy Convex Hull）被严格定义为包络所有这些已知相点的下凸多面体包络面（Lower Convex Hull）！

图 6-1：无机材料热力学能量凸包（Convex Hull）与多相热力学竞争分解机制示意图

图 6-2：热力学能量凸包（Convex Hull）相图、分解驱动力（ΔE_hull）与可合成亚稳态区间定量验证

2.2 凸包距（Energy Above Convex Hull, ΔEhull\Delta E_{\text{hull}}）的严密数学形式化 ​

对于任意一个通过生成模型预测的新晶体候选 xx，其化学组分为 cx\mathbf{c}_x，能量为 E(x)E(x)。 它在已知相图竞争体系中的凸包距（ΔEhull\Delta E_{\text{hull}}），被严格定义为：该晶体能量与由现有竞争相任意线性组合生成相同组分所能达到的理论最低能量之差！

这在数学上是一个严格的线性规划（Linear Programming, LP）对偶问题：

ΔEhull(x)=E(x)−min⁡α∑i=1MαiE(pi)\Delta E_{\text{hull}}(x) = E(x) - \min_{\boldsymbol{\alpha}} \sum_{i=1}^M \alpha_i E(p_i)

约束条件（组分守恒与质量守恒）：

∑i=1Mαici=cx,∑i=1Mαi=1,αi≥0  (∀i∈{1,…,M})\sum_{i=1}^M \alpha_i \mathbf{c}_i = \mathbf{c}_x, \quad \sum_{i=1}^M \alpha_i = 1, \quad \alpha_i \ge 0 \; (\forall i \in \{1,\dots,M\})

凸包距的物理判决准则： ​

  • ΔEhull=0 meV/atom\Delta E_{\text{hull}} = 0 \text{ meV/atom}： 新晶体落在凸包表面或凸包顶点，它比该组分下人类已知的任何分解产物都要稳定，构成自然界真正的基态绝对稳定新相（Thermodynamically Ground State）！
  • 0<ΔEhull≤25 meV/atom0 < \Delta E_{\text{hull}} \le 25 \text{ meV/atom}： 新晶体具有极其微弱的热力学分解驱动力。然而，由于常温下的热运动涨落能量为 kBTroom≈25.7 meVk_B T_{\text{room}} \approx 25.7 \text{ meV}，且固态重构相变通常伴随着巨大的动力学势垒（Activation Barrier），这类材料在现实中极其容易作为**长寿命亚稳态材料（Metastable Materials）**被成功合成（如金刚石相对于石墨即为亚稳态）；
  • 25<ΔEhull≤70 meV/atom25 < \Delta E_{\text{hull}} \le 70 \text{ meV/atom}： 处于亚稳态上限，需通过特殊实验外场（如超高压高温烧结、外延薄膜晶格失配应力诱导）方可合成；
  • ΔEhull>100 meV/atom\Delta E_{\text{hull}} > 100 \text{ meV/atom}： 热力学分解驱动力过大，在合成过程中几乎瞬间自发相分离为其他已知二元/多元晶相，属于毫无现实价值的“数字幻觉”。

3. 动力学稳定性：全布里渊区声子色散谱与虚频检验 ​

热力学凸包仅仅排除了化学相分离的风险。但在微观晶格层面，晶体还必须抵御晶格微观振动的自发相变崩溃。这一关卡由晶格动力学（Lattice Dynamics）与声子谱守卫。

3.1 简谐近似与动力学矩阵（Dynamical Matrix） ​

在玻恩-奥本海默近似下，晶体势能面展开为原子小位移 uli\mathbf{u}_{li} 的泰勒级数。在简谐近似（Harmonic Approximation）下，二阶原子力常数张量（Interatomic Force Constants, IFC）为：

Φαβ(li,l′j)=∂2E∂uliα∂ul′jβ\Phi_{\alpha\beta}(li, l'j) = \frac{\partial^2 E}{\partial u_{li\alpha} \partial u_{l'j\beta}}

引入平面波声子本征解，在倒空间波矢 q\mathbf{q} 处构建动力学矩阵 D(q)D(\mathbf{q})：

Diα,jβ(q)=1MiMj∑l′Φαβ(0i,l′j)exp⁡(iq⋅(Rl′−R0))D_{i\alpha, j\beta}(\mathbf{q}) = \frac{1}{\sqrt{M_i M_j}} \sum_{l'} \Phi_{\alpha\beta}(0i, l'j) \exp\left( i \mathbf{q} \cdot (\mathbf{R}_{l'} - \mathbf{R}_0) \right)

求解本征方程得到声子色散关系的平方本征频率：

D(q)eqν=ωqν2eqν,ν=1,…,3ND(\mathbf{q}) \mathbf{e}_{\mathbf{q}\nu} = \omega_{\mathbf{q}\nu}^2 \mathbf{e}_{\mathbf{q}\nu}, \quad \nu = 1, \dots, 3N

3.2 虚频的物理本质与软模相变（Soft Mode Collapse） ​

若动力学矩阵的某一阶本征值 ωqν2<0\omega_{\mathbf{q}\nu}^2 < 0，其频率为纯虚数：

ω=i∣ω∣⟹u(t)∝exp⁡(−iωt)=exp⁡(∣ω∣t)\omega = i |\omega| \quad \Longrightarrow \quad u(t) \propto \exp(-i \omega t) = \exp( |\omega| t )

  • 在真实的物理世界中，这意味着微观原子在偏离平衡位置时，不仅感受不到拉回原位的弹性恢复力，反而受到顺向加速的失稳排斥力！
  • 极微小的热涨落位移都会随着时间指数发散，诱发晶格自发发生软模畸变或宏观晶体崩塌崩溃。
  • 检验标准：只有通过密度泛函微扰理论（DFPT）扫描晶体全布里渊区高对称路径（如 Γ→X→M→R→Γ\Gamma \to X \to M \to R \to \Gamma），确认整条声子色散曲线处处非负（ωqν≥0\omega_{\mathbf{q}\nu} \ge 0），该材料才具备在现实世界中稳定站立的物理资格！

4. 机械力学与弹性稳定性：Born 稳定性判据 ​

除了微观振动，晶体在宏观弹性形变下也必须满足热力学正定性。对于任意微小的应变张量 ε\boldsymbol{\varepsilon}，体系的形变应变能必须严格为正：

U=12∑i,j=16Cijεiεj>0U = \frac{1}{2} \sum_{i,j=1}^6 C_{ij} \varepsilon_i \varepsilon_j > 0

这要求 6×66 \times 6 的弹性刚度矩阵 CC 必须是正定矩阵。例如，对于最常见的立方晶系（Cubic Crystals），著名的 Born 稳定性判据 为：

C11−C12>0,C11+2C12>0,C44>0C_{11} - C_{12} > 0, \quad C_{11} + 2C_{12} > 0, \quad C_{44} > 0

若生成模型给出的晶胞弹性常数违背上述不等式，晶体在受到哪怕一微克的机械剪切力时都会瞬间碎裂。

5. 工业级代码实战：线性规划求解能量凸包距 ΔEhull\Delta E_{\text{hull}} ​

下面我们使用 Python 与 scipy.optimize.linprog，完全独立实现计算材料学核心的相图分解与能量凸包距（ΔEhull\Delta E_{\text{hull}}）线性规划求解器：

python
"""
文件名: convex_hull_stability_solver.py
功能: 基于线性规划 (Linear Programming) 实现多组分化学相图的能量凸包距计算，
      定量评估新生成材料的热力学绝对稳定性与亚稳态分解驱动力。
"""

import numpy as np
from scipy.optimize import linprog


class PhaseDiagramConvexHull:
    """
    热力学相图凸包 (Convex Hull) 求解引擎
    """
    def __init__(self, elements: list):
        self.elements = elements
        self.num_elements = len(elements)
        self.known_phases = []

    def add_known_phase(self, name: str, composition_dict: dict, energy_per_atom: float):
        """
        向历史已知数据库中录入稳定/竞争晶相
        composition_dict: 例如 {"Li": 1, "Fe": 1, "O": 2}
        energy_per_atom: 归一化单原子形成能 (eV/atom)
        """
        # 将化学计量比归一化为组分分数向量
        total_atoms = sum(composition_dict.values())
        comp_vector = np.zeros(self.num_elements)
        for i, el in enumerate(self.elements):
            comp_vector[i] = composition_dict.get(el, 0.0) / total_atoms

        self.known_phases.append({
            "name": name,
            "comp": comp_vector,
            "energy": energy_per_atom
        })

    def calculate_energy_above_hull(self, cand_comp_dict: dict, cand_energy_per_atom: float) -> dict:
        """
        利用对偶单纯形线性规划求解候选晶体的能量凸包距 ΔE_hull
        """
        total_atoms = sum(cand_comp_dict.values())
        c_target = np.zeros(self.num_elements)
        for i, el in enumerate(self.elements):
            c_target[i] = cand_comp_dict.get(el, 0.0) / total_atoms

        M = len(self.known_phases)
        # 目标函数系数: 极小化基底各相能量的线性组合 c^T * x
        c_obj = np.array([p["energy"] for p in self.known_phases])

        # 等式约束矩阵 A_eq: 组分原子守恒 + 组分和归一化为 1
        # A_eq @ alpha = b_eq
        A_eq = np.zeros((self.num_elements + 1, M))
        for j, p in enumerate(self.known_phases):
            A_eq[:self.num_elements, j] = p["comp"]
            A_eq[self.num_elements, j] = 1.0

        b_eq = np.zeros(self.num_elements + 1)
        b_eq[:self.num_elements] = c_target
        b_eq[self.num_elements] = 1.0

        # 变量边界: 各相摩尔混合分数 alpha_i >= 0
        bounds = [(0, None) for _ in range(M)]

        # 线性规划求解
        res = linprog(c_obj, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs")

        if not res.success:
            raise RuntimeError("线性规划求解失败，可能目标组分超出了相图边界！")

        hull_energy_at_comp = res.fun
        delta_e_hull = cand_energy_per_atom - hull_energy_at_comp

        # 提取优势分解产物产率
        decomposition_products = []
        for j, alpha in enumerate(res.x):
            if alpha > 1e-4:
                decomposition_products.append((self.known_phases[j]["name"], alpha))

        return {
            "delta_e_hull_eV": delta_e_hull,
            "delta_e_hull_meV": delta_e_hull * 1000.0,
            "is_ground_state": bool(delta_e_hull <= 1e-4),
            "is_synthesizable_metastable": bool(delta_e_hull <= 0.030),  # < 30 meV/atom
            "competing_hull_energy": hull_energy_at_comp,
            "decomposition_pathway": decomposition_products
        }


# =====================================================================
# 单元验证模块: 模拟真实三元锂电池固态电解质体系相图稳定性检验
# =====================================================================
if __name__ == "__main__":
    print(">>> [STEP 1] 构建 Li-La-Zr-O (LLZO 固态电解质母系) 热力学相图环境...")
    elements = ["Li", "La", "Zr", "O"]
    pd = PhaseDiagramConvexHull(elements)

    # 录入单质与常见基态二元/三元已知稳定相 (DFT形成能数据)
    pd.add_known_phase("Li_metal", {"Li": 1}, 0.0)
    pd.add_known_phase("La_metal", {"La": 1}, 0.0)
    pd.add_known_phase("Zr_metal", {"Zr": 1}, 0.0)
    pd.add_known_phase("O2_gas",    {"O": 2},  0.0)

    pd.add_known_phase("Li2O",     {"Li": 2, "O": 1}, -2.05)
    pd.add_known_phase("La2O3",    {"La": 2, "O": 3}, -3.85)
    pd.add_known_phase("ZrO2",     {"Zr": 1, "O": 2}, -3.78)
    pd.add_known_phase("Li6Zr2O7", {"Li": 6, "Zr": 2, "O": 7}, -3.12)
    pd.add_known_phase("La2Zr2O7", {"La": 2, "Zr": 2, "O": 7}, -3.92)

    print(">>> [STEP 2] 检验一个通过生成模型预测的石榴石结构 Li7La3Zr2O12 (LLZO) 候选物...")
    # 测试案例 1: 完美基态结晶 (能量极低, 应落在凸包上)
    candidate_comp = {"Li": 7, "La": 3, "Zr": 2, "O": 12}
    result_stable = pd.calculate_energy_above_hull(candidate_comp, cand_energy_per_atom=-3.45)

    print(f"    候选 1 (优化构型): 凸包距 ΔE_hull = {result_stable['delta_e_hull_meV']:.2f} meV/atom")
    print(f"    - 是否热力学基态: {result_stable['is_ground_state']}")
    print(f"    - 是否具备高可合成性: {result_stable['is_synthesizable_metastable']}")
    print(f"    - 理论竞争分解路径: {result_stable['decomposition_pathway']}")

    # 测试案例 2: 虚假高能量结构 (能量偏高, 应自发分解)
    print("\n>>> [STEP 3] 检验一个能量偏高、结构存在局域畸变的次级生成候选...")
    result_unstable = pd.calculate_energy_above_hull(candidate_comp, cand_energy_per_atom=-3.25)
    print(f"    候选 2 (畸变构型): 凸包距 ΔE_hull = {result_unstable['delta_e_hull_meV']:.2f} meV/atom")
    print(f"    - 是否热力学基态: {result_unstable['is_ground_state']}")
    print(f"    - 是否具备高可合成性: {result_unstable['is_synthesizable_metastable']}")
    print(f"    - 理论分解产物组合: {result_unstable['decomposition_pathway']}")

    assert result_stable['delta_e_hull_meV'] < result_unstable['delta_e_hull_meV'], "凸包距排序物理逻辑颠倒！"
    assert result_unstable['delta_e_hull_meV'] > 100.0, "不稳定相未被正确识别出高凸包驱动力！"
    print("\n>>> [SUCCESS] 热力学能量凸包线性规划验证全部通过！")

6. 本章小结 ​

本章为无机固态材料的逆向生成建立了不容逾越的“物理质量质检站”：

  1. 破除形成能迷信：深刻论证了为何负形成能无法证明稳定性，确立了以**能量凸包距（ΔEhull\Delta E_{\text{hull}}）**作为固态化学合成可行性的黄金标尺；
  2. 量化亚稳态生存边界：明确了 ΔEhull<25∼50 meV/atom\Delta E_{\text{hull}} < 25 \sim 50 \text{ meV/atom} 的热力学生存红线；
  3. 微观动力学与宏观力学双重锁死：解析了全布里渊区声子色散谱虚频（ω2<0\omega^2 < 0）与 Born 弹性常数正定性的决定性意义；
  4. 交付工业级相图 LP 算法：亲手构建了高维组分相图分解与凸包距求解器。

在完成了无机固态晶体的物理大审判之后，有机小分子领域又将迎来怎样的“现实引力”？如何确保生成的分子能被湿法化学家在实验室中合成出来？

进入 第 07 章：现实引力（分子可合成性量化与计算机辅助逆合成规划 CASP）！

《原子智能》· 纸质出版预备版 · PolyAI Team 著

Copyright © 2024-2026 PolyAI Team. 保留所有权利。