《分子模拟教程(精品).ppt》由会员分享，可在线阅读，更多相关《分子模拟教程(精品).ppt（76页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第三章第三章 分子模拟方法分子模拟方法 1.蒙特卡罗蒙特卡罗(Monte Carlo)方法基础方法基础 2.分子动力学分子动力学(Molecular Dynamics)方法基础方法基础 3.程序讲解程序讲解 本章具体讲解内容本章具体讲解内容掌握分子模拟方法的必备知识：掌握分子模拟方法的必备知识：掌握分子模拟方法的必备知识：掌握分子模拟方法的必备知识：n n编程技能编程技能编程技能编程技能(Fortran or C/C+)(Fortran or C/C+)n n统计物理学（统计力学）统计物理学（统计力学）统计物理学（统计力学）统计物理学（统计力学）：统计物理学基础；统计物理学基础；统计物理学基础
2、；统计物理学基础；系综原理；系综原理；系综原理；系综原理；非平衡统计力学基础；非平衡统计力学基础；非平衡统计力学基础；非平衡统计力学基础；涨落理论涨落理论涨落理论涨落理论n n分子热力学分子热力学分子热力学分子热力学 ：分子间相互作用理论；分子间相互作用理论；分子间相互作用理论；分子间相互作用理论；分布函数理论分布函数理论分布函数理论分布函数理论n n气体分子运动论气体分子运动论气体分子运动论气体分子运动论n n其它其它其它其它分子模拟的目的分子模拟的目的:为什么要进行分子模拟？为什么要进行分子模拟？1.将分子聚集体的性质与如下方面相联系：将分子聚集体的性质与如下方面相联系：n分子的微观相互作
3、用分子的微观相互作用n分子聚集体的结构分子聚集体的结构n分子的动力学过程分子的动力学过程2.分子模拟对实验进行补充，使我们能够：分子模拟对实验进行补充，使我们能够：n预测现有或新材料的性质预测现有或新材料的性质n在分子水平研究宏观现象在分子水平研究宏观现象n获得实验无法或难以发现的东西获得实验无法或难以发现的东西什么是计算机分子模拟方法？什么是计算机分子模拟方法？什么是计算机分子模拟方法？什么是计算机分子模拟方法？分子模拟的定义：分子模拟的定义：分子模拟的定义：分子模拟的定义：统计力学基本原理出发，将一定数量的分子输入计算统计力学基本原理出发，将一定数量的分子输入计算统计力学基本原理出发，将一
4、定数量的分子输入计算统计力学基本原理出发，将一定数量的分子输入计算机内进行分子微观结构的测定和宏观性质的计算。机内进行分子微观结构的测定和宏观性质的计算。机内进行分子微观结构的测定和宏观性质的计算。机内进行分子微观结构的测定和宏观性质的计算。q 按照获得微观态的方法不同，分子模拟分为：按照获得微观态的方法不同，分子模拟分为：(1)蒙特卡罗方法蒙特卡罗方法(Monte Carlo,MC)(2)分子动力学方法分子动力学方法(Molecular Dynamics,MD)(3)混合方法混合方法(hybrid method，HM)计算机分子模拟的发展历史：计算机分子模拟的发展历史：计算机分子模拟的发展历
5、史：计算机分子模拟的发展历史：1.1.蒙特卡罗方法（蒙特卡罗方法（蒙特卡罗方法（蒙特卡罗方法（MCMC）1953 1953 Metropolis,Metropolis,UlamUlam,RosenbluthRosenbluth and Tell and Tell，Los Alamos National LabLos Alamos National LabMonte Carlo simulation of hard sphere.Monte Carlo simulation of hard sphere.2.2.分子动力学方法（分子动力学方法（分子动力学方法（分子动力学方法（MDMD）19571
6、957Alder and Alder and WainwrigthWainwrigth，Livermore LabLivermore LabMolecular dynamics simulation of hard Molecular dynamics simulation of hard spheres.spheres.微观与宏观微观与宏观微观与宏观微观与宏观分分子子模模拟拟在在微微观观尺尺度度与与实实验验室室的的宏宏观观世世界界之之间间起起着着桥梁的作用：桥梁的作用：u 给定分子间的相互作用给定分子间的相互作用u“准确准确”预测预测研究体系的性质研究体系的性质MCMC与与与与MDMD的区别
7、：的区别：的区别：的区别：n nMC:MC:构型平均，不包含动力学部分；构型平均，不包含动力学部分；构型平均，不包含动力学部分；构型平均，不包含动力学部分；利用概率行走产生微观态。利用概率行走产生微观态。利用概率行走产生微观态。利用概率行走产生微观态。n nMD:MD:时间平均，产生动力学性质；时间平均，产生动力学性质；时间平均，产生动力学性质；时间平均，产生动力学性质；利用运动轨线随时间的变化来产生一系列微观态。利用运动轨线随时间的变化来产生一系列微观态。利用运动轨线随时间的变化来产生一系列微观态。利用运动轨线随时间的变化来产生一系列微观态。计算机分子模拟的发展历史（续）：计算机分子模拟的发
8、展历史（续）：计算机分子模拟的发展历史（续）：计算机分子模拟的发展历史（续）：从上个世纪九十年代初期以来，计算机模拟技术得到了从上个世纪九十年代初期以来，计算机模拟技术得到了飞速发展，主要基于三个方面的发展：飞速发展，主要基于三个方面的发展：v分子力场的发展（基石）分子力场的发展（基石）（Amber，OPLS、Compass）原子间的键长、键角、分子间的内聚能等原子间的键长、键角、分子间的内聚能等v模拟算法（途径）模拟算法（途径）v计算机硬件（工具）计算机硬件（工具）HPCx计算机分子模拟的特点：计算机分子模拟的特点：计算机分子模拟的特点：计算机分子模拟的特点：n n原子水平的模拟原子水平的模
9、拟原子水平的模拟原子水平的模拟n n计算机实验计算机实验计算机实验计算机实验n n检验理论、筛选实验检验理论、筛选实验检验理论、筛选实验检验理论、筛选实验n n科学研究中的第三种方法科学研究中的第三种方法科学研究中的第三种方法科学研究中的第三种方法分子模拟中涉及的几个基本概念：分子模拟中涉及的几个基本概念：n 模拟计算盒子或模拟胞腔模拟计算盒子或模拟胞腔Simulation box(cell)装装有有一一定定数数目目流流体体分分子子的的研研究究对对象象，它它是是我我们们要要研研究究的的宏观体系的缩微模型。宏观体系的缩微模型。立方形胞腔立方形胞腔n 周期性边界条件周期性边界条件(Periodic
10、 boundary condition,PBC)u在小体系中，边界效应总是很显著。在小体系中，边界效应总是很显著。u在在包包含含1000个个原原子子的的简简单单立立方方晶晶体体中中488个原子处于边界上。个原子处于边界上。u在在包包含含1000000个个原原子子的的简简单单立立方方晶晶体体中中仍然有仍然有 6%的原子在边界上。的原子在边界上。在模拟中，考虑具有真实边界的对象，不切合实际：在模拟中，考虑具有真实边界的对象，不切合实际：增强了有限尺寸效应增强了有限尺寸效应 人为造成的边界会影响流体的性质人为造成的边界会影响流体的性质当某个粒子运动出模拟盒子的某一边界时，另外一个影像当某个粒子运动出
11、模拟盒子的某一边界时，另外一个影像粒子从另一对立边界进入到此盒子中。粒子从另一对立边界进入到此盒子中。n 周期性边界条件周期性边界条件(Periodic boundary condition,PBC)u本本体体体体系系的的近近似似：中中心心盒盒子子在在X，Y和和Z方方向向无无限扩展；限扩展；u消消除除人人为为形形成成边边界界的的表表面效应；面效应；u保保证证中中心心盒盒子子中中的的粒粒子子数恒定。数恒定。u只只需需要要跟跟踪踪中中心心盒盒子子中中各粒子的运动各粒子的运动。n 周期性边界条件的算法：周期性边界条件的算法：yxLo采用数学采用数学函数：函数：FLOOR(r/L):返回不超过返回不超
12、过r/L的最大整数的最大整数FLOOR(4.8)has the value 4.FLOOR(-5.6)has the value-6.xL/2o-L/2采用数学采用数学函数：函数：r/L0,ANINT(r/L)=AINT(r/L+0.5)r/L 0,ANINT(r/L)=AINT(r/L-0.5)n 周期性边界条件的算法：周期性边界条件的算法：yn 最小影像转化原理最小影像转化原理(Minimum image convention)定义：定义：定义：定义：中中中中心心心心元元元元胞胞胞胞中中中中的的的的一一一一个个个个粒粒粒粒子子子子只只只只与与与与此此此此元元元元胞胞胞胞中中中中的的的的其其
13、其其它它它它N N1 1个个个个粒粒粒粒子子子子，或它们的最近邻影像发生相互作用。或它们的最近邻影像发生相互作用。或它们的最近邻影像发生相互作用。或它们的最近邻影像发生相互作用。适用条件：适用条件：适用条件：适用条件：粒粒粒粒子子子子间间间间相相相相互互互互作作作作用用用用势势势势能能能能的的的的截截截截断断断断距距距距离离离离必必必必须须须须不不不不大大大大于于于于模模模模拟拟拟拟中中中中心心心心元元元元胞长度的一半。胞长度的一半。胞长度的一半。胞长度的一半。此两粒子与中心粒子的距此两粒子与中心粒子的距离相等，但是：离相等，但是：黑色球发生作用黑色球发生作用绿色球不发生作用绿色球不发生作用此
14、两粒子是与中心原子此两粒子是与中心原子相互作用的最近邻影像相互作用的最近邻影像n最小影像转化原理的算法：最小影像转化原理的算法：采用数学采用数学函数：函数：r/L0,ANINT(r/L)=AINT(r/L+0.5)r/L 0,ANINT(r/L)=AINT(r/L-0.5)n 截断势能截断势能(Truncating the Potential)n n本本本本体体体体体体体体系系系系采采采采用用用用周周周周期期期期性性性性边边边边界界界界条条条条件件件件描描描描述：述：述：述：不不不不可可可可能能能能将将将将所所所所有有有有粒粒粒粒子子子子与与与与它它它它们们们们影影影影像像像像粒粒粒粒子子子子
15、间间间间的的的的相相相相互互互互作作作作用用用用全全全全都都都都计计计计算。算。算。算。必必必必须须须须在在在在不不不不大大大大于于于于中中中中心心心心盒盒盒盒子子子子长长长长度度度度的的的的一一一一半半半半处处处处进进进进行行行行截截截截断断断断，以以以以便便便便与与与与最小影像转化原理最小影像转化原理一致。一致。一致。一致。n n粒粒粒粒子子子子间间间间的的的的相相相相互互互互作作作作用用用用主主主主要要要要来来来来自自自自于于于于截截截截断断断断范范范范围围围围内内内内，而而而而范范范范围围围围外外外外的的的的贡贡贡贡献献献献很很很很小小小小，可忽略不计。可忽略不计。可忽略不计。可忽略不
16、计。截断范围内截断范围内的的相互作用相互作用u截断势能函数的形式：截断势能函数的形式：简单截断势能函数简单截断势能函数(Truncated Potential)：缺点：缺点：rc:截断距离或半径截断距离或半径 势能在截断处不连续势能在截断处不连续,当一对分子穿越边界时，总能量不守恒。当一对分子穿越边界时，总能量不守恒。分子间力在截断处为无穷大，分子间力在截断处为无穷大，MD运动过程不稳定。运动过程不稳定。忽略截断半径之外的所有作用忽略截断半径之外的所有作用位移截断势能函数位移截断势能函数(Shifted and Truncated Potential)：缺点：缺点：分子间力分子间力仍然仍然在截
17、断处不连续。在截断处不连续。优点：优点：势能在截断处连续，但不影响分子间力的大小势能在截断处连续，但不影响分子间力的大小 分子间力在截断处不为无穷大分子间力在截断处不为无穷大u截断势能函数的形式：截断势能函数的形式：常用于常用于MC和和MD模拟中模拟中位移力截断势能函数位移力截断势能函数(Shifted-Force Potential)：常用于常用于MD模拟中模拟中优点：优点：势能和分子间力均在截断处连续势能和分子间力均在截断处连续u截断势能函数的形式：截断势能函数的形式：u截断势能函数的对比：截断势能函数的对比：位移力截断势能位移力截断势能简单截断势能函数简单截断势能函数一、一、Monte
18、Carlo模拟方法基础：模拟方法基础：亦称统计模拟或随机抽样方法，亦称统计模拟或随机抽样方法，statistical simulation method 利用利用随机数随机数进行数值模拟的方法进行数值模拟的方法Monte Carlo名字的由来：名字的由来：是由是由Metropolis在二次世界大战期间提出的：在二次世界大战期间提出的：Manhattan计划，计划，研究与原子弹有关的中子输运过程；研究与原子弹有关的中子输运过程；Nicholas Metropolis(1915-1999)Monte-Carlo,Monaco投硬币，掷骰子投硬币，掷骰子Monte CarloMonte CarloM
19、onte Carlo方法计算方法计算方法计算方法计算方法计算方法计算PiPiPi值值值值值值n n随机数的定义和特性随机数的定义和特性随机数的定义和特性随机数的定义和特性什么是随机数？什么是随机数？单个的数字不是随机数单个的数字不是随机数;是指一个数列，其中的每一个体称为随机数，其值与是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；数列中的其它数无关；在一个均匀分布的随机数中，每一个体出现的概率是在一个均匀分布的随机数中，每一个体出现的概率是均等的；均等的；v例如：在例如：在0,1区间上均匀分布的随机数序列中，区间上均匀分布的随机数序列中，0.00001与与0.5出现的机会均
20、等出现的机会均等随机数应具有的基本特性随机数应具有的基本特性v随机数序列应是独立的、互不相关的随机数序列应是独立的、互不相关的(uncorrelated)：即序列中的任一子序列应与其它的子序列无关；即序列中的任一子序列应与其它的子序列无关；v长的周期长的周期(long period)：v均匀分布的随机数应满足均匀性均匀分布的随机数应满足均匀性(Uniformity)：随机数序列应是均匀的、无偏的，即：如果两个子区间的随机数序列应是均匀的、无偏的，即：如果两个子区间的“面面积积”相等，则落于这两个子区间内的随机数的个数应相等。相等，则落于这两个子区间内的随机数的个数应相等。例如：对例如：对0,1
21、)区间均匀分布的随机数，如果产生了足够多的随机数，区间均匀分布的随机数，如果产生了足够多的随机数，而有一半的随机数落于区间而有一半的随机数落于区间0,0.1不满足均匀性不满足均匀性如果均匀性不满足，则会出现序列中的多组随机数相关的情如果均匀性不满足，则会出现序列中的多组随机数相关的情况况均匀性与互不相关的特性是有联系的均匀性与互不相关的特性是有联系的实际应用中，随机数都是用数学方法计算出来的，这些算法具实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；有周期性，即当序列达到一定长度后会重复；v有效性（有效性（Efficiency):模拟结果可靠模拟结
22、果可靠 模拟产生的样本容量大模拟产生的样本容量大所需的随机数的数量大所需的随机数的数量大随机数的产生必须快速、有效，最好能随机数的产生必须快速、有效，最好能够进行并行计算。够进行并行计算。n随机数与随机数发生器随机数与随机数发生器得到一个可能的随机数序列，是在计算机上实现得到一个可能的随机数序列，是在计算机上实现Monte Carlo方法的关键方法的关键随机数的产生方法：随机数的产生方法：0,1区区间间上上均均匀匀分分布布的的随随机机数数是是Monte Carlo模模拟拟的的基基础础，服服从从任任意意分分布布的的随随机机数数序序列列可可以以用用0,1区区间间均均匀匀分分布布的的随随机机数数序列
23、作适当的变换或舍选后求得。序列作适当的变换或舍选后求得。(0,1)2 -1(-1,1)利用随机数表，如利用随机数表，如Tippett于于1972年发表的随机数表年发表的随机数表；占用太多的计算机内存占用太多的计算机内存采用物理方法，如利用电子线路的热噪声等；采用物理方法，如利用电子线路的热噪声等；昂贵而且不便重复昂贵而且不便重复:伪随机数（伪随机数（Pseudo-Random Number)递推到一定次数后，出现周期性的重复现象。递推到一定次数后，出现周期性的重复现象。利用数学递推公式利用数学递推公式一一旦旦公公式式和和初初值值定定下下来来，整整个个随随机机数数序序列列便便被被确确定定下下来来
24、，而而且且每每一一个个随随机机数数只只被被它它前前面面的的那那个个数数唯唯一一确确定定，因因此此这这类类随随机机数数并不是真正的随机数。并不是真正的随机数。n nMonte CarloMonte Carlo方法基本思想方法基本思想方法基本思想方法基本思想当当所所求求的的问问题题是是某某种种事事件件出出现现的的概概率率，或或是是某某个个随随机机变变量量的的期期望望值值时时，它它们们可可以以通通过过某某种种“随随机机试试验验”的的方方法法，得得到到这这种种事事件件出出现现的的频频率率和和概概率率，或或者者得得到到这这个个随随机机变变量量的的统统计计平平均均值值，并用它们作为问题的解。并用它们作为问
25、题的解。n nMonte CarloMonte Carlo方法解决的问题方法解决的问题方法解决的问题方法解决的问题问题本身是确定性问题，要求我们去寻找一个随机过程，使问题本身是确定性问题，要求我们去寻找一个随机过程，使该随机过程的统计平均就是所求问题的解。该随机过程的统计平均就是所求问题的解。问题本身就是随机过程，我们可以根据问题本身的实际物理问题本身就是随机过程，我们可以根据问题本身的实际物理过程来进行计算机模拟和跟踪，并采用统计方法求得问题的过程来进行计算机模拟和跟踪，并采用统计方法求得问题的解。解。n nMonte CarloMonte Carlo方法的特点方法的特点方法的特点方法的特点
26、计算的收敛性和收敛速度均与问题的维数无关，适合解决计算的收敛性和收敛速度均与问题的维数无关，适合解决高维问题。高维问题。对问题的适应能力强。对问题的适应能力强。收敛速度仅为样本数的收敛速度仅为样本数的-1/2次，因而计算耗时大。次，因而计算耗时大。n nMonte CarloMonte Carlo方法的应用举例：方法的应用举例：方法的应用举例：方法的应用举例：计算积分：计算积分：常用的积分方法求解：常用的积分方法求解：将将积积分分区区域域a,b均均匀匀地地划划分分成成N各各分分区区间间，则则积积分分结结果果可可近近似似地表示成：地表示成：x=(b-a)/N简单的简单的Monte Carlo积分
27、方法求解：积分方法求解：n n利用均匀分布的随机数发生器，从利用均匀分布的随机数发生器，从利用均匀分布的随机数发生器，从利用均匀分布的随机数发生器，从 a,ba,b 区间产生一系列随机区间产生一系列随机区间产生一系列随机区间产生一系列随机数数数数x xi i，i i=1,=1,2,2,.,.,N N其中其中 X为均匀分布，并且为均匀分布，并且 X a,bn n近似求解近似求解近似求解近似求解Eg(XEg(X):):n n近似求解积分近似求解积分近似求解积分近似求解积分:随机抽样随机抽样n n当当当当我我我我们们们们用用用用简简简简单单单单Monte Monte CarloCarlo计计计计算算
28、算算积积积积分分分分时时时时，若若若若该该该该函函函函数数数数为为为为常常常常数数数数函函函函数，数，数，数，g g(x(x)=constant)=constant，则取样数不管多少，准确度为则取样数不管多少，准确度为则取样数不管多少，准确度为则取样数不管多少，准确度为100100。n n如如如如果果果果在在在在积积积积分分分分区区区区间间间间内内内内，g g(x(x)为为为为一一一一平平平平滑滑滑滑函函函函数数数数，则则则则简简简简单单单单Monte Monte CarloCarlo方方方方法法法法较较较较为为为为准准准准确确确确，反反反反之之之之，如如如如果果果果g g(x(x)的的的的变
29、变变变动动动动很很很很剧剧剧剧烈烈烈烈，则则则则简简简简单单单单Monte CarloMonte Carlo方法的误差会变大。方法的误差会变大。方法的误差会变大。方法的误差会变大。说明：说明：说明：说明：重要性重要性Monte Carlo抽样方法抽样方法uu在在在在 g g(x(x)变变变变化化化化剧剧剧剧烈烈烈烈时时时时，如如如如果果果果以以以以Monte Monte CarloCarlo方方方方法法法法取取取取样样样样，最最最最好好好好依据依据依据依据g(xg(x)的大小来决定取样率。的大小来决定取样率。的大小来决定取样率。的大小来决定取样率。uu当当当当|g g(x(x)|)|的的的的值
30、值值值较较较较大大大大时时时时，对对对对 g g(x)dx(x)dx的的的的贡贡贡贡献献献献也也也也较较较较大大大大，如如如如果果果果没没没没被被被被选中，则结果的误差极大。选中，则结果的误差极大。选中，则结果的误差极大。选中，则结果的误差极大。uu解解解解决决决决方方方方式式式式：改改改改变变变变 x x被被被被选选选选中中中中的的的的机机机机率率率率，让让让让|g g(x(x)|)|值值值值较较较较大大大大的的的的点点点点被被被被选中的机率增加。选中的机率增加。选中的机率增加。选中的机率增加。uu采采采采用用用用权权权权重重重重分分分分布布布布函函函函数数数数(Weight Weight
31、distribution distribution functionfunction)w(xw(x):决决决决定每个定每个定每个定每个x x被选中的机率。被选中的机率。被选中的机率。被选中的机率。重要性抽样的定义：根据一定的分布形式进行的随机抽样。重要性抽样的定义：根据一定的分布形式进行的随机抽样。重要性抽样的定义：根据一定的分布形式进行的随机抽样。重要性抽样的定义：根据一定的分布形式进行的随机抽样。uuw w(x(x)必须归一化必须归一化必须归一化必须归一化，即在积分区间内即在积分区间内即在积分区间内即在积分区间内 w(x)dxw(x)dx=1=1。uu由由由由于于于于 x x 的的的的选选
32、选选取取取取已已已已被被被被 w(xw(x)扭扭扭扭曲曲曲曲，所所所所以以以以计计计计算算算算积积积积分分分分时时时时要要要要把把把把这这这这部部部部分分分分 还还还还 回去回去回去回去：若一共取：若一共取：若一共取：若一共取样了样了样了样了N N个个个个x x，则积分值为：则积分值为：则积分值为：则积分值为：重要性重要性Monte Carlo抽样方法抽样方法 Metropolis Monte Carlo方法方法我们所模拟的系统最终要达到的平衡分布是我们所模拟的系统最终要达到的平衡分布是我们所模拟的系统最终要达到的平衡分布是我们所模拟的系统最终要达到的平衡分布是BoltzmanBoltzman
33、分布：分布：分布：分布：Boltzmann Boltzmann 概率分布函数：概率分布函数：概率分布函数：概率分布函数：我我我我们们们们如如如如果果果果能能能能够够够够产产产产生生生生这这这这种种种种分分分分布布布布，我我我我们们们们就就就就能能能能够够够够计计计计算算算算系系系系统统统统的的的的大大大大多多多多数数数数性性性性质质质质，但但但但这这这这是是是是不不不不可可可可能能能能的的的的，因因因因为为为为我我我我们们们们不不不不知知知知道道道道Z Z的的的的值值值值，但但但但是是是是对对对对于于于于任意两个状态，我们有：任意两个状态，我们有：任意两个状态，我们有：任意两个状态，我们有：可
34、可以以在在相相空空间间中中构构造造一一个个马马尔尔科科夫夫链链，使使相相空空间间中中的的样样本本点随着链的增长逐步趋近于点随着链的增长逐步趋近于Boltzman分布。分布。一个序列一个序列x0,x1,x2,xn,如果对任何如果对任何n都有：都有：则此序列是一个则此序列是一个Markov链。链。要求：要求：任任何何一一次次实实验验的的结结果果依依赖赖于于前前一一次次的的试试验验，并并且且近近依依赖赖于于前一次的试验。前一次的试验。马尔科夫马尔科夫(Markov)链链可可以以证证明明：通通过过构构造造Markov链链，体体系系中中最最终终的的平平衡衡分分布布就是就是Boltzman分布分布Metr
35、opolis平衡条件平衡条件(Detailed balance condition)：平衡条件：平衡条件：系系统统处处于于状状态态X的的概概率率正正比比于其于其Boltzman因子：因子：如果如果 是是对对称的称的:Metropolis Monte Carlo方法的算法：方法的算法：a)给出一个初始状态，并计算系统的能量给出一个初始状态，并计算系统的能量Eoldb)随机产生一个新状态，并计算新系统的能量随机产生一个新状态，并计算新系统的能量Enewc)如果如果E(Enew Eold)0,则计算则计算Boltzman因子：因子：e)在在(0,1)区间上区间上产生一个均匀分布的随机数产生一个均匀分
36、布的随机数 ；f)如果如果则接受新状态并回到则接受新状态并回到b)g)否则保留原值并回到否则保留原值并回到b)1.1.正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法(Canonical MC Simulation)(Canonical MC Simulation)具有确定的粒子数具有确定的粒子数N、温度、温度T和体积和体积V对于含对于含对于含对于含N N个粒子的系统，位型（构型）的配分函数：个粒子的系统，位型（构型）的配分函数：个粒子的系统，位型（构型）的配分函数：个粒子的系统，位型（构型）的配分函数：某个特定构型的发生概率为某个特定构型的发生
37、概率为某个特定构型的发生概率为某个特定构型的发生概率为 P PNVTNVT(r(rN N)1.1.正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法正则系综蒙特卡罗模拟方法(Canonical MC Simulation)(Canonical MC Simulation)Monte Carlo模拟中模拟中任一物理量的计算：任一物理量的计算：位型积分位型积分概率密度概率密度系统处于位型系统处于位型rN的概率密度的概率密度a)给出一个初始状态，并计算系统的能量给出一个初始状态，并计算系统的能量Uoldb)随机产生一个新状态，并计算新系统的能量随机产生一个新状态，并计算新系统的
38、能量Unewc)如果如果U(Unew Uold)0,则计算则计算Boltzman因子：因子：e)在在(0,1)区间上区间上产生一个均匀分布的随机数产生一个均匀分布的随机数 ；f)如果如果则接受新状态并回到则接受新状态并回到b)g)否则保留原值并回到否则保留原值并回到b)正则系综正则系综MC模拟算法的组织：模拟算法的组织：正则系综正则系综MC模拟算法的流程图：模拟算法的流程图：给定每个分子的初始位置，给定每个分子的初始位置，ri(0)随机选取一个分子，并随机移动到新的位置随机选取一个分子，并随机移动到新的位置计算移动前后的系统能量变化计算移动前后的系统能量变化U拒绝移动拒绝移动U 0?Exp(U
39、)(0,1)?统计系统的热力学性质及其它物理量统计系统的热力学性质及其它物理量统计性质不变？统计性质不变？打印结果，结束打印结果，结束YesNoNoNo接受移动接受移动YesYes大大约约循循环环107到到108次次Monte Carlo模拟中几个热力学量的计算：模拟中几个热力学量的计算：N个粒子系统中的总势能：个粒子系统中的总势能：假设采用截断势能函数：假设采用截断势能函数：Uc：截断范围内的总势能；：截断范围内的总势能；Ulrc：截断半径外对截断半径外对势能的长程校正势能的长程校正(Long-range correction)对于对于LJ流体：流体：含含N个粒子系统中的压力：个粒子系统中的
40、压力：Wc：截断范围内的总维里项：截断范围内的总维里项(Virial)；Plrc：截断半径外对压力的长程校正：截断半径外对压力的长程校正Read simulationparametersStartInitialize positionsof all particlesNew simulation?Read oldconfigurationMonte CarloloopStopyesnoMonte Carlo loop SubroutineStartStopTrial moveSatisfy Metropolisrule?Accept thetrial moveUpdate energyand
41、virialSample thepressureEnd ofsimulation?yesnoyesnoMain program正则系综正则系综MC模拟程序基本结构：模拟程序基本结构：正则系综正则系综MC模拟程序模拟程序F11讲解（讲解（LJ,NVT）：）：*READ INPUT DATA*初始状态：初始状态：READ(*,(A)Title !运行作业题目运行作业题目READ(*,*)NStep !运行步数运行步数READ(*,*)Iprint !打印步数打印步数READ(*,*)Isave !保存步数保存步数READ(*,*)Iratio !调整步数调整步数READ(*,(A)CNFile !
42、位型文件位型文件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 ！随机移动的最大距离！随机移动的最大距离DensLJ=De
43、nsDens =Dens/(Sigma*3)IF(Rcut.GT.0.5)Stop Cut-Off Too Large 模拟盒子的边长为模拟盒子的边长为1*Read Initial Configuration*Call ReadCN（CNFile）正则系综正则系综MC模拟程序模拟程序F11讲解（讲解（LJ,NVT）：）：初始位型：初始位型：参阅程序参阅程序F23Call FCC需要自己给定所有粒子初始位置需要自己给定所有粒子初始位置u面心立方面心立方(face-centered cubic,FCC)：正则系综正则系综MC模拟程序讲解（模拟程序讲解（LJ,NVT）：）：长程校正：长程校正：Sr3
44、=(Sigma/Rcut)*3Sr9=Sr3*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=(V+Vlrc)/Real(N)Ws=(W+Wlrc)
45、/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=Acatma+1.0 Elseif(Exp(-DeltVb).Gt.Ranf(Dummy)Then V =V+DeltV W =W+DeltW RX(i)=RXInew Acatma=Acatma+1.0
46、 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+RZIJ*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.
47、0最小影像原理最小影像原理算法：算法：Metropolis算法算法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 Of I In The New C
48、onfiguration*Call Energy(RXInew,RYInew,RZInew,I,Rcut,Sigma,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)的计算：的计算：计算表达式：
49、计算表达式：物理意义：物理意义：流流体体中中距距一一个个分分子子为为 r 处处出出现现另另一一个个分分子子的的几几率率密密度度，它它反反映映流体中短程有序的特点流体中短程有序的特点。rr+r在在模模拟拟中中通通常常在在盒盒子子长长度度一一半半的的范范围内，考察围内，考察g(r)随距离的变化。随距离的变化。分子模拟中径向分布函数分子模拟中径向分布函数g(r)的算法：的算法：第一步：计算球壳层间的距离，并初始化一些变量：第一步：计算球壳层间的距离，并初始化一些变量：Ngr:g(r)的统计次数的统计次数Delg:球壳层间的距离球壳层间的距离nhist:球壳层的层数球壳层的层数模拟开始时需要给定：模拟
50、开始时需要给定：nhist 以及统计以及统计g(r)的步数间隔。的步数间隔。分子模拟中径向分布函数分子模拟中径向分布函数g(r)的算法：的算法：第二步：统计第二步：统计g(r):分子模拟中径向分布函数分子模拟中径向分布函数g(r)的算法：的算法：第三步：系综统计平均计算第三步：系综统计平均计算g(r):2.2.巨正则系综蒙特卡罗模拟方法巨正则系综蒙特卡罗模拟方法巨正则系综蒙特卡罗模拟方法巨正则系综蒙特卡罗模拟方法 (Grand Canonical MC Simulation)(Grand Canonical MC Simulation)n n恒定恒定恒定恒定 V,T,V,T,和和和和 ，体系的