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

第8章 可极化静电势分子力场(4)

通过对方程(2-40)进行能量最小化,就会得到每一个变量Aq,B的解。但是如果这样,就等于没有对变量进行边界条件的限制,此时通过最小化得到一系列的线性方程组:

如果对浮动电荷变量Aq进行限制,使每一个分子都保持电中性,可以通过Lagrange乘因子的方法来处理,或者对一约化系列的非限制性变量q"进行转换,并可以找到矩阵C,使得Cq"q,在此时方程的解为:

此方法在计算长程相互作用时,用简单的库仑相互作用表示,对短程采用屏蔽库仑相互作用。当BB"r小于BB"K时,这种屏蔽相互作用采用的是有效距离BB"K而不是真实距离BB"r。

实际上,截断半径大小等于两种原子类型的半径之和,每一种原子类型的截断半径都可以通过拟合得到。

通过在分子周围步设网格点,采用了最小二乘法拟合从头计算静电势的方法来拟合参数。对分子的每一个构象,都进行了一系列的静电微扰v(r):

网格点的位置,KΦ是由于外势的微扰,在网格点k处的静电势变化,如果上式用矩阵表示为ΦRq,AB,qq,那么。参数的拟合就归结为对误差函数2的最小化问题,对所有的微扰求和,假设变量q是非限制的;拟合的一个重要条件就是矩阵J必须是正的有限值,以至于以上方程能找到最小值,而不至于造成“极化塌陷”,由于通过网格点拟合静电势不够稳定,因此可以通过矩阵的奇异值分解法,或者通过误差函数、退火和共轭梯度法来拟合参数。

2002年,Kaminski等人提出了相对于以上模型相对简单的键增加电荷与诱导偶极相结合的静电势模型[105],同时偶极位点不仅包括原子,而且扩展到氧的孤对电子,拟合精密的从头计算结果,建立了第一代完全的蛋白质极化力场,应用于气态的多肽构象研究。首先,他们把关于位点i的诱导偶极矩i能量表达为:

方程的二次项是我们熟悉的诱导偶极矩的自身能量;i是位点i的偶极极化度;一次项系数i表示位点i的内禀电场的负值,也就是说,在没有其它位点电场的情况下,电场依然存在,实质上是对一个孤立分子引进了一个永久偶极矩。那么,在该模型下的极化静电势能的表达式就可以表达为:

键增加电荷,在键相邻的原子上分配一定的电荷,同时在外势变化时候电荷可以在键相连的原子间转移,电荷变化量相等,同时也可以保证分子呈电中性;ijklJ,是位点i,j和k,l键增加电荷之间相互作用的耦合标量元。

是位点i,j键增加电荷与在位点k处偶极相互作用的耦合矢量元。

表示位点i和j的偶极之间的相互作用的耦合二级张量。

通过对方程(2-54)总能量的最小化,从而求得偶极矩,那么就要求:

对以上方程可以通过自洽场方法,或者矩阵对角化方法来求解,文章中提到了如何解决在诱导偶极模型的极化塌陷问题,他们指出:虽然目前没有使用静电屏蔽作用,但是,Lennard-Jones势能项可以防止其发生。偶极矩可以赋以虚质量和动能,利用Lagrangian方法对空间坐标进行积分。

2.3.4极化力场小结

以上我们介绍了三种应用于蛋白质体系的极化力场,它们都毫无疑问地为我们提供了一些更加有效的模型来描述和计算生物蛋白质分子内和其在溶液中的静电极化,解释和计算了通常固定电荷力场所不能描述和计算的现象和物理量。该模型把诱导偶极表现为一个无质量虚电荷位点(Drude位点),通过一个假想弹簧与原子位点相连,连接两个位点的向量由极化位点的极化率决定,该方法在计算方法和计算效率上较诱导点偶极模型有显着提高,但该模型的物理本质不明朗,其模拟结果同诱导偶极模型相当。那么,各种不同的模型既有共同之处,也有各自不同的特点,下面我们就以氮甲基乙酰胺分子为例加以讨论:

1.诱导偶极模型利用分子偶极矩来拟合原子偶极极化度参数,参数拟合相对简单,而且能有效地利用实验极化度数据,该模型能更好地符合实验值,如果采用各向同性极化度,需要12个原子固定电荷参数和12个原子各向同性极化度参数,相对多极模型,它的参数较少,但是该模型无法描述不同方向极化差异;使用各向异性的偶极极化度,就会需要12个原子固定电荷参数和36个原子各向同性极化度参数,这时计算量与其偶极位点个数的9倍成正比,计算变得更加复杂。Ponder等人的诱导多极模型,虽然能够更好地描述分子的极化,尤其在存在氢键和共轭的体系,但是计算量会大大增加,对于像蛋白质这样的大分子体系来说,计算实在存在一定困难。

2.无论是诱导偶极或者是多极模型,他们都没有考虑到在外场的作用下,或者是分子自身构象变化时,分子内部和分子间的电荷转移,实际上这种影响在生物体系中非常重要。

3.Rick的浮动电荷模型,采用Mulliken孤立原子的电负性和硬度,屏蔽库仑势采用原子Slater轨道库仑重叠积分的方法,该方法对限制性体系(刚性,如键长固定)的计算是合理的,在模拟初始时计算一次,之后就采用不变屏蔽势,但是对于非限制性体系来说,计算量很大,很难适用;该方法只考虑了沿着化学键方向的电荷极化,对非共面或分叉氢键等都不能体现有效合理的极化。但总体来说,该模型的提出是第一次把浮动电荷模型应用到力场中,具有里程碑的意义。

4.CHARMM浮动电荷模型,它在模型上采用了Rick模型,不同之处在于,采用了半经验的屏蔽库仑势,比轨道重叠积分方法简单,大大地减少了计算时间;另外,他们还采用了线性响应方法来拟合电负性和硬度参数。该方法也只考虑了沿着化学键方向的电荷极化,对非共面或分叉氢键等都不能体现有效合理的极化。但是,他们提出了利用分子或体系极化度来拟合硬度参数的方案,如果设置非原子中心点电荷,如孤对电子和π电荷,利用各向异性分子极化度就可以解决以上弊端。

5.浮动电荷和诱导偶极相结合的静电势模型,从原理上讲,它既考虑到了外势变化体系内部的电荷转移,又体现了各向异性的静电极化,利用线性响应模型确定参数。为我们提供了一个很好的极化力场静电势模型。但是,该模型的参数确定和静电势的计算都很复杂,对于一个NMA分子而言,需要12个电负性和12个硬度参数,还需要36个各向异性极化度参数,计算量与浮动电荷点的个数的平方成正比,与偶极点的个数的9倍成正比,对于蛋白质体系来说运算量太大。所以说,我们应该平衡考虑计算的精确度和计算量,目前也没有足够的证据证明该模型的计算结果相对其它模型有明显的优势。

6.Kiminski等人的浮动电荷和诱导偶极模型,在Stern等人的模型的基础上,采用了键增加电荷,也就是说,只考虑相临成键原子间的电荷转移,从而简化了模型;另一方面,他们又增加了氧原子的孤对电子区域的诱导偶极位点,可以更好的描述氢键作用。但是目前他们仍旧只是采用了固定的键增加电荷,浮动电荷点也只是放置在原子位置上,没有放置孤对电子和电荷位点,尽管使用了氧原子的孤对电子区域的诱导偶极位点,但是不能忽略电荷的转移极化;另外,只考虑成键原子间的电荷转移也是不太合适的。