JCP：基于第一性原理分子动力学热力学积分的离子溶剂化自由能计算

本文提出了基于第一性原理分子动力学热力学积分（TI）计算离子溶剂化能的方法，该方法具有坚实的理论基础。计算离子溶剂化能的TI过程包括热力学循环，该热力学循环包括空腔形成、离子添加和空腔消除三个步骤。

首先，采用一个外加势场，在排斥势的作用下，在水溶液中创建一个空腔。以空腔半径R_c为反应坐标，改变球型外加限制势的半径，以0.4 Å为增量，逐步创建半径为3 Å的空腔。然后以离子的质量为反应坐标，通过一个质量百分比因子以10%的增量从0%变化到100%，将一个离子逐渐加入到空腔中。对于这个子过程（离子加入体系），在每个AIMD步骤计算了两个DFT体系：一个有阳离子，另一个没有阳离子。体系的总能量和原子受力使用这两个AIMD对应的质量的百分比因子进行加权平均。使用这些加权平均的原子受力来确定AIMD轨迹。由于空腔的存在，这两个体系在DFT SCF中都很容易收敛，解决了非物理的原子构型在第一性原理计算中自洽收敛会非常困难的问题。第三，在含有阳离子的水溶液中，通过逐渐减小R_c（势阱的截断半径）来消除空腔，通过每一步 0.4 Å 的改变量逐步将限制势半径减小到 -1 Å 。这三个步骤相加产生了H₂O（l）和M ⁿ⁺· H₂O（l）之间的自由能差。由于整个过程在平均静电势为零的周期性超胞中进行，因此本文对带电体系进行了特殊的修正。

