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

第2章 前言(2)

为简单地描述由原子组成的分子体系的非谐振动,通常根据分子的各种结构单元、势函数分解成键长伸缩能Estr、键角弯曲能Ebend、二面角扭曲能Etors、范德华作用能Evdw、静电作用能Eelec和以上各项的耦合相互作用能Ecorss等,总能量EMM可表示为:

对于更精确的力场势函数,计算光谱等性质,常常加入耦合相互作用项Ecross,其可以包括键长伸缩-键角弯曲交叉项、二面角扭曲-键长伸缩交叉项等,同时也可加入氢键函数项等其他修正项。图1-1形象地体现了各种结构单元对应的能量项,下面将对上述各势能项分别讨论。

1.2.1键长的伸缩振动

分子中最强烈相互作用的原子形成化学键,通常情况下,化学键的长度在其平衡位置附近呈小幅度的振动。描述此种作用的势能项被称为键伸缩振动势能项。键伸缩势能函数的最简单的表示就是在平衡位置处进行Taylor展开,其截止到二级展开项的表达式为:

0R为A-B键的平衡几何键长,由于力场方法通常得到的是相对势能,所以E(0)可以作为常数,而从势能函数中去掉。在求一阶微商时,上式中第二项也为零,因而,化学键A-B的伸缩振动能可以用谐振势能函数表示:

如果化学体系的张力很大(strainedandcrowded),用谐振势能函数不足以模拟键的伸缩振动,用谐振函数计算得到的结果与实验值也有很大的偏差。这就需要对谐振函数进行改进,最直接的方法就是加入高阶非谐振动项以校正其误差:

加入高阶项虽然能在一定程度上更好地描述键的伸缩振动,但同时也增加了需要拟合的参数,而且,键长伸缩振动势能函数的高阶展开不能带来正确的极限行为,立方非谐振常数k3通常是负值,如果Taylor展开在立方项截断,则键长较大时能量值趋向于负数,而四次项常数k4通常为正值;如果在四次项截断,则键长较大时能量趋向正数。利用这样的函数求势能的极小值,如果初始结构不好就得不到正确的最优几何构型。化学键伸缩趋近无限大的正确极限是化学键的断裂,也就是能量趋近于离解能。

能够比较正确描述这一极限行为的函数为Morse势能函数[20]:

Morse势能函数并不经常使用,因为它需要三个参数来描写一个化学键,但是Morse势能函数能够非常好地描述从平衡位置到离解状态较大范围内化学键的行为。图1-2描述了各种势能函数在平衡键长附近的势能曲线。所以,从总体来说,Morse势能函数在较大的键长变化范围内都能很好地描述键的伸缩势能变化。

谐振模型计算结果同实验相比较示意图

键伸缩势能函数中的平衡键长参数一般用高分辨率X射线结晶学方法测定,精度可达0.0001nm,也可由电子衍射和微波光谱测定。力常数可由振动光谱得到。优化力场参数时,一般选择一定数量的模型化合物,根据原子和键所处的化学环境来进一步的修正,使得力场参数不但能重现分子的结构和能量,而且能够较好地模拟分子的振动光谱,然后再计算其他没有作为模型化合物的分子性质来验证这套参数的可靠性。

1.2.2键角的弯曲振动

与键长伸缩振动类似,胡克定律或是简谐振动也可以用来描述键角的弯曲振动,即:

每一个角度对势能的贡献也是利用常数和参考值来描述,如果从平衡位置弯曲一个角度需要较小的能量,则力常数相对来说也较小。张力比较大的有机化合物与无张力的有机化合物之间,键角的差值最大为15,在这种情况下,非谐振动现象不明显[21,22]。当键角过分小时,如含有三元环和四元环的化合物,非谐振现象就很明显,这时谐振模型就要做相应的修正。

与键伸缩振动类似,加入高阶项是常用的修正方法。

1.2.3二面角扭转相互作用

键的伸缩和键角的弯曲通常被称为“硬自由度”,如果要是它们偏离平衡位置很远往往需要非常大的能量。通常,分子中结构和能量的变化是由二面角扭转和非键势能项共同作用的结果。分子中连续存在相键连的四个原子A-B-C-D,二面角(torsionalangle)是指A-B和C-D键组成的角ω(如图1.3所示),角度ω的取值范围在[0,360]或是[-180,180]。二面角扭转势能与键长的伸缩振动和键角的弯曲振动有两点显着的不同。首先是能量函数是以ω为周期的周期性变化函数,也就是当ω旋转360后的能量和0时的能量相同。其次是二面角的扭转能一般很低,因而偏离最低能量构象的结构很容易发生。在对构象的研究中,二面角扭转势能就是非常重要的。通常用二面角ω的余弦函数的Fourier展开来描述其能量,给出周期性变化的描述。

