蒙特卡罗方法.ppt

1、蒙特卡罗方法,5.1 基本思想和一般过程,5.1.1 基本思想 蒙特卡罗方法亦称统计模拟方法 利用随机数进行数值模拟的方法 在物理问题中，通常要研究一些随机变量的分布规律，比如分布函数的形式，分布参数的估值等。利用解析方法来求解非常困难，利用实验研究随机变量的分布需要反复进行大量的实验，不仅耗资巨大，而且有些物理量或物理行为根本就难以测量或无法测量。这时就可以根据待求随机问题的变化规律和物理现象本身的统计规律人为地构造出一个合适的概率模型来模拟实际的物理过程，并按照该模型进行大量的统计实验，使它的某些统计参量正好是待求问题的解。这种方法称为蒙特卡罗方法 。,5.1 基本思想和一般过程,Buff

2、on投针实验 1768年，法国数学家Comte de Buffon利用投针实验估计值,长度为 l的针随机地落在相距为dl 的一组水平线之间，求针与线相交的概率？,针相对于平行线的位置可以用一个随机向量表示 随机向量平均分布在区间0,d)0,). 其概率密度函数为1/d. 针与平行线相交的概率为,(5.1),5.1 基本思想和一般过程,5.1.2 马尔科夫(Markov)过程 设一个系统的状态序列为x0,x1,xn,如果对任何n都有概率 p(xn|xn-1 ,x1,x0) = p(xn|xn-1) 则称此序列为马尔科夫序列或马尔科夫链。 设将粒子的位置空间(一维，二维或三维)划分为许多大小相等的

3、网格。如果粒子跳到新的位置后就忘记了原来的格点位置，这就是随机行走。如果粒子跳到新的位置后还记得原来的格点位置，并回避与走过的路径交叉或重复，就叫做自回避行走。,(5.2),5.1 基本思想和一般过程,马尔科夫链是一种随机行走的运动状态。任何一次运动的结果仅仅依赖于前一次运动的结果，因此序列x0, x1, , xn发生的概率可以因式分解为 p(x0 x1 xn-1xn) = p(xn | xn-1) p(x1 | x0)p(x0) 初始状态的绝对概率p(x0)=1。因此将这些条件概率称之为单步跃迁概率，或转移概率。 pij=p(xj | xi) = p(xi xj) 从 MC模拟的应用来说，马

4、尔科夫链的最重要的性质是，存在着各个状态的一个不变分布而最终状态应当遵从某一特定分布。若温度不变，则产生的状态就相应于热平衡分布,(5.3),(5.4),(5.5),5.1 基本思想和一般过程,不变分布定义 若一个概率分布uk满足以下条件 uk0，k=1，2，. 则称概率分布uk是不变的或定常的。其中跃迁概率pij形成一个随机矩阵。,(5.6),(5.7),5.1 基本思想和一般过程,遍历态定义 一个状态如果是非周期性的并且是经久的，具有一个有限的平均返回时间，则称这个状态是遍历态。一个由遍历态组成的马尔科夫链叫做遍历的马尔科夫链。 绝对概率取向定理 一个简约的非周期性链，当且仅当它是一个遍历

5、链时，这个链具有一个不变分布。这时对于一切的k都有uk0，并且不论其初始分布如何，各个绝对概率都趋于uk。 这个条件保证了，不管初始分布如何，马尔科夫链中的状态最终会遵从某个惟一的分布。这个定律在计算物理学的模拟领域中具有基础重要性。,5.1 基本思想和一般过程,5.1.3 蒙特卡罗的一般过程 Monte Carlo的计算过程就是用数学方法在计算机上实现对随机变量的模拟，以求得对问题的近似解。 用蒙特卡罗方法解题的一般过程可归结为以下三个步骤： (1) 构造或描述问题的概率过程 对于本身就只有随机性质的问题，如粒子输运问题，主要是正确地描述和模拟这个概率过程。对于本来不是随机性质的确定性问题，

