摘 要：運(yùn)用分子動力學(xué)方法，采用Hamid模型構(gòu)建初始結(jié)構(gòu)，模擬了水化硅酸鈣的結(jié)構(gòu)，并分析比較了不同狀態(tài)水化硅酸鈣結(jié)構(gòu)的異同，且通過與試驗(yàn)數(shù)據(jù)的比較驗(yàn)證了模擬結(jié)果的可靠性。模擬結(jié)果表明：晶態(tài)水化硅酸鈣的結(jié)構(gòu)長、短程都有序，而其無定形態(tài)只存在近程有序的結(jié)構(gòu)特點(diǎn)；無定形水化硅酸鈣的基本結(jié)構(gòu)單元同其晶態(tài)一樣仍為硅氧四面體，且硅氧四面體鏈以Q2形式連接；模擬得到的原子間距離、配位數(shù)等結(jié)構(gòu)參數(shù)基本與實(shí)驗(yàn)值相符合。
關(guān)鍵詞：水化硅酸鈣；結(jié)構(gòu)；分子動力學(xué)
0前言
水化硅酸鈣是決定水泥性能的關(guān)鍵組分，因此其組成、結(jié)構(gòu)及性質(zhì)自上世紀(jì)50年代Grudemo[1]和Taylor[2]的開創(chuàng)性工作以來一直是水泥科學(xué)研究中的重要內(nèi)容，各國學(xué)者都進(jìn)行了深入研究，提出了一系列結(jié)構(gòu)模型，如類Tobermorite和Jennite模型[3]、固溶體模型[4]、固溶體模型[5]和納米結(jié)構(gòu)、中介結(jié)構(gòu)假說[6]。同時，水化硅酸鈣結(jié)構(gòu)的研究方法和手段也變得越來越豐富，如運(yùn)用X-射線衍射根據(jù)其彌散的衍射峰位置確定其結(jié)構(gòu)大致類似于托貝莫來石[3]，運(yùn)用核磁共振測定水化硅酸鈣結(jié)構(gòu)中硅氧四面體的聚合狀態(tài)[7]，運(yùn)用X-射線光電子能譜根據(jù)O1s結(jié)合能的變化來區(qū)別橋氧與非橋氧[8]，運(yùn)用場發(fā)射高分辨透射電鏡觀察無定形水化硅酸鈣中的納米結(jié)晶區(qū)[6]。尤其是近年來分子動力學(xué)（molecular dynamics, MD）的發(fā)展，為在原子水平上探討水化硅酸鈣的精細(xì)結(jié)構(gòu)提供了一種嶄新的方法。上世紀(jì)90年代中期，F(xiàn)aucon[9]就運(yùn)用分子動力學(xué)方法成功地模擬了水化硅酸鈣結(jié)構(gòu)不穩(wěn)定地原因、硅氧四面體鏈的斷裂機(jī)理和陽離子（Al）取代硅后對結(jié)構(gòu)的影響，但緊隨其后未見任何水化硅酸鈣的分子動力學(xué)模擬成果見諸于報道。值得注意的是，對用作模型參照的晶態(tài)托貝莫來石、羥基硅鈣石結(jié)構(gòu)也從晶體學(xué)角度上詳細(xì)地研究了它們的對稱性、晶胞參數(shù)、鍵長等
結(jié)構(gòu)性質(zhì)[10,11]，這為構(gòu)建C–S–H分子動力學(xué)模擬的初始結(jié)構(gòu)提供了精確的數(shù)據(jù)依據(jù)。模擬過程中，所使用的分子動力學(xué)軟件為DL_POLY(http://www.cse.clrc.ac.uk/msi/software/DL_POLY），本課題組對其定義勢函數(shù)的源代碼部分進(jìn)行了適當(dāng)?shù)男薷摹?
1 模擬實(shí)驗(yàn)
1.1 初始結(jié)構(gòu)創(chuàng)建
經(jīng)反復(fù)試驗(yàn)本文選擇了Hamid模型[11]來構(gòu)建初始結(jié)構(gòu)。其晶體學(xué)參數(shù)如下所示?；瘜W(xué)式：Ca4Si6O14(OH)4·2H2O；晶系：單斜；晶胞參數(shù)：am=0.669 nm，bm=0.739 nm，cm=2.277 nm，γ＝123.49°；空間點(diǎn)群：P21；模擬尺寸設(shè)為5am×4bm×cm，即am方向尺寸為3.695 nm，bm方向尺寸為2.956 nm，cm方向尺寸為6.831 nm；各原子對應(yīng)的離子電荷采用經(jīng)驗(yàn)電荷（formal charges），如Si的電荷?。?。為了考慮Si–H，O–H，H–H的兩體作用勢和Si–O–H，H–O–H的三體作用勢，必須將OH和H2O拆為單個原子。查得O–H鍵長0.96 nm，∠Si–O–H為109.5°，∠H–O–H為104.5°，于是本文在設(shè)定O的位置就是OH/H2O的位置的前提下并根據(jù)鍵長鍵角求得H的位置，如表1所示。這樣處理可能存在如下問題：一是不一定符合原始的對稱關(guān)系；二是這樣生硬地將H放在這一位置可能使之與另一原子相距很近，進(jìn)而導(dǎo)致排斥強(qiáng)烈，使結(jié)構(gòu)不穩(wěn)定。因此，為了消除這樣處理對結(jié)構(gòu)造成的潛在不利影響，在選定相互作用勢及參數(shù)后本文還對初始結(jié)構(gòu)進(jìn)行了優(yōu)化處理，處理方法如下：在正則系綜（NVT）中，300 K下平衡5000步（時間步長為0.1 fs）；然后以NVT系綜的輸出作為微正則系綜（NVE）的初始輸入，在相同溫度下以相同步長再平衡5000步。
1.2 相互作用勢設(shè)定
該過程的目的就是選擇適當(dāng)?shù)南嗷プ饔脛萁?jīng)驗(yàn)函數(shù)及其參數(shù)，并設(shè)定分子動力學(xué)模擬控制參數(shù)。其方法為以本文提出的初始結(jié)構(gòu)作為輸入，選定勢函數(shù)及設(shè)定參數(shù)后在NVE下運(yùn)行1000步，其中平衡900步。因?yàn)镹VE系綜在合適的初始條件下能夠很快達(dá)到平衡（通常為數(shù)百步），若系統(tǒng)能夠達(dá)到平衡且數(shù)據(jù)不溢出（沒有太高的能量作用），則勢函數(shù)和參數(shù)是合理的。
（1）兩體作用勢經(jīng)驗(yàn)函數(shù)選擇及參數(shù)設(shè)定
由于DL_POLY單獨(dú)定義庫侖勢，因此兩體作用勢就是van de Waals力短程作用勢。經(jīng)過反復(fù)試驗(yàn)，本文最終選擇了BHM（Born-Huggins-Meyer）勢，其函數(shù)表達(dá)式如式（1），各參數(shù)見表2。
為了體現(xiàn)水的作用（他們有可能離解為H+和OH-），Si–H，H–H，O–H的兩體作用勢除了BHM勢還需加上如式（2）所示的額外項(xiàng)[12,13]。Si–H只需加一項(xiàng)，而H–H則需加兩項(xiàng)，O–H更是需加三項(xiàng)，各參數(shù)見表2。兩體作用勢的截?cái)喟霃饺?.76 nm。
（2）三體作用勢經(jīng)驗(yàn)函數(shù)選擇及參數(shù)設(shè)定
三體作用勢是一種比van de Waals力作用勢更短程、更弱的作用勢，其對總能量的貢獻(xiàn)通常弱于1％。本文所選三體作用勢經(jīng)驗(yàn)函數(shù)的表達(dá)式如式（3），其參數(shù)如表3所示，截?cái)喟霃饺?.345 nm。
?。?）庫侖勢算法選取及參數(shù)設(shè)定
根據(jù)本文選擇的平行六面體周期性邊界條件，相應(yīng)地庫侖勢算法選用Ewald求和算法。DL_POLY在指定長程作用力（靜電作用力）截?cái)嗪颓蠛途龋ㄏ鄬φ`差）的前提下，可由其內(nèi)部函數(shù)估算出求和精確性的控制參數(shù)，包括收斂參數(shù)α，靜電作用力在正空間的截?cái)喟霃?I>rcut和最大倒易矢量。
1.3 結(jié)構(gòu)模擬
采用NVE系綜，速度和位置算法采用LF（Leap-Frog）算法，邊界條件選為平行六面體周期性邊界條件，步長設(shè)為1 fs，速度、位置、力等軌跡數(shù)據(jù)和結(jié)構(gòu)能量、總能量、溫度等統(tǒng)計(jì)數(shù)據(jù)每10步存取一次，總共平衡20 000步。為了節(jié)省模擬時間，還采用了復(fù)合時間步（multiple timestep）算法[14]。
2 模擬結(jié)果與討論
2.1 偏徑向分布函數(shù)
圖1為Ca-O，Si-O，O-O，Si-Si的偏徑向分布函數(shù)（partial radial distribution function, RDF），圖中Crystal對應(yīng)的結(jié)構(gòu)為優(yōu)化處理過的初始結(jié)構(gòu)，Liquid對應(yīng)的結(jié)構(gòu)為以優(yōu)化過的初始結(jié)構(gòu)為輸入制備得到的高溫結(jié)構(gòu)，Amorphism對應(yīng)的結(jié)構(gòu)為由高溫結(jié)構(gòu)冷卻得到的常溫結(jié)構(gòu)。
由圖1知，Crystal對應(yīng)結(jié)構(gòu)的各RDF的各層次峰布滿整個半徑r取值范圍，且在較大半徑r范圍內(nèi)其值并不趨于1，說明該結(jié)構(gòu)不但近程有序而且遠(yuǎn)程也是有序的；Amorphism和Liquid對應(yīng)結(jié)構(gòu)的各RDF只存在較明顯的第一近鄰峰和明顯寬化的第二近鄰峰，其他各峰均消失，且在較大半徑r范圍內(nèi)其值趨于1，表明Amorphism和Liquid對應(yīng)結(jié)構(gòu)處于近程有序、遠(yuǎn)程無序的結(jié)構(gòu)狀態(tài)。比較Amorphism和Liquid對應(yīng)結(jié)構(gòu)的各RDF還發(fā)現(xiàn)，前者的第一近鄰峰都要比后者的尖銳，即前者峰的半高寬要比后者的小，表明前者的近程有序程度要比后者的明顯。另外仔細(xì)研究Amorphism的Si-O的RDF發(fā)現(xiàn)，其第一峰不僅比Crystal的低而且還寬，另外在其第二近鄰峰中的4 Å附近出現(xiàn)4了峰的劈裂現(xiàn)象（如圖1(c)中的小圖所示），這是明顯的無定形態(tài)特征。綜上所述，Crystal對應(yīng)的結(jié)構(gòu)就是晶態(tài)的C-S-H，而Liquid對應(yīng)的結(jié)構(gòu)則是其高溫液態(tài)相，Amorphism對應(yīng)的結(jié)構(gòu)則是無定形的C-S-H。
2.2 配位數(shù)
配位數(shù)（coordination numbers, CN）在一定程度上也能反應(yīng)模擬對象的局域結(jié)構(gòu)。比較圖2可知，晶態(tài)C-S-H對應(yīng)的各CN曲線在半徑r取值范圍內(nèi)總存在著幾個明顯的臺階，說明它不但有確切的最近鄰配位，而且還有次近鄰甚至更遠(yuǎn)近鄰的配位，這與晶體的近、遠(yuǎn)程都有序的結(jié)構(gòu)特征是一致的。液態(tài)結(jié)構(gòu)和無定形結(jié)構(gòu)對應(yīng)的CN曲線只是在較小半徑r范圍內(nèi)才有臺階，而在較大范圍內(nèi)曲線幾乎呈直線上升，說明他們只存在相距很近的近鄰配位，而在相距較遠(yuǎn)的范圍內(nèi)原子個數(shù)與相距半徑r幾乎呈正比，即不存在特定的長程配位關(guān)系。這與液體和無定形體的近程有序而遠(yuǎn)程無序的結(jié)構(gòu)特點(diǎn)是一致的。
為了進(jìn)一步考察三者的近程結(jié)構(gòu)，由圖2（c）中的小圖知，無論是晶體還是無定形體、液體都在CN為4處出現(xiàn)了較寬范圍內(nèi)的長臺階，表明Si的O最近鄰配位在三者中都是4配位，說明三者的基本結(jié)構(gòu)單元均為硅氧四面體。另外由圖2（d）中的小圖知，無定形體在CN等于2處出現(xiàn)了較寬范圍內(nèi)的長臺階且在較大半徑r范圍內(nèi)并無其他臺階出現(xiàn)，表明無定形水化硅酸鈣中Si與Si之間存在2配位關(guān)系，即硅氧四面體鏈以Q2形式連接，這與Wieker的NMR測試結(jié)果非常一致。
2.3 結(jié)構(gòu)參數(shù)
表4為晶態(tài)C-S-H和無定形C-S-H的鍵長或原子間距離、配位數(shù)等結(jié)構(gòu)信息。鍵長或原子間距離取RDF第一近鄰峰的峰值對應(yīng)的半徑r。配位數(shù)取CN曲線的第一拐點(diǎn)對應(yīng)的值。
注：圓括號中的數(shù)字為計(jì)算配位數(shù)時的截?cái)喟霃?；a為從文獻(xiàn)[14]中得到的0.9 nm托貝莫來石數(shù)據(jù)；b為從文獻(xiàn)[14]中得到的1.1 nm托貝莫來石數(shù)據(jù)；c為鍵長平均值。
由表4知，除了Ca的O配位數(shù)外其他數(shù)據(jù)都與實(shí)驗(yàn)值基本符合，因此在考慮到模擬體系尺寸有限、邊界效應(yīng)和統(tǒng)計(jì)誤差的影響的前提下，可以認(rèn)為本文采用的相互作用勢函數(shù)和初始結(jié)構(gòu)是合理的，模擬結(jié)果基本能重現(xiàn)實(shí)際的微觀結(jié)構(gòu)。至于Ca的O配位數(shù)與實(shí)驗(yàn)值有較大差別，其原因在于MD和實(shí)驗(yàn)的統(tǒng)計(jì)方法不同。
3 結(jié) 論
以Hamid的托貝莫來石結(jié)構(gòu)模型為基礎(chǔ)，將H2O、OH分解為單個原子按照鍵長、鍵角要求置于某一位置，從而構(gòu)建得到初始結(jié)構(gòu)。由MD模擬得到的偏徑向分布函數(shù)知，晶態(tài)C–S–H不僅長程有序而且短程也有序，而無定形C–S–H確實(shí)只存在近程有序的結(jié)構(gòu)特點(diǎn)。由各原子對的配位情況知無定形C-S-H的基本結(jié)構(gòu)單元同晶態(tài)C-S-H一樣仍為硅氧四面體，且硅氧四面體鏈以Q2形式連接。在采用合理的相互作用勢和初始結(jié)構(gòu)的前提下，MD模擬得到了與實(shí)驗(yàn)值基本一致的原子間距離、配位數(shù)等結(jié)構(gòu)參數(shù)。
參考文獻(xiàn)
[1] A. Grudemo. Discussion of the structure of cement hydration compounds[A]. Proceeding of the third international symposium on the chemistry of cement [C]. London: 1952. 247–253.
[2] H. F. W. Taylor. Proposed structure for calcium silicate hydrates gel [J]. J Am Ceram Soc, 1986, 69(6): 464–467.
[3] H. F. W. Taylor. Cement Chemistry [M]. 2nd edition, London: Thomas telford, 1997. 128–134.
[4] M. W. Grutzeck, T. J. Larosa, S. Kwan. Characteristics of C-S-H gels [A]. Proceeding of the 10th ICCC [C], Gothenburg, Sweden: 1997. 2ii067.
[5] I. G. Richardson, G. W. Groves. The composition and structure of C-S-H in hardened slag cement pastes [A]. Proceedings of the 10th ICCC [C], Gothenburg, Sweden, 1997. 2ii068.
[6] D. Viehland, J. F. Li, L. J. Yuan, et al. Mesostructure of calcium silicate hydrate (C-S-H) gels in portland cement paste: Short-range ordering, nanocrystallinity, and local compositional order [J]. J Am Ceram Soc, 1996, 79(7): 1 731–1 744.
[7] X. Cong, R. J. Kirkpatrick. 29Si MAS NMR study of the structure of calcium silicate hydrate [J]. Adv Cem Based Mater, 1996, 3(3): 144–156
[8] L. Black, A. Stumm, K. Garbev, et al. X-ray photoelectron spectroscopy of aluminium-substituted tobermorite [J]. Cem Concr Res, 2005, 34(1): 51–55
[9] P. Faucon, J. M. Delaye, J Yirlet, et al. Study of the structural properties of the C-S-H (I) by molecular dynamics simulation [J]. Cement and Concrete Research, 1997, 27(10): 1 581–1 590.
[10] S. Merlino, E. Bonaccoris, T. Armbruster. Tobermorites: their real structure and order-disorder (OD) character [J]. American Mineralogist, 1999, 84: 1 613–1 621.
[11] S. a. Hamid. The crystal structure of the 11 Å natural tobermorite Ca2.25[Si3O7.5(OH)1.5].1H2O [J]. Zeitschrift für Kristallographie, 1981, 154: 189–198.
[12] B. P. Feuston, S. H. Garofalini. Oligomerization in silica sols [J]. J Phys Chem, 1990, 94: 5 351–5 356.
[13] F. H. Stillinger, A. Rahman. Revised centred force potentials for water [J]. J Chem Phys, 1978, 68(2): 666–670.
[14] W. Smith, T. R. Forester, I. L. Todorov, et al. The DL_POLY user manual. http://www.cse.clrc.ac.uk/msi/software/DL_POLY
