朱后禹, 徐 靜, 匙玉華, 趙聯明, 郭文躍
(中國石油大學(華東) 理學院, 山東 青島 266580)
計算材料學(computational materials science)是材料科學與計算機交叉學科。該學科利用現代高速計算機技術,結合物理學、化學、生物學、材料學等相關學科的基礎知識,可以對材料從微觀到宏觀多尺度的結構、性能等進行理論計算[1]。
計算機模擬技術可以模擬超高溫、超高壓等極端環境下的材料使用性能,模擬材料在使用條件下的性能演變規律和機理,進而實現材料使用性能的改善。因此,在現代材料學領域中,計算機實驗已經成為與傳統實驗室的實驗具有同等地位的研究手段[2]。 對于原子層次上的問題的求解,一般基于密度泛函理論(density functional theory,DFT)[3]的第一性原理(first principles)計算方法[4]。在現代石油工業中,加氫脫硫過程是一個重要環節,為此本文嘗試對加氫脫硫反應從理論上進行實驗設計。借助于Materials Studio程序包(Accelrys公司)中的DMol3模塊進行計算。Materials Studio軟件是一個跨越各種尺度的研究平臺,不僅可以解析電子結構,還可以預測材料的宏觀結構性能。實驗采用了基于密度泛函理論的第一性原理計算方法,依托超算中心集群資源,通過模擬計算給出了噻吩在MoP表面的吸附幾何結構、能量,以及脫硫反應的機理。基于學生所掌握的基本理論知識,選取難度適中的實驗內容、適合計算材料學實驗教學。
由于石油燃料的輸運過程的增加以及環境標準的提高,全球的石油領域工作人員一直在努力尋找一種可以更為清潔的開采使用能源的方法,脫硫反應便是其中的一個重點。石油煉化過程中,一個重要的脫硫方法便是加氫脫硫[5],它是石油化合物精煉的一個尤為重要的催化反應[6]。在脂肪族硫化物的加氫脫硫反應中,含鈷和鎳的二硫化鉬(MoS2)催化劑的引入對其有很好的效果[7],但是對噻吩類物質效果不佳[8]。過渡金屬磷化物在加氫脫硫反應中的催化作用在近幾十年已經得到了廣泛的研究[9]。研究發現,磷化鉬(MoP)比傳統的硫化鉬在加氫脫硫反應中起到了更好的催化性能[10-11]。本實驗應用Materials Studio軟件,通過基于周期性密度泛函理論對噻吩在磷化鉬不同表面上的吸附量化計算模擬,加深了學生對于吸附能、吸附位點、過渡態等概念的理解。
模型建立是模擬計算中最基礎的一步,所構建的模型是否合理決定了實驗是否可以正確進行。模型建立分為如下步驟:
(1) 首先在可視界面建模時從軟件自帶結構庫中導出MoP晶胞;
(2) 利用軟件建模工具進行切面、建立超晶胞、建立真空層等操作,綜合考慮模擬真實性與計算精度、效率,超晶胞大小取3層厚度;
(3) 增加噻吩分子在MoP表面,此時需要考慮噻吩在MoP表面可能的各種位置,從而確定其表面最穩定的位置。
本次實驗的所建的初始模型結構見圖1。

