分子模拟教程(径向分布函数)ppt课件.ppt

《分子模拟教程(径向分布函数)ppt课件.ppt》由会员分享，可在线阅读，更多相关《分子模拟教程(径向分布函数)ppt课件.ppt（76页珍藏版）》请在三一办公上搜索。

1、第三章 分子模拟方法,1. 蒙特卡罗(Monte Carlo)方法基础 2. 分子动力学(Molecular Dynamics)方法基础 3. 程序讲解,本章具体讲解内容,掌握分子模拟方法的必备知识：,编程技能 (Fortran or C/C+) 统计物理学（统计力学）： 统计物理学基础； 系综原理； 非平衡统计力学基础； 涨落理论分子热力学 ：分子间相互作用理论； 分布函数理论气体分子运动论其它,分子模拟的目的:,为什么要进行分子模拟？,将分子聚集体的性质与如下方面相联系：分子的微观相互作用分子聚集体的结构分子的动力学过程,分子模拟对实验进行补充，使我们能够：预测现有或新材料的性质在分子水平

2、研究宏观现象获得实验无法或难以发现的东西,什么是计算机分子模拟方法？,分子模拟的定义： 统计力学基本原理出发，将一定数量的分子输入计算机内进行分子微观结构的测定和宏观性质的计算。,按照获得微观态的方法不同，分子模拟分为：蒙特卡罗方法 (Monte Carlo, MC)分子动力学方法 (Molecular Dynamics, MD) (3) 混合方法 (hybrid method，HM),计算机分子模拟的发展历史：,1. 蒙特卡罗方法（MC）1953 Metropolis, Ulam, Rosenbluth and Tell， Los Alamos National LabMonte Carlo

3、 simulation of hard sphere.2. 分子动力学方法（MD）1957Alder and Wainwrigth，Livermore LabMolecular dynamics simulation of hard spheres.,微观与宏观,分子模拟在微观尺度与实验室的宏观世界之间起着桥梁的作用：,给定分子间的相互作用 “准确”预测研究体系的性质,MC与MD的区别：,MC: 构型平均，不包含动力学部分； 利用概率行走产生微观态。MD:时间平均，产生动力学性质；利用运动轨线随时间的变化来产生一系列微观态。,计算机分子模拟的发展历史（续）：,从上个世纪九十年代初期以来，计算机

4、模拟技术得到了飞速发展，主要基于三个方面的发展：分子力场的发展（基石） （Amber，OPLS、Compass）原子间的键长、键角、分子间的内聚能等模拟算法（途径）计算机硬件（工具）,HPCx,计算机分子模拟的特点：,原子水平的模拟计算机实验检验理论、筛选实验科学研究中的第三种方法,分子模拟中涉及的几个基本概念：,模拟计算盒子或模拟胞腔,Simulation box (cell),装有一定数目流体分子的研究对象，它是我们要研究的宏观体系的缩微模型。,立方形胞腔,周期性边界条件(Periodic boundary condition, PBC),在小体系中，边界效应总是很显著。在包含1000个原

5、子的简单立方晶体中488个原子处于边界上。在包含1000000个原子的简单立方晶体中仍然有 6%的原子在边界上。,在模拟中，考虑具有真实边界的对象，不切合实际： 增强了有限尺寸效应 人为造成的边界会影响流体的性质,当某个粒子运动出模拟盒子的某一边界时，另外一个影像粒子从另一对立边界进入到此盒子中。,周期性边界条件(Periodic boundary condition, PBC),本体体系的近似：中心盒子在X，Y和Z方向无限扩展；消除人为形成边界的表面效应；保证中心盒子中的粒子数恒定。只需要跟踪中心盒子中各粒子的运动。,周期性边界条件的算法：,采用数学函数：,FLOOR(r/L): 返回不超过

6、r/L的最大整数,FLOOR (4.8) has the value 4.FLOOR (-5.6) has the value -6.,采用数学函数：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),周期性边界条件的算法：,y,最小影像转化原理(Minimum image convention),定义：,中心元胞中的一个粒子只与此元胞中的其它N1个粒子，或它们的最近邻影像发生相互作用。,适用条件：,粒子间相互作用势能的截断距离必须不大于模拟中心元胞长度的一半。,此两粒子与中心粒子的距离相等，但是：黑色球发生作用

