王世偉朱朋哲李瑞

1)(北京科技大學機械工程學院,北京 100083)

2)(北京交通大學機械與電子控制工程學院,北京 100044)

1 引 言

碳納米管[1]由于具有優(yōu)異的電學、力學和熱學性能,被期望應用于納米器件和納機電系統(tǒng)中[2,3].由于表面界面效應,納米尺度下摩擦和能量耗散更為顯著,這引起了研究者的重視.碳納米管在制備和化學改性過程中,表面通常帶有多種基團,這些表面基團對其納米尺度摩擦特性均存在重要影響.由于碳納米管和石墨烯在應用中具有很多相似之處,表面基團對它們摩擦的影響均受到了研究者的關注.Wal等[4]發(fā)現(xiàn)表面氟化納米管摩擦系數(shù)低至0.002—0.07,指出表面改性的碳納米管作為固體潤滑劑具有廣闊的應用前景.石墨烯氟化后情況則有所不同,Kwon等[5]的實驗和理論結果表明石墨烯在表面氟化后摩擦力增大.Ko等[6]指出石墨烯表面氫化、氟化和氧化后均具有較大的摩擦力.Dong等[7]認為氫化石墨烯摩擦力升高的原因是氫化導致表面的原子粗糙度增大.Li等[8]研究了碳納米管端帽與基底垂直接觸的情況,結果表明碳納米管所受的摩擦力隨界面間羥基含量的增加而增大.Wang等[9]采用密度泛函理論(density functional theory,DFT)研究了氧化石墨烯的層間滑移行為,結果表明界面間氫鍵的形成大大增加了界面間的靜摩擦力.Chen等[10]指出摩擦力隨滑動速度的變化取決于界面的化學性質(zhì),表面接枝諸如—OH,—COOH和—NH2等基團后能夠形成氫鍵網(wǎng)絡,摩擦力隨著滑動速度的增加而減小.Zheng等[11]的結果表明若石墨烯在Ge(111)基體上排列適當,即使在氟化或氧化之后,石墨烯的摩擦力也能夠保持在極小的水平.

摩擦過程總是伴隨著能量的轉(zhuǎn)移和耗散,一些學者對此做了深入研究.對于納米碳材料,如Eckstei等[12]研究了摻雜碳納米管的能量耗散;Kim等[13]指出隨著溫度升高,納米石墨諧振器的能量耗散增加;Hu等[14]綜述了原子尺度摩擦的模型,討論了能量耗散的機理,并分析了多層石墨烯能量耗散的具體過程.一些研究者也探討了能量耗散的影響因素.Wang等[15]采用一維PrandtlTomlinson(P-T)模型,研究了原子尺度摩擦的能量轉(zhuǎn)移和耗散,指出在黏滑運動過程中積累的能量并不總是完全耗散.考慮到聲子能量耗散的事實,動摩擦不僅取決于界面性質(zhì),也取決于固體中體積原子的性質(zhì).Kajita等[16]通過實驗和模擬比較了不同同位素單晶金剛石的動力學摩擦,觀察到13C金剛石的摩擦力比12C金剛石的約小3%,這表明13C金剛石摩擦力更小的原因在于原子質(zhì)量對聲子能量耗散的抑制.Cannara等[17]測量了表面接枝氫、氘的單晶金剛石和硅的摩擦力,結果表明接枝氫的表面具有更高的摩擦力,并指出表面振動頻率的改變將影響摩擦力.

上述研究工作解釋了界面基團影響納米碳材料摩擦特性的原因,分析了影響能量耗散的因素,但目前的工作中對碳納米管接枝羥基后摩擦特性變化的研究較少,并且未從聲子耗散的角度深入分析羥基影響界面能量耗散的過程.本文采用分子動力學模擬研究了界面接枝羥基后碳納米管在硅基底上摩擦力、能量和聲子態(tài)密度的變化,從聲子耗散的角度解釋了界面羥基影響碳納米管摩擦的機理.

