分子动力学方法——分子的瞬时频率模式，用于评估经典水力场和刚性水力场的性能

《The Journal of Physical Chemistry B》：Molecular Dynamics-Instantaneous Frequencies of Molecules Method Modes to Address the Performance of Classical and Rigid Water Force Fields

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

分子力场（FFs）的准确参数化对于从分子动力学模拟中获得可靠的描述符至关重要。对FFs进行基准测试需要能够直接或间接检测化学环境变化引起的相互作用势能变化的方法。本研究使用“分子动力学-分子瞬时频率”（MD-IFM）方法作为FF基准测试的工具，该方法允许我们从分子动力学模拟中计算振动光谱可观测量，包括中心频率（ω0）、频率波动幅度（σ）和光谱扩散相关时间（τc）。我们通过评估常见水分子力场（TIP3P、SPC/E、TIP4P/Ew、TIP5P和OPC）在高阶MB-pol势下的振动描述符的有效性，以及这些描述符在解析溶剂化和动力学方面的应用来展示该方法的优势。虽然经典力场在预测水的振动频率方面表现良好，但它们系统性地低估了局部环境的异质性（σ），并且与MB-pol相比显示出更快的光谱扩散动力学（更低的τc值）。分析表明，这些不足与对氢键网络的简化表示有关。此外，我们的结果验证了MD-IFM作为一种高效且稳健的FF验证方法的有效性。

引言
准确描述分子系统的结构和动态行为对于理解其宏观性质（如密度、沸点和熔点、导电性等）至关重要。(1,2) 因此，分子结构、溶剂化动力学和扩散动力学的表征非常重要。例如，电池电解质中的缓慢扩散动力学和强分子间相互作用通常会导致高粘度系统，从而具有低导电性和较差的传输性能。(3) 另一个例子是酶，它们需要活性位点中氨基酸的精确排列来结合底物并稳定过渡态，因为底物方向的微小变化会显著改变催化过程的反应动力学。(4,5) 此外，复杂溶液（如深共晶溶剂和离子液体等）的多组分性质可能导致纳米异质性的存在，从而影响溶质的溶解能力。(6) 实验上，广泛使用光谱技术来研究复杂系统的分子结构（如散射技术）和动力学（如时间分辨技术）。(7?22) 然而，这些实验技术得到的信号和指标通常难以解释。因此，通常将计算建模与实验结合使用，以揭示系统的分子尺度行为。(23?57) 为了生成预测性的分子动力学模拟，需要准确描述势能面（PES）。为此，PES可以通过从头算方法、半经验方法或经典力场（FFs）构建。

从头算和半经验方法能提供更准确的系统描述，但计算成本比经典力场更高。例如，即使使用超级计算机，这些计算也只能适用于几百个原子，并且最佳情况下的时间尺度限制在皮秒级别。(58) 相比之下，经典力场给出的分子系统表示更快且计算效率更高，因此在现代的大规模分子模拟中得到广泛应用。(59) 然而，这些方法正确表示PES复杂性的能力在许多情况下是有限的。因此，需要一种能够对FFs进行基准测试的方法，以评估分子表示的适用性，从而确保可以自信地使用FFs生成分子动力学模拟，并从中提取信息以解释实验结果。

存在不同的方法来基准测试计算方法；一些例子包括计算相变或溶剂化的能量，并将其与实验或更高层次的理论进行比较。(60?64) 通常，基准测试会检查FF产生的PES的形状（即最小值和凹度）。(65?67) 然而，这些特性也可以从由于分子间相互作用而在分子中观察到的PES波动中得出。(68) 通过分子动力学（MD）模拟，可以经典地探索PES的形状和波动。因此，基准测试方法可以通过其化学环境的动态波动间接检查PES的形状。为此，能够表征PES形状和波动的可观测量是必不可少的。

