馬 穎

(湘潭大學材料與光電物理學院，低維材料及其應(yīng)用技術(shù)教育部重點實驗室，湘潭 411105)

非晶態(tài)石英的變電荷分子動力學模擬*

馬 穎

(湘潭大學材料與光電物理學院，低維材料及其應(yīng)用技術(shù)教育部重點實驗室，湘潭 411105)

(2010年1月13日收到;2010年4月15日收到修改稿)

基于迭代變電荷方法，用分子動力學模擬了非晶態(tài)石英的結(jié)構(gòu)與振動特征.首先利用熔化-猝火方法得到了非晶石英的平衡結(jié)構(gòu).在此基礎(chǔ)上獲得了體系不同原子對之間的對關(guān)聯(lián)函數(shù)、鍵角分布函數(shù)和振動頻譜等，結(jié)果與實驗數(shù)據(jù)均符合較好.變電荷方法的計算結(jié)果表明，非晶石英體系內(nèi)粒子的電荷與石英晶體內(nèi)粒子電荷顯著不同，并且出現(xiàn)了較大的漲落.

分子動力學，變電荷，非晶石英

PACS:61.43.Fs，61.43.Bn，63.50.Lm

1.引 言

石英(SiO2)也許是人類歷史上認識最早、應(yīng)用最廣泛的材料之一.常溫常壓下的石英晶體，即 α-石英，屬于三方晶系，在現(xiàn)代工業(yè)中具有廣泛的用途.隨著溫度、壓強等條件的變化，石英還能存在于多種不同的晶體結(jié)構(gòu)中.非晶態(tài)石英，俗稱玻璃，同樣也是重要的工業(yè)材料，具有優(yōu)良的光學、機械、熱穩(wěn)定性能和良好的抗輻照性能，廣泛用于微電子、化工、光學等領(lǐng)域［1］.因此，人們從實驗和理論兩方面對晶態(tài)和非晶態(tài)石英展開了大量的研究，從微觀機制上理解其宏觀性質(zhì)［2—5］.

近年來，隨著計算機軟硬件的飛速發(fā)展，計算機模擬的手段日益成為材料科學研究中的重要工具.現(xiàn)有的實驗手段還很難得到材料在原子尺度上的直觀圖像，而計算機模擬的方法則無疑能有效地彌補這一缺陷.廣泛采用的模擬方法包括第一性原理計算、分子動力學模擬和蒙特卡羅方法等.例如，Liu等［6］運用第一性原理方法計算了各種SiO2晶體的穩(wěn)定性.與第一性原理方法相比，分子動力學不但能夠提供材料在原子尺度上的結(jié)構(gòu)信息，而且由于不考慮體系的電子結(jié)構(gòu)，計算效率要高得多，因此也被廣泛地用于研究晶態(tài)和非晶態(tài)石英.例如，分子動力學模擬表明，石英晶體的 α—β相變?yōu)橛行颉獰o序型相變［7］.非晶態(tài)石英的結(jié)構(gòu)、表面、振動譜等也能由分子動力學模擬描述［8—10］.

遺憾的是，目前對晶態(tài)或非晶態(tài)石英的分子動力學模擬大都采用了所謂的固定電荷假設(shè)，亦即模擬體系內(nèi)粒子所帶的電荷在模擬過程中保持不變，且同種粒子所帶電荷相同.然而，理論和實驗證據(jù)都表明，粒子所帶電荷與粒子所處的環(huán)境密切相關(guān)［11］.如果所模擬的體系為晶體，那么在通常條件下，晶體內(nèi)原子在平衡位置附近做近簡諧振動.因此，原子周圍環(huán)境基本保持不變，且同類原子所處環(huán)境基本相同，固定電荷假設(shè)仍然適用［12］.對于非晶態(tài)物質(zhì)而言，體系內(nèi)不再存在長程有序性，因此不同位置的同類原子所處的環(huán)境不再相同.在模擬過程中，原子所處的環(huán)境也有可能發(fā)生較大的改變.可見，對于非晶體系，固定電荷假設(shè)不再適用.此外，對于缺陷、表面、界面等非周期性體系，固定電荷假設(shè)同樣不再適用［13］.為了模擬上述體系，人們發(fā)展了包括電荷平衡(charge equilibration，Qeq)法［14］、電荷漲落(fluctuation charge，F(xiàn)Q)法［15］、迭代變電荷 (iterative fluctuation charge，IFC)法［16］等的變電荷分子動力學方法.這些方法已經(jīng)在對分子團［17］、晶界［18］等非周期性體系的分子動力學模擬中得到了成功的應(yīng)用.對非晶石英的變電荷分子動力學的系統(tǒng)研究則尚未見報道.本文利用迭代變電荷方法對非晶態(tài)石英進行了分子動力學模擬.模擬得到的非晶態(tài)石英的結(jié)構(gòu)性質(zhì)等與實驗結(jié)果符合較好.特別是對體系粒子電荷的分析表明，在非晶態(tài)石英中粒子電荷與晶態(tài)石英顯著不同，且在模擬過程中粒子電荷有較大的漲落.

