熔融氯化鈉晶體介電常數(shù)的模擬計算方法

陳 新 陳 衛(wèi) 張現(xiàn)仁 任 瑛 李糧生

(1.中國科學(xué)院過程工程研究所 多相復(fù)雜系統(tǒng)國家重點實驗室,北京 100190;2.北京化工大學(xué) 有機無機復(fù)合材料國家重點實驗室,北京 100029;3.北京化工大學(xué) 北京市能源環(huán)境催化重點實驗室,北京 100029;4.中國科學(xué)院大學(xué),北京 100049;5.電磁散射重點實驗室,北京 100854)

引 言

隨著微波技術(shù)的不斷發(fā)展,電介質(zhì)材料的介電常數(shù)作為描述材料電磁特性的重要參數(shù),反映了材料對外界電磁場的響應(yīng),目前已廣泛應(yīng)用于飛機隱身材料[1-2]、雷達探測[3-4]以及5G 通訊[5-6]等重要領(lǐng)域。介電弛豫光譜法[7-9]由于具有非侵入性、檢測迅速等優(yōu)點,常用于測定電介質(zhì)材料的介電常數(shù),但該方法往往受限于實際應(yīng)用中的復(fù)雜環(huán)境以及缺乏較為精確的參數(shù)模型。

近年來,隨著計算機技術(shù)的不斷發(fā)展,分子動力學(xué)模擬方法為學(xué)者在原子或分子尺度上理解材料的介電特性提供了重要幫助。除了可以有效預(yù)測實際工作情況下電介質(zhì)的介電特性,目前分子模擬方法已經(jīng)將頻率拓展到太赫茲(THz)光譜。Kirkwood[10]提出了介電常數(shù)的經(jīng)典統(tǒng)計理論,這是利用計算機模擬介電常數(shù)的基石;Neumann 等[11-13]基于線性響應(yīng)理論,利用分子動力學(xué)模擬將晶體的偶極矩漲落與周期性邊界條件下的頻率相關(guān)介電常數(shù)聯(lián)系起來,得到頻率相關(guān)介電常數(shù)的計算表達式;Caillol等[14-16]以離子-溶劑混合物體系為研究對象,根據(jù)電流和微觀極化的相關(guān)函數(shù),計算與頻率相關(guān)的介電常數(shù)和電導(dǎo)率,得到相對復(fù)雜的復(fù)介電常數(shù)表達式,該表達式包含偶極子-偶極子、電流-偶極子、電流-電流相關(guān)函數(shù)的貢獻;Thomas 等[17]提出了利用從頭計算分子動力學(xué)的方法計算二丁醇的磁偶極矩;Chen 等[18]利用機器學(xué)習(xí)結(jié)合分子動力學(xué)的方法提高了介電常數(shù)的計算精度;Yang 等[19]研究了高分子溶液中濃度變化導(dǎo)致的聚合物結(jié)構(gòu)轉(zhuǎn)變對體系介電常數(shù)的影響。目前,對電介質(zhì)材料介電常數(shù)的研究主要集中于固體材料和聚合物體系方面,但是針對體系中存在含凈電荷的離子、分子或熔融鹽晶體的計算機模擬介電常數(shù)的相關(guān)研究卻很少。這是由于傳統(tǒng)的偶極矩漲落方法無法適用于體系中存在含自由電荷的粒子的情況,其原因是:體系內(nèi)存在導(dǎo)電介質(zhì),無法維持靜電場平衡;周期性邊界對離子偶極矩計算會產(chǎn)生數(shù)值跳躍,這是因為偶極矩定義為粒子所帶電荷與位移矢量的乘積,當體系中施加周期性邊界條件時,對于類似水的這種極性分子,其具有完整的分子信息,周期性邊界條件并不會影響分子的偶極矩,然而對于含有自由電荷的體系,自由電荷可跨越周期性邊界條件,導(dǎo)致體系的偶極矩值不能收斂[20]。針對上述問題,Sega 等[21-22]提出了利用電導(dǎo)率計算NaCl 溶液介電頻譜的方法;Schr?der 等[23-25]致力于計算離子液體介電常數(shù)的研究,建立了計算離子液體介電常數(shù)和電導(dǎo)率的理論基礎(chǔ),研究分析了極化的微觀機制對離子液體廣義介電常數(shù)或磁化率的貢獻。