振动可观测量，如中心频率（ω0）、频率波动幅度（σ）和来自振动光谱扩散的去相关时间（τc），是描述PES特性的有价值指标。前两个指标提供了关于分子及其所经历的分子环境的平均结构的信息，而后者通过振动频率相关性的丧失描述了系统的动态。(54,70) 实验上，这些量可以通过结合线性和非线性光谱技术（如傅里叶变换红外光谱（FTIR）和二维红外光谱（2DIR）或三脉冲光子回波方法轻松获得。(69) 因此，振动光谱及其可观测量使得实验和分子模拟之间可以进行一一对应的比较。

频率图(71)已被广泛用于计算模拟振动可观测量。然而，频率图通常需要大量的参数化，并且它们在不同溶剂和更异质环境中的适用性并不总是有保证的。(72?74) 其他方法(38,54,75?77)可用于将分子坐标转换为振动可观测量，从而计算光谱线形、峰位移或频率-波动相关函数（FFCF）。(78) 然而，这些方法往往会增加计算成本。在此，我们使用“分子动力学-分子瞬时频率”（MD-IFM)(58,79)作为FFs的基准测试工具。MD-IFM是一种无需参数化的方法，用于从分子动力学模拟中计算上述振动光谱可观测量。除了避免复杂的参数化外，MD-IFM的计算成本较低，它结合了经典分子动力学和频率计算，并采用了半经验级别的理论GFN2-xTB。(80) 相比之下，使用DFT进行相同计算的速度要慢一千多倍。(79,80) 此外，MD-IFM已被证明能够使用这种计算框架准确再现N-甲基乙酰胺的酰胺I模式以及氯仿中的C-H伸缩振动。(58,81) 鉴于MD-IFM的准确性依赖于对PES的适当描述，它成为评估水分子力场性能的独特工具。

涉及水的分子模拟很常见，为此开发了许多FFs，包括常见的经典水模型，如3位点模型（TIP3P (82) 和 SPC/E (83)），这些模型为氧原子和氢原子分配了固定电荷。这些经典FFs在计算效率和准确性之间取得了良好的平衡(84)，但没有考虑对外部电场（极化）响应的电荷分布变化。此外，更复杂的FFs（如TIP4P/Ew (85,86)、OPC (87) 和 TIP5P (88)）在“虚拟”位点引入了额外电荷，以更好地捕捉静电性质和多极效应，并调节交换相互作用以提高交换和色散能量计算的准确性。(89) 其他FFs（如AMOEBA (90,91)）使用诱导点偶极子和多极表示来引入极化效应。例如，CHARMM Drude模型(92)也为原子添加了Drude粒子以达到相同的目的。TIP4P/PQ (93) 和 POL3 (94)也是特定于水的可极化模型，广泛应用于具有高电子环境的系统，如离子溶剂或水-金属界面。(95,96) 虽然可极化模型考虑了电荷分布对外部电场的动态响应，并提供了比固定电荷模型更准确的化学模拟，但由于计算成本较高，它们并未被广泛采用。

最近，MB-pol被开发出来以准确表示水。MB-pol是一种特别重要的力场(97?99)，可能是最接近通用模型的选择，因为它能够在气相和凝聚相中准确再现水的性质(100)，而不需要依赖单独的模型（例如用于冰模拟的TIP4P/Ice (101)）。MB-pol的多体势能特性使其能够准确描述许多结构和热力学性质，其中一些具有亚化学精度。(102) 尽管神经网络势能是一个有前景的参考点(103,104)，并且从头算方法总体上能获得更好的结果，但由于其对水簇振动光谱、液态水的红外和拉曼光谱、H2O和HDO水探针的红外带以及2DIR光谱中OD伸缩模式的光谱扩散过程的广泛验证，我们选择MB-pol作为评估其他水模型的金标准。(47,50,52,105?107) 然而，像可极化FFs一样，MB-pol也仅适用于相对较短的分子动力学模拟，这使得描述水中的大分子系统（如蛋白质）在数百纳秒内的模拟变得极其昂贵。(108,109)

在这项工作中，我们的目标是展示MD-IFM方法作为一种工具，用于研究在由不同FFs（TIP3P (82)、TIP4P/Ew (85,86)、TIP5P (88)、SPC/E (83) 和 OPC (87)）表示的显式水分子浴中的HDO分子的上述振动可观测量，同时使用MB-pol作为描述水溶剂化和动力学的金标准。