图1. 使用 PWmat 设置的TI 模型的原理图。(a)热力学过程的分解。Mⁿ⁺（vacuum）、H₂O(l) 和 Mⁿ⁺·H₂O(l) 分别代表真空、液态水和液态带电水中的阳离子。TI 用于创建空腔 (ΔG_1c)、添加离子 (ΔG' ₂) 和消除空腔 (ΔG_1e)。ΔG' ₂ 和ΔG₂之间的区别在于ΔG' ₂ 考虑了离子的自由能。空腔是位于WB和CCWB中心的球形区域。(b) 原子结构示意图。绿色、红色、粉色和灰色球体代表Li离子、O原子、H原子和势阱引起的空穴。空腔尺寸由势阱的截止半径决定。

上述已经提到，TI方法将一个给定的热力学过程分为三个步骤。每个步骤中热力学量的变化通过热力学量相对于反应坐标参数的偏导数积分得到。本文计算的离子溶解的溶剂化自由能，对应的反应过程可表示为：

所以溶剂化能ΔG_s可以通过下式计算：

其中，⟨⟩表示玻尔兹曼因子的统计平均值。在模拟中，它表示 AIMD 轨迹的平均值。另外，F代表总受力。在区间积分中，使用了从-1.0到3.0 Å的10个点。空腔消除是类似的过程。

由于周期性带电体系的电势是发散的，因此在第一性原理计算软件中都加入了均匀的背景电荷来保证体系实际上是电中性的（凝胶模型）；因此有离子和没离子情况下的电势零点不一样，不能直接做对比，需要找一个对齐的势能。类似的原理在计算缺陷能级时也会涉及。因此为了处理带电体系，需要引入静电能量校正。该静电能量校正包括两项：一项是电势校正项，另一项是镜像相互作用项。为了处理带电体系，本研究使用电位对齐的方法来消除带有中心带电阳离子的水盒子（CCWB）和不带阳离子的中性水盒子（WB）之间的电子电位差。此外，还通过合理的估计考虑了镜像相互作用项。同时使用 DFT-D2 近似考虑 vdW（范德华）相互作用。

首先，为了确定该电位偏移，使用表面溶剂计算确定真空能级设置为零时的实际平均电位。在表面溶剂计算中，真空能级可以移至零，因此系统电势不是任意的。如果正确电位表示为 V（中性），则校正项为

这里，求和是针对氧原子的。另外，“a”是采样区域的半径，r是距采样区域中心的氧原子距离。对于表面计算，采样中心放置在表面区域的中心。对于 CCWB，采样中心位于 (0, 0, 0)。V_r(O)表示PWmat代码提供的每个O原子的势值。

在周期性边界条件下，镜像相互作用项对于静电能量也至关重要。镜像交互项是由离子及其镜像电荷在周期性结构中引起的。根据Suo等人的镜像电荷相互作用校正研究，原子单位的镜相互作用项（E_q）为

其中E(Mⁿ⁺(gas))是气相中Mⁿ⁺的总能，T是温度（设置为300 K）。S是Mⁿ⁺的熵，对于压力为P的理想气体，其熵可以计算为

其中，−RT ln 24.46 表示离子从 1 原子气态转移到 1 mol L⁻¹ 溶液期间标准态差。求得Li⁺、Na⁺、K⁺、Be²⁺、Mg²⁺ 和 Ca²⁺ 的 G[Mⁿ⁺(gas)]，进而求得离子溶剂化能及其标准差，最终计算出，Li⁺、Na⁺、K⁺、Be²⁺、Mg²⁺ 和 Ca²⁺的离子溶剂化能分别为 -4.800、-3.974、-2.854、-25.705、-19.680 和 -16.650 eV，标准差为 0.023、0.023、0.035、0.052、0.041 和 0.051 eV。溶剂化能计算结果与实验数据吻合较好。

图 2 比较了本文提出的模型和五种主流理论模型获得的离子溶剂化能的模拟结果。目前存在五种主流理论模型：（1）水力场MD，主要采用经典力场MD的显式溶剂模型，或结合准化学理论的DFT相互作用势与MD；

（2）杂化溶剂模型，内部有显式离子水簇，外部有隐式水模型（偶极连续溶剂）；

（3）基于玻恩模型或平均球近似（MSE）的近似模型；

（4）隐式溶剂模型，其中利用泊松模型构建偶极连续介质溶剂模型；

（5）基于准化学理论的簇-连续介质模型。一般来说，‍水力场MD方法得到的计算结果相对实验结果偏正。相比之下，近似模型和隐式溶剂模型由于采用刚性拟合方法，产生的结果相对实验结果偏负。此外，虽然混合溶剂模型是一种良好且有效的方法，但它依赖于显式的溶剂簇大小。所提出的模型与之前提出的混合溶剂模型具有相似的性质。然而，当前提出的基于 TI 的模型不依赖于团簇大小，因此更可靠。

此外，还计算了阴离子（F⁻）的溶剂化能以验证所提出模型的普适性，计算结果（-4.480 eV）与实验数据相差约0.5 eV，与阳离子结果相似。

图2. 使用Exp（实验），水力场MD，混合溶剂模型，近似模型，隐式溶剂模型，和簇-连续介质模型等不同模型获得的 Li⁺、Na⁺、K⁺、Be²⁺、Mg²⁺ 和 Ca²⁺的溶剂化能

这里，M(s) 代表块体金属的吉布斯自由能。因此，中性块体金属的溶剂化电压 (ΔV) 可确定为

其中 eU = -4.44 eV 是标准氢电极电势，因为 ΔV 是使用 U 作为参考来定义的。为了获得溶剂化电压，固体金属自由能计算如下：

其中，E(M)表示块体的基态能量，ZPE表示零点能量，TS表示熵校正，N表示晶胞中的原子数，δμ (M)表示单个金属元素的密度泛函形成能系统误差的修正项。

图3显示了Li、Na、K、Be、Mg和Ca金属的电极电位及其标准差，金属的电极电位分别为-2.89、-2.65、-2.75、-2.18、-2.45和-3.13 V，标准差分别为0.023、0.023、0.035、0.026、0.020 和 0.025 V。相应的实验值为-3.04、-2.74、-2.91、-1.97、-2.35和-2.87 V。在溶剂化电压结果中，Ca的计算结果与实验结果的差别最大（0.26 V）  。

图3. 第一性原理MD热力学积分方法计算得到的溶剂化电压与实验的对比

综上所述，本文提出了一种使用分子动力学热力学积分的方法计算离子溶剂化能。计算的Li⁺、Na⁺、K⁺、Be²⁺、Mg²⁺ 和 Ca²⁺ 的溶剂化能与实验结果接近，并且在实验误差范围之内。还获得了这些离子的溶剂化电压。此外，与之前提出的研究离子溶解问题的三种方法中的水团簇模型相比，该模型不受团簇大小依赖以及缺乏严格的准谐近似所带来的不确定性的影响。因此，本文提出的方法为计算离子溶剂化能提供了可靠且稳健的方法。以上所有的DFT计算模拟均来自于软件PWmat。