蒙地卡羅法與表面科學

蒙地卡羅法與表面科學

---

什麼是蒙地卡羅法（Monte Carlo method）呢？好像是與賭博有很大的關係吧！是的，的確如此。它是利用機率學與統計理論來處理、解決平日可能接觸的科學或工程問題。尤其是用在目前正 蓬勃發展中的表面科學，更可發揮得淋漓盡緻了。所以本文先從簡單的幾個數學例子入門，進而探索到複雜的多原子世界（表面科學）。

如何將分子動力學（molecular dynamics）的理論模擬應用在目前正蓬勃發展中的表面科學上，是極富有挑戰性的研究工作。在本刊去年十二月號＜分子動力模擬與表面科學＞一文，曾做 了簡單的介紹；而這對於表面科學的探索，則有賴於結合物理及化學的理論與概念。現在如欲帶入數學的理念，那「蒙地卡羅法」即是最佳的選擇了。

蒙地卡羅法

這方法是在第二次世界大戰末期時，由當時著名的物理學家和數學家，如馮紐曼（vonNeumann）、鄂拉姆（Ulam）、 費米及梅卓普立斯（Metropolis）所發展出來的。當時這方法是用來研究在可引起核分裂的材料中中子的擴散運動，所以它可說是核子武器發展中的附帶 產物。因為這方法和賭博的機率、亂數（random number；見附註）的產生有關係，因此以最著名的賭場蒙地卡羅命名。但因為統計上「量」的需要，此方法一直到1980年後電腦發展一日千里，才真正廣泛地應用於科學的研究而大出風頭。

首先先舉著名「布丰針」（Buffon needle）問題來說明這方法。這是一位法國自然科學家布丰（comte de Buffon）於1777年提出，如果將一根長度L的針「隨意」丟到一平面上，而這平面上是由很多條等間距d（d＞L）且平行的直線組成（見圖一），則這 根針和其中一條直線相交的機會是多少？答案是2L/πd。

證明也簡單：先求出針長L在與平行線垂直的投影分量L sin A，所以針與直線相交的機是（L∣sinA│）／d。因為針丟擲的可能角度是屬於平均分布，也就是說各種角度都有相同的機會，所以可利用sinA絕對值的 平均值，來求出我們所要的機率P，其計算方法是一簡單的積分

【瀏覽原件】

若我們已知針長L及兩平行線的間距是d，則經由多次此重複此實驗，可得到機率P，而且也可以經由此公式推算出π的正確數值。當然用電腦來模擬、計算，是最合適不過了。所以可寫一小程式來測試此公式，也可經此運算建立起對蒙地卡羅法的初步認識。表一顯示，丟擲次數從1,000增至100,000,000，利用以上公式所求得的π值也隨著次數的增加而更接近正確值（3.1415926…）。

接下來的例子，說明類似射擊的問題。假設有一正方形和一個內接圓（見圖二），我們可想像正方形好似一可能射擊到的範圍，而這 內接圓的面積即代表所擊中的靶面。現在假定正方形邊長是2，則圓形的半徑是1。同樣地，利用電腦產生一組亂數（x,y），這x和y的值都必須是 [－1,1]之間的實數，若x²+y²小於1，則記擊中一次，經由多次重複此實驗，得到一「命中率」（命中數目除以全部射擊數目）Q。此Q值即是

【瀏覽原件】

所以π值即等於命中率乘以4（見表二），π值隨著射擊數目增加而更趨近正確值。

到現在為止，或許已對蒙地卡羅法有了初步的認識。但何謂蒙地卡羅法呢？一般而言，是指藉由亂數產生而完成的計算方法。換言之，也就是利用機率、統計的方法來解決機率或非機率的問題。其應用非常廣泛，包括物理、化學、數學、經濟、工程等，漸漸成為科學研究上不可或缺的有效工具。

高維積分上的應用

現在考慮一維積分的問題（見圖三），其型式如

【瀏覽原件】

其中f(x)是x的任意函數，此積分值即代表曲線下的面積。若利用我們前面已介紹過的方法來處理這問題，則可假設此積分型式可改寫成不連續的加和型式

【瀏覽原件】

其中x_i為分布於a與b之間的隨機數列（亂數），共取N次，然後代入函數f，就如同抽樣檢查一樣，得到不同的N個f(x_i)值，最後除以樣本總數N，得到一平均值【瀏覽原件】，面積（積分值）即是此平均值乘以橫軸上下限的差（b－a）（見圖三）。