6、比如计算定积分、解线性方程组及偏微分方程边值问题等，要用蒙特卡罗方法求解，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求的问题的解。,5.1 基本思想和一般过程,(2) 实现从已知概率分布的抽样 有了明确的概率过程后，为了实现过程的数字模拟，必须实现从已知概率分布的随机数的抽样，进行大量的计算随机模拟实验，从中获得随机变量的大量试验值。各种概率模型具有不同的概率分布，因此产生已知概率分布的随机变量，是实现Monte Carlo方法的关键步骤。最简单、最基本、最重要的一个概率分布是(0，1)上的均匀分布 (或称矩形分布)。随机数就是具有这种均匀分布的随机变量。对于其他复杂概率模型的概

7、率分布可以用数学方法在此基础上产生。因此，随机数是Monte Carlo模拟的基本工具。,5.1 基本思想和一般过程,(3) 建立各种统计量的估计，得到问题的求解 一般说来，构造了概率模型并能从中抽样后，即可对所得抽样值集合进行统计处理，从而产生待求数字特征的估计量，给出问题的求解及解的精度估计。,5.2 随机数和伪随机数,由单位矩形分布中所产生的简单子样称为随机数序列，其中的每一个体称为随机数。而所谓伪随机数，则是指用数学递推公式所产生的随机数。由于这种办法属于半经验性质，因此只能近似地具备随机数性质。判断产生伪随机数的某种方法的好坏，首先是它均匀性和独立性是否能较好地具备；其次是它的费用大

8、小，在电子计算机上产生伪随机数，费用即指所用计算机的时间。由蒙特卡罗方法解题一般过程可知，随机数的产生是实现Monte Carlo方法的关键步骤，是Monte Carlo模拟的基本工具。,5.2 随机数和伪随机数,5.2.1 随机数 矩形分布也常称为均匀分布，其中最基本的是单位矩形分布，其分布密度函数如下： 利用某些物理现象可以在电子计算机上产生随机数，但其产生的随机数序列无法重复实现，使程序无法进行复算，结果无法验证。同时需要增添随机数发生器和电路联系等附加设备，费用昂贵。因此，在Monte Carlo方法中一般不采用随机数，而采用伪随机数。,(5.8),5.2 随机数和伪随机数,5.2.2

9、 伪随机数 伪随机数是用数学方法产生的随机数，在给定初值下，由以下的递推公式 确定 (n=1，2，)。 由此产生的随机数并不相互独立，可通过适当地选取递推公式来近似满足独立性要求；另一方面，在电子计算机表示中在(0，1)之间的随机数是有限的，当产生的随机数出现下述情况：,(5.9),(5.10),5.2 随机数和伪随机数,在随机数序列中就出现周期性循环的现象，这与随机数的要求也是相违背的，但Monte Carlo法中计算所用的个数也是有限的，只要其个数不超过随机数产生周期性的个数即可。正是基于这种原因，将数学方法产生的随机数称为伪随机数。,5.3 随机抽样,5.3.1 简单抽样 如果有足够多满

10、足使用要求的伪随机数，就可以有效地运用蒙特卡罗方法给出高维定积分的近似值。由均匀分布中选取随机数的积分方法称为简单抽样蒙特卡罗方法。 考虑一维定积分 其近似值为,(5.11),(5.12),5.3.1 简单抽样,这种积分数值近似解的统计误差以 的比例减小。这就是说，近似值与所求积分真值之间的偏差随n趋于无穷大而变为零。由此看出，积分 的值可以由n个函数值h(x1)、 h(x2) 、h(xn)的平均值近似给出，这里的n是一个很大的数。上述方法就相当于把n个坐标x1、x2xn均匀分布于区间a,b且认为每个点是等权重的，然后在n个离散子集中通过抽样被积函数值而给出积分的近似值，这样一种积分方法又被称

11、为“确定性抽样”。,5.3.1 简单抽样,代替在区间a,b 上等间隔的选取子区间，我们可以在均匀分布于区间a,b 上的m个随机选取的坐标x1、x2xm中进行抽样，然后求平均，从而得到一个与上述所求等效的积分近似值。这种方法要求，在所研究的区间上应给出大量的分别由m个无相关随机数组成的链(即马尔科夫链)。 将被积函数在这些点的响应取值看作一个统计样本集：,(5.13),非权重数值积分方法示意图,确定性抽样,随机抽样,5.3.1 简单抽样,由随机抽样法得到的积分数值近似解的统计误差以 的比例衰减。这种认为自变数在所研究区间均匀分布，并随机选择被积函数值的蒙持卡罗方法就叫做简单抽样法或非权重随机抽样

