CAD+MLIP：高效计算固体振动自由能与热力学性质的技术实践

1. 项目概述：当机器学习势函数遇见固体热力学

在计算材料科学领域，预测一个材料在特定温度和压力下的稳定性，本质上是在计算它的自由能。自由能是决定材料相变、合成可行性乃至最终性能的“总指挥”。然而，精确计算自由能，尤其是其中的振动熵贡献，一直是个老大难问题。传统上，热力学积分（TI）被视为“金标准”，但其计算成本高得吓人，动辄需要数百万个分子动力学（MD）步，且实现路径复杂，对新手极不友好。准谐近似（QHA）虽然便宜，但它本质上是个“零温”近似，对于非谐效应显著的材料（比如某些高温相或软模材料），预测结果常常失灵。

近年来，机器学习原子间势（MLIP）的崛起，让我们能以接近第一性原理计算的精度，实现比其快几个数量级的分子动力学模拟。这为解决自由能计算的“采样”难题提供了新武器。但光有快枪还不够，我们还需要一个能高效利用MD采样数据、准确提取振动信息的方法。协方差原子位移（CAD）方法，正是这样一把“好刀”。它绕开了复杂的积分路径，直接从MD轨迹中原子位移的统计涨落（协方差矩阵）出发，反推出一个“有效”的力常数矩阵，进而得到包含非谐效应的振动频率和所有衍生的热力学性质。

简单来说，这个项目的核心就是：用机器学习势函数驱动高效的分子动力学模拟，再用CAD方法从模拟数据中“榨取”出振动熵和自由能。它瞄准的是那些需要大量统计采样或涉及复杂非谐效应的材料体系，比如合金、高温相变材料或者像锂这样的轻元素金属。下面，我就结合自己在锂金属体系上的实操经验，拆解这套方法的工作流、关键细节、避坑指南以及它和QHA、TDEP等方法的对比。

2. 技术路线解析：为什么是CAD+MLIP？

2.1 方法选型的逻辑：在精度与效率的钢丝上行走

计算材料的热力学性质，本质是在精度、效率和通用性之间做权衡。下图清晰地展示了各种原子间相互作用模型的帕累托前沿：

精度-效率权衡与CAD方法定位

CAD方法的核心输入是高质量的MD轨迹。因此，势函数的选择直接决定了整个计算的可行性与可靠性。DFT虽然精度高，但其O(N³)的计算复杂度使得产生一条收敛的、足够长的MD轨迹成本巨大，难以进行系统的温度-压力扫描。经验势函数速度最快，但其拟合精度往往有限，且泛化能力差，难以保证在训练域外（如高压、高温）的可靠性。

MLIP恰恰卡在了这个甜蜜点上。它通过神经网络学习DFT数据，能以接近DFT的精度，实现O(N)甚至更优复杂度的力场评估。这意味着，我们可以用单台GPU工作站，在几十分钟内完成一个数千原子体系、数万步的NVT模拟，并获得可靠的原子受力。这种“降维打击”般的能力，使得在保持精度的前提下进行大规模的CAD计算成为可能。我选择的NequIP势函数，因其基于等变神经网络架构，能隐式捕捉长程相互作用，对力（决定动力学）的预测尤其准确，这对于振动性质的计算至关重要。

注意：选择MLIP时，切勿只看重其报告的能量/力均方根误差（RMSE）。务必针对你的目标性质（如晶格常数、弹性常数、声子谱）进行系统性的基准测试。一个在能量上误差很小的势，可能在声子软化或相变能垒上表现不佳。

2.2 CAD方法原理：从原子舞步到热力学量

CAD方法的巧妙之处在于其深刻的统计物理内涵。它不直接求解复杂的运动方程，而是认为，在平衡态MD模拟中，原子围绕其平衡位置的“舞蹈”（位移涨落）已经编码了全部的有效相互作用信息。

1. 核心公式推导设想一个包含N个原子的超胞。在MD模拟中，我们记录每个原子i在µ方向（x, y, z）偏离其理想晶格位置的位移 ( x_{i\mu} )。计算整个轨迹中所有原子位移的协方差矩阵 ( \mathbf{\Sigma}U )，其矩阵元为： [ \Sigma{i\mu, j\nu} = \langle x_{i\mu} x_{j\nu} \rangle ] 这里 (\langle \cdot \rangle) 表示系综平均。这个3N×3N的矩阵，捕捉了所有原子运动之间的关联。