這種方法可推廣至任意高維空間的積分。現再以二維空間的積分（雙重積分）為例：

【瀏覽原件】

可以轉換為

【瀏覽原件】

以上欲求x_i分布於a與b之間，先求出一亂數s（0＜s＜1），再乘以（b－a）再加a即可。同理欲求得y_i分布於c與d之間，只要將s乘以（d－c）再加c即可。

雖然在維數較低的空間積分，傳統的數值積分比蒙地卡羅法快速而準確，但是在高維積分時，蒙地卡羅積分法是唯一的選擇。在運算過程中我們唯一要注意的是，盡量減少計算誤差；而此高維積分只要取足樣本數目（N夠大），則結果必然可達到所要求的準確境界，和維數並沒關係。

氣體分子的亂步行為

一般而言，氣體分子的速率是非常的快，在常溫常壓下，空氣的分子大約是以每秒500公尺的速率行走。但問題是空氣分子並非能任意地走得很遠，平均每秒中發生10⁹次 碰撞。而在任一次碰撞後，分子的運動方向會改變，所以在一秒鐘內分子實際只走了0.1公分的直線距離。這整個過程中，可用「亂步」（random walks）來描述，這也就是所謂的擴散（diffusion）運動。所謂「亂步」，即是指運動的粒子，在沒有外力的作用下，任意行走，每走一步後，即忘 記其原先走過的軌跡（喪失記憶），下一步更沒有一定的方向可遵循，就如同一個醉漢走在十字路口上，不知下一步會踏向何方。

位能及作用力的介入

從統計力學上來考慮，若有一系統包括了N個粒子並局限於體積V內，系統的溫度為T，這也就是所謂的「正則系集」（canonicalensemble；見附文）。

以我們有興趣的N個質點的系統為例，一物理量的平均值可寫成

【瀏覽原件】 （1）

式（1）的是全部的作用位能，而各種可能的狀態是以波茲曼 （Boltzmann）分布。從一個粒子的運動著手，先任意給一距離，若新的位置所造成的新組態，其位能比原先的位能還低，則此粒子可移至新的位置。若新 位能比舊位能還高，則將此位能差dU代入波茲曼因子exp（－dU/KT），另外再產生一介於0與1之間的亂數ζ，比較這兩個數值：（1）exp（－dU /KT）＞ζ，則此粒子移至新位置（2）exp（－dU/KT）≦ζ，則此粒子退至舊位置如此的步驟，一直延續下去，則整個系統內粒子如何運動，便十分清 楚顯現出來。

蒙地卡羅法雖然在探討分子系統的組態（configmation）空間上是個非常有用的工具，但為了獲得與時間相關的訊息，則分子動力學的方法是不可或缺的。本刊去年十二月號裡曾介紹過此理論，它只不過是解N個粒子的牛頓運動方程而已，所不同於蒙地卡羅法的地方是將新舊組態的取捨法，代之以與時間自然演變的數值積分法。

分子吸附於固體表面上

分子吸附於一固體的表面層上，它的運動行為完全取決於所有的相互作用位能（也就是全部的作用力），所以如何取得正確的作用位能是解決問題的關鍵所在。一般而言，作用位能包括：

一、吸附原子與吸附原子的作用位能；

二、吸附原子與表面層原子的作用位能；

三、吸附原子與固體內部原子的作用位能；

四、表面層與固體內部原子間的作用位能；

五、表面層間原子互相作用位能；

六、固體內部原子間的作用位能。

由此可見，欲求得正確的總位能值，幾乎是不可能。唯有做一些假設化簡困難度，才有可能估計出與實驗結果相符的數值來。至於固 體內要取多少個原子納入計算才算合理？這就決定於總位能值的差別。我們可從小的範圍（取較少的原子數目）開始計算總位能值，然後一直增加原子的數目（擴大 範圍），直到總位能值不再明顯增加，就可停止。此時的範圍與數目即是我們所要的，這種方法也就叫做「收斂性測驗」（convergence test。此外，我們也要注意，若粒子是雙原子分子的話，分子的內能必須加上（即是加上電子能、振動能和轉動能），才不致有誤差。

