ZYM+20111224-分子动力学.ppt

1、分子动力学Molecular Dynamics,甄云梅 20111215,内容,分子动力学 1.简介 2.牛顿运动方程的数值解法 3.积分步长的选取 4.简化单位 5.分子动力计算流程 6.边界条件 7.系综,分子动力学-简介,什么是分子动力学？ 分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。,分子动力学-简介,基本原理,摘自计算化学从理论到分子模拟陈敏伯,基本原理,摘自计算化学从理论到分子模拟陈敏伯,Verlet法（最常见） 如果

2、将粒子的位置按泰勒级数展开： 如果将式子中的t换为-t，则上式变为： 将上面两式相加得到：,分子动力学-牛顿运动方程的数值解法,分子动力学-牛顿运动方程式的数值解法,Verlet法,如果两式相减，则可得到： 即可以由t+t和t-t的位置得到t时刻的速度. 缺点：此式含有1/t,因为t很小，因此非常容易导致误差，为了矫正这个缺点，Verlet又发展出了Leap frog法（跳蛙法）。,分子动力学-牛顿运动方程的数值解法,需要连续记录两个时刻得位置。,Verlet法,Leap frog法,Leap frog法,分子动力学-牛顿运动方程的数值解法,对两式分别求导得到：,分子动力学-牛顿运动方程的数值

3、解法,需要知道上一个时刻得位置ri(t)，中间时刻的速度vi(t-(1/2)h) 和力Fi(t),Leap frog法,预测矫正法Gear Predictor-Corrector,预测矫正法,分子动力学-牛顿运动方程的数值解法,上式,预测矫正法,分子动力学-牛顿运动方程的数值解法,分子动力学-牛顿运动方程的数值解法,分子动力学-积分步长的选取,分子动力学计算的基本思想是赋予分子体系初始运动状态之后利用分子的自然运动在相空间中抽取样本进行统计计算，积分步长就是抽样的间隔。 分子动力学中，最重要的工作为如何选取合适的积分步长，在节省时间的同时也保证计算的精确性。 原则: 积分步长小于系统中最快运动

4、周期的十分之一。 太长的步长会造成分子间的激烈碰撞，体系数据溢出；太短的步长会降低模拟过程搜索相空间的能力。,分子动力学-积分步长的选取,例1,例2,分子动力学-积分步长的选取,分子动力学-积分步长的选取,分子动力学-简化单位,研究分子或原子系统时，如果采用国际单位制，原子质量以g为单位，则通常的原子质量约为10-22g级别；若位置以cm为单位，则通常的量纲为10-8cm；同样积分步长用s做单位通常在10-1310-16s。这些量纲非常小，实验中很容易引起误差，因此实际计算时通常采用简化单位。,分子动力学-简化单位,和为LJ函数的参数值。 原子的质量单位简化为原子的质量数，长度单位简化为。,例

5、如：氩（39.5）的熔点为84K，=0.24kcal/mol,=3.504，系统的密度为1.784g/cm3,设计算时的步长10-12s，压力1atm，则各物理量可以简化为：,分子动力学-简化单位,如果系统中含有不止一种原子，则选取其中一种原子的LJ参数作为简化单位的标准，例如水分子中有H和O两种原子，可以选取O的LJ参数进行简化。,分子动力学-简化单位,分子动力学-分子动力计算流程,系统容量 在进行分子动力学计算前，必须先估计运算的可行性，一般分子动力学计算的容量为数千个原子级别，如果系统过大则超出计算的范围。 对于N个原子的系统，每计算一步需要计算1/2N(N+1)组远程作用力，这是运算过

6、程中最耗时的部分，若原子增加一倍，计算时间则为1/22N(2N+1)为原来的四倍多。,分子动力学-分子动力计算流程,系统容量,分子动力学-分子动力计算流程,研究对象所处的时间范围,通常分子动力研究所选取的积分步长为飞秒fs（1fs=10-15s），若以目前一般的个人工作站或较强的个人电脑从事1000个原子系统的计算，累积100万步即研究10-9s（1ns）的 时间范围，需要两星期的时间?。 因此从实际的角度来讲，分子动力学适合研究反应或运动时间小于1ns的体系，而不适合较慢的反应或运动。例如蛋白质折叠在10-3s（1ms）级别,则需要非常长的时间。,研究对象所处的时间范围,执行分子动力学计算时

7、，将一定数目的分子放在一定形状的盒子中，并使它的密度和实验密度相符合，再选定实验的温度，即可以着手计算。 通常是将分子随机的放置在盒子中，或是按照其结晶位置排列，这将作为分子动力模拟的初始位置。 计算开始时，还必须知道原子的初始速度。系统中所有原子速度产生的动能总和应当满足玻尔兹曼分布：,计算过程,计算过程,分子动力学-分子动力计算流程,原子数目,热力学温度,玻尔兹曼常数,原子速度可以按此式随机产生，只要总动能符合3/2NTkB即可。,由初始位置和速度开始，计算每一步新的速度和坐标，再由新的速度计算系统的温度： 如果得到的速度过大或过小，则会使得到的温度大于或小于设定的温度很多，这时需要校正速