接下来是关键一步：将这个位移协方差矩阵与系统的有效力常数矩阵联系起来。在谐振子近似下，可以通过统计力学严格推导出，质量加权的协方差矩阵 (\tilde{\mathbf{\Sigma}})（即 (\tilde{\Sigma}{i\mu, j\nu} = \Sigma{i\mu, j\nu} / \sqrt{m_i m_j})）与质量约化的力常数矩阵 (\tilde{\mathbf{C}})（即 (\tilde{C}{i\mu, j\nu} = C{i\mu, j\nu} / \sqrt{m_i m_j})）满足以下关系： [ \tilde{\mathbf{C}} = \frac{1}{\beta} \tilde{\mathbf{\Sigma}}^{-1}, \quad \beta = \frac{1}{k_B T} ] 这个公式是CAD方法的基石。它告诉我们，通过计算MD轨迹中原子位移的协方差并求逆，再乘以温度因子，我们就可以得到系统在当前温度下的“有效”力常数矩阵。这个“有效”力常数已经包含了所有非谐效应（原子间的非简谐相互作用）对平均结构的平均影响。

2. 从矩阵到物理量得到有效力常数矩阵 (\tilde{\mathbf{C}}) 后，对其进行对角化，即可得到3N个本征值 (\lambda_k) 和对应的本征向量。每个本征值对应一个正则模的频率（扣除平移和旋转的零频模后）： [ \omega_k = \sqrt{\lambda_k} ] 这些 (\omega_k) 就是包含了非谐效应的“有效”振动频率。一旦有了频率，所有热力学量都可以通过标准的谐振子统计公式计算：

• 振动熵： [ s_{\text{vib}} = k_B \sum_{k} \left[ \frac{\beta \hbar \omega_k}{e^{\beta \hbar \omega_k} - 1} - \ln(1 - e^{-\beta \hbar \omega_k}) \right] ]
• 振动自由能（亥姆霍兹）： [ f_{\text{vib}} = \frac{1}{N} \sum_{k} \left[ \frac{1}{2} \hbar \omega_k + \frac{1}{\beta} \ln(1 - e^{-\beta \hbar \omega_k}) \right] ]
• 声子色散与态密度：利用本征向量，可以将这些频率映射到倒空间的第一布里渊区，从而绘制出有限温度下的声子色散关系图和态密度。

实操心得：公式 ( \tilde{\mathbf{C}} = \beta^{-1} \tilde{\mathbf{\Sigma}}^{-1} ) 在推导时假设了势能面是谐性的。对于真实体系，CAD方法实际上构建了一个“最佳”的谐性势，使得在该谐性势下产生的位移涨落与真实MD模拟的涨落一致。因此，它提取的是包含非谐效应的有效谐性频率，这是其优于传统零温谐性计算的核心。

3. 完整工作流与实操要点

基于CAD方法计算材料振动性质，可以梳理为下图所示的清晰工作流。这个流程高度模块化，每个环节都有需要特别注意的细节。

flowchart TD A[确定目标材料与相] --> B[获取/训练MLIP] B --> C{MLIP精度<br>是否达标?} C -- 否 --> D[在目标P/T范围内<br>重新训练或微调MLIP] D --> B C -- 是 --> E[执行NVT/NPT MD模拟<br>（存储轨迹快照）] E --> F[CAD后处理] subgraph F [CAD后处理核心步骤] F1[计算并对称化<br>位移协方差矩阵] F2[对角化得到<br>有效频率与声子谱] F3[代入谐振子公式<br>计算热力学性质] end F --> G[输出: 熵/自由能/声子谱等]

3.1 第一步：势函数的准备与验证

这是整个计算的基石，决不能马虎。以我使用的锂金属NequIP势函数为例：

1. 数据准备：训练数据应覆盖你感兴趣的温度和压力范围。对于锂，我包含了从50K到500K、0到5 GPa范围内的多种原子构型，通过主动学习（Active Learning）不断补充MD采样中遇到的新构型，确保势函数在相关相空间有良好覆盖。
2. 训练与测试：将数据集按8:1:1分为训练集、验证集和测试集。除了监控能量和力的RMSE，我特别关注以下性质的测试：
   • 静态性质：平衡晶格常数、弹性常数、空位形成能。
   • 动态性质：在几个关键温度点（如100K, 300K）计算声子谱，与DFT计算结果或实验数据对比。声子谱的吻合度是预测振动熵可靠性的直接证据。
   • 相稳定性：计算不同晶体结构（如BCC, FCC）在0K下的能量差，与DFT结果对比。
