分子动力学模拟.ppt

《分子动力学模拟.ppt》由会员分享，可在线阅读，更多相关《分子动力学模拟.ppt（63页珍藏版）》请在三一办公上搜索。

1、中国石油大学,分子动力学模拟,一、系综理论,本章主要内容,二、分子动力学方法,分子动力学模拟,三、模拟细节,四、参量的计算,五、液态水的MD模拟,六、误差分析,七、分子动力学模拟方法的应用,分子动力学模拟（molecular dynamics simulation，简称MD）方法首先是由Alder和Wainwright提出的，现已逐渐成为预测系统特性、验证理论和改进模型的计算工具。MD方法的基本思想是把物质看成由原子和分子组成的粒子系统，从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解出运动方程而得到系统中全体粒子在相

2、空间中的轨道，然后统计得到系统的热力学参数、结构和输运特性等。也就是由体系的微观性质来求算其宏观性质。属于微观尺度的模拟技术。,一、系综理论,分子动力学模拟,相空间或称相宇，是经典的系综理论中为了描述系统的微观状态而引入的一个基本概念。,1.1 Phase space,空间中的每一点，代表系统一种可能的运动状态，系统的不同运动状态，可以用相空间中不同的点来代表。从经典力学可知，系统的运动状态随时间变化，因此相空间中代表系统运动状态的点将在相空间中移动，形成轨道。,分子动力学模拟,处在一定宏观条件下，大量性质和结构完全相同的系统的集合。又称统计系综。系综是采用统计方法描述系统的统计规律性时引入的

3、一个最基本的概念。,1.2 系综(Ensemble),相空间中的不同点代表了系统的不同运动状态，也可以把它们看作是性质和结构相同而初始条件不同的系统所处的状态。这些系统的集合就是系综。统计系综的各个系统反映了实际系统在不同时刻的面貌。,分子动力学模拟,1.3 准各态历经假说,一个力学体系在长时间的运动中，它的代表点可以无限接近相空间中的任何点，系综平均等于长时间的时间平均。,系综的方法可用于由彼此存在相互作用的大数量粒子所组成的系统，按实际系统所处的不同宏观条件，应当采用不同的系综（微正则系综、正则系综、巨正则系综）及其对应的不同形式的分布函数。,统计方法的特点是由微观量求宏观量。,分子动力学

4、模拟,1.4 常用到的系综,正则系综（NVT）：一个粒子数为N,体积为V,温度为T和总动量为守恒量的系综,在这个系综中系统的粒子数(N),体积(V)和温度(T)都保持不变,并且总动量为零。,微正则系综（NVE）：它是孤立的、保守的系统的统计系综。在这种系综中，体系与外界不交换能量，体系的粒子数守恒，体系的体积也不发生变化。,巨正则系综：由与温度恒定的大热源和化学势恒定的大粒子源接触，具有确定体积的系统构成的统计系综。,等温等压系综（NPT）,分子动力学模拟,1.5 统计系综（Statistical ensembles）,、(T,P N)constant pressure（等温等压系综）,、(E

5、,V,N)microcanonical（微正则系综）,、(T,V,N)canonical,constant volume（正则系综）,、(T,V,)grand canonical,分子动力学模拟,2.1 Newtonian mechanics,In the MD method,every new distribution is derived from the previous one by using the interactions between the particles.The interactions depends on the position of the particles

6、.,In that potential the particle feels a force,According to Newtons second law,二、分子动力学方法,分子动力学模拟,2.2 Potential energy functions,对一个由N个原子构成的简单系统，其势能项由下式给出,式中右端第一项是外场（如电场、磁场、声场等）对系统的作用；第二项是两体势即系统中每两个粒子间的相互作用；第三项是三体势，表示系统中每三个粒子间的相互作用,有效两体势,它包含多体效应，可很好地反映系统粒子间的相互作用。,分子动力学模拟,简单粒子系统,采用Lennard-Joans势、硬球势、软

7、球势、方阱势等模型,复杂粒子系统,采用多中心的位置-位置（site-site）相互作用模型，但其中心间的相互作用仍然采用简单的位能模型。,带电分子,库仑（Coulomb）相互作用。,处于外场中,静电相互作用是长程相互作用,分子动力学模拟,1、Lennard-Joans势,下面仅对简单系统的相互作用模型给予简介,分子动力学模拟,2.硬球势（Hard-Sphere）,3.软球势（Soft-Sphere）,通常，v 是为整数的参数。,4.方阱势（Square-Well）,分子动力学模拟,2.3 Calculations of force,velocity,position,The initial d

8、istribution of the Molecular dynamics simulation is generated in a random distribution.,In simulation:,Each particle is also assigned an initial velocity vi,分子动力学模拟,The force causes an acceleration,Which in turn modifies the initial velocity vi as,And modifies the initial position ri as,The MD simul