2 模型和方法

本文的模型由單壁碳納米管SWCNT(10,10)和硅基底組成,如圖1所示.模型(a)為理想模型,包括水平取向的碳納米管和硅基底;在模型(b)中,僅硅基底接枝了羥基;在模型(c)中,碳納米管和硅基底均接枝了羥基.羥基和碳原子、硅原子之間以化學鍵的形式結合,其結合力比物理吸附力強,可以防止羥基脫附[18].碳納米管的直徑為1.356 nm,硅基底在x,y方向上的尺寸分別為8.01 nm和7.98 nm,硅基底底部原子固定以近似模擬大塊硅片.碳納米管羥基含量指羥基與碳納米管中碳原子數(shù)目之比,硅基底羥基含量指羥基與硅基底中最上層硅原子數(shù)目之比.模型在x[100],y[010]方向均采用周期性邊界條件.

圖1 碳納米管和硅基底組成的模型 (a)理想模型;(b)硅基底接枝羥基;(c)碳納米管和硅基底均接枝羥基Fig.1.Illustration of simulation models:(a)Ideal model:carbon nanotube on Si substrate;(b)Si substrate grafted with hydroxyl groups;(c)Carbon nanotube and Si substrate both grafted with hydroxyl groups.

碳納米管原子間的C—C成鍵作用采用AIREBO勢[19]描述,硅基底的原子間相互作用采用TERSOFF勢函數(shù)[20]描述,由于AIREBO勢函數(shù)中未考慮O原子,Argyris等[21,22]指出可采用OPLS_AA力場[23]描述羥基改性后的體系中界面上Si—O—H,C—O—H之間的相互作用.碳納米管與硅基底間僅存在范德瓦耳斯力作用,采用12-6 Lennard-Jones勢[24]描述,當碳納米管和硅基底均接枝羥基時會形成氫鍵,氫鍵作用由DREILDING力場[25]描述.

模擬系綜為正則系綜(NVT),體系溫度控制在300 K,采用廣義Langevin方程控制溫度,動力學方程的數(shù)值積分采用Velocity-Verlet算法,計算步長為1 fs.采用LAMMPS[26]實現(xiàn)分子動力學模擬過程.在模擬中首先進行能量最小化,隨后使體系弛豫100 ps至穩(wěn)定,之后給碳納米管一個恒定的速度使其在硅基底上沿運動方向勻速運動.

3 結果與討論

3.1 界面接枝羥基對碳納米管摩擦行為的影響

界面接枝羥基對碳納米管所受的摩擦力產(chǎn)生了重要影響.圖2顯示了接枝不同比例的羥基時碳納米管所受的側(cè)向力.為了更清晰地對比顯示各條曲線,圖2中對側(cè)向力數(shù)據(jù)進行了曲線擬合.可以看出,與模型(a)(即圖中Ideal對應的曲線)和模型(b)(即圖中0%/10%對應的曲線)相比,模型(c)(即圖中5%/5%,5%/10%,10%/10%,10%/20%對應的曲線)中界面均接枝羥基后,碳納米管所受的側(cè)向力顯著增大,并且側(cè)向力的起伏隨著接枝羥基數(shù)目的增加而增大.接枝不同比例的羥基后界面的平均摩擦力及平均氫鍵數(shù)目如圖3所示,可以看出平均摩擦力隨羥基數(shù)目的增加而逐步變大;界面間的平均氫鍵數(shù)量與摩擦力呈現(xiàn)相同趨勢,表明界面氫鍵數(shù)量與摩擦力緊密相關.Chen等[10]指出氫鍵作用對界面摩擦力的影響明顯,本文的結果與其結論一致.

圖2 碳納米管所受的側(cè)向力Fig.2.Lateral force on carbon nanotube.