圖1 MoP表面模型
本次實驗的模擬計算包括模型優化和計算。優化計算得到的體相MoP的晶格參數:a=b=0.326 6 nm和c=0.319 8 nm,與實驗結果(a=b=0.322 3 nm,c=0.319 1 nm)符合較好。后續計算也采用此參數以保證MoP能量保持在最低值。交換相關能以廣義梯度近似(GGA)以及Perdew和Wang提出的(PW91)泛函計算[12-13]采用密度函數半核贗勢(DSPP)方法描述金屬原子的核內電子,用雙極化基組加極化函數(DNP)處理價電子波函數,在DSPP優化結構基礎上,采用全電子進行單點能量計算來校正磷化鉬表面和吸附物質的電子性質[14-16]。所有的計算都用自旋極化進行。能量、位移和梯度的收斂標準分別為5.442×10-4eV,5×10-4nm和1.088 eV/nm。本次實驗選取噻吩(T)在(001)表面上穩定吸附的T-bridge-fcc和(010)表面上穩定吸附的T-hollow-top構型分別開展計算。反應過渡態計算采用CI-NEB方法。
在MoP(001)上,噻吩有平躺吸附和垂直式兩種類型的分子吸附,平躺吸附包括T-bridge-hollow(T-bridge-fcc和T-bridge-hcp)、T-hollow(T-fcc和T-hcp)和T-cross-bridge構型,圖2示出了典型的吸附構型。為了研究在本結構中表面硫原子所起到的作用,計算了兩種類型的表面硫原子,一種是S原子共享表面一個Mo原子,記為SS;另一種是C原子共享表面一個Mo原子,記為SC。在這種T-bridge-fcc構型中,噻吩分子中心在Mo2—Mo4橋位上方,噻吩S原子位于fcc位的上方,且分別與C2和C5共用Mo2和Mo4,另外2個C原子共用Mo3原子。在MoP表面,與金屬表面比較,分子碳環平面傾斜了11°(圖2(a))。在S-修飾的表面上,表面原子S和噻吩的共吸附結構見圖2(b)和圖2(c)。與在清潔表面上不同,形成一個新的S—Mo鍵,在硫修飾的表面上,噻吩中靠近表面S原子的鉬原子的鍵長被拉伸了6.6%和3.1%,C2—S鍵縮短約為0.9%。同時觀察到由表面硫原子引起的吸附物的移動。C—C和C—S鍵略微改變(C—C鍵0.2%,C—S鍵2.1%),噻吩在磷化鉬(001)面上的吸附能,SS結構中減小8.5%和SC結構中減少21.2%。

圖2 在清潔和S-修飾的MoP(001)面上噻吩(T)的T-bridge-fcc吸附結構(C為黑色、S黃色、H白色、Mo淺灰色、P紫色)
在MoP(010)表面上,類似于MoP(001)的吸附情況見圖3。噻吩環上的碳原子和噻吩下面的4個Mo原子的標記如圖3(a)所示。為了研究表面S原子的影響,計算了兩種類型的表面S原子,一種是S原子共享表面一個Mo原子,記為SS,另一種是C原子共享表面一個Mo原子,記為SC。在T-hollow-top構型(見圖3(c))中,S原子位于Mo4原子的top位,C3位于Mo1和Mo2間的橋位之外,其余C原子位于鄰近的Mo的top位。在清潔的磷化鉬(010)上,噻吩的吸附形成了一個新的S—Mo鍵(0.248 2 nm)和4個新的C—Mo鍵(約0.232 0 nm),引起了C—S(約0.011 nm),C2—C3和C4—C5(約0.09 5 nm)鍵的較大程度的拉伸和C3—C4鍵(0.004 9 nm)輕微的伸展。此種吸附的吸附能為-2.55 eV。
在S-修飾的表面上,由于表面S原子和噻吩S和C原子共享Mo原子,結果T-hollow-top有兩種穩定的共吸附構型,即T-hollow-top_SS和T-hollow-top_SC(見圖3(b)和(c))。與清潔表面的情況相比,S修飾的表面因為形成了新的S—Mo鍵,涉及的噻吩的S/C—Mo鍵被拉伸0.001 4/0.002 4 nm,而噻吩的骨架鍵(C—C和C—S鍵)在0.000 1~0.002 3 nm的范圍內變化。吸附能分別為-2.27和-2.28 eV,反映出表面硫原子減弱了噻吩在磷化鉬(010)面上的吸附。