蒙地卡羅法用 在吸附原子或分子在表面上的運動，是非常的合適，尤其是擴散運動。若吸附原子或分子非常多，以致於構成吸附層原子三度空間的成長增加，這種過程叫「晶格成 長」（crystal growth，見圖四），而結果稱為「薄膜」（thin film）。這薄膜可包括數十至數千層原子層（atomic layer），也有人稱之巨大薄膜為「超晶格」（superlattice）；目前很多科學家與工程師都專注於這方面的研究。在理論方面，我想要探討這種 現象與機制，蒙地卡羅法則是最佳的電腦數值模擬方法了。

更進一層了解表面科學的奧祕

蒙地卡羅法是一種機率與統計結合在一起的數學方法，若能應用在各種科學與工程方面，所得到的結果將會令人非常滿意。尤其在表面科學急速發展的今日，這種藉由電腦運作的模擬方法就顯得非常重要。

此外，與蒙地卡羅法並行的「分子動力學」理論也是很適當的選擇。所以結合這兩種方法，可對固體表面做最正確、詳盡的模擬；其結果與實驗相互印證後，對於表面科學的奧祕，能更增進一層的了解。當然，筆者對於蒙地卡羅法了解非常有限，而這方法是極為博大精深，應用的範圍與程度亦是非常的廣泛而深入。如何將這方法發揮得淋漓盡緻，則待有興趣的讀者去努力發揚了。

參考資料

1. Metropolis, N., A. W. Rosenbluth, M. N.Rosenbluth and A. H. Teller, 1953, Journal of Chemical Physics, 21：1087.

2. Kalos, M. H. and P. A. Whitlock, 1986,Monte Carlo Methods, Volume I: Basics, John Wiley ＆ Sons, New York.

3. Allen, M. P. and D. J. Tildesley, 1987, Computer Simulation of Liquids, Clarendon Press,Oxford.

4. Ciccotti, G., D. Frenkel and I. R. McDonald, Simulation of Liquids and Solids, Elsevier Science Publishers B. V., New York.

5. Georgiev, N., A. Milchev, M. Paunov and B. Dunweg, 1992, Surface Science, 264：455.

6. Zangwill, A., 1988, Physics at Surface, Cambridge University Press, Cambridge.

7. Polanyi, J. C., R. J. Williams S. F. O'shea, 1991, Journal of Chemical Physics, 94：978.

相克東任職於中央研究院物理研究所

註：所謂「亂數」，也叫「隨機數」，就是具有一定的分布函數（distribution function）且各數間是相互獨立的數列（number sequence）。分布函數種類很多，最基本的是均勻分布。通常若未加指明，則隨機數是指均勻分布的亂數而言。在本文中所舉醉漢的例子及未受外力分子的 擴散運動，都是屬於高斯分布（Gaussian distribution）。早在1939年，統計學家就已利用轉盤來產生十幾萬個亂數。但是為了重複實驗及分析，如何重複產生相同之亂數數列的方法，便 成為統計學上極其重要的事件。這方法必須是有系統性及固定性，並且必須兼顧到所產生出來的數列，能具有良好的隨機性（也可稱之為雜亂性）。現在我們也可利 用一簡單的代數式，將它寫成一小副程式，利用高速電腦於頃刻之間便能產生數以萬計的亂數。有興趣的讀者，可參考數值方法（numerical method）方面的書籍。

梅卓普立斯理論

如果一個系統包括N個質點，並且每一個質點的位置都被標明為（x,y,z）及其動量皆表示為（P_x,P_y,P_z），所以在相空間（phasespace）內的「維數」（dimension）即可定義為d＝（3＋3）N。如此則需要d個數值來標示這N個質點所在空間內的任一點。為了簡化這記號，我們可用一大寫的符號X≡(x₁,x₂,x₃,...,x_d)來代表在這相空間中任一點的位置。

任何一點，我們可選取一試驗數值X_t。而此新的試驗數值是否被接受，則取決於X的分布情形。現在我們即定義P（X）是此分布函數。

假設現在是X_n，則下一步的組態X_n+1，如何決定呢？第一步先建立試驗值：X_t=X_n+δ，其中δ是一微小變化值，而其大小取決於我們所要探討的問題而定。

第二步計算出 r=P(X_t)/P(X_n)的比值如r≧1，則此新試驗值X_t被接受。

X_n+1=X_t=X_n+δ,r≧1

如r＜1，則此新試驗值X_t被接受的機率是r，詳細接受與否，則視第三步而定。

第三步，如果r＜1的情況，我們則要產生一亂數ζ，其值介於[0,1]。所以若ζ＞r，則新的試驗值X_t不被接受；但若是ζ≦1，則新X_t值亦可被接受。