对于结构多样的含硝基化合物，该方法能够提供稳健且易于解释的熔点预测结果

《Energetic Materials Frontiers》：Robustly interpretable melting-point prediction for structurally diverse nitro-containing compounds

　　冯晓宇|陈芳|宋思伟|王瑞辉|王毅|张庆华
西北工业大学航天学院，西安，710072，中国

**摘要**
对于含硝基的柔性有机分子而言，准确预测其熔点对于优化其在制药、能源材料和化学加工中的性能至关重要。通过使用九种代表性的含硝基化合物来研究分子间相互作用以及振动熵效应对熔点的影响，我们开发了一个双描述符框架，该框架整合了来自晶体结构的焓描述符和来自从头算分子动力学（AIMD）的熵描述符。通过比较线性多元回归（MLR）模型和非线性KNN模型，并结合SHAP分析，发现KNN模型（整合了分子间相互作用（BEC_max1-3、BEC_total）、分子范德华体积（VW）和分子振动特征（RMSD_ave、E_b和ΔE））具有高准确性（训练集：R2 = 0.83，MAE = 12.48°C；测试集：R2 = 0.79，MAE = 15.29°C）和可解释性。机制分析表明，熔化过程主要受焓和熵的共同作用支配。这一框架为新型含硝基柔性有机化合物的合理设计建立了结构-性质联系。

**1. 引言**
含硝基有机化合物在医学、化学工程和国防领域得到广泛应用。其中，三硝基甲苯（TNT）因其独特的80°C附近的可逆固液相变而特别值得关注，这使其成为几乎不可替代的熔融载体。因此，熔点是一个关键参数，决定了含硝基有机化合物的功能性质。这一基本的热力学参数不仅在化合物的初步鉴定和纯度分析中起着关键作用，也是获取沸点、水溶性、液体粘度等性质的重要依据。然而，实验测定的熔点仍然耗时且效率低下，无法从庞大的化学库中筛选出具有所需性质的目标化合物。因此，开发预测这种分子性质的工具是非常必要的。
近年来，出现了多种预测方法，包括基团贡献方法（GCM）、分子动力学（MD）模拟和定量结构-性质关系（QSPR）模型。值得注意的是，混合QSPR-机器学习（ML）框架受到了关注，这些框架利用了先进的分子描述符（如SMILES、Mordred、RDKit、MOLMAPs等）和ML算法（如图卷积网络（GCN）、深度神经网络（DNN）、核岭回归（KRR）、随机森林（RF）、SISSO和XGBoost）。尽管最先进的模型能够实现平均绝对误差（MAE）< 30°C（最佳MAE为24.7°C），但描述符的选择对性能的影响比算法复杂性更为重要。当前的描述符主要基于量子化学计算或固态理论编码分子结构特征。然而，这些描述符往往缺乏直观的物理化学可解释性，无法捕捉动态的分子内/间相互作用（图1a）。例如，尽管晶体堆积参数（如二聚体结合能和晶格振动模式）与熔化转变有直接的热力学相关性，但它们仍然被忽视。整合这些参数可以提高预测准确性，并提供关于相变现象的机制洞察。

**2. 方法**
**2.1. 模型分子的设计与合成**
为了构建高质量的熔点数据集，分子选择经过精心设计，以平衡结构相似性和可控的变化。4-氨基-3,5-二硝基-1H-吡唑（化合物1）作为核心骨架，确保了焓和静电基线的一致性。通过一步亲核取代在N1位置进行系统修改（补充信息，第1节；图2a，顶部），生成了具有不同长度的烷基链（C2–C4）和极性官能团（如硝酸酯（?ONO2）和叠氮化物（?N3）的衍生物。这种方法将扰动局限在侧链身份上，同时保留了共有的硝基-芳香骨架，从而分离了取代基性质对热行为的影响。

**2.2. 焓描述符的提取**
从热力学的角度来看，熔点的本质是固液相变过程中焓和熵的协同效应的结果（Tm=ΔHm/ΔSm）。这种协同作用不仅体现在分子间相互作用对熔化焓的显著调节上，还深刻反映在分子排列顺序变化对熵变化的决定性影响上。因此，熔点本质上是分子间相互作用能量（主导ΔH）和分子构型自由度（主导ΔS）的综合热力学表现。这为理解分子结构与热力学性质之间的相关性提供了理论基础。为了构建参数化的焓描述符，使用实验解析的晶体结构（表S1–S8）通过簇结合能（BEC）方法分析了分子间相互作用（图2c，右侧；图S2和表S9–S21）。

