含復(fù)雜近鄰的二維正方格子鍵滲流的蒙特卡羅模擬*

尋之朋 郝大鵬

(中國(guó)礦業(yè)大學(xué)材料與物理學(xué)院,徐州 221116)

1 引言

滲流模型蘊(yùn)含著豐富的物理,是統(tǒng)計(jì)物理學(xué)領(lǐng)域的基礎(chǔ)模型之一[1-6].一方面,研究滲流模型對(duì)理解自然界中的相變和臨界現(xiàn)象、分形、標(biāo)度理論等均有著重要的科學(xué)價(jià)值.另一方面,滲流具有廣泛的應(yīng)用,如流體流經(jīng)多孔物質(zhì)、森林火災(zāi)、流行病的傳播(如最近爆發(fā)的新冠肺炎)等.因此,對(duì)滲流的深入研究具有重要的理論價(jià)值和現(xiàn)實(shí)意義.

建立各種格點(diǎn)滲流模型是研究滲流常用的方法之一.在指定格點(diǎn)模型上,每一個(gè)格點(diǎn)(或鍵)是獨(dú)立的并以一定的概率p被占據(jù)(不被占據(jù)的概率為 1-p),被占據(jù)的格點(diǎn)(或鍵)可以形成毗連的團(tuán)簇.當(dāng)占據(jù)概率p從 0 逐漸增大到某一臨界值時(shí),格點(diǎn)模型上將開始出現(xiàn)大到能夠貫穿整個(gè)系統(tǒng)的團(tuán)簇,并且這時(shí)的系統(tǒng)通常表現(xiàn)出連續(xù)相變的特性,稱這樣的系統(tǒng)為滲流,并且稱此時(shí)的臨界概率為滲流閾值,用pc表示.滲流閾值pc是滲流的核心參數(shù)之一,只有準(zhǔn)確地確定滲流閾值,才能更好地研究其臨界行為以及確定臨界指數(shù).

無(wú)論從解析角度還是采用蒙特卡羅模擬等數(shù)值方法,許多格點(diǎn)模型被廣泛研究以獲得滲流閾值pc以及相應(yīng)的臨界指數(shù).在眾多格點(diǎn)模型中,含復(fù)雜近鄰的格點(diǎn)滲流模型由于多方面的應(yīng)用而成為該領(lǐng)域的熱點(diǎn)研究課題之一.例如,該類模型點(diǎn)滲流可以鏡像到格點(diǎn)模型上擴(kuò)展形狀(如圓形、球形)的吸附滲流過(guò)程[7,8].鍵滲流不僅與小世界網(wǎng)絡(luò)有著密切的聯(lián)系[9],還可以從滲流的角度去研究分析流行病的傳播過(guò)程[10,11].其次對(duì)這類滲流模型的深入研究,既能豐富相關(guān)的滲流理論,又對(duì)研究滲流閾值與格點(diǎn)模型的結(jié)構(gòu),尤其是格點(diǎn)配位數(shù)之間的關(guān)系具有重要的指導(dǎo)作用.此外,這類格點(diǎn)滲流模型的結(jié)構(gòu)介于離散滲流與連續(xù)滲流之間,對(duì)其進(jìn)行深入探討也有助于揭示離散滲流和連續(xù)滲流之間轉(zhuǎn)變的物理機(jī)制.