2.模擬方法

分子動力學模擬的關(guān)鍵在于勢函數(shù).對SiO2體系，較早采用的勢函數(shù)包括 Born-Mayer-Huggins (BMH)勢和修正的BMH勢等.它們都能較好地模擬非晶石英的結(jié)構(gòu)特征［8］.在 BMH勢和修正的BMH勢中，Si原子和O原子的電荷分別取為+4和-2，亦即認為 Si—O鍵是完全離子性的.事實上，Si—O鍵還有部分共價性.第一性原理計算表明，在原子團中，每個 Si—O鍵之間的電荷轉(zhuǎn)移大約為0.7 e［19］.由于Si和O的配位數(shù)分別為4和2，對應(yīng)的離子電荷則分別為 +2.8和 -1.4.當然，在晶態(tài)或非晶態(tài)SiO2塊體材料中，Si和O的電荷與在原子團中的電荷不同.此外，在分子動力學模擬中，往往還適當調(diào)節(jié)Si和O的電荷以得到最佳的模擬結(jié)果.一般而言，在塊體(晶態(tài)或非晶態(tài))SiO2中，Si和O的電荷分別位于 +1.4—+2.8和 -0.7—-1.4之間.比如 Tsuneyuki等［19］提出的 TTAM勢，van Beest等［20］提出的BKS勢等，都屬于這一類的勢函數(shù).然而，BMH勢，TTAM勢和BKS勢等都有一個明顯的缺陷，就是將帶電離子作為點電荷來處理，利用點電荷間的庫侖定律直接計算二者的相互作用.事實上，由量子力學基本原理可知，用“電子云”來描述電子繞核運動更為合適;而兩帶電離子間的相互作用，則應(yīng)該計算對應(yīng)的電子云密度的庫侖積分［17］.這顯然會導(dǎo)致計算量的大幅增加.最近，Ma和 Garofalini［7］提出了一種新的勢函數(shù)，利用經(jīng)驗公式來計算庫侖積分，從而較好地解決了這一問題，并且分子動力學模擬結(jié)果表明，該勢函數(shù)能更精確地給出SiO2的晶體結(jié)構(gòu)、彈性常數(shù)和聲子色散譜，并成功地模擬了石英晶體的 α—β相變.因此，在本文中，我們采用了 Ma和 Garofalini［7］提出的勢函數(shù)，該勢函數(shù)的具體形式及參數(shù)詳見文獻［7］.此外，對非晶態(tài)材料的分子動力學模擬，有必要采用變電荷方法.在本文中，我們采用了所謂的迭代變電荷方法來得到隨環(huán)境而變化的離子電荷［16］.比之于其他變電荷方法，該方法有計算效率高、程序簡便等優(yōu)點，詳見文獻［16］.為了得到非晶態(tài)石英的平衡結(jié)構(gòu)，我們采取了常用的熔化-淬火法［8］.首先，隨機生成尺寸為3 nm×3 nm×3 nm，凈電荷為零的SiO2體系.體系初始密度為2.2 g/cm3，即石英玻璃密度的實驗值.由于初始結(jié)構(gòu)中Si原子和O原子的位置是隨機生成，并不能代表非晶態(tài)石英的真實結(jié)構(gòu)，因此將體系升溫至5000 K，并保持30 ps以保證初始結(jié)構(gòu)成為完全無序的非晶態(tài).隨后，以降溫速率為500 K/20 ps將體系溫度降至300 K.在降溫過程中，體系溫度由Nose-Hoover熱浴控制［21，22］.當體系溫度達到300 K時，采用了 Parrinello-Rahman的恒壓分子動力學方法［23］，體系壓強維持在 1 atm (101.325 kPa)，模擬時間為30 ps，從而得到非晶態(tài)石英的最終平衡結(jié)構(gòu).在整個熔化-淬火模擬過程中時間步長均取為1 fs.