3. 部署与效率测试：在目标计算平台（如带GPU的服务器）上测试单步力和能量的评估速度。对于CAD，我们需要进行数万步的MD，因此每秒的步数（steps/second）至关重要。确保你的MLIP接口与MD引擎（如LAMMPS, ASE）兼容良好。

3.2 第二步：分子动力学模拟设置

MD模拟为CAD提供原材料——平衡轨迹。设置不当会导致结果不收敛或完全错误。

1. 系综选择：通常使用NVT（恒温恒容）系综进行生产模拟。但为了确定该温度下的平衡体积，需要先进行NPT（恒温恒压）模拟。我的流程是：先运行NPT模拟（~20 ps）使体系弛豫到平衡体积，取最后一段的平均体积作为该温度下的平衡体积，再用此体积进行NVT模拟用于CAD分析。
2. 超胞大小与模拟时长：这是收敛性测试的核心。我以BCC锂在300K为例：
   • 系统大小：测试了从64到2000个原子不等的超胞。如图2a所示，振动熵和自由能在约1000个原子时基本收敛。但在高通量筛选中，使用500个原子的超胞，其自由能误差仅~0.04 meV/atom，在可接受范围内，能大幅节省计算资源。
   • 模拟时长：测试了从5000到40000个MD步。如图2b所示，约20000个步长（使用2 fs时间步，即40 ps）后，熵和自由能收敛。这比热力学积分所需的数百万步少了两个数量级。
3. 参数设置：
   • 时间步长：对于锂这样的轻元素，我通常使用1-2 fs。过大的步长（如5 fs）可能导致能量漂移，尽管有研究指出对声子谱计算可能有益，但需谨慎测试。
   • 热浴：使用Nosé-Hoover链或Langevin热浴来精确控温。确保热浴参数设置合理，避免温度振荡。
   • 初始速度：使用不同的随机种子生成初始速度，重复运行2-3次模拟，以评估统计误差。如图2所示，对于稳定相，不同初始速度带来的结果差异（标准误差）非常小。
   • 采样间隔：MD轨迹中相邻快照是高度相关的。为了获得独立的样本，需要设置采样间隔（stride）。我通常每10-20步保存一次快照，并确保总采样帧数（用于计算协方差矩阵的帧数）足够（通常需要数千帧）。

3.3 第三步：CAD后处理核心算法与实现

拿到MD轨迹后，真正的CAD计算是后处理过程。我通常用Python编写脚本，流程如下：

1. 轨迹处理与位移计算：

   import numpy as np from ase.io import read from ase.geometry import get_layers # 读取轨迹，假设已经去除了平衡部分 traj = read('md_trajectory.xyz', index=':') # 读取所有帧 # 获取理想晶格位置 (通常取第一帧或平均位置作为参考) reference_positions = traj[0].get_positions() # 或者使用平均位置（更稳定） # avg_positions = np.mean([atoms.get_positions() for atoms in traj], axis=0) num_atoms = len(traj[0]) num_frames = len(traj) displacements = np.zeros((num_frames, num_atoms, 3)) for i, atoms in enumerate(traj): # 将原子映射回原胞，消除周期性边界条件引起的跳跃 # 这里假设体系是完整的晶体，没有扩散 # 对于复杂情况，需要更精细的Wigner-Seitz原胞映射 current_pos = atoms.get_positions() # 简单差值（适用于小位移，无扩散） disp = current_pos - reference_positions # 考虑周期性边界条件修正 (PBC wrap) cell = atoms.get_cell() disp = disp - np.round(disp / cell.lengths()) * cell.lengths() displacements[i] = disp

