螢石的晶體化學(xué)基因特性量化計算與分子動力學(xué)模擬

史新章 曹 釗，2，3，4 張金山 王介良 曹永丹 董紅娟

（1.內(nèi)蒙古科技大學(xué)礦業(yè)研究院，內(nèi)蒙古包頭014010；2.廣東省資源綜合利用研究所，廣東廣州510650；3.稀有金屬分離與綜合利用國家重點實驗室，廣東廣州510650；4.內(nèi)蒙古自治區(qū)礦業(yè)工程重點實驗室，內(nèi)蒙古包頭014010）

氟（F）是新能源、新材料、國防、光電、化工等領(lǐng)域 必不可少的化學(xué)材料［1］。螢石（CaF2）作為氟的主要來源，其含氟量高達48.67%，是地殼中最重要的含氟礦物相，已被列為我國的戰(zhàn)略性礦產(chǎn)資源［2-3］。

浮選是礦物分離的主要方法之一，浮選藥劑在礦物顆粒暴露面上吸附的穩(wěn)定性差異決定了礦物分離的效果。螢石常與稀土礦、鎢礦、方解石等礦物伴生，脂肪酸類捕收劑是其浮選分離的常用藥劑［4］。由于螢石與稀土礦、方解石、白鎢礦等礦物的可浮性相近，其浮選分離是行業(yè)內(nèi)難題［5］。因此，研究螢石的晶體化學(xué)、晶面性質(zhì)以及藥劑在螢石不同晶面上的吸附構(gòu)型和吸附能等礦物基因特性對研發(fā)新型浮選藥劑、實現(xiàn)螢石高效浮選分離有重要的指導(dǎo)意義。

分子模擬方法可以從原子角度去研究礦物晶體和晶面性質(zhì)，預(yù)測藥劑分子在礦物表面的吸附構(gòu)型和作用方式，定量計算藥劑在礦物表面作用過程中的能量及電荷變化等基因特性［6-8］。與傳統(tǒng)的浮選試驗相比，分子模擬技術(shù)具有耗時短、成本低、易操作等優(yōu)點，已經(jīng)成為礦物加工與分選領(lǐng)域工作者研究吸附機理、設(shè)計開發(fā)新型浮選藥劑分子的重要手段之一［9］。Gao等［5，10］基于密度函數(shù)理論（DFT）計算了幾種螢石晶面的表面能和表面斷裂鍵密度等礦物基因特性，發(fā)現(xiàn)表面能與螢石的斷裂鍵密度呈線性相關(guān)。Zhu等［6］憑借量化計算進一步探究了在石英（101）面上新型石英捕收劑α-BLA（α溴代月桂酸）的作用機制，觀察到Ca（OH）+的活化影響是因為其在石英表面與α-BLA間發(fā)揮了“橋聯(lián)”作用。Maldonado［11］等利用第一性原理計算模擬預(yù)測了螢石的表面穩(wěn)定性，發(fā)現(xiàn)螢石（111）面表面能最低，是螢石最穩(wěn)定的晶面。前人研究成果表明，模擬結(jié)果與實際微浮選試驗結(jié)果基本一致，證明了分子模擬技術(shù)研究礦物晶體化學(xué)基因特性的有效性與實用性。

螢石的浮選分離通常在水溶液中進行，水在NaOL水解和螢石顆粒的溶解中起重要作用。而大多數(shù)模擬研究在真空系統(tǒng)中進行，水溶液對液固界面的影響被簡化甚至忽略，導(dǎo)致模擬的結(jié)果不具說服力。本文利用DFT計算了螢石晶體的（111）面、（110）面、（311）面和（100）面4個常見解理面和暴露面的斷裂鍵密度，對螢石的晶體化學(xué)、晶面性質(zhì)等基因特性進行表征；通過在水溶液系統(tǒng)中的分子動力學(xué)（MD）模擬，研究并闡述了NaOL在不同螢石晶面上的吸附行為以及相關(guān)的吸附構(gòu)型和相互作用能，結(jié)合相關(guān)的量化計算以期探明NaOL在螢石晶體不同暴露面上的吸附機理。