本文基于Neumann 介電理論公式[11-13],進一步推導(dǎo)了適用于離子溶液或熔融鹽晶體體系的介電常數(shù)解析表達式;以熔融氯化鈉(NaCl)晶體為研究算例,利用分子動力學(xué)方法測定了不同溫度下NaCl晶體的介電常數(shù)和頻率相關(guān)的介電譜,研究了NaCl晶體熔化過程中介電特性的變化;最后利用經(jīng)典外加電場法計算了熔融NaCl 晶體的靜態(tài)介電常數(shù),從而對本文的方法進行了驗證。

1 介電常數(shù)計算及分子模擬方法

1.1 介電常數(shù)計算

基于線性響應(yīng)理論,Neumann 等[11-13]將分子動力學(xué)模擬得到的晶體偶極矩漲落與周期性邊界條件下的頻率相關(guān)介電常數(shù)聯(lián)系起來,得到與頻率ω相關(guān)的經(jīng)典介電常數(shù)ε的計算表達式(式(1))。

式中:V為模擬體系的體積,kB為玻爾茲曼常數(shù),T為模擬體系溫度。

式中,φ(t)是總偶極矩Mtot的自相關(guān)函數(shù)。

考慮一個由N個帶電粒子組成、體積為V、溫度為T、粒子位移矢量為ri以及所帶電荷量為Qi(i=1,2,…,N)的體系,根據(jù)Schr?der 等[26]的研究,可以將體系的總偶極矩Mtot劃分為平動偶極矩MJ和質(zhì)心修正偶極矩MD,如公式(4)所示。

考慮到模擬系統(tǒng)中存在Nmol個分子,分子i包含ni個電荷,每個自由離子也被視為ni=1 的分子,因此分別表示分子i上的第j個電荷的電荷量和位置向量,則總偶極矩可以表示為

公式(5)中,等號右邊第一項表示質(zhì)心矯正的分子偶極矩μi,可以表示為

公式(5)中,等號右邊第二項是平動偶極矩MJ,

將公式(4)和(8)帶入公式(1)的傅里葉-拉普拉斯變換中,可以得到頻率相關(guān)介電常數(shù)的表達式(式(9))。

由于公式(9)考慮了體系中分子偶極貢獻與離子電流貢獻,故可廣泛作為體系的介電常數(shù)解析表達式。對于本文考慮的熔融NaCl 晶體體系,隨著溫度升高,NaCl 晶體結(jié)構(gòu)被破壞,體系中含有自由移動的離子電流,離子i的坐標位置等于其質(zhì)心坐標位置,那么體系中分子偶極矩μi=0,即體系中不存在分子偶極矩貢獻,故僅考慮體系中離子電流貢獻,其頻率相關(guān)介電常數(shù)的解析表達式為

當頻率小于103Hz 時,體系的介電常數(shù)稱為靜態(tài)介電常數(shù),利用電流法計算靜態(tài)介電常數(shù)εstatic,其表達式為

1.2 分子模擬方法

利用分子建模腳本Moltemplate[27]搭建了包含有1 000 個原子的NaCl 晶體的初始構(gòu)型,立方體模擬盒在3 個方向上的大小均為2.8 nm。對模擬盒的3 個方向都應(yīng)用周期性邊界條件,選用經(jīng)典原子間相互作用勢函數(shù)Born-Mayer -Huggins(BMH)勢[28-30]描述熔融NaCl 晶體離子之間的相互作用勢。

2月8日，徐云天讓弟弟給父親打電話，謊稱得了重病。徐河來了，見小兒子根本沒病，他轉(zhuǎn)身就要走。徐云天鼓足勇氣勸道：“爸，弟弟今年就要中考了，就算你跟我媽的緣份到頭了，也要替弟弟著想??！”

式中:uij為NaCl 晶體離子之間的相互作用勢,Aij為BMH 勢排斥參數(shù),Cij和Dij為范德華參數(shù),rij為粒子i和j之間的距離,ρij為離子對ij的長度參數(shù),qi和qj分別為i和j的電荷量,ε0為真空介電常數(shù)。