2. 构建与对称化协方差矩阵：

   # 重塑位移数组: (num_frames, 3N) disp_flat = displacements.reshape(num_frames, -1) # shape: (F, 3N) # 计算协方差矩阵 covariance_matrix = np.cov(disp_flat, rowvar=False, bias=True) # shape: (3N, 3N) # **关键步骤：对称化** # 未经对称化的协方差矩阵可能不满足晶体的空间群对称性，导致声子谱出现非物理的劈裂。 # 对称化操作将矩阵元投影到其对称性允许的空间。 def symmetrize_covariance_matrix(cov_matrix, atoms, symprec=1e-5): from ase.spacegroup import Spacegroup # 获取空间群操作 sg = Spacegroup(atoms.get_chemical_symbols()[0]) # 简单示例，实际需根据结构确定 # 获取对称操作（旋转矩阵和平移向量） sym_ops = sg.get_symop() # 这是一个简化的示意。实际实现需要： # 1. 将协方差矩阵按原子索引分块。 # 2. 对每一对原子块 (i,j)，应用所有对称操作，找到等价原子对。 # 3. 将所有等价矩阵块取平均。 # 此过程较复杂，通常可利用现有软件库（如phonopy的对称性例程）或自行实现。 # 对称化后，矩阵的本征值/向量将更平滑，声子谱更合理。 return symmetrized_matrix # 应用对称化 cov_sym = symmetrize_covariance_matrix(covariance_matrix, traj[0])

3. 计算有效频率与热力学量：

   # 质量加权 masses = traj[0].get_masses() # 形状 (N,) mass_vector = np.repeat(np.sqrt(masses), 3) # 形状 (3N,) # 创建质量加权矩阵 M_ij = 1/sqrt(m_i * m_j) mass_matrix = np.outer(1/mass_vector, 1/mass_vector) # 形状 (3N, 3N) cov_mass_weighted = cov_sym * mass_matrix # 求逆得到有效力常数矩阵 (公式2) beta = 1.0 / (kB * temperature) # temperature 为模拟温度 effective_force_constants = (1.0 / beta) * np.linalg.pinv(cov_mass_weighted) # 使用伪逆更稳定 # 对角化 eigenvalues, eigenvectors = np.linalg.eigh(effective_force_constants) # 剔除3个平移和3个旋转零频模（对于非线形分子） # 通常通过设置一个小的频率阈值来实现，如 1e-5 THz threshold = 1e-5 nonzero_indices = np.where(np.sqrt(np.abs(eigenvalues)) / (2*np.pi) > threshold)[0] omega_k = np.sqrt(np.abs(eigenvalues[nonzero_indices])) # 角频率，取绝对值防止负值 # 注意：负的本征值（对应虚频）可能出现在亚稳或不稳定相中，需谨慎分析。 # 计算振动熵 (公式3) beta_hbar_omega = beta * hbar * omega_k exp_term = np.exp(beta_hbar_omega) # 避免除零，对极小的omega_k做处理 with np.errstate(divide='ignore', invalid='ignore'): term1 = beta_hbar_omega / (exp_term - 1) term1 = np.nan_to_num(term1, nan=0.0) # 当omega_k->0, term1->1 term2 = np.log(1 - np.exp(-beta_hbar_omega)) term2 = np.nan_to_num(term2, nan=0.0) # 当omega_k->0, term2->0 s_vib_per_mode = kB * (term1 - term2) s_vib_total = np.sum(s_vib_per_mode) / num_atoms # 每原子的熵

3.4 第四步：自由能构建与相图预测

单个相的熵和自由能计算出来后，最终目标是构建自由能面，预测相变。

1. 总自由能构成：对于单质体系，每原子的总亥姆霍兹自由能 ( f(v, T) ) 为： [ f(v, T) = u_{\text{static}}(v) + f_{\text{vib}}(v, T) + f_{\text{el}}(v, T) ]

   • ( u_{\text{static}} )：静态晶格能。即0K下完美晶体的势能。这部分强烈建议使用更高精度的DFT单点能计算，而不是MLIP的能量。因为自由能差往往在meV/atom量级，DFT对静态能量的计算更可靠，且计算成本相对MD模拟可忽略。
   • ( f_{\text{vib}} )：振动自由能。即由CAD计算得到的部分（公式5）。
   • ( f_{\text{el}} )：电子自由能。在金属中，温度升高会激发电子，贡献额外的熵和自由能。可通过DFT计算有限温度下的电子态密度，并用费米-狄拉克分布积分得到。对于锂，在300K以下，此项贡献很小（约0.04 ( k_B )），但不可忽略。