方法
MD-IFM方法首先通过经典分子动力学探索配置空间，然后使用一种已被证明能够足够准确计算振动频率的半经验方法逐一处理大量溶剂配置的振动频率。为了保留PES探索的结果，对FF描述的化学环境内的分子探针进行几何优化。所有选择（如簇大小、样本长度等）都是在指标收敛后仔细挑选的，以便将获得的结果归因于FF的性能而非其他来源。

对于每个选定的FF，使用AMBER24软件进行了分子动力学（MD）模拟。(110) 使用PACKMOL软件构建了边长为40 ?的水盒子。(111) 在NVE系综中从0 K加热到300 K之前进行了能量最小化；然后在NPT系综中进行了1 ns的平衡步骤，随后进行了12个60 ps的NVE系综的生产步骤。生产运行的最后50 ps用于后续的频率计算，每50 fs提取一个快照，总共12,000个快照。所有阶段的牛顿方程积分时间步长设置为2 fs。使用SHAKE来约束涉及氢原子的键。对于MB-pol，每48 fs提取一个快照。

频率计算
对于频率计算，通过选择固定水分子周围前三层溶剂壳层内的所有水分子来提取一个簇，该水分子作为分子探针，这是基于径向分布函数（图S1），从而得到一个包含30个分子的簇。我们还计算了每个系统中选定探针的第一层溶剂壳层中的水分子数量；为此，径向分布函数的第一个最小值被用作距离截止值。仅优化中心分子探针（其中一个氢原子被氘替换），并使用GFN2-xTB算法进行频率计算，如先前报道的那样(58)，同时保留了FF描述的簇中其余水分子的采样几何结构以及这种环境对HDO分子振动频率的影响。

HDO探针的每个振动模式的ω0值是通过频率的强度加权平均值获得的。频率分布也是使用得到的12,000个频率值构建的，σ则是从频率的强度加权标准差中获得的。FFCF是通过频率随时间的自相关性计算得出的（C(t) = t)δω(0)），其中δω(t)表示时间t时的频率与轨迹上平均频率之间的差异，遵循Kubo形式主义。（69）FFCF采用双指数衰减模型进行建模（C(t)=A1e?t/τ1+A2e?t/τc+y0）（54），以获得每个正常模式的τc作为最长特征时间。不确定性计算的标准误差直接来自拟合参数，或者根据SE=s√N计算得出，s是标准差或加权标准差（对于中心频率）。报告的不确定性值是基于95%置信区间的1.96倍标准误差。

结果与讨论
通过MD模拟评估了水（HDO）分子平均观察到的分子环境及其动态，使用了不同的FFs，这些FFs基于中心频率（ω0）、频率波动幅度（σ）和振动光谱扩散的去相关时间（τc）。这些值是为HDO分子的三种振动模式计算的：HOD弯曲、OD伸展和OH伸展。MD-IFM方法对每种振动模式的振动可观测量得出的结果在表1中呈现，每个单独的指标在以下部分中分别讨论。

表1. 通过MD-IFM方法获得的ω0、σ和τc的一般结果，用于评估的FFs和HDO振动模式
力场 HOD弯曲 OD伸展 OH伸展
ω0 (cm–1) σ (cm–1) τc (ps) ω0 (cm–1) σ (cm–1) τc (ps) ω0 (cm–1) σ (cm–1) τc (ps)
TIP3P 1350.2 ± 0.5 26.0 ± 0.3 0.23 ± 0.0 252 0 ± 27 0.39 ± 0.0 346 8 ± 29 0.43 ± 0.0
SPC/E 1339.8 ± 0.5 27.4 ± 0.5 0.86 ± 0.0 247 3 ± 28 0.80 ± 0.0 340 6 ± 21 15 ± 7 0.65 ± 0.0
TIP4P/Ew 1348.5 ± 0.5 25.5 ± 0.6 0.40 ± 0.0 249 8 ± 27 0.80 ± 0.0 341 8 ± 21 0.83 ± 0.0
OPC 1341.2 ± 0.5 24 ± 10 4.4 ± 0.0 255 16 ± 20 0.64 ± 0.0 350 7 ± 28 40.8 ± 0.0
TIP5P 1346.1 ± 0.5 25.8 ± 0.6 0.95 ± 0.0 248 8 ± 28 69 ± 0.0 342 6 ± 21 10 1.09 ± 0.0
MB-pol 1340.9 ± 0.6 28.9 ± 0.7 0.51 ± 0.0 247 28 51 2 ± 0.0 340 3 ± 21 15 8.13 ± 0.0