12、法。 随机抽样法的真正优势表现在对较高维积分的近似求解，诸如在多体动力学和统计力学中所遇到的问题。蒙待卡罗方法对较高维体系的积分误差仍是 ，而这时梯形定则给出的误差变为1/m2/D，这里D为维数。,5.3.1 简单抽样,将其推广到多维的情况 式中(x1i,x2i,xki)=xi 表示从由s个无相关抽样矢量构成的子集中抽取的第i个矢量，其每一个矢量都是由k个随机选取且等权重的矢量分量组成。每一矢量都代表在由所求积分表述的多维物体上的抽样点。,(5.14),5.3.1 简单抽样,原子尺度多体问题的模拟是高维蒙特卡罗积分方法的典型应用领域。然而，必须借助统计力学方法才能从大量的微观数据中获得与宏观性

13、质有关的数据。例如，在离散能量情况下，特性量q的正则系综平均NVT可写为: =r1rN, p1pN表示相空间组态或相空间中的点，NVT ()表示相空间分布函数。若用概率密度wNVT ()表示，则：,(5.15),(5.16),5.3.1 简单抽样,对于蒙特卡罗积分法，采用多种S个离散抽样矢量i，上式可以表示为,(5.17),5.3.2 重要抽样,描述多体系统状态的大多数平均量都可以用积分形式表示。为了选择合适的积分方法，我们必须考虑这些平均量的某些性质。首先，平均量应是多维积分因为它们与N个粒子的独立坐标和动量矢量(对于原子气就有6N个变量)有依赖关系。其次，被积函数值可以变更n个数量级，即明

14、显石同于配分函数中的玻耳兹曼因子。这就是说某些组态对积分有较大的贡献、而另一些组态的重要性则可以忽略。第二，所选用的方法通常要能进行陡变曲线函数的积分，例如在微正则系综中所见到的相空间密度函数。,5.3.2 重要抽样,把任意陡的被积函数变换成非常平滑的函数且调整积分区间的想法是至要抽样法的基本思想。换句话由简单抽样法扩展为重要抽样法，其一个最主要的改进应当是使用了权重被积函数。这就是说，所使用的伪随机数是从非均勾分布中选取的。这种操作方法允许我们把精力集中于在空间区域对函数值的计算与评价，使其对积给出恰当的贡献。引入权重函数g(x)，则对积分J得估算可以写成：,(5.16),5.3.2 重要抽

15、样,把变量x换为y(x) 则积分变为 关于这个积分的蒙特卡罗近似求解前面已讨论过，即随机抽样在区间y(xa)，y(xb)上均匀分布的n个点，然后对相应于h(x)/g(x)的离散值取平均，即：,(5.17),(5.18),5.3.2 重要抽样,即： 权重函数引入到蒙特卡罗方法的过程： 函数h(x)在a和b之间的积分J可以表示为 假定h(x)=exp(-x/2)，根据其级数展开式1-x/2+x3/2-x3/6+1-x/2，则原来的基本表达式可以表示为：,(5.19),(5.20),5.3.2 重要抽样,利用 和 则,(5.21),(5.22),(5.23),(5.24),5.3.2 重要抽样,变量

16、代换修改了积分区间，同时也改变了被积函数的表达形式。因此，对同一积分我们可以通过下列步骤来完成，即：由对y选择平均分布的随机数代替对原变量x的随机值的选取，用求权重函数h(x)g(x)的积分代替对原始函数h(x)的积分。如果y的取值是一个恰当的分布，则对于接近于1的函数h(x)g(x)的积分将会相对简单一些。 以权重函数代替状态函数为基础，Metropolis蒙特卡罗方法将用于处理化学组分恒定的正则系综和化学组分变化的巨正则系综。在这些经典情况下，通常假定其权重函数具有玻耳兹曼标准型的农达形式。 这就体现了非量子粒子的概率分朽指数地依赖于其能量的基本原则。,5.4 Metropolis蒙特卡罗