圖3 接枝不同比例羥基時,界面的平均摩擦力和平均氫鍵數(shù)目Fig.3.Average friction force and average hydrogen bond numbers between interfaces when grafted with different hydroxyl group ratio.

本文也研究了界面接枝羥基后,碳納米管的手性角和直徑對摩擦的影響.圖4(a)顯示了直徑近似、不同手性角時單位長度碳納米管所受的平均摩擦力,其中SWCNT(11,9),SWCNT(12,8),SWCNT(13,7),SWCNT(14,6),SWCNT(15,0)的摩擦力類似,但SWCNT(10,10)的摩擦力比其他的摩擦力大,其原因在于SWCNT(10,10)為扶手椅型,碳原子在軸向規(guī)則排列,與理想的硅基底表面易形成摩爾紋結構,界面間相互作用更強.不同直徑時單位長度碳納米管所受的平均摩擦力如圖4(b)所示.可以看出,SWCNT(7,7)所受的平均摩擦力最小,SWCNT(10,10),SWCNT(15,15)和SWCNT(20,20)所受的摩擦力類似,SWCNT(25,25)所受的摩擦力最大.即直徑對摩擦力有一定的影響,摩擦力會隨直徑的增加而增大.其原因在于隨著碳納米管直徑的增大,碳納米管底部受基底的吸引而形成平臺,導致其與硅基底間的接觸面積增大,進而引起摩擦力的增加.

圖4 單位長度碳納米管所受的平均摩擦力 (a)不同手性角條件下;(b)不同直徑條件下Fig.4.Average friction force on carbon nanotube per unit length:(a)Different chiral angles;(b)different diameters.

3.2 界面接枝羥基對體系能量耗散的影響

為了觀察體系的能量耗散情況,我們統(tǒng)計了接枝不同比例的羥基時體系總能量的變化,如圖5所示.可以看出,當碳納米管接枝羥基的比例相同時,體系總能量基本相同,硅基底接枝羥基的比例對體系總能量影響較小.當碳納米管接枝羥基的比例不同時,總能量則存在較大差異,如碳納米管接枝羥基比例為10%與5%時總能量差值達470 eV.總體來看,體系總能量的變化趨勢與平均摩擦力具有相似性,即隨著接枝羥基比例的增加而增大.對于某一具體的體系,其總能量在模擬過程中基本保持穩(wěn)定.

圖5 接枝不同比例的羥基時體系的總能量Fig.5.Total energy of the simulation system when grafted with different hydroxyl group ratio.

圖6 接枝不同比例的羥基時碳納米管的勢能Fig.6.Potential energy of carbon nanotube when grafted with different hydroxyl group ratio.

接枝不同比例的羥基時,碳納米管的勢能變化如圖6所示,與理想模型(即圖中的Ideal模型)以及僅硅基底接枝羥基的模型(即羥基比例為0%/10%的模型)相比,界面均接枝羥基后碳納米管的勢能明顯增大.碳納米管的勢能隨其接枝羥基比例的增加而增大,但當碳納米管接枝羥基比例相同、硅基底接枝羥基比例不同時,碳納米管的勢能隨基底接枝羥基比例的增大而減小;當界面接枝羥基比例為10%/10%和10%/20%時,這一趨勢表現(xiàn)的更為明顯.其原因在于界面間的相互作用隨基底接枝羥基數(shù)目的增加而增強,如界面氫鍵數(shù)量增加.本文對不同手性角、不同直徑碳納米管的勢能和總能量也進行了統(tǒng)計分析,其規(guī)律與SWCNT(10,10)相同.