中心频率（ω0）考虑了化学环境及其平均结构的影响，使我们能够评估FF所做的整体PES描述。HOD弯曲模式的ω0值范围在约1339到约1350 cm–1之间（表1），显示出所有FFs的相似性能。然而，这些值与报道的实验值1458 cm–1相比被低估了。（32）此外，SPC/E和OPC得出的值与MB-pol非常接近（图1A），而其他FFs的值相差很小，介于0.4%到0.7%之间。对于OD伸展模式，ω0值的范围更分散，大约在约2470 cm–1到约2551 cm–1之间（表1）。虽然OD伸展在预测的中心频率上比HOD弯曲模式有更大的分散，但值仍然在相对较小的范围内（即，误差<3.3%）。然而，在这种情况下，所有经典FFs都与MB-pol有显著偏差，除了SPC/E。后者的中心频率误差仅为0.1%，而其他FFs的误差高达3.3%（图1A），但相对于实验值约2500 cm–1，（112）TIP4P/Ew产生了最接近的结果。OH伸展模式显示出与OD伸展预期的频率偏移，其中心频率范围在约3403 cm–1到约3507 cm–1之间，但与OD伸展模式的偏差相似。相对于MB-pol的百分比差异再次在0.1%到3.0%之间变化，SPC/E显示出最低的偏差，OPC显示出最大的偏差（图1A），这次MB-pol的结果最接近实验值约3400 cm–1。（55,112）值得注意的是，SPC/E的出色性能可能与其有效的极化校正有关，（83）该校正考虑了水分子的更真实的偶极矩，从而可能更好地表示水的平均结构，因此中心频率更为准确。

图1
（A）使用经典FFs预测的ω0与MB-pol相比的百分比差异。差异大于零表示MB-pol得出的值较低，反之亦然。穿过x轴的误差条表示值之间的差异可以忽略不计。
（B）使用经典FFs预测的ω0与使用MB-pol预测的ω0的相关性图。黑线代表y = x函数，垂直虚线是MB-pol FF的振动模式位置。

通过直接比较不同经典FFs对HDO不同模式的预测中心频率与MB-pol（图1B），可以更好地了解它们性能的差异。在这种表示中，较大的误差表现为距离对角线较远。结果显示，所有FFs的偏差都很小，其中OPC的偏差较小但更明显。此外，OD和OH伸展频率的预测不如弯曲模式的预测精确。总体而言，所有FFs产生的HDO振动中心频率与MB-pol相比具有合理的准确性。

如前所述，OD和OH伸展模式的中心频率不如同一分子的弯曲模式的预测精确。这种差异可以归因于每种水振动模式对化学环境变化的敏感性不同，或者等效地，水分子的氢键状态。（54,55）例如，氢键数量较少或处于气相的水分子的振动频率比处于四面体氢键网络中的水分子表现出蓝移。（113）因此，FF对分子相互作用的精确描述，即氢键，在这种情况下，应该导致计算出的中心频率有较小的偏差。这也为评估FF的质量提供了一个可靠的指标，与黄金标准进行比较。为了评估溶剂壳对每种FF的振动模式频率的影响，计算了每个HDO模式的中心频率作为第一溶剂壳中水分子数量的函数（所有FF的选择值为3、4和5，SPC/E也为6，因为其数量与3相当，图S2）。