**2.3. 熵描述符的处理**
准确的熵确定长期以来一直是热力学建模中的一个重大挑战。为此，我们将量子化学计算与分子动力学模拟相结合。从量子化学优化的分子几何结构中得到的范德华体积（VW）用作熵描述符，以描述静态分子体积效应（图2d，左侧）。在柔性分子系统中，体积膨胀主要来源于空间位阻较大的基团（例如支链或取代环）。这些结构特征降低了晶体堆积效率，但在液态下允许完全的构象自由度，从而增加了熔化熵（ΔS）。此外，从头算分子动力学（AIMD）模拟用于推导五个与熵相关的描述符，以表征分子振动和构型变化（图2d）。在B97-3c理论水平上，对结构优化的分子进行了AIMD模拟；模拟在298.15 K的温度下进行，持续时间为2.0 × 10^5 fs。通过对结构变形和能量波动的系统性分析，量化了原子尺度上的熵相关现象。平均RMSD（RMSD_ave）是在最后1.0 × 10^4 fs的轨迹段（1.9 × 10^5–2.0 × 10^5 fs）上使用Visual Molecular Dynamics（VMD）软件计算的。为了确保对结构和动态性质研究的准确性，还获得了另一个描述符RMSD_ave（ALIGN），以消除由于相对于初始构象的平移和旋转运动可能引入的噪声和不确定性。

