书城自然科学分子模拟力场方法与应用
8392100000015

第15章 水分子力场(1)

分子力场的各种概念都能够利用液态水分子的各种经验模型来检验。尽管水分子很小,但是作为所研究的各种力场模型来说,水分子是很好的范例,它的许多可观察的性质,比如结构、动力学性质和热力学性质都可以简单地用计算机模拟,并且和已知的实验比较,从而确定经验力场模型的准确性。液态水分子体系也是一种富有挑战性的模拟体系,到目前为止曾提出过大量的经验模型。用来模拟液态水的大多数力场模型都是利用有效的对势能函数,没有清晰地考虑三体作用以及极化作用。在简单的分子间点作用模型中,每个水分子都是刚体模型,分子间相互作用利用库仑作用和vanderWaals作用描述,非刚体的水分子模型允许分子构象的改变,目前逐渐发展的模型则清晰地考虑了多体相互作用和极化作用。

6.1简单的水分子模型

常见的简单水分子模型[207~215]从早期的ST2和ST4模型逐渐发展到简单点电荷力场(SimplePointCharge,SPC),以及可转移的分子间势能函数(TransferableIntermolecularPotential,TIPnP),这些模型都是利用刚体水分子结构,即固定的键长和键角。水分子间的势能面利用vanderWaals相互作用,以及三点、四点或五点间的静电相互作用来描述。通常使用的三点、四点和五点的简单水分子力场模型画在图6-1中,相应的力场参数列在表6-1中。

Berendsen等[212]提出的SPC模型、Jorgensen等[215-216]提出的TIP3P模型以及改进的SPC/E模型都是利用三点计算静电相互作用,其中氢原子上带有部分正电荷,与其相平衡的氧原子上带有部分的负电荷。不同水分子间的vdw相互作用通过Lennard-Jones势函数计算,但是只考虑了不同水分子中氧原子之间的相互作用,没有考虑氢原子之间以及氧原子和氢原子之间的vdw相互作用。TIP3P、SPC和SPC/E模型在水分子的几何构型、部分电荷和vdw参数上都略有不同。最早的四点模型是由Bernal和Fowler在1933年提出的[207],简称BF模型。虽然该模型目前已经很少使用,但由于其奠定了水分子模型的发展基础,因而还很重要。最常使用的四点水分子力场是由Jorgensen提出的TIP4P模型,该模型把负电荷从氧原子移到沿氢原子方向的角HOH的平分线上(图6-1),氧原子上没有电荷分布。五点模型最早是由Stillinger和Rahaman在1974年提出的ST2模型,其主要特点是考虑了水分子的四面体结构,氢原子上带有部分正电荷,负电荷放在氧原子的两个孤对电子上,而且利用开关函数调整静电相互作用。最近,Jorgensen等人又提出了TIP5P模型[218],虽然它的几何结构类似于ST2模型(图6-1),但是氢原子和氧原子孤对电子上的部分电荷以及vdw参数都存在一定的差异(详见表6-1)。下面简单地介绍这些水分子力场模型的势能表达形式。

早期水分子ST2和ST4模型描述的势能面包括水分子间的vanderWaals相互作用和静电相互作用。Vdw作用只考虑氧原子,静电相互作用考虑了每个水分子上两个氢原子和氧原子两个孤对电子上带有的部分电荷,并且只计算不同水分子之间的作用。ST2和ST4模型的特点是在静电相互作用中乘以一个开关函数S(r):

其中1和2表示两个不同的水分子,RL=2.0160,RU=3.1287。需要强调的是,ST2模型与真实水分子中氢键的强度和取向是一致的。ST4模型和ST2模型的区别在于:

(1)氢原子的部分电荷qH从0.2357增加到0.2457;(2)孤对电子-氧原子-孤对电子(lp-O-lp)之间的夹角从理想的四面体角度109.47°增加到134.47°。

简单点电荷力场(SPC)将水分子视为刚体分子,其键长和键角的值固定。类似于ST模型,SPC力场也将二个水分子之间的相互作用分为vdw和库仑作用,而且仅O-O原子间有vdw参数。水分子中每个原子均带有部分电荷,不但计算了分子间库仑静电相互作用,还计算了分子内静电作用。SPC力场的形式为:

利用SPC力场可以计算水分子体系的各种热力学性质,但是利用该力场计算的扩散系数与实验值有较大的偏差(误差在200%~300%)。SPC/E力场(extendedSPC)是对SPC的改进,增加了氧原子和氢原子的部分电荷。SPC模型中氧原子电荷qO为-0.82,而SPC/E模型中氧原子电荷为-0.8472,部分电荷的增加使SPC/E模型给出较大的偶极矩、介电常数以及扩散系数等性质。2002年Glttli等人又对SPC模型进行三点改进,提出了SPC/A和SPC/L模型,这三点改进分别是:(1)vdw相互作用不仅仅作用在氧原子上,同时也作用在氢原子上;(2)增加了氧原子和氢原子的部分电荷;(3)改变了水分子的几何构型(HOH键角和OH键长)。改进的SPC/A和SPC/L模型的水分子势能函数为:

其中C6(i,j)和C12(i,j)分别为原子i和原子j间Lennard-Jones势能函数的吸引和排斥作用项系数。SPC/A模型与SPC模型的区别在于SPC/A模型在氢原子间存在一个吸引vdw作用,氧原子间的Lennard-Jones势的排斥系数较小,并且在原子上有较小的部分电荷;SPC/L模型与SPC模型的区别在于SPC/L模型的键长较大,键角较小,原子上的部分电荷较小,同时在氢原子间也存在一个吸引vdw作用。

