第 05 章：无机晶体生成实战：从 CDVAE 到 MatterGen 与 GNoME

“在有机分子中，原子通过共价键手拉手围成小巧的分子孤岛；
而在无机固态晶体中，原子被晶胞矩阵以空间群的严苛秩序，在三维无限空间中复制成浩瀚的晶格星阵。”

1. 引言：有机小分子与固态晶体的物理鸿沟

许多在计算机视觉或自然语言处理领域功成名就的深度学习学者，初次迈入晶体材料设计领域时，往往会经历巨大的挫败。他们试图把训练分子图（如生成 SMILES 字符串或二维分子的邻接矩阵）的成熟模型直接照搬到钛酸钡、氮化镓或高熵合金上，结果生成的结构无一例外全部崩溃。

为什么用于小分子的生成模型无法直接照搬到固态晶体？其内在根源在于两大截然不同的物理本性：

       有机小分子 (Molecular Matter)              固态无机晶体 (Periodic Solid Crystals)
  ─────────────────────────────────────     ─────────────────────────────────────────
  1. 局域离散共价键 (明确的键长与键级)        1. 离域非共价相互作用 (离子键、金属电子海、配位多面体)
  2. 孤立有限系统 (边界为真空，尺寸固定)      2. 无限周期性格子 (周期性边界条件 PBC，体积无限)
  3. 欧氏空间刚体平移旋转                   3. 230 个空间群对称性约束与晶胞应变张量变换

如果模型不能显式对晶胞基矢矩阵的应变形变、原子在三维环面上的周期性折叠以及无限镜像原子间的相互作用进行物理保真建模，生成稳定晶体便无从谈起。

本章将系统拆解无机晶体生成的数理表征三元组，并深度剖析近年来该领域的三大里程碑巨作——CDVAE、DeepMind GNoME 与 2025 年微软发表于 Nature 的 MatterGen！

2. 晶体生成的数学三元组： $x = (L, \mathbf{S}, \mathbf{Z})$

在数学形式化层面，任意一块理想单晶固态材料，均由其最小重复单元——**晶胞（Unit Cell）**唯一确定。一个晶胞由以下三元组构成：

x = (L, \mathbf{S}, \mathbf{Z})

2.1 晶胞基矢矩阵 $L \in \mathbb{R}^{3 \times 3}$

定义晶胞在三维实空间中的三个基矢为 $\mathbf{a}, \mathbf{b}, \mathbf{c} \in \mathbb{R}^3$ ：

L = [\mathbf{a}, \mathbf{b}, \mathbf{c}] = \begin{bmatrix} a_x & b_x & c_x \\ a_y & b_y & c_y \\ a_z & b_z & c_z \end{bmatrix}

晶胞体积：由标量三重积严格决定： $V_{\text{cell}} = \det(L) = \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) > 0$ ；
晶格参数（Lattice Parameters）：包含 3 个边长与 3 个夹角：
$a = \|\mathbf{a}\|, \quad b = \|\mathbf{b}\|, \quad c = \|\mathbf{c}\|$
$\alpha = \angle(\mathbf{b}, \mathbf{c}), \quad \beta = \angle(\mathbf{a}, \mathbf{c}), \quad \gamma = \angle(\mathbf{a}, \mathbf{b})$
旋转约化：任意对晶胞施加刚体旋转 $R \in SO(3)$ ，体系的本质物理性质保持不变。通常通过下三角或上三角 QR 分解消除该旋转自由度。

2.2 分数坐标 $\mathbf{S} \in [0, 1)^{N \times 3}$ 与三维环面拓扑

每个原子在晶胞中的实际笛卡尔坐标 $\mathbf{r}_i \in \mathbb{R}^3$ 是连续实数，但我们在生成模型中必须将其解耦为分数坐标（Fractional Coordinates）：

\mathbf{r}_i = L \mathbf{s}_i = s_{i,1}\mathbf{a} + s_{i,2}\mathbf{b} + s_{i,3}\mathbf{c}, \quad \text{其中 } \mathbf{s}_i \in [0, 1)^3