1 螢石結(jié)構(gòu)與研究晶面確定

1.1 螢石晶胞結(jié)構(gòu)

螢石的晶胞結(jié)構(gòu)和參數(shù)由美國礦物學(xué)家晶體結(jié)構(gòu)數(shù)據(jù)庫獲得［12］，其晶胞參數(shù)：a=b=c=0.546 31 nm，α=β=γ=90°，空間群為Fm3 m，Z=4，屬于立方晶系。螢石的晶胞結(jié)構(gòu)如圖1所示。Ca2+在螢石晶體結(jié)構(gòu)中以立方緊密方式堆積，四面體的空隙位置被F-占據(jù)。Ca2+與附近的8個F-結(jié)合成8配位，形成Ca-F8立方體，F(xiàn)-與附近的4個Ca2+結(jié)合成4配位，形成正四面體［13］。Ca—F鍵的鍵長為0.236 6 nm。Ca—F鍵的離子鍵分?jǐn)?shù)為89.14%，化學(xué)鍵的離子鍵分?jǐn)?shù)大于50%，可認(rèn)為該化學(xué)鍵為離子鍵，因此Ca—F鍵為強離子鍵［5］。

1.2 研究晶面的確定

根據(jù)螢石晶體的X射線衍射（XRD）譜圖（圖2）確定主要強峰，找出這些主要強峰對應(yīng)的晶面。衍射峰的前 4個強度對應(yīng)的（111）、（220）、（311）和（400）面更容易通過機械粉碎和研磨獲得。

附著能（AE）模型是基于Hartman和Perdok提出的PBC理論［14-16］。運用PBC理論可以確定晶體薄片的附著能Eatt，晶面的附著能越低，其生長速率越小，因而越容易成為晶體形貌中最常見的晶面［17］。本文基于Material Studio 2017軟件的Morphology模塊計算了螢石晶體的形貌，計算結(jié)果如圖3和表1所示。

由圖3和表1可知，{111}晶面族的晶面附著能最小，晶面占比為100%，結(jié)合XRD測試結(jié)果，最終確定（111）、（110）、（311）和（100）面為可研究晶面（（110）面，即（220）面的平行晶面；（100）面，即（400）與（200）面的平行晶面）。

2 模擬方法

2.1 計算模塊的選擇

利用Materials Studio 2017軟件中的CASTEP和Forcite模塊進行相關(guān)計算。CASTEP模塊用于密度泛函計算，包括螢石晶體結(jié)構(gòu)與表面模型的優(yōu)化和能量計算等任務(wù)。Forcite模塊用于油酸鈉與螢石不同晶面的分子動力學(xué)模擬，該模擬在模塊中的COMPASSⅡ（原子模擬研究的凝聚相優(yōu)化分子勢）力場下進行。COMPASSⅡ力場包括常見的有機分子和無機材料，包括金屬、金屬鹵化物、鋁硅酸鹽和金屬氧化物等［18］，這使其可以合理地應(yīng)用于構(gòu)建的模型。

2.2 晶胞優(yōu)化

通過MS2017中的Build Crysatls模塊對螢石晶胞建模，然后采用CASTEP模塊對其進行幾何優(yōu)化，其中交換相關(guān)泛函選擇廣義梯度近似Wu-Cohen（WC）方法［19］。幾何優(yōu)化的收斂參數(shù)設(shè)置為：最大原子間作用力為0.3 eV/nm，最大原子位移為1×10-4nm，各個原子的能量收斂誤差小于1×10-5eV，最大原子內(nèi)應(yīng)力0.05 GPa。SCF收斂精度設(shè)置為Fine，其收斂誤差不超過1×10-6eV，贗勢采用OTFG Ultrasoft，平面波截斷能Energy cutoff設(shè)置為394.6 eV，k點取樣密度為6×6×6，k點取樣間隔separation在幾何優(yōu)化過程中設(shè)置為4×10-3/nm。上述參數(shù)在后續(xù)計算中若無特殊情況不做更改。

