流體熵相關(guān)性質(zhì)的Monte Carle摸擬新方法

1、第 4 頁 共 4 頁2001.12.28 流體熵相關(guān)性質(zhì)的Monte Carle摸擬新方法 韓 敏 (無機(jī)化學(xué)專業(yè) MG0124006) 摘要提出了用于計算實際體系的熵性質(zhì)的Monte Carle多級分子取樣摸擬方法.應(yīng)用這一方法,對硬球流體的化學(xué)勢及Helmholtz自由能進(jìn)行了計算,得到了初步的結(jié)果.計算化學(xué)勢時不存在通常實驗粒子方法所遇到的高密度問題.該方法特別適合規(guī)律性的系統(tǒng)研究,較之普通模擬方法要有效得多.模擬得到的硬球體系無限稀釋組分的超額化學(xué)勢與對比直徑的關(guān)系,在相變區(qū)域為一條雙凹曲線;無論是在相變區(qū)還是在單相區(qū),Canahan-Starling公式對這一關(guān)系的描述均有較大偏差

2、.關(guān)鍵詞 Monte Carle 分子模擬 硬球 化學(xué)勢 Helmholtz自由能 熵 計算機(jī)分子模擬是伴隨著計算機(jī)的出現(xiàn)而發(fā)展起來的一門技術(shù),在科學(xué)研究中的地位正日益提高,在一些領(lǐng)域如相平衡的研究中,常常需要計算如自由能熵化學(xué)勢等與體系的配分函數(shù) 性質(zhì)有關(guān)的量.這些性質(zhì)稱為熵相關(guān)性質(zhì).通過分子模擬, 體系的系蹤平均性質(zhì)如內(nèi)能壓力對配分函數(shù)直接作計算可以得到, 但要計算熵相關(guān)性質(zhì),目前倘有一定的難度.許多研究工作者為此作了大量的工作,概括起來有以下幾種方法.第一方法是設(shè)計一個可逆途徑,使所研究的狀態(tài)與以知熵相關(guān)性質(zhì)的狀態(tài)相連,兩個狀態(tài)的自由能之差可由 (A/ V) N,T=-P (1)(A/T

3、)/( (1/T)N,V=E (2) 積分求得1,該方法非常直觀,但需要進(jìn)行多個分子模擬以求得可逆途徑上的壓強(qiáng)或能量,并且如果以可逆途徑通過相變區(qū)域時,這種方法就不適用了.第二種方法是利用“傘形取樣”2或“接收比”3等方法來求算實際體系與某一具有相似勢能的參考體系的自由能之差.這種方法目前已廣泛應(yīng)用于非理想流體混合物性質(zhì)的計算,但若參考體系與實際體系差別較大,仍需通過多個模擬,才能保證計算的精度.第三種方法是直接根據(jù)模擬得到的構(gòu)型采用試驗粒子的方法求算化學(xué)勢4,5,該方法的一個缺陷是難以應(yīng)用于高密度流體,雖然已有一些旨在克服高密度問題的努力,如Mezei提出的孔基巨正則系蹤方法6,但都不能從根

4、本上解決這一問題. 針對分子模擬中熵相關(guān)性質(zhì)計算量大的難題,提出采用多級取樣模擬的方法,使得在單一分子模擬過程中完成所設(shè)計途徑上各狀態(tài)點的模擬,從而求得體系的熵相關(guān)性質(zhì).1.理論 在正則系蹤中,氦姆霍茲自由能與配分函數(shù)有如下關(guān)系: A=-KTlnQn,v,t-(3)可以設(shè)置一個描述體系某一特征的參數(shù)l ,使得l =1時代表實際體系, l =0代表某一參考體系 ,若l在0,1區(qū)間取n個值l i(i=1,2.n)使l i=0, l i+1=1,任意相鄰的參數(shù)l i和l i+1,所描述的兩個體系自由能之差可以表示為:A(l i+1)-A(l i)=-KTlnQn,v,t(l i+1)/Qn,v,t(

5、l i)-(4)將配分函數(shù)的比值表示成系蹤平均的形式,可以得到 A(1)-A(0)=-KTln<exp-b(u(l i+1)u(l i)> l i-(5)上面的<> l i表示在參數(shù)為li時進(jìn)行模擬的系蹤平均,由上式可以看出,實際體系與參考體系自由能之差可以通過l1 . l 2 . l 3 . l n等系蹤平均計算得到。顯然,通過n個分子來計算自由能之差是一項艱苦的工作.為減輕工作量,采用多級取樣模擬的方法,讓參數(shù)l在模擬過程中逐級變化，從而得到各個參數(shù)下體系的系蹤平均.應(yīng)用多級取樣模擬方法,可以推算體系的化學(xué)勢.自由能.混合物的超額熱力學(xué)性質(zhì)等.對化學(xué)勢的計算,可以表

6、示成等容下N+1個粒子體系與N個粒子體系的自由能之差的計算,可通過變化N+1個粒子尺寸的多級取樣模擬完成.對自由能,可采用所有粒子尺寸的多級取樣模擬,得到實際體系與理想氣體或某一已知自由能的參考流體的自由能之差.對于混合物的超額自由能,則可通過NPT正則系蹤里的多級取樣模擬進(jìn)行計算.2.模擬過程本文對硬球流體的化學(xué)勢和亥姆霍茲自由能進(jìn)行了計算.模擬的粒子數(shù)為54,采用立方周期性邊界條件,初始構(gòu)型為面心立方晶格,在開始多級取樣模擬前經(jīng)2-100000構(gòu)型達(dá)到平衡.下面對化學(xué)勢和自由能的計算作簡要的說明.2.1化學(xué)勢的模擬對于硬球體系化學(xué)勢的模擬,參數(shù)可選為一個粒子與其他粒子的直徑比,即l =-d

7、/ d0,采用多級取樣法計算當(dāng)該粒子直徑為d和0時體系自由能之差.該粒子直徑為0時的粒子體系與等容下N-1個粒子體系的自由能之差可嚴(yán)格由 AN(0)-AN-1=Ukl-KTln(1-rpd03/6)-(6) 得到.這里, AN(0), AN-1分別表示中心粒子直徑為 0的N粒子體系和N-1個離子體系的亥姆霍茲體系的自由能.因此,硬球流體 的超額化學(xué)勢為: mex = m-mid=AN(1)-AN(0)-KTln(1-rpd03/6)-(7)再求算化學(xué)勢的多級取樣模擬中,本文人為地選定標(biāo)號為1的粒子為中心粒子,并 逐級變化其直徑.由于計算自由能之差的被積函數(shù)只于中心粒子的作用有關(guān),為提高該函數(shù)系

8、蹤平均精度,采用歐圍齊的優(yōu)先取樣法7,使體系中各粒子 的選擇移動幾率和其與中心粒子的距離有關(guān),本文定義的幾率函數(shù)形式為: W(ri)=1/(ri2+d02)-(8)則粒子的移動幾率為 WI=W(ri)/W(rj)-(9)其中為粒子與中心粒子的距離.對于硬球流體,計算自由能之差的被積函數(shù)可簡化為: 1 minri>li+1d0 EXP-b(U (li+1)-U(li)= -(10) 0 minri/>li+1d0若定義為在參數(shù)的這級模擬中各粒子與中心粒子的最近距離大于的構(gòu)型數(shù)可得(EXP-b(U(li+1)-U(li) li=NK(i+1)/NK(i)-(11)對于硬球體系,在每一級

9、的取樣模擬中,若參數(shù)為,理論上可以計算得到所有的,則和所描述體系的自由能之差,可以利用所有的各級模擬所獲得的結(jié)果,按式(1 2)計算得出,以提高統(tǒng)計精度.AN(lI+1)-AN(lI)=-KTlnNK(i+1)/NK(i)-(12)3. 2亥姆霍茲自由能的模擬硬球體系亥姆霍茲自由能的模擬與化學(xué)勢的模擬不一致.不同之處為在這里代表第i級模擬中粒子之間的最近距離大于ljd0的構(gòu)型數(shù).另外,在模擬中無須采用owicki優(yōu)先取樣法.在l=0時體系的自由能即為理想氣體的自由能,由此即可得到各個值時的自由能.自由能是一種廣延性質(zhì),與體系大小成正比.因此,為保證一定的統(tǒng)計精度,體系越大,在間劃分的級數(shù)就越多

10、.若總的模擬構(gòu)型一定,則各級的模擬構(gòu)型就越少,相應(yīng)地就會造成大的誤差漲落.為獲得實用的模擬結(jié)果,需要合適的數(shù)據(jù)處理方法.本文采用三次平滑技術(shù),得到了較光滑的曲線.本文模擬了對比密度r*=rd3在0-1.2范圍內(nèi)的亥姆霍茲自由能,d0即為r*=1.2時的粒子直徑,在范圍內(nèi)取值,使2的取值間隔為0.001d02.在每級模擬中產(chǎn)生1k個構(gòu)型,其中前1/3的構(gòu)型用于平衡而舍去,剩余構(gòu)型用于計算.由于l在變化范圍中存在相變區(qū)域1,為保證模擬的順利進(jìn)行,本文使由1逐級減小至0的途徑進(jìn)行模擬.3.結(jié)果與討論3.1 化學(xué)勢的模擬化學(xué)勢的 大小如何,總能選擇適當(dāng)?shù)募墧?shù)n以獲得足夠精確的結(jié)果 .圖1表示了本文模擬

11、的結(jié)果,與近似結(jié)果8進(jìn)行 了比較,偏差較大.如果增大模擬的粒子數(shù)和構(gòu)型數(shù),可使精度進(jìn)一步提高.在本文所模擬 的 密度范圍r*=0.851.00內(nèi),通常所用的粒子實驗法已難以應(yīng)用,而采用本法可得到較合理的結(jié)果。圖1中在相變區(qū)(r*>0.944)隨密度的增加，可以看出硬球流體化學(xué)勢先增大后減小的趨勢。Pade近似在該區(qū)域已不適用，故與模擬結(jié)果有較大偏差。 多級取樣模擬的一個優(yōu)越之處是可以同時求得所設(shè)計途徑上各狀態(tài)點的熵相關(guān)性質(zhì)，以便對其作系統(tǒng)的研究。在求算圖1中各密度下的化學(xué)勢的過程中，同時可得到對應(yīng)密度下半徑為的無限稀釋組分的化學(xué)勢。圖2表示了本文模擬的無限稀釋組分超額化學(xué)勢與對比直徑的關(guān)

12、系曲線，同時還給出了相關(guān)密度下由CarnahanStarling公式9計算的結(jié)果。在相變區(qū)，模擬的結(jié)果表明，無限稀釋組分超額化學(xué)勢與對比直徑的關(guān)系為一條具有雙拐點的雙凹曲線，這一特征是CarnahanStarling公式無法描述的，由此可知，根據(jù)這類曲線的形狀可以鑒定體系是否分相。3. 2亥姆霍茲自由能圖3是模擬得到的硬球流體的氦姆霍茲自由能與密度的關(guān)系以及由公式計算得到的結(jié)果?？梢钥闯觯?.90的密度范圍內(nèi)多級取樣模擬結(jié)果與公式的計算結(jié)果相差不大。若在多級取樣模擬中增加各級構(gòu)型數(shù)，就會得到精確的結(jié)果。圖3還給出了Frankel計算得到的硬球面心立方晶格的自由能10，與本文模擬的結(jié)果相近。圖

13、2，3所得的結(jié)果在多級取樣模擬中通過一次普通長度的模擬就可得到。因此，對有規(guī)律性系統(tǒng)的研究比普通模擬方法要有效得多，一方面是由于在多級取樣模擬中，相鄰級的平衡構(gòu)型相差較小，逐級變化參數(shù)時，體系可達(dá)到平衡；另一方面是由于不同參數(shù)下平衡構(gòu)型所包含的信息得到了充分利用，從而提高了所求性質(zhì)的統(tǒng)計精度。 參 考 文 獻(xiàn)1 Hoover W G,Ree F H.J. Chem.Phys.1967,47:48732 Torrie G M,Valleall J.P. J.Comput.Phys.1997,23:1873 Bennett C H. Comput.Phys.1987,22:2454 Fixman M.J. Chem.Phys.1983,78:42235 Deitrick G L. J.Chem.Phys.1989,90:2