环面拓扑（3-Torus $\mathbb{T}^3$ ）：在周期性边界条件（PBC）下，分数坐标 $0.0$ 与 $1.0$ 在物理上是完全重合的同一个点！因此，原子位移的扩散加噪不能采用普通的欧氏高斯噪声，而必须采用定义在环面上的卷绕正态分布（Wrapped Normal Distribution）。
最小镜像约定（Minimum Image Convention）：两原子 $i$ 与 $j$ 在无限晶格中的真实最近邻位移向量，必须遍历跨越周期性晶胞边界的全部格矢平移：
$\mathbf{r}_{ij} = L \left( (\mathbf{s}_i - \mathbf{s}_j) - \text{round}(\mathbf{s}_i - \mathbf{s}_j) \right)$

2.3 原子化学序数 $\mathbf{Z} \in \{1, \dots, 118\}^N$

代表晶胞内 $N$ 个原子的化学元素种类（原子序数），以离散类别或 One-Hot 编码形式嵌入。

3. 里程碑 1：晶体扩散变分自编码器（CDVAE, ICLR 2021）

2021 年，MIT 的 Tian Xie 与 Jeffrey Grossman 团队提出了CDVAE（Crystal Diffusion Variational Autoencoder），正式拉开了扩散模型在晶体物理学中应用的帷幕。

图 5-1：三维周期性晶体晶胞基矢与原子分数坐标协同去噪扩散生成机制示意图

CDVAE 的核心解耦设计：

编码器（Encoder）：使用具备周期性图卷积（Periodic GNN）的网络，将周期性晶体结构压缩为一个全局潜在向量 $z \in \mathbb{R}^d$ ；
两阶段生成解耦：
- 第一步：解码器首先从潜变量 $z$ 预测晶胞基矢矩阵 $L$ 与晶胞内包含的原子数量 $N$ ；
- 第二步：在固定的晶胞 $L$ 内部，从纯高斯白噪声初始化的原子分数坐标 $\mathbf{s}^{(T)}$ 出发，利用等变力场预测网络逐步执行退火朗之万动力学去噪（Annealed Langevin Dynamics）：
  $\mathbf{s}^{(t-1)} = \mathbf{s}^{(t)} + \frac{\alpha_t}{2} \mathbf{f}_\theta(\mathbf{s}^{(t)}, L, \mathbf{Z}, t) + \sqrt{\alpha_t} \, \boldsymbol{\epsilon}_t \pmod 1$
这在物理上极其优美：它模拟了高温等离子体在特定晶胞体积限制下，原子相互排斥、自发寻找势能极小值点并最终“凝固”结晶的微观热力学全过程！

4. 里程碑 2：MatterGen（Nature 2025，微软 AI4Science）

尽管 CDVAE 开创了先河，但它将晶胞矩阵 $L$ 作为确定性回归输出，未能实现晶胞基矢与原子坐标的全局概率协同演化。2025 年 1 月，微软研究院 AI4Science 团队在国际顶刊 Nature 正式发表了重磅突破——MatterGen（Generative Model for Inorganic Materials Design）。

MatterGen 的革命性跨越：