8、度，一般所允许的温度范围为： 如果计算的温度超过允许范围，则将所有原子的速度乘以一个校正因子： 使系统温度重新调整为设定温度：,计算过程,分子动力学-分子动力计算流程,在到达热平衡之前的轨迹信息和速度不需要保存，因为那时的物理意义不严格，真正对我们有意义的是达到平衡后的系统信息。 通常的分子动力学模拟需要累积数百万步以供分析，因而会产生极大的存储。以1000个原子为例，储存每步的位置和速度约要10kB空间，则10000步需要100MB。 因为分子动力学计算的步长很短，每一步移动的距离也很小，通常每隔1020步存储一次来节省硬盘空间。,计算过程,分子动力学-分子动力计算流程,计算过程,分子动力学

9、-分子动力计算流程,分子动力学-分子动力计算的初始设定,执行分子动力学计算时，必须选择恰当的初始条件，如恰当的初始温度、初始速度、初始位置和积分步长等，如果初始条件选择不当则会浪费很长的时间才能达到热平衡，或者说根本达不到平衡。 初始结构 初始结构的选取以越接近模拟体系的理想结构越好。如要模拟体心立方的固体系统时选取了面心立方作为初始结构，便是不明智之举。 通常在模拟前作一些实验，降低系统的局部高能区，使其接近于能量最低点的理想结构，避免计算中由于偏离理想结构太多而导致失败。,分子动力学-分子动力计算的初始设定,初始结构,对于均匀相的液体，通常选用它的固相晶体结构作为初始模型，如果有XRD实验

10、得到的结构数据，则可以直接采用；如果没有，通常采用均匀分布的面心立方结构排列这些分子或原子。 对于复杂液体，通常参考XRD或NMR实验得到的结构位置，或者采用分子力学能量最小化的结果作为起点。,对于小分子可以随意指定其方向，如水和二氧化 碳分子等；,对于大分子则考虑取向时以不发生重叠为宜。,对于分子系统，通常将分子的质心放在面心立方 的格点上，分子的取向视系统而定。,分子动力学-分子动力计算的初始设定,初始结构,通常系统初始速度的产生，是先在-1+1之间选取一个呈高斯分布的随机数，再将其乘以系统所有粒子的平均速度： 即可以得到符合麦克斯韦-玻尔兹曼分布的粒子速度： 通常在实验前还要检查各方向动

11、量是否为0，避免计算过程中系统移动而导致能量不稳定。,初始速度,初始速度,分子动力学-分子动力计算的初始设定,分子动力学-边界条件,分子动力学-边界条件,周期性边界条件,周期性边界条件使人们能够基于少量粒子的模拟的同时考虑周围溶液和溶质分子对研究体系的影响。 把研究的体系看作一个特定形状空间包围的区域，如果采用了周期边界条件，则这个基本单元会沿着所有的方向进行周期性扩展以形成一个无限的周期排列。如图,周期性边界条件,摘自计算机药物分子设计许筱杰,如图中的每个单元周围有8个临接映射单元；如果是三维周期，则有26个。 映射单元中粒子的坐标通过基本单元中粒子坐标的平移操作得到。当基本单元中的粒子离开

12、这个单元进入一侧的映射单元时，其映射粒子会从基本单元的另一侧进入基本单元。,分子动力学-边界条件,周期性边界条件,摘自计算机药物分子设计许筱杰,分子动力学-边界条件,周期性边界条件,摘自计算机药物分子设计许筱杰,非周期性边界条件,分子动力学-边界条件,非周期性边界条件,当采用非周期性边界条件考虑溶剂效应时，可以在研究的重要部位加入非周期的水分子环境，比如球形水环境以及层状水环境。 在实际的分子力学或分子动力学模拟中，为了加快计算的速度，可以把体系划分为柔性区、约束区和完全刚性区。 可以在活性位点的重心为中心分别定义半径为R1和R2的两个球形范围。在模拟中，半径为R1的球体 内的受体、配体和水分

