最近，TIP4P-FQ模型计算了气态、液态和固态冰Ih三种不同相态的偶极矩，在较低的温度条件下计算的介电常数与实验值还是存在一定的偏差。Rick等人表明如果一个力场模型能够很好地预测溶质周围溶剂的结构，它就能够适应各种不同温度的计算。目前没有一个极化力场模型能够真正给出水分子体系中各部分确切的电荷分布，这对于计算偶极矩和介电常数是非常重要的，因为这两者都依赖于分子体系中电子密度的分配方式。Jorgensen等人也表示，在较高的温度和压力条件下的模拟计算中，对结果进一步改善的首要任务就是利用更多的静电相互作用点以及更加精确的极化。ABEEM-7P模型恰好满足了这两点的要求，利用了七点电荷(每个水分子中有三个原子，二个化学键和二个氧原子的孤对电子)计算静电相互作用。而且原子-键电负性均衡方法可以精确地计算水分子体系的电荷，最重要的一点是体系的电荷能够随外势环境的改变而改变。本节中我们继续利用ABEEM-7P模型计算在269~348K温度范围内液态水分子体系的各种性质，如电荷分布、偶极矩、介电常数、汽化热和径向分布函数等性质，详细地分析了这些性质和温度的依赖关系。
图8-7给出了这些电荷的分布情况。ABEEM-7P模型计算的O原子、H原子、O-H化学键和孤对电子的电荷在298K和其它温度下有一定的相似性。在不同的温度下，ABEEM-7P模型计算的O原子和H原子都带有部分正电荷，与其相平衡的O-H化学键和孤对电子带有部分的负电荷。如260K，O原子的平均电荷为0.1130；273K，H原子的平均电荷为0.3815；310K，O-H化学键的平均电荷为-0.1516；348K，O原子孤对电子的平均电荷为-0.2776。另外，从O原子、H原子、O-H化学键和孤对电子的电荷我们也能够大致地判断出参与形成氢键的受束缚的H原子和孤对电子以及没有参与形成氢键的自由的H原子和孤对电子。比如，和气态水分子团簇(H2O)n(n=1~6)的电荷相类似，在260K液态水体系中某一单体水分子的H原子电荷为0.3894和0.3963，显然这两个H原子是都受到束缚的，都参与了氢键的形成；在273K某单体水分子中H原子的电荷分别为0.4685和0.2716，则带有0.4685个正电荷的H原子一定是受到束缚的，参与了氢键的形成，而另一个带有0.2716个正电荷的H原子则是自由的，没有参与氢键的形成。在不同的温度下，如果能够详尽地分析液态水分子体系中受到束缚的或是自由的H原子和孤对电子，对理解水分子体系不同温度下的复杂氢键网络结构有着巨大的帮助作用。对氢键网络的深入了解，无疑将是破解水分子体系各种特殊性质的关键所在。
在260K、273K、298K、310K和348K温度下，ABEEM-7P模型计算的O原子、H原子、O-H化学键和孤对电子电荷主要存在两点差别。(1)不同温度下，相应原子、化学键和孤对电子的平均电荷存在一定的差别，尽管这种差别并不是很大。比如，在260K、273K、298K、310K和348K温度下，O原子平均电荷分别为：0.1130、0.1133、0.1132、0.1131和0.1121，与气态单体水中O原子电荷(0.1125)相差的绝对值分别为：
0.0005、0.0008、0.0007、0.0006和0.0004。其中，348K时O原子上的平均电荷最小，并且与气态单体水中O原子的电荷最为接近。H原子、O-H化学键和孤对电子的结果与O原子类似，260K、273K、298K、310K和348K温度下，液态水中H原子的平均电荷和气态单体水中H原子电荷(0.2897)相差的绝对值为：0.0915、0.0918、0.0916、0.0901和0.0833，液态水中O-H键的平均电荷与单体水中O-H键电荷(-0.1552)相差的绝对值为0.0036、0.0038、0.0036、0.0036和0.0030，孤对电子和单体水(-0.1908)相差的绝对值分别为0.0953、0.0956、0.0954、0.0939和0.0868。这样的结果是合理的，因为随着温度的升高，液态水的体积增大，密度减小，不同水分子之间的长程相互作用减小，则在较高的温度下原子所带电荷的绝对值要比低温时所带的电荷小。(2)从O原子、H原子、O-H化学键和孤对电子在不同温度下电荷的分布情况来看(图8-7)，在较低的温度时，O原子(图8-7a)和H原子(图8-7b)的电荷分布有两个较为明显的峰，并且随着温度的升高，这两个峰逐渐变得不明显。这说明，在较低的温度时，液态水分子体系中明显含有两种类型的氢原子和氧原子，一种是参与形成氢键的受束缚的原子，另一种是没有参与形成氢键的自由的原子，此时单体水分子的排列较为有序。随着温度的升高，单体水的排列逐渐变得杂乱无章，不能够清晰地分辨出参与形成和没有参与形成氢键的原子，在氧原子和氢原子的电荷分布上就反映不出两个明显的峰。同时，由于水分子体系中氢键的存在，在较低的温度时不同水分子之间的长程和短程相互作用都很大，当温度不断升高时，尽管短程相互作用仍然也很大，但是长程相互作用在不断削弱。总之，利用ABEEM-7P模型计算的水分子体系的O原子、H原子、O-H化学键和孤对电子电荷在不同的温度下能够定量地给出不同的数值，尤其是利用了氢键相互作用区的klp，H(Rlp，H)参数清晰地考虑了水分子间近程的相互作用，为进一步从分子水平上了解液态水体系的结构奠定了基础。
ABEEM-7P模型计算的260K、273K、310K和348K温度下液态水体系中平均偶极矩和介电常数列在表8-5中，并且与已有的TIP5P、TIP4P-FQ模型的计算值和实验值进行了比较。从表8-5中我们可以看出，在不同的温度下，ABEEM-7P模型计算的液态水体系中单体水的偶极矩分布范围是不一样的。在260K、273K、298K、310K和348K，单体水的偶极矩的分布范围分别2.39~3.38D、1.98~3.36D、2.30~3.39D、2.08~3.32D和2.18~3.27D，在相应的温度下，最大值和最小值之间分别相差0.99D、1.38D、1.09D、1.24D和1.09D。单体水偶极矩分布范围最大的是在273K，最小的是在260K。那么，这样的偶极矩分布是否和液态水的一些特殊性质有关呢？这一问题目前我们虽然不能给出详细的解答，但是我们可以粗略地分析一下，273K也就是0C，与液态水的最大密度温度4C最接近，此时溶液的密度最大，体积最小，体系中单体水的极化相应就大。在更低的温度下(260K)，溶液中水分子的排列是有规则的，每个水分子周围的环境比较相似，受到其它水分子的作用相差不大，因而偶极矩分布的范围最小。
到目前为止，不同温度下液态水的单体平均偶极矩在实验上还没有测得，只是利用一种归纳模型得到过冰Ih晶格的偶极矩为3.09D。对于液态水分子体系，ABEEM-7P模型计算的单体水平均偶极矩随着温度的升高逐渐减小，如在260K、273K、298K、310K和348K，ABEEM-7P模型计算的平均偶极矩分别为(表8-6)2.83D、2.82D、2.80D、2.78D和2.74D，与TIP4P-FQ模型取得了一致的结果。但是，ABEEM-7P模型计算的平均单体偶极矩普遍比TIP4P-FQ大，如在260K，ABEEM-7P的平均偶极矩为2.832D，TIP4P-FQ为2.805D，两者之间相差0.027D；在310K，ABEEM-7P的计算值为2.762D，TIP4P-FQ为2.606D，相差0.156D。从260K到310K，ABEEM-7P模型的平均偶极矩改变了0.05D，这个数值远远小于TIP4P-FQ模型在该温度范围内的偶极矩改变值(0.2D)。在260K、273K、298K、310K和348K温度下，ABEEM-7P模型得到的偶极矩分布画在图8-8中。随着温度的升高，偶极矩的分布逐渐变窄，在260K、273K、298K、310K和348K，半峰宽分别为0.489、0.470、0.422、0.404和0.350。
ABEEM-7P模型计算的偶极矩分布随着温度升高逐渐变窄，以及平均偶极矩随着温度降低逐渐变大，都说明ABEEM-7P模型在较低的温度给出更加极化的水分子环境，这和已有的结论相一致。
图8-9中给出了ABEEM-7P模型计算的在260~348K温度范围内的介电常数，以及与已有的实验值、其它水分子经验力场模型进行的比较。在260K、273K、298K、310K和348K，ABEEM-7P模型计算的介电常数分别为92、83、76、73和63，在273K、298K、310K和348K，实验测得的介电常数分别为87.9、78.4、74.2和62.6，在相应的实验温度下，我们的计算结果与实验值偏差的绝对值分别为4.9、2.4、1.2和0.4。可见，在较低的温度时，ABEEM-7P结果与实验值相差较大，大概偏差6%。
但是，随着温度的升高差别逐渐减小，在348K我们的计算结果几乎能够模拟实验测量值，偏差只有0.6%。TIP4P-FQ模型在260K、273K、298K和310K计算的介电常数分别为105、97、79和78，虽然在298K它给出了与实验值最为接近的结果，但是从整体上来说都比实验值大，TIP5P模型计算的298K和310K的介电常数为91.8和81.5，也比相应的实验值大。图4.10给出了更加清晰的比较，随着温度的升高，ABEEM-7P的介电常数与实验值符合得非常好，但在较低的温度与实验值还有一定的差别。TIP4P-FQ模型在高温时的介电常数过大，在低温时的介电常数过小，表现出对温度的较大依赖作用。TIP5P模型则在整个温度范围都过高地估计了介电常数。
8.4金属离子水溶液性质
在自然界中，金属离子起着广泛和不同的作用，例如，它们作为反离子，维持浓度梯度，是生物体系酶作用的重要元素。在分子水平上理解简单金属离子的水化是化学和生化方面的基础问题。对于分析水溶液中的金属离子的各种不同的特性，各种实验技术如NMR、X-ray衍射、Raman光谱等，以及计算机模拟如分子动力学模拟和MonteCarlo模拟都是强有利的工具。从头计算的分子动力学（abinitiomoleculardynamics，AIMD）模拟和Car-Parrinello方法[321~324]都是非常受欢迎的，因为它们能较精确地确定分子间力。例如，White等对一个钠离子和53个水分子进行了第一原则的分子动力学(first-principlesmoleculardynamics，FPMD)模拟[323]；Lyubartsev等利用Car-Parrinello方法模拟了一个锂离子和32个水分子[324]；Marx等用Car-Parrinello方法模拟一个Be2+离子和31个水分子，总的模拟时间是1ps。混合的量子力学和分子力学（hybridquantummechanics/molecularmechanics，QM/MM）方法也发展起来模拟各种体系，这种相对精确和有效的方法已经被用于研究金属离子水溶液体系。但是这些方法的耗时性限制了模拟时间和所研究体系的大小，这样的模拟就不能很好地预测一些动力学性质(需要较长时间的模拟)。经典的(classical)计算机模拟可以对较大体系进行较长时间的模拟，但其需要一个非常精确的分子间相互作用势能函数。对于金属离子水溶液的计算机模拟研究，大多数都是利用固定的点电荷模型（相当少的模拟采用浮动电荷模型），而固定电荷模型的模拟不能很好地描述体系静电环境随着时间的变化。本文将用基于ABEEM/MM的浮动电荷模型对金属离子水溶液进行分子动力学(MD)模拟，即ABEEM/MM-MD模拟。
8.4.1自由能计算自由能是热力学中最重要的物理量。自由能可以是亥姆霍兹（Helmholtz）自由能或吉布斯（Gibbs）自由能。Helmholtz自由能要在NVT系综下测得，Gibbs自由能要在NPT系综下测得。这里主要讨论NVT系综计算Helmholtz自由能。根据统计力学，NVT系综下，Helmholtz自由能与配分函数Q（partitionfunction）相关，即：
AkTlnQB（8-4）自由能差（自由能的变化）可以如下计算：
理论上，配分函数可以直接由经典力学的配分函数定义来计算，而且自由能差的计算只取决于起始态的能量差。
其中，H(pN，rN)是系统的哈密顿量，为动能kineticK和势能potentialE的加和。在实际计算中，配分函数转化为系综的平均，而且动能项可以消去，则哈密顿量差最终变为势能的差。实际上，自由能的变化小于2kTB时，上式的结果非常精确，但若自由能差超过此范围，则导致统计的不准确性，即若状态1与状态2的相空间(phasespace)的重叠性高，则可利用这种方法计算，反之，若重叠性很低，则不易得到正确的自由能变化值。为了克服重叠性低的障碍，可利用分段的方法。设状态i介于状态1和状态2之间，则依此，可将状态1到状态2的变化分成许多自由能变化不大的阶段，分别计算各个阶段的自由能差的变化，再将其相加得到两态的总的自由能差。这是计算自由能差的第一种方法——热力学微扰法（thermodynamicperturbation），即是将状态1改变为状态2的势能变化分成许多中间态，每个状态的能量可以写成参数的函数（由0变到1）：