為了探討界面接枝羥基影響摩擦和能量耗散的深層原因,我們對體系聲子態(tài)密度的變化進行了分析.體系的振動信息是聲子態(tài)密度研究的基礎,計算時使用速度自相關函數(shù)法提取體系振動信息并進一步做頻譜分析.這種方法由Dickey等[27]在1969年通過理論推導提出,其基本原理為:體系在穩(wěn)態(tài)振動時,其速度自相關函數(shù)和振動譜互為Fourier變換/反變換.本文選取相應特征時間點并輸出當時1 ps范圍內(nèi)的原子坐標與速度信息,利用速度自相關函數(shù)法研究體系的振動并提取振動信息做頻譜分析,比較不同模型及不同狀態(tài)的振動態(tài)密度分布.

接枝不同比例的羥基時,對碳納米管的所有原子進行統(tǒng)計得到其聲子態(tài)密度如圖7所示.可以看出,接枝不同比例的羥基,其聲子態(tài)密度存在較大的差異,理想模型和羥基比例為0%/10%時聲子態(tài)密度基本相同(如圖7(a)、圖7(b)所示),羥基比例為5%/5%和5%/10%(如圖7(c),(d)所示),10%/10%和10%/20%時的聲子態(tài)密度比較相似(如圖7(e),(f)所示),其原因是碳納米管接枝羥基后羥基的振動在碳納米管的聲子態(tài)密度中開始體現(xiàn),碳納米管接枝羥基比例相同則呈現(xiàn)出類似的結果.

當碳納米管不接枝羥基和接枝羥基比例為0%/10%時,如圖7(a),(b)所示,聲子態(tài)密度圖中存在9,18和51 THz三個振動峰.因為此時碳納米管均未接枝羥基,所以可以確定這三個振動峰為碳納米管的固有振動峰.Dresselhaus等[28]指出SWCNT的拉曼光譜中環(huán)呼吸振動峰(radial breathing mode,RBM)、無規(guī)振動峰(D峰)、剪切振動峰(G峰)的波數(shù)分別為180,1350和1582 cm?1,根據(jù)k=f/c(式中k為波數(shù),f為頻率,c為光速),計算得到對應的頻率分別為5.40,40.47和47.43 THz,其他譜線在SWCNT的拉曼光譜中也會出現(xiàn),如環(huán)呼吸振動的倍頻峰(2RBM)[29].拉曼光譜中D峰是缺陷和無序度的反映,G峰是有序度的反映,由于本文計算使用的為理想的SWCNT,應該不存在D峰,或者D峰非常微弱.對比模擬的振動峰,可以看出RBM,G峰出現(xiàn)了大約4 THz的頻移;D峰峰值不明顯,18 THz時為2RBM.這一結果表明計算得到的聲子態(tài)密度與拉曼光譜存在良好的對應性.

圖7 接枝不同比例的羥基時,碳納米管的聲子態(tài)密度圖譜 (a)Ideal;(b)0%/10%;(c)5%/5%;(d)5%/10%;(e)10%/10%;(f)10%/20%Fig.7.Phonon density of states of carbon nanotube when grafted with different hydroxyl groups ratio:(a)Ideal;(b)0%/10%;(c)5%/5%;(d)5%/10%;(e)10%/10%;(f)10%/20%.

羥基的引入對碳納米管晶格振動產(chǎn)生了調(diào)制,對碳納米管振動模態(tài)存在較大影響.如圖7(c)—(f)所示,引入羥基以后碳納米管聲子態(tài)密度出現(xiàn)了37,43,67,86 THz四個新的振動峰,這四個振動峰即為羥基的振動峰,并且羥基振動峰成為主要振動峰.隨著羥基比例的增加,碳納米管固有頻率處峰值降低.其原因在于碳納米管接枝羥基后,碳原子和氧原子形成共價鍵,并且氧原子質(zhì)量較大,其振動對碳原子的振動擾動較大,使碳納米管固有的晶格振動受到抑制,因此碳納米管固有頻率處峰值降低.由于氫原子質(zhì)量比較小,且氫原子不和碳原子直接接觸,因此對碳納米管的振動影響較小.隨著羥基比例的增加,37和43 THz頻率處峰值降低,因此這兩個峰值是由氧原子振動引起的,原因是羥基增加以后,碳原子和氧原子形成的共價鍵增多,它們之間相互抑制,導致氧原子振動峰值降低.當碳納米管接枝羥基比例由5%變?yōu)?0%時,氫原子數(shù)目變化較小,67和86 THz處峰值隨羥基比例的變化不明顯,可以推斷此處的振動峰由氫原子引起.

