分子模拟教程课件

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

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

3、d Tell， Los Alamos National Lab Monte Carlo simulation of hard sphere. 2. 分子动力学方法（MD）1957Alder and Wainwrigth， Livermore Lab Molecular dynamics simulation of hard spheres.,5,PPT学习交流,微观与宏观,分子模拟在微观尺度与实验室的宏观世界之间起着桥梁的作用：,给定分子间的相互作用 “准确”预测研究体系的性质,6,PPT学习交流,MC与MD的区别：,MC: 构型平均，不包含动力学部分； 利用概率行走产生微观态。 MD: 时间

4、平均，产生动力学性质； 利用运动轨线随时间的变化来产生一系列微观态。,7,PPT学习交流,计算机分子模拟的发展历史（续）：,从上个世纪九十年代初期以来，计算机模拟技术得到了飞速发展，主要基于三个方面的发展： 分子力场的发展（基石） （Amber，OPLS、Compass） 原子间的键长、键角、分子间的内聚能等 模拟算法（途径） 计算机硬件（工具）,HPCx,8,PPT学习交流,计算机分子模拟的特点：,原子水平的模拟 计算机实验 检验理论、筛选实验 科学研究中的第三种方法,9,PPT学习交流,分子模拟中涉及的几个基本概念：,模拟计算盒子或模拟胞腔,Simulation box (cell),装有

5、一定数目流体分子的研究对象，它是我们要研究的宏观体系的缩微模型。,立方形胞腔,10,PPT学习交流,周期性边界条件(Periodic boundary condition, PBC),在小体系中，边界效应总是很显著。 在包含1000个原子的简单立方晶体中488个原子处于边界上。 在包含个原子的简单立方晶体中仍然有 6%的原子在边界上。,在模拟中，考虑具有真实边界的对象，不切合实际： 增强了有限尺寸效应 人为造成的边界会影响流体的性质,11,PPT学习交流,当某个粒子运动出模拟盒子的某一边界时，另外一个影像粒子从另一对立边界进入到此盒子中。,周期性边界条件(Periodic boundary c

6、ondition, PBC),本体体系的近似：中心盒子在X，Y和Z方向无限扩展； 消除人为形成边界的表面效应； 保证中心盒子中的粒子数恒定。 只需要跟踪中心盒子中各粒子的运动。,12,PPT学习交流,周期性边界条件的算法：,采用数学函数：,FLOOR(r/L): 返回不超过r/L的最大整数,FLOOR (4.8) has the value 4. FLOOR (-5.6) has the value -6.,13,PPT学习交流,采用数学函数：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),周期性边界条件的算

7、法：,y,14,PPT学习交流,最小影像转化原理(Minimum image convention),定义：,中心元胞中的一个粒子只与此元胞中的其它N1个粒子，或它们的最近邻影像发生相互作用。,适用条件：,粒子间相互作用势能的截断距离必须不大于模拟中心元胞长度的一半。,此两粒子与中心粒子的距离相等，但是： 黑色球发生作用绿色球不发生作用,此两粒子是与中心原子相互作用的最近邻影像,15,PPT学习交流,最小影像转化原理的算法：,采用数学函数：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),16,PPT学习交流,

8、截断势能(Truncating the Potential),本体体系采用周期性边界条件描述： 不可能将所有粒子与它们影像粒子间的相互作用全都计算。 必须在不大于中心盒子长度的一半处进行截断，以便与最小影像转化原理一致。 粒子间的相互作用主要来自于截断范围内，而范围外的贡献很小，可忽略不计。,截断范围内的相互作用,17,PPT学习交流,截断势能函数的形式：,简单截断势能函数(Truncated Potential)：,缺点：,rc: 截断距离或半径,势能在截断处不连续,当一对分子穿越边界时，总能量不守恒。 分子间力在截断处为无穷大，MD运动过程不稳定。,忽略截断半径之外的所有作用,18,PPT

