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

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

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

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

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

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

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

1 螢石結構與研究晶面確定

1.1 螢石晶胞結構

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

1.2 研究晶面的確定

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

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

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

2 模擬方法

2.1 計算模塊的選擇

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

2.2 晶胞優(yōu)化

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

2.3 礦物表面模型與DFT計算

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

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

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

2.4 NaOL溶液模型

借助MS的Amorphous Cell模塊構建NaOL溶液模型，如圖4所示。力場設置為COMPASSⅡ，采用力場分配電荷平衡法（force-field-assigned）計算原子電荷。NaOL溶液模型的尺寸設置為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 分子動力學模擬細節(jié)

將優(yōu)化好的螢石單胞切割，建立（111）面、（110）面、（311）面和（100）面，使用supercell將UV值設置為10×10、7×10、6×10和10×10的表面單元，算法采用smart，收斂水平為10-4kcal/mol，然后對每個晶面進行幾何優(yōu)化以確定底部原子。在求和方法里將van der Waals設置為Atom-based，electrostatic設置為Ewald，Atom-based求和的截止距離為1.55 nm，采用力場分配電荷平衡法計算原子電荷。為結合NaOL溶液箱，在優(yōu)化表面的基礎上，建立真空層厚度為0的螢石表面模型。按照與之前相同的參數(shù)設置，再次對由礦物表面和NaOL溶液箱組成的復合模型進行幾何優(yōu)化。為獲得藥劑與礦物表面吸附平衡狀態(tài)下的最佳位置，進行MD模擬以使能量收斂。在MD模擬之前，（111）面、（110）面、（311）面和（100）面的礦物-試劑復合物的最終模型尺寸分別設置為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設置為Velocity Scale，時間步長為1 fs，模擬時間為200 ps，使得整個系統(tǒng)的能量達到平衡。接著將thermostat設置為Nose，Q比率為0.01，溫度為298 K，時間步長為1 fs，總模擬時間為1 000 ps。礦物表面與藥劑相互作用的相對親和力由相互作用能表示，由下式計算［18］：

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

3 結果與討論

3.1 螢石表面晶體特性分析

礦物的晶面性質受其晶體化學的影響，特別是對于細小的礦物顆粒。隨著礦物顆粒粒徑的減小會暴露更多的晶面，圖5為螢石不同晶面的結構及其斷裂鍵情況。由晶體的周期性可知在（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位點吸附結合，礦物晶面上的Ca位點越多其活性越強，與藥劑的反應程度越高［5］，由此可得知（111）和（100）面活性較強。

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

3.2 相互作用能比較

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

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

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

3.3 吸附構型性質

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

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

為了進一步闡明吸附構型，計算了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ī)律。上述結果表明NaOL在螢石（111）、（311）和（100）面發(fā)生穩(wěn)定吸附；螢石（110）面附近沒有C和O原子分布，而且C和O原子的濃度及其與螢石（110）面間的距離在吸附前后沒有發(fā)生較明顯的變化，結合相互作用能與吸附構型的分析，表明NaOL幾乎不與（110）面發(fā)生作用。

4 結論

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

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

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