9、ation can describe systems that evolve in time.The new positions are derived from the Newtonian law of motions and therefore deterministic.,分子动力学模拟,2.4 Equations of motion,为了在计算机上解运动方程，必须为微分方程建立一个有限差分格式，从差分方程中再导出位置和速度的递推关系式。这些算法是一步一步执行的，先算t 时刻的位置和速度，然后在此基础上计算t+1时刻的位置和速度。,微分方程最为直接的离散化格式来自泰勒展开:,1、有限差分

10、方法-预测校正法,分子动力学模拟,2、有限差分方法-Verlet算法,、Verlet算法的一般形式,为了用数值方法求解微分方程，,采用有限差分方法离散化得到差分格式,两式相加得,分子动力学模拟,上式提供了一个方法，从粒子在前两步（t和t-h）时刻的位置以及t时刻的作用力来得到粒子在t+h时刻的位置。,速度对轨道计算没有关系，但对估算动能及速度相关函数（用来研究粒子的输运特性）非常有用。速度形式为,分子动力学模拟,、Verlet 蛙跳格式（leap-frog）,The velocity varies so one has to choose a reasonable average value

11、to be used.The velocity at the middle of the step ought to be a good compromise,分子动力学模拟,、Verlet速度相加形式,分子动力学模拟,分子动力学模拟,3.1 初始化位型,模拟开始时首先要初始化位型分布，既首先要给定微元中分子的初始位置和初始速度。,、简单立方晶格,、体心立方晶格,、面心立方晶格,三、模拟细节,分子动力学模拟,分子动力学模拟,初始速度,模拟时，各粒子的初始速度按麦克斯韦速度分布取样。,Maxwells distribution law of velocity,指平衡状态下理想气体分子速度分布的统

12、计规律。在平衡状态下，分布在任一速率区间内的分子数与总分子数的比率为,分子动力学模拟,3.2 Periodic boundary conditions,通常选取小数体系（几十个到数千个分子）作为研究对象，但由于位于表面和内部的分子受力差别较大，将会产生表面效应。,为了消除此效应，同时建造出一个准无穷大体积，使其可以更精确地代表宏观系统，必须引入周期性边界条件。,把小数体系看作一个中心元胞，此中心元胞被与中心元胞相同的其它元胞所包围，当某分子离开中心元胞时，该分子将在整个点格上移动以致它从中心元胞对面重新进入这个元胞。,周期性边界条件：,分子动力学模拟,分子动力学模拟,3.3 Calculati

13、on of interactions,i,j 两个离子的相互作用势,体系总势能可表示为：,在计算势能时必须考虑相互作用的力程，对粒子数为N的模拟系统，原则上任何两个粒子间都有相互作用，在计算体系势能时须进行N（N-1）/2次运算。,为了提高计算效率，实际模拟过程中通常在一个方便的力程上将位势截断，用以减少计算势能所消耗的时间。,分子动力学模拟,对势能的最大贡献来自于粒子的近邻区域，位势截断常用的方法是球形截断，截断半径一般取2.5或3.6，对截断距离之外分子间相互作用能按平均密度近似的方法进行校正。,分子动力学模拟,Cell method,Neighbor list method,1、近邻表-

14、Verlet近邻表,2、近邻表-网格近邻表,分子动力学模拟,3.4 标度与趋衡,正则系综的MD模拟，模拟过程中必须保证体系温度恒定，所以在每一时间步要对分子的速度进行标度使其满足正则系综的条件，标度因子为，其中T是指定的体系温度，T*是每一时间步的瞬时动能温度。,标度：,分子动力学模拟,一种按上述方法建立的系统不会具有所要的能量，而且很可能这个状态不对应一个平衡态，而所期望的体系参量是在平衡态通过统计平均得到的，所以为使系统达到平衡需要一个趋衡阶段。,趋衡：,对于正则系综，体系趋衡是通过增减体系能量来实现的。如果在进行足够多的运算步数以后，系统能持续给出确定的平均动能和平均势能的数值，那么就认

15、为平衡已经建立。,分子动力学模拟,3.5 The MD units,这里我们以能量（焦耳，J）和粒子间距离（米，m）为基本单位来对其它量进行约化，使其成为无量纲量。仅举例如下：粒子数密度：温度：能量：压强：时间：力：,分子动力学模拟,A simulation run produces the raw data in form of a very large disk file,The file contains the complete state of the system at each step.,The disk file is processed after the simulati

16、on is finished.It contains at least all the positions and velocities of all particles.This information is sufficient to calculate all the properties of the system.However,it is more economical to calculate properties during the simulation and store them in the file rather than reading the file and c