9、学习交流,位移截断势能函数(Shifted and Truncated Potential)：,缺点：,分子间力仍然在截断处不连续。,优点：,势能在截断处连续，但不影响分子间力的大小 分子间力在截断处不为无穷大,截断势能函数的形式：,常用于MC和MD模拟中,19,PPT学习交流,位移力截断势能函数(Shifted-Force Potential)：,常用于MD模拟中,优点：,势能和分子间力均在截断处连续,截断势能函数的形式：,20,PPT学习交流,截断势能函数的对比：,位移力截断势能,简单截断势能函数,21,PPT学习交流,一、Monte Carlo模拟方法基础：,亦称统计模拟或随机抽样方法，

10、statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,投硬币，掷骰子,22,PPT学习交流,Monte Carlo方法计算Pi值,23,PPT学习交流,随机数的定义和特性,什么是随机数？,单个的数字不是随机数;,是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；,在一个均匀分布的随机数中，每一个体出

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

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

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

14、试验”的方法，得到这种事件出现的频率和概率，或者得到这个随机变量的统计平均值，并用它们作为问题的解。,Monte Carlo方法解决的问题,问题本身是确定性问题，要求我们去寻找一个随机过程，使该随机过程的统计平均就是所求问题的解。,问题本身就是随机过程，我们可以根据问题本身的实际物理过程来进行计算机模拟和跟踪，并采用统计方法求得问题的解。,29,PPT学习交流,Monte Carlo方法的特点,计算的收敛性和收敛速度均与问题的维数无关，适合解决高维问题。,对问题的适应能力强。,收敛速度仅为样本数的-1/2次，因而计算耗时大。,30,PPT学习交流,Monte Carlo方法的应用举例：,计算积

15、分：,常用的积分方法求解：,将积分区域a,b均匀地划分成N各分区间，则积分结果可近似地表示成：,x = (b-a)/N,31,PPT学习交流,简单的Monte Carlo积分方法求解：,利用均匀分布的随机数发生器，从a,b区间产生一系列随机数xi，i=1, 2, ., N,其中 X为均匀分布，并且 Xa,b,近似求解Eg(X):,近似求解积分:,随机抽样,32,PPT学习交流,当我们用简单Monte Carlo计算积分时，若该函数为常数函数，g(x)=constant，则取样数不管多少，准确度为100。 如果在积分区间内，g(x)为一平滑函数，则简单Monte Carlo方法较为准确，反之，如

16、果g(x)的变动很剧烈，则简单Monte Carlo方法的误差会变大。,说明：,33,PPT学习交流,重要性Monte Carlo抽样方法,在 g(x) 变化剧烈时，如果以Monte Carlo方法取样，最好依据g(x)的大小来决定取样率。 当|g(x)|的值较大时，对g(x)dx的贡献也较大，如果没被选中，则结果的误差极大。 解决方式：改变 x被选中的机率，让|g(x)| 值较大的点被选中的机率增加。 采用权重分布函数(Weight distribution function) w(x) :决定每个x被选中的机率。,重要性抽样的定义：根据一定的分布形式进行的随机抽样。,34,PPT学习交流,

17、w(x)必须归一化，即在积分区间内w(x)dx=1。 由于 x 的选取已被 w(x) 扭曲，所以计算积分时要把这部分还回去：若一共取样了N个x，则积分值为：,重要性Monte Carlo抽样方法,35,PPT学习交流,Metropolis Monte Carlo方法,我们所模拟的系统最终要达到的平衡分布是Boltzman分布：,Boltzmann 概率分布函数：,我们如果能够产生这种分布，我们就能够计算系统的大多数性质，但这是不可能的，因为我们不知道Z的值，但是对于任意两个状态，我们有：,可以在相空间中构造一个马尔科夫链，使相空间中的样本点随着链的增长逐步趋近于Boltzman分布。,36,P

18、PT学习交流,一个序列x0, x1, x2, ,xn,如果对任何n都有：,则此序列是一个Markov链。,要求：,任何一次实验的结果依赖于前一次的试验，并且近依赖于前一次的试验。,马尔科夫(Markov)链,可以证明：通过构造Markov链，体系中最终的平衡分布就是Boltzman分布,37,PPT学习交流,Metropolis平衡条件 (Detailed balance condition)：,平衡条件：,系统处于状态X的概率正比于其Boltzman因子：,如果是对称的:,38,PPT学习交流,Metropolis Monte Carlo方法的算法：,给出一个初始状态，并计算系统的能量Eol

