第六章高分子科學中的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。鏈所占的格子分數