通过自适应力匹配方法利用电子结构预测环己烯的性质

《The Journal of Physical Chemistry C》：Predicting Properties of Cyclohexene with Electronic Structure Methods through Adaptive Force Matching

　　高分辨率图像
下载MS PowerPoint幻灯片

电子结构计算是确定分子能量和结构的强大工具；然而，其计算成本限制了它们在大型系统中的应用，并阻碍了需要集合平均值的有限温度性质的计算。我们证明了自适应力匹配（AFM）方法可以用来弥合这一差距，仅基于电子结构信息就能预测分子的宏观性质。作为概念验证，我们将研究重点放在环己烯上。通过将AFM模型与使用MP2和B3LYP-D3方法计算的参考力进行拟合，我们成功预测了16种不同的性质。对于所有有可用实验值的性质，AFM模型都取得了非常好的吻合度。对于其他性质，AFM模型的预测填补了缺失实验数据的空白。虽然分别拟合纯相和水合相的模型能够非常吻合地预测环己烯的溶解度，但对于某些性质（如环己烯和水之间的界面张力），则需要一个混合相模型。这些结果表明，使用基于AFM的分子动力学直接从电子结构信息预测广泛的宏观性质是可行的且准确的，为计算材料科学开辟了新的途径，并加速了新型化学物质的发现和设计。

**I. 引言**
电子结构计算已成为确定小分子构象能量和结构信息的强大工具。(1?3) 然而，对于大型系统来说，精确的电子结构计算可能会变得极其昂贵。即使是计算效率更高的密度泛函理论（DFT）方法，通常也仅适用于包含不超过几百个原子的系统。(4?6) 电子结构方法面临的另一个挑战是，大多数材料性质（如液体的自由能、密度和沸点）都是在有限温度下的集合性质。为了准确确定这些性质，需要执行集合平均值，而这不仅仅依赖于关于少数几种构象的信息。因此，需要一种更高效的模型来预测这些宏观性质。

应对这一挑战的一种方法是将量子力学（QM）信息映射到分子力学力场中。(7?10) 这些力场随后可以与分子动力学或蒙特卡罗采样等方法结合使用，以计算有限温度下的宏观性质，从而实现集合性质的预测。尽管从电子结构计算中推导出力场有着悠久的历史，(11?16) 最大的进展是在水建模方面取得的，(17,18) 例如MB-Pol (19) 和各种基于自适应力匹配（AFM）的方法(20) 在预测水的宏观性质方面表现出良好的准确性。(21?24) 相比之下，为其他材料（如碳氢化合物）开发纯基于电子结构的模型仍然较为有限。尽管最近的AFM模型在预测密度和汽化热等特定性质方面显示出潜力，(25,26) 但尚不清楚这些仅通过电子结构计算开发的模型是否能够准确预测给定材料的广泛物理性质。通过这种方法全面预测多种物理性质的能力尚未完全确立。

在这项概念验证研究中，我们选择了环己烯作为研究对象。作为尼龙工业合成中的中间体，环己烯具有重要的实际意义。其分子结构包含一个双键，导致相对紧张的构象。此外，环己烯的特点是分子间相互作用较弱，其相互作用能量主要由色散力主导。从基于电子结构的模型的角度来看，这类色散作用系统面临两个挑战。第一个挑战是众所周知的DFT难以准确模拟动态相关性，(27?29) 这导致了色散现象。另一个挑战是这类弱结合系统的凝聚能较小；纯环己烯的凝聚能的1 kcal误差可能会转化为较大的百分比误差。这种误差预计会影响密度预测，从而导致密度的显著误差。例如，在基于AFM的CO2模型中，由于基组叠加误差引入的微小误差导致了密度的显著误差。(30) 虽然AFM在模拟极性系统（如水、(22,23,31) 水合离子、(32,33) 和水合小分子,(34,35) 方面取得了良好的成功，但对于像环己烯这样的弱相互作用化合物来说，它是评估我们方法性能的一个好案例。

事实上，已经为纯相和水相中的环己烯创建了AFM模型。(25) 然而，尚未对这些模型在广泛性质上的性能进行全面评估。在这项研究中，我们将评估纯相模型预测一系列物理性质的能力，包括密度（ρ）、汽化焓（ΔHvap）、表面张力（γ）、粘度（η）、扩散常数（DP）、等温压缩性（κT）、汽化自由能（ΔGvap）、临界温度（TC）、压力（PC）和密度（ρC）以及沸点（Tb）。