對(duì)含復(fù)雜近鄰的格點(diǎn)滲流模型的研究可以追溯到由Dalton,Domb 和Sykes[12-14]提出的“等效近鄰模型”.之后,研究人員在該方向開展了廣泛的研究,如菱形結(jié)構(gòu)緊密區(qū)域的長(zhǎng)程點(diǎn)滲流[15],考慮最近鄰和次近鄰的體心立方格子的點(diǎn)滲流和鍵滲流[16],含第四近鄰格點(diǎn)的面心立方格子的點(diǎn)滲流[17]等.d’Iribarne 等[18-20]從圖論的角度系統(tǒng)計(jì)算了長(zhǎng)程關(guān)聯(lián)(最大計(jì)算到了第10 近鄰格點(diǎn))的全部11 種阿基米德格子的點(diǎn)滲流.Malarz 和Galam[21]引入“復(fù)雜近鄰”的概念,即不同近鄰格點(diǎn)的各種組合而并不局限于緊密近鄰格點(diǎn),隨后該課題組對(duì)二維正方格子、三維簡(jiǎn)單立方格子、四維超立方格子的點(diǎn)滲流開展了一系列的研究工作[22-26].最近,基于多種有效算法的發(fā)展,含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流也得到了廣泛的研究[27-30].

滲流閾值pc與格點(diǎn)模型結(jié)構(gòu)(如配位數(shù)z)之間的依賴關(guān)系一直是探討的焦點(diǎn)之一.例如文獻(xiàn)[23]指出一些三維格點(diǎn)模型的點(diǎn)滲流閾值與格點(diǎn)配位數(shù)之間滿足冪律關(guān)系pc～z-a,指數(shù)a=0.790(26) .Galam 和Mauger[31],van der Marck[32],以及其他研究人員也分析研究了許多其他體系,得出類似的冪律關(guān)系[29,30].Domb[12]指出在相同的格點(diǎn)單元形狀下,含擴(kuò)展近鄰的格點(diǎn)模型點(diǎn)滲流在配位數(shù)z較大時(shí)的漸近行為與連續(xù)滲流閾值ηc有關(guān);該論點(diǎn)得到了其他學(xué)者的贊同和推進(jìn)[7,8,20],最近更深入的研究[33]揭示出該漸近關(guān)系為pc～2dηc/z.對(duì)于鍵滲流,該漸近關(guān)系趨于Bethe 格子的結(jié)果pc～1/(z-1);以及當(dāng)z為有限值時(shí)滿足有限尺寸修正zpc-1～a1z-x,其中在二維及三維情況下x=(d-1)/d[34,35].這些理論分析結(jié)果均得到了大量數(shù)值模擬工作的驗(yàn)證[27,28,33,35].

但是,以上漸近行為及有限尺寸修正行為僅對(duì)含緊密近鄰的格點(diǎn)滲流模型是有效的,而對(duì)于含復(fù)雜(非緊密)近鄰格點(diǎn)滲流模型的分析卻遇到了困難.其主要原因是這類格點(diǎn)滲流模型表現(xiàn)出所謂的“簡(jiǎn)并”現(xiàn)象,即一個(gè)配位數(shù)z對(duì)應(yīng)于多個(gè)格點(diǎn)模型結(jié)構(gòu),如圖1 所示的含復(fù)雜近鄰格點(diǎn)的二維正方格子(圖中“?”代表中心格點(diǎn),①代表最近鄰或第1近鄰格點(diǎn),②代表次近鄰或第2 近鄰格點(diǎn),······),“第1 近鄰+第2 近鄰”與“第1 近鄰+第3 近鄰”所對(duì)應(yīng)的配位數(shù)是相同的,均為8;但是這兩種結(jié)構(gòu)的滲流閾值卻是不相同的.因此,對(duì)于這類格點(diǎn)滲流模型,僅僅考慮配位數(shù)一個(gè)參量是不夠的,如何綜合考慮格點(diǎn)模型的拓?fù)浣Y(jié)構(gòu)來(lái)消除“簡(jiǎn)并”是解決該問(wèn)題的關(guān)鍵所在.

圖1 正方格子中心點(diǎn)(“ ? ”)的近鄰格點(diǎn),到第 5 近鄰Fig.1.Neighbors of a central site (“ ? ”) on a square lattice,up to 5 th nearest neighbors.