2. 从亥姆霍兹自由能到吉布斯自由能：实验相变通常在恒定压力下发生，因此我们需要吉布斯自由能 ( g(P, T) )。 [ g(P, T) = \min_v [ f(v, T) + P v ] ] 实际操作中，我采用“等温线法”：

   • 在目标温度T下，选取一系列体积 ( {v_i} )。
   • 对每个 ( v_i )，运行NVT-MD模拟，并用CAD计算 ( f_{\text{vib}}(v_i, T) )。
   • 加上该体积下的 ( u_{\text{static}}(v_i) ) 和 ( f_{\text{el}}(v_i, T) )，得到 ( f(v_i, T) )。
   • 用状态方程（如Birch-Murnaghan）拟合 ( f(v, T) ) 曲线，找到给定压力P下使 ( f + Pv ) 最小的体积 ( v_{eq} )，对应的最小值即为 ( g(P, T) )。
   • 如图4a所示，通过此方法找到的平衡体积与NPT-MD模拟直接得到的平均体积高度一致，验证了流程的自洽性。

3. 相变温度预测：对竞争相（如锂的FCC和BCC）重复上述过程，得到它们随温度变化的吉布斯自由能曲线 ( g_{\text{FCC}}(T) ) 和 ( g_{\text{BCC}}(T) )。相变温度 ( T_c ) 由两条曲线的交点决定。如图4f所示，CAD、TDEP和QHA均预测锂的FCC→BCC马氏体相变发生在~150K，高于实验值80K。这揭示了当前方法在预测微小自由能差（~1 meV/atom）时面临的普遍挑战。

4. 关键问题、排查技巧与经验实录

在实际操作中，你会遇到各种问题。下面是我踩过坑后总结的排查清单和应对策略。

4.1 收敛性判断：你的模拟跑够了吗？

这是CAD计算中最常见的问题。不收敛的结果毫无意义。

• 症状：熵或自由能值随模拟步长或系统大小剧烈波动，没有达到平台。
• 诊断与解决：
  1. 系统大小：如图2a，逐步增大超胞（如 4x4x4, 6x6x6, 8x8x8...），绘制目标性质（如 ( s_{\text{vib}} )）随原子数的变化图。当变化小于你的误差容限（例如，对于自由能，1 meV/atom是可接受的化学精度阈值）时，认为收敛。
  2. 模拟时间：如图2b，固定系统大小，增加模拟总步长。分析最后一段轨迹（如后20 ps）计算的性质是否稳定。务必丢弃初始弛豫阶段的数据（通常前5-10 ps），只从平衡后的轨迹开始采样。
  3. 采样独立性：计算原子位移的自相关函数，确定衰减时间。采样间隔应大于此衰减时间，以确保帧间独立。如果资源有限，可以尝试对轨迹进行块平均（block averaging）来估计误差。
  4. 快速检查：一个经验法则是，对于简单金属/半导体，~1000原子、~40 ps的NVT模拟通常能给出收敛的熵值。对于更软或非谐性更强的材料，可能需要更长的模拟时间。

4.2 声子谱出现虚频：是物理现象还是计算假象？

在CAD计算中，对角化有效力常数矩阵有时会得到负的本征值，对应虚频（imaginary frequency）。

• 可能原因1：统计采样不足。这是最常见的原因，尤其在小体系或短时模拟中。原子位移的协方差矩阵未能准确捕捉所有振动模式的相关性。
  • 解决：增加模拟时间或系统大小。对称化协方差矩阵对此有巨大改善。如图2所示，对称化对熵值影响不大，但对获得平滑、物理的声子色散曲线至关重要。
• 可能原因2：模拟失稳或相变。在接近相变温度或压力下，当前相可能亚稳态，MD模拟中可能出现局部结构起伏或预相变，导致动力学不稳定。
  • 解决：检查MD轨迹的均方根位移（RMSD）和能量漂移。如果RMSD异常大或能量持续上升，说明模拟失稳。尝试从不同初始速度重复模拟，如果结果发散（如图4f中接近熔点时），则表明该条件下相不稳定。
• 可能原因3：势函数精度不足。MLIP在训练域外预测不力，给出了不准确的力，导致非物理的动力学行为。
  • 解决：用DFT计算相同超胞的声子谱进行对比。如果DFT没有虚频而CAD有，很可能是势函数问题。需要在当前条件下补充训练数据，重新训练或微调势函数。

核心经验：不要一看到虚频就认为是方法失败。对于某些亚稳相或相变路径上的鞍点，虚频是物理真实的，指示着失稳模式。关键在于结合模拟稳定性分析和势函数验证来综合判断。

4.3 低温下的失败：经典模拟的量子困境