圖3 噻吩(T)在清潔和S-修飾的MoP(010)表面上的穩定吸附構型(C原子黑色、S原子黃色和H原子白色)
在MoP(001)表面,T-bridge-fcc脫硫的相關勢能面(PES)以及所涉及的中間體的結構見圖4。第一步,在TS1-A中,C2—S鍵斷裂,反應能壘僅為0.18 eV,斷裂距離為0.216 7 nm, C2—C3鍵經歷了一個大的移動稍微偏向橋位,隨后使C2原子移動到鄰近的fcc位點,C3原子移到原來在IS中的C3和C4共享的表面Mo原子的頂部,形成了一個相當穩定的末態(FS),thiolate-A過程放熱2.05 eV, C—S鍵斷裂,組建一個有Mo原子的六元環,并且C—S鍵中插入了1個金屬原子;第二步, TS2-A過程,C—S鍵的活化是由于C5—S鍵的伸縮振動,系統穩定1.31 eV,在FS((C4H4+S)-A)中,C5原子占據fcc位點,并且解離的S原子朝著相反的方向移向鄰近的hcp位點,從而導致大的C5—S分離(0.482 6 nm)。在TS中,噻吩的S原子和C4—C5基位于IS和FS之間的中間位置,即噻吩的S原子位于橋位和C5原子位于表面Mo原子的頂位,且C5—S距離為0.262 0 nm。該過程活化能為0.44 eV。

圖4 MoP(001)表面上噻吩bridge-fcc構型的直接脫硫路徑(T-bridge-fcc和吸附于無限遠處2個H原子的總能量作為能量參考(eV))
圖5給出了T-bridge-fcc構型的初始加氫路徑以及相關結構。定義噻吩2號位(圖5上面的路徑)和噻吩3號位(圖5下面的路徑)的加氫路徑分別為M和N。在過程M中,初始狀態是鄰近C2、C3的hcp位點的氫和bridge-fcc的位點的噻吩原子共同吸附。在過程TS-M中,氫原子向C2運動,H—Mo距離0.175 5 nm,C2—C3鍵偏向C2尾端的H原子, C2—S鍵斷裂,斷裂距離0.222 1 nm,C2—Mo鍵略拉長,C2與進攻的氫原子的間距是0.226 1 nm。過渡態后,C2—C3—C4角得到較大擴展,C2—C3鍵位于鄰近fcc位點的頂點,且攻擊的H與C2成鍵,最終C2—S分離(分離距離0.312 8 nm)達到穩定末態。此過程放熱1.35 eV,能壘是0.63 eV。在IS-N構型里,攻擊的H在C3原子的fcc位點上,與IS-M相比有點不穩定。在過渡態TS-N里,H到達與C3相鄰的鉬原子的頂點,H—Mo鍵距離0.182 1 nm。C3原子已經調整到便于攻擊H原子的位置。這種結構預測1.48 eV的加氫能壘,在TS-N之后,3-MHT-bridge-fcc形成,該過程吸熱0.82 eV。很明顯,清潔表面上,噻吩2號位的加氫更具有優勢。

圖5 MoP(001)表面上噻吩T-bridge-fcc構型的初始加氫路徑(相應參數標注同于圖4)
在清潔MoP(010)表面,T-hollow-top構型的脫硫反應發生了2個連續的C—S鍵斷裂(見圖6)。C2—S鍵斷裂產生thiolate(C4H4S)-A,能壘為0.74 eV,反應能-0.42 eV。TS1-A中,S位移導致C2—S鍵斷裂,C5—S鍵略微收縮(0.006 6 nm)。然后,穩定的中間體thiolate-A成為第一個C—S鍵斷裂的終態(FS),S和C2保持在近鄰bridge位,C2—S距離又被拉長。下一個C—S鍵斷裂通過TS2產生穩定的FS((C4H4+S)-A),該系統比thiolate-A系統能量低0.32 eV。在TS2-A中,C5—S鍵已斷裂,斷裂距離0.248 3 nm,預測活化能壘0.55 eV。FS中,C4H4已各自移動到四重hollow位點,末端兩個碳原子位于top和bridge處,解離的S位于bridge位。

圖6 MoP(010)表面上噻吩T-hollow-top構型的直接脫硫路徑(T-hollow-top加上無限遠處吸附的2個H原子總能量作為能量參考(eV))
對于T-hollow-top在清潔的表面上(見圖7),在位置2和3處的加氫能壘為0.97 eV和0.98 eV,反應能各為0.46 eV和0.56 eV。