為了更深入地探討滲流閾值與格點(diǎn)模型的幾何或拓?fù)浣Y(jié)構(gòu)之間的關(guān)聯(lián),本文基于圖論[36]、有機(jī)化學(xué)[37,38],以及最近Malarz[39]關(guān)于點(diǎn)滲流中模型參數(shù)識(shí)別的基本思想,引入新的標(biāo)量參數(shù)ξ=其中zi和ri分別為第i近鄰格點(diǎn)的配位數(shù)及其到中心格點(diǎn)的距離,綜合地考慮格點(diǎn)滲流模型的拓?fù)浣Y(jié)構(gòu)來(lái)消除這種“簡(jiǎn)并”,并將該標(biāo)量參數(shù)應(yīng)用于研究分析含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流.同時(shí),基于一種高效的單團(tuán)簇生長(zhǎng)算法對(duì)含復(fù)雜近鄰的二維正方格子鍵滲流進(jìn)行了大量的蒙特卡羅模擬研究,得到了這些格點(diǎn)滲流模型高精度的鍵滲流閾值.所得到的滲流閾值與新的標(biāo)量參數(shù)很好地滿足冪律關(guān)系pc∝ξ-γ,數(shù)值擬合得指數(shù)γ≈1 .

2 理論基礎(chǔ)

團(tuán)簇的統(tǒng)計(jì)分布是研究滲流的一個(gè)中心物理量,通過(guò)該分布的特性即可確定滲流閾值及相應(yīng)的臨界指數(shù).首先定義團(tuán)簇的尺寸為一個(gè)團(tuán)簇中所含格點(diǎn)(或鍵)的個(gè)數(shù),用s表示.則團(tuán)簇的分布通常用ns(p) 表示,它是指在占據(jù)概率p下,尺寸為s個(gè)格點(diǎn)(或鍵)的團(tuán)簇的數(shù)目.在實(shí)際使用時(shí),將對(duì)團(tuán)簇的分布進(jìn)行歸一化處理.在滲流閾值pc處,ns(pc)滿足關(guān)系式

其中,τ為Fisher 指數(shù),Ω為冪律關(guān)系的一階修正指數(shù).τ和Ω均被認(rèn)為具有普適性,只與系統(tǒng)的維度有關(guān).A0和B0為與具體的模型有關(guān)的兩個(gè)參數(shù),是非普適的.由(1)式可得出模型上每個(gè)格點(diǎn)(或鍵)屬于尺寸大于或等于s的團(tuán)簇的概率為

其中,A1=A0/(τ -2),B1=(τ -2)B0/(τ+Ω-2) .(2)式兩邊同時(shí)乘以sτ-2得

可以看出,當(dāng)s較大時(shí),sτ-2P≥s隨s-Ω的變化將呈線性關(guān)系.利用該線性關(guān)系可以確定滲流閾值,因?yàn)楫?dāng)p/=pc時(shí),關(guān)系是非線性的.

當(dāng)p/=pc時(shí),需要引入標(biāo)度函數(shù).在標(biāo)度極限s →∞和p→pc下,(p-pc)sσ趨于常數(shù),P≥s滿足

式中σ為另一普適指數(shù).對(duì)標(biāo)度函數(shù)f(x) 作泰勒展開

其中,C2=B2f′(0) .這里假設(shè)f(0)=1,將(5)式代入(4)式可得

其中D2=A2C2.(6)式表明,當(dāng)s較大時(shí),sτ-2P≥s在p→pc時(shí)將趨于常數(shù),而在p遠(yuǎn)離pc時(shí)將偏離該常數(shù).這提供了確定滲流閾值的另一個(gè)途徑.

3 數(shù)值模擬結(jié)果與分析