如图4e所示，在低温区（<100K），用经典MD结合CAD公式（5）计算的吉布斯自由能 ( g = f_{\text{vib}} + Pv ) 与实验和QHA结果出现显著偏差，甚至出现非单调行为，这违反了热力学第二定律。

• 根源：对于锂这样的轻元素，在低温下，核量子效应（NQE）变得显著。经典MD基于玻尔兹曼分布和牛顿力学，完全忽略了量子隧穿、零点能等效应。在10K时，锂的热德布罗意波长可达~2 Å，与原子位移相当，经典描述完全失效。
• 解决方案：
  1. 路径积分分子动力学（PIMD）：这是最严格的处理NQE的方法，但计算成本比经典MD高1-2个数量级。可以将PIMD与CAD结合，用路径积分采样的轨迹来构建协方差矩阵。
  2. 量子校正：对经典CAD结果进行事后校正。例如，使用Wigner量子校正或更简单的谐性零点能校正。对于自由能，一个常见的近似是：( F_{\text{quantum}} \approx F_{\text{classical}} + F_{\text{ZPE}} )，其中零点能校正 ( F_{\text{ZPE}} = \frac{1}{2}\sum_k \hbar \omega_k )。
  3. 使用量子热力学公式：如图4e所示，如果使用公式 ( f = u - Ts )（其中u是MD模拟的平均内能，s是CAD熵），在高温区与实验符合得更好。这是因为该公式在一定程度上隐含了部分非谐效应，但在低温区仍会失败。
  • 我的建议：对于含有氢、锂等轻元素或低温下的研究，必须考虑核量子效应。如果计算资源允许，PIMD+CAD是最佳选择。否则，需要明确说明低温结果的局限性，并采用量子校正。

4.4 与其他方法的对比：CAD vs. QHA vs. TDEP

为了让你更清楚何时选择CAD，这里做一个直接对比：

固体振动性质计算方法对比

个人体会：CAD的最大优势在于其**“无脑”式的后处理**。一旦有了收敛的MD轨迹，计算就是几个矩阵操作，几乎没有调参的烦恼。而TDEP需要进行力匹配优化，其结果可能依赖于初始猜力和优化算法。对于锂这样的金属，两者结果非常接近（图4c,d,f），说明在满足收敛条件的前提下，CAD是TDEP的一个极佳替代。但对于某些非常复杂的体系，TDEP通过拟合可能对噪声更稳健，这点需要根据具体体系测试。

5. 总结与展望：CAD方法的定位与最佳实践

经过在锂金属体系上的完整实践，我认为CAD方法结合MLIP，为计算固体有限温度振动性质提供了一个在效率、精度和易用性上取得优异平衡的方案。它绝非万能，但在其适用范围内非常强大。

最佳实践路线图：

1. 明确问题：你的材料是否强非谐？是否涉及高温相变？如果是，CAD或TDEP比QHA更合适。
2. 准备势函数：投入足够时间验证MLIP。声子谱测试是必选项。确保训练数据覆盖你的目标相空间。
3. 系统收敛性测试：这是绝对不能跳过的一步。针对你的具体体系和温度，进行系统大小和模拟时长的扫描，确定可靠的参数。
4. 执行与监控：运行MD时，实时监控温度、压力、能量和RMSD。一旦发现异常（如熔融、结构坍塌），立即分析原因。
5. 后处理与对称化：计算CAD时，务必对协方差矩阵进行对称化，这是获得物理声子谱的关键。
6. 交叉验证：如果可能，用QHA或TDEP在个别点上进行计算对比。与已有的实验数据（如比热、声子谱）对比。
7. 警惕量子效应：对于轻元素或低温（T < Debye温度），必须考虑核量子效应的校正，或直接使用PIMD。

最后一点心得：CAD计算出的熵本身可以作为一个亚稳性的描述符。如果一个相的CAD熵值异常高，或者在不同重复模拟中发散，即使MD轨迹没有显示明显的结构相变，也可能暗示着该相在模拟条件下是动力学不稳定的，或者存在一个能垒很低的其他相变路径。这为我们分析材料的相对稳定性提供了一个除能量和结构外的动态视角。

这套方法的价值在于其通用性。一旦流程打通，它可以像流水线一样应用于高通量的材料筛选，快速评估成千上万种材料在服役温度下的振动自由能，从而加速新型热电材料、超导材料或结构合金的发现。计算成本的降低，使得我们能够去探索那些传统方法难以触及的复杂相变和动力学过程。