3.結(jié)果與討論

為了驗證模擬結(jié)果的正確性，我們首先考察了平衡態(tài)非晶石英的結(jié)構(gòu)性質(zhì).圖1給出了在300 K下表征非晶態(tài)石英結(jié)構(gòu)的Si—O，O—O，Si—Si對關(guān)聯(lián)函數(shù).很明顯可以看出，在非晶態(tài)石英中僅存在短程有序性，而長程有序性消失.上述對關(guān)聯(lián)函數(shù)中第一峰的位置則分別對應(yīng)于非晶態(tài)石英中 Si—O，O—O，Si—Si的平衡距離，詳見表1.表中同時列出了Soules［24］，F(xiàn)euston和Garofalini［8］利用固定電荷方法得到的模擬結(jié)果，以及實驗結(jié)果［25，26］.不難看出，本文的模擬結(jié)果與實驗結(jié)果符合最好.

圖1 在300 K下表征非晶態(tài)石英結(jié)構(gòu)的 Si—O，O—O，Si—Si對關(guān)聯(lián)函數(shù)

表1 各種原子間距的比較

除了第一峰位置以外，Si—O—Si鍵和O—Si—O鍵的鍵角分布同樣是表征非晶態(tài)石英結(jié)構(gòu)的重要參數(shù).圖2給出了模擬得到的Si—O—Si鍵和O—Si—O鍵鍵角分布函數(shù).由圖可知，O—Si—O鍵集中分布于100°—120°的區(qū)域，峰值位于109°.這一角度與晶態(tài)石英中的 O—Si—O鍵角基本符合，表明在非晶石英中，Si仍然位于四面體中心，而與之相鄰的4個O則位于四面體頂點處.Si—O—Si鍵角的分布則寬泛得多，在110°—180°之間，且無明顯峰值.計算可得平均Si—O—Si鍵角為141°.這些特征與Feuston和Garofalini［8］的固定電荷分子動力學模擬結(jié)果基本一致，但平均 Si—O—Si鍵角比Feuston和 Garofalini的結(jié)果小.這是由于二者的電荷分布不同，從而導(dǎo)致了結(jié)構(gòu)上的差異.當然，兩種方法的結(jié)果均與實驗結(jié)果基本一致［8］.綜合對關(guān)聯(lián)函數(shù)與鍵角分布函數(shù)的結(jié)果可知，非晶態(tài)石英實際是由共角的四面體連接而成.在晶態(tài)石英中，這些四面體的排布規(guī)則，形成了有序的周期性晶體結(jié)構(gòu).而在非晶石英中，這種有序的周期性消失，但仍然保留了基本的四面體結(jié)構(gòu)單元.

圖2 模擬得到的非晶石英中Si—O—Si鍵和O—Si—O鍵的鍵角分布

我們進一步計算了非晶石英的振動頻譜.原子的振動頻譜直接依賴于原子間的相互作用力.因此，振動頻譜的計算可以直接反映分子動力學模擬結(jié)果.實驗上，非晶態(tài)石英的振動譜有三個主要的特征峰，其頻率分別為400，800和1100 cm-1，對應(yīng)于O在垂直于 Si—O—Si平面、在 Si—O—Si平面內(nèi)沿 Si—O—Si鍵角角平分線方向和在 Si—O—Si平面內(nèi)垂直于 Si—O—Si鍵角角平分線方向的振動.此外，通常還在200 cm-1左右出現(xiàn)一個峰值.在分子動力學中，振動頻譜可以通過對速度自關(guān)聯(lián)函數(shù)的Fourier變換得到［10］.圖3給出了非晶石英的振動頻譜，其峰值分別位于 187，354，770和 1104 cm-1，均與實驗數(shù)據(jù)符合較好.這無疑進一步證明了本文所采用的勢函數(shù)和變電荷分子動力學方法的正確性.

圖3 模擬得到的非晶石英的振動頻譜