聲子態(tài)密度曲線所覆蓋區(qū)域代表了所激發(fā)聲子數(shù)目的多少,聲子態(tài)密度大,體系耗散的內(nèi)能多.曲線峰值及平穩(wěn)過渡段的升高和降低代表了模型激發(fā)出該頻率聲子數(shù)目的增加與減少.碳納米管接枝羥基后,從聲子態(tài)密度的變化可以觀察能量耗散方式的變化.如圖7(c)—(f)所示,隨界面接枝羥基的增加,碳納米管原有振動峰的峰值降低,羥基振動峰逐漸成為主要振動峰,即碳納米管中的能量耗散方式由碳納米管的振動逐漸向羥基的振動轉(zhuǎn)變.當碳納米管接枝羥基的比例達到10%時,通過羥基耗散能量成為了碳納米管中能量耗散的主要方式.

圖8 接枝不同比例羥基時,整個體系的聲子態(tài)密度圖譜 (a)Ideal;(b)0%/10%;(c)5%/5%;(d)5%/10%;(e)10%/10%;(f)10%/20%Fig.8.Phonon density of states of the simulation system when grafted with different hydroxyl group ratio:(a)Ideal;(b)0%/10%;(c)5%/5%;(d)5%/10%;(e)10%/10%;(f)10%/20%.

圖9 體系接枝羥基比例為5%/10%時,不同階段的聲子態(tài)密度圖譜,其中(a)—(c)為接枝羥基碳納米管的聲子態(tài)密度,(d)—(f)為整個體系的聲子態(tài)密度 (a)碳納米管未接觸硅基底前的自由狀態(tài);(b)弛豫階段;(c)碳納米管在硅基底上運動;(d)碳納米管未接觸硅基底前的自由狀態(tài);(e)弛豫階段;(f)碳納米管在硅基底上運動Fig.9.Phonon density of states at different stages when the hydroxyl group ratio is 5%/10%(a)—(c)are the phonon density of states of carbon nanotube grafted with hydroxy groups,(d)—(f)are the phonon state density of the whole simulation system:(a)The free state before carbon nanotube contacts with Si substrate;(b)after relaxation;(c)carbon nanotube moves on the Si substrate;(d)the free state before carbon nanotube contacts with Si substrate;(e)after relaxation;(f)carbon nanotube moves on the Si substrate.

體系的聲子態(tài)密度等于體系各部分聲子態(tài)密度之和,硅基底及其接枝羥基的聲子振動模態(tài)在整個體系的聲子態(tài)密度中可以體現(xiàn),因此此處給出了整個體系聲子態(tài)密度的變化,如圖8所示.接枝不同比例的羥基時整個體系圖譜均存在碳納米管的三個固有振動峰,如圖8(a)所示.Hart等[30]指出硅的拉曼光譜單峰在520 cm?1附近,轉(zhuǎn)換成頻率為14.90 THz附近,因此頻率約18 THz的峰值為碳納米管的2RBM與硅振動峰的疊加.因為體系中硅的數(shù)量很多,此頻率處的振動峰最為顯著.圖8(c)—(f)中其他峰為羥基的振動峰,頻率與圖7中類似.從圖8(c)到圖8(f),體系接枝羥基的數(shù)目逐漸增加,羥基引起振動峰的峰值逐漸增大.圖8(c)和圖8(d)相比,硅基底接枝羥基數(shù)目增加,37和86 THz的振動峰有明顯的增大,其他峰值變化不大;而圖8(e)和圖8(f)相比,隨著羥基數(shù)目增加,后者除37和86 THz的振動峰有增加外其他頻率振動峰均明顯降低,峰值的降低是由于接枝羥基數(shù)目較多而導致碳納米管和硅基底的晶格振動受到了抑制.從圖8(a)—(f)可明顯看出,隨著接枝羥基比例的增加,羥基振動占整體能量耗散的比例越來越大,當界面接枝羥基比例為10%/20%時,能量耗散的主要方式由碳納米管和硅基底的振動轉(zhuǎn)變?yōu)榻缑媪u基的振動.