圖7 MoP(010)表面上噻吩T-hollow-top構型的初始加氫路徑
在本實驗的實施過程中,前期布置學生調研相關文獻,調動學生主觀能動性; 課上以學生為主體討論實驗步驟,分析實驗結果;課后上交完整的實驗報告。由于本實驗依托于先進的超算中心計算平臺,使得學生在掌握理論和實驗知識的同時,也掌握了Windows和Linux雙系統操作能力。該實驗鞏固了學生的專業基礎知識,促進了理論與實踐結合,提高了學生的科研素養。教學實踐表明,學生更愿意接受這種參與性強、內容具有前沿性的研究型實驗。
[1] 張躍,谷景華,尚家香,等. 計算材料學基礎[M]. 北京: 北京航空航天大學出版社,2007.
[2] 堅增運,劉翠霞,呂志剛,計算材料學[M]. 北京: 化學工業出版社,2012.
[3] Hohenberg P,Kohon W. Inhomogeneous electron gas[J]. Phys Rev,1964(136):B864-B871.
[4] Segall M D,Lindan P J D,Probert M J,et al. First-principles simulation: ideas illustrations and theCASTEP code[J]. J Phys: Cond Matt,2002(14):2717-2743.
[5] Rodriguez J A, Kim J Y, Hanson J C, et al. Physical and Chemical Properties of MoP, Ni2P, and MoNiP Hydrodesulfurization Catalysts: Time‐Resolved X‐Ray Diffraction, Density Functional, and Hydrodesulfurization Activity Studies[J]. Cheminform, 2003, 34(39):1-8.
[6] Wang Huamin, Roel Prins. HDS of benzothiophene and dihydrobenzothiophene over sulfided Mo/γ-Al2O3[J] Appl Catal A: Gen, 2008, 350(2):191-196.
[7] Krebs E, Silvi B, Daudin A, et al. A DFT study of the origin of the HDS/HydO selectivity on Co(Ni)MoS active phases[J]. J Catal, 2008, 260(2):276-287.
[8] Yang R T, Hernández-Maldonado A J, Yang F H. Desulfurization of transportation fuels with zeolites under ambient conditions[J]. Science, 2003, 34(40):79-81.
[9] Oyama S T. Novel catalysts for advanced hydroprocessing: transition metal phosphides[J]. J Catal, 2003, 216(1):343-352.
[10] Phillips D C, Sawhill S J, Self R, et al. Synthesis, Characterization, and Hydrodesulfurization Properties of Silica-Supported Molybdenum Phosphide Catalysts[J]. J Catal, 2002, 207(2):266-273.
[11] Wu Z, Sun F, Wu W, et al. On the surface sites Of MOP/SiO2 catalyst under sulfiding conditions: IR spectroscopy and catalytic reactivity studies[J]. J Catal, 2004, 222(1):41-52.
[12] Perdew J P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas[M]// Density-functional approximation for the correlation energy of the inhomogeneous electron gas. New York: NRC Research Press,1986:8822-8824.
[13] Perdew J P, Wang Y. Accurate and simple analytic representation of the electron-gas correlation energy[J]. Phys Rev B: Condens.Matter, 1992, 45(23):13244.
[14] Feng Z, Liang C, Wu W, et al. Carbon Monoxide Adsorption on Molybdenum Phosphides: Fourier Transform Infrared Spectroscopic and Density Functional Theory Studies[J]. J Phys Chem B, 2003, 107(49):13698-13702.
[15] Jun Ren, Chunfang Huo, Xiaodong Wen, et al. Thiophene Adsorption and Activation on MoP(001), γ-Mo2N(100), and Ni2P(001):Density Functional Theory Studies[J]. J Phys Chem B, 2006,110(45):22563.
[16] Li G, Zhu H, Zhao L, et al. Theoretical Survey of Thiophene Hydrodesulfurization Mechanism on Clean and Single-Sulfur-Atom-Modified MoP(001)[J]. J Phys Chem.C, 2016, 120(40):23009-23023.