單團(tuán)簇生長(zhǎng)算法的基本規(guī)則[29,30]如下:在格點(diǎn)模型上任取一個(gè)格點(diǎn)作為種子,以該種子為起點(diǎn)開始生長(zhǎng)單獨(dú)的團(tuán)簇.在周期邊界條件下,格點(diǎn)模型上任意一個(gè)格點(diǎn)都可以作為種子,本文在數(shù)值模擬中選取為坐標(biāo)原點(diǎn).從種子格點(diǎn)開始,團(tuán)簇以一定的概率p占據(jù)(不被占據(jù)的概率為 1-p)與之連接的鍵從而向近鄰格點(diǎn)生長(zhǎng).在系統(tǒng)尺寸為L(zhǎng)×L=16384 × 16384 的正方格子上,針對(duì)每一個(gè)格點(diǎn)滲流模型,獨(dú)立地生長(zhǎng)不少于 5×108個(gè)系綜(即獨(dú)立的團(tuán)簇).從統(tǒng)計(jì)上來(lái)講,這些團(tuán)簇具有不同的尺寸,尺寸較小的團(tuán)簇將很快停止生長(zhǎng),而有的團(tuán)簇則可以生長(zhǎng)到很大的尺寸.因此在實(shí)際的數(shù)值模擬中,需要設(shè)定團(tuán)簇生長(zhǎng)的上臨界截?cái)嚅撝党叽鐏?lái)避免出現(xiàn)邊界環(huán)繞,同時(shí)控制程序運(yùn)行時(shí)間,當(dāng)團(tuán)簇的生長(zhǎng)達(dá)到預(yù)設(shè)的上臨界截?cái)嚅撝党叽鐣r(shí)停止生長(zhǎng).然后,將尺寸落在區(qū)間 [ 2n,2n+1-1] (n=0,1,2,···)內(nèi)的團(tuán)簇的數(shù)量記錄到結(jié)果向量的第n個(gè)值中.對(duì)于達(dá)到上臨界截?cái)嚅撝党叽缛栽谏L(zhǎng)的團(tuán)簇,它將被統(tǒng)計(jì)到結(jié)果向量的最后一個(gè)值中.對(duì)于本文所模擬的所有格點(diǎn)模型,上臨界截?cái)嚅撝稻O(shè)置為 216,這也就意味著結(jié)果向量的長(zhǎng)度為17.采用單團(tuán)簇生長(zhǎng)算法在數(shù)值模擬中不需要存儲(chǔ)記錄大量的中間迭代數(shù)據(jù),僅需要對(duì)結(jié)果向量進(jìn)行處理即可,這是該算法優(yōu)于其他方法的一個(gè)方面.此外,新的單團(tuán)簇生長(zhǎng)算法[29,30]在原有的基礎(chǔ)上采用了有效的技術(shù)手段對(duì)算法進(jìn)行改進(jìn),成功地避免了每次單個(gè)團(tuán)簇生長(zhǎng)之后將整個(gè)格點(diǎn)模型信息清除.優(yōu)化的基本思想如下:將每個(gè)格點(diǎn)的初始值設(shè)置為零,對(duì)于第n個(gè)團(tuán)簇(對(duì)應(yīng)于統(tǒng)計(jì)系綜次數(shù)),晶格模型上任何值小于n的格點(diǎn)都被視為未被占據(jù).當(dāng)一個(gè)格點(diǎn)在新團(tuán)簇生長(zhǎng)的過(guò)程中被占據(jù)時(shí),其值賦為n.這個(gè)過(guò)程節(jié)省了大量時(shí)間,因此我們可以模擬非常大的晶格尺寸,并且不必在每個(gè)團(tuán)簇(實(shí)際上許多團(tuán)簇都很小)生長(zhǎng)之后清除整個(gè)晶格的信息,使計(jì)算效率得到了顯著的提高.在文獻(xiàn)[29]中,我們給出了該優(yōu)化算法的程序偽代碼.

對(duì)于普適的標(biāo)度指數(shù),目前已得到二維情況下的解析 解τ=187/91,σ=36/91和Ω=72/91,本文在數(shù)值分析中直接使用這些普適參數(shù)值.于是根據(jù)蒙特卡羅模擬數(shù)據(jù),相應(yīng)的物理量,如sτ-2P≥s,便可易于計(jì)算得到.

