分子力学与分子动力学

1、分子力学和分子动力学Molecular Mechanics and Molecular Dynamics甄云梅20111210内容简介 力场1 分子力场的势函数形式2分子力场分类3 力场的参数化能量最小化热力学性质分子力学简介分子力学简介分子力学，又叫力场方法(force field method),是基于经典牛顿力学方程的一种计算分子的平衡结构和能量的方法。与量子力学不同，它求解的是Newton方程，而不是薛定谭方程。基本假设: 1、The Born-Oppenheimer Approximation对势能面的经 验性拟合。量子力学中的薛定谭方程为：(非相对论和无时间依赖的情分子力学简介基本

2、假设基于Born-Oppenheimer近似，其物理模型可描述为：原子核的 质量是电子质量的103J05倍，电子速度远远大于原子核的运动 速度，每当核的分布形式发生微小变更，电子立刻调整其运动 状态以适应新的核场。这意味着，在任一确定的核分布形式下，电子都有相应的运动 状态；同时核间的相对运动可视为所有电子运动的平均结果。 所以电子的波函数只依赖于原子核的位置，而不是他们的动能。 于是这个近似认为，电子的运动与原子核的运动可以分开处理, 可以将上式分解为/电子运动/ 的波函数厂电子运动方& (Hel + VNN )y/e (R, r) = EeIy/e (R, r) 核运动方程:(7； + E

3、el )y/N (R) = E屮“ (&核运动的波函数A电子运动方程：(比严VnQ%(R,门二Em(R，门 核运动方程:忆+岛必(& = Eg(R)方程中的能量Ee(势能面)仅仅是原子核坐标有关。相应的，方程 所表示的为在核势能面E(R)的核运动方程。直接求解方程，釆用的是从头算或者是半经验，这样的量化计算都是把电子的波函数和能量处理成原子核坐标的函数。由于量子化学求解电子波函数和势能面耗时巨大，常常将势能面进行经验性的拟合, 成为力场，由此构成分子力学的基础。将方程用牛顿运动方程代替，势能面釆用力场拟合，就构成了分子 动力学的基础。分子力学简介基本假设 2、简单作用模型对体系相互作用的贡献来

4、自诸如键伸缩、键角的开合、单键的旋转等等。即使使用类似Hooke定律这样的简单函数，也 能令力场运转良好。口 3、力场的可移植性（关键属性）仅在少数情况下通过测试的一套函数，可以用来解决更广范围内的问题。进一步讲，从小分子得来的数据可以用来研究 类似高分子的大分子。简单分子力场分子力场是分子力学的核心。分子力学的基本理论就是一个分子力场由分子内相互作用和分子间相互作用两大部分构成, 即力场的势能包括成键和非键相互作用，所有的势能的总和即 为分子的构象能。E(= Z 耳 + E耳(&)b ondednonbon ded分子力学简介简单分子力场由于分子力学是经验的计算方法，不同的分子力学方法会釆用

5、不同的势能函数(Potential Energy Function, PEF)表达式，而且力场参数值也会不同。般将分子的PEF分解成五分子力学简介分子力学简介部分:v(r”)=工与a -+ 工与(Q _ 0加2 + Y y(l + cos(nG-/)bonds 厶angles 厶torsions 厶分子力学简介分子力学简介分子力学简介EboMrorate- cilong -bond=F+ Fbond -stretchcnigle-bend分子力学简介分子力学简介键伸缩能键弯曲能二面角扭转能non-bonded=Evcm-der-Waalselectrostatic范徳华作用能静电作用能分子力学

6、简介简单分子力场这样的一个简单的力场方程如何用来计算分子的构象能?Angles CCC x 1 C-C-H x 10 H-C-H x 7Non-bonded H-H x 21 H-C x 6V(/-：=务/+齢U+Ixtndxfifigtesv工-y(l + COS(H6t-/)+壹（4昇尸-（如+潭） tl%G4亦必.分子力学简介简单分子力场 General Features Force fields are empiricalThere is no “correct” form of a force fieldForce fields are evaluated based solely

7、on their performance. Force field are parameterized for specific propertiesStructural properties.Energy.Spectra Force field definitionFunctional form (usually a compromise between accuracy and computational efficiency) parameters (transferability assumed). Atom typesAtomic number (e.g., C5 N, O5 H).

8、Hybridization (e.g., SP3, SP2, SP).Environment (e.g.，cyclopropane, cyclobutane).分子力学分子力场的势函数形式键伸缩能分子力学分子力场的势函数形式分子力学分子力场的势函数形式键伸缩能0WVOVb= X撷(IMBindungenVa= E討川-%尸Winkcl952 寺ExtraplanarwinkelVo =52心【1 + cos(呼-d)DihcdralwinkclVk/W-分子力学分子力场的势函数形式键伸缩能%= EQ刚/(4低0好)Paare(i j)VvdW = 龙c12(ij)/r； 一 C如gAPaarc

9、(ij)E$ = De exp(-A(/ - /0) -12De： depth of the potential energy minimum和/2DeM： messo)： frequency of the bond vibration (w= 7k/jj) l0: the reference value of the bondOptimum莫斯函数(Morse Function)Es = De exp(-A(/ - Zo) -12TRIPOS, Cherm-X, CHARMm和AMBER釆用谐振 子函数形式CVFF, DRIEDING和UFF既支持莫斯函数也支持谐振 子模型MM2和MMX用二

10、阶泰勒展开的莫斯函数MM3, CFF和MMFF94用三阶泰勒展开的莫斯函数分子力学分子力场的势函数形式键角弯折能分子力学分子力场的势函数形式键角弯折能键角弯折能Angle Bending谐振子模型1 9EB = -kfe-e平衡键角键角键角弯折力常数谐振子模型在偏离平衡位置不大的情况下（10以内）可以取得很好的结果。采用谐振子的力场包括:TRIPOS, CHEM-X, CHARMm,AMBER以及CVFF等分子力学分子力场的势函数形式二面角扭转能二面角扭转能 Torsion Rotationet =+ COS（Ht /）分子力学分子力场的势函数形式二面角扭转能分子力学分子力场的势函数形式二面角

11、扭转能匕为势垒高度(barrier height),定量描述了二面角旋转的难易程度； N为多重度(multiplicity),指键从0。至U360。旋转过程中能量极小点的个数; 丫为相因子(phase factor),指单键旋转通过能量极小值时二面角的数值。 co为扭转角度(torsion angle)大部分力场如AMBER，TRIPOS, CHEM-X, CHARMm, COSMIC, DREIDING和CVFF采 用较简单的势函数形式第二代力场如MM2, MM3, CFF及MMFF94 采用傅里叶级数形式由于二面角的扭转对总能量的贡献小于键长和键角的贡献, 一般情况下二面角的改变要比键长和

12、键角的变化自由得多。 因此在一些处理大分子的力场中常保持键长、键角不变， 只考虑二面角及其他的作用而优化整个分子的构象和能量。分子力学分子力场的势函数形式静电相互作用范德华相互作用能Lennard-Jones势函数w=（-r-（）nr n rA6F1Jr.12ruRepulsion regimevan der Waals attraction regimeOptimum energyI1r为原子对间的距离；w为势阱深度， 为势能参数，因原子的种类各异。正的部分为排斥势，负的部分为吸引势取6,加取12时，叫做LJ 6-12势函数，用于AMBER, CVFF, CHARMm,DREIDING, U

13、FF以及TRIPOS等力场静电相互作用 Electrostatic Contributions点电荷法：通过经验规则或者量化计算确定每个原子上的部分电荷(partial charge),两个原子之间的静电作用用库仑公式来计算。分子力学分子力场的势函数形式静电相互作用分子力学分子力场的势函数形式静电相互作用偶极矩法：根据某些规则计算出每个化学键的偶极矩，通过计算偶极偶极相互作用來描述静电相互作用。dipoleVdipole = K 巴(cos 力 - 3 cos cos a) rij匕切。/產分子间或分子内偶极偶极相互作用的能量 严是两个偶极的偶极矩 /是两个偶极矩间的角度门和 是连接两个偶极向

14、量间的夹角点电荷法VS偶极矩法两种方法在处理有机小分子体系的时候效率相似，但 是当用来处理带电生物大分子体系时，偶极矩方法显 得过于耗时。 MM2, MM3和MMX用键偶极矩法计算静电相互作用其它力场釆用点电荷方法计算点电荷方法的问题在于如何把电荷分配到原子上量子化学计算法电荷可以由多极矩、热力学性质、静电势拟合得来经验规则法分子力学分子力场的势函数形式交叉相互作用项Crossing Terms键伸缩-键伸缩相互作用键伸缩-键角弯折相互作用键伸缩-二面角旋转相互作用键角弯折键角弯折相互作用键角弯折-二面角旋转相互作用应用 TRIPOS, CHEM-X, AMBER, DREIDING, UFF

15、和 COSMIC力场中没有相互作用项 MM2和MMFF94只支持键伸缩-键角弯折相互作用项 MM3力场支持键伸缩-键角弯折、键角弯折-键角弯折、 键伸缩.二面角旋转相互作用项 CVFF和CFF91都支持能量是相对的由不同的方法计算得到的能量的绝对值是毫无意义的。只 有当它与同体系的其他构象计算得到的能量相比较时才有 意义。-比较不同程序计算得到的能量值无意义-用同一种程序时，比较不同分子的能量值无意义分子力学分子力场分类传统力场分子力学分子力场分类传统力场传统力场 AMBER力场 Kollman group, 1984最初仅为蛋白质和核酸体系提供相应的原子类型和力场参数 1990,发展了适用于

16、多糖模拟的力场参数(Homan 1990) 1995, 2000,加入了适用于有机小分了的原了类型和参数 2torsions 厶E =- -)2 + y (0+ 工 才l + cos(6 刃anglesbonds 厶cne/x 厶f/ 、12/ 6QiQj +r17乙丿r4亦/ IJ+工 CHARMm力场(Chemistry at Harvard Macromolecular mechanics) Karplus group, 1983适用于各种分子性质的计算和模拟对于从孤立的小分子到溶剂化的大生物体系的多种模 拟体系都可以给出较好的结果但不适合于有机金属配合物e+2X(力-龙。+工詈7 +Z

17、ij &儿匕0冷 ijE二工心(厂-汀+工觴(0_%)2+工卄对心cos(M)b&06 J ?丁 人口 B :乙 T12- CVFF力场(Consistent Valence Force Field) Dauber-Osguthorpe group, 1988适用范围包括有机小分子和蛋白质体系扩展后可用于某些无机体系的模拟，如硅酸盐、铝 硅酸盐、磷铝化合物主要用于预测分子的结构和结合自由能E =$1 -严fY 弘(0 - 尸 +乙1 + cos(M) + 工叱b0(fZ+工工必)+工工伤卫-仇)0-%)b b90 夕+工工坨2-%)(%) +工尸则cos0(&-%)-%) +工工 加b &0Z

18、 / MMX力场 Allinger group, 1989包括MM2和MM3主要针对有机小分子函数形式比较复杂，包含交叉项也可用于生物大分子体系，但是速度会比较慢E =Es + 工 Eb + 工 Ey + 工Eop + 工 Esb + 工 Ey + 工 E“+工 Eq分子力学分子力场分类第二代力场第二代力场 CFF力场 CFF91主要适用于模拟有机小分子、蛋白质以及小分子蛋 白质之间的相互作用 CFF95除了适用于蛋白质和有机小分子体系，还可用于有机 高分子体系的模拟，如聚碳酸酯及多糖PCFF在CFF91的基础上，还适用于聚碳酸酯、三聚氧胺甲 醛树脂、多糖、核酸、分子筛等其他无机和有机材料体系

19、的 模拟。同时还提供20种金属离子的参数分子力学分子力场分类第二代力场 COMPASS力场（Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies用于原子水平模拟研究的凝聚态优化的分子力场） Sun group, 1994第一个能同时预测气态和凝聚态性质的从头力场适用于常见的有机小分子、无机小分子和高分子最新的力场也适合于金属、金属氧化物及金属卤化物的模 拟可较为准确的预测晶体的多种性质分子力学分子力场分类第二代力场 MMFF94力场(Merk Molecular Force Field)

20、Hagler, 1996目前最为准确的力场之一定义了非常完备的原子类型既适用于有机小分子，也适用于大分子体系，如蛋白质E+工Eb +工Et十工Eop +工Esb十工Ev +工E“+工分子力学分子力场分类分子力场的选择通用力场基于规则的力场(Rule-based Force Fields)较好的可移植性，适用范围广精度差 ESFF力场 (Extensive Systematic Force Field)用于有机分子、无机分子及有机金属化合物的结构预测 UFF力场(Universal Force Field)覆盖了周期表中所有元素，应用最广泛的通用力场 Dreiding 力场适用于有机分子、生物大

21、分子及位于主族的无机分子分子力场的选择蛋白质分子的模拟：首选AMBER力场和CHARMm力场，也可 用CFF力场、CVFF力场和MMFF94力场核酸分子的模拟：采用AMBER力场、CHARMm力场、 MMFF94力场或用户自定义的力场小分子蛋白质复合物体系的模拟：首选CHARMm力场和MMFF94力场,也可用CVFF力场和CFF力场高分子的模拟：首选COMPASS力场，也可用PCFF力场和CFF95力场分子力学力场的参数化分子力学力场的性能即它的计算结果的准确性和可靠性主要取决于势能函数和结构参数。这些有关力常数，结构参数的“本征值”的置定过程称为力场的参数化。参数化的过程要在大量的热力学、光

22、谱学实验数据的基础上进行，有时也需要由量子化学计算的结果提供数据。各类键长、键角的“本征值” 一般取自晶体学、电子衍射或其 他的谱学数据，键伸缩和角变力常数主要由振动光谱数据确定, 扭转力常数经常要从分子内旋转位垒来推算。分子力学能量最小化分子力学最重要的内容是根据适合的力场计算分子各种可能构象的势能，势 能最低的构象为最为稳定的构象。寻找势能最低点的过程称为能量最小化， 所得到的构象称为几何优化构象。分子的几何优化构象是计算分子性能的基 础。由于数学上只能保证求得局部极小值，即实现局部优化，而不一定能求得全 局最小值。所以得到的是在这一构象附近的一相对稳定的构象。这里介绍的 能量极小化方法只

23、是求局部极小的方法。粗结构分子几何优化能量极小构象= F = 0; 4odXi dx-分子力学常用的优化方法最陡下降算法共辄梯度算法Steepest Descents SD 一级微商算法Conjugate Gradients - CONJ二级微商算法 Newton-Raphson Method能量极小化算法-最陡下降法(SD)最速下降法求能量极小点的原理如卜一：很据能虽 梯度的负值-V(g)指出的方向进行一维搜索达到极 小,然后樽求笫二次的TE(qh再沿它指出的方向逬 行一维搜索到极小.继续卜。去直到符合预设的能量收 敛标准.于是在势能面负梯(91也【=T(N523)以二维问题为例(图2 5.

24、2-1)：设起始状态点在三护 必) 度自线方向上任意点的位置为接着求得沿着梯度方向上任意点的能暈即得到图2.521(a).其中.满足n严&(/)的点就是图2.5.2.1(a)中的C点.这样就完成了最速下降法的第1步. 城后从C点出发作最速下降法的第2步一不断进行下去完成全部最速F降法(图2.5.2-Hb).0)辰速下降袪的一址搜索(b)收图2.5.2-1加速下降法的 维搜索及K收敛路径可以证明最速卜降法的能鼠优化搜索路径中的每一步的方向都垂貢于前一步的方向.:能量极小化算法-共轨梯度法(CONJ)基本思想是把最速下降法与二次函数共純方向结合起来、利用L2知位置处的梯度构 造一组共轨方向、并沿着

25、这一组方向进行一维搜索，求出目标函数的极小点存储量少，计算方便。IT. r峠1川：川川小丫出比|”叫“ 卄小|川呷严叫“叫共轨梯度法是一个典型的共轨方向 法，它的每一个搜索方向是互相共轨的, 而这些搜索方向d仅仅是负梯度方向与 上一次迭代的搜索方向的组合，因此，能量极小化算法Newton-Raphson法利用函数的一次及二次导数求最小值的方法中，最简单的是牛顿-拉森法。一元多次函数U（x）对丑点的二次泰勒展开式为tU（R）=U（工k） +（H Hk）U（Nk） +（H 工女）2/&k）/2 +（3-23） 将式（3-23）对工微分得Ux）的一次导数为：LT（H）=ir（九）+ （ 乂一打）b（

26、竝）（3-24）在极小值点，x=x*函数的一次微分为零，厅（）=0。将此关系代 入式（3-24），得土17（工）=0 = /（顶）+ （工XV ）LTw（Xk） 一职-齡35）若函数为多元函数.则可得类似的关系式，即：新=赤一0（九）一1（佥）(3-26)式中，沪（冠）t为Hessian矩阵的反矩阵以函数/(x,y) = x +2y为例,其Hessian矩阵为:ay3 v3.t(3-27)其反矩阵为:1/201/4(3-28)故依式(3-26)t若由(9, 0,9.0)开始寻找函数f9y)的最小值点，得:9.09.0、4/20厂(0 1/4(3-29)此点即此函数的最小值点。牛顿拉森法利用一次

27、及二次导数，起始 位置离最小值点愈近则效果愈好。此方法因为需要求矩阵的反矩阵，因此 需要较多的储存空间。此外，He韵加矩阵的本征值必须全为正数，否则会 导致移到能量较高点的情形发生。牛顿拉森法因其计算量较大，一般适用 于较小的分子，通常应用于求约含:能量极小化算法比较最陡下降法：方向支化丸，收欽优化辐度丸共辘梯度法收飲快，易陷入局部势阱，对初始结构偏寓不大 Newton-Raphson法计算量较丸，当撤南小肘收欽快构象A. B的焙差III分子力学热力学性质构象能分子在任一构彖的势能称为分子的构象能。対多原子分 子而言，势能为分于键长、键角.扭转角等内坐标的函 数。不同构象所对应的势能则构成分子

28、的多维势能面。点!陰如潜何 种构瓠 对“势依照简谐振子的假设为：3N-6jUvib = S 伽（豆 + MiT 1 t 1f H(AB)=ASE(A)+Zkl4/AB)+ALUAB)+ZLU(AB 一般温度状况下，分子异构物的平均移N转动能相同;分子转动能分子振动能异构物间构象能之差分子移动能AH(A-B)=ASE(A-B)+ AUvib(A-B)分子生成焙分子的生成熔(enthalpy of formarion)为分子的熔与其组成原子的熔的差值。依照热力学,气态分子的标准生成沁98为例如在压力为lbar (lbar=105Pa),温度为298K条件下的反应:xC(石墨，s)+/Hg(g)亠C丹出加(g)其物种数目的变化