图2中呈现的中心频率结果显示，经典FFs再现了振动频率随第一溶剂壳中水分子数量的预期变化（溶剂色变效应），其中水分子数量越多，形成氢键的概率越高。在HOD弯曲模式的情况下，TIP4P/Ew、OPC、TIP5P和MB-pol FFs显示中心频率随水分子数量连续减少，而SPC/E和TIP3P则总体上减少但不是单调的。此外，中心频率的范围较小（约3 cm–1），这与弯曲模式对化学环境的敏感性较低一致。（32）对于HDO的两种伸展模式，中心频率显示出预期的较大变化范围（约25 cm–1），这与这些模式对氢键环境的较高敏感性一致。值得注意的是，对于HDO伸展，所有FFs都一致地再现了中心频率随水分子数量增加的预期红移（图2B,C），表明大多数FF正确地表示了水分子的氢键状态对中心频率的影响。

图2
（A）HOD弯曲模式，（B）OD伸展模式，（C）OH伸展模式的第一溶剂壳中水分子数量的ω0值。

频率波动幅度（σ）
频率分布及其形状有助于描述MD模拟中HDO分子观察到的构象空间的异质性。由于σ代表频率分布的标准差，较高的σ值表示频率值的方差较大，这与探针化学环境的可变性密切相关。这些值可以在一定程度上与IR带的线宽进行比较，或者根据Kubo形式主义从FFCF的预指数因子得出。（69）
结果显示，HOD弯曲模式的分布接近单一高斯分布（图S3和表S1），表明这种水模式的频率波动具有随机行为，而伸展模式的分布（图S4和S5以及表S2和S3）不能用单一高斯表示，表明振动动态更为复杂。（114）具体来说，HOD弯曲模式的频率波动幅度（σ）在所有FFs中范围约为24到29 cm–1（表1）。这些σ值与MB-pol的值（图3A）显著不同。值得注意的是，MB-pol呈现了HOD弯曲最宽的频率分布，而OPC显示最窄的分布，其余的分布在中间。与MB-pol相比，SPC/E的描述最接近，差异仅为?5.1%。所有其他FFs的偏差都很大，最高达到?16.4%。在OD伸展模式的情况下，σ值要大得多，范围从63 cm–1到82 cm–1。然而，相对于MB-pol的百分比差异与弯曲模式相似（即?0.3%到?23.8%）。此外，SPC/E和TIP5P与MB-pol的差异不大，而OPC由于频率分布较窄，偏差最大。最后，OH伸展模式显示出最宽的分布，σ = 85 – 115 cm–1，与MB-pol的百分比差异范围从?0.3%到?26.5%，并且与OD伸展模式的模式类似；即SPC/E、TIP4P/Ew和TIP5P最接近MB-pol，而OPC由于其狭窄的分布具有最大的差异。

图3
（A）使用经典FFs预测的σ与MB-pol相比的百分比差异。差异大于零表示MB-pol得出的值较低，反之亦然。穿过x轴的误差条表示值之间的差异可以忽略不计。
（B）使用经典FFs预测的σ与使用MB-pol预测的σ的相关性图。黑线代表y = x函数，垂直虚线是振动模式的位置。

图3B中的相关图总结了所有经典FFs与MB-pol相比的分布宽度结果。从这些结果可以看出，所有FFs预测的频率分布都比MB-pol更窄。虽然这种低估在两种HDO伸展模式中更为明显，但与MB-pol的误差在统计上对于所有模式来说都是相似的（图3A），除了TIP3P和OPC的结果有统计学上的差异。