本文使用 S Q-a,b,···來(lái)表示含第a近鄰,第b近鄰,······的二維正方格子.基于理論公式(6),圖2為在占據(jù)概率分別為 0.250365 ,0.250367 ,0.250368,0.250369 及 0.250371 時(shí) S Q-1,2 格子 鍵滲流sτ-2P≥s隨sσ的變化曲線.可以看出:當(dāng)團(tuán)簇尺寸s較小時(shí),由于有限尺寸效應(yīng),sτ-2P≥s隨sσ急劇下降;當(dāng)s較大時(shí),sτ-2P≥s隨sσ的變化逐漸表現(xiàn)出線性行為,并且隨著占據(jù)概率p趨于滲流閾值pc,曲線逐漸趨于水平,即趨于常數(shù).根據(jù)圖2 中線性部分的性質(zhì),對(duì)每一個(gè)占據(jù)概率下曲線的線性部分計(jì)算斜率,滲流閾值的中心值可以通過(guò)作圖

來(lái)進(jìn)行確定,如圖2 中的插圖所示.通過(guò)插圖中的線性部分與橫軸的交點(diǎn),可以確定滲流閾值的中心值為pc=0.2503683 .

圖2 SQ-1,2 格子鍵滲流在不同占據(jù)概率p 下sτ-2P≥s隨 s σ 的變化曲線,其中 τ =187/91 ,σ =36/91 .插圖表示主圖中所示結(jié)果線性部分的斜率隨占據(jù)概率p 的變化關(guān)系Fig.2.Plot of s τ-2P≥s vs.s σ with τ =187/91 and σ=36/91 for the bond percolation of the S Q-1,2 lattice under different values of p.The inset indicates the slope of the linear portions of the curves shown in the main figure as a function of p.

另一方面,基于理論公式(3),在占據(jù)概率分別為p=0.250365 ,0.250367 ,0.250368 ,0.250369,以 及 0.250371 時(shí) S Q-1,2 格子鍵滲流sτ-2P≥s隨s-Ω的變化曲線如圖3 所示.由圖3 可以很明顯地看出,當(dāng)占據(jù)概率p偏離滲流閾值pc時(shí),曲線表現(xiàn)出明顯的偏離線性行為;而當(dāng)占據(jù)概率p逐漸趨于滲流閾值pc時(shí),曲線表現(xiàn)出越來(lái)越好的線性行為.于是,通過(guò)圖3 可以確定鍵滲流閾值的范圍為0.250368＜pc＜0.250369.

圖3 SQ-1,2 格子鍵滲流在不同占據(jù)概率p 下sτ-2P≥s隨 s -Ω 的變化曲線,其中 τ =187/91 ,Ω=72/91Fig.3.Plot of s τ-2P≥s vs.s -Ω with τ =187/91 and Ω=72/91 for the bond percolation of the S Q-1,2 lattice under different values of p.

綜上兩種方法,可以確定 S Q-1,2 格子鍵滲流閾值為pc=0.2503683(7),其中括號(hào)里面的數(shù)字表示閾值末位的誤差.對(duì)于本文所模擬的其他格點(diǎn)模型,其分析過(guò)程與上述類似;因此不再贅述每一個(gè)格點(diǎn)模型具體的結(jié)果圖形及分析過(guò)程,而是直接給出鍵滲流閾值的數(shù)值處理結(jié)果,并將結(jié)果總結(jié)在表1 中.

表1 含復(fù)雜近鄰格點(diǎn)的二維正方格子的鍵滲流閾值Table 1.Bond percolation thresholds of square lattice with complex neighborhoods.

