分子模拟教程(径向分布函数)

第三章分子模拟方法1.蒙特卡罗(MonteCarlo)方法基础2.分子动力学(MolecularDynamics)方法基础3.程序讲解本章具体讲解内容掌握分子模拟方法的必备知识：编程技能(FortranorC/C++)统计物理学（统计力学）：统计物理学基础；系综原理；非平衡统计力学基础；涨落理论分子热力学：分子间相互作用理论；分布函数理论气体分子运动论其它分子模拟的目的:为什么要进行分子模拟？1.将分子聚集体的性质与如下方面相联系：分子的微观相互作用分子聚集体的结构分子的动力学过程2.分子模拟对实验进行补充，使我们能够：预测现有或新材料的性质在分子水平研究宏观现象获得实验无法或难以发现的东西什么是计算机分子模拟方法？分子模拟的定义：统计力学基本原理出发，将一定数量的分子输入计算机内进行分子微观结构的测定和宏观性质的计算。按照获得微观态的方法不同，分子模拟分为：(1)蒙特卡罗方法(MonteCarlo,MC)(2)分子动力学方法(MolecularDynamics,MD)(3)混合方法(hybridmethod，HM)计算机分子模拟的发展历史：1.蒙特卡罗方法（MC）1953Metropolis,Ulam,RosenbluthandTell，LosAlamosNationalLabMonteCarlosimulationofhardsphere.2.分子动力学方法（MD）1957AlderandWainwrigth，LivermoreLabMoleculardynamicssimulationofhardspheres.微观与宏观分子模拟在微观尺度与实验室的宏观世界之间起着桥梁的作用：给定分子间的相互作用“准确”预测研究体系的性质MC与MD的区别：MC:–构型平均，不包含动力学部分；–利用概率行走产生微观态。MD:–时间平均，产生动力学性质；–利用运动轨线随时间的变化来产生一系列微观态。计算机分子模拟的发展历史（续）：从上个世纪九十年代初期以来，计算机模拟技术得到了飞速发展，主要基于三个方面的发展：分子力场的发展（基石）（Amber，OPLS、Compass）原子间的键长、键角、分子间的内聚能等模拟算法（途径）计算机硬件（工具）HPCx计算机分子模拟的特点：原子水平的模拟计算机实验检验理论、筛选实验科学研究中的第三种方法分子模拟中涉及的几个基本概念：模拟计算盒子或模拟胞腔Simulationbox(cell)装有一定数目流体分子的研究对象，它是我们要研究的宏观体系的缩微模型。立方形胞腔周期性边界条件(Periodicboundarycondition,PBC)在小体系中，边界效应总是很显著。在包含1000个原子的简单立方晶体中－488个原子处于边界上。在包含1000000个原子的简单立方晶体中－仍然有6%的原子在边界上。在模拟中，考虑具有真实边界的对象，不切合实际：•增强了有限尺寸效应•人为造成的边界会影响流体的性质当某个粒子运动出模拟盒子的某一边界时，另外一个影像粒子从另一对立边界进入到此盒子中。周期性边界条件(Periodicboundarycondition,PBC)本体体系的近似：中心盒子在X，Y和Z方向无限扩展；消除人为形成边界的表面效应；保证中心盒子中的粒子数恒定。只需要跟踪中心盒子中各粒子的运动。周期性边界条件的算法：LrLrrxxx0LrrrxxxyxLo)]LFLOOR(Dble[Lrrr采用数学函数：FLOOR(r/L):返回不超过r/L的最大整数FLOOR(4.8)hasthevalue4.FLOOR(-5.6)hasthevalue-6.xL/2o-L/22/LrrLrxxx2/LrrLrxxx)LANINT(Lrrr采用数学函数：r/L0,ANINT(r/L)=AINT(r/L+0.5)r/L0,ANINT(r/L)=AINT(r/L-0.5)周期性边界条件的算法：y最小影像转化原理(Minimumimageconvention)定义：中心元胞中的一个粒子只与此元胞中的其它N－1个粒子，或它们的最近邻影像发生相互作用。适用条件：粒子间相互作用势能的截断距离必须不大于模拟中心元胞长度的一半。此两粒子与中心粒子的距离相等，但是：黑色球发生作用绿色球不发生作用此两粒子是与中心原子相互作用的最近邻影像最小影像转化原理的算法：2/,,,LrrLrxijxijxij采用数学函数：2/,,,LrrLrxijxijxij)LANINT(Lrrrr/L0,ANINT(r/L)=AINT(r/L+0.5)r/L0,ANINT(r/L)=AINT(r/L-0.5)截断势能(TruncatingthePotential)本体体系采用周期性边界条件描述：–不可能将所有粒子与它们影像粒子间的相互作用全都计算。–必须在不大于中心盒子长度的一半处进行截断，以便与最小影像转化原理一致。粒子间的相互作用主要来自于截断范围内，而范围外的贡献很小，可忽略不计。截断范围内的相互作用截断势能函数的形式：①简单截断势能函数(TruncatedPotential)：缺点：rc:截断距离或半径•势能在截断处不连续,当一对分子穿越边界时，总能量不守恒。•分子间力在截断处为无穷大，MD运动过程不稳定。cccrrrrrUrU0)()(忽略截断半径之外的所有作用②位移截断势能函数(ShiftedandTruncatedPotential)：cccsrrrrrUrUrU0)()()(缺点：分子间力仍然在截断处不连续。优点：•势能在截断处连续，但不影响分子间力的大小•分子间力在截断处不为无穷大截断势能函数的形式：常用于MC和MD模拟中③位移－力截断势能函数(Shifted-ForcePotential)：常用于MD模拟中优点：•势能和分子间力均在截断处连续截断势能函数的形式：cccrrcsfrrrrrrdrrdUrUrUrUc0)()()()()(截断势能函数的对比：位移－力截断势能简单截断势能函数一、MonteCarlo模拟方法基础：亦称统计模拟或随机抽样方法，statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来：•是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；NicholasMetropolis(1915-1999)Monte-Carlo,Monaco投硬币，掷骰子MonteCarlo方法计算Pi值随机数的定义和特性什么是随机数？单个的数字不是随机数;是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；在一个均匀分布的随机数中，每一个体出现的概率是均等的；例如：在[0,1]区间上均匀分布的随机数序列中，0.00001与0.5出现的机会均等随机数应具有的基本特性随机数序列应是独立的、互不相关的(uncorrelated)：即序列中的任一子序列应与其它的子序列无关；长的周期(longperiod)：均匀分布的随机数应满足均匀性(Uniformity)：随机数序列应是均匀的、无偏的，即：如果两个子区间的“面积”相等，则落于这两个子区间内的随机数的个数应相等。例如：对[0,1)区间均匀分布的随机数，如果产生了足够多的随机数，而有一半的随机数落于区间[0,0.1]不满足均匀性如果均匀性不满足，则会出现序列中的多组随机数相关的情况均匀性与互不相关的特性是有联系的实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；有效性（Efficiency):模拟结果可靠模拟产生的样本容量大所需的随机数的数量大随机数的产生必须快速、有效，最好能够进行并行计算。随机数与随机数发生器①得到一个可能的随机数序列，是在计算机上实现MonteCarlo方法的关键②随机数的产生方法：[0,1]区间上均匀分布的随机数是MonteCarlo模拟的基础，服从任意分布的随机数序列可以用[0,1]区间均匀分布的随机数序列作适当的变换或舍选后求得。(0,1)2-1(-1,1)利用随机数表，如Tippett于1972年发表的随机数表；占用太多的计算机内存采用物理方法，如利用电子线路的热噪声等；昂贵而且不便重复:伪随机数（Pseudo-RandomNumber)递推到一定次数后，出现周期性的重复现象。利用数学递推公式一旦公式和初值定下来，整个随机数序列便被确定下来，而且每一个随机数只被它前面的那个数唯一确定，因此这类随机数并不是真正的随机数。MonteCarlo方法基本思想当所求的问题是某种事件出现的概率，或是某个随机变量的期望值时，它们可以通过某种“随机试验”的方法，得到这种事件出现的频率和概率，或者得到这个随机变量的统计平均值，并用它们作为问题的解。MonteCarlo方法解决的问题•问题本身是确定性问题，要求我们去寻找一个随机过程，使该随机过程的统计平均就是所求问题的解。•问题本身就是随机过程，我们可以根据问题本身的实际物理过程来进行计算机模拟和跟踪，并采用统计方法求得问题的解。MonteCarlo方法的特点•计算的收敛性和收敛速度均与问题的维数无关，适合解决高维问题。•对问题的适应能力强。•收敛速度仅为样本数的-1/2次，因而计算耗时大。MonteCarlo方法的应用举例：计算积分：baIdxg(x)NiiNNiiNxfNabxxfI11limΔlim①常用的积分方法求解：将积分区域[a,b]均匀地划分成N各分区间，则积分结果可近似地表示成：Δx=(b-a)/N,Δ1xi-xi②简单的MonteCarlo积分方法求解：利用均匀分布的随机数发生器，从[a,b]区间产生一系列随机数xi，i=1,2,...,Ndxa-b1g(x)a)-(bdxg(x)Ibabag(x)a)-(bI其中X为均匀分布，并且X[a,b]N1i)limig(xN1g(x)N近似求解E[g(X)]:g(x)a)-(bI近似求解积分:随机抽样当我们用简单MonteCarlo计算积分时，若该函数为常数函数，g(x)=constant，则取样数不管多少，准确度为100％。如果在积分区间内，g(x)为一平滑函数，则简单MonteCarlo方法较为准确，反之，如果g(x)的变动很剧烈，则简单MonteCarlo方法的误差会变大。说明：)exp()(2xxf)100exp()(2xxf②重要性MonteCarlo抽样方法在g(x)变化剧烈时，如果以MonteCarlo方法取样，最好依据g(x)的大小来决定取样率。当|g(x)|的值较大时，对∫g(x)dx的贡献也较大，如果没被选中，则结果的误差极大。解决方式：改变x被选中的机率，让|g(x)|值较大的点被选中的机率增加。采用权重分布函数(Weightdistributionfunction)w(x):决定每个x被选中的机率。重要性抽样的定义：根据一定的分布形式进行的随机抽样。w(x)必须归一化，即在积分区间内∫w(x)dx=1。由于x的选取已被w(x)扭曲，所以计算积分时要把这部分[还]回去：若一共取样了N个x，则积分值为：NiiixwxgNdxxwxgxwdxxgI1)()(1)()()()(②重要性MonteCarlo抽样方法③MetropolisMonteCarlo方法我们所模拟的系统最终要达到的平衡分布是Boltzman分布：Boltzmann概率分布函数：jkTEkTEijieepZepkTEii我们如果能够产生这种分布，我们就能够计算系统的大多数性质，但这是不可能的，因为我们不知道Z的值，但是对于任意两个状态