此外,圖8(c)—(f)也顯示,隨著體系接枝羥基數(shù)目的增加,羥基振動對應的聲子態(tài)密度峰逐漸增大,聲子態(tài)密度曲線所覆蓋區(qū)域變大,能量耗散增加.由圖5、圖6可以看出體系的總能量隨著體系接枝羥基比例的增加而增大,碳納米管的勢能也隨碳納米管上接枝羥基比例的增大而增大,但隨著硅基底上接枝羥基比例的增大而減小.綜合比較圖5、圖6、圖8可以看出,總體來看,體系接枝羥基的比例和總能量直接相關,總能量大時體系的能量耗散也大.但接枝羥基比例為10%/20%和10%/10%時情況則有所不同,10%/20%體系接枝的羥基比例比10%/10%體系的更高,但圖8(f)中碳納米管對應的振動峰降低,原因是由于界面作用增強,碳納米管的勢能降低,這說明碳納米管的能量耗散情況與體系整體的能量耗散存在差異.

我們以羥基比例為5%/10%的體系為例,探討接觸和剪切運動對能量耗散的影響,如圖9所示.與碳納米管未接觸硅基底前(圖9(a))相比,碳納米管在硅基底上運動時,碳納米管聲子態(tài)密度(圖9(c))的峰值均有明顯的降低,且低頻振動峰消失;碳納米管在硅基底上運動時,與碳納米管未接觸硅基底前(圖9(d))相比,整個體系的聲子態(tài)密度圖譜(圖9(f))中氫原子引起的86 THz振動峰增大,同樣低頻振動峰也消失.其原因可歸結為當碳納米管和硅基底接觸后,由于界面羥基間的相互作用,原子的振動狀態(tài)發(fā)生了變化,整體表現(xiàn)為碳納米管上羥基的振動減弱,硅基底上氫原子的振動增加.弛豫階段的聲子態(tài)密度(圖9(b),(e))與運動階段的差別不大,這表明碳納米管和硅基底從未接觸到接觸對能量耗散的影響比剪切運動更大.

4 結 論

本文采用分子動力學模擬研究了界面接枝羥基對碳納米管摩擦和能量耗散方式的影響,得到了以下結論:

1)界面間引入羥基后,由于界面間氫鍵的形成與斷裂,碳納米管所受的摩擦力明顯增加,界面接枝羥基比例越高,平均摩擦力越大.手性角對摩擦力有一定影響,扶手椅型碳納米管的摩擦力比手性型和鋸齒型大.直徑的增加會導致摩擦力增大,原因是碳納米管底部與硅基底間的接觸面積增大.

2)體系的總能量及碳納米管的勢能隨羥基數(shù)目增加而逐步增加.碳納米管接枝羥基比例相同的情況下,系統(tǒng)的總能量基本相同.碳納米管接枝羥基比例保持不變時,碳納米管的勢能隨硅基底接枝羥基比例升高而降低,接枝羥基比例較高時這一現(xiàn)象更為明顯,原因是硅基底和碳納米管之間的相互作用變大.總體來看,體系總能量大時體系的能量耗散也大,但通過碳納米管的振動耗散的能量也可能降低.