此外，本文采用的變電荷分子動力學方法的計算結(jié)果表明，在非晶石英中，Si原子和O原子的電荷均不再是常數(shù)，而是出現(xiàn)了較大的漲落.圖4和5分別給出了體系內(nèi)Si原子平均電荷和O原子平均電荷隨模擬時間的演化.為便于比較，在圖4和5中，我們還給出了用變電荷分子動力學方法得到的相同溫度下 α-石英體系內(nèi)原子平均電荷隨模擬時間的演化.由圖可知，第一，非晶態(tài)與晶態(tài)石英內(nèi)的原子電荷顯著不同，這顯然是因為二者具有不同的結(jié)構(gòu).在非晶態(tài)石英中，部分Si—O鍵長僅有1.4左右，對應(yīng)的 Si原子電荷為1.78，同時部分 Si—O鍵長增至1.8左右，對應(yīng)的 Si原子電荷為1.46.而在晶態(tài)石英中Si—O鍵長為1.61，Si原子電荷在1.59左右.這表明隨著Si—O鍵長的減小，Si原子電荷增加較快，而隨著Si—O鍵長的增加，Si原子電荷減小較慢.因此，雖然在非晶石英中的平均Si—O鍵長仍為1.61，但是Si原子的平均電荷有所增加.類似地，O原子平均電荷則有所降低.原子電荷與鍵長間的這種非線性關(guān)系同樣存在于H2O等體系中［27］.第二，在同一種體系內(nèi)，粒子所帶電荷隨時間而變.這是由于粒子在平衡位置附近振動，對應(yīng)的其所帶電荷也在平衡電荷附近上下漲落.第三，非晶態(tài)石英中的粒子電荷漲落較大.這是由于非晶態(tài)石英結(jié)構(gòu)較不規(guī)則，粒子所處的環(huán)境變化較大所致.上述這些特征，顯然是任何固定電荷方法所無法描 述 的.雖 然 Feuston 和 Garofalini［8］以 及Garofalini［10］采用的固定電荷分子動力學模擬也能夠給出較為精確的非晶石英的結(jié)構(gòu)與振動特征，但是我們的模擬結(jié)果表明，變電荷方法還能夠正確地描述體系內(nèi)粒子電荷的變化.

圖4 非晶石英和α-石英體系內(nèi)Si原子平均電荷隨模擬時間的演化

圖5 非晶石英和α-石英體系內(nèi)O原子平均電荷隨模擬時間的演化

4.結(jié) 論

本文運用變電荷分子動力學方法模擬了非晶態(tài)石英的對關(guān)聯(lián)函數(shù)、鍵角分布函數(shù)和振動頻譜等.此外，還模擬了體系內(nèi)粒子電荷的演化，這是變電荷方法最顯著的特征.在非周期性體系中，這種演化往往顯得較為重要.例如，有文獻表明SiC中晶界薄膜的電荷轉(zhuǎn)移與薄膜平衡厚度相關(guān)［18］.本文對非晶態(tài)石英的模擬結(jié)果表明，變電荷方法不但能夠較為精確地描述該體系的結(jié)構(gòu)與振動特性，而且還揭示了非晶石英與晶態(tài)石英內(nèi)原子電荷的不同.本文的成功模擬也表明，本文采用的勢函數(shù)與變電荷方法還能進一步推廣，如用于模擬非晶或晶態(tài)石英體系內(nèi)的缺陷、表面、界面、熔融態(tài)石英等非周期性體系.

［1］Zhou Z X，Wang H L，Shen Y Q，Liu D J，Liu H，He S Y，Yang D Z 2008 Acta Phys.Sin.57 592(in Chinese)［周忠祥、王宏利、申艷青、劉大軍、劉 海、何世禹、楊德莊2008物理學報57 592］

［2］Yan F P，Wang L，Wei H，F(xiàn)u Y J，Jian W，Zheng K，Mao X Q，Li J，Liu L S，Peng J，Jian S S 2009 Acta Phys.Sin.58 1793(in Chinese)［延鳳平、王 琳、魏 淮、傅永軍、簡偉、鄭 凱、毛向橋、李 堅、劉利松、彭 健、簡水生2009物理學報58 1793］