優(yōu)化后的晶胞參數(shù)為：a=b=c=0.551 98 nm，α=β=γ=90°，空間群為Fm3m。將優(yōu)化后的晶胞參數(shù)與實驗值對比發(fā)現(xiàn)誤差小于3%，表明參數(shù)設(shè)置合理，該模型可以用于后續(xù)計算。

2.3 礦物表面模型與DFT計算

將優(yōu)化后的晶胞通過supercell建立2×1的超晶胞，沿著不同的方向?qū)ζ淝懈睿纬桑?11）面、（110）面、（311）面以及（100）面，在不同晶面上方建立真空層得到各自的表面結(jié)構(gòu)。隨后在CASTEP模塊中選擇WC泛函計算螢石不同晶面的表面能。表面層和真空層厚度會對表面結(jié)構(gòu)的穩(wěn)定性造成影響，其穩(wěn)定性可由表面能表征，越低的表面能代表著表面結(jié)構(gòu)穩(wěn)定性越好［6］。當(dāng)表面層和真空層厚度繼續(xù)增加時，表面能會趨于穩(wěn)定，所以對其進行收斂性測試為后續(xù)計算構(gòu)建最穩(wěn)定的表面模型［20］。

分別采用式（1）和式（2）計算不同晶面的斷裂鍵密度和表面能［21］。

其中，Db和Nb分別為晶面斷裂鍵密度和某晶面單胞區(qū)域內(nèi)的未飽和鍵（斷裂鍵）個數(shù)；S為該晶面上單胞區(qū)域的面積；Eslab和Ebulk分別為表面結(jié)構(gòu)和單胞的總能量；Nslab和Nbulk分別為表面結(jié)構(gòu)和單胞的總原子數(shù)；A為表面結(jié)構(gòu)沿Z軸方向的面積；2表示表面結(jié)構(gòu)沿Z軸方向有上下2個表面。

2.4 NaOL溶液模型

借助MS的Amorphous Cell模塊構(gòu)建NaOL溶液模型，如圖4所示。力場設(shè)置為COMPASSⅡ，采用力場分配電荷平衡法（force-field-assigned）計算原子電荷。NaOL溶液模型的尺寸設(shè)置為3.9 nm×3.9 nm×3.6 nm、3.8 nm×3.9 nm×3.6 nm、4.1 nm×3.9 nm×3.4 nm、3.9 nm×3.9 nm×3.5 nm。

2.5 分子動力學(xué)模擬細(xì)節(jié)

將優(yōu)化好的螢石單胞切割，建立（111）面、（110）面、（311）面和（100）面，使用supercell將UV值設(shè)置為10×10、7×10、6×10和10×10的表面單元，算法采用smart，收斂水平為10-4kcal/mol，然后對每個晶面進行幾何優(yōu)化以確定底部原子。在求和方法里將van der Waals設(shè)置為Atom-based，electrostatic設(shè)置為Ewald，Atom-based求和的截止距離為1.55 nm，采用力場分配電荷平衡法計算原子電荷。為結(jié)合NaOL溶液箱，在優(yōu)化表面的基礎(chǔ)上，建立真空層厚度為0的螢石表面模型。按照與之前相同的參數(shù)設(shè)置，再次對由礦物表面和NaOL溶液箱組成的復(fù)合模型進行幾何優(yōu)化。為獲得藥劑與礦物表面吸附平衡狀態(tài)下的最佳位置，進行MD模擬以使能量收斂。在MD模擬之前，（111）面、（110）面、（311）面和（100）面的礦物-試劑復(fù)合物的最終模型尺寸分別設(shè)置為3.9 nm×3.9 nm×14.7 nm（總原子數(shù)：7 380）、3.8 nm×3.9 nm×14.4 nm（總原子數(shù)：7 680）、4.1 nm×3.9 nm×13.9 nm（總原子數(shù)：7 740）和3.9 nm×3.9 nm×14.2 nm（總原子數(shù)：7 980）。

