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

第7章 可极化静电势分子力场(3)

这里,(0)αE是原子α的基态能量,iααjr表示分子i中原子α与分子j中原子β之间的距离,()αβiααjJr是库仑相互作用,()iααjVr表示任何外加的非库仑相互作用。

依据电负性均衡原理,AQ为分子中各原子电负性相等时的电荷,是分子电负性。

通过上述方程可以得到一系列的线性方程组,同时保证整个分子电中性的限制。那么,这种限制可分成两种情况:

(1)整个体系保持电中性,体系所有原子的化学势都相等,分子间有电荷转移。

(2)体系中每个分子保持电中性,每个分子中的所有原子的化学势都相等,分子间没有电荷转移。

为了保持电中性,采用了Lagrangian乘因子方法,对整个体系保持电中性的方程表达为1L,对每个分子保持电中性的方程则为2L。

这里,αm是原子的质量,αM是虚电荷质量,单位为时间平方除以电荷平方(时间平方/电荷平方),为Lagrange乘因子。原子核的运动符合Newton运动方程,电荷随时间变化方程为:

在算法上,两种均衡方法没有区别,电荷所受的作用力为平均电负性和该点的瞬时电负性的差值。如果它的瞬时电负性高于平均电负性,那么该原子电荷就受到使其电荷减少的力,直到二者电负性相等。在实际应用中,他们选择了第二种电荷限制方案,也就是说,不允许分子间的电荷转移。这是比较合理的,因为在较远距离,两个原子之间的隧道效应很小。电荷的质量QM是一个虚质量,所以应该选择得足够小以保证满足在原子核的位置变化后,电子能迅速地调整到原子核势场相适应的位置,在分子力学计算中服从Born-Oppenheimer近似前提,在TIP4P-FQ模型中QM等于6.0105,SPC-FQ模型为6.9105,但是也要保证运动方程的解不是太小,电子与核之间没有热量交换,电荷的自由度保持在0K左右。根据分子动力学知识,这一点可以通过Noséthermostatting方法来保证电荷在较低的温度(6K左右)。

根据Rappé和Goddard的电荷均衡方案(QEQ)[134],分子内的库仑相互作用项()ijJr等于以两个原子为中心的Slater轨道的库仑积分:

n1为主量子数,为指数,iA是归一化因子,对于r0时,()iiJr就是在该力场中,当两个原子距离较近时利用屏蔽静电相互作用势,只有当两原子间距离超过2.5时,ijJ才等于ij1r。最后根据电负性均衡原理,得到一系列的方程,加上电中性限制,就可求解每个原子的电荷,从而利用力场的函数方程计算体系静电相互作用能量。

2004年,Patal等人在以往的浮动电荷力场方法的基础上,建立了CHARMM浮动电荷力场[108~109],该模型基本采用了Rick等人的模型,但是在确定模型参数电负性和硬度时,采用了Banks等人的线性响应模型,拟合从头计算数据得到。Stern等人后来发展的浮动电荷和诱导偶极相结合模型也利用该方法拟合参数,这将在下一节详细论述。

CHARMM浮动电荷力场决定了原子电负性和硬度参数时采用线性响应方法,就是利用类似水分子的虚偶极子作为模型分子的一外加静电场,对该分子进行扰动,将这样的虚偶极子放置在分子的周围,如在氢键位置、垂直于芳香环平面的上下(此时电子云充当氢键受体)和一些其它任意位置。有外场时的静电势能为:

方程(2-34)中的最后一项为外势能贡献项,iΦ表示原子i的外静电势。那么,在没有外场时其静电能量为:

这里0iQ是在无外场时原子i的电荷。ij是原子的硬度,硬度参数可以通过组合方法获得,这里,ijR是原子i和原子j的距离,局域屏蔽库仑势当ijR大于2.5时等于1r,这种相互作用只应用于具有(1,2)、(1,3)和(1,4)关系的原子间作用势。一旦超出1-4关系就通过简单的库仑相互作用势计算。对不同分子间的静电相互作用采用简单的库仑作用势。这种屏蔽作用模型比利用原子中心的Slate轨道的库仑重叠积分方法计算要节约计算时间。

在有无外场条件下达到平衡时的条件分别是:

由上式就可以得到以下两种情况下的两个方程:

用方程(2-38a)减去(2-38b)就会得到试探外场对模型分子电荷分布的响应表达式:于是就得到了原子硬度和由于外场引起的电荷响应之间的关系式。

对于一个误差函数(2-40),表示参数化模型和实验结果的差值,该模型采用了密度泛函理论计算的电荷响应结果作为标准数据。密度泛函方法计算中,利用B3LYP交换相关能函数形式和6-31++G(2d,p)基组,拟合静电势电荷方法计算在外势下分子的电荷响应。

通过最小化以上误差函数,就会得到优化的模型参数。在参数化的过程中,首先从简单的烷烃碳原子和氢原子类型开始,然后到更复杂化学环境下的原子类型。

该方法又设计了另外一种确定硬度参数的方法。我们知道,偶极极化度张量可以从硬度矩阵元导出,其方程为:

原子相对几何中心的坐标,这一方程可以拟合量子化学计算和实验的分子极化度来拟合一系列的原子硬度参数,但是由于该模型采用的是原子中心的浮动电荷,所以无法用来表示各向异性的极化度,如果对芳香环结构的分子加入非共面的点电荷,对氧和氮原子加入孤对电子点电荷,那么就可以用该方法更加简单地计算模型参数,同时也可以在参数化的过程中充分考虑到实验的结果,使得参数更适应实际应用体系。

2.3.3浮动电荷与诱导偶极相结合模型

单独的浮动电荷模型和诱导偶极模型都在一定程度上对极化力场的发展起到了重要的作用。但是目前的每一种单独使用的模型都有各自的不足,尽管浮动电荷模型计算效率高,但是浮动电荷模型存在一定的局限性,如对于一个平面分子,不能体现非平面方向的极化,又如像水溶液中与氧原子形成的分叉氢键,也很难准确定性的描述;而诱导偶极模型又不能体现分子中各原子随外势变化而引起的电荷转移,而且计算量大,该模型的计算量与其偶极位点个数的9倍成正比。在这种情况下,1999年,Stern等人发展了两者相结合的方法,采用与之前的浮动电荷模型相同的线性反应方法获取参数,该模型利用该方法得到令人满意的三体能和丙氨酸二肽、四肽的构象能,因而他们又利用此前的参数拟合方法,对组合方法进行了研究。

总的静电相互作用能的表达式为:

表示原子A的电负性、硬度和浮动点电荷;0A和0BE是原子A所受的外势和原子B所受的外场;Bα和B是原子B的极化度张量和诱导偶极矩;AAAq12Jq是任意原子A关于自身浮动电荷Aq的自身电荷极化能量项;q是原子A的浮动电荷与外势的相互作用能;AAA"A"qJq是原子A和A"的浮动电荷之间的静电相互作用能。

原子B的诱导偶极与外场相互作用势能项;BBB"B"T是原子B和B"的诱导偶极之间的相互作用能。

原子A的浮动电荷与B原子的诱导偶极之间的相互作用能。

定义现在,就可以定义一个ABN3N维的向量q和v和一个ABN3N乘以ABN3N维的矩阵J,AN表示浮动电荷的数目,BN表示偶极点的数目。