19、d,随机产生一个新状态，并计算新系统的能量Enew,如果E(EnewEold)0, 则接受新状态并回到 b),如果E(EnewEold)0, 则计算Boltzman因子：,在(0, 1)区间上产生一个均匀分布的随机数 ；,如果,则接受新状态并回到b),否则保留原值并回到b),39,PPT学习交流,1. 正则系综蒙特卡罗模拟方法(Canonical MC Simulation),具有确定的粒子数N、温度T和体积V,40,PPT学习交流,对于含N个粒子的系统，位型（构型）的配分函数：,某个特定构型的发生概率为 PNVT(rN),1. 正则系综蒙特卡罗模拟方法(Canonical MC Simula

20、tion),41,PPT学习交流,Monte Carlo模拟中任一物理量的计算：,位型积分,概率密度,系统处于位型rN的概率密度,42,PPT学习交流,给出一个初始状态，并计算系统的能量Uold,随机产生一个新状态，并计算新系统的能量Unew,如果U(UnewUold)0, 则接受新状态并回到 b),如果U(UnewUold)0, 则计算Boltzman因子：,在(0, 1)区间上产生一个均匀分布的随机数 ；,如果,则接受新状态并回到b),否则保留原值并回到b),正则系综MC模拟算法的组织：,43,PPT学习交流,正则系综MC模拟算法的流程图：,给定每个分子的初始位置，ri(0),随机选取一个

21、分子，并随机移动到新的位置,计算移动前后的系统能量变化U,拒绝移动,U0 ?,Exp(U) (0,1)?,统计系统的热力学性质及其它物理量,统计性质不变？,打印结果，结束,Yes,No,No,No,接受移动,Yes,Yes,大约循环107到108次,44,PPT学习交流,Monte Carlo模拟中几个热力学量的计算：,N个粒子系统中的总势能：,假设采用截断势能函数：,Uc：截断范围内的总势能； Ulrc：截断半径外对势能的长程校正(Long-range correction),对于LJ流体：,45,PPT学习交流,含N个粒子系统中的压力：,Wc：截断范围内的总维里项(Virial)；,Plr

22、c：截断半径外对压力的长程校正,46,PPT学习交流,Read simulation parameters,Start,Initialize positions of all particles,New simulation?,Read old configuration,Monte Carlo loop,Stop,yes,no,Monte Carlo loop Subroutine,Start,Stop,Trial move,Satisfy Metropolis rule?,Accept the trial move,Update energy and virial,Sample the p