晶格矩阵、原子坐标与元素类型的完全协同联合扩散（Joint Diffusion）： MatterGen 彻底摒弃了先生成晶胞再填充原子的两阶段割裂，而是建立了一个统一的多通道随机微分方程（SDE），在连续扩散空间中同时对晶格应变张量 $\boldsymbol{\varepsilon}_t \in \mathbb{R}^{3 \times 3}$ 、原子分数坐标 $\mathbf{s}_t \in \mathbb{T}^3$ 和连续松弛的元素嵌入向量 $\mathbf{z}_t$ 执行逆向时间去噪：
$\mathrm{d} \begin{bmatrix} \boldsymbol{\varepsilon}_t \\ \mathbf{s}_t \\ \mathbf{z}_t \end{bmatrix} = \mathbf{F}\left( \boldsymbol{\varepsilon}_t, \mathbf{s}_t, \mathbf{z}_t, t \right) \mathrm{d}t + \mathbf{G}(t) \mathrm{d}\mathbf{w}$
模块化属性适配器（Property Adapter Modules）：通过在冻结的预训练基础扩散骨干上插入微型适配器网络，MatterGen 首次实现了极其精准的多物理属性定制生成：
- 能够直接指定目标空间群（如立方晶系 $Fm\bar{3}m$ 、六方晶系 $P6_3/mmc$ ）；
- 能够直接输入目标弹性体积模量（ $K_V$ ）、磁性自旋序（铁磁/反铁磁）与电子带隙；
物理基态命中率飞跃：生成的结构中，处于热力学能量凸包之上或极度接近凸包（ $\Delta E_{\text{hull}} < 50 \text{ meV/atom}$ ）的高稳定材料比例提升了数倍，彻底解决了此前晶体生成模型常常输出“原子剧烈碰撞或晶格爆裂”的难题。

5. 里程碑 3：GNoME（Nature 2023，DeepMind）——主动学习数据飞轮

如果说 MatterGen 代表了生成算法的巅峰，那么 DeepMind 发表在 Nature 上的 GNoME（Graph Networks for Materials Exploration） 则展现了规模化（Scaling Laws）与主动学习数据飞轮的恐怖力量。

GNoME 团队依靠这一套周而复始自我进化的飞轮机制，将当时全球人类已知的稳定晶体数量直接扩充了 一个数量级（新增 220 万种全新稳定固态晶体结构），彻底改写了固态物理学的相图版图！

6. 工业级代码实战：带周期性边界条件（PBC）的晶格特征计算与扩散去噪

图 5-2：周期性边界条件（PBC）下晶体扩散生成前后的径向分布函数 g(r) 与原子硬球碰撞消除验证

python

"""
文件名: crystal_periodic_diffusion.py
功能: 实现晶体无限周期性边界条件 (PBC) 的张量级最小镜像几何计算，
      并构建周期性晶胞去噪预测网络的最小可运行原型。
"""

import math
import torch
import torch.nn as nn


class PeriodicCrystalGeometry:
    """
    无机固态晶体几何核心运算器:
    处理晶胞基矢矩阵 L、体积计算及周期性边界条件 (PBC) 下的距离矩阵
    """
    @staticmethod
    def get_cell_volume(L: torch.Tensor) -> torch.Tensor:
        """
        计算晶胞体积 V = det(L)
        L 维度: (batch_size, 3, 3), 基矢按列排布 [a, b, c]
        """
        return torch.det(L)

    @staticmethod
    def pbc_displacement(frac_coords_i: torch.Tensor,
                         frac_coords_j: torch.Tensor,
                         lattice: torch.Tensor) -> torch.Tensor:
        """
        基于最小镜像法则 (Minimum Image Convention) 计算两点在三维周期性环面上的真实实空间位移向量
        d_ij = L * ((s_i - s_j) - round(s_i - s_j))
        """
        # 分数坐标差值: (N, N, 3)
        delta_frac = frac_coords_i.unsqueeze(1) - frac_coords_j.unsqueeze(0)
        # 环面折叠到 [-0.5, 0.5) 区间
        delta_frac_pbc = delta_frac - torch.round(delta_frac)

        # 投影至三维真实欧氏空间: r = L * s
        # lattice: (3, 3)
        real_displacement = torch.matmul(delta_frac_pbc, lattice.T)
        return real_displacement


