Evdw是表示围绕原子核的电子云之间的相互作用,因而在上面的计算中假设原子是球形的。这种假设不满足两种情况。首先是当相互作用的原子中有氢原子,氢原子是单电子原子,该电子总是与邻近原子有一定的键连,所以氢原子核周围的电子分布不是球形的,而是偏向邻近的原子的。解决这种各向异性的一个方法就是在计算Evdw时,把位置移向键的方向。在MM2[27]和MM3方法[28]中采用的协调因子为0.92,也就是代入Evdw的距离是0.92倍的X-H键距离。氢原子周围的电子密度主要取决于X原子的性质,与电负性大的原子(如氧原子和氮原子)相连的氢原子的vdw半径较小,而与电负性小的原子(如碳原子)相连的氢原子的vdw半径较大。大部分力场中都包含了许多不同类型的氢原子,如与碳相连的氢,与氧相连的氢等。其次是含有孤对电子的原子,如氧原子和氮原子。孤对电子要比键上的电子弥散大,表现为原子在孤对电子方向要“大”一些。一些力场在孤对电子的方向放上虚原子来模拟其行为,类似于其他原子,孤对电子有自己的vdw参数,只是这些参数值要远远小于氢原子,这样就可以使氧原子或氮原子在孤对电子方向略有“膨胀”。
氢键的相互作用应该给予特殊的考虑,氢键是由氢原子与电负性较大的原子如氧原子、氮原子或是氧原子、氮原子的孤对电子形成的。一般氢键的键能在56kcal/mol(1cal=4.1868J,1kcal=4186.75J),而单键的能量一般在60110kcal/mol,vdw相互作用能在0.10.2kcal/mol。氢键的能量主要来于自带正电荷的氢原子与带负电荷的其它原子或孤对电子之间的静电吸引。这种特殊的稳定作用可以通过较大的势阱深度和较小的vdw半径来模拟,而且氢键的vdw作用方程也略有不同,通常使用改进的Lennard-Jones方程:
在一些力场中,还在和氢键距离有关的相关项中加入角度部分(1cosXHY)或是(1cosXHY)4来计算氢键的vdw相互作用。目前对氢键的处理倾向于不特殊考虑其vdw参数或方程形式,而仅从静电相互作用方面来解释。
VanderWaals相互作用的参数可以通过多种途径获得。在早期的力场中,一般是通过分析晶体堆积获得参数,这样的vdw参数可以得到非常好的实验几何结构以及热动力学性质如升华能。近期的力场是通过模拟液态体系获得vdw参数,通过优化这些参数可以计算更广泛的液态体系的各种动力学和热力学性质。
1.2.6静电相互作用非键相互作用的另一部分就是静电相互作用,静电作用是由于分子内部电荷的极化,产生的带有正电荷和负电荷部分的相互作用。计算静电相互作用的一种方法是将电荷分配在原子上,另一种方法是将键上分配一定的偶极矩[27~28],虽然这两种方法给出相似的结果,但并不完全相同,只有在计算分子间长程相互作用时,这两种方法才给出相同的结果。
点电荷的相互作用可以用库仑势能函数(CoulombPotential)表示:
其中,是介电常数,原子电荷QA,QB通常为拟合参数。由于氢键主要由缺少电子的氢原子和电负性较大的原子(如氧原子和氮原子)之间的静电吸引作用,合理的选择部分电荷则可以比较准确地模拟这种相互作用。有效的介电常数可以模拟周围分子(溶液)的极化作用,值为1对应的是真空介质,较大的值削弱了电荷-电荷之间的长程相互作用。尽管没有理论上的证明,通常值取为1~4。在一些力场中,介电常数是依赖于距离变化的0ABr,库仑相互作用则表示为:2AB0ABQQ/(r)。
类似于vdw相互作用,传统的静电相互作用仅包括两体作用。然而对于极化体系,三体相互作用则非常重要,大概占两体作用的10%~20%[29]。三体相互作用是由于第三个原子的极化作用导致的两个相互作用的原子电荷的改变,“多体”效应可以用原子的极化来模拟[30]。由于极化作用大大地增加了计算时间,因而其应用受到了很大的限制。对于极化体系忽略极化作用可能是现代分子力学中一个主要的缺陷,但是一些力场在选取电荷参数时考虑了平均极化作用,使得原子的点电荷要略大于孤立分子时原子上的电荷。
电负性大的元素比电负性小的元素更容易吸引电子,在分子中产生不均衡的电荷分布。这种电荷分布可以用多种方式来表达,一种普遍的方法就是在整个分子内进行部分电荷重排,这些电荷可用来描述分子的静电相互作用。如果电荷被限制在原子核周围,通常称之为部分电荷模型或净电荷模型。在两个分子(或是同一分子不同部分)之间的电荷相互作用可用库仑定律来描述:
其中NA和NB是两个分子中点电荷的个数。如果点电荷只分布在原子上,那么原子点电荷参数可以通过四种方法来获取:第一种,通过半经验方法,如利用电负性均衡方法[119~120,134~135,140,175];第二种,通过拟合量子化学计算分子波函数布居分析方法,如Mulliken布居分析[33]、AIM(AtomsinMolcule)方法[36]、Lwdin方法[35]和自然轨道布居分析方法(Naturalpopulationanalysis,NPA)[36]等;第三种,拟合从头计算分子静电势电荷方法[37];第四种,拟合实验偶极矩方法[31,32]。目前在力场中常用的是拟合静电势电荷方法,该方法可以得到比较合理的静电势能。当然,在实际中常常考虑多种因素来拟合点电荷参数,如静电势电荷、分子偶极矩、分子团簇的结合能和分子的热力学数据等。
另一种计算庞大体系的分子间静电相互作用的方法是中心多极展开法(CentralMultipleExpansion)[38]。这种方法是将整个分子视为一个整体,利用分子的电次极计算静电相互作用。分子的各电次极分别为:零次极为分子所带的电荷,二次极为分子的偶极(Dipole),四次极为分子的四极(Quadrupole),八次极为分子的八极(Octopole)等。通常最重要的是分子最低不为零的次极。非中性的离子或分子,如Na+、Cl-、NH4+、CH3CO2-等均带有电荷。中性的分子如H2O、NH3、CH3Cl等最低的电次极为偶极,如CO2、N2等分子最低的电次极为四极。CH4和CCl4等分子最低的电次极为八极。各种电次极可以用适当的电荷排列表示,如偶极可表示为两个带相反电量的电荷排在适当的距离。下面简单地介绍中心多极展开法。
设一个分子中带有q1和q2的电荷,分别位于-z1和z2的位置,如图1-8所示。
的四极,依此类推。由上式可知,P点的势能可展开为分子多极的多项式,称之为中心多极展开法。此展开式中,仅第一不为零的次极项与坐标无关。若总电荷为零,则偶极不变;若总电荷与偶极均为零,则四极与选择的位置无关。通常在计算中选取电荷分布的中心为原点。
利用中心多极展开法计算静电相互作用的优点在于该方法具有较高的效率。以苯分子为例,若以原子所带点电荷模型计算两个苯分子的静电作用,需要计算144个相互作用项。但苯为无电荷且无偶极矩的分子,最低不为零的次极矩为四极矩。
因此,以中心多极展开法计算两个苯分子间的静电作用仅需要计算一个四极矩-四极矩作用项。当两个分子间的距离与分子的尺度相当,则不能引用中心多极展开法。
只有当分子间的距离大于各个分子的中心至其最远电荷距离的和时,才能应用中心多极展开法。
顾名思义,非键相互作用项对成键势能项中的原子间的势能应当没有贡献,如两个直接相连的原子间无Evdw和Eelec相互作用,它们之间的相互作用主要用Estr来描述;对于分子CH3(CH2)50CH3,两端的H原子之间的相互作用等同于不同分子中H原子间的相互作用,因此描述这类分子的力场函数中就应该包含Evdw/Eelec项。但是在实际中,非键相互作用项的应用范围不尽相同,大多数力场中当间隔三个或更多化学键时就包含了原子对间的Evdw/Eelec项,这就意味着在描述A-B-C-D型分子的力场中对于A和D原子对间的相互作用既包括Etors,又应该包括Evdw/Eelec项。或者在一些力场中有1,3位的vdw相互作用(也就是A-B-C分子中A和C间的vdw作用),这样的力场称为Urey-Bradley。在这种情况下,扭曲三个原子相连的分子所需要的能量包括了Evdw和Ebend项。在现代分子力场中,Estr计算分子中所有相互键连的原子对间(1,2)的伸缩振动,Ebend计算间隔一个化学键的原子对间(1,3)的弯曲振动。类似的,Etors计算1,4间的扭转运动,Evdw/Eelec描述1,4之间或是间隔更多化学键的原子对间的相互作用。
1.2.7各种相互作用的耦合
方程(1-4)中的前五项在通常的力场能量表达式中是普遍存在的。最后一项Ecross,表示前面各项的相互耦合势能项。以水分子为例,平衡构型的键长为0.958,键角为104.5。如果角度被压缩到90,由电子结构方法优化得到的最低能量构象的键长为0.968,略长于平衡几何键长。类似的,如果角度变大,则最低能量构象的键长就小于平衡几何键长。可以定性地理解为,角度变小氢原子就会靠得更近,这将导致氢原子之间的排斥变大,从而增大了键长。如果在力场能量表达式中只有前五项,就不能很好地模拟键长和键角之间的耦合,这样就可以考虑加入一个额外的即依赖于键长又依赖于键角的相互耦合项。
Ecross通常可以写为单个坐标的Taylor展开式。在所有的耦合相互作用项中最重要的就是键长和键角的耦合,对于A-B-C分子,键长和键角的耦合可以表示为:
通常情况下,耦合相互作用项中的常数并不是依赖所有相关联的原子类型。
比如,对于A-B-C分子,原则上stretch/bend耦合作用项的参数应该包括A,B,C三个原子,但是通常只考虑中心原子,也就是kABC=kB,或者不依赖原子类型选取一个普遍适用的常数。必须注意到,如果分子的几何构型远离平衡几何结构,上面的耦合相互作用不是很稳定。比如,化学键被无限拉大,如果小于0则Estr/bend将会趋向-。如果键长的伸缩振动势能用简谐振子描述就不会产生这样的问题,但是如果用Morse势能函数,在几何构型优化和模拟中就要特别注意避免键长过大。
1.3分子力场的参数化
力场方法的一个最基本的假设就是结构单元可以在不同的分子中相互转移,这就给力场的参数化提供了可靠的依据。也就是说,只要确定力场中原子或基团的类型,从实验或量子力学计算结果中选取恰当而又充分的相关数据,就可以得到可以相互转移的力场参数。利用电子和X光衍射可以得到键伸缩势能函数的平衡距离;利用红外和拉曼振动光谱可以得到相关的力常数;同样,势能函数和其参数的质量也需要用实验来验证,例如汽化热,溶解热,分子密度和扩散系数等。从头计算方法的发展和应用对力场的参数化起到非常重要的作用,利用从头计算方法获取参数有很多优点:(1)可以为参数化提供比实验更加充足的数据;(2)能够使力场应用于更加广阔的领域和分子类型;(3)使其参数更具有可转移性和一致性。力场的参数化是一个非常重要而又烦琐的工作,当增加一些新的参数来计算新的目标分子的性质时,需要大量的数据来拟合力场参数。对于烷烃分子来说,如果采用一般的力场函数,需要28个不同的参数。在MM2力场中,对30个不同原子类型就需要3722不同的参数。对键长伸缩振动和键角弯曲振动势能项参数来说,在不同的力场中具有一定的转移性。但对非常敏感的非键相互作用和二面角扭转势能项,由于它们之间存在着较强的耦合,所以要花费大量的时间优化非键相互作用和二面角扭转作用的参数。
分子光谱很难利用力场方法模拟,可以适当地考虑耦合相互作用项Ecross,则可以大大改善其计算结果。一些力场的参数化过程充分地利用分子模拟方法得到的热力学数据,如OPLS(OptimizedParametersforLiquidSimulations)力场[39],从而使其非键相互作用参数能够更加准确地模拟体系的分子间相互作用和相关液态性质。