**3. 结果与讨论**
**3.1. 描述符相关性分析**
为了确定熵和焓描述符分类的合理性，我们系统地分析了不同描述符之间的相关性。相关矩阵通过热图可视化，颜色梯度从红色（强正相关，1）到蓝色（强负相关，?1），白色表示可忽略，揭示了不同描述符之间的相关性（图3a）。以VW和BEC_var之间的线为边界，热图的左下角和右上角部分颜色更深，不同类型描述符之间的相关性有显著差异，这意味着左下角的描述符在性质上更为相似，而右上角的描述符则更为相似，即左下角的描述符可以被认为是与焓相关的，而右上角则由所谓的与熵相关的描述符占据。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图3。(a) 不同描述符之间相关性的热图。(b) 分子簇（化合物1）中的结合能和描述焓的特征值。(c) 熔点与焓因子之间的关系。五个局部焓描述符（BEC_max1至BEC_max5）表现出强相关性（R = 0.81至0.98），尤其是在BEC_max1-2和BEC_max1-3之间（R = 0.98），表明在表征分子相互作用方面存在功能冗余。相反，BEC_var与其他焓参数表现出明显的负相关性（R = ?0.62至?0.80），反映了其在晶体稳定性相对于二聚体结合能参数的比较效应。与熵相关的描述符表现出补偿关系：VW与RMSD_ave（R = 0.80）和RMSD_ave（ALIGN）（R = 0.92）呈强正相关，反映了分子体积的增加如何增强构象灵活性——这通过增加的均方根偏差和放大的构象无序得到证明。相反，VW与E_b（R = ?0.81）之间的显著负相关性揭示了一个热力学权衡：体积膨胀可能会限制分子的振动自由度并减少振动能量。这项综合分析为选择预测含硝基柔性有机化合物熔点的多变量模型中的描述符建立了基本标准。3.2. 焓描述符的影响在分子簇中，很明显，簇内的相邻分子与中心分子（绿色区域）之间存在显著的弱相互作用。在每个相邻分子中，只有局部原子位点对整体相互作用有实质性贡献（红色区域）（图3b，顶部和图S2）。46,47 二聚体结合能计算量化了这些差异：化合物1显示出六个高能量的二聚体对，范围从?33.0到?46.0 kJ·mol?1，而其他构型的相互作用较弱（≤12 kJ·mol?1），甚至存在排斥作用（3.4 kJ·mol?1，图S2）。通过二聚体结合能的数值表示可以有效地突出这些相互作用的具体差异（图3b，中间）。初步分析显示熔点（Tm）与七个焓描述符之间存在中等但统计上显著的相关性，包括五种局部二聚体结合能组合：BEC_max1、BEC_max1-2、BEC_max1-3、BEC_max1-4、BEC_max1-5，以及两个全局结合能参数：BEC_total和BEC_var（图3b，底部）。在评估的参数中，平均总结合能（BEC_total）与Tm的相关性最强（R2 = 0.291，图3c），突显了其在晶体结构稳定性中的潜在作用。然而，BEC_total的有限预测能力表明，仅凭焓描述符不足以准确预测熔点。此外，在五种局部二聚体结合能组合中，BEC_max1-2（R2 = 0.225）和BEC_max1-3（R2 = 0.219）显示出最强的相关性。这表明熔点行为主要由最强的二聚体相互作用的协同作用决定，而不是累积的结合效应。这种位点特异性优势与晶体热稳定性中的“最弱环节”假设一致，即关键分子相互作用的破坏可以触发晶格崩塌。尽管这项研究的结果不够确凿，但它们仍然支持了提出的双描述符策略。这强调了对于分子结构存在显著差异的化合物，熔点不能忽视由烷基侧链等柔性结构引入的熵贡献。重要的是，焓描述符的参数化揭示了几个不同的特征：（1）分子内引入烷基链结构通常会增加分子体积，同时减弱了分子间相互作用；（2）侧链中含有极性官能团的化合物通常表现出增强的分子间相互作用，这在化合物3、5、7和9中很明显。这与我们先前知识中的结构-性能关系一致。3.3. 熵描述符的影响烷基侧链的延伸和极性官能团的引入与范德华体积（VW）的膨胀（1134.5–1872.3 Bohr3）在硝基化合物系列中呈线性相关（图4a，左侧；表S22）。分子体积的增加与熔点的降低之间的这种相关性证实了引入较长侧链或环状结构导致的体积增加破坏了分子的有序排列，显著增强了系统的无序性（熵增加），从而导致熔点降低。下载：下载高分辨率图像（968KB）下载：下载全尺寸图像图4。(a) 分子体积和分子振动模式（化合物1）以及描述熵的特征值。(b) 熔点与熵因子之间的关系。硝基化合物侧链的结构多样性在模拟条件下诱导出不同的热振动模式。在298.15K时，从头算分子动力学模拟显示这些化合物表现出平移、旋转和振动运动，柔性侧链的振动幅度大于刚性的吡唑核心（图4a，右侧；图S3；补充视频1）。观察表明，在热刺激下，分子结构的熵贡献主要来自柔性侧链。结构和能量波动是使用在1.9 × 10?–2.0 × 10? fs内采样的平衡态配置来量化的（图4a，右侧）。以下是与本文相关的补充数据：下载：下载视频（1MB）多媒体组件18. 通过RMSD_ave（ALIGN）和RMSD_ave进行的结构分析揭示了0.029–0.169 ?和2.813–4.697 ?的构象位移（表S22），而平衡态能量方差（E_var）量化了从1.23 × 10??到8.15 × 10??的热力学波动（表S27）。这些指标表明，延长的侧链放大了平衡态的结构重组和能量分散。能量松弛动力学，其特征是E_b（?0.75至?0.12）和ΔE（Emax-Eequ，其中Emax代表最大能量，Eequ代表平衡能量）（0.020–0.039 Hartree）在不稳定到平衡的转变期间（表S27），反映了更灵活的结构在达到平衡时表现出更明显的能量变化，近似于更大的熵变化。鉴于该数据集的特定分子结构特征，研究表明熵效应在降低熔点方面也起着重要作用，而不仅仅是纯焓效应。利用熵描述符集，我们建立了第二个熔点模型。与基于焓描述符的初始熔点模型相比，这个新模型表现更好。对于特征参数，基于结构变化的描述符（如RMSD_ave（ALIGN）和RMSD_ave）的拟合结果不太令人满意，这表明仅基于结构变化的程度来近似熵效应是相对粗糙的（图4b）。相比之下，进一步分析熔点（Tm）和熵描述符表明Tm与分子体积VW之间存在强线性相关性（R2 = 0.795）。同时，作为熵描述符的E_b也表现出一定的稳健性（R2 = 0.759），而能量差ΔE显示出中等的预测能力（R2 = 0.686），突出了其在量化熵贡献方面的次要作用。同样，参数化熵描述符揭示了几个显著特征：（1）仅依赖结构指标来近似熵有显著局限性；（2）分子结构中的较长侧链通常会导致结构和能量在热效应下的更明显变化，包括结构灵活性的逐渐增加，以及柔性系统中E_b和ΔE的增加，这与熵和能量重新分配之间的关系中的热力学原理一致。3.4. Tm预测模型性能为了研究分子描述符对含硝基柔性化合物熔点的影响，并确定最佳描述符组合，本研究建立了多元线性回归（MLR）和非线性k最近邻（KNN）模型。为了严格评估模型的泛化能力，我们实施了以下数据集划分策略：本研究中新合成的九种硝基化合物和我们之前研究中的三种硝基化合物构成了训练集，而从文献中收集的包含不同功能领域的柔性硝基化合物，如能量材料（ADNNMP、BOM、AMDNP、MDNNMP）和药物分子（甲硝唑、MTZ和奥尼唑、RON）构成了独立的测试集（图5c，底部）。这种策略允许通过结合结构相似和构象多样的化合物来有针对性地评估模型的预测和泛化能力。下载：下载高分辨率图像（776KB）下载：下载全尺寸图像图5。(a) 熔点预测模型的性能。(b) 熔点预测模型的SHAP分析。(c) 测试集的分子结构。在特征选择阶段，相关性分析揭示了某些描述符之间的强线性依赖性（例如，描述符BEC_max1-4和BEC_max1-5之间的相关系数R2 = 0.99），表明存在冗余（图3a）。数据通过z分数标准化进行了标准化，以提高模型的可解释性和泛化能力，并在独立的外部验证集（测试集）上进行了后续评估。为了确定平衡预测能力和泛化性能的最佳描述符组合，通过在训练集上对所有可能的描述符组合进行MLP和KNN模型的彻底训练，并使用最小平均绝对误差（MAE）作为选择标准来确定表现最佳的模型。线性模型（MLR）和非线性模型（KNN）被比较以确定它们的最佳特征组合和性能（图5a）。MLR模型使用了VW、RMSD_ave和E_b（训练：R2 = 0.75，MAE = 16.74°C；测试：R2 = 0.70，MAE = 19.04°C），而KNN模型使用了BEC_max1-3、BEC_total、VW、RMSD_ave、E_b和ΔE（训练：R2 = 0.83，MAE = 12.48°C；测试：R2 = 0.79，MAE = 15.29°C）。线性和非线性模型都强调了三个熵描述符VW、RMSD_ave和E_b的重要性，这与前一节中建立的这些熵描述符与熔点之间的强线性相关性一致。此外，在非线性模型中包含焓描述符（如BEC）显著提高了模型性能，强调了它们对熔点预测的非线性贡献。SHAP分析（图5b）显示VW和RMSD_ave与熔点呈线性反比关系，而E_b则呈直接线性关系。这一发现与分析结果一致，其中柔性基团被证明引入了更松散和更无序的分子结构，通常导致熔点降低。相比之下，BEC_max1-3和BEC_total主要以非线性方式贡献。尽管基于BEC的焓描述符没有显示出强线性相关性，但它们在非线性模型中的包含至关重要，因为它们有效地捕捉了预测熔点时分子间弱相互作用的显著作用。此外，ΔE在KNN模型中显示出线性趋势，但其包含并没有显著改善线性模型，表明其作用更多是协同的而不是独立的线性。基于平均绝对SHAP值的分析，发现描述符ΔE、E_b和RMSD_ave在最佳模型中具有最高的SHAP值（分别为6.36、6.02和4.79），表明代表分子振动能量的描述符在预测熔点中起主导作用。此外，VW（平均绝对SHAP值 = 2.63），反映了柔性基团的立体阻碍，以及分子间相互作用描述符BEC_total（平均绝对SHAP值 = 2.28）和BEC_max1-3（平均绝对SHAP值 = 2.53）对熔点预测提供了重要的补充贡献（图5b）。4. 结论为了建立一个化学上可解释且高精度的含硝基有机化合物熔点预测框架，我们构建了一个包含九种定制设计模型化合物的数据集，并提取了基于热力学的描述符。这些描述符包括与焓相关的描述符：局部参数（例如，最大二聚体结合能BEC_max1）和全局参数（例如，结合能方差BEC_var），以及通过结构动力学捕捉分子灵活性的熵描述符（例如，均方根偏差RMSD_ave）和能量松弛特性（例如，能量松弛因子E_b、能量方差E_var）。通过比较线性MLR和非线性KNN模型，确定KNN模型为最佳模型，它使用了描述符组合BEC_max1-3、BEC_total、VW、RMSD_ave、E_b和ΔE。其训练性能为R2 = 0.83和MAE = 12.48°C，测试性能为R2 = 0.79和MAE = 15.29°C。SHAP分析进一步表明，基于灵活性的熵描述符有效地捕捉了控制熔点行为的主导因素，而基于相互作用的焓描述符提供了重要的非线性补充贡献。该模型是专门为柔性硝基芳香/杂环化合物开发的。对于具有相似分子灵活性的硝基杂环化合物，当前模型可能具有合理的适用性。然而，由于其稳定性控制熵变化机制的根本差异，将其直接应用于刚性笼状能量分子可能是不可靠的。本研究建立的二元描述符框架具有内在的通用性。通过未来的改进，特别是通过调整或补充描述符集以捕捉这些系统的独特热力学和结构特征，其应用可以扩展到其他化学类别，如硝基杂环化合物或笼状分子。这项工作强调了物理化学驱动的描述符在预测建模中的价值。将数据集扩展到更广泛的化学领域，有望提高模型的精确度及其在工业和制药领域的适用性。

CRediT作者贡献声明：
冯晓宇：撰写原始稿件、数据可视化、数据管理。
陈芳：撰写原始稿件、数据可视化、实验研究、形式分析、数据管理。
宋思伟：模型验证、方法论研究、概念构建。
王瑞辉：数据可视化、软件开发、实验研究。
王毅：撰写审稿稿件、资源协调、资金筹集、概念构建。
张庆华：撰写审稿稿件、项目监督、项目管理、资金筹集。