7、绿色球不发生作用,此两粒子是与中心原子相互作用的最近邻影像,最小影像转化原理的算法：,采用数学函数：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),截断势能(Truncating the Potential),本体体系采用周期性边界条件描述：不可能将所有粒子与它们影像粒子间的相互作用全都计算。 必须在不大于中心盒子长度的一半处进行截断，以便与最小影像转化原理一致。粒子间的相互作用主要来自于截断范围内，而范围外的贡献很小，可忽略不计。,截断范围内的相互作用,截断势能函数的形式：,简单截断势能函数(Truncat

8、ed Potential)：,缺点：,rc: 截断距离或半径,势能在截断处不连续,当一对分子穿越边界时，总能量不守恒。 分子间力在截断处为无穷大，MD运动过程不稳定。,忽略截断半径之外的所有作用,位移截断势能函数(Shifted and Truncated Potential)：,缺点：,分子间力仍然在截断处不连续。,优点：,势能在截断处连续，但不影响分子间力的大小 分子间力在截断处不为无穷大,截断势能函数的形式：,常用于MC和MD模拟中,位移力截断势能函数(Shifted-Force Potential)：,常用于MD模拟中,优点：,势能和分子间力均在截断处连续,截断势能函数的形式：,截断势

9、能函数的对比：,位移力截断势能,简单截断势能函数,一、Monte Carlo模拟方法基础：,亦称统计模拟或随机抽样方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,投硬币，掷骰子,Monte Carlo方法计算Pi值,随机数的定义和特性,什么是随机数？,单个的数字不是随机数;,是指一个数列，其中的每一个体

10、称为随机数，其值与数列中的其它数无关；,在一个均匀分布的随机数中，每一个体出现的概率是均等的；,例如：在0,1区间上均匀分布的随机数序列中，0.00001与0.5出现的机会均等,随机数应具有的基本特性,随机数序列应是独立的、互不相关的(uncorrelated)：,即序列中的任一子序列应与其它的子序列无关；,长的周期(long period)：,均匀分布的随机数应满足均匀性(Uniformity)：,随机数序列应是均匀的、无偏的，即：如果两个子区间的“面积”相等，则落于这两个子区间内的随机数的个数应相等。,例如：对0,1)区间均匀分布的随机数，如果产生了足够多的随机数，而有一半的随机数落于区间

11、0,0.1不满足均匀性,如果均匀性不满足，则会出现序列中的多组随机数相关的情况均匀性与互不相关的特性是有联系的,实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；,有效性（Efficiency):,模拟结果可靠,模拟产生的样本容量大,所需的随机数的数量大,随机数的产生必须快速、有效，最好能够进行并行计算。,随机数与随机数发生器,得到一个可能的随机数序列，是在计算机上实现Monte Carlo方法的关键,随机数的产生方法：,0,1区间上均匀分布的随机数是Monte Carlo模拟的基础，服从任意分布的随机数序列可以用0,1区间均匀分布的随机数序列作适当

12、的变换或舍选后求得。,(0,1) 2 -1(-1,1),利用随机数表，如Tippett于1972年发表的随机数表； 占用太多的计算机内存采用物理方法，如利用电子线路的热噪声等； 昂贵而且不便重复,:,伪随机数（Pseudo-Random Number),递推到一定次数后，出现周期性的重复现象。,利用数学递推公式,一旦公式和初值定下来，整个随机数序列便被确定下来，而且每一个随机数只被它前面的那个数唯一确定，因此这类随机数并不是真正的随机数。,Monte Carlo方法基本思想,当所求的问题是某种事件出现的概率，或是某个随机变量的期望值时，它们可以通过某种“随机试验”的方法，得到这种事件出现的频率

13、和概率，或者得到这个随机变量的统计平均值，并用它们作为问题的解。,Monte Carlo方法解决的问题,问题本身是确定性问题，要求我们去寻找一个随机过程，使该随机过程的统计平均就是所求问题的解。,问题本身就是随机过程，我们可以根据问题本身的实际物理过程来进行计算机模拟和跟踪，并采用统计方法求得问题的解。,Monte Carlo方法的特点,计算的收敛性和收敛速度均与问题的维数无关，适合解决高维问题。,对问题的适应能力强。,收敛速度仅为样本数的-1/2次，因而计算耗时大。,Monte Carlo方法的应用举例：,计算积分：,常用的积分方法求解：,将积分区域a,b均匀地划分成N各分区间，则积分结果可

