第六章 分子动力学模拟 Molecular Dynamics –MD 6.1引言分子动力学模拟方法是在牛顿力学的理论框架下，根据体系内分子之间的相互作用势，获得每个原子随时间运动的轨迹，通过系综平均，可以得到感兴趣的与结构和动力学性质有关的物理量，如：平均原子坐标，平均能量、平均温度及原子运动的自相关函数等。

这些物理量是通过对每个原子的运动轨迹，即微观量求平均而得到的宏观量，因此可以与实验观测量进行比较。

用计算机模拟方法在向空间采样方法有两种： （1） 随机采样 MC （2） 确定性方法MD以上讲过的MC （Monte Carlo ）采样方法就是随机方法，与随机方法不同，确定性方法是按照动力学规律使系统在相空间运动。

分子动力学模型就是一种确定性方法。

它的基本出发点是从一个完全确定的物理模型出发，通过解牛顿运动方程而得到原子运动的轨迹。

我们感兴趣的可测量的客观物理量可以通过相空间的采样求系综平均而得到。

在多态历经假设成立的情况下，系综平均与长时间平均是相同的。

⎰∞→∞==τττ01))(),((limdt t p t q A A A系综其中q,p 为t 的函数。

A 表示系综平均，∞A 表示无穷长时间平均。

因模拟时间总是有限的。

对耦分子体系，当模拟时间大于分子的弛豫时间时，有限观测时间可以变成为无穷长的。

当弛豫模拟〉τt ，模拟t 可认为∞，因物理上的∞是不可能的。

6.2基本原理 1.动力学方程基本动力学方程包括在经典力学（CM ）框架下的牛顿方程和在量子动力学（QM ）框架下的薛定谔方程。

在常温下，经典的牛顿方程对研究生物分子体系的结构和动力学性质已经足够了，因为这时体系的量子效应并不十分重要。

但是，对研究包含隧道效应的反应时间问题时，量子效应十分明显，这时就必须用QM 方程来模拟体系的量子动力学性质。

QM:含时薛定谔方程为),(),(t r i t r H t→∂∂→∧-=ψψ （2.1）其中∧H 为哈密顿算符，),(t r →ψ为波函数，→r 表示一系列原子坐标，即),,(21→→→→=N r r r r 。

（关于量子动力学模拟的有关内容下面有专门章节去讲，这里就不赘述了） CM:经典的牛顿运动方程不同坐标系下有不同的表达形式。

在直角坐标系下，牛顿方程可写为)(22t F m i dt r d ii →=→（2.2）其中i m 为原子质量，力→i F 通过对势函数求一阶导数而得到，即→→→→∂∂→-=iN r r r r U i F ),,(21 （2.3）其中→i r 为直角坐标，),,(21→→→N r r r U 为N 个粒子系统的相互作用势函数（U为半经验势）。

若用广义坐标i q 来描写经典力学体系，运动方程的形式就变为拉格朗日方程（lagrangian Equation ）0..=-∂∂∂∂∙iiq L q Ldt d (2.4)其中拉格朗日函数为体系动能与势能之差，即)(),(),(q U q q K q q L -=∙∙(2.5)上式中广义坐标{}M q q q q ,,21≡，M 为体系的自有度数。

拉氏方程常用来研究体系与热浴的耦合以及键角约束等。

因这些问题选用了广义坐标比较方便。

经典力学的另一种表达形式是哈密顿正则方程，ii i i q H P P H q ∂∂-=∂∂=∙∙（2.6）其中系统哈密顿量H 可表为动能和事能之和，H 为系统的特征函数，包含了完整的系统力学行为的全部信息。

)(),(),(q U q p K q p H +=（2.7）其中广义动量i P 为 []∙∙∂-∂=∂∂=iii q q U q p K q L P )(),(,i=1,2,……M （2.8）在实际的MD 模拟中，直角坐标系下的牛顿运动方程的积分具有简单的形式，因为势函数U 只与直角坐标有关。

2.运动方程中的力MD 模拟有一个基本的假设，即原子之间的半经验相互作用势函数可以有足够精确的精度来描写体系原子之间的相互作用，并可应用到象蛋白质分子这样复杂的体系中去。

原子相互作用势可表为直角坐标→r 及一系列力参数S 的函数。

),.....,;,.......,(),(2121M N s s s r r r U s r U →→→→=（2.9）势函数U 可以分解为不同类型相互作用之和，例如可以把相互作用势按相互作用原子的个数进行分解∑∑∑∑∑→>>>→→→→>>→→→>→→→+++++=alli N lk j i l k j i kj i k j i ji j i ii s r Us r r r r Us r r r U s r r U s r U s r U ),(........),,,,(),,,(),(),(),()()4()3(,)2()1(（2.10）下面我们分项讨论以上不同类型的相互作用势及其所对应的相互作用力：01第一项为一体相互作用以体相互作用常用于研究带电离子i q 在电场→E 中的力)(t E q F i i →→=另外，原子的位置约束势pr U 可以写为谐和势的形势∑⎥⎥⎦⎤⎢⎢⎣⎡-=→→→→ii i pr i pr pr r t r K r k r U 200)(21),,(（2.11） 其中Pr iK 为位置约束力常数，0i r →为约束参考位置。

相应的约束力为：))(()(0→→→--=i i pri i r t r K t F （2.12） 在MD 模拟中，原子位置约束可用来建立比较合理的体系。

例如，当蛋白质分子周围放入水分子时，由于一开始溶质与溶剂原子之间可能存在不合理的接触，若用这一初始构型进行EM 或MD 模拟，势必会造成蛋白质分子的变形。

所以，可对溶质原子进行位置约束而让溶剂分子充分运动。

然后再把溶质原子的位置约束去掉与水分子一起作充分MD 模拟。

这样作可以获得比较合理的蛋白质构象。

02第二项为二体相互作用“对”相互作用包括共价键相互作用、非键范德华和库仑相互作用以及NMR 实验中测定的原子之间的距离约束等。

最简单的共价键能为谐振函数[]∑-=→nnn b nb b b t bK b k r U 200)(21),,(（2.13）其中b n K为键能力常数，0nb 为平衡键长。

原子i 和j 之间的键长→=ij n r t b )(，且→→→-=j i ij r r r 。

作用在共价键原子上的力为：⎥⎦⎤⎢⎣⎡--=→→→)()())(()(0t r t r b t b K t F ij ij n n b n i （2.14）)(t F F i j →→-=键能的另一种模型为Morse 势 （2.15）常用的非键相互作用能的最简单形式为 （2.16）相应的作用于原子i 上的非键相互作用力为[)()(20)(6)(12)(4)(t r t r ij r ji t R A t R B i ij ij ij ijij ijt r q q t F →→⎥⎥⎦⎤+-=επε（2.17）对范德华相互作用中排斥项12/ij ij r B 还可以用指数函数描述[]0/)(exp ijij ij R t r B - 其中0ij R 为附加的参数。

相应的排斥力为)()()()exp()(00t r t r R t r R B i ij ij ijij ijijt F →→-=在计算（2.16）式的非键相互作用求和时，原则上要对所有原子对求和。

但在实际应用时，还必须排除近邻原子对的相互作用。

因为第一、二近邻的原子对之间的距离太近，若不排除掉这些非键对，非键能就会变得很大。

对于第三近邻的原子，不同力场有不同的处理方法。

如ECEpp 力场。

就排除掉与第三近邻原子之间的非键作用力。

对AMBER 力场，第三个相邻原子队之间的非键相互作用要根据标准值进行调节（）对GROMOS 力场，第三相邻原子队之间的范德华参数和比标准值要小。

GROMOS 力场，芳香环的第三近邻非键对要排除掉，关于排它原子的定义见左03第三项为三体相互作用例如势函数中的键角弯曲能项包括i,j,k 三个原子，键角能为： （2.18） 其中baK为键角相互作用力参数，0θ为平衡键角。

键角n θ的值可用键向量→ij r 和→kj r 来表示)arccos(→→→→=kjij kj ij r r r r n θ（2.19） 图作用在原子i 上的力为（2.20） 其中 （2.21）同理可以算出作用在原子k 上的力→k F ，则作用在原子j 上的力为))()((t F t F F k i j →→→+-= 图04第四项为四体相互作用正常二面角和非正常二面角由四个原子i,j,k,l 组成，二面角ϕ由原子（i,j,k ）和原子（j,k,l ）组成的两个平面来定义，其能量项为0=n δ或π6,5,4,3,2,1=n m（2.22）其中da n K 和id n K 为力参数，正常二面角n ϕ的值可在（0,2π）范围内取值，而非正常二面角的值被约束在0n ϕ附近，非正常二面角的引入是为了保持分子的手性（chirality ）或者是使侧链的芳香环保持在一个平面内。

（关于正常二面角和非正常二面角的定义见后） 对正常二面角，作用在原子i 上的力为： （2.23）对非正常二面角，作用在原子i 上的力为： （2.24）由二面角的定义，n ϕ 可用键向量表示 （2.25）其中 sign 为符号传送函数。

⎪⎪⎩⎪⎪⎨⎧-=无符号传送a a x a sign ),( 000=<≥a x x05多体相互作用若考虑多体相互作用对非键相互作用的贡献，需要考虑分子的极化，为三极矩、四极矩等。

在实际模拟时，N 体相互作用项可用对相互作用去近似，具体是通过调整相互作用参数而实现的，称为有效相互作用。

例如液相中的极性分子，其偶极矩会受到周围溶剂分子的影响而产生诱导偶极，平均偶极将会增大。

最简单的形式是多项同性点偶极。

偶极矩在i 位点上的诱导偶极矩（induced dipole ）为 （2.26）其中i α为极化率，→i E 为原子处的电场，ij T 为场张量50243ijij j i ij rr r r T πε-=→→ （2.27）方程（2.26）可用迭代法求解对于生物大分子体系，即使仅考虑对相互作用，非键对数正比于2N ，N 为原子数。

显然，原子个数越多，非键相互作用的计量数就会越大。

通常用截断半径（cutoff ）来减少长程非键对数。

对长程静电相互作用的处理方法，下面用专门一章来讨论。

6.3 边界条件计算机模拟只能模拟有限的系统，模拟的系统的粒子数远远小于阿佛伽德罗常数23010*02.6=N 为了减小有限尺度对模拟结果的影响，正确处理边界条件就显得十分必要。

常用的边界条件的处理方法有下面三种： 1. 真空边界条件（Vacuum boundary Condition ）真空边界条件是最简单的边界条件，与气相（gas phase ）类似，即环境为压强0=p 的气体。