選用經(jīng)典分子模擬軟件Lammps 進行分子動力學(xué)模擬[31],選取隨機數(shù)種子為825 577,使體系中粒子的速度在給定的溫度下滿足玻爾茲曼分布;之后運行100 ps 的能量最小化過程,防止不合理的初始結(jié)構(gòu)導(dǎo)致模擬過程崩潰。采用NPT 系綜,分別使用Nosé-Hoover 熱浴法和Berendsen 耦合方法使體系維持給定的溫度和1 個標準大氣壓。使用Ewlad 求和法計算長程庫侖相互作用,截斷半徑為12 ?,采用Velocity-Verlet 算法求解運動方程的積分,時間步長設(shè)置為1 fs;每次進行MD 模擬過程時首先運行2 ns 以達到平衡狀態(tài),之后進行1 ns 的數(shù)據(jù)采樣過程,每隔0.005 ps 進行一次體系平均,從而計算靜態(tài)介電常數(shù)、頻率相關(guān)介電譜以及進行后處理。使用可視化軟件VMD[32]對軌跡文件進行觀察。本研究的模擬溫度區(qū)間為300 ～1 500 K,調(diào)節(jié)不同溫度大小是為了揭示NaCl 晶體的相變演化過程,研究相變過程中晶體材料介電特性的變化。NaCl 晶體的初始構(gòu)型由VMD 繪制[32],如圖1 所示。

圖1 NaCl 晶體的初始結(jié)構(gòu)模型Fig.1 Initial structure model of NaCl crystals

2 結(jié)果與討論

2.1 NaCl 晶體相變演化

在1 個標準大氣壓下,分別在300、500、700、900、1 100、1 180、1 200、1 220、1 240、1 300、1 500 K的溫度下計算體系的體積V,結(jié)果如圖2 所示?？梢钥闯?體系體積隨著溫度的升高而增加。當溫度升高至1 240K 時,體積開始快速膨脹,這表明固態(tài)NaCl 開始熔化為液態(tài),這與實驗測得的NaCl 的熔點(1 074 K)[33]基本吻合。當溫度達到1 300 K 時,體系已經(jīng)完全成為液態(tài)。

圖2 NaCl 的體積-溫度相圖Fig.2 Volume-temperature phase diagram of NaCl

此外,徑向分布函數(shù)(radial distribution function,RDF)[34]也可以作為分析體系結(jié)構(gòu)變化的重要參數(shù)用來研究物質(zhì)的有序性,其定義為一個系統(tǒng)在距離一個標記粒子的r處找到其他粒子的概率。圖3 為不同溫度下體系中Na-Cl 的RDF 圖,其中g(shù)(r)Na-Cl表示鈉離子周圍距離半徑r到r+dr處發(fā)現(xiàn)氯離子的概率。當溫度由1 180 K 升高至1 240 K時,RDF 在不同位置呈現(xiàn)出4 個離散峰,并且隨著溫度升高,g(r)Na-Cl幾乎沒有增加,最大峰值大約在3.4 左右,表示此時熔融態(tài)NaCl 晶體中仍有較強的約束作用,結(jié)構(gòu)只存在輕微擾動。隨著溫度升高到1 300 K 時,原來有序的熔融NaCl 晶體結(jié)構(gòu)被破壞,相互作用力減弱,分子結(jié)構(gòu)變得雜亂無序,RDF 峰值發(fā)生了顯著變化,并且僅保留兩個較為明顯的離散峰,表明體系已經(jīng)轉(zhuǎn)變?yōu)橐簯B(tài)。

圖3 不同溫度下Na-Cl 的徑向分布函數(shù)Fig.3 Radial distribution function of Na-Cl at different temperatures

2.2 介電頻譜

根據(jù)上文推導(dǎo)的介電常數(shù)的計算方法,為了獲得不同溫度下熔融NaCl 晶體體系的頻率相關(guān)性,首先需要獲得體系中離子電流-電流隨時間的變化規(guī)律。由于時間相關(guān)函數(shù)本身就是信號與自身之間相關(guān)性程度的度量,因此深入分析相關(guān)函數(shù)的行為有助于理解頻率相關(guān)介電譜的特性[21]。圖4 給出了1 180 ～1 300 K 的溫度下熔融NaCl 晶體的電流相關(guān)函數(shù)。可以看出,NaCl 電流自相關(guān)函數(shù)在約10-1ps處呈現(xiàn)極大的振蕩行為,隨后電流自相關(guān)函數(shù)逐漸衰減為零,表明該模擬體系在給定的時間尺度內(nèi)達到收斂。值得注意的是,當溫度升高至1 300 K 時,不僅電流相關(guān)函數(shù)迅速降低,而且峰值出現(xiàn)的頻率位置也發(fā)生了顯著變化,這可能表明隨著體系相態(tài)結(jié)構(gòu)的改變,共振頻率發(fā)生了變化。由于時間相關(guān)函數(shù)十分依賴于較長的模擬時間尺度,并且在分子模擬過程中,由于原子軌跡上的每個點只提供一個值,難以對每個粒子殘留物進行平均,這導(dǎo)致即使在長時間模擬下,無窮遠處的自相關(guān)函數(shù)仍然在零附近不停振蕩。此外,使用傅里葉-拉普拉斯變換的方法計算解析介電頻譜時,離散傅里葉變換算法將從分子動力學(xué)模擬中獲得的有限數(shù)據(jù)點進行周期性擴展,這并不是公式中簡單的連續(xù)函數(shù)積分。為了減少數(shù)學(xué)算法差異,往往需要更多的數(shù)據(jù)點,這增加了模擬的負擔,而離散傅里葉變換算法中使用的矩形積分方法將會引入顯著的誤差。因此,本文利用非線性最小二乘法將模擬數(shù)據(jù)擬合成特定的連續(xù)函數(shù)[35],并直接解析計算這些函數(shù)的傅里葉-拉普拉斯變換,這大大平滑了介電譜中的統(tǒng)計粗糙度。