14、近似地表示成：,x = (b-a)/N,简单的Monte Carlo积分方法求解：,利用均匀分布的随机数发生器，从a,b区间产生一系列随机数xi，i=1, 2, ., N,其中 X为均匀分布，并且 Xa,b,近似求解Eg(X):,近似求解积分:,随机抽样,当我们用简单Monte Carlo计算积分时，若该函数为常数函数，g(x)=constant，则取样数不管多少，准确度为100。如果在积分区间内，g(x)为一平滑函数，则简单Monte Carlo方法较为准确，反之，如果g(x)的变动很剧烈，则简单Monte Carlo方法的误差会变大。,说明：,重要性Monte Carlo抽样方法,在 g(

15、x) 变化剧烈时，如果以Monte Carlo方法取样，最好依据g(x)的大小来决定取样率。当|g(x)|的值较大时，对g(x)dx的贡献也较大，如果没被选中，则结果的误差极大。解决方式：改变 x被选中的机率，让|g(x)| 值较大的点被选中的机率增加。采用权重分布函数(Weight distribution function) w(x) :决定每个x被选中的机率。,重要性抽样的定义：根据一定的分布形式进行的随机抽样。,w(x)必须归一化，即在积分区间内w(x)dx=1。由于 x 的选取已被 w(x) 扭曲，所以计算积分时要把这部分还回去：若一共取样了N个x，则积分值为：,重要性Monte C

16、arlo抽样方法,Metropolis Monte Carlo方法,我们所模拟的系统最终要达到的平衡分布是Boltzman分布：,Boltzmann 概率分布函数：,我们如果能够产生这种分布，我们就能够计算系统的大多数性质，但这是不可能的，因为我们不知道Z的值，但是对于任意两个状态，我们有：,可以在相空间中构造一个马尔科夫链，使相空间中的样本点随着链的增长逐步趋近于Boltzman分布。,一个序列x0, x1, x2, ,xn,如果对任何n都有：,则此序列是一个Markov链。,要求：,任何一次实验的结果依赖于前一次的试验，并且近依赖于前一次的试验。,马尔科夫(Markov)链,可以证明：通过

17、构造Markov链，体系中最终的平衡分布就是Boltzman分布,Metropolis平衡条件 (Detailed balance condition)：,平衡条件：,系统处于状态X的概率正比于其Boltzman因子：,如果是对称的:,Metropolis Monte Carlo方法的算法：,给出一个初始状态，并计算系统的能量Eold,随机产生一个新状态，并计算新系统的能量Enew,如果E(EnewEold)0, 则接受新状态并回到 b),如果E(EnewEold)0, 则计算Boltzman因子：,在(0, 1)区间上产生一个均匀分布的随机数 ；,如果,则接受新状态并回到b),否则保留原值并

18、回到b),1. 正则系综蒙特卡罗模拟方法(Canonical MC Simulation),具有确定的粒子数N、温度T和体积V,对于含N个粒子的系统，位型（构型）的配分函数：,某个特定构型的发生概率为 PNVT(rN),1. 正则系综蒙特卡罗模拟方法(Canonical MC Simulation),Monte Carlo模拟中任一物理量的计算：,位型积分,概率密度,系统处于位型rN的概率密度,给出一个初始状态，并计算系统的能量Uold,随机产生一个新状态，并计算新系统的能量Unew,如果U(UnewUold)0, 则接受新状态并回到 b),如果U(UnewUold)0, 则计算Boltzma

19、n因子：,在(0, 1)区间上产生一个均匀分布的随机数 ；,如果,则接受新状态并回到b),否则保留原值并回到b),正则系综MC模拟算法的组织：,正则系综MC模拟算法的流程图：,给定每个分子的初始位置，ri(0),随机选取一个分子，并随机移动到新的位置,计算移动前后的系统能量变化U,拒绝移动,U0 ?,Exp(U) (0,1)?,统计系统的热力学性质及其它物理量,统计性质不变？,打印结果，结束,Yes,No,No,No,接受移动,Yes,Yes,大约循环107到108次,Monte Carlo模拟中几个热力学量的计算：,N个粒子系统中的总势能：,假设采用截断势能函数：,Uc：截断范围内的总势能；