class CrystalDenoisingNetwork(nn.Module):
    """
    晶体协同扩散去噪网络核心骨架:
    输入带噪原子分数坐标 s_t 与时间步 t，预测逆向去噪得分力场向量
    """
    def __init__(self, num_species: int = 16, hidden_dim: int = 64):
        super().__init__()
        self.atom_embedding = nn.Embedding(num_species, hidden_dim)

        # 时间步正弦位置编码
        self.time_mlp = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim),
            nn.SiLU(),
            nn.Linear(hidden_dim, hidden_dim)
        )

        # 周期性等变距离消息传递层
        self.message_layer = nn.Sequential(
            nn.Linear(hidden_dim * 2 + 1, hidden_dim),
            nn.SiLU(),
            nn.Linear(hidden_dim, hidden_dim)
        )

        # 分数坐标去噪输出头 (预测位移受力向量 delta_s)
        self.coord_denoiser = nn.Sequential(
            nn.Linear(hidden_dim, hidden_dim),
            nn.SiLU(),
            nn.Linear(hidden_dim, 3)
        )

    def get_sinusoidal_time_embedding(self, timesteps: torch.Tensor, dim: int) -> torch.Tensor:
        """标量扩散时间步的傅里叶/正弦嵌入"""
        half_dim = dim // 2
        freqs = torch.exp(-math.log(10000) * torch.arange(0, half_dim, dtype=torch.float32, device=timesteps.device) / half_dim)
        args = timesteps[:, None].float() * freqs[None, :]
        return torch.cat([torch.sin(args), torch.cos(args)], dim=-1)

    def forward(self, frac_coords: torch.Tensor, atom_types: torch.Tensor,
                lattice: torch.Tensor, t: torch.Tensor) -> torch.Tensor:
        """
        前向去噪推断:
        frac_coords: (N, 3) - 当前带噪分数坐标
        atom_types: (N,)    - 原子类型索引
        lattice: (3, 3)     - 晶胞矩阵
        t: 标量扩散时间步
        """
        N = frac_coords.shape[0]
        h = self.atom_embedding(atom_types)  # (N, hidden_dim)

        # 融入时间步嵌入
        t_embed = self.get_sinusoidal_time_embedding(t.view(1), h.shape[-1])
        t_feat = self.time_mlp(t_embed).expand(N, -1)
        h = h + t_feat

        # 计算周期性边界条件下的所有原子间实空间距离
        disp = PeriodicCrystalGeometry.pbc_displacement(frac_coords, frac_coords, lattice)
        dist = torch.norm(disp, dim=-1, keepdim=True)  # (N, N, 1)

        # 消息聚合: 聚合邻居原子的周期性几何信息
        h_i = h.unsqueeze(1).expand(-1, N, -1)
        h_j = h.unsqueeze(0).expand(N, -1, -1)
        edge_features = torch.cat([h_i, h_j, dist], dim=-1)  # (N, N, 2*hidden + 1)

        messages = self.message_layer(edge_features)
        # 消除自身对自身的消息贡献
        eye_mask = torch.eye(N, dtype=torch.bool, device=h.device).unsqueeze(-1)
        messages = torch.where(eye_mask, torch.zeros_like(messages), messages)

        # 邻域求和池化
        h_updated = h + messages.sum(dim=1)

        # 预测去噪分数坐标修正量
        delta_frac_coords = self.coord_denoiser(h_updated)
        return delta_frac_coords