nV通常指势能垒的高度,但是不确切,因为在展开中,其他项和非键项也对绕中心键图1-3中的2-3旋转起重要作用。n表示多重度,它的值表示在区间[0,360]内该势能函数给出的极小点个数;是相位因子。当n=1时,旋转周期是360;n=2时,旋转周期是180,极小点在0和180;n=3时,周期是120,极小点在60,-60和180,极大值在120,其他类似。

1.2.4非共面弯曲振动

实际分子中,键角的弯曲振动除了平面内振动外,还包括平面外振动。在力场的应用中发现,一些体系通过一般的势能函数不能正确地得到分子结构,如甲基酮等羰基化合物分子,只用键伸缩和键角弯曲势能项得不到实际的平面结构,尽管平面结构使键角弯曲项能量升高,但是π其键离域能更低而使之最终处于相对能量更低的平面结构。这种情况下,就需要在力场中加入非共面弯曲势能项(或非正常二面角扭转势能项)来保证其得到合理平面结构。如图1-4所示,在环丁酮分子中定义ω非正常二面角扭转势能函数项,ω表示原子1与5-3-2平面所成角,用以上方程形式,以保持其二面角在0或者180。在一些力场中还加入了键长伸缩振动-键角弯曲振动的交叉项。目前主要的方法有三种:一是,利用非正常二面角扭转势能项,用方程(1-12)表示:

这种模型定义更容易理解,首先定义一个平面,如图1-5所示,该平面由中心原子和其他两个原子确定,第四个原子与该平面所成的角为,如果四个原子在一个平面内,则角为0°。模型三中的h表示第四个原子与上述平面的距离。通过以上的校正可以保证力场得到合理的分子结构和能量。

1.2.5VanderWaals相互作用

VanderWaals相互作用能描述除以上键连情况外的非键连的两个原子间的排斥和吸引作用,和静电作用一起称为非键相互作用。在原子间距离较大时,Evdw为零;而在原子间距离较小时,Evdw表示出较大的排斥作用。在量子力学中,这是由于两个原子的电子云相互重叠,带负电荷的电子间相互排斥而引起的。在一定的距离时,这两个电子云之间由于电子相关效应又存在着较弱的相互吸引作用,这种吸引作用也可以解释为诱导的偶极-偶极相互作用。即使分子(或分子的某一部分)不存在永久偶极矩,电子的瞬时运动也会产生不均匀的分布,偶极矩将会在邻近分子(或者分子的不同部分)产生诱导偶极矩,这样就产生了吸引作用。从理论上说,这种吸引作用随着两个原子之间距离6次方的倒数的变化而变化。事实上,诱导偶极-偶极相互作用仅仅是从诱导的四极-偶极、四极-四极等相互作用的贡献中获得的,并且随着r-8,r-10等的变化而变化。r-6仅仅是在较远距离时的渐近行为,因而这种相互作用力也可称之为“色散力”或“London力”[17]。Vdw相互作用力是稀有气体之间的主要相互作用(这就是为什么氩可以成为液体和固体的原因),并且它也是非极化分子如烷烃之间的主要相互作用。

在较小的距离时,Evdw表现为较大的正值;当两个原子恰好相互“接触”时,Evdw表现出较小的负值;当距离再增大时,Evdw趋近于零。满足这种关系的一般方程形式为:

从理论上说,排斥作用的具体表达形式是不可能获得的,但是,它仅要求当r趋向于无限大时,排斥作用趋于零,并且它要比r-6更快地趋于零。能够满足上面要求的比较通用的势能形式是Lennard-Jones(LJ)函数[23]。其中排斥作用项用r-12表示:

Lennard-Jones势能函数中有两个可调参数:碰撞直径(在这个距离时相互作用能为零)和势阱深度,图1-6给出了这些参数的图形表示。Lennard-Jones能量表达也可以用rm表示,在这个距离时,势能函数有一个极小值,并且能量对原子核间距离的一阶微分等于零(也就是:/0ELJr)。因而,可以很容易地得到表达式:6mr2。

Lennard-Jones势能函数也可以表示为:

Lennard-Jones势能函数的特点是有一个相互吸引部分,它随着r-6的变化而变化,而相互排斥部分随着r-12的变化而变化,这两个部分画在图1-7中。吸引相互作用部分r-6类似于Drude模型中色散能的理论处理。虽然排斥部分r-12没有理论上的讨论,尤其是量子力学中没有建议这样的指数项,但是对于稀有气体指数12给出非常合理的计算(对于碳氢体系指数12给出的势能函数太陡)。目前,Lennard-Jones势能函数被广泛的使用,尤其是计算大分子体系,原因之一就是从r-6很容易得到r-12。在计算排斥相互作用时也可以用到不同的指数表达形式,比如在一些力场中9或者10就可以给出不太陡的势能函数。因而,vanderWaals相互作用的一般形式可以写为:

从电子结构的理论上讲,排斥作用是由于电子波函数的重叠,而且随着逐渐远离原子核,电子密度呈指数形式降低(氢原子的确切波函数就是e指数形式),这似乎也给出了一些关于排斥部分选择e指数形式的证明。“e指数-r-6”的Evdw相互作用形式,也称之为“Buckingham”或“Hill”形式的vdw函数:

其中,r0和的物理意义和前面相同,是自由选择的参数。如果选择=12,则描述长程相互作用时就等同于Lennard-Jones势函数;如果=13.772,就给出平衡位置时Lennard-Jones势能函数的平衡力常数,参数也可以利用拟合常数。(1-20)式给出的e指数-6形式的势函数在处理较近原子之间的相互作用时会产生一些问题,比如在什么距离该势函数会发生“翻转”(turnover)。当r趋近于零时,指数部分趋近于常数A,而r-6项趋近于-。因而当优化几何结构时,两个原子距离很近就会导致原子核的聚变(fusion)。

Evdw的第三种表达方法是用Morse势能函数表示,在长程相互作用时,Morse势能函数不包含r-6项,但是正如前面所讲,它实际上包含了r-8和r-10项。描述Evdw的Morse势能函数中D和远比Estr中的小,而且r0要更长一些。

三种描述Evdw的方程主要区别是在近程相互作用的排斥部分,对于传统力场,Lennard-Jones和e指数-6势函数都过高地估计了近程排斥作用,而且还存在近程“反转(inverting)”问题,从化学角度考虑这些问题,它们似乎都无关紧要,因为超过100kcal/mol的能量足以破坏化学键,在实际计算中也从不在这部分取样。目前,对于极化力场势能函数,在形成氢键原子间静电吸引项在近程极化升高,甚至导致“极化塌陷”,这恰好与e指数-6形式的vdw势函数行为相反,因而在极化力场中有很好的应用价值[103~104]。势函数的相互吸引部分是分子间相互作用最基本的组成部分,这三种Evdw的表达式基本一致。Morse和e指数-6势函数形式都给出较好的Evdw,u36890矼可能是由于这两种方法有三个可调参数,而在Lennard-Jones势能函数中只含有两个参数。

尽管指数方程存在一些缺陷,但大多数的力场都采用Lennard-Jones形式的势能函数,因为Lennard-Jones势能函数不但可以给出相对合理的结果,并且计算简单。

多原子体系的vdw相互作用不可避免地包含了不同原子类型的相互作用。比如利用两点模型计算两个CO分子之间的Lennard-Jones势能,vdw参数不仅包含了C-C和O-O之间的作用,同时也需要包含C-O之间的相互作用。如果一个体系中有N个不同类型的原子,则需要N(N1)/2组参数来计算不同原子的相互作用。vdw参数的确定是一个很困难也很耗时的过程,通常假定不同类型原子间的参数是由相同类型原子间的参数拟合而得的。Lorentz-Berthelot拟合原则给出A-B间相互作用的碰撞直径AB等于A-A和B-B间碰撞直径的算术平均值,而势阱深度等于几何平均值:

Lorentz-Berthelot的拟合原则应用于相似的分子时给出非常好的计算结果。它的主要不足是利用几何平均过高地估计了势阱深度。在一些力场中也可以利用几何平均来代替算术平均计算碰撞直径,或是利用算术平均代替几何平均计算势阱深度。

如果分子中两个原子相距三个化学键(A-B-C-D中A和D原子,也称之为1,4原子),它们的vdw和静电相互作用不同于其他的非键相互作用,要进行特殊的处理。这些原子的相互作用对于中心键的旋转势能以及扭转势能都能有很大的贡献。1,4非键相互作用通常需要乘以经验校正因子,比如,1984年的AMBER力场[24]中建议的vdw和静电相互作用的校正因子为2.0(1995年的AMBER力场中静电部分的校正因子为1/1.2)。对于1,4非键原子的校正原因主要是:vdw中排斥部分采用r-12的指数形式(相对于e指数形式,它过于陡峭),对1,4原子有显着的影响。另外,当两个1,4原子距离很近时,电荷通过相连的键发生了重新分布,这将部分降低相互作用能,因为如果具有相同距离的两个原子处于不同的分子中,电荷是不可能发生这样重新分布的。