蒙特卡洛方法在高分子材料中的应用课件

第六章高分子科学中的MonteCarlo方法

MonteCarlo方法——一个十分独特的名字Monte-Carlo,MonacoMonteCarlo原为地中海沿岸Monaco的一个城市的地名，气候温和，景色怡人，人口不到一万，是世界闻名的大赌场。将MonteCarlo作为一种计算方法的命名固然已经赋予了新的内容。然而，顾名思义，MonteCarlo方法的随机抽样特征在它的命名上得到了反映。第六章高分子科学中的MonteCarlo方法Mon1MC方法的发展归功于核武器早期工作期间LosAlamos（美国国家实验室中子散射研究中心）的一批科学家。vonNeumann,Metropolis,Ulam和Kahn等人在电子计算机上对中子行为进行随机抽样模拟，通过对大量中子行为的观察推断出所要求算的参数。LosAlamos小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MC在各种问题中的应用。学术界一般将Metropolis和Ulam在1949年发表的论文作为MonteCarlo方法诞生的标志。MC方法的发展归功于核武器早期工作期间LosAlamos（26.1MonteCarlo方法的基本思想MonteCarlo方法在数学上称其为随机模拟(randomsimulation)方法、随机抽样(randomsampling)技术或统计试验(statisticaltesting)方法．它的最基本思想是：为了求解数学、物理及化学等问题，建立一个概率模型或随机过程，使它的参数等于问题的解；当所解的问题本身属随机性问题时，则可采用直接模拟法，即根据实际物理情况的概率法则来构造MonteCarlo模型；然后通过对模型或过程的观察抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。在高分子科学中的MonteCarlo模拟主要采用直接模拟方法。6.1MonteCarlo方法的基本思想MonteC3MonteCarlo方法的突出特点是，它的解是由试验得到的，而不是计算出来的。其程序结构简单，解题时受问题条件限制的影响较小，具有广泛的适应性。但不能解决精确度要求很高的问题。蒙特卡洛方法需要大量的随机数，计算量很大，人工计算需耗费大量的时间，利用计算机可大大减少计算时间，增加试验次数以提高计算精度，因此，蒙特卡洛方法的广泛应用与计算机技术的发展是不可分割的。MonteCarlo方法的突出特点是，它的解是由试验得到的4设所要求的量x是随机变量ξ的数学期望E(ξ)，那么用MonteCarlo方法来近似确定x的方法是对ξ进行N次重复抽样，产生相互独立的ξ值的序列ξl，ξ2，…，ξN，并计算其算术平均值：根据Kolmogorov的大数定理则有：即当N充分大时，成立的概率等于1，亦即可以用作为所求量x的估算值。设所要求的量x是随机变量ξ的数学期望E(ξ)，那么用Mont5例6-1用统计试验方法求圆周率π考虑边长为1的正方形，以其一角为圆心和边长为半径，在正方形内画一条1／4圆弧，如图所示。在正方形内等概率地产生n个随机点(xi，yi)，i=l，2，3…，n，设n个随机点中有k个点落在四分之一圆弧内，显然，当n→∞时有以下关系成立：因而，圆周率π的估值为：例6-1用统计试验方法求圆周率π在正方形内等概率地产生n6判断随机点(xi，yi)是否位于圆内的判别式为：用一对(0，1)随机数Ul，U2分别模拟随机变量的取值xi和yi，当时，则计数器k值增1。这个判别式就是蒙特卡洛方法的概率模型。当试验次数n足够大时，所得的估值的精度也随之提高。判断随机点(xi，yi)是否位于圆内的判别式为：用一对(0，7例6-2.蒲丰氏问题ComtedeBuffon(1707-1788)FrenchNeedleexperiment,1777例6-2.蒲丰氏问题ComtedeBuffon(178Buffon投针问题：平面上画很多平行线，间距为a。向此平面投掷长为l（

l＜a）的针,求此针与任一平行线相交的概率p。可以证明求出π值其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。Buffon投针问题：平面上画很多平行线，间距为a。向此平面9一些人进行了实验，其结果列于下表：一些人进行了实验，其结果列于下表：106.2MonteCarlo方法与高分子科学MonteCarlo模拟与高分子科学结下了不解之缘是由于高分子科学本身的特点所决定的，因为在高分子科学中存在着大量可供进行MonteCarlo直接模拟的随机性问题。如：由于聚合反应本身的随机性特点，高分子系综内各个成员之间存在着与其生成机理密切相关的特定分布，即体系中所生成的高分子链并非具有相同的分子量，而是存在着所谓的分子量分布问题；在多元聚合中，多元共聚物不仅具有分子量分布，而且导致了不同种单元在高分子链上的排列问题，即所谓的序列分布；在多官能团的聚合反应中的支化和凝胶化问题；高分子链的热降解和辐射降解等等，无一不是随机性问题。6.2MonteCarlo方法与高分子科学Monte11MonteCarlo方法在现代高分子科学中的应用主要具有以下特征:由于高分子凝聚态物理的发展，高分子体系的MonteCarlo研究从对单链的研究转向对高浓度多链体系的研究。由静态平衡态问题向动态和非平衡态问题发展也是当前高分子MonteCarlo模拟的重要特征。高分子链的分子运动学，尤其是高浓度多链体系的分子运动问题是当前研究的重要方面。人们对共混和嵌段共聚物的界面、高分子和液晶的界面、高分子链的吸附、晶态和非晶态的界面性质和相互扩散问题开展了MonteCarlo模拟研究。高分子MonteCarlo方法的新算法也是值得研究的。MonteCarlo方法在现代高分子科学中的应用主要具有以126.3随机数与伪随机数产生均匀分布随机数的方法可以采用物理方法和数学方法。最简单的产生随机数的物理方法是掷骰子游戏；采用电学噪声的变化也可产生随机数。但物理方法产生随机数的“费用”很高，且速度慢。因此，实际应用的随机数一般均在计算机上采用数学方法来产生。用数学方法产生的随机数一般均采用某种确定性的表达式来实现，因此其并非真正的随机，故通常称其为“伪随机数”。用数学方法产生伪随机数的优点是因为它借助于迭代公式，所以特别适合于计算机。而且其产生的速度快、费用低。目前，多数的计算机均附带有“随机数发生器”。6.3随机数与伪随机数产生均匀分布随机数的方法可以采用物13用数学迭代方法产生的随机数存在两个问题：1、递推公式和初始值a1、a2、…、ak确定后，整个随机数序列便被唯一确定下来。即任意一个随机数被前面的随机数唯一确定了，不满足随机数相互独立的要求。2、既然随机数序列是用递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限多的，因此这样的随机数序列就不可能不出现重复地无限继续下去。这种随机数序列出现周期性的循环现象是与随机数的要求相矛盾的。对第一个问题不能从本质上改变，但只要递推公式选得好随机数的相互独立性是可近似满足；第二个问题，则不是本质的，因为用MonteCarlo方法解任何问题时，所用随机数个数总是有限的，只要保证不超过伪随机数序列出现循环现象的长度即可。用数学迭代方法产生的随机数存在两个问题：1、递推公式和初始值14用数学迭代方法产生随机数均存在周期现象，随着迭代过程的不同，其效果也各不相同。一般满足下列要求的产生方法才可被认为是好的：(1)随机性和统计独立性要好；(2)容易在计算机上实现；(3)省时，存贮量小；(4)伪随机数的周期长。用数学迭代方法产生随机数均存在周期现象，随着迭代过程的不同，15乘同余法乘同余法由Lehmer首先提出。由于采用乘同余法具有在计算机上容易实现、快速等优点，因此乘同余法已被广泛采用。乘同余法的迭代公式为，作为[0，1]区间上均匀分布的伪随机数序列。（给出初始值x0及参数λ、M）当周期很大时，可用乘同余法乘同余法由Lehmer首先提出。由于采用乘同余法具有16一个简单的例子一个简单的例子17上面的例子中，第一个随机数生成器的周期长度是10，而后两个的周期长度只有它的一半。我们自然希望随机数的周期越长越好，这样得到的分布就更接近于真实的均匀分布。上面的例子中，第一个随机数生成器的周期长度是10，而后两个18表：乘同余法的参数及周期

表：乘同余法的参数及周期19MonteCarlo方法的核心就是随机数的使用，因此计算机模拟结果的优劣将强烈地依赖于伪随机数的质量。

伪随机数的均匀性伪随机数的独立性对于已经产生的随机数质量的检验主要是：MonteCarlo方法的核心就是随机数的使用，因此计算机20

伪随机数的均匀性检验可用xn的矩来判别，均匀性好的随机数序列在N→∞时应满足下列要求：一阶矩二阶矩三阶矩四阶矩伪随机数独立性检验一般采用χ2检验。

伪随机数的均匀性检验可用xn的矩来判别，均匀性好的随机数序21随机变量的抽样：前面讨论了[0，1]均匀分布的伪随机数的产生，然而在实际应用中概率分布的形式是多种多样的。并满足：产生[0，1]随机数r，如果条件满足，则认为事件Ai发生。一、从随机事件中抽样：假设随机事件的出现概率分别为Pi

(i＝1，2，…n)。为了对随机事件Ai进行抽样，首先需构造累积概率：随机变量的抽样：并满足：产生[0，1]随机数r，如果条件一、22例6-3.掷骰子点数的抽样掷骰子点数X=n的概率为：选取随机数ξ，如则在等概率的情况下，可使用如下更简单的方法：其中［］表示取整数。例6-3.掷骰子点数的抽样掷骰子点数X=n的概率为：23二、连续型分布的抽样：连续型分布的一般形式如下：这里f(t)为分布的概率密度函数。如果分布函数的反函数存在，则连续型分布的一般抽样方法是通过其反函数直接抽样：这里r是[0，1]均匀分布的随机数，F-1为F(x)的反函数。二、连续型分布的抽样：这里r是[0，1]均匀分布的随机数，F24在［a,b］上均匀分布的分布函数为：例6-4.在［a,b］上均匀分布的抽样在［a,b］上均匀分布的分布函数为：例6-4.在［a,b25其抽样方法为：这里r是[0,1]区间均匀分布的随机数。其抽样方法为：这里r是[0,1]区间均匀分布的随机数。26MonteCarlo方法的估值精度ε与试验次数N的平方根成反比，若精度提高10倍，则试验次数N要增加100倍。收敛速度慢是蒙特卡洛方法的主要缺点。蒙特卡洛方法的精度估算有概率性质，它并不断言精度一定好于ε，而只是表明，所算精度以接近于1的概率不超过某一界限，这是蒙特卡洛方法与其它确定性误差计算的根本区别之处。

MonteCarlo方法的估值精度ε与试验次数N的平方根成27例6-5：中子扩散问题原子核反应堆的壁是铅制的，对中子起屏蔽作用。中子从反应堆内侧进入壁内与铅原子发生碰撞。求出穿透铅壁中子数的百分比，被吸收入铅壁中子数的百分比，以及重新返回反应堆中子数的百分比。入口铅墙（长为3d）d例6-5：中子扩散问题原子核反应堆的壁是铅制的，对中子起屏蔽28解：设壁厚为常量3d，中子是垂直进入壁内的，并设每个中子在壁内每次走过d（平均自由程）才与铅原子碰撞，碰撞后以随机的方向弹射，再走过d的距离，和第二个铅原子碰撞，如此继续下去。最后，有三种情况（1）中子穿透铅壁；（2）被铅壁吸收（假定经过8次碰撞后，没有穿透或返回，则认为被吸收；（3）重新返回反应堆。解：设壁厚为常量3d，中子是垂直进入壁内的，并设每个中子在壁29现在研究对中子运动的模拟：假设一个中子在壁内处于与壁内侧距离为x的位置上与铅原子碰撞，然后以θ角的方向弹射，那么θ是[0，2π]之间的均匀分布的随机数。中子经过弹射后，与壁内侧的距离x变为：x+dcos(2πy)若（1）x>3d

则中子穿透铅壁（2）x<0则中子返回反应堆（3）0≤x≤3d

则继续下一次碰撞，重复这个过程直至中子脱离铅壁或8次碰撞后被吸收为止。对5000个中子进行模拟的结果为：穿透26.3%；吸收22%；返回51.7%现在研究对中子运动的模拟：对5000个中子进行模拟的结果为：306.4MonteCarlo方法在聚合物研究中的应用示例聚乙烯分子结构的模拟共聚物序列分布的MonteCarlo算法高分子无规行走(randomwalks)链的模拟研究高浓度多链体系动力学的“空格扩散算法”6.4MonteCarlo方法在聚合物研究中的应用示31聚乙烯分子结构的模拟聚乙烯分子是由重复的(—CH2

—)单体组成的长分子链，由于碳原子有四个共价键，其空间构型如图所示。Flory根据统计力学理论，导出柔性分子链的非晶态结构取无规线团的构象，各高分子链之间可以相互贯通，它们可以缠结。该无规线团模型在70年代利用中子小角散射技术得到了证实。聚乙烯分子结构的模拟聚乙烯分子是由重复的(—CH2—32聚乙烯分子的空间构型在平面上的投影，可以近似地看成如图所示的结构。模拟程序先定义八个方向，并给出每个方向对应的数值，如图。当分子链段方向为3时，其后面分子链的可能取向方向为2、3或4，它们在聚乙烯中是等概率的。至于下面分子链向哪个方向运动，可由计算机产生的随机数来决定。这样就可模拟出聚乙烯的无规线团状分子结构。聚乙烯分子的空间构型在平面上的投影，可以近似地看成如图所示的33共聚物序列分布的MonteCarlo算法共聚反应的MonteCarlo研究开展得较早，所涉及的主要问题是组成和序列分布问题，其主要目的是通过共聚产物的序列分布来获得单体的活性比和鉴别不同的反应机理。共聚反应的MonteCarlo算法比较简单，因此我们只是简要地介绍其基本算法。共聚物序列分布的MonteCarlo算法共聚反应的Mont34具有末端效应两元共聚反应：末端效应是指只有端点上的单体单元对聚合反应的速率常数有影响。对于两元共聚反应的四种增长反应可记为：由此还可定义活性比，具有末端效应两元共聚反应：由此还可定义活性比，35假定，各速率常数与链长无关，而且引发和终止过程可忽略(一般当高分子的链长很长时均可认为引发和终止过程的影响可忽略)，则由—M1*到—M1*的转变概率为：这里[M1]表示投料浓度，而由—M1*转变为—M2*的转变概率为：相应地有：假定，各速率常数与链长无关，而且引发和终止过程可忽略(一般当36MonteCarlo模拟程序可由如下步骤构成:(1)设增长链的第一个单元为—Mi(i＝1，2)，根据[M1]、[M2]、[—M1*]、[—M2*]的浓度可计算活性比r1，r2和转变概率pij；(2)产生一个单位区间内均匀分布的随机数ξ。(3)因pi1+pi2＝l，故若ξ<pi1则在增长链上加上一个Ml单体，并认为其生成了—M1*

。若认为单体Ml和M2的浓度在增长过程中一直保持恒定，则转回步骤(2)继续进行模拟。但若认为单体浓度是可变的，则由于Mi单体消耗了一个分子故必须重新计算浓度[Mi](i＝l，2)，然后再回到步骤(1)继续进行模拟。MonteCarlo模拟程序可由如下步骤构成:(1)设增37(4)若上式不满足，即ξ＞pi1，则表明发生—M1*到—M2*的增长反应。因此，在增长链上加上一个M2单体，并认为增长链的端基己转变为—M2*

。对于恒定浓度的情况，转回步骤(2)，而对于非恒定浓度的情况，则计算变化后的浓度再转回步骤(1)；(5)上述步骤一直重复，直至达到所需的链长或所需的单体转化率。在模拟过程中可统计各感兴趣的量，如链上Ml和M2单体的组成和序列分布等。必须指出，上述过程只模拟了一根链的情况。为了获得较高的统计精度，可重复多条链后进行平均。(4)若上式不满足，即ξ＞pi1，则表明发生—M1*到—M38采用上述MonteCarlo算法，Motoc等模拟了丙烯酸甲酯(M1)和氯丙烯(M2)的共聚反应。他们设r1=0.08，r2=5.1，模拟所得的共聚物组成和三元组的百分数与实验值的比较见表。结果表明，对于该体系可近似地认为只存在末端效应。采用上述MonteCarlo算法，Motoc等模拟了丙烯酸39MonteCarlo模拟结果与实验结果的比较

MonteCarlo模拟结果与实验结果的比较40高分子无规行走(randomwalks)链的模拟无规行走(randomwalks，简称RW)链模型用来研究柔性高分子链在稀溶液中的大尺度性质。RW链即所谓的自由连接链(freelyjointedchain，简称FJC)。其基本特征是：链中两相邻键的夹角(键角)可任意选择，每个键的内旋转角也可任意取值，链中非直接键接的链单元与链单元之间不存在任何相互作用。高分子无规行走(randomwalks)链的模拟无规行走(41蒙特卡洛方法在高分子材料中的应用ppt课件42以计算均方末端距＜R2＞为例，RW链的简单抽样法计算过程如下：(1)键向量r1的起点放在坐标原点，并令k＝1；(2)对于其后的键向量通过产生在半径为a(通常令a＝1)的球面上均匀分布的随机点，并以其作为rk向量的终点以及rk+1向量的起点；(3)计算末端距向量Rk＝Rk-1+rk；(4)如果k=n，则将令Rk=R，并求R2，如果k＜n，则将k+l替代k，并返回到(2)。对于RW链的最基本特征是：<R2>RW∝n高分子物理：<R2>RW=n

l2以计算均方末端距＜R2＞为例，RW链的简单抽样法计算过程如下43格子链模型格子链(latticechain)模型的基本做法是将空间离散化，即链单元只能取空间中某些人为规定的格点(latticesite)。显然，格子链在细节上与真实链有较大的差别，但高分子链的许多统计性质(大尺度行为)并不依赖于链模型的细节。格子链模型格子链(latticechain)模型的基本做法44方格子模型是把空间离散化为一个立方点阵，即链单元的空间坐标只能在这个点阵空间所定义的格点上取值。为了便于计算，通常格子的边长取为1。因其空间维数的不同，人们给予方格子链以不同的名称。在两维空间里，人们一般称其为方格子链(squarelatticechain)；在三维空间中，称其为立方格子链(cubiclatticechain)。方格子模型是把空间离散化为一个立方点阵，即链单元的空间坐标只45一、RW链的抽样：采用直接抽样法生成RW链的方法十分简单，以两维方格子链为例，可由如下几个步骤构成：(1)将第一个链节固定在坐标原点上，并设格子的边长为1；

(2)产生(0，3)整数序列的随机数；(3)由随机数的数值按预先规定的规则，决定下一个链节所达的格点。设有一链单元的坐标为(x，y)，则由伪随机数的具体数值来决定下一个链单元的坐标位置；一、RW链的抽样：采用直接抽样法生成RW链的方法十分简单，以46当随机数为0时，x+1→x；当随机数为l时，y+1→y；当随机数为2时，x-1→x；当随机数为3时，y-1→y；当随机数为0时，x+1→x；47(4)重复(2)，(3)直至所需链长n，记录所需的结果，诸如最后一个链单元的坐标位置rn等；(5)重复(1)~(4)直至达到所需的高分子链的分子数(样本容量)M，由所生成的链的样本，可计算链构象统计的特征量，如均方末端距。这里，为样本中第l条链的末端距平方，(4)重复(2)，(3)直至所需链长n，记录所需的结果，诸48研究高浓度多链体系动力学的“空格扩散算法”为了使得格子链模型能更为有效地推广到研究高浓度直至‘熔体”的高分子体系的动力学问题，陆建明和杨玉良在Larson等提出的键长涨落模型的基础上提出了适合于研究高浓度多链体系动力学的“空格扩散算法”。按照键长涨落模型，模型链的键长允许取两个数值，即方格子的边长(一般取为1)和格子的对角线(边长的√2倍)。由于键长的可涨落性，因此每个格点的配位数分别为8(两维)和18(三维)。研究高浓度多链体系动力学的“空格扩散算法”为了使得格子链模型49下图给出了典型的微松弛模式和禁阻运动模式。考虑到主要想模拟高浓度多链体系，因此采用空格作为算法的运动主体，具体算法可归结为以下几个步骤。下图给出了典型的微松弛模式和禁阻运动模式。考虑到主要想模拟高50在具有周期边界条件的元胞中规则地按所需浓度排入所需链长的链，设链长为n，链的总数为N条，而少量的空格在排布中尽可能分布均匀。本文元胞大小取为LX×LY=44×44，高分子链长取为n=21，高分子数为N=88。链所占的格子分数