MD模擬使用NVT系綜（恒定體積和恒定溫度），thermostat設(shè)置為Velocity Scale，時間步長為1 fs，模擬時間為200 ps，使得整個系統(tǒng)的能量達到平衡。接著將thermostat設(shè)置為Nose，Q比率為0.01，溫度為298 K，時間步長為1 fs，總模擬時間為1 000 ps。礦物表面與藥劑相互作用的相對親和力由相互作用能表示，由下式計算［18］：

式中，Ecom為優(yōu)化后復(fù)合模型的總能量；Erea+water+ion和Esur+water+ion分別表示沒有礦物表面的系統(tǒng)和沒有試劑的系統(tǒng)能量；Ewater+ion是鈉離子和水分子的能量。相互作用能的值（ΔE）表示吸附系統(tǒng)的穩(wěn)定性。相互作用能值越負(fù)表明礦物表面和藥劑之間的相互作用越強，吸附過程可以自發(fā)地進行，而作用能為正則意味著藥劑與礦物表面的相互作用較弱，吸附過程不能自發(fā)地進行?；贛D模擬的結(jié)果，使用Forcite模塊進一步計算了沿晶體表面Z軸NaOL中代表性原子（C和O原子）的濃度分布，以闡明吸附構(gòu)型。

3 結(jié)果與討論

3.1 螢石表面晶體特性分析

礦物的晶面性質(zhì)受其晶體化學(xué)的影響，特別是對于細(xì)小的礦物顆粒。隨著礦物顆粒粒徑的減小會暴露更多的晶面，圖5為螢石不同晶面的結(jié)構(gòu)及其斷裂鍵情況。由晶體的周期性可知在（111）面上存在不飽和的Ca原子，其具有4個Ca—F斷裂鍵。除了Ca—F斷裂鍵以外，在（110）面上還觀察到有2個不飽和的F原子，其具有2個F—Ca斷裂鍵。在（311）面上觀察到1個不飽和的F原子，其具有1個F—Ca斷裂鍵。（100）面的斷裂鍵情況與（111）面相似，也具有4個Ca—F斷裂鍵。NaOL主要是與螢石晶面的Ca位點吸附結(jié)合，礦物晶面上的Ca位點越多其活性越強，與藥劑的反應(yīng)程度越高［5］，由此可得知（111）和（100）面活性較強。

表2和表3為螢石不同晶面的表面能和斷裂鍵密度的計算結(jié)果。由表2和表3可知，螢石不同晶面的斷裂鍵密度（Db）按以下順序遞減：（111）面＞（100）＞（311）面＞（110）面。通過分析表面能、斷裂鍵密度和面間距等礦物基因特性可知，雖然（100）和（111）面都具有100%的Ca—F斷裂鍵，但是（100）面的面間距小于（111）面，其表面能大于（111）面，因此（111）面為螢石的最穩(wěn)定解理面；（100）面由于其表面能較高不易暴露，（110）面比（311）面更容易暴露，上述結(jié)論與前人的研究結(jié)果一致［5，11］。

3.2 相互作用能比較

圖6為NaOL在不同螢石晶面上吸附前后的復(fù)合模型。NaOL捕收劑被溶解成油酸陰離子（去質(zhì)子化的RCOO—基團）和鈉離子（Na+）。由圖6可知，在經(jīng)過MD模擬后，幾乎所有的油酸陰離子都吸附在（111）、（311）和（100）面上，這些晶面的水膜被破壞，其中（111）面的吸附效果最為顯著，而螢石（110）面幾乎沒有油酸陰離子吸附。