圖4 不同溫度下NaCl 的電流相關(guān)函數(shù)Fig.4 Current correlation function of NaCl at different temperatures

由于電流自相關(guān)函數(shù)的行為類似于阻尼振蕩,通常在幾皮秒內(nèi)迅速振蕩,隨后逐漸衰減到零,因此可以使用多指數(shù)余弦相關(guān)函數(shù)來擬合相應(yīng)的模擬結(jié)果。

表1 電流相關(guān)函數(shù)的擬合參數(shù)Table 1 Fitting parameters of the current correlation function

圖5 不同溫度下NaCl 電流相關(guān)函數(shù)的擬合結(jié)果Fig.5 Fitting results of the current correlation function of NaCl at different temperatures

將電流擬合函數(shù)代入公式(10),利用Python腳本可以得到熔融NaCl 晶體的頻率相關(guān)介電譜。圖6 為對數(shù)坐標系中熔融NaCl 晶體在不同溫度下的頻率相關(guān)介電譜。在1 180 ～1 240 K 的溫度范圍內(nèi),隨著頻率的增加,大約在3 THz 處外場與NaCl 晶體的偶極子振蕩頻率發(fā)生共振,介電常數(shù)的實部εreal開始加速增大,達到峰值后快速減小,演化為接近于無外場作用下的極化狀態(tài)。同樣地,在達到共振頻率時,虛部介電損耗εimag從0 逐漸增大到峰值,隨后又減小到0。當溫度進一步升高到1 300 K 時,NaCl 開始熔化發(fā)生相變,與NaCl為晶體狀態(tài)時相比,熔融NaCl 晶體的介電譜發(fā)生了顯著變化。在NaCl 熔化以后,由于體系中存在不同組分,各組分的介電響應(yīng)相互重疊,介電譜的共振吸收峰的帶寬變寬,波形發(fā)生了明顯變化,導(dǎo)致體系的介電頻譜十分復(fù)雜。可以發(fā)現(xiàn),此時介電損耗峰值降低到約1.5,而實部處的共振峰值幾乎消失。這是因為NaCl 晶體結(jié)構(gòu)隨體系溫度改變而發(fā)生了變化,體系中離子組分隨頻率增加發(fā)生離子極化,分子的固有偶極矩的轉(zhuǎn)向極化不再響應(yīng)。

圖6 不同溫度下NaCl 的頻率相關(guān)介電譜Fig.6 Frequency dependent dielectric spectra of NaCl at different temperatures

2.3 外加電場分子動力學(xué)模擬

為了進一步驗證上述推導(dǎo)的介電常數(shù)計算的準確性,本文使用外加電場法計算靜態(tài)介電常數(shù)。所謂靜態(tài)介電常數(shù)是指當頻率很低或接近零時,介電常數(shù)是與電場頻率無關(guān)的量。利用外加電場法計算體系的靜態(tài)介電常數(shù)εstatic,其表達式為

式中:E為外電場強度,ME為施加電場后介質(zhì)的偶極矩,V為模擬體系的體積。因此,可以利用分子動力學(xué)模擬的方法在平衡態(tài)的NaCl 分子構(gòu)型的Z方向上施加靜電場EZ,電場強度以0.01 V/? 為間隔,從0.01 V/? 逐漸增大至0.1 V/?。然而當T=1 100 K 時,如果繼續(xù)向體系中施加0.1 V/? 的電場強度,此時NaCl 晶體的偶極矩大小受電場強度的影響會發(fā)生劇烈波動,這會無法保證電場強度與誘導(dǎo)偶極矩之間高度一致的線性關(guān)系,從而影響靜態(tài)介電常數(shù)的估算。因此本文選擇從0.001 V/? 逐漸增大至0.01 V/?。根據(jù)電場誘導(dǎo)偶極矩和電場強度之間的線性關(guān)系,利用最小二乘法擬合可得到不同溫度下NaCl 晶體的靜態(tài)介電常數(shù)。

圖7 為誘導(dǎo)偶極矩和電場強度之間的線性擬合結(jié)果?？梢园l(fā)現(xiàn),隨著體系施加的電場強度的增加,誘導(dǎo)偶極矩增大。在300 ～1 100 K 的溫度范圍內(nèi),電場強度和誘發(fā)偶極矩之間保持著良好的線性關(guān)系。隨著溫度進一步升高至1 300 K,NaCl 晶體開始熔化,采樣點的波動明顯增大,這可能與NaCl 晶體結(jié)構(gòu)被破壞有關(guān)。

圖7 不同溫度下電場誘導(dǎo)偶極矩MZ和電場強度EZ之間的線性關(guān)系Fig.7 Linear relationship between the electric field induced dipole moment and the electric field intensity at different temperatures

分別采用電流法和外加電場法計算NaCl 的靜態(tài)介電常數(shù)與溫度的相關(guān)性,結(jié)果如圖8 所示?？梢钥闯?當NaCl 為固態(tài)晶體時,靜態(tài)介電常數(shù)隨著溫度的升高而逐漸增加;當溫度升高到固液臨界點(T=1 100 K)時,靜態(tài)介電常數(shù)達到最大;隨著溫度進一步升高,NaCl 發(fā)生固-液相變過程,體系的靜態(tài)介電常數(shù)迅速下降。電流法和外加電場法的計算結(jié)果有著良好的一致性,電流法能夠有效反映NaCl 體系相變過程中靜態(tài)介電常數(shù)的變化規(guī)律,廣泛適用于熔融鹽晶體以及溶液中含有自由電荷的體系。根據(jù)文獻[36]的實驗測定結(jié)果,室溫下NaCl 晶體的靜態(tài)介電常數(shù)在5.45 ～5.9 之間。本文用電流法計算得到的靜態(tài)介電常數(shù)為4.5(T=300 K),小于實驗觀測值,這可能與描述NaCl 晶體的力場參數(shù)有關(guān)。在今后的研究中,本文將采用機器學(xué)習(xí)結(jié)合量子力學(xué)第一性原理的計算方式來提高力場的精度。此外,由公式(11)可知,對于電流法,如果要準確計算靜態(tài)介電常數(shù),需要從0 積分到無窮大,而實際上本文首先用最小二乘法將電流關(guān)聯(lián)函數(shù)擬合成一個解析函數(shù),然后對解析函數(shù)進行積分,這樣會引入誤差,導(dǎo)致計算結(jié)果偏低。

圖8 采用電流法和外加電場法計算得到的NaCl靜態(tài)介電常數(shù)與溫度的相關(guān)性Fig.8 Correlation between the static dielectric constant of NaCl and the temperature calculated by the current method and the external electric field method

3 結(jié)論

(1)基于經(jīng)典介電理論公式,進一步推導(dǎo)了適用于熔融鹽晶體以及溶液中含有自由電荷的體系的介電常數(shù)表達式,該公式可以用于計算熔融NaCl 晶體的介電性質(zhì)。

(2)利用電流擬合相關(guān)函數(shù)分析了不同溫度下NaCl 的頻率相關(guān)介電譜,結(jié)果表明:隨著溫度升高(300 ～1 100 K),在共振吸收峰處,介電常數(shù)的實部先增大,達到峰值后快速減小;當溫度升高至1 300 K 時,NaCl 晶體發(fā)生相變,介電損耗顯著下降,在共振頻率處介電常數(shù)的實部峰值幾乎消失。

(3)分別利用電流法和外加電場法計算了不同溫度下NaCl 的靜態(tài)介電常數(shù),結(jié)果顯示:隨著體系溫度升高(300 ～1 100 K),NaCl 晶體的介電常數(shù)呈上升趨勢;當溫度高于相變點時,NaCl 的靜態(tài)介電常數(shù)快速下降。電流法和外加電場法的計算結(jié)果有著良好的一致性,表明電流法可以有效描述NaCl 的相變過程。