表1 第2 列給出了每一個(gè)格點(diǎn)模型的配位數(shù)z.可以看出存在明顯的“簡(jiǎn)并”現(xiàn)象,即許多不同的格點(diǎn)模型具有相同的配位數(shù)z,但鍵滲流閾值卻不相同.圖4 給出了這些格點(diǎn)模型鍵滲流閾值pc隨格點(diǎn)配位數(shù)z變化的對(duì)數(shù)-對(duì)數(shù)曲線,也可以明顯地看出“簡(jiǎn)并”行為.這些結(jié)果說(shuō)明,對(duì)于含復(fù)雜(非緊密)近鄰的格點(diǎn)滲流模型,滲流閾值不僅僅取決于格點(diǎn)配位數(shù),而需要對(duì)格點(diǎn)模型的其他幾何參數(shù)進(jìn)行更深入的分析.

圖4 表1 中所示格點(diǎn)模型 pc 隨z 變化的對(duì)數(shù)-對(duì)數(shù)曲線Fig.4.The log-log plot of pc vs.z for the lattices shown in Table 1.

以上分析表明需要采取有效的技術(shù)手段,引入新的參數(shù)來(lái)消除“簡(jiǎn)并”.最近,Malarz[39]綜合考慮第i近鄰格點(diǎn)的配位數(shù)zi以及其到中心格點(diǎn)距離的平方引入了新的標(biāo)量參數(shù)

來(lái)有效消除“簡(jiǎn)并”,并將其應(yīng)用于分析含復(fù)雜近鄰格點(diǎn)的二維三角格子的點(diǎn)滲流,取得了成功.實(shí)際上,這種消除簡(jiǎn)并的方式在圖論中分析樹結(jié)構(gòu)的拓?fù)洳蛔冃訹36],有機(jī)化學(xué)中分析分子拓?fù)渲笖?shù)[37,38]等方面也有成功的應(yīng)用.本文嘗試采用這種方法來(lái)分析含復(fù)雜近鄰二維正方格子的鍵滲流.首先,根據(jù)(8)式計(jì)算得到每一個(gè)格點(diǎn)模型對(duì)應(yīng)的標(biāo)量指數(shù)ξ,如表1 中第3 列所示.對(duì)應(yīng)于不同近鄰格點(diǎn)的zi和ri2參見(jiàn)表2.從計(jì)算結(jié)果可以看出,引入新的參數(shù)后,在很大程度上消除了“簡(jiǎn)并”現(xiàn)象.然后,作出滲流閾值pc隨參數(shù)ξ變化的對(duì)數(shù)-對(duì)數(shù)曲線,如圖5 所示.結(jié)果顯示,滲流閾值pc隨ξ的變化呈現(xiàn)出很好的冪律關(guān)系pc∝ξ-γ,并且線性擬合得到冪律指數(shù)γ≈1 .可見(jiàn),(8)式消除簡(jiǎn)并的方法也能夠成功地應(yīng)用于含復(fù)雜近鄰格點(diǎn)模型的鍵滲流.

表2 正方格子不同近鄰格點(diǎn)的相關(guān)參數(shù)Table 2.Parameters of different nearest neighbors on square lattice.

圖5 表1 中所示格點(diǎn)模型 pc 隨ξ 變化的對(duì)數(shù)-對(duì)數(shù)曲線Fig.5.The log-log plot of pc vs.ξ for the lattices shown in Table 1.

4 結(jié)論

綜上,為了深入分析滲流閾值與格點(diǎn)模型結(jié)構(gòu)之間的關(guān)系,本文基于高效的單團(tuán)簇生長(zhǎng)算法,對(duì)含復(fù)雜近鄰的二維正方格子鍵滲流進(jìn)行了大量的蒙特卡羅模擬研究.得出主要結(jié)論如下:

1) 單團(tuán)簇生長(zhǎng)算法能夠有效地模擬含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流,并獲得了高精度的鍵滲流閾值.對(duì)于這類格點(diǎn)模型,由于模型中存在交叉(非相交)鍵,這在很大程度上限制了之前許多算法或方法的有效應(yīng)用,而單團(tuán)簇生長(zhǎng)算法則提供了模擬此類格點(diǎn)模型的有效途徑.