3)界面接枝羥基后,接枝羥基的碳納米管和整個體系的聲子態(tài)密度均出現(xiàn)羥基對應的振動峰.隨接枝羥基數(shù)目的增加,羥基對應的振動峰峰值增大,碳納米管和硅基底的固有振動峰降低,界面間能量耗散的主要途徑由碳納米管和硅基底的振動轉(zhuǎn)變?yōu)榻缑骈g羥基的振動.

感謝密歇根大學彭慶博士和北京科技大學謝璐博士的討論.

[1]Iijima S 1991Nature354 56

[2]Scarselli M,Castrucci P,de Crescenzi M 2012J.Phys.:Condens.Matter24 313202

[3]Liew K M,Wong C H,He X Q,Tan M J,Meguid S A 2004Phys.Rev.B69 1738

[4]van der Wal R L,Miyoshi K,Street K,Tomasek A,Peng H,Liu Y,Margrave J,Khabashesku V 2005Wear259 738

[5]Kwon S,Ko J H,Jeon K J,Kim Y H,Park J Y 2012Nano Lett.12 6043

[6]Ko J H,Kwon S,Byun I S,Jin S C,Park B H,Kim Y H,Park J Y 2013Tribol.Lett.50 137

[7]Dong Y L,Wu X W,Martini A 2013Nanotechnology24 375701

[8]Li R,Mi J X 2017Acta Phys.Sin.66 046101(in Chinese)[李瑞,密俊霞 2017物理學報 66 046101]

[9]Wang L F,Ma T B,Hu Y Z,Wang H 2012Phys.Rev.B86 125436

[10]Chen J,Ratera I,Park J Y,Salmeron M 2006Phys Rev.Lett.96 236102

[11]Zheng X,Lei G,Yao Q,Li Q,Miao Z,Xie X,Qiao S,Wang G,Ma T,Di Z,Luo J,Wang X 2016Nat.Commun.7 13204

[12]Eckstein K H,Hartleb H,Achsnich M M,Sch?ppler F,Hertel T 2017ACS Nano11 10401

[13]Kim S Y,Park H S 2009Appl.Phys.Lett.94 101918

[14]Hu Y Z,Ma T B,Wang H 2013Friction1 24

[15]Wang Z J,Ma T B,Hu Y Z,Xu L,Wang H 2015Friction3 170

[16]Kajita S,Tohyama M,Washizu H,Ohmori T,Watanabe H,Shikata S 2015Tribology Online10 156

[17]Cannara R J,Brukman M J,Cimatu K,Sumant A V,Baldelli S,Carpick R W 2007Sci.318 780

[18]Sun Y,Yang S,Chen Y,Ding C,Cheng W,Wang X 2015Environ.Sci.Technol.49 4255

[19]Brenner D W,Shenderova O A,Harrison J A,Stuart S J,Ni B,Sinnott S B 2002J.Phys.:Condes.Matter14 783

[20]Tersof fJ 1988Phys.Rev.B37 6991

[21]Argyris D,Tummala N R,Striolo A,Cole D R 2008J.Phys.Chem.C112 13587

[22]Hughes Z E,Shearer C J,Shapter J,Gale J D 2012J.Phys.Chem.C116 24943

[23]Damm W,Frontera A,Tirado-Rives J,Jorgensen W L 1997J.Comput.Chem.18 1955

[24]Ruof fR S,Hickman A P 1993J.Phys.Chem.97 2494

[25]Mayo S L,Olafson B D,Goddard W A III 1990J.Phys.Chem.94 8897

[26]Plimpton S 1995J.Comput.Phys7 1

[27]Dickey J M,Paskin A 1969Phys.Rev.188 1407

[28]Dresselhaus M S,Dresselhaus G,Saito R,Jorio A 2005Phys.Rep.409 47

[29]Yin Y,Vamivakas A N,Walsh A G,Cronin S B,Unlü M S,Goldberg B B,Swan A K 2007Phys.Rev.Lett.98 037404

[30]Hart T R,Aggarwal R L,Lax B 1970Phys.Rev.B1 638