计算材料学ppt课件.ppt

计算材料学Computational Materials Science,第二讲 蒙特卡罗（Monte Carlo）方法,主讲：张晖电话：13865606861Email: ,本堂课 主要内容,Monte Carlo模拟发展简介Monte Carlo模拟基本原理Monte Carlo模拟典型算法Monte Carlo模拟典型应用,蒙特卡洛法是什么？,蒙特卡洛(Monte Carlo)方法，是在简单的理论准则基础上，采用反复随即抽样的方法，解决复杂系统的问题。其实质是一种概率和统计的问题。 蒙特卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。,MC的基本思想,MC基本思想很早以前就被人们所发现和利用。17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。但要真正实现随机抽样是很困难的，甚至几乎是不可能的。 高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。,确定性系统,随机性系统,模拟,自然界,Monte-Carlo模拟，即随机模拟（重复“试验”）,重复试验,计算机模拟,Monte Carlo方法：,亦称统计模拟方法，statistical simulation method 利用随机数进行数值模拟的方法,Monte Carlo名字的由来：,是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；,Monte Carlo是摩纳哥（monaco)的首都，该城以赌博闻名,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,Monte Carlo方法简史,简单地介绍一下Monte Carlo方法的发展历史,1、Buffon投针实验：,18世纪，法国数学家Comte de Buffon利用投针实验估计的值,1777年法国科学家布丰提出的一种计算圆周率的方法随机投针法，即著名的布丰投针问题。这一方法的步骤是： 1) 取一张白纸，在上面画上许多条间距为d的平行线。 2) 取一根长度为l（ld） 的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m 3）计算针与直线相交的概率 布丰本人证明了，这个概率是: p=2l/(d) ，为圆周率 ：利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料,实验者 年代 投掷次数 相交次数 圆周率估计值 沃尔夫 1850 5000 2531 3.1596 史密斯 1855 3204 1219 3.1554 德摩根 1880 600 383 3.137 福克斯 1884 1030 489 3.1595 拉泽里尼 1901 3408 1808 3.1415929 赖纳 1925 2520 859 3.1795 布丰投针实验是第一个用几何形式表达概率问题的例子，他首次使用随机实验处理确定性数学问题，为概率论和蒙特卡罗方法的发展起到一定的推动作用。,Monte Carlo方法之随机数的产生 许多计算机系统都有随机数生成函数,F90: call random_seed call random_number(a)2、ISEED=RTC() X=RAN(ISEED)Y=RAN(ISEED）,Matlab: x=rand(N) 产生元素在(0, 1)间随机分布的N*N矩阵 s= rand(state,0) 重设该生成函数到初始状态,注意：上述随机数序列均具周期性，如上页random子程序的周期约230。,实例一、计算值,计算过程：1、构造或描述问题的概率过程2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果计算均值,Monte Carlo方法之典型算法与应用,考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点，若有M个点落于“图形”内，则该“图形”的面积近似为M/N。,用该方法计算的基本思路是：1、根据圆面积的公式：s=R2,当R=1时，S=。2、由于圆的方程是：x2+y2=1(x2为x的平方的意思),因此1/4圆面积为x轴、y轴和上述方程所包围的部分。3、如果在1*1的正方形中均匀地落入随机点，则落入1/4圆中的点的概率就是1/4圆的面积。其4倍，就是圆面积。由于半径为1，该面积的值为的值。,REAL R,R1,R2,PIISEED=RTC()N0=0N=300000DO I=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R1.0)N0=N0+1END DOPI=4.0*N0/NWRITE(*,*)PIEND,MC的优点,MC与传统数学方法相比，具有直观性强，简便易行的优点，该方法能处理一些其他方法无法解决的负责问题，并且容易在计算机上实现，在很大程度上可以代替许多大型、难以实现的复杂实验和社会行为。无污染、无危险、能摆脱实验误差。,Monte Carlo方法利用随机抽样的方法来求解物理问题;,数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;,Monte Carlo方法并非只能用来解决包含随机过程的问题: 例如:用Monte Carlo方法计算定积分. 对这样的问题可将其转换成相关的随机过程, 然后用Monte Carlo方法进行求解,注意以下两点：,实例二 定积分计算,事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。下面考虑一个简单的定积分,! 计算x*2在(0,1)上积分计算过程：1、构造或描述问题的概率过程:产生服从分布f（x）的随机变量Xi( )（i=1,2, ,N）2、从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果计算均值( ),面积的计算,f (x),x,辛普逊方法,I = Sn,11,REAL YY=0N=300000ISEED=RTC()DO I=1,N X=RAN(ISEED) Y=Y+X*2/N END DOWRITE(*,*)YEND,limx2dx (dx 0),MC的应用,自然现象的模拟：,数值分析：,数学问题，求积分，求逆矩阵，解线性代数等,经济学模拟：库存问题，随机服务系统中排队问题,人口问题：人口的出生，传染病的蔓延；乃至动物的生态竞争,金属学：扩散、组织长大、相变过程,蒙特-卡洛模拟的意义,能研究不同边界、不同材料的影响 理论不可能、实验耗费太大 用于实验设计无污染 反应堆防护 核弹爆炸能摆脱实验误差 作理论和实验的桥梁,Monte Carlo模拟的步骤：,根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数； （即构造或描述问题的概率过程）,从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果； 有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是实现MC方法的关键步骤，其中最基本的是（0，1）均匀分布。,对模拟结果进行分析总结，预言物理系统的某些特性。,4. 模拟结果的检验,Monte Carlo方法另一个重要问题：随机数,随机数：由单位矩阵分布中所产生的简单子样称为随机数序列，其中的每一个个体称为随机数。 但真正的随机数的不适合电子计算机上使用，因为它需要很大的存储量。利用某些物理现象可以在电子计算机上产生随机数，需要增添随机数发生器和电路联系等附加设备。,伪随机数：,是有数学递推公式所产生的随机数。（近似的具备随机数的性质。） An+1=T(A);An+1= An+k+1伪随机的优点和缺点：判断伪随机数好坏的方法： 1、它能够有较好的均匀性和独立性； 2、它的费用大小，即指所消耗计算机的时间； 3、容量要求尽可能大。,随机数产生的办法,产生均匀分布随机数的几种方法； （1）物理方法；（2）数学方法。伪随机数产生方法：加同余法乘同余法乘加同余法取中方法逆变换法合成法筛选法 。,关于随机数的几点注意,注1 由于均匀分布的随机数的产生总是采用某个确定的模型进行的，从理论上讲，总会有周期现象出现的。初值确定后，所有随机数也随之确定，并不满足真正随机数的要求。因此通常把由数学方法产生的随机数成为伪随机数。 但其周期又相当长，在实际应用中几乎不可能出现。因此，这种由计算机产生的伪随机数可以当作真正的随机数来处理。注2 应对所产生的伪随机数作各种统计检验，如独立性检验，分布检验，功率谱检验等等。,应用之二 生日问题,MC模拟,假设有n个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50的几率人群中至少有2个人生日相同。,n 理论几率 模拟几率0.117 0.1100.411 0.4120.527 0.5200.706 0.6920.941 0.93650 0.986 0.987,INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC()DO J=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1 DO I=1,365 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR=1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DOEND DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END,MC在材料学领域的应用 随机行走,背景,如，布朗运动最简单、无限制随机行走(Unrestricted randon walk,RW ),平均平方端-端位移:,自然科学和社会生活中很多现象都与随机运动有关,可以模拟的内容？扩散；分子运动；。,如图所示，第i个分子在经过N步随机行走后距原点距离为R，对n个分子每步的位移平方求和后取平均值就得到了所有分子距原点的方均距离：,！Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XX(1:1000,1:1000)REAL XXM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XX(J,I)：X*X ,J:第n个原子,I:第几步跳跃!XXM(I): THE MEAN OF XX WRITE(*,*) “原子个数JMAX，跳动次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第n个原子跳跃X=0DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN0.5)THEN X=X+1ELSE X=X-1 END IF XX(J,I)=X*XEND DOEND DOOPEN(1,FILE=“f:DIF1.DAT)DO I=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX !WRITE(1,*) I, XXM(I)END DOCLOSE(1)END,!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,Y,XY(1:1000,1:1000)REAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第几个原子,I:第几步跳跃!XYM(I): THE MEAN OF XY WRITE(*,*) “原子个数JMAX，跳跃次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几个原子X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THEN x=xy=y-1 END IFIF(RN.LT.0.5.AND.RN.GE.0.25)THEN x=xy=y+1 END IFIF(RN.LT.0.75.AND.RN.GE.0.5)THEN x=x-1y=y END IFIF(RN.GE.0.75)THEN x=x+1y=y END IFXY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=“f:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOCLOSE(1)END,!Monte Carlo Simulation of Two Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREAL XYM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XY(J,I)：X*Y ,J:第n个实验,I:第几步跳跃!XYM(I): THE MEAN OF XY WRITE(*,*) “原子个数JMAX，跳动次数 IMAXREAD(*,*) JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DO J=1,JMAX !第n个原子X=0 !Y=0 !DO I=1,IMAX !第几步跳跃RN=4*RAN(ISEED)+1X=X+XN(RN) Y=Y+YN(RN)XY(J,I)=X*X+Y*YEND DOEND DOOPEN(1,FILE=C:DIF2.DAT)DO I=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I)/JMAX !WRITE(1,*) I, XYM(I)END DOEND ！做三维空间随机行走？,实际环境是复杂的:可变，缺陷，风， 季节，电场等。,空间中有随机性缺陷,不退行走,自回避行走,例如模拟高分子的位形 用随机行走方法模拟高分子位形是用随机行走的轨迹代表高分子的位形，行走过的位置代表的是构成分子的原子或官能团，因此，无限制随机行走忽略了体斥效应。 不退行走就是禁止在每一步行走后立即倒退，可以解决刚走的一步与上一步重叠的问题。但不退行走没有完全解决高分子的体斥效应。 自回避行走就是所有已走过的位置不能再走，这样就完全解决了体斥效应问题。,气体分子的随机行走 假设在一个空旷封闭的房间中心滴上一滴香水，挥发的香水分子随机的与空气中的粒子发生碰撞，最终会扩散到整个房间。这里我们将通过计算机模拟400个分子在二维平面内的随机运动,讨论熵变,现在来讨论一下熵在我们这个气体扩散模型中是什么意思。在刚开始的时候，我们的系统处在一个非常有序的状态：所有的分子都处于原点。这时系统的熵为零，系统不存在任何的无序度。随着时间的推移，分子开始不断的向外扩散，这时系统出现了无序状态，熵开始逐渐增加。,系统的熵和我们预测的一样在随时间增大，但当所有分子布满整个边界以内区域时，系统的熵开始趋于稳定。从这一结果中我们可以了解：当所有的分子随机的布满整个区域时，虽然当我们跟踪某一个确定分子时，它还是在区域内到处乱窜，但每个小区域内分子的密度却不会再变化了。所以，一旦气体分子扩散到整个区域以后，不管我们再等上多少时间，系统的熵都不会再有太大的起伏。换句话说，让系统自动回到开始的状态，即所有分子都在原点的状态，已经不可能了。,第一个程序：画一个圆晶粒,USE MSFLIB !INTEGER XR,YR !在区域中画一个圆！ READ(*,*) XR,YRPARAMETER XR=400,YR=400 INTEGER R,S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，Y0Y0=YR/2R=190 !MIN(X0-10,Y0-10) 圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J) IER=SETPIXEL(I,J)END DOEND DOEND,画晶界,!画一个圆USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNS XN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2 ! 圆心位置X0，Y0Y0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IF IER=SETPIXEL(I,J)END DOEND DOEND,第一性原理根据原子核和电子互相作用的原理及其基本运动规律，运用量子力学原理，从具体要求出发，经过一些近似处理后直接求解薛定谔方程的算法，习惯上称为第一性原理1。第一性原理，英文First Principle，是一个计算物理或计算化学专业名词，广义的第一性原理计算指的是一切基于量子力学原理的计算。我们知道物质由分子组成，分子由原子组成，原子由原子核和电子组成。量子力学计算就是根据原子核和电子的相互作用原理去计算分子结构和分子能量（或离子），然后就能计算物质的各种性质。从头算（ab initio）是狭义的第一性原理计算，它是指不使用经验参数，只用电子质量，光速，质子中子质量等少数实验数据去做量子计算。但是这个计算很慢，所以就加入一些经验参数，可以大大加快计算速度，当然也会不可避免的牺牲计算结果精度。,那为什么使用“第一性原理”这个字眼呢？据说这是，来源于“第一推动力2”这个宗教词汇。第一推动力是牛顿创立的，因为牛顿第一定律说明了物质在不受外力的作用下保持静止或匀速直线运动。如果宇宙诞生之初万事万物应该是静止的，后来却都在运动，是怎么动起来的呢？牛顿相信这是由于上帝推了一把，并且牛顿晚年致力于神学研究。现代科学认为宇宙起源于大爆炸，那么大爆炸也是有原因的吧。所有这些说不清的东西，都归结为宇宙“第一推动力”问题。第一推动一定由某种原理决定。这个可以成为“第一原理”。但是为什么称量子力学计算为第一性原理计算？大概是因为这种计算能够从根本上计算出来分子结构和物质的性质，这样的理论很接近于反映宇宙本质的原理，就称为第一性原理了。第一性原理就是从头计算，不需要任何参数，只需要一些基本的物理常量，就可以得到体系基态的基本性质的原理。,材料的显微组织结构与材料的性能具有密切的联系.而晶粒的平均尺寸是材料性能的一个重要指标,晶粒尺寸的变化对材料的塑性、韧性、硬度、强度、耐磨性等力学性能及机械加工性能具有重大的影响.如超塑性变形中,在温度场和应力场的共同作用下晶粒会发生长大现象.它会大大地降低材料的超塑性性能,使材料无法充满型腔,产生裂痕或发生断裂现象. 因此,掌握和控制材料制备和材料加工中晶粒长大过程一直是学者们致力解决的热点问题,对材料的晶粒长大过程有了较深刻的了解,并且能提出一些控制晶粒长大的有效办法.但这些仍然无法对该问题作出满意的解答.一方面,由于实际材料往往无法满足理论模型的假设条件而导致理论预测与实验结果的偏离;另一方面,由于实验条件的制约,很难获得充足而准确的数据,因此它无法正确反映晶粒长大的特征,近年来,由于计算机技术的突飞猛进,为利用计算机技术来模拟材料的晶粒长大过程研究开辟了一条新途径,成为解决晶粒长大理论模型和实验模型遇到的困难的必由之路. 当前,晶粒长大的计算机模拟的基本方法有两种:一是数值模拟,它是在理论模型建立的晶粒长大速度方程的基础上,利用数值方法对方程进行求解,获得稳态晶粒尺寸分布函数及其相关变量之间的数学关系;二是图像模拟,它是将设计一定几何阵列来构造初始组织,依据晶粒长大的物理和数学规律通过程序操作初始组织发生演化,生动地再现了材料的显微组织的演变过程. 由于数值模拟以理论模型为基础,它严重地依赖于理论模型的合理性.而图像模拟能够综合地反映晶粒长大的动力学和拓扑学关系,是晶粒长大计算机模拟的重点.常用的图像模拟方法又有各种确定性方法和蒙特卡罗方法及元胞自动机等随机性方法.,Finite difference approximation of differential equations,A differential equation can be approximated by a finite difference scheme. For example,Forward time,Backward space,Central space,Forward time-,Central space,FICK 第二定律,Fick 第二定律稳态扩散解,REAL C0(0:1000+1),C(0:1000+1) !C(DISTANCE)C0=0.1C0(0)=0.8 !BOUNDARY CONDITIONC0(1001)=0.1 !BOUNDARY CONDITIONC=C0 OPEN(1,FILE=F:DIF.dat) DO JT=1,50000 ！Time DO IX=1,1000 ! distanceC(IX)=C0(IX)+0.45*(C0(IX+1)+C0(IX-1)-2.0*C0(IX) C0(IX)=C(IX) IF(JT=50000) WRITE(1,*) IX,C(IX) END DO C0=C END DOEND,