23、ressure,End of simulation?,yes,no,yes,no,Main program,正则系综MC模拟程序基本结构：,47,PPT学习交流,正则系综MC模拟程序F11讲解（LJ, NVT）：,* READ INPUT DATA *,初始状态：,READ (*,(A) Title ! 运行作业题目 READ (*,*) NStep ! 运行步数 READ (*,*) Iprint ! 打印步数 READ (*,*) Isave ! 保存步数 READ (*,*) Iratio ! 调整步数 READ (*,(A) CNFile ! 位型文件 READ (*,*) Dens

24、! 对比密度 READ (*,*) Temp ! 对比温度 READ (*,*) Rcut ! 对比截断半径,48,PPT学习交流,无因次量：,49,PPT学习交流,正则系综MC模拟程序F11讲解（LJ, NVT）：,量纲变换：,Beta = 1.0 / Temp Sigma = ( Dens / Real ( N ) ) * ( 1.0 / 3.0 ) Rmin = 0.70 * Sigma ！判断粒子发生重叠时的距离 Rcut = Rcut * Sigma ！截断半径 DRmax = 0.15 * Sigma ！随机移动的最大距离 DensLJ = Dens Dens = Dens / (

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

26、Sr3 * 3 Vlrc12 = 8.0 * Pi * DensLJ * Real ( N ) * Sr9 / 9.0 Vlrc6 = - 8.0 * Pi * DensLJ * Real ( N ) * Sr3 / 3.0 Vlrc = Vlrc12 + Vlrc6 Wlrc12 = 4.0 * Vlrc12 Wlrc6 = 2.0 * Vlrc6 Wlrc = Wlrc12 + Wlrc6,52,PPT学习交流,算法：能量求和：,Call Sumup ( Rcut, Rmin, Sigma, Ovrlap, V, W ) If ( Ovrlap ) Stop Overlap In Init

27、ial Configuration Vs = ( V + Vlrc ) / Real ( N ) Ws = ( W + Wlrc ) / Real ( N ) Ps = Dens * Temp + W + Wlrc Ps = Ps * Sigma * 3,53,PPT学习交流,* Check For Acceptance * DeltV = Vnew - Vold DeltW = Wnew - Wold DeltVb = Beta * DeltV If ( DeltVb .Lt. 75.0 ) Then If ( DeltV .Le. 0.0 ) Then V = V + DeltV W =

28、W + DeltW RX(i) = RXInew . Acatma = Acatma + 1.0 Elseif ( Exp ( - DeltVb ) .Gt. Ranf ( Dummy ) ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew Acatma = Acatma + 1.0 Endif Endif Acm = Acm + 1.0,算法：Metropolis算法,54,PPT学习交流,算法：势能的计算,DO 100 J = 1, N IF ( I .NE. J ) Then RXIJ = RXI - RX(J) RXIJ = RXIJ -

29、 Anint ( RXIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ IF ( RIJSQ .LT. RcutSQ ) THEN SR2 = SIGSQ / RIJSQ SR6 = SR2 * SR2 * SR2 VIJ = SR6 * ( SR6 - 1.0 ) WIJ = SR6 * ( SR6 - 0.5 ) V = V + VIJ W = W + WIJ ENDIF ENDIF 100 Continue V = 4.0 * V W = 48.0 * W / 3.0,最小影像原理,55,PPT学习交流,算法：Metropolis算

30、法,RXIOld = RX(I) . * Calculate The Energy Of I In The Old Configuration * Call 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

31、Of I In The New Configuration * Call Energy ( RXInew, RYInew, RZInew, I, Rcut, Sigma,Vnew, Wnew ),随机移动,56,PPT学习交流,* Adjust Maximum Displacement * Ratio = Acatma / Real ( N * Iratio ) If ( Ratio .Gt. 0.5 ) Then DRmax = DRmax * 1.05 Else DRmax = DRmax * 0.95 Endif Acatma = 0.0 Endif If ( Mod ( Step, I

32、print ) .Eq. 0 ) Then,算法：步长的调整,57,PPT学习交流,分子模拟中径向分布函数g(r)的计算：,计算表达式：,物理意义： 流体中距一个分子为 r 处出现另一个分子的几率密度，它反映流体中短程有序的特点。,r,r+r,在模拟中通常在盒子长度一半的范围内，考察g(r)随距离的变化。,58,PPT学习交流,分子模拟中径向分布函数g(r)的算法：,第一步：计算球壳层间的距离，并初始化一些变量：,Ngr: g(r)的统计次数,Delg: 球壳层间的距离,nhist: 球壳层的层数,模拟开始时需要给定： nhist 以及统计g(r)的步数间隔。,59,PPT学习交流,分子模拟中

33、径向分布函数g(r)的算法：,第二步：统计g(r):,60,PPT学习交流,分子模拟中径向分布函数g(r)的算法：,第三步：系综统计平均计算g(r):,61,PPT学习交流,62,PPT学习交流,2. 巨正则系综蒙特卡罗模拟方法 (Grand Canonical MC Simulation),恒定 V, T, 和 ，体系的粒子数发生波动； 可用于预测 EOS-type的性质，但主要是用来模拟吸附过程； 系统的微观态的分布与正则系综类似； 建立的随机过程须增添粒子数的随机增减； 模拟是一个等化学势面上的Markov过程。,63,PPT学习交流,巨正则系综配分函数,对于原子系统，位型（构型）的配分函数：,其中， s 为标度坐标，r = V1/3 。概率密度为：,64,PPT学习