背景导读

离子溶解是电化学中一个重要的过程。电化学反应中许多重要的参数，例如电化学还原电位、无限稀释活度系数、亨利定律溶解常数和离子溶解度等，都与离子的溶剂化能有关。然而，由于测量技术和数据处理的困难，离子溶剂化能的实验测量经常存在不确定性，这也导致许多小分子以及带电粒子的实验数据的缺失。因此，发展一种精确的理论计算方法来确定液相中的离子和小分子的溶剂化能对于电化学的发展至关重要。

文章简介

计算自由能在理论上的一种合理的方法是使用热力学积分（TI）。作为一种计算“炼金术”方法，TI需要在初态和末态之间插入一系列连续变化构型，而这些构型不一定具有实际的物理意义。例如在计算离子的溶剂化自由能的过程中，初态是没有离子的状态，末态是离子存在的状态，那么在这之间需要插入一系列“离子并不完全存在”的状态。对于经典力场而言，由于离子与溶剂的相互作用（主要指静电相互作用和LJ相互作用）是显式存在的，因此可以通过调节这些相互作用大小来模拟离子的存在状态。然而，想在第一性原理计算中引入类似的过程，不仅在流程上更加麻烦，还很容易造成由非物理的原子结构带来的自洽收敛的困难。因此，计算“炼金术”方法尚未广泛应用于第一性原理计算。近日，中国科学院半导体所汪林望教授团队和北京航空航天大学张千帆教授团队提出了一种使用基于第一原理计算的TI方法来计算离子溶剂化能的通用方法, 该方法通过构造合适的热力学循环并加入静电修正的方式，巧妙的避免了引入非物理的原子构型而造成的收敛问题。该方法也适用于其他小分子。相关研究成果发表在了国际著名期刊《The Journal of Chemical Physics》 上，论文第一作者为北京航空航天大学材料科学与工程学院在读博士生蔺超，本文计算部分采用PWmat完成。与TI相关的分子动力学模拟则通过修改后的特殊ASE接口实现。

主要内容

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

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

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

01TI方法计算溶剂化能的总流程

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

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

02空腔的产生、消除和离子添加的溶剂化自由能变化

①空腔产生的自由能变化 

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

②离子添加的溶剂化自由能 

      其中，α1和α2分别设置为0和1，并且〈〉表示使用玻尔兹曼因子的统计平均值，可以用AIMD轨迹平均值代替。

③静电能量校正 

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

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

理想情况下，V（阳离子，电荷）应取在远离阳离子的位置，V（中性）应取在表面计算的中心处， 然而，与带电缺陷的晶体计算不同，水溶液有很大的波动。因此，为了获得平均电势，对一段时间内和特定区域上的氧原子电势进行平均。具体来说，

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

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

03溶剂化能和溶剂化电压

①溶剂化能 

由以上热力学积分的三个步骤的结果，得到H2O（l）和Mn+ ·H2O（l）之间的自由能差 ，即G(Mn+⋅H2O(l)) − G(H2O(l))。再结合文中第二个等式，需要计算G(Mn+(gas))来获得溶剂化能。G(Mn+(gas)) 使用以下方程确定：

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

其中，−RT ln 24.46 表示离子从 1 原子气态转移到 1 mol L−1 溶液期间标准态差。求得Li+、Na+、K+、Be2+、Mg2+ 和 Ca2+ 的 G[Mn+(gas)]，进而求得离子溶剂化能及其标准差，最终计算出，Li+、Na+、K+、Be2+、Mg2+ 和 Ca2+的离子溶剂化能分别为 -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+、Be2+、Mg2+ 和 Ca2+的溶剂化能

②溶剂化电压

结果还用于计算溶剂化电压，即电解金属所需的电压。溶剂化电压的测定使用

这里，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+、Be2+、Mg2+ 和 Ca2+ 的溶剂化能与实验结果接近，并且在实验误差范围之内。还获得了这些离子的溶剂化电压。此外，与之前提出的研究离子溶解问题的三种方法中的水团簇模型相比，该模型不受团簇大小依赖以及缺乏严格的准谐近似所带来的不确定性的影响。因此，本文提出的方法为计算离子溶剂化能提供了可靠且稳健的方法。以上所有的DFT计算模拟均来自于软件PWmat。