17、方法,Metropolis蒙特卡罗算法是一种权重方法(或称重要随机抽样方法)。在热平衡情况下，对于用极限概率分布ens()表征的特定统计系综，利用这种方法可生成系统状态。在化学组分恒定的系统(正则、微正则和等温等压系综)中，其概率分布是系统哈密顿量的函数；而在化学组分变化的系统(巨正则系综)中，概率分布则是其化学势的函数。 对于任意的两个系统状态i和j，可以用转移概率ij联系起来，这里ij表示系统从状态i运动到状态j的概率。系统第i个、第j个状态的概率密度函数分别由ens(i) 和 ens(j)给出。,5.4 Metropolis蒙特卡罗方法,对于Metropolis蒙特卡罗解法，转移矩阵的矩

18、阵元可写为： 式中，表示Markov链的对称随机矩阵。由在 ens(i) ens(j)情况下计算矩阵元的定则可知， Metropolis算法,(5.25),5.4 Metropolis蒙特卡罗方法,只需要知道比例系数ens(i) /ens(j)，而没必要知道系统配分函数ZNVT。利用矩阵的对称性有： 和 方程式(5.26)是每一个状态i必定存在于相空间某一区域的自然结论。方程式(5.27)表明了微观可逆性原理。这一点可用来证明马尔科夫链的极限分布等于概率密度分布函数ens() 。,(5.26),(5.27),5.4 Metropolis蒙特卡罗方法,Metropolis算法可由随机地(或系统地

19、)选择并更新试探状态i来执行计算任务，并按照下式抽样随机矩阵： 式中，N是粒子可能的位置数，并根据ens(i) /ens(j) 。计算组态能量值的变化。,(5.28),5.4 Metropolis蒙特卡罗方法,正则和微正则系综的Metropolis方法 在Metropolis蒙特卡罗方法中，若对于系综的N个原子，都安排一个初始位置就可以计算这种组态的哈密顿量。显然，其系综组态变化依赖于所给定的宏观特性值；例如在正则和微正则系统中，其基本主要组分数是给定的。 在蒙特卡罗模型中，通过任意地或系统地改变粒子位置可以给出新的组态。例如，把原子的位置由i变到j就代表了系统状态的变化，并记为ij，在相空间

20、它表示系统由i点运动到j点。根据所考虑的粒子之间的具体相互作用，这一位移将使系统能量由H(i)变为H(j)。然而原子的这一任意位移并不影响系统的组成。,5.4 Metropolis蒙特卡罗方法,在这一变化过程中系统的哈密顿量变化 H(ij) = H(j) - H(i) 如果新的能量值比变化前的能量值小，则该位移使得系统处于较低能量状态。这时，变更即被接受且其位移原子保留在新的位置上。 然而，如果新的能量值大于前一个能量值，则该变更仅以一定概率ens(i) /ens(j)被接受，在正则系统中这个概率由玻耳兹曼因子给出,(5.29),(5.30),5.4 Metropolis蒙特卡罗方法,当H(i

21、j) 取正值时，接受组态变化的概率可表示为 根据Metropolis等人(1953年)的理论，若假定在0和1之间有一个随机数，则新的组态可按照下列定则确定： 如果新的组态被驳回，亦即没有被接受，则把原位置记为新位置并记数一次，然后采用其他随机选定的原子重复上述过程。,(5.31),(5.32),5.4 Metropolis蒙特卡罗方法,巨正则系综，其基本组分的浓度是没有限定的。这样，就可任意选择一个原子并通过改变一次原子种类而得到一个新的组态j。显然，这一过程将影响到系统的化学组成。 同样，组态的变更将以一定的概率被接受。在巨正则系综，因其组成变化而引起的能量变化U(ij) 、也可由计算得到。

22、对新的组态可按照下列标准进行判断：如果得到的新的能量小于原来的能量，则组态的这一变更被接受，如果新能量大于原来的能量，则这一变更将以某一概率被接受。 在巨正则系综，当能量变化U(ij) 为正值时，其组成改变被接受的概率可出下式给出：,(5.33),5.4 Metropolis蒙特卡罗方法,式中， U表示混合能及混合物化学势改变量之和。同理，根据方程式(5.32)确定新的组态是否被接受。如果新组态不被接受，则把原组态记作新的组态，而后再任意地或系统地选择其他原子重复上述过程。,5.5 自旋蒙特卡罗模型,在固体物理学和材料科学等领域，当处理离散空间晶格和局域相互作用等问题时，必须引进考虑其他自由度