藥劑在礦物表面上的相互作用能越負(fù)，藥劑與礦物表面的親和力越高，表4為NaOL與不同螢石晶面的相互作用能。

由表4可知，NaOL與螢石（100）面的相互作用能和（111）面接近，但比（311）面的要負(fù)得多，NaOL幾乎不與（110）面發(fā)生吸附作用。結(jié)合前文中的有關(guān)結(jié)論，NaOL更傾向于與最穩(wěn)定的（111）面吸附結(jié)合。

3.3 吸附構(gòu)型性質(zhì)

NaOL在不同螢石晶面上的吸附構(gòu)型如圖7所示。

由圖7可知：幾乎所有的NaOL以雙核雙配位構(gòu)型與（111）面上的Ca位點結(jié)合，這是油酸鹽與含鈣礦物相互作用的最穩(wěn)定構(gòu)型［5］；螢石（110）面有50%的F—Ca斷裂鍵，大量水分子通過與F原子形成氫鍵吸附在（110）面導(dǎo)致其被較厚的水膜覆蓋，高親水性使得NaOL在（110）面的吸附變得困難；在NaOL與（311）面的配合物結(jié)構(gòu)中發(fā)現(xiàn)了單核單配位、雙核單配位和雙核雙配位3種吸附構(gòu)型，由于其表面有17%的F—Ca斷裂鍵通過與水分作用形成了氫鍵，水分子比（111）面更接近（311）面；NaOL與（100）面的配合物結(jié)構(gòu)中發(fā)現(xiàn)單核單配位和雙核雙配位2種吸附構(gòu)型，其中雙核雙配位構(gòu)型占主導(dǎo)地位。

為了進一步闡明吸附構(gòu)型，計算了NaOL中代表性原子（C和O原子）沿晶面Z軸的濃度分布，以定量說明NaOL在螢石不同晶面上的分布情況，如圖8所示。

由圖8可知：NaOL在吸附后，C和O原子與螢石晶面間的距離變化明顯，并在螢石（111）、（311）和（100）面上有規(guī)律地分布；其中C原子分布在距螢石（111）面約0.27～2 nm的范圍內(nèi)，而O原子分布在距（111）面約0.20～0.4 nm的范圍內(nèi)，最強的原子鍵長為0.20 nm，在（311）和（100）面也發(fā)現(xiàn)了同樣的規(guī)律。上述結(jié)果表明NaOL在螢石（111）、（311）和（100）面發(fā)生穩(wěn)定吸附；螢石（110）面附近沒有C和O原子分布，而且C和O原子的濃度及其與螢石（110）面間的距離在吸附前后沒有發(fā)生較明顯的變化，結(jié)合相互作用能與吸附構(gòu)型的分析，表明NaOL幾乎不與（110）面發(fā)生作用。

4 結(jié)論

（1）通過X射線衍射分析與MS晶體形貌計算，確定了螢石的4個常見暴露晶面，分別為（111）、（110）、（311）和（100）面；通過對螢石不同晶面的表面能比較，發(fā)現(xiàn)（111）面的表面能最小，表明（111）面是螢石最穩(wěn)定的解理面。

（2）通過DFT計算表面斷裂鍵分析，螢石晶面的各向異性與其原子種類與排布有關(guān)。（111）和（100）面具有100%的Ca—F斷裂鍵密度，對NaOL有很強的親和力，吸附作用較強；（110）面和（311）面隨著螢石被粉碎并研磨成更細(xì)小的顆粒時會暴露出來，由于其分別含有50%和17%的F—Ca斷裂鍵而親水，從而降低對NaOL的親和力，吸附作用較弱。

（3）通過相互作用能以及原子濃度分布分析，OL-陰離子與（111）面的相互作用能絕對值最大，以雙核雙配位構(gòu)型優(yōu)先與其相互作用并產(chǎn)生穩(wěn)定吸附；OL-陰離子與（110）面的相互作用能絕對值最小，NaOL幾乎不與其發(fā)生作用。