［3］Jia W Y，Pei L W 1984 Acta Phys.Sin.33 1092(in Chinese)［賈惟義、裴力偉1984物理學報33 1092］

［4］Gao S J，Ouyang S X 2003 Acta Phys.Sin.52 1292(in Chinese)［高祀建、歐陽世翕2003物理學報52 1292］

［5］Wang J F，Li H L，Tang R M，Cha J X，He S A 1982 Acta Phys.Sin.31 1423(in Chinese)［王積方、李華麗、唐汝明、查濟璇、何壽安1982物理學報31 1423］

［6］Liu F，Garofalini S H，King-Smith D，Vanderbilt D 1993 Phys. Rev.B 49 12528

［7］Ma Y，Garofalini S H 2006 Phys.Rev.B 73 174109

［8］Feuston B P，Garofalini S H 1988 J.Chem.Phys.89 5818

［9］Feuston B P，Garofalini S H 1989 J.Chem.Phys.91 564

［10］Garofalini S H 1982 J.Chem.Phys.76 3189

［11］Mortier W J，Ghosh S K，Shankar S 1986 J.Am.Chem.Soc. 108 4315

［12］Ma Y，Garofalini S H 2008 J.Chem.Phys.128 084505

［13］Nakano A 1997 Comp.Phys.Commun.104 59

［14］Rappe A K，Goddard W A 1991 J.Phys.Chem.95 3358

［15］Rick S W，Stuart S J，Berne B J 1994 J.Chem.Phys.101 6141

［16］Ma Y，Garofalini S H 2006 J.Chem.Phys.124 234102

［17］Louwen J N，Vogt E T C 1998 J.Mol.Catal.A 134 63

［18］Ma Y，Chen S D，Xie G F 2009 Acta Phys.Sin.58 7792(in Chinese)［馬 穎、陳尚達、謝國鋒2009物理學報58 7792］

［19］Tsuneyuki S，Tsukada M，Aoki H，Matsui Y 1988 Phys.Rev. Lett.61 869

［20］van Beest W H，Kramer G，Santen R A 1990 Phys.Rev.Lett. 64 1955

［21］Nose S 1984 Mol.Phys.52 255

［22］Hoover W G 1985 Phys.Rev.A 31 1695

［23］Parrinello M，Rahman A 1980 Phys.Rev.Lett.45 1196

［24］Soules T F 1979 J.Chem.Phys.71 4570

［25］Mozzi R L，Warren B E 1969 J.Appl.Crystallogr.2 164

［26］Narten A H 1972 J.Chem.Phys.56 1905

［27］Guillot B，Guissani Y 2001 J.Chem.Phys.114 6720

PACS:61.43.Fs，61.43.Bn，63.50.Lm

Variable charge molecular dynamics simulation of vitreous silica*

Ma Ying
(Key Laboratory of Low Dimensional Materials and Application Technology of Ministry of Education，F(xiàn)aculty of Materials，Optoelectronics and Physics，Xiangtan University，Xiangtan 411105，China)

13 January 2010;revised manuscript

15 April 2010)

Variable charge molecular dynamics simulations of vitreous silica have been performed based on the iterative fluctuation charge model.The vitreous silica was formed using the standard melt-quench process.The pair distribution function，angle distribution and frequency spectrum were then obtained.The results are in good agreement with experimental data.More importantly，our results show that the atomic charge in vitreous silica is significantly different from those in the crystalline silica，and larger fluctuations are observed.

molecular dynamics，variable charge，vitreous silica

*國家自然科學基金(批準號:10702059)、高等學校博士學科點專項科研基金 (新教師基金)(批準號:20070530009)、教育部留學回國人員科研啟動基金(批準號:2008890)、國家博士后科學基金(批準號:20090451102，201003515)和湘潭大學材料科學與工程博士后流動站資助的課題.

*Project supported by the National Natural Science Foundation of China(Grant No.10702059)，the Research Fund for the Youth Scholars of the Specialized Research Fund for the Doctoral Program of Higher Education of China(Grant No.20070530009)，the Scientific Research Starting Foundation for the Returned Overseas Chinese Scholars of Ministry of Education，China(Grant No.2008890)，the China Postdoctoral Science Foundation(Grant Nos.20090451102，201003515)，and the Post-Doctoral Station of Materials Science and Engineering at Xiangtan University，China.