23、的模型和方法。在这些方法中，能量变化通常是由于粒子自旋的翻转引起的，而不是由于粒子的位移或粒子交换。这对于通过在晶格格点配置特性粒子而预测多体系统综合特性来说，具有重要意义。 在这些模型中，最简单的就是1/2自旋伊辛(Ising)模型，它可作为磁性材料或二元合金系的一个粗糙的近似方法。q态波茨(Potts)模型是原始伊辛模型的推广或考虑多自旋之后的改进形式。在原始伊辛模型中，自旋带有一定的随意性，不是一个很明确的物理量。在介观尺度的计算材料学领域，基于波茨自旋模型的模拟方法具有特别重要的实用价值。,5.5 自旋蒙特卡罗模型,通过存在于任意晶格上的粒子自旋的随机翻转，以及采用Metropolis

24、蒙特卡罗抽样方法对引起的能量变化加权平均，我们可计算得到这些自旋晶格模型的热力学性质及其演化。 1/2自旋伊辛模型 在1/2自旋伊辛模型中，由在规则晶格格点上的原子或分子之间的对相互作用能求和，就可计算出系统内能。如果用磁学的语言表述，则伊辛模型是由存在相互作用且与外部磁场也有互作用的自旋自由度的离散集构成。 原始的伊辛自旋模型仅考虑最近邻相互作用。若在伊辛模型中考虑到长程相互作用，则有时被称为扩展伊辛模型。,5.5 自旋蒙特卡罗模型,原始的伊辛模型是研究固体中分子磁矩铁磁性有序化的较为恰当的方法。假定在格座i的自族变量Si有两个不同的状态亦即“自旋向上” Si +1或“自旋向下” Si -1

25、。这较好的描述了1/2自旋的情况，尽管自旋被赋予了经典自由度以及没有结出量子角动量的对易规则。将伊辛模型扩展成为一个真正的量子方法，是在海森堡模型中实现的。因为每个格座的自旋可假设只有两个状态，所以伊辛模型适宜于研究具有周期结构的二元台金的原子组态。由此可见，自旋就表示了各个格座的占有情况。,5.5 自旋蒙特卡罗模型,对铁磁有序化的经典伊辛系综，由最近邻相互作用能推导的哈密顿量可表示为： 式中，J为有效相互作用能；B代表某个强度热力学场变量，例如磁场；表示对所有最近邻进行一对一对的求和；Si为自旋变量，其指向为格座i处局域磁矩的方向。 对于含有组分A、B、C的三元合金型伊辛系统，其哈密顿量HA

26、BC可写为：,(5.34),5.5 自旋蒙特卡罗模型,式中， 表示k球内AB对的数目；Vij表示组分i和j的原子之间的相互作用能。对于每一种原子，比如A，其总原子数的表达式为：,(5.35),(5.36),5.5 自旋蒙特卡罗模型,海森堡自旋模型 伊辛模型只能应用于自旋矢量与磁场引起的量子化方向平行或反平行的倩况。这就表明，伊辛模型的哈密顿量只有用于描述在自旋空间具有大的各向异性的磁体才是合效的。 然而，在实际系统中，自族偏离量子化轴的涨落总是在一定程度上必然存在。对此，海森堡提出了一个适合的哈密顿量，即：,(5.37),5.5 自旋蒙特卡罗模型,式中，x1、x2及x3是自旋空间的笛卡尔轴；J

27、i表示各向异性作用能；B代表外场。若J0，则又返回到经典伊辛模型，也就是方程式(5.34)。 海森堡模型和伊辛模型最本质的区别就在于，前者的自旋算符是不对易的，因而，海森堡模型不是一个经典自旋模型，而是一个量子力学自旋模型。,5.5 自旋蒙特卡罗模型,q态波茨自旋模型 在介观尺度关于相变问题的预测研究方面，q态波茨模型具有特别的意义和作用。 波茨模型的基本思想是：代替在伊辛模型中使用的二重自旋变量，而采用广义自旋变量Si，用其表示q个可能状态中的一个态；同时只计不同近邻情况下的相互作用。对此，哈密顿量可由下式给出，,(5.38),5.5 自旋蒙特卡罗模型,由式(5.38)可以看出，对于型哈密顿