20、Ulrc：截断半径外对势能的长程校正(Long-range correction),对于LJ流体：,含N个粒子系统中的压力：,Wc：截断范围内的总维里项(Virial)；,Plrc：截断半径外对压力的长程校正,Read simulationparameters,Start,Initialize positionsof all particles,New simulation?,Read oldconfiguration,Monte Carloloop,Stop,yes,no,Monte Carlo loop Subroutine,Start,Stop,Trial move,Satisfy Me

21、tropolisrule?,Accept thetrial move,Update energyand virial,Sample thepressure,End ofsimulation?,yes,no,yes,no,Main program,正则系综MC模拟程序基本结构：,正则系综MC模拟程序F11讲解（LJ, NVT）：,* READ INPUT DATA *,初始状态：,READ (*,(A) Title ! 运行作业题目READ (*,*) NStep ! 运行步数READ (*,*) Iprint ! 打印步数READ (*,*) Isave ! 保存步数READ (*,*) Ir

22、atio ! 调整步数READ (*,(A) CNFile ! 位型文件READ (*,*) Dens ! 对比密度READ (*,*) Temp ! 对比温度 READ (*,*) Rcut ! 对比截断半径,无因次量：,正则系综MC模拟程序F11讲解（LJ, NVT）：,量纲变换：,Beta = 1.0 / TempSigma = ( Dens / Real ( N ) ) * ( 1.0 / 3.0 )Rmin = 0.70 * Sigma ！判断粒子发生重叠时的距离Rcut = Rcut * Sigma ！截断半径DRmax = 0.15 * Sigma ！随机移动的最大距离DensL

23、J = DensDens = Dens / ( Sigma * 3 )IF ( Rcut .GT. 0.5 ) Stop Cut-Off Too Large ,模拟盒子的边长为1,* Read Initial Configuration *Call ReadCN（CNFile）,正则系综MC模拟程序F11讲解（LJ, NVT）：,初始位型：,参阅程序F23,Call FCC,需要自己给定所有粒子初始位置,面心立方 (face-centered cubic, FCC)：,正则系综MC模拟程序讲解（LJ, NVT）：,长程校正：,Sr3 = ( Sigma / Rcut ) * 3Sr9 = Sr

24、3 * 3Vlrc12 = 8.0 * Pi * DensLJ * Real ( N ) * Sr9 / 9.0Vlrc6 = - 8.0 * Pi * DensLJ * Real ( N ) * Sr3 / 3.0Vlrc = Vlrc12 + Vlrc6Wlrc12 = 4.0 * Vlrc12Wlrc6 = 2.0 * Vlrc6Wlrc = Wlrc12 + Wlrc6,算法：能量求和：,Call Sumup ( Rcut, Rmin, Sigma, Ovrlap, V, W )If ( Ovrlap ) Stop Overlap In Initial Configuration Vs

25、 = ( V + Vlrc ) / Real ( N )Ws = ( W + Wlrc ) / Real ( N )Ps = Dens * Temp + W + WlrcPs = Ps * Sigma * 3,* Check For Acceptance *DeltV = Vnew - VoldDeltW = Wnew - WoldDeltVb = Beta * DeltVIf ( DeltVb .Lt. 75.0 ) Then If ( DeltV .Le. 0.0 ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew . Acatma = Ac

26、atma + 1.0 Elseif ( Exp ( - DeltVb ) .Gt. Ranf ( Dummy ) ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew Acatma = Acatma + 1.0 EndifEndifAcm = Acm + 1.0,算法：Metropolis算法,算法：势能的计算,DO 100 J = 1, NIF ( I .NE. J ) ThenRXIJ = RXI - RX(J) RXIJ = RXIJ - Anint ( RXIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + R

27、ZIJ * RZIJIF ( RIJSQ .LT. RcutSQ ) THENSR2 = SIGSQ / RIJSQSR6 = SR2 * SR2 * SR2VIJ = SR6 * ( SR6 - 1.0 )WIJ = SR6 * ( SR6 - 0.5 )V = V + VIJW = W + WIJ ENDIFENDIF100 ContinueV = 4.0 * VW = 48.0 * W / 3.0,最小影像原理,算法：Metropolis算法,RXIOld = RX(I) .* Calculate The Energy Of I In The Old Configuration *Cal

28、l Energy ( RXIold, RYIold, RZIold, I, Rcut, Sigma, Vold, Wold )C * Move I And Pickup The Central Image *RXInew = RXIold + ( 2.0 * Ranf ( DUMMY ) - 1.0 ) * DRmax. RXInew = RXInew - Anint ( RXInew )* Calculate The Energy Of I In The New Configuration *Call Energy ( RXInew, RYInew, RZInew, I, Rcut, Sig

29、ma,Vnew, Wnew ),随机移动,* Adjust Maximum Displacement *Ratio = Acatma / Real ( N * Iratio )If ( Ratio .Gt. 0.5 ) Then DRmax = DRmax * 1.05Else DRmax = DRmax * 0.95EndifAcatma = 0.0EndifIf ( Mod ( Step, Iprint ) .Eq. 0 ) Then,算法：步长的调整,分子模拟中径向分布函数g(r)的计算：,计算表达式：,物理意义：流体中距一个分子为 r 处出现另一个分子的几率密度，它反映流体中短程有序的

30、特点。,r,r+r,在模拟中通常在盒子长度一半的范围内，考察g(r)随距离的变化。,分子模拟中径向分布函数g(r)的算法：,第一步：计算球壳层间的距离，并初始化一些变量：,Ngr: g(r)的统计次数,Delg: 球壳层间的距离,nhist: 球壳层的层数,模拟开始时需要给定： nhist 以及统计g(r)的步数间隔。,分子模拟中径向分布函数g(r)的算法：,第二步：统计g(r):,分子模拟中径向分布函数g(r)的算法：,第三步：系综统计平均计算g(r):,2. 巨正则系综蒙特卡罗模拟方法 (Grand Canonical MC Simulation),恒定 V, T, 和 ，体系的粒子数发生

31、波动；可用于预测 EOS-type的性质，但主要是用来模拟吸附过程；系统的微观态的分布与正则系综类似；建立的随机过程须增添粒子数的随机增减；模拟是一个等化学势面上的Markov过程。,巨正则系综配分函数,对于原子系统，位型（构型）的配分函数：,其中， s 为标度坐标，r = V1/3 。概率密度为：,Metropolis GCMC algorithm,产生巨正则系综的马尔可夫链的过程涉及到三种不同的随机移动： 随机移动一个粒子，算法与正则系综相同； 随机添加一个粒子 随机删除一个粒子,三种随机移动方式的接受概率：,采用逸度（Fugacity）：,Panagiotopoulos in 1987

32、introduced a clever way to simulate coexisting phases without an interface,Gibbs 系综 MC (GEMC),Two simulation volumes,Thermodynamic contact without physical contact,MC simulation includes moves that couple the two simulation volumes,Particle exchange equilibrates chemical potential,Volume exchange eq

33、uilibrates pressure,the coupled moves enforce mass and volume balance,尤其适用于研究纯流体或混合物的相平衡问题；此方法不能用于涉及到非常稠密流体的相平衡问题；此方法能同时获得共存相的各自密度及其组成；此方法避免了共存相界面的问题。临界点附近由于波动很大，很难模拟，所以临界点常采用其它方式确定。,GEMC模拟方法的特点:,GEMC 的配分函数,对于原子系统，位型（构型）的配分函数:,N=N1+N2V=V1+V2 constant T,GEMC 模拟算法：,随机选择一个粒子进行移动 (NVT).改变每个模拟盒子的体积，但总体积保

34、持不变 (NTP).盒子间交换粒子 (VT)。,Phase Behavior of Alkanes,Panagiotopoulos group,Alkane Mixtures,Siepmann group:乙烷庚烷Ethane + n-heptane,True Random Numbershttp:/www.random.org/http:/www.fourmilab.ch/hotbits/http:/ Number Generatorshttp:/random.mat.sbg.ac.at/http:/www.math.utah.edu/alfeld/Random/Random.htmlhttp:/,随机发生器的资源,Monte Carlo方法的资源,http:/csep1.phy.ornl.gov/mc/mc.htmlhttp:/www.chem.unl.edu/zeng/joy/mclab/mcintro.htmlhttp:/csep1.phy.ornl.gov/rn/rn.htmlhttp:/,