除了研究纯环己烯的性质外，我们还探讨了环己烯在水合环境中的行为。值得注意的是，AFM模型不依赖于组合规则，这意味着不能从纯环己烯和纯水的参数推导出水合环境中的环己烯参数。相反，必须开发一个单独的模型来描述水合环境中的环己烯行为。这种水合模型能够可靠地预测水合自由能（ΔGhyd）、水合焓（ΔHhyd）和扩散常数（DH）。通过确定ΔGhyd和ΔGvap，我们可以预测环己烯的溶解度。然而，虽然可以通过分别拟合纯相和水合相的模型来预测溶解度，但这种方法无法确定环己烯-水混合物的所有性质。例如，确定环己烯和水之间的界面张力（γint）需要一个能够描述两相的单一模型。为了解决这个问题，我们通过同时拟合纯相和水合相的环己烯来开发了一个混合相模型。该模型使我们能够预测环己烯-水混合物的γint，并且我们还报告了它在模拟纯环己烯液体方面的性能。

**II. 计算细节**

**A. 模型开发**
AFM已被证明是将高级电子结构势能表面映射到高效力场模型的强大工具。(33,36,37) AFM通常独立于其他非键合参数来拟合色散力。通过将相对较弱但重要的色散相互作用与库仑力和短程排斥力解耦，我们可以开发出更稳健的模型。在拟合色散力之后，AFM过程包括三个主要步骤：采样步骤、量子力学（QM）或QM/分子力学（MM）步骤以及力匹配（FM）步骤。

在采样步骤中，通常会执行分子动力学（MD）模拟，以生成代表感兴趣构象空间的训练集。QM或QM/MM步骤涉及使用参考电子结构方法计算力，这为确定力场参数提供了必要的数据。实际参数的确定在FM步骤中进行，通常包括多个阶段。例如，通常先拟合分子间参数，然后再拟合分子内参数，这减少了拟合强耦合参数的数值挑战。

对于环己烯，之前已经为纯相和水合相分别开发了AFM模型。(25) 纯相模型是在不使用偏电荷的情况下，根据使用MP2方法计算的参考力进行拟合的，色散相互作用使用对称自适应扰动理论（SAPT）(38) 以MP4的质量进行计算。对于像纯环己烯这样的弱相互作用系统，MP2的准确性很重要，而水合相模型则是根据使用计算要求较低的计算方法B3LYP-D3(BJ) 计算的力进行拟合的。(39?42) AFM的计算成本主要集中在QM或QM/MM计算上。使用AMD EPYC 9454 CPU上的ORCA软件，每个MP2计算需要45分钟，该CPU使用12个核心。使用相同的硬件，每个B3LYP-D3(BJ)计算大约需要15分钟。为了生成每个纯相和水合模型的所有参考力，大约需要进行800次这样的单点力计算。由于这些模型的开发过程已在之前描述过，我们将在支持信息中仅简要总结所使用的程序。

虽然使用纯相和水合相的两个单独模型可以确定包括环己烯在水中的溶解度在内的广泛性质，但它不允许模拟水中的环己烯液滴或计算两种液体之间的界面张力。为了便于混合物的模拟，必须开发一个能够在一次模拟中描述两相的单一模型。这需要根据纯相和水合相的参考力来拟合模型。值得注意的是，环己烯液体和水合环己烯的AFM模型使用不同的偏电荷。由于环己烯是非极性的，之前开发的纯环己烯模型在拟合时没有使用偏电荷，而水合环己烯模型则使用了偏电荷。为了模拟与液态水接触的环己烯，环己烯必须使用一组偏电荷。虽然有一种可能性是拟合一组适用于纯环己烯和水合环己烯的电荷，但我们选择了更简单的方法：在混合相模型中也使用水合环己烯模型的偏电荷。我们认为，对于非极性系统来说，实际偏电荷的重要性较低，正如我们在不使用偏电荷的情况下成功开发这些模型所证明的那样。

另一个复杂因素是，纯相模型和水合模型采用了不同的原子类型方案，后者使用了更多的原子类型，如图1所示。由于排斥力和库仑项高度耦合，具有不同偏电荷的两个原子应该使用不同的排斥项进行拟合。在之前的纯相模型开发中放弃使用偏电荷，使其能够仅使用两种碳原子类型和两种氢原子类型，而水合模型使用了三种碳原子类型和三种氢原子类型。混合模型将使用与水合模型相同的原子类型。

**图1**
图1. (a) CHEn和(b) CHEh模型所使用的原子类型。

**高分辨率图像**
下载MS PowerPoint幻灯片

为了更好地理解混合相模型的性能，我们在图2中展示了各个相模型和混合相模型的分子间力、扭矩以及原子力散点图。对于各个相模型，纯相QM参考值是根据纯环己烯（CHE）模型CHEEn进行拟合的，而水合相QM/MM参考值是根据水合模型CHEh进行拟合的。混合相模型CHEm则是根据两相的参考计算进行拟合的。完美的拟合是对角线上的直线。从图2中可以清楚地看到，比较单个模型和混合模型时，拟合质量没有明显差异。

**图2**
图2. 比较模型性能与所拟合的参考QM力和扭矩的散点图。(a) 显示了CHEm模型（红色）和CHEh及CHEEn模型（黑色）的分子力，(b) 分子扭矩，以及(c) 原子力。