17、alculating them afterwards.,The potentials and forces certainly belong to this class because they are needed at each time step anyway in order to calculate the new position and velocities of the particles.,四、参量的计算,分子动力学模拟,1内能,模拟体系的内能由两部分组成,2比热,涨落是体系中热力学量在不同时刻出现的波动。系统比热可以通过能量的涨落导出，表示如下,分子动力学模拟,3.压强,体

18、系的压强可通过维里定理导出：,分子动力学模拟,4扩散系数,在模拟过程中，体系的迁移系数可以通过平衡相关函数导出。由Green-Kubo方法，计算扩散系数的公式如下：,由于受计算机硬件的限制，在模拟过程中不可能对时间积分至无穷。为了选择合适的积分上限，在模拟过程中应首先计算各水分子的速度自相关函数，从水体系归一化的速度自相关函数可以看出在什么长度的时间间隔时，速度自相关函数逐渐趋零，也就是说在什么时间范围内速度基本上不再相关，那么求扩散系数所选取的积分区间就可以确定。,为速度自相关函数（VACF）,分子动力学模拟,5、粘滞系数,利用平衡分子动力学模拟方法计算平衡系统的粘滞系数的基础是Green-

19、Kobo理论。根据该理论，粘滞系数由下式给定,分子动力学模拟,6、径向分布函数-g(r),径向分布函数是与时间无关的关联性的量度，精确地说就是与中心粒子间距为rr+dr 的体积元内分子的平均数密度。,在实际模拟过程中我们以元胞边长的1/2为半径作圆球，将其分割为厚度相等的100个球壳，计算分子间距R满足（ij，且i,j 依次取值）的分子对个数N(i)，i=1,2,.100.,分子动力学模拟,分子动力学模拟方法的基本步骤,1.确定研究对象，如确定所研究的体系是属于微正则系综、正则系综、巨正则系综呢还是属于等温等压系综等等；2.建立适当的位能模型，如假定所研究的体系位势形式满足硬球势、软球势或勒纳

20、德琼斯势等等；3.初始化位形和速度；4.建立体系粒子的运动方程并进行求解，求解时需对位势进行截断，并应用周期性边界条件和最小影像数约定，每一时间步对速度进行标度以使体系达到趋衡；5.通过统计平均来求解一些我们感兴趣的热力学参量等。,五、液态水的MD模拟,分子动力学模拟,1模拟参数,本文将水系统视为正则系综进行MD模拟，模拟所用分子数为256，水分子初始位型按面心立方晶格分布，晶格所占空间体积由水的密度导出。各分子起始速度按Maxwell分布取样。模拟温度分别取=273K、293K、313K、333K、353K、373K。模拟过程中为保证体系温度恒定，在每一时间步对分子的速度进行标度；计算体系位

21、能时采用球形截断法，截断半径为2.5，对截断距离之外分子间相互作用能按平均密度近似的方法进行校正；运动方程的求解采用Velocity Verlet Algorithm算法；时间步长为210-15s，每次模拟的总步数为105步，其中前50000步使体系趋衡，后50000步用来统计各种平衡性质和结构性质；模拟过程中各参量均采用约化单位。,分子动力学模拟,计算体系的内能有着重要的意义，不仅是可以通过能量求出比热、表面张力等热力学量，更为重要的是内能（或势能）的演变是我们在模拟过程中判断体系平衡的依据。体系是否达到平衡直接关系到热力学量的真实性、稳定性和准确性。,2水的内能和比热,分子动力学模拟,表1

22、 水的势能模拟结果与文献值和实验值的比较,表2 水的比热的模拟结果和实验结果的比较,分子动力学模拟,3水的扩散系数和粘滞系数,表3 水的扩散系数模拟结果与文献值和实验值的比较,分子动力学模拟,4水的径向分布函数,径向分布函数是对物质结构的反映，在研究物质特别是液态时有着重要的作用。,分子动力学模拟,1、位能模型,位能模型是对体系中分子或离子之间相互作用势的反映，位能模型选择的正确与否决定了模拟的成败。但是许多情况下由于研究对象的复杂性，使对问题的认识还不够深入，导致所选择的位能模型和实际情况存在一定差距；同时受计算机硬件的限制，也不可能选择更接近实际的复杂模型，如多体相互作用模型，从而导致误差

23、出现。,表1 关于纯水结构的理论模型,六、误差分析,分子动力学模拟,2、模拟分子数,分子动力学模拟的主要任务就是要从微观来求算体系的宏观性质，所以模拟过程中参与模拟的分子数越多，模拟结果就越真实，越接近体系的宏观属性。例如Tanaka等人采用1.36104个水分子对水体系进行了模拟，其扩散系数更加接近实验值。但也应考虑到：模拟过程中由于分子受力和势能的计算需要涉及微元内所有的分子对，所以计算量很大，非常消耗机时，一般占总计算量的80%左右。模拟体系分子数扩大一倍，总模拟时间并不仅增加一倍，而是呈级数形式上升。,要想解决大数粒子体系的模拟问题，还须引入一些模拟技巧，如Verlet Neighbo