28、量当近邻格座上是相同自旋的粒子时其交换相互作用能为零；而当近邻为不同自旋的粒子时则有非零的交换相互作用能。,5.6 蒙特卡罗方法的误差,存在于蒙特卡罗近似方法中的主要误差来源于随机数产生器和统计处理方法。由数字计算机提供的随机数，不但不是真正的互不相关，同时还表现出周期性。因此，在进行大量的使用之前，对所用序列数的周期性进行检验是很有必要的。 对数值蒙特卡罗实验的统计不确定度、可以根据关于大数的中心极限定理来讨论。假设在区间0，1的原始积分J表示为：,(5.39),5.6 蒙特卡罗方法的误差,其随机求积分的近似公式为 由统计学知识可以得到： 式中， 表示h的方差，它是h与其在积分区域平均值之间

29、偏差的一种度量。,(5.40),(5.41),5.6 蒙特卡罗方法的误差,上述式(5.41)揭示了蒙持卡罗近似法中的两个重要的统计问题：第一，误差随着n的增大以n-1/2的速度减小；第二，近似方法的精度随着 的减小而提高。这就是说，当积分函数尽可能的光滑时我们可以获得最佳的积分近似值。 对于低维积分的近似计算，蒙特卡罗积分法的有效性比经典的确定性方法的有效性要差。例如，梯形求积分法对于n步计算所给出的误差为1/n2，这一结果好于由蒙特卡罗近似法给出的误差n-1/2 。然而，蒙特卡罗积分法的主要优势表现在对较高维积分的求解。在较高维体系蒙持卡罗法的积分误差仍然是n-1/2的规模，而梯形方法产生的

30、误差变成了1/n2/D (D为维数)。,5.7 MC方法在材料科学中的应用,晶粒长大是纯金属、合金、陶瓷等多晶体材料中最普遍的现象， 对材料性能具有很重要的影响。 1.晶粒长大过程的模拟 莫春立等采用蒙特卡罗方法对晶粒长大过程进行了模拟。 首先将晶粒结构表示成二维的随机数， 假设随机数在0 到35 之间， 这里的每个数对应着一个晶粒的方位即晶向。以相同晶向并且相邻的区域表示一个晶粒， 在不同晶向的部分即是晶粒边界见晶粒结构图5.1。,图5.1 晶粒结构图,晶粒长大过程的模拟,模拟过程如下： 对计算的区域格点进行剖分， 以确定晶粒结构； 按顺序对每个格点标定序号(从1 到晶粒总数目N ) ； 对

31、每个格点随机地给定一个晶向数目； 计算格点间的相互作用， 每个格点位向的变化导致整体晶粒的长大。具体变化过程如下： a. 随机选取一个初始格点； b. 若此点属于晶界， 那么它可以随机转变为它相邻的位向； c. 计算转变前后的能量变化E (晶向变化产生的能量变化 E )； d. 若E 小于等于零， 则新晶向被接受. 如果E大于零， 则新晶向按一个概率W 被接受，,晶粒长大过程的模拟,如果系统大小为N，那么N个再定向的尝试为1个蒙特卡罗步长。 对于一个小段晶界， 其移动速度v 是由局部自由能G 或化学势来驱动的， 即 在纯金属中， G 和相同，有,(5.42),(5.43),(5.44),晶粒长

32、大过程的模拟,由晶向变化所引起的界面能量变化是晶粒长大的驱动力。对一个晶格点来说， 其界面能与和它相邻的格点的相互作用有关。界面能是晶粒晶向差的函数 式中，Si，Sj为晶粒晶向数，由1Q间选取。矩阵Mij为 其中J 为常数并与晶界能成比例关系， D为delta 函数。,(5.44),(5.45),晶粒长大过程的模拟,可以用二维四边形或六角形点阵进行模拟见图5.2。六角形点阵进行模拟时要考虑周围的六个最紧邻格点，用四边形点阵模拟时要考虑格点周围的4 个最紧邻格点，以及4 个次紧邻格点，并注意边界条件对模拟结果的影响， 为减少边界效应， 一般采用周期性边界条件。晶向数目Q 的选取一般可为16 或3