两种模型的均方根误差（RMSE）在表1中报告，同时提供了拟合力和参考力之间线性回归的斜率和R2值。分子力和扭矩的RMSE值反映了模型再现参考力分子间组分的能力。有趣的是，CHEm模型在再现分子力和扭矩方面的表现比使用单独模型大约好2%。使用AFM时，环己烯-水和水环己烯-环己烯相互作用的非键合短程参数是独立的。因此，将纯相参考数据添加到训练集中不会影响非键合的环己烯-水参数，因为CHEm模型使用与CHEh相同的偏电荷。因此，两个模型得出了相同的环己烯-水相互作用结果。观察到的分子间RMSE略有改善，这归因于在纯相中分子间相互作用的再现性得到了提高。这种改善表明，即使是在水环境中拟合的部分电荷，也能使模型性能略优于不使用部分电荷的模型。表1显示了模型与QM参考值之间力或扭矩的RMSE、斜率和R2值，分别对应于将两种相拟合到各自的模型中，以及将两种相拟合到CHEm模型中。实验和预测的环己烯性质

| 性质 | 实验值 | 预测值 | ΔGhyd | ΔHhyd | ΔHn/a | 密度 | ΔHvap | ΔGvap | DPn/a | κT | γ | η | Tb | ΔGsol | γint |
| ---- | ------ | ------- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- | ------- | ------ |
| (H) | 3.2 ± 0.3 kJ/mol | 3.8 ± 0.3 kJ/mol | -32.20 kJ/mol | (0.99 ± 0.01) × 10–5 cm2/s | 811. g/L | 33.47 kJ/mol | 18.911 kJ/mol | (2.4 ± 0.2) × 10–5 cm2/s | 1.0366 × 10–3 MPa–1 | 26.17 mN/m | 0.625 mPa s | 560.45 K | 44.3 bar | 276.7 g/L |
| (M) | 3.8 ± 0.3 kJ/mol | -31.6 ± 1.0 kJ/mol | (0.99 ± 0.07) × 10–5 cm2/s | 839 ± 2 g/L | 32.9 ± 0.1 kJ/mol | 17.4 ± 0.3 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 22.6 ± 2.0 mN/m | 0.49 ± 0.01 mPa s | 513 K | 49 bar | 283.25 g/L |
| (N) | 3.2 ± 0.3 kJ/mol | -31.8 ± 0.8 kJ/mol | (0.99 ± 0.01) × 10–5 cm2/s | 851 ± 1 g/L | 34.5 ± 0.3 kJ/mol | 18.2 ± 0.2 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 25.4 ± 0.7 mN/m | 0.66 ± 0.02 mPa s | 344 K | 303.9 g/L | 341 K |
| (H,N) | 1.51 kJ/mol | 3.8 ± 0.3 kJ/mol | -32.20 kJ/mol | (0.99 ± 0.01) × 10–5 cm2/s | 811. g/L | 32.9 ± 0.1 kJ/mol | 17.4 ± 0.3 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 26.17 mN/m | 0.49 ± 0.01 mPa s | 560.45 K | 44.3 bar | 283.25 g/L |
| (M,N) | 1.51 kJ/mol | 3.8 ± 0.3 kJ/mol | -31.8 ± 0.8 kJ/mol | (0.99 ± 0.07) × 10–5 cm2/s | 851 ± 1 g/L | 34.5 ± 0.3 kJ/mol | 18.2 ± 0.2 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 25.4 ± 0.7 mN/m | 0.66 ± 0.02 mPa s | 344 K | 303.9 g/L | 341 K |
| (B) | 1.51 kJ/mol | 3.8 ± 0.3 kJ/mol | -31.6 ± 1.0 kJ/mol | (0.99 ± 0.01) × 10–5 cm2/s | 811. g/L | 32.9 ± 0.1 kJ/mol | 17.4 ± 0.3 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 22.6 ± 2.0 mN/m | 0.49 ± 0.01 mPa s | 560.45 K | 44.3 bar | 283.25 g/L |
| (C) | 1.51 kJ/mol | 3.8 ± 0.3 kJ/mol | -31.8 ± 0.8 kJ/mol | (0.99 ± 0.01) × 10–5 cm2/s | 811. g/L | 32.9 ± 0.1 kJ/mol | 17.4 ± 0.3 kJ/mol | (1.9 ± 0.3) × 10–5 cm2/s | (1.0 ± 0.1) × 10–3 MPa–1 | 25.4 ± 0.7 mN/m | 0.66 ± 0.02 mPa s | 560.45 K | 44.3 bar | 283.25 g/L |

注：部分数据在中文中可能具有多重含义，因此在翻译时进行了适当的解释和调整。例如，“ΔGhyd”在某些上下文中表示“水化自由能变化”，而在其他上下文中可能表示“总自由能变化”。此外，一些单位（如kJ/mol和cm2/s）在中文中直接保留，因为它们在国际上被广泛接受。