# =====================================================================
# 单元验证模块: 验证 PBC 最小镜像距离与晶体扩散去噪前向传播
# =====================================================================
if __name__ == "__main__":
    torch.manual_seed(42)

    print(">>> [STEP 1] 构建一个立方钙钛矿晶胞 (如钛酸锶 SrTiO3 模型)...")
    # 立方晶胞: a = b = c = 3.905 Å, 夹角 90°
    lattice = torch.eye(3) * 3.905
    vol = PeriodicCrystalGeometry.get_cell_volume(lattice.unsqueeze(0)).item()
    print(f"    晶胞矩阵已加载, 真实晶胞体积: {vol:.4f} Å³")

    # 钙钛矿 5 原子分数坐标: Sr(0,0,0), Ti(0.5,0.5,0.5), O(0.5,0.5,0), O(0.5,0,0.5), O(0,0.5,0.5)
    frac_coords = torch.tensor([
        [0.0, 0.0, 0.0],
        [0.5, 0.5, 0.5],
        [0.5, 0.5, 0.0],
        [0.5, 0.0, 0.5],
        [0.0, 0.5, 0.5]
    ])
    atom_types = torch.tensor([0, 1, 2, 2, 2])  # 假定类别编码

    print(">>> [STEP 2] 检验周期性边界条件 (PBC) 几何跨边界连通性...")
    # 检验 corner(0,0,0) 与相邻晶胞内原子的最小距离
    disp_matrix = PeriodicCrystalGeometry.pbc_displacement(frac_coords, frac_coords, lattice)
    dist_matrix = torch.norm(disp_matrix, dim=-1)
    print(f"    Sr 原子到 Ti 原子的三维对角最近邻真实距离: {dist_matrix[0, 1].item():.4f} Å")
    expected_dist = 3.905 * math.sqrt(3) / 2.0
    assert abs(dist_matrix[0, 1].item() - expected_dist) < 1e-4, "PBC 最小镜像计算距离错误！"

    print(">>> [STEP 3] 执行晶体扩散去噪网络前向单步计算...")
    denoiser = CrystalDenoisingNetwork(num_species=8, hidden_dim=32)
    t = torch.tensor([150])  # 中间扩散时间步

    # 注入微小扰动噪声
    noisy_coords = frac_coords + torch.randn_like(frac_coords) * 0.05
    predicted_corrections = denoiser(noisy_coords, atom_types, lattice, t)

    print(f"    输出去噪位移场尺寸: {predicted_corrections.shape}")
    assert predicted_corrections.shape == (5, 3), "去噪网络输出张量维度异常！"
    print(">>> [SUCCESS] 晶体无限周期性几何与扩散去噪模块单元验证全部通过！")

7. 本章小结

本章深入跨越了有机小分子与固态晶体材料之间的物理鸿沟：

确立了晶体数理三元组：形式化阐明了晶胞基矢矩阵 $L$ 、三维环面分数坐标 $\mathbf{S} \in \mathbb{T}^3$ 与元素种类 $\mathbf{Z}$ 在周期性边界条件（PBC）下的严密处理法则；
复盘了晶体生成三大里程碑：
- CDVAE 实现了潜空间晶胞与朗之万动力学原子去噪的两阶段解耦；
- 2025 年 Nature 的 MatterGen 实现了应变张量、坐标与元素的全联合协同扩散与多物理属性精准适配；
- DeepMind GNoME 演示了主动学习等变图力场飞轮扩充 220 万新晶体的工业级规模化实践。
编写了可运行的 PBC 几何张量引擎：亲手构建了最小镜像约定距离计算与晶体扩散去噪网络。

然而，计算机生成出的晶体无论在几何上多么规整美观，它在真实物理世界中是否真的能够稳定存在？这需要直面第一道生死玄关——进入 第 06 章：现实引力（热力学凸包稳定性与动力学声子谱检验）！

第 05 章：无机晶体生成实战：从 CDVAE 到 MatterGen 与 GNoME ​

1. 引言：有机小分子与固态晶体的物理鸿沟 ​

2. 晶体生成的数学三元组：x=(L,S,Z)x = (L, \mathbf{S}, \mathbf{Z})x=(L,S,Z) ​

2.1 晶胞基矢矩阵 L∈R3×3L \in \mathbb{R}^{3 \times 3}L∈R3×3 ​

2.2 分数坐标 S∈[0,1)N×3\mathbf{S} \in [0, 1)^{N \times 3}S∈[0,1)N×3 与三维环面拓扑 ​

2.3 原子化学序数 Z∈{1,…,118}N\mathbf{Z} \in \{1, \dots, 118\}^NZ∈{1,…,118}N ​

3. 里程碑 1：晶体扩散变分自编码器（CDVAE, ICLR 2021） ​

CDVAE 的核心解耦设计： ​

4. 里程碑 2：MatterGen（Nature 2025，微软 AI4Science） ​

MatterGen 的革命性跨越： ​

5. 里程碑 3：GNoME（Nature 2023，DeepMind）——主动学习数据飞轮 ​

6. 工业级代码实战：带周期性边界条件（PBC）的晶格特征计算与扩散去噪 ​

7. 本章小结 ​