所有FFs在频率分布上的一致结果是OPC和TIP3P倾向于为HOD伸展模式产生更窄的分布。分析HDO分子的第一溶剂壳中的水分子数量（图S2）显示，OPC FF产生的环境对于水分子的氢键环境稍微不那么异质。然而，TIP3P的情况并非如此，因为它具有更宽的氢键状态分布（图S2）。因此，从OPC-MD计算出的振动频率分布具有较小的波动，因为HDO分子在模拟中经历的化学环境较少多样（图S2）。相比之下，TIP3P的情况不能简单地用相同的氢键论点来解释（图S2），但它可能与HDO分子在模拟中观察到的较少异质的环境有关。虽然像TIP3P和TIP4P/Ew这样的模型使用接近实验气相值的O–H键长，但OPC使用了显著更短的键长0.8724 ?和更窄的H–O–H角103.6°。这种水分子结构旨在优化体介电性质的分子四极矩（87），但可能会改变探针周围的局部静电环境。因此，OPC是一种针对宏观整体性质优化的力场（FF），它在描述局部溶剂化环境的微观异质性时牺牲了精度。通过评估HDO第一溶剂化层中水分子数量变化对频率分布的影响（使用与ω0分析相同的值），也研究了溶剂化结构对频率分布形状的影响。在HOD弯曲模式中，潜在的分布不仅可以用单一的高斯函数很好地表示，而且它们的最大值也大致位于相同的频率上（见图4A和表S4）。这一观察结果与中心频率结果以及之前报道的水弯曲模式的缺乏溶剂化色变现象一致。（32）此外，它们的不对称系数（ACs）作为偏度的代表，相对较小，这与它们的高斯形状相符。相比之下，HDO的拉伸模式的潜在分布具有明显的红移（见图4、表S6和S8），这直接与拉伸模式对第一溶剂化层中水分子数量的敏感性有关。另外，所有子分布都是不对称的，并且表现出负偏度，这一点通过它们的ACs可以看出来（见表S5、S7和S9）。有趣的是，HDO的拉伸模式的ACs明显大于其弯曲模式的ACs。此外，在所有情况下，ACs与MB-pol相当，除了OPC，其ACs较小。再次强调，频率子分布中的大不对称性反映了氢键状态对HDO分子的影响。当HDO分子的第一溶剂化层中有三个水分子时，这种效应在OD和OH拉伸模式的分布中可以清楚地看到（见图4）。在这种情况下，所有频率分布的偏度都明显高于那些对应于更多溶剂化水分子的情况。从分子角度来看，这种大的不对称性可能是由于少数水分子在空间上的稀疏分布造成的，因为这个溶剂化层中的微小几何变化应该会导致显著的频率波动。这里观察到的频率分布的不对称性与之前的研究结果一致（32,113），这些研究将这种不对称性归因于处于不同氢键状态的水分子群体。值得注意的是，所有FF和MB-pol都能一致地再现分布中的红移和不对称性，这强烈验证了它们的存在。

图4：使用TIP3P（A）、SPC/E（B）、TIP4P/Ew（C）、OPC（D）、TIP5P（E）和MB-pol（F）进行的模拟得到的HOD弯曲模式、OD拉伸模式和OH拉伸模式的频率分布。图中还展示了与第一溶剂化层中不同数量的水分子相关的子分布的贡献。

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

与之前的可观测量（ω0和σ）不同，振动光谱扩散及其相关特征时间取决于分子环境的动态。因此，这一指标提供了分子动力学（MD）模拟中结构变化的时间尺度估计，并且可以通过2DIR光谱使用CLS等方法实验测量得到。（78,115）换句话说，光谱扩散的动态描述了液体中发生的潜在分子过程的时间演变，这些过程导致状态记忆随时间的丧失，或者等效地，描述了振动频率的去相关化动态。（69）在这里，光谱扩散代表了由于溶剂化层的时间演变而引起的水振动频率波动的自相关。对于水来说，已经报告称振动频率相关性的丧失与飞秒时间尺度（50–400 fs）内的快速氢键断裂和形成以及皮秒时间尺度（0.7–1.2 ps）内的缓慢氢键网络重组有关。（54,55,116）这一分子现象被MB-pol成功再现。（50）在这里，光谱扩散的动态是通过FFCF计算的，其中每个FFCF都用双指数衰减函数来描述（见图S6–S8和表S10–S12），以获得特征时间。所有经典的FF（表1）都定性地再现了水的预期双指数衰减行为（R2 > 0.97，见图S6–S8和表S10–S12），这与水光谱扩散的快速和慢速成分一致。（13?22）然而，这里的分析仅关注皮秒成分，因为它可以直接与从2DIR光谱实验获得的实验值进行比较。图5A显示了氢键网络重组的τc值与MB-pol相比的百分比差异，对于不同的经典FF和HOD振动模式。对于HOD弯曲模式，OPC得到的τc与MB-pol最为接近，误差为?13.0%，而其他FF的差异较大，范围从?54.9%到86.3%。同样，两种拉伸模式的皮秒去相关时间与MB-pol相比的偏差范围从?18.7%到?67.9%。然而，在拉伸情况下，所有FF的氢键网络重排的特征时间都观察到比MB-pol小。此外，TIP4P和TIP5P始终比其他FF得到更接近的值。值得注意的是，在MD模拟过程中，HDO分子的两个氢原子被同等对待，因此预期OD和OH拉伸模式会表现出相同的动态行为。然而，并非所有FF都是如此，只有TIP3P、TIP4P/Ew和MB-pol的τc值相当，这表明它们的构象空间探索可能不够充分。考虑到它们使用了250 ps的动态进行构象空间采样，这本应足以描述水光谱扩散的皮秒时间尺度，这一结果令人惊讶。

