21k()(1)k（8-8）其中，k表示任意一个中间态的势能项，2k表示状态2的势能项，1k表示状态1的势能项。当自由能由ii的取样方法为向前取样法（forwardsampling），当自由能由ii的取样方法为向后取样法（backwardsampling），当同时计算ii与ii的取样方法为双向取样法（double-wisesampling）。计算自由能差的另两种方法为热力学积分法（thermodynamicintegration）和缓慢生长法（theslowgrowthmethod）。热力学积分法即是自由能的变化可以用积分式表示：
缓慢生长法即是将值“缓慢”改变，每一步的()iH都非常接近()iH，理论上，无论采用哪一种方法，所得的自由能变化均应相同，在实际计算中可以根据需要而选择某一种方法。
本文中，钠离子和镁离子的水合自由能的计算是利用自由能微扰和分子动力学模拟。计算经过两个阶段：生长阶段（growth）和带电阶段（charging）[223]。首先，离子电荷设为0，改变离子的vanderWaals参数，根据：
其中，final和final为钠离子或镁离子的终态的vanderWaals参数。生长阶段的最后结构，即1的最后结构为带电阶段的最初结构。带电阶段的模拟计算根据：
finalq()（8-12）在式（8-12）中，finalq为钠离子或镁离子的终态的电荷（+1或+2）。锂离子和钾离子的水合自由能计算以钠离子的带电阶段的最后结构为起点，而铍离子和钙离子的水合自由能计算以镁离子的带电阶段的最后结构为起点。vanderWaals参数根据式（8-8）线性地由钠离子变化到锂离子和钾离子，或由镁离子变化到铍离子和钙离子。
在这些模拟中，对于每个值，50ps的平衡期和150ps的产物期。最后Helmholtz溶剂化自由能由得到的轨迹计算得到，利用自由能微扰方法：
AkTEkTBBlnexp/（8-13）式中的E和E分别为体系的初态（利用初态的力场参数）的势能和微扰态（利用微扰态的力场参数）的势能。
8.4.2不同温度下锂离子水溶液的性质各种实验方法可以测得锂离子水溶液的结构性质，例如，测得锂离子水合数的变动范围为4~6。从头计算确定锂离子是一种四配位的结构，但从头计算只能研究小的团簇，对于溶液性质的研究无能为力。因此，各种计算机模拟方法发展起来，用于模拟凝聚态的特点，但计算机模拟的首要问题是确定一个好的势能模型。基于ABEEM/MM的势能模型计算锂离子水分子团簇的各项性质取得了较好的结果，在这一部分，我们要在锂离子水溶液的分子动力学模拟中进一步测试这个势能模型。
ABEEM/MM-MD模拟计算的248K、268K、298K、328K、348K和368K温度下的锂-氧径向分布gLiO(r)展示在图8-10中，从中可以看出，锂离子水溶液的结构受温度的影响情况。
图8-10不同温度下的Li-O的径向分布从图中我们可以观察到，随着温度的升高，gLiO(r)的第一峰和第二峰都有一定程度的降低。在高温时，gLiO(r)的第一峰有所降低，但是降低的程度并不是很大，这表明即使在较高的温度时，锂离子与其第一近邻水分子之间仍存在较强的相互作用。总体上看，径向分布函数随温度（248~368K）的变化并不很明显，这与其它经验力场的模拟结果以及实验得到的结论是一致的。积分离子氧的径向分布函数gLiO(r)到第一个最小值rmin1，根据min1第一水合层的水分子数ncoor，这些值列于表8-7中。随温度的变化，锂离子的水合结构主要以4配位和5配位的结构为主。
表8-7锂离子水溶液在不同温度下的各项性质T(K)ncoorDH2O(×10-5cm2/s)DLi+(×10-5cm2/s)Residencetime(ps)MSDMSDVAC2484.720.670.250.2598.52684.671.140.500.4670.62984.541.921.131.1647.03284.474.231.651.7129.03484.615.312.592.5021.53684.636.343.083.0919.0转移(transport)是指分子从一个区域流向另一个区域的性质，比如溶液处于一种溶质分布不平衡的状态，这样溶质就要发生扩散直至浓度达到均一的状态。如果体系中有一定的热力学梯度，则能量就要流动直到温度处处相等，同样，动量梯度会产生黏度(viscosity)。转移的存在意味着体系并没有达到平衡。作为一种动力学性质，扩散系数（diffusioncoefficient）无论在实验和模拟过程中都可以被非常准确地获得，因而扩散系数是评价经验力场模型的一个非常重要和基本的参数，此外，扩散系数不但能够反映出分子间势能函数的短程相互作用，还能够表现出长程相互作。扩散系数是和均方位移(meansquaredistance)rt)r0)2相关的，Einstein方程给出均方位移等于2Dt。在三维体系中，均方位移表示为：
上式中，只有当t才是严格成立的，()irt是指分子i的质量中心在时刻t的位置向量，平均值是指要计算体系中所有的分子。除了利用均方位移可以计算扩散系数，速度自相关函数(velocityautocorrelationfunction)C(t)v也可以用来计算扩散系数，根据：
其中，kB是Boltzmann常数，T是绝对温度，m是粒子的质量。因为分子动力学模拟可以提供特定时间的数据资料，这就使得某一时刻的性质与其前面t时刻或后面t时刻的这一性质或不同性质有一定的相关性。相关方程可以表示为：
C(t)x(t)y(t)xy（8-15）当x和y表示不同的性质时，相关函数为交叉相关函数(cross-correlationfunction)；当x和y表示同一性质时，此函数为自相关函数(autocorrelationfunction)。自相关函数表示体系保持它以前“记忆”的程度（换句话说，体系多长时间“失去”它的记忆）。
速度自相关函数即为一个简单的例子，它表示t时刻的速度与0时刻的速度的相关性。
当时间长时，粒子的速度与其初速度已无关联性，故相关函数值趋近于零，反之，时间短时粒子速度与其初速度的关联性强，相关函数值较大。因为粒子的速度为向量，故相关函数值为负值，表示与初速度的方向相反。相关函数为统计力学中最重要的一种函数，可以由不同形式的相关函数计算各种与时间有关的物理量的平均值。
速度自相关函数除以(0)(0)iivv即得到归一化的速度自相关函数C(t)v：
其中，()ivt和(0)iv为t时刻和0时刻的速度。归一化的速度自相关函数在0时刻的值为1，随着时间t的增加逐渐趋近于0。失去相关性的时间经常被称为相关时间(correlationtimeorrelaxationtime)，在分子动力学模拟中，模拟时间一定要大大长于相关时间以减少计算的不确定性。
在本文中，水的扩散系数利用均方位移来计算得到，锂离子的扩散系数利用均方位移和速度自相关函数来计算得到。计算的不同温度下的水的扩散系数和离子的扩散系数列于表8-7中。从表中数据可以看出，无论是水的扩散系数还是离子的扩散系数都是随温度升高而增大，这是与实验测量相一致的，而且利用均方位移和速度自相关函数计算得到的锂离子的扩散系数具有很好的一致性。为了更好地对比我们计算的结果和实验及其它经验模型的模拟结果，图8-11和图8-12中给出了我们模拟的不同温度下的水的扩散系数和离子的扩散系数，以及实验数据和其它经验势能模型的模拟结果。从图中可以看出，我们计算的水的扩散系数和锂离子的扩散系数与实验数据符合得比较好，而略大于1.7M的LiCl水溶液的模拟结果，主要是由于利用不同的势能模型。我们的结果与实验数据的吻合证明了ABEEM/MM-MD模拟计算的锂离子水溶液的动力学性质的准确可靠性。
是Heaviside单位函数，如果一个水分子i在半径为r的区域范围内，其值为1，否则为0，Nr是在这个区域内的水分子的个数，并且我们允许水分子短暂离开区域r的最长时间不超过2ps。如果水分子在离子周围具有较长的停留时间，2ps的离开时间并不能影响统计计算结果。这样，根据式（8-17）就可以得到第一水合层水分子在离子周围的停留时间（如果选第一水合层的半径为r，那么，Nr即为第一水合层水分子的个数）。
而且，我们可以根据下面的近似式来计算大于10ps的第一水合层水分子的停留时间这些计算结果列于表8-7中，我们可以看出：随温度的升高，第一水合层水分子的停留时间呈现降低的趋势，这是由于温度的升高，水分子的运动加快，不能长时间束缚在锂离子周围。我们将不同温度下的第一水合层水分子的停留时间拟合成Aexp(E/kT)的函数形式，其中，E为第一水合层水分子和外层水分子发生一次交换所需要的能量。基于ABEEM/MM的势能模型的模拟给出E为10.91±0.12KJ/mol。
Egorov等利用指数势能函数、Dang和Heinzinger的Lennard-Jones势能函数进行分子动力学模拟，分别给出内外层水分子发生一次交换所需要的能量E为12.9±0.2KJ/mol、11.5±0.2KJ/mol和11.7±0.2KJ/mol。不同温度下的第一水合层水分子的停留时间以及拟合函数曲线（我们的计算结果和Egorov实验组的模拟结果）展示于图8-13中。
图8-13不同温度下的锂离子第一水合层水分子的停留时间以及拟合函数曲线锂离子周围水分子的动力学性质明显不同于外层水体的动力学性质，因此，我们研究了水分子质心的跳跃的长度(jumplength)和角度(angle)，以及离子的这些动力学性质，列于表8-8中。jumplength是指在1fs时间内粒子移动的平均距离，jumpangle是指在两个连续的时间步长内粒子转移的平均角度。从表中数据可以看出，在我们研究的温度范围内，温度对跳跃的长度和角度有很小的影响；但我们可以明显地看到这样的顺序：锂离子的跳长和跳角明显地大于第一水合层的水分子，更大于外层的水分子。虽然跳跃的长度和角度受温度影响很小，但锂离子周围的水分子的动力学性质明显不同于外层水分子的动力学性质，其证明了锂离子对第一水合层水分子的动力学性质有较大的影响。对比以前的报道，我们统计计算的跳长有些小，而跳角有些大，这是由于运用不同的势能模型以及统计的不确定性的原因。
表8-8随温度变化的跳长和跳角Temperature(K)248268298328348368(Li+)(degree)5.334.995.265.374.965.78(H2O-inner)(degree)2.061.951.762.131.952.27(H2O-outer)(degree)1.131.161.141.131.121.17distance(Li+)(×10-38.148.877.849.069.289.33distance(H2O-inner)(×10-35.085.745.795.635.945.94distance(H2O-outer)(×10-34.855.465.205.446.065.338.4.3室温下金属离子水溶液的性质
金属离子水溶液的结构性质主要包括离子-水分子的径向分布函数，角度分布函数（包括离子周围的水分子的取向角分布函数和氧-离子-氧角度分布函数），以及在不同金属离子水溶液中离子对周围水分子结构的影响（即ABEEM-7P水分子结构的变化）。下面我们对这三个方面分别进行讨论。
围绕一个离子的总的结构特点可以用离子-氧径向分布函数gIO(r)和离子-氢的径向分布函数gIH(r)来描述。计算得到的单价碱金属离子（Li+、Na+和K+），二价碱土金属离子（Be2+、Mg2+和Ca2+），以及过渡金属离子（Fe2+和Fe3+）-水的径向分布函数分别展示在图8-14和图8-15中，以及它们的连续积分数。从离子-水的径向分布函数图可以看出，峰位的趋向是与离子半径密切相关的。
如图8-14所示，对于锂离子-氧的径向分布函数，第一个最大峰位于1.95埃，第一个水合层很好地分离于第二个水合层，而且两层之间有很深的最小值区域，其证实了相对稳定的第一个水合层结构。对于钠离子-氧的径向分布函数和钾离子-氧的径向分布函数，它们的最大峰分别位于2.31埃和2.77埃，第一个水合层并没有很好地分离于第二个水合层，两层之间没有很深的最小值区域，说明在模拟过程中两层之间有很明显的水交换反应发生。对于二价的碱土金属离子，由于离子和水分子之间存在较强的静电吸引力，模拟过程中没有发现第一水合层水分子和第二水合层水分子之间的交换现象，这是与单价的碱金属离子不同的。在Be2+、Mg2+和Ca2+离子-氧的径向分布函数中，第一个峰很好地分离于第二个峰，明显地定义了第一水合层，那些尖而高的第一个峰分别位于1.67埃、2.00埃和2.38埃处。对于Be2+离子，我们模拟的结果与X-ray衍射实验测得的Be2+-O距离（1.67埃）、以及Car-Parrinello模拟[125]得到的Be2+-O距离（1.66±0.06埃）有很好的一致性；对于Mg2+离子和Ca2+离子，我们得到的径向分布函数的峰位也与实验数据有很好的吻合，如Neilson和Enderby给出Mg2+-O距离和Ca2+-O距离的实验值分别为2.06埃和2.46埃，而Hewish等进行中子衍射实验获得Ca2+-O峰位于2.4埃处。