13、子是完全柔性的； 而半径为R2的球体内的受体和水分子则 是加上谐振子回复力约束的；半径为R2 的球体完全是刚性的。 基于这样的策略，在有效考察体系 重要部位的同时可以大大减小计算的消耗。,摘自计算机药物分子设计许筱杰,分子动力学-系综,分子动力学-系综,平衡分子动力学模拟，总是在一定的系综下进行的。,研究扩散系数,研究相变,研究溶液特性,温度调控机制可以使系统的温度维持在给定值，也可以根据外界环境的温度使系统温度发生涨落。 一个合理的温控机制能够产生正确的统计系综，即调温后各粒子位形发生的概率可以满足统计力学法则。,系综的控温,系综的控温,分子动力学-系综,Berendsen温控机制(Bere

14、ndsen外部热浴法),系综的控温,分子动力学-系综,基本思想：假设系统和一个恒温的外部热浴耦合在一起，通过热浴吸收和释放能量来调节系统的温度，使之与恒温热浴保持一致。,摘自计算化学从理论到分子模拟陈敏伯,摘自计算化学从理论到分子模拟陈敏伯,系统常压条件 压力诱导下的相变 系统的压力可以通过改变体积来调节。,系综的控压,系综的控压,分子动力学-系综,Berendsen控压机制(Berendsen外部热浴法),系综的控压,分子动力学-系综,摘自计算机药物分子设计许筱杰,分子动力学-计算结果分析,分子动力学-计算结果分析,运动轨迹分析,热力学特性计算,径向分布函数计算,与时间有关的物理量计算,自由

15、能的计算,直接动画显示 随时间的变化曲线 统计分布曲线,分子动力学-计算结果分析,运动轨迹分析,运动轨迹分析,系统的温度 可以根据系统的平均动能得到系统的平均温度。,热力学特性计算,分子动力学-计算结果分析,热力学特性计算,动能,温度,系统的自由度,动能的统计涨落 系统动能的统计涨落(fluctuation)定义为： 通过统计涨落可以计算出系统的定容比热容（heat capacity at constant volume）Cv：,分子动力学-计算结果分析,热力学特性计算,原子总数,系统动能,玻尔兹曼常数,分子动力学-计算结果分析,径向分布函数计算,系统中参考分子周围rr+dr范围内的分子数目为

16、dN，则定义径向分布函数g(r)为：,为系统的密度,若系统的总分子数为N，可得：,所以径向分布函数可以解释为系统的区域密度和平均密度之比。,因为：,参考分子附近的区域密度不同于系统密度，但与参考分子距离远 时区域密度应与平均密度相同，即当r值大时径向分布函数的值应当接近于1.,径向分布函数计算,分子动力计算径向分布函数的方法为：,分子动力学-计算结果分析,径向分布函数计算,与时间有关的物理量计算,分子动力学-计算结果分析,与时间有关的物理量计算,分子动力计算除了计算系统的平均值以外，最重要的是计算系统的各种动态特性。,均方位移,分子动力计算当中，分子由起始位置不停的移动，每一瞬间的位置都不相同

17、，如果以 表示时间t时粒子i的位置，则粒子位移平方的平均值称为均方位移（mean square displacement，MSD），即：,根据统计原理，只要分子数目足够多，计算的时间足够长，系统的任何一个瞬间都可以作为起点，计算的平均值相等，因此计算MSD时常采用多时间起点：,设一共运行n步，各步位移 ,将此轨迹分为相等数目的两部分: 计算MSD时每次均取n/2组数据的平均:,分子动力学-计算结果分析,与时间有关的物理量计算,根据爱因斯坦扩散定律，有：,式中D为粒子的扩散系数，因此当模拟的时间很长时，均方位移对时间曲线的斜率即为6D，由此可求出粒子的扩散速率。,分子动力学-计算结果分析,与时间

18、有关的物理量计算,分子动力学-计算结果分析,与时间有关的物理量计算,相关函数,设A，B为系统两个不同的物理量（位置、速度等），则A或B的自身相关函数(auto correlation function)可表示为：,分子动力学-计算结果分析,与时间有关的物理量计算,分子动力学-计算结果分析,与时间有关的物理量计算,选取的最大轨迹桢差。,对于相关函数的计算，通常有两种方法：,1.直接计算法 缺点：占用的存储空间大； 2.快速傅里叶变换法: 复杂，但更有效,相关函数可以算什么？,分子动力学-计算结果分析,与时间有关的物理量计算,速度自相关函数可以计算粒子的扩散系数，关系为,对C和H的VACF积分，得

19、到：,表明在系统中H原子运动的速度要快于C原子。,1,2,与原子运动相关的频率可以由速度相关函数的傅里叶变换得到：,为频率 Re表示取括号中的实数部分,频率的图谱密度,根据轨迹将上式转化为：,分子动力学-计算结果分析,与时间有关的物理量计算,3,分子偶极的相关函数可以用来计算红外光谱，光谱的线状函数为:,红外光谱仪检测到的频率位置和光谱计算得到光谱位置应该相同，因此可以对比验证所用力场参数的正确性。,可以计算拉曼散射光谱，拉曼光谱的函数为：,4,式中， 表示 归一化的转动相关函数（rotational correlation function）,分子动力学-计算结果分析,与时间有关的物理量计算,5,分子动力学-计算结果分析,与时间有关的物理量计算,自由能的计算,分子动力学-计算结果分析,自由能的计算,根据统计力学，分子体系的自由能可由下式计算，,自由能微扰（FEP）方法,对于生物大分子或溶液中的分子，由于其柔性较大，直接计算体系的自由能较为困难，常用自由能微扰方法计算两状态的自由能差来