图5：（A）使用经典FF预测的τc与MB-pol相比的百分比差异。差异大于零表示MB-pol得到的值较低，反之亦然。穿过x轴的误差条表示值之间的差异可以忽略不计。（B）使用经典FF预测的τc与使用MB-pol预测的τc的相关图。黑线代表y = x函数，垂直虚线代表振动模式。

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

通过MB-pol和经典FF之间的相关图（图5B）进一步比较了光谱扩散动态。该图显示计算出的值有显著的散布，其中大多数FF得到的τc值小于MB-pol的值。这一结果表明，经典FF倾向于在HDO分子周围的化学环境中产生更快的动态。这种现象可能与使用完全刚性水分子的模型有关，这些模型考虑了旋转和平移自由度，但忽略了更精细的微观结构运动。像4位点和5位点这样的更复杂模型明确地将负电荷与氧中心分离，与3位点模型相比，能够更真实地表示四面体氢键网络。因此，3位点模型中氧原子上的电荷集中可能会人为地使势能表面变平，这一特征可能会影响局部结构和动态的描述，导致性能不如MB-pol。总体而言，所有非极化的经典FF都无法定量预测光谱扩散，表现为低估了去相关时间，且与MB-pol产生的值相比误差超过20%。然而，即使时间尺度不是很准确，结果也表明TIP4P/Ew和OPC模型更准确地捕捉了氢键动态，尽管在后一种情况下，这可能是通过补偿误差实现的，因为OPC没有正确描述频率波动（见前一节）。

最后，还必须讨论经典近似在构象采样中的影响。核量子效应（NQEs），如零点能和隧穿，在液态水动力学中起着不可忽视的作用。（117）特别是，NQEs通常会软化氢键网络，允许比经典模型预测的更大的质子离域。在我们的MD-IFM结果背景下，由于包括了量子核效应，频率波动幅度可能会受到影响，预计会增大。此外，τc值也可能改变，因为量子波动促进了氢键网络的重新配置，降低了氢键断裂和形成的有效障碍。（117）

结论

这项研究表明，MD-IFM方法是评估分子FF的强大工具。特别是，MD-IFM成功地用于评估一系列常见的非极化水模型，以MB-pol作为黄金标准，并提供了对FF表示液态水结构和动态准确性的独特见解。为此评估的FF包括：TIP3P、SPC/E、TIP4P/Ew、TIP5P和OPC。所有这些FF在预测平均振动频率方面都表现良好，表明它们总体上正确描述了水分子经历的平均局部结构。基于被探测水分子的第一溶剂化层中水分子数量的基础分布分析证实，每个力场都一致地再现了由于氢键状态变化而预期的振动溶剂化色变。然而，对频率分布的分析揭示了水氢键网络固有异质性的系统性低估。这种误差在OPC模型中尤为明显。在动态方面，所有经典FF都表现出比MB-pol更快的振动频率去相关动态，这表明经典模型在描述氢键重排方面存在不足。虽然这项研究表明FF之间没有明确的优劣之分，但更复杂的4位点和5位点模型在描述化学环境动态方面往往表现更好。总体而言，我们强调MD-IFM方法是通过非昂贵的振动可观测量计算来评估经典力场的宝贵工具，这些可观测量考虑了化学环境结构和动态。