24、ur List和Link List等方法。,分子动力学模拟,表2 不同方法每时间步所消耗CPU时间的比较,分子动力学模拟,3、求解运动方程的算法,求解分子运动方程的方法有两种，一种是Verlet Algorithm，它是对二阶微分运动方程直接求解的一种方法；另一种则是Gear预测校正算法（Gear Predictor-Corrector Algorithm），它是通过泰勒级数展开的方法实现的。,在时间步长较小的情况下，通过高阶Gear算法得到的结果更为精确，但是模拟者对大步长更感兴趣，因为相同的总模拟步数，单步长越大，模拟体系演化的时间越长，就越接近实际情况，模拟结果就越真实，在这种情况下Ve

25、rlet Algorithm算法则更具吸引力。,分子动力学模拟,4、总时间步、时间步长,一般来说模拟总时间步越多，单步长越大，模拟结果就接近实际情况，但是由于受计算机硬件的限制，不可能对一个足够长的时间过程进行模拟，同时单步长较大，对体系演变反映的跨度超长，计算误差也就越大；如果单时间步长选择太短，小于体系分子位置的演变，也会造成误差的增加，所以模拟时要结合所选取的求解分子运动方程的算法，合理选择步长以减小误差。,求解中的另一个重要问题是时间步长的选择：步长h连同实现的MD步数，决定了相空间有多大一部分被抽样，看起来越大越好，这样可以对更大的部分抽样，但是由h所决定的时间标尺和实际体系发生变化

26、的时间标尺并不一定完全相同，所以时间步长的选取应根据研究对象的具体情况而定,分子动力学模拟,5、位能截断,为提高计算效率，实际模拟过程应进行位势截断，对截断距离之外分子间相互作用能按平均密度近似的方法进行校正，虽然进行了校正，但难以避免误差的出现，特别是在分子间存在长程相互作用时，进行位势截断会导致严重的失真，在这种情况下则必须采用Ewald求和的方法。,上述的误差分析可以看出，误差产生的原因具有多样性。为了减少误差、提高精度，单纯从某一个方面入手是盲目和片面的，在误差产生的过程中它们是相互影响、相互制约的。所以实际模拟过程中，必须认真分析个方面的情况，在精度和误差之间作出平衡，统筹兼顾，以达

27、到最佳的模拟效果。,分子动力学模拟,Adsorption of methane on zeolite,分子筛是具有均匀的微孔、其孔径与一般分子大小相当为一类吸附剂或薄膜类物质。在化学化工中经常被用作催化剂或者催化剂的载体。,七、分子动力学模拟方法的应用,分子动力学模拟,Layer-Cell Modeler easily creates simulation of water-benzene interface,在体系内部物理性质和化学性质完全均一的一部分称为“相”（Phase)。相与相之间在指定的条件下有明显的界面，在界面上，从宏观的角度看，性质的改变是飞跃式的。,分子动力学模拟,MgO Su

28、bstrate Sputtered by Ar Plasma,分子动力学模拟,Aqueous LiBr subject to external electric field,分子动力学模拟,分子动力学模拟,石墨烯/硅纳米线复合材料自卷曲过程,平行放置,相互靠近,完成第一层包覆,剩余石墨烯完成包覆,形成稳定结构,分子动力学模拟,实验测定与 MD 模拟结果的比较W.Linert,and F.Renz,J.Chem.Inf.Comput.Sci.,33,776(1993),Experimentally Determined,MD-predicted,(1,2,4)-4-(1,1-e二甲乙基)-2-烃

29、基环戊羰基酰胺晶体,分子动力学模拟,MD预测的顺磁性和反磁性冰晶体结构,Ferromagnetic ice,Antiferromagnetic ice,O.A Karim&A.D.J.Haymet,J.Chem.Phys.,89,6889(1988),分子动力学模拟,Ru-Al 合金断裂过程动态模拟,C.S.Becquart,D.Kim,J.A,Rifkin,and P.C.Clapp,Mat.Sci.Engin.,A170,87(1993),断裂点周围的损坏区域,分子动力学模拟,Molecular dynamics simulations of the adsorption of industrial silane molecules at a zinc oxide surface,分子动力学模拟,分子动力学模拟,液相条件下单分子在金属表面的吸附,水溶液中A(7)、C(11)、E(15)和H(21)分子在Fe 和FeCO3表面的平衡吸附构型,Fe,FeCO3,氨基上的N原子易被极性水分子极化成为(NH3)+,分子动力学模拟,结束,分子动力学模拟,