33、6或更多，Q 过小则晶粒长大不规则。,1,1,1,1,1,1,?,晶粒长大过程的模拟,模拟中采用六角形格点， 其范围为300300 个格点元素， 周期性边界条件， 对材料的晶粒长大过程进行模拟。 微观结构在长大过程中的拓扑学上的变化以及晶粒尺寸d与蒙特卡罗步长MCS 之间的关系见图5.3。可以看到， 不同时刻对应的晶粒尺寸及分布， 即微观组织每一步演化的过程， 这比单纯用平均晶粒直径来表示晶粒长大动力学更直观， 更能揭示演变过程的实质。经多次模拟利用回归方法可以建立模拟平均晶粒尺寸d与蒙特卡罗步长之间关系,(5.46),(a) 100MCS (b) 200MCS (c) 400MCS 图5.3

34、 微观结构与蒙特卡罗步长MCS之间的关系,蒙特卡罗方法,蒙特卡罗方法先将系统用一个哈密顿量描述，并选择一个对问题合适的系综。然后用同这个系综相联系的分布函数和配分函数、就可以计算所有的可观察量了MC方法的基本想法是，对主要的贡献抽样以得到可观察量的估值。 蒙特卡罗方法只给出同系统的位形特性有关的信息，这与分子动力学方法相反，后者导出真实的动力学特性。分子动力学方法给出对时间的依赖关系和系统的坐标和动量等方面的信息。通过选择一个适当的系综，MC方法可以计算在粒子数、体积和温度不变时的可观察量。MC方法一个重大的优点是，很容易实现许多种系综。并且在一次模拟计算过程中还可以改变系综！,蒙特卡罗方法,

35、蒙特卡罗方法 (1)在相空间中规定一个初始点x0. (2)产生一个新状态x (3)计算跃迁概率 w(x，x) (4)产生一个均匀随机数R(0,1) (5)如果跃迁概率W小于随机数R，那么把老状态算作一 个新状态并回到第2步 (6)否则接受新状态并回到第2步。,蒙特卡罗方法,在上面的算法中第1步类似于分子动力学中的初始化。马尔可夫链会失去对初始态的记忆，因此初始状态精确地是什么在很大的程度上并不重要。但是仍然得小心，因为初始状态可能会落到与问题不相干的那一部分相空间之中。在第2步，我们随机地选取一个新状态或新位形例如、在一个Lennand-Jones系统的蒙特卡罗模拟中，将随机地选一个原子并把它