可转移的分子间势能函数(TIPnP)也将水分子视为刚性分子,即键长和键角的值恒定。TIPnP力场因其形式简单,常与其它力场合并应用于生化系统的计算。TIPnP力场的基本形式为:

其中,a和b为不同水分子上的原子。从图6-1可以看到,TIP3P模型中单体水不具有四面体结构,这导致了TIP3P力场模拟的液态水结构与实验值不一致,因而Jorgensen将TIP3P模型改进为TIP4P模型,大大改善了对液态水结构的模拟。TIP5P模型的建立由于充分考虑了水分子中氧原子的孤对电子,因此更容易解释水分子体系中氢键的形成,从而给出了更加准确的最大密度温度。

另外一种简单的水分子力场模型是中心力场模型(CentralForce,CF)[116]。中心力场模型没有利用vdw和库仑静电相互作用描述水分子体系的势能面,而是分别给出了O-O、O-H和H-H原子之间的势能函数表达形式,计算了其它水分子力场不能给出的水分子内部的振动模式(intramolecularmodesofvibration)和液态水的自离解能力(thecapacityofself-dissociationintheliquidwater)。其表达形式为:

其中r为分子内或分子间原子对的距离。利用中心力场模型计算液态水的结构能够清晰地给出水分子内氢原子-氢原子之间的距离,也能够给出水分子形成氢键时氧原子-氧原子之间的距离,与X-ray衍射实验所观测的结果相吻合。

6.2极化电荷模型

简单的水分子模型虽然能够给出大范围内纯液态水的各种性质,但是这样的模型都是利用固定的电荷参数,不能描述原子(或分子)周围电场改变而产生的静电相互作用的改变。在固定电荷力场中,电荷的分布必须反映出一定状态下的平均场电荷,因而不能转移应用到其它热力学状态和非均相体系。由于非均相体系中具有较强的电场梯度,如离子或溶质-溶液表面的存在,所以能够详细地描述极化效应或多体作用的力场模型可以给出更加正确的非均相体系的各种性质。比如在这样的模型中孤立水分子的偶极矩更接近于气态水的偶极矩而不是“有效”的液态水偶极矩。相对于固定电荷力场模型,极化力场(PolarizationForceField)就是一个“有效”的势能函数,它能够处理电子密度对周围环境的响应。对于分子来说,这种“响应”主要包括两种类型:一个原子周围的局域电子密度的改变,以及电子密度在不同原子之间或原子与化学键间的转移。处理前一种效应的模型可以称之为偶极极化模型(DipolePolarizableModels)[220~250],处理后一种效应的模型称为浮动电荷模型(FluctuatingChargeModels)[250~266],另外还有一种是结合了偶极极化和浮动电荷的水分子力场模型[155]。

偶极极化模型是在20世纪70年代后期建立的处理多体效应的模型。通常认为在经验力场中加入偶极相互作用对于改善计算机模拟溶液体系是必要的,尤其是对于多相化学体系[242,243,249]。偶极极化模型的参数是利用气态单体水、水分子团簇以及液态水的性质拟合得到的,一般认为,满足这两种热力学状态(或密度)就应该能够描述其它的状态。极化模拟方法基于诱导偶极[223,224]来计算局域电子密度的改变,水分子体系是第一个应用极化模型来描述的系统[225,233,237,239,240,242]。

Stillinger和David在1978年首先提出了偶极极化模型[220],并且计算了水分子体系以及水分子离解产物(H,OH)的优化几何结构和偶极矩。该模型是对中心力场模型的改进,一方面保持了中心力场模型允许分子振动和离解的性质,另一方面包括了中心力场模型没有的光学振动上的电子极化。Stillinger和David的偶极极化模型假定水分子体系是氢和氧粒子通过一定的力聚集在一起的集合,液态水中水分子的离解产生了离子片段,形成溶剂化物:

HOHOH2(6-6)假定粒子氢裸露出质子,而粒子氧带有两个单位的负电荷,这样相互作用能就可以写为两部分之和:

IIIUUU(6-7)其中,第一项包含了体系中所有粒子对间的相互作用能之和:

第二项是一个非加和势能项,它的表达形式是通过极化粒子的经典静电学得到的。

由于氢离子仅有一个质子,氢氢之间的相互作用可以写为:

其中e代表质子电荷,类似的UOH和UOO满足一定的渐近限制(即r):

和氢离子不同,离子氧的电子壳层有一定的空间伸展方向,所以在较小的r处UOH和UOO将偏离正常的Coulombic相互作用形式。UOH应该表现出一定的共价键性质,而UOO体现出电子云的重叠排斥。这样在氧离子上加上一个极化率参数,为了计算非离解水分子的气态性质[251],=1.4443。在经典的静电学中,原子i诱导产生的偶极矩

i是由极化率和电场Ei决定的:

对于一组给定的粒子,方程(6-11)和(6-12)是线性相关的,并且唯一地决定了电场Ei和诱导偶极矩

i。根据这种经典静电学得到的极化能量可以简单地写为:

由于氧原子核周围电子云的空间伸展,方程(6-11)、(6-12)和(6-14)必须进行一定的修改。Stillinger和David在模型中加入空间离域响应(spatialdelocalization),但是保持了线性极化响应的特点(linearpolarizationresponse),首先利用修改的矢量场Gi代替了Ei:

其中,K(r)只是一个数量校正方程,实质上它只在氧原子半径周围不为零。类比于方程(6-11),诱导偶极矩可以表示为:

iiμG(6-16)另外的校正是加入了极化能量,参照(6-14)式: