單晶鋅電極在應(yīng)變電荷耦合下的溶解機(jī)制-跨尺度計(jì)算與機(jī)器學(xué)習(xí)

蔣紅星 孫俁哲 孫 升 張統(tǒng)一

（上海大學(xué)材料基因組工程研究院，上海 200444）

20世紀(jì)以來(lái)，金屬材料的使用量急劇增長(zhǎng)，其開(kāi)采冶煉加工過(guò)程消耗了大量的能源，每年因?yàn)殡娀瘜W(xué)腐蝕造成了大量金屬材料的浪費(fèi)，所以金屬的腐蝕溶解問(wèn)題得到了廣泛的關(guān)注。金屬電極的腐蝕受應(yīng)力應(yīng)變和電化學(xué)環(huán)境的影響，因此有必要研究應(yīng)變-電荷耦合場(chǎng)下金屬電極的腐蝕特性。Gutman等［1］研究發(fā)現(xiàn)，外加力場(chǎng)會(huì)導(dǎo)致金屬材料電化學(xué)勢(shì)的改變，在金屬電極發(fā)生彈性變形時(shí)，隨著拉伸／壓縮應(yīng)力的增加，金屬的開(kāi)路電壓呈線性降低的趨勢(shì)。這是由于力學(xué)化學(xué)作用使應(yīng)變能轉(zhuǎn)變?yōu)殡娀瘜W(xué)能，即金屬的表面活性增大，腐蝕傾向也增大。國(guó)內(nèi)關(guān)于應(yīng)變腐蝕的研究報(bào)道［2-3］也很多，如 Zhang等［3］研究了拉應(yīng)力作用下X80鋼在酸性溶液中的電化學(xué)行為，發(fā)現(xiàn)隨著拉應(yīng)力的增加，電極的開(kāi)路電位降低，反應(yīng)活性增加。Mavrikakis等［4］通過(guò)密度泛函理論（density functional theory，DFT）研究了應(yīng)變對(duì)金屬電極表面反應(yīng)活性的影響，計(jì)算中對(duì)Ru（0001）晶面施加面內(nèi)應(yīng)變，發(fā)現(xiàn)表面應(yīng)變與O原子／CO分子的吸附能和活化能壘之間存在關(guān)聯(lián)性，相對(duì)于無(wú)應(yīng)變表面，拉伸應(yīng)變表面降低了對(duì)O原子／CO分子的吸附能以及CO分解過(guò)程的能壘，并將這種效應(yīng)歸因于金屬d能帶中心隨應(yīng)變的變化。

鋅具有蘊(yùn)藏量豐富、低毒、無(wú)污染、高能量密度以及低成本等特點(diǎn)，而被廣泛用作電池負(fù)極材料。其中可充電鋅離子電池（zinc-ion battery，ZIB）由于其良好的環(huán)境兼容性、極高容量密度、低廉的價(jià)格和很高的安全性，被認(rèn)為是今后大規(guī)模能量?jī)?chǔ)存的方向之一［5-6］。鍍鋅是一種犧牲陽(yáng)極保護(hù)陰極的電化學(xué)防腐蝕方法，其工藝成熟、成本較低，被廣泛應(yīng)用于制造業(yè)［7］。

但目前的試驗(yàn)方法無(wú)法從微觀層面研究腐蝕過(guò)程中原子的溶解和電荷的轉(zhuǎn)移過(guò)程，而跨尺度計(jì)算則可以通過(guò)原子計(jì)算結(jié)合宏觀溶液計(jì)算較好地模擬鋅的電化學(xué)溶解。Schwarz等［8］通過(guò)使用耦合密度泛函理論／連續(xù)介質(zhì)溶液方法，在寬電壓范圍內(nèi)的密集計(jì)算，揭示了吸附在金屬晶體臺(tái)階邊緣上水的部分氧化機(jī)制。本文選擇密度泛函／連續(xù)介質(zhì)溶液方法，對(duì)鋅電極的溶解腐蝕過(guò)程進(jìn)行了詳細(xì)研究，并結(jié)合符號(hào)回歸算法預(yù)測(cè)了鋅電極發(fā)生溶解腐蝕的臨界電極表面應(yīng)變-電荷組合。

1 數(shù)學(xué)建模與方法

圖1 電極／溶液模型和鋅電極／溶質(zhì)鋅原子結(jié)構(gòu)示意圖Fig.1 Schematic of the electrode／solution model and zinc electrode／solute zinc atoms

計(jì)算系統(tǒng)在整個(gè)電化學(xué)回路中的位置如圖1所示，電極為Zn，并浸入含NaF的電解質(zhì)溶液中。電極使用DFT計(jì)算，電解質(zhì)溶液采用連續(xù)介質(zhì)模型計(jì)算［9］，計(jì)算采用超胞方法。電極模型采用了22個(gè)Zn原子，其中20個(gè)原子形成一個(gè)含有5層原子的平板電極（每層有2×2個(gè)Zn原子），另2個(gè)鋅溶質(zhì)原子分布在平板的兩側(cè)，形成對(duì)稱結(jié)構(gòu)，如圖1中的插圖所示。Zn電極與電解質(zhì)溶液接觸的表面方向?yàn)槊芘帕骄Ц竦模?001＞方向，［0001］方向被定義為 z軸，Zn電極最外層的原子面則位于z＝0。電極兩側(cè)填充電解質(zhì)溶液，超胞沿z向的長(zhǎng)度被設(shè)定為44.8?，電極厚度為10.28?。

本文使用開(kāi)源軟件 JDFTx［10-13］進(jìn)行密度泛函／連續(xù)介質(zhì)溶液耦合計(jì)算，電解質(zhì)溶液模型采用SaLSA模型，該模型具有非局域介電響應(yīng)特性［9］，對(duì)電解質(zhì)溶液中強(qiáng)極化和帶電溶質(zhì)的計(jì)算的準(zhǔn)確度較高。使用濃度為1 M的隱性Na+和F-離子為鋅電極提供徳拜屏蔽效應(yīng)。DFT計(jì)算采用平面波展開(kāi)，平面波矢空間的K點(diǎn)網(wǎng)格劃分為12×12×2，截?cái)嗄芰吭O(shè)定為707 eV，并采用周期性邊界條件，交換關(guān)聯(lián)函數(shù)使用廣義梯度近似（general gradient approximation，GGA）［14］，鋅 原 子 使 用Ultrasoft贗勢(shì)方法［15］。電極的電子態(tài)密度計(jì)算使用線性四面體插值方法［16］。

計(jì)算前先進(jìn)行結(jié)構(gòu)優(yōu)化，計(jì)算步驟如下：（1）將建好的Zn電極模型去除2個(gè)溶質(zhì)原子后放入溶液中進(jìn)行晶格參數(shù)優(yōu)化（即整個(gè)結(jié)構(gòu)處于能量最小化），得到Zn電極HCP結(jié)構(gòu)的平均晶格參數(shù)為：；（2）將兩個(gè)溶質(zhì)Zn原子放入第1步的結(jié)構(gòu)兩側(cè)進(jìn)行能量最小化（第1步中的結(jié)構(gòu)保持固定，僅松弛2個(gè)溶質(zhì)原子的位置），溶質(zhì)原子會(huì)吸附在鋅電極表面3個(gè)相鄰Zn原子的幾何中心上方。在之后一系列計(jì)算中電極的晶格參數(shù)和溶質(zhì)原子在平行于電極表面內(nèi)的位置都保持不變，改變?nèi)苜|(zhì)Zn原子與電極表面的距離（從 z＝1．55?至 z＝14．65?）進(jìn)行一系列能量計(jì)算。在密度泛函計(jì)算中，金屬表面的應(yīng)變通常是通過(guò)改變表面方向的晶格參數(shù)來(lái)實(shí)現(xiàn)的［4，17-19］，本文也采用該方法，即通過(guò)同時(shí)改變晶格a和b方向的參數(shù)來(lái)施加面內(nèi)應(yīng)變。

DFT計(jì)算時(shí)將凈電荷施加在鋅電極上，因此凈電荷位于電極和2個(gè)孤立Zn原子中。表面電荷由電極最外表面原子的電荷密度相對(duì)于電極內(nèi)部原子的電荷密度的變化來(lái)定義，通常使用吉布斯突變表面模型，定義為表面每單位表面積的凈電荷。為了區(qū)分這兩類凈電荷，將整個(gè)計(jì)算系統(tǒng)的凈電荷簡(jiǎn)稱為凈電荷，而表面的凈電荷稱為表面電荷。凈電荷是由施加的電位引起的，而表面電荷是由自由表面引起的。表面電荷隨凈電荷的變化而變化。

密度泛函／連續(xù)介質(zhì)溶液計(jì)算中分別使用自洽場(chǎng)（SCF）和Polak-Ribière共軛梯度法進(jìn)行電子和電解質(zhì)的弛豫。在兩次連續(xù)的SCF自洽循環(huán)中，能量收斂的相對(duì)準(zhǔn)確度被設(shè)定為10-7Hartree（1 Hartree≈27.211 4 eV）。采用 L?wdin電荷布居分析［21-22］計(jì)算電極和溶質(zhì)Zn原子的電荷分布。

由于計(jì)算Zn溶解／吸附的臨界（凈電荷、電極面內(nèi)應(yīng)變）組合需花費(fèi)較多時(shí)間和成本，因此基于機(jī)器學(xué)習(xí)中的符號(hào)回歸算法（symbolic regression）［23-24］，找出擬合訓(xùn)練數(shù)據(jù)的函數(shù)關(guān)系式，然后利用函數(shù)關(guān)系式，求得具有給定性質(zhì)的參數(shù)組合，進(jìn)而加速尋找臨界凈電荷-電極表面應(yīng)變組合。

2 結(jié)果和討論

為了計(jì)算方便，采用22個(gè)鋅原子的凈電荷總值表征系統(tǒng)的凈電荷。圖2（a）和圖2（b）分別是在凈電荷為3.5｜e｜和4.0｜e｜時(shí)，不同電極面內(nèi)應(yīng)變下能量隨電極-溶質(zhì)原子距離的變化（圖中實(shí)線表示使鋅溶質(zhì)原子自動(dòng)溶解的應(yīng)變）。能量Ead的定義為：

式中：E［Zn（slab）+Zn（z→∞）］表示計(jì)算系統(tǒng)的能量，此時(shí)溶質(zhì)鋅原子遠(yuǎn)離電極表面，Zn（z→∞）］由 Zn（z＝14.7?）近似，在 z≥7?時(shí)溶質(zhì)鋅原子和Zn電極之間的相互作用可忽略不計(jì)。

如圖2（a）所示，在凈電荷為 3.5｜e｜、未施加應(yīng)變的條件下，在2.46?＜z＜4.5?范圍內(nèi)，系統(tǒng)能量處于1 eV的平臺(tái)，即在此范圍內(nèi)溶質(zhì)原子的溶解／吸附近似處于平衡狀態(tài)。在電極面內(nèi)施加5%拉伸應(yīng)變時(shí)，在1.6?＜z＜4.5?范圍內(nèi)出現(xiàn)了一段能量峰值，即從2.46?處的能量低點(diǎn)（電極表面的吸附位置）到4.5?處的能量高點(diǎn)，約0.17 eV的能量勢(shì)壘，該勢(shì)壘會(huì)阻礙溶質(zhì)原子溶解進(jìn)入電解質(zhì)溶液。在電極面內(nèi)施加 -5% 應(yīng)變（壓縮應(yīng)變）時(shí)，吸附能單調(diào)遞減，說(shuō)明表面溶質(zhì)原子產(chǎn)生了自發(fā)溶解。z＞4.5?時(shí)，不同應(yīng)變下系統(tǒng)能量均不斷降低，當(dāng)z＞6.9?時(shí)基本穩(wěn)定，表明溶質(zhì)原子在遠(yuǎn)離電極表面時(shí)與電極間的相互作用微弱，而在4.5?＜z＜6.9?時(shí)能量明顯降低，說(shuō)明溶液中的溶質(zhì)原子將難以跨越能壘而傳遞到電極上。如上所述，（3.5｜e｜，0）為一組臨界組合（電荷，應(yīng)變）。

如圖2（b）所示，在凈電荷為 4｜e｜、未施加應(yīng)變的條件下，系統(tǒng)能量隨溶質(zhì)原子-電極表面間距增加而單調(diào)減小，說(shuō)明表面溶質(zhì)原子產(chǎn)生了自發(fā)溶解，且這一過(guò)程不可逆。在電極面內(nèi)施加5%和15%拉伸應(yīng)變時(shí)，在2.46?＜z＜4.5?范圍內(nèi)，能量下降坡度隨著拉伸應(yīng)變的增加逐漸減緩，此處能量坡度定義為溶質(zhì)原子從z＝4.5?進(jìn)入到吸附位置z＝2.46?時(shí)的能量勢(shì)壘，該值在0%、5%和15%應(yīng)變下分別為1.02、0.63和0.11 eV。在施加20%拉伸應(yīng)變時(shí)，能量坡度為負(fù)值，在z＝2.46?時(shí)最低為0.61 eV，在z＝4.5?時(shí)最高為0.82 eV，能量勢(shì)壘為0.21 eV，阻礙溶質(zhì)原子進(jìn)入電解質(zhì)溶液。在z＞4.5?范圍內(nèi)，能量變化趨勢(shì)與圖2（a）的相似，4.5?＜z＜6.9?范圍內(nèi)的能量坡度阻止了溶質(zhì)原子吸附到電極表面。

圖2 不同凈電荷下能量隨電極面內(nèi)應(yīng)變的變化Fig.2 Energy profiles of zinc atoms on a zinc electrode under different surface strains at different net charges

圖3（a、b）和圖 3（c、d）分別是在凈電荷為3.5｜e｜和 4.0｜e｜時(shí)，不同電極面內(nèi)應(yīng)變下溶質(zhì) Zn原子和電極最外表面鋅原子所帶電荷隨電極-溶質(zhì)原子距離的變化。圖3（a）中，溶質(zhì)原子的最大電荷均出現(xiàn)在z＝4.4?時(shí)，在0%、5%、-5%應(yīng)變下該數(shù)值分別為 1.65｜e｜、1.53｜e｜和1.76｜e｜。在［1.6，4.4］?范圍內(nèi)，不同應(yīng)變下每個(gè)溶質(zhì)Zn原子所帶電荷均隨著電極-溶質(zhì)原子間距的增加而增加，如在5%拉伸應(yīng)變下，溶質(zhì)原子所帶電荷從1.11｜e｜增加到1.53｜e｜。圖3（b）中，在［1.6，4.4］?距離范圍內(nèi)，不同應(yīng)變下電極最外表面每個(gè)鋅原子所帶電荷均隨著電極-溶質(zhì)原子間距的增加單調(diào)減少，在0%應(yīng)變下，鋅原子所帶電荷從0.32｜e｜減小到 0.07｜e｜，5%應(yīng)變下從0.19｜e｜減小到 0.08｜e｜，在 -5%應(yīng)變下從 0.16｜e｜減小到0.03｜e｜。當(dāng)z＞6.9?時(shí)，在壓縮應(yīng)變和無(wú)應(yīng)變條件下，電極最外表面鋅原子所帶電荷接近，約0.09｜e｜，在拉伸應(yīng)變下，約 0.11｜e｜。

圖3（c）中，隨著拉伸應(yīng)變的增加，溶質(zhì)鋅原子所帶電荷降低，如在z＝6.9?，電極面內(nèi)拉伸應(yīng)變?yōu)?%、5%、15%和20%時(shí)，溶質(zhì)鋅原子所帶電荷分別為 1.89｜e｜、1.81｜e｜、1.63｜e｜和 1.53｜e｜。圖3（d）中，在［1.6，3.4］?范圍內(nèi)，0%應(yīng)變下，電極最外表面鋅原子所帶電荷最高，如在z＝2．56 ?時(shí)為0.23｜e｜，5%、15%、20%應(yīng)變下分別對(duì)應(yīng)為0.17｜e｜、0.18｜e｜、0.19｜e｜，該結(jié)果與圖 3（b）中相似，即應(yīng)變導(dǎo)致在［1.6，3.4］?范圍內(nèi)電極最外表面原子所帶電荷下降，但拉伸應(yīng)變的增加導(dǎo)致電荷升高，壓縮應(yīng)變的增加導(dǎo)致電荷下降。在z＞4．4?時(shí)，隨著拉伸應(yīng)變的增加，最外表面原子所帶電荷升高且隨著溶質(zhì)鋅原子與電極表面距離的增加基本不變，如在 z＝6.9?時(shí)，0%、5%、15% 和 20%下分別為 0.07｜e｜、0.10｜e｜、0.12｜e｜、0.13｜e｜。結(jié)合圖2與圖3（a）和3（c）中的實(shí)線部分，可見(jiàn)在3.5｜e｜凈電荷和-5%應(yīng)變、4.0｜e｜凈電荷和 0%應(yīng)變、4.0｜e｜凈電荷和 5%應(yīng)變、4.0｜e｜凈電荷和15%應(yīng)變條件下，溶質(zhì)鋅原子溶解進(jìn)入電解質(zhì)溶液所帶電荷分別為1.45｜e｜、1.89｜e｜、1.81｜e｜、1.63｜e｜，比較能量曲線發(fā)現(xiàn)，溶質(zhì)原子所帶電荷越接近經(jīng)典理論的+2價(jià)態(tài)，溶質(zhì)Zn原子越易溶解。

為了尋找凈電荷4｜e｜下溶質(zhì)鋅原子溶解的臨界應(yīng)變值，使用SR進(jìn)行符號(hào)回歸訓(xùn)練，圖4顯示了在凈電荷為4｜e｜時(shí)，不同溶質(zhì)原子-電極表面距離和拉伸應(yīng)變條件下系統(tǒng)能量的符號(hào)回歸數(shù)值模型。將訓(xùn)練過(guò)程中的4組數(shù)據(jù)（共40個(gè)，圖4中黑色數(shù)據(jù)點(diǎn)）均作為訓(xùn)練集，每個(gè)數(shù)據(jù)包含了溶質(zhì)原子-電極表面間距與電極面內(nèi)應(yīng)變2個(gè)特征量和1個(gè)系統(tǒng)能量輸出量。使用 gplearn（https：／／gplearn.readthedocs.io）中的符號(hào)回歸模塊對(duì)歸一化后的數(shù)據(jù)集進(jìn)行訓(xùn)練，符號(hào)回歸模型使用的參數(shù)有：population_size＝5000，迭代次數(shù)（generations）為 113，計(jì)算精度（stopping_criteria）設(shè)置為0.008，交叉概率（p_crossover）為0.7，子樹(shù)變異概率p_subtree_mutation＝0.1，hoist變異概率p_hoist_mutation＝0.05，點(diǎn)變異概率p_point_mutation＝0.1，最大樣本抽取 max_samples＝0.9，過(guò)程打印 verbose＝1，稀疏系數(shù) parsimony_coefficient＝0.000 12。根據(jù)圖4中符號(hào)回歸訓(xùn)練模型，在［0%，20%］應(yīng)變范圍內(nèi)等間距取100個(gè)應(yīng)變值并代入圖4的符號(hào)回歸模型，比較發(fā)現(xiàn)在應(yīng)變?yōu)?6.75%時(shí)，溶質(zhì)原子在電極表面的溶解／吸附最接近于臨界狀態(tài)。

圖3 在凈電荷為3.5｜e｜和4.0｜e｜時(shí)，不同電極面內(nèi)應(yīng)變下溶質(zhì)Zn原子和電極最外表面鋅原子所帶電荷隨電極-溶質(zhì)原子距離的變化Fig.3 Charge curves of solute zinc atoms and zinc atoms on the outermost surface of electrode as a function of distance from electrode to solute atoms under different electrode surface strains when the net charges are 3.5｜e｜and 4.0｜e｜

圖4 凈電荷為4｜e｜時(shí)的符號(hào)回歸訓(xùn)練模型Fig.4 Symbolic regression trainingmodel under net charge of 4｜e｜

為了驗(yàn)證SR預(yù)測(cè)的臨界組合（電荷，應(yīng)變）的準(zhǔn)確性，進(jìn)行了凈電荷為4｜e｜、電極面內(nèi)應(yīng)變?yōu)?6.75%條件下的計(jì)算。如圖5所示，實(shí)線為符號(hào)回歸模型預(yù)測(cè)值，方塊為采用JDFTx計(jì)算的電極表面應(yīng)變?yōu)?6.75%時(shí)的數(shù)值，可見(jiàn)計(jì)算值與符號(hào)回歸模型預(yù)測(cè)值基本吻合，說(shuō)明了符號(hào)回歸模型的預(yù)測(cè)精度較高。在2.46?＜z＜4.5?范圍內(nèi)，模型預(yù)測(cè)值和計(jì)算值都為0.61 eV左右，沒(méi)有出現(xiàn)能量勢(shì)壘，即溶質(zhì)原子在吸附位置處于臨界狀態(tài)。在z＞6.9?時(shí)，兩者均逐漸穩(wěn)定，說(shuō)明符號(hào)回歸模型可以預(yù)測(cè)溶質(zhì)原子進(jìn)入電解質(zhì)液后與電極表面微弱的相互作用。結(jié)合JDFTx計(jì)算和符號(hào)回歸的預(yù)測(cè)值，得到了兩組溶質(zhì)原子溶解／吸附臨界組合（電荷，應(yīng)變），分別為（3.5｜e｜，0%）和（4｜e｜，16.75%）。從圖中還能看出，當(dāng)電荷小于6.9?時(shí)，另外兩組以電極表面應(yīng)變?yōu)樽宰兞康那€的應(yīng)變值存在差異，而在z＞6.9?后逐漸與電極表面應(yīng)變16.75%時(shí)的曲線重合，并趨于穩(wěn)定，與上述解釋相符。

圖5 符號(hào)回歸算法預(yù)測(cè)曲線和JDFTx計(jì)算值Fig.5 Curves predicted by using symbol regression algorithm and values calculated by JDFTx

3 結(jié)束語(yǔ)

本文使用密度泛函／連續(xù)介質(zhì)溶液模型耦合計(jì)算研究了應(yīng)變-電荷耦合場(chǎng)下Zn電極的溶解機(jī)制。在凈電荷3.5｜e｜下分別在電極表面內(nèi)施加了5%（拉伸）和-5%（壓縮）應(yīng)變，鋅電極表面溶質(zhì)原子的近似溶解／吸附平衡狀態(tài)被打破，施加5%應(yīng)變抑制了溶質(zhì)Zn原子的溶解，而施加-5%應(yīng)變則促進(jìn)了溶質(zhì)Zn原子的溶解。此外，在系統(tǒng)凈電荷為4.0｜e｜以及施加不同拉伸應(yīng)變的條件下，拉伸應(yīng)變?cè)酱笠种谱饔迷矫黠@，直至出現(xiàn)能量勢(shì)壘阻止了溶質(zhì)Zn原子的溶解，且拉伸應(yīng)變的增加導(dǎo)致溶質(zhì)Zn原子進(jìn)入電解質(zhì)溶液所帶電荷也下降。結(jié)合JDFTx計(jì)算值和符號(hào)回歸模型，可以預(yù)測(cè)固定凈電荷下Zn電極表面溶質(zhì)原子溶解／吸附平衡狀態(tài)的臨界應(yīng)變。但是本文計(jì)算得出的拉伸應(yīng)變抑制溶解的結(jié)果與通常認(rèn)知不符，可能的原因是通常試驗(yàn)采用的是宏觀多晶材料，而本文計(jì)算的是一種完美單晶金屬材料。