36、随机地移到某一定半径内的另一个位置上。在第3步列第5步要作出决定，是接受还是拒绝一次蒙特卡罗变动为计么我们需要随机数R(0,1)呢?因为概率P(RW(x，x)等于W(x，x)。该算法是任何蒙特卡罗模拟(不论是什么系综)的基础。,蒙特卡罗方法,蒙特卡罗模拟进行的是哪种平均呢?乍看之下也许会认为，MC方法进行的是系综平均。但是现在已经清楚，在一次 MC模拟过程中，进行的实际上是时间平均，与分子动力学中的平均相似。只不过在蒙特卞罗方法中轨道是穿过位形空间，而分子动力学中的轨道则是穿过位置-动量空间。 迄今一直暗中假定，对跃迁概率的选取是满足遍历性限制的。马尔可夫过程中的遍历性指的是，任何一个状态可以

37、从任何其他状态到达。更强的说法是，任何状态必须可以从任何其他状态经过有限步跃迁到达。如果不能满足，那么系统的状态就会分属几个遍历类，各类之间禁止任何跃迁。在刚球高密度堆集的系统中就发生这种情况。这种系统中可能没有从六角密堆集状态到面心立方密堆集状态的跃迁。,蒙特卡罗方法（微正则系统）,考虑一个结定的体积V内具有固定粒子数N的守恒系统，微正则分布由一个函数表示，因此配分函数为 其中E是系统的固定能量。系统的位形只能是其哈密顿量被约束为E的那些位形。 用配分函数计算物理量的方法如下： 对于任何一个可观察量A，有一个依赖于系统状态的函数A(x)与之相联系。可观察量A等于系踪平均,(5.46),(5.

38、47),蒙特卡罗方法（微正则系统）,在微正则系综中，一切状态预先具有相等的权重，如(5.47)式中的函数所示计算上式右端积分的蒙特卡罗方法的总想法是，对系统能到达的那部分相空间抽样然后求和。因此必须建立一个算法，使得系统以遍历方式在相空间中的恒定能量曲面上运动。在微正则系综MC方法中，由于一切系统有相等的权重，系统是按随机行走方式在这个曲面上运动(图5.3)。如果这个随机行走是简单随机行走，而不是一个比方说自回避随机行走，那么就定义了一个马尔可夫过程。,图5.3 相空间中恒定能量曲面上的随机行走,蒙特卡罗方法（微正则系统）,假设通过某种方式产生出一个状态x，它满足H(x)E。一旦处于恒定能量曲

39、面上之后，抽样算法必须能够产生更多仍然在这个曲面上的状态。我们把恒定能量曲面约束放松一点，允许偏离这个曲面值，使得能量处于E- H(x) E+这个区域中。在配分函数中引进一个额外的自由度就可以做到这一点，我们管这个自由度叫一个妖精，它具有能量ED：,(5.48),蒙特卡罗方法（微正则系统）,妖精的作用类似于动能项在分子动力学中所起的作用。它通过在系统附近穿行和输运能量而改变位形，从而使系统在固定能量曲面上产生一次随机行走。不过，我们必须限制妖精的能量，否则它将会吸收掉全部能量！例如，规定妖精的能量只能为正值。,蒙特卡罗方法（微正则系统）,上面概述的过程用算法可表示如下： NVE蒙特卡罗方法 (

40、1)建立一个状态x使得H(x)E. (2)设定妖精的能量ED(例如ED0) (3)选择系统的一个部分 (4)改变系统的局部状态，使得xx (5)计算所产生的能量变化即HH(x)- H(x) (6)如果能量降低了，接受这个改变，今ED ED-H并把x当作新位形回到第3步 (7)否则, 只有当妖精携带足够的能量即ED-H 0时，才接受位形的改变。这时ED ED-H并把x 当作新位形。 (8)回到第3步。,蒙特卡罗方法（微正则系统）,这个算法的第6步和第7步保证系统将弛豫到热平衡。此外，第7步还保证了妖精的能量恒为正。 在概念上我们可以把妖精想象成一只温度计。随着妖精先后与系统的不同部分相接触，它可

41、以保温度计那样取得或交出能量。初始时，妖精有任意的分布。系统起着热库的作用并使妖精热化(thermalize)最终能量服从玻耳兹曼分办，这使得可以计算温度：,(5.49),蒙特卡罗方法（微正则系统）,用Ising模型作为进行蒙特卡罗模拟的一个实际例子。 Ising模型的定义如下：令GLd为一个d维点格。在每个座点i上有一个自旋Si它可以取值+1或-1。各个自旋通过一个交换耦合能J相互作用。此外，还允许存在一个外磁场Q。系统的哈密顿量为 上式等号右边第一项只对最近邻求和。符号代表单个自旋的磁矩如果交换常量J为正，这个哈密顿量是铁磁体的一个模型，即各个自旋倾向于同方向排列若J为负，交换作用是反铁磁

42、性的，自旋倾向于交错反向排列。下面假定相互作用是铁磁性的，J0,(5.50),蒙特卡罗方法（微正则系统）,Ising模型呈现出相变。它有一个临界点Tc在临界点上发生二阶相变。当温度T高于Tc时，序参量即磁化强度m(向上的自族数减去向下的自旋数再除以总自旋个数)在无外磁场时为0。当温度T低于Tc时，存在着一个二重简并的自发磁化。这个模型的相因如图5.4所示。,图5.4 三维Ising模型的相图,M是磁化强度， T是温度， Tc是临界点,蒙特卡罗方法（微正则系统）,为了计算例如一个三维Ising模型的磁化强度。可以用微正则系综蒙特卡罗方法。磁化强度将是能量的函数。但是，有了妖精能量的分布之后，也可得到磁化强度和温度的函数关系。为简单起见，我们令外场为0。 令E是恒定的能量值，并设已建立了具有这一能量值的一个自旋位形s=