陳 浩 顏克鳳* 李小森 陳朝陽 張 郁 徐純剛

(1.中國科學(xué)院 廣州能源研究所， 廣州 510640； 2.中國科學(xué)院 廣州能源研究所 中國科學(xué)院天然氣水合物重點實驗室，廣州 510640； 3.中國科學(xué)院 廣州天然氣水合物研究中心， 廣州 510640；4.廣東省新能源和可再生能源研究開發(fā)與應(yīng)用重點實驗室， 廣州 510640；5.中國科學(xué)院大學(xué)， 北京 100049)

引 言

天然氣水合物是以CH4、C2H6、CO2等氣體分子為客體，以H2O分子為主體，在低溫高壓下形成的一種籠型晶體化合物，其中主體分子以氫鍵相互連接形成不同的籠型孔穴，客體分子占據(jù)孔穴。天然氣水合物的晶體結(jié)構(gòu)主要包括I型(SI)、II型(SII)和H型(SH)。I型水合物晶胞是由2個512小孔穴和6個51262大孔穴總共46個H2O分子形成的體心立方結(jié)構(gòu)；II型水合物晶胞是面心立方結(jié)構(gòu)，包含136個H2O分子，由8個51264大孔穴和16個512小孔穴組成；H型水合物晶胞是簡單的立方結(jié)構(gòu)，包含34個H2O分子，晶胞中有3種不同的孔穴，即3個512孔穴、2個435663孔穴和1個51268孔穴。在自然界中天然氣水合物主要以SI晶體結(jié)構(gòu)存在，客體分子主要為CH4。由于具有儲量大、熱值高、無污染等優(yōu)點，天然氣水合物作為21世紀(jì)的一種新能源，其開采研究對于我國后續(xù)能源的可持續(xù)發(fā)展具有重要意義。目前，天然氣水合物的開采方法主要有熱激法、降壓法和化學(xué)試劑法，各個方法均有一定的局限，因此開采方法的優(yōu)化，以及新型開采方法的研發(fā)是天然氣水合物開采技術(shù)研究的重點。CO2置換開采水合物方法是一種新型的開采方法，在從水合物中置換出天然氣的同時注入CO2從而生成CO2水合物。該方法在開采出水合物資源的同時存儲CO2，可有效地減少碳排放，降低溫室效應(yīng)[1]。同時，CO2水合物的形成可有效防止水合物分解引發(fā)的地質(zhì)問題。因此，CO2置換開采天然氣水合物技術(shù)被認(rèn)為是可用于水合物開采和減少碳排放的具有巨大應(yīng)用潛力的一種手段。

實驗和理論研究表明，注入CO2置換天然氣水合物中的CH4是可行的[2]。從熱力學(xué)上看，CO2水合物體系的平衡條件要比CH4水合物體系的平衡條件溫和[3]。從動力學(xué)上看，當(dāng)溫度為250 K時，CO2置換CH4水合物過程的吉布斯自由能ΔG是負(fù)值，這意味著置換反應(yīng)是自發(fā)進行的[4]。Yuan等[5]利用三維反應(yīng)器對CO2置換CH4水合物過程的有利條件進行研究，發(fā)現(xiàn)具有游離氣的水合物儲層適合置換開采。對CO2置換CH4水合物過程的微觀分析和組分測量表明，置換主要發(fā)生在水合物相中[6]，CO2水合物和CH4水合物同時存在于水合物中，置換過程不存在CH4水合物的解離[7]。有部分研究表明CH4水合物中只有51264籠中的CH4被置換，因此置換效率較低[8]。為提高置換效率，探究過程中的置換機理成為CO2置換CH4水合物技術(shù)研究的重點。目前，宏觀實驗和微觀測試的研究能夠獲得置換過程的熱力學(xué)和動力學(xué)特性，但還不足以解釋置換機理。計算機模擬是在微觀尺度上揭示置換機理的一種有效方法。近年來，研究者們將分子動力學(xué)(molecular dynamics，MD)模擬用于研究CO2置換CH4水合物過程的置換機理，模擬結(jié)果表明，置換過程可以分為3個步驟[9]:首先打破水合物籠子，然后CH4分子跑出籠子，同時CO2分子進入空籠并占據(jù)籠子進行置換。Liu等[10]利用MD模擬發(fā)現(xiàn)置換出的CH4分子在形成氣泡的過程中并沒有出現(xiàn)漆膜，而且CO2分子的存在能夠促使CH4分子形成氣泡，防止CH4重新生成水合物，從分子角度闡述了CO2的置換機理。MD模擬能從分子尺度解析置換的動態(tài)過程，量子力學(xué)(quantum mechanics，QM)模擬則能從量子尺度，通過從頭算(Ab initio calculation)以及密度泛函理論(density functional theory，DFT)的方法對水合物主客體分子間的電子效應(yīng)進行分析，從量子層面解釋置換機理[11]。Geng等[12]利用MD 模擬和QM模擬對CH4、CO2以及CH4- CO2水合物的穩(wěn)定性進行研究，通過分子尺度的徑向分布函數(shù)、均方位移、擴散系數(shù)以及量子尺度的結(jié)合能的對比分析，發(fā)現(xiàn)CH4- CO2混合水合物的穩(wěn)定性最高，此結(jié)果與Ota等[13]的實驗結(jié)果相一致。Liu等[4]利用從頭算方法對CO2取代CH4占據(jù)水合物籠子的結(jié)合能變化及熱力學(xué)特性進行研究，發(fā)現(xiàn)CO2占據(jù)51262籠在熱力學(xué)上的有利性。

由上述研究可知，水合物中主、客體分子對水合物的穩(wěn)定性具有重要的影響。為更深入地闡述CO2置換CH4水合物過程的本質(zhì)，本文采用DFT方法對水合物中主、客體雙分子聚體之間的相互作用進行模擬，通過對稱性匹配微擾理論(symmetry adapted perturbation theory，SAPT)進行能量分析，解析主、客體相互作用中靜電、色散、誘導(dǎo)和交換4種相互作用各自的貢獻。同時計算了波函數(shù)信息，分析了約化密度梯度函數(shù) (reduced density gradient，RDG)、獨立梯度模型(independent gradient model，IGM)、靜電勢等參數(shù)，探討主、客體對水合物籠子穩(wěn)定性影響的作用機理，為深層次地闡述置換機理提供理論依據(jù)。

1 模型及計算方法

本文利用DFT方法對水合物中主、客體分子聚體，即H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體進行研究，針對它們的幾何結(jié)構(gòu)、相互作用能量、相互作用方式進行計算模擬。模擬通過幾何優(yōu)化確定雙分子聚體的最優(yōu)結(jié)構(gòu)，然后計算單點能，采用式(1)計算雙分子聚體的相互作用能。

ΔEa-b=Etotal-Ea-Eb

(1)

式中，ΔEa-b表示雙分子聚體中分子之間的相互作用能，a、b分別表示聚體中的分子；Etotal表示聚體的總能量；Ea和Eb分別表示雙分子聚體中各個分子的能量。模擬結(jié)構(gòu)優(yōu)化及能量計算都在Gaussian16程序C.01版本[14]上進行。

本文首次采用高階SAPT[21]對水合物中主、客體雙分子聚體進行能量分解，評估雙分子聚體相互作用中靜電、色散、誘導(dǎo)和交換4種相互作用對主、客體相互作用的貢獻，以及影響置換過程的最主要的相互作用方式。整個SAPT計算都在PSI4程序[22]上進行。同時在幾何優(yōu)化后的結(jié)構(gòu)的基礎(chǔ)上，使用Multiwfn程序[23]計算波函數(shù)信息。通過RDG分析(也稱為non-covalent interaction，NCI分析)[24]進行弱相互作用的可視化研究。采用IGM[25]對弱相互作用的強度進行分析。通過靜電勢分析獲得靜電作用在CO2置換CH4水合物過程中的影響方式。

2 結(jié)果與討論

2.1 主、客體雙分子聚體結(jié)構(gòu)優(yōu)化和相互作用能

水合物的晶體結(jié)構(gòu)主要由主體分子和客體分子組成的籠子構(gòu)建而成，其中主體分子與客體分子間的相互作用決定籠子的穩(wěn)定性。在CO2置換CH4水合物過程中，水合物籠子的主體H2O分子、客體CH4分子和置換的CO2分子之間的相互作用影響著整個置換過程。本文采用不同的泛函在6- 311+G(2d,2p)基組下進行結(jié)構(gòu)優(yōu)化，然后采用CCSD(T)/aug- CC- pVTZ水平計算相互作用能。模擬優(yōu)化得到的H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體的幾何構(gòu)型如圖1所示，雙分子聚體相互作用能計算結(jié)果如表1所示。從表1中可以看出，經(jīng)B3LYP- D3泛函優(yōu)化后計算得到的雙分子聚體的能量相比于M062X- D3、WB97XD、B97D這3個泛函的計算結(jié)果，更接近于基準(zhǔn)MP2函數(shù)的計算結(jié)果，表明B3LYP- D3泛函更適合用于描述H2O- H2O、H2O- CH4、H2O- CO2之間的相互作用。由于基準(zhǔn)MP2函數(shù)計算時間較長[26]，因此在B3LYP- D3/6- 311+G(2d,2p)水平下進行CO2置換CH4水合物過程的幾何結(jié)構(gòu)優(yōu)化。

圖1 雙分子聚體的優(yōu)化結(jié)構(gòu)Fig.1 Optimized structures of dimers

表1 H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體相互作用能Table 1 Interaction energies of H2O- H2O, H2O- CH4,H2O- CO2 dimers

在B3LYP- D3/6- 311+G(2d,2p)水平下進行的3種雙分子聚體結(jié)構(gòu)優(yōu)化的基礎(chǔ)上，分別利用B3LYP、B3LYP- D3、M062X、M062X- D3、WB97XD、B97D、MP2等方法在6- 311+G(2d,2p)基組下進行雙分子聚體相互作用能的計算，并與CCSD(T)/aug- CC- pVTZ水平下的計算結(jié)果進行對比，如表2所示。從表2可以看出，相比于其他計算水平，M062X/6- 311+G(2d,2p)水平下計算獲得的分子間相互作用能與CCSD(T)/aug- CC- pVTZ水平下的計算結(jié)果最接近，表明其計算精度較其他方法的要高。CCSD(T)/aug- CC- pVTZ水平主要用于計算小分子單體，對于大分子體系計算較為困難[27]，

表2 不同計算水平下的雙分子聚體相互作用能Table 2 Interaction energies of dimers for different calculation levels

因此，M062X/6- 311+G(2d,2p)水平是CO2置換CH4水合物過程的QM研究的較優(yōu)計算結(jié)合能參數(shù)。

從圖1、表1和表2可以看出，在H2O- H2O聚體中，不同H2O分子中的H原子與O原子由氫鍵連接形成相互作用，相互作用能為-21.05 kJ/mol。文獻[28]對氫鍵進行了系統(tǒng)分類，指出弱氫鍵相互作用能小于-10.46 kJ/mol時，相互作用由靜電和色散作用主導(dǎo)；“weak-to-medium”氫鍵相互作用能在-10.46～-58.60 kJ/mol時，相互作用由靜電作用主導(dǎo)；中氫鍵相互作用能在-46.04～-62.79 kJ/mol之間時，相互作用由靜電作用主導(dǎo)；強氫鍵相互作用能大于-62.79 kJ/mol時，相互作用由靜電和感應(yīng)作用主導(dǎo)。因此H2O分子中H原子與相鄰的H2O分子中O原子間的相互作用介于弱氫鍵相互作用與中氫鍵相互作用之間。在H2O- CH4聚體中，H2O分子中的H與CH4中的C形成微弱的吸引作用，相互作用能為-4.06 kJ/mol。在H2O- CO2聚體中，H2O中的O原子與CO2中的C原子相互作用，相互作用能為-12.18 kJ/mol。可以看出H2O分子之間的相互作用明顯高于H2O- CH4和H2O- CO2，因此在水合物中H2O分子與客體分子之間相互作用偏弱，更傾向于與其他H2O分子相互作用形成籠子。并且與H2O- CH4相比，H2O- CO2具有更高的相互作用能，當(dāng)形成籠型水合物時，則更利于籠子的聚集，得到的CO2水合物較CH4水合物更加穩(wěn)定[4]。

2.2 主、客體雙分子聚體相互作用中各項相互作用的貢獻

在化學(xué)體系中，弱相互作用為分子間或分子內(nèi)的非化學(xué)鍵作用，包括范德華作用、pi- pi堆積作用、氫鍵、二氫鍵、鹵鍵等。SAPT利用片段間的微擾將弱相互作用分解為靜電能(Eelst)、交換能(Eexch)、色散能(Edis)和誘導(dǎo)能(Eind)4部分，從而從能量角度更深入地理解弱相互作用的本質(zhì)。本文通過SAPT對CO2置換CH4水合物過程中主、客體雙分子聚體的能量進行分析，解析主、客體相互作用中靜電、色散、誘導(dǎo)和交換4種相互作用的貢獻。H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體在不同SAPT計算水平下的相互作用能列于表3。將其與CCSD/aug- CC- pVTZ水平下的相互作用能比較發(fā)現(xiàn)，在SAPT下，計算級別的提升并不一定會伴隨計算精度的提升，并且加入δMP2項來考慮誘導(dǎo)和色散間的高階耦合作用，加入(CCD)項以優(yōu)化色散項并不會使得相互作用能的計算值更接近CCSD/aug- CC- pVTZ水平下的計算值。通過對比分析發(fā)現(xiàn)SAPT2+/aug- CC- pVDZ、SAPT2+/jun- CC- pVTZ和SAPT2/aug- CC- pVTZ 水平下計算的相互作用能更接近CCSD/aug- CC- pVTZ水平下的計算結(jié)果，具有較高的計算精度，因此它們在分析H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體相互作用時更具優(yōu)勢。在Parker等[21]有關(guān)于效率計算的研究中指出SAPT2+/jun- CC- pVTZ水平耗時是SAPT2+/aug- CC- pVDZ和SAPT2/aug- CC- pVTZ水平的5倍以上。當(dāng)進行小分子聚體計算時，耗時差異較小，當(dāng)進行水合物籠體的相關(guān)研究時，耗時差異則較大。而在實際的水合物籠體計算中，SAPT2/aug- CC- pVTZ水平計算過程的臨時文件占用幾T的硬盤儲存空間，存在因硬盤空間不足而使得計算終止的問題。因此，考慮到計算精度與計算時長，水合物在SAPT分析中較適合采用SAPT2+/aug- CC- pVDZ水平計算，同時SAPT2+/aug- CC- pVDZ水平是文獻[21]提出的SAPT計算中的“銀”標(biāo)準(zhǔn)，能夠進行精準(zhǔn)的能量分解研究。

表3 不同SAPT計算水平下雙分子聚體相互作用能Table 3 Interaction energies of dimers for different SAPT calculation levels

在SAPT方法中，弱相互作用能可以被分解為具有不同物理意義的組分，即ΔE=Eelst+Eexch+Eind+Edis。其中，Eelst主要反映主客體之間的經(jīng)典靜電相互作用；Eexch表示在主客體中由電子的費米子行為引起的單體波函數(shù)以及反對稱的重疊交換排斥作用；Eind描繪了誘導(dǎo)作用，包括極化以及兩種單體之間的電荷轉(zhuǎn)移；Edis是由主客體的電子之間的庫侖相關(guān)性引起的色散貢獻。分子占據(jù)籠子時，主、客體相互作用能的差異明顯，為了闡明水合物中主、客體相互作用的電子性質(zhì)并探究影響主客體間穩(wěn)定性的最主要的相互作用形式，我們對H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體在SAPT2+/aug- CC- pVDZ計算水平下進行SAPT分析，結(jié)果見表4。從表4中可以發(fā)現(xiàn)相互作用能的負(fù)值是由于靜電、誘導(dǎo)、色散3項的總和大于交換排斥項的貢獻產(chǎn)生的，因此總的相互作用主要為吸引作用。為了考察各吸引作用對總吸引作用的貢獻，更直觀深入地理解主客體相互作用的性質(zhì)，將表4中的數(shù)據(jù)轉(zhuǎn)化為每個吸引項的貢獻分?jǐn)?shù)C，如式(2)所示。

表4 SAPT2+/aug- CC- pVDZ水平下計算的雙分子聚體相互作用能物理分量Table 4 The physical components of the interaction energy of the dimers calculated at the SAPT2+/aug- CC- pVDZ level

(2)

其中EX表示Eelst、Eind或Edis。SAPT衍生的各相互作用對總吸引作用的貢獻分?jǐn)?shù)如圖2所示。從圖中可以看出，在H2O- H2O聚體中靜電作用占比為64.25%，即吸引作用主要是由靜電作用引起的，而誘導(dǎo)和色散作用的影響較小，占比分別為19.53%和16.21%。在H2O- CH4聚體中靜電作用占比為36.00%，低于色散作用的45.20%，即吸引作用主要是由色散作用引起的，而誘導(dǎo)作用占比僅為18.80%，對總的吸引作用影響較小。在H2O- CO2聚體中色散作用占比為29.62%，相比于H2O- CH4聚體中高的色散貢獻下降15.58%，靜電作用占比為59.38%，較H2O- CH4聚體增加23.38%。由此發(fā)現(xiàn)在雙分子聚體中，隨著相互作用能的增加，靜電作用有上升趨勢，色散作用則表現(xiàn)為下降，即靜電作用占比與相互作用能大小成正比，色散作用占比與相互作用能大小成反比。因此，可以推斷客體分子CO2和CH4與主體H2O分子作用時靜電和色散相互作用占比的差異導(dǎo)致了兩種客體分子形成水合物時結(jié)構(gòu)和能量的差別。

圖2 SAPT衍生的各相互作用對總吸引作用的貢獻分?jǐn)?shù)Fig.2 The contribution percentage of each interaction derived from SAPT to the total attraction

綜上分析可知，在3種聚體的能量分解中H2O- H2O聚體的吸引作用主要由靜電作用貢獻，色散和誘導(dǎo)作用的影響較小。在H2O- CH4、H2O- CO2聚體的對比中發(fā)現(xiàn)，與H2O相互作用的分子由CH4轉(zhuǎn)變?yōu)镃O2時色散作用的影響減弱，靜電作用的影響加強，表明了靜電作用在CO2置換CH4中的重要地位。由此推斷，增加靜電作用是提升置換效率的有效方法，例如加入與H2O分子有較高靜電作用的試劑分子，可以提高置換效率。

2.3 主、客體分子的靜電勢分析

靜電勢V(r)(ESP)是由原子核和電子在周圍空間每一點r處產(chǎn)生的分子間或分子內(nèi)真實的物理特性，其值能夠通過衍射實驗以及計算的方法確定，計算公式為

(3)

式中，ZA表示位于RA處原子核A所帶的電荷數(shù)，ρ(r)表示分子的電子密度，r′表示距離r單位電子電荷增量距離。在任意點r上的V(r)反映了分子中每個核和每個電子在該點的貢獻[29]。在任何給定的區(qū)域，V(r)是正還是負(fù)取決于原子核的正電荷或電子的負(fù)電荷在特定位置的主導(dǎo)性。V(r)是正的區(qū)域為親核性，V(r)是負(fù)的區(qū)域為親電性。計算過程中當(dāng)電子密度為0.001 e/Bohr3(0.006 784 e/?3)時，在分子表面繪制靜電勢等值面，如圖3所示。V(r)的最正值和最負(fù)值分別用Vmax和Vmin表示。圖3(a)為主體H2O分子的靜電勢等值面圖，從圖中可以看出H2O分子中H原子靜電勢的表面電位為正，且Vmax值為189.12 kJ/mol和189.20 kJ/mol，O原子靜電勢的表面電位為負(fù),Vmin值為-140.44 kJ/mol。由此可以得出H2O- H2O聚體中H原子與O原子之間的靜電勢差達到329.64 kJ/mol。圖3(b)為客體CH4分子的靜電勢等值面圖，可以看出H原子具有正的表面電荷，Vmax值為35.37 kJ/mol，Vmin點位于CH4的兩個H原子之間。由圖3(a)和圖3(b)得出，H2O分子的Vmax與CH4分子的Vmin之差以及CH4分子的Vmax與H2O分子的Vmin之差分別為197.45 kJ/mol和175.81 kJ/mol，因此H2O分子的Vmax位點更利于同CH4分子的Vmin位點相互作用，即H2O分子中的H與CH4分子中的C形成相互作用，結(jié)果與圖1一致。圖3(c)為客體CO2分子的靜電勢等值面圖，由圖發(fā)現(xiàn)碳原子表面點位為正，Vmax值為118.08 kJ/mol，氧原子表面靜電勢為負(fù)，Vmin值為-48.14 kJ/mol。由圖3(a)和圖3(c)得出，在H2O- CO2聚體中H2O分子的Vmax與CO2分子的Vmin之差以及CO2分子的Vmax與H2O分子的Vmin之差分別為237.26 kJ/mol和258.52 kJ/mol，因此H2O分子的O原子更利于與CO2分子的C原子形成相互作用，這與圖1中的幾何結(jié)構(gòu)優(yōu)化結(jié)果相同。在H2O- CO2聚體與H2O- CH4聚體的對比中，H2O- CO2聚體中兩種分子的靜電勢差值(258.52 kJ/mol)高于H2O- CH4聚體的相應(yīng)值(197.45 kJ/mol)，但低于H2O- H2O聚體的329.64 kJ/mol，因此H2O與CO2相互作用的強度高于H2O與CH4分子，而低于H2O分子與H2O分子間的相互作用，這與相互作用能計算結(jié)果(表1)的趨勢相匹配。圖4為 H2O- H2O、H2O- CH4、H2O- CO2雙分子聚體的靜電勢表面穿透圖，穿透位置與Vmax、Vmin的最大差值一致，對于H2O- H2O、H2O- CH4、H2O- CO2這3種聚體，其分子間靜電勢表面相互穿透距離分別為1.16 ?、0.45 ?、0.98 ?，相互作用強度與相互穿透距離的大小具有一致性，進一步驗證了靜電作用在H2O- H2O、H2O- CH4、H2O- CO2相互作用中的關(guān)鍵影響。

圖3 靜電勢等值面圖Fig.3 Isosurface diagrams of electrostatic potential

圖4 靜電勢表面穿透圖Fig.4 Surface penetration diagrams of electrostatic potential

2.4 RDG及IGM分析相互作用可視化

RDG分析基于電子密度及其降低的梯度，可定性鑒定和表征各種強度的弱相互作用。IGM理論則基于電子密度，能夠識別和量化由于相互作用而導(dǎo)致的電子密度凈梯度衰減，定義分子間相互作用的區(qū)域，從而直接表征分子間的相互作用。本文通過計算波函數(shù)信息，對主、客體分子進行RDG和IGM分析，定向研究主、客體分子對水合物籠子穩(wěn)定性的影響機理。RDG分析是通過生成雙分子聚體真實空間相應(yīng)區(qū)域的梯度等值面來定位它們，通過這些區(qū)域的sign(λ2)ρ值來區(qū)分相互作用的種類，其中sign(λ2)函數(shù)是電子密度Hessian矩陣的第二大本征值λ2的符號，用以反映相互作用類型；弱相互作用臨界點的ρ的數(shù)值和鍵的強度具有正相關(guān)性。sign(λ2)ρ值從0.05→0→-0.05逐漸減小，相互作用變化為強互斥作用(位阻作用)→范德華作用→強吸引作用(氫鍵、鹵鍵等)，在RDG圖中表現(xiàn)為顏色從紅色→綠色→藍色的變化，如圖5所示。當(dāng)sign(λ2)ρ>0時表現(xiàn)為互斥作用，sign(λ2)ρ<0時表現(xiàn)為吸引作用。從圖5中可以看出H2O- H2O聚體中H2O分子之間形成了較強的吸引作用，即氫鍵相互作用，而H2O- CH4、H2O- CO2聚體中分子間呈現(xiàn)出范德華相互作用。H2O- CO2聚體的范德華作用形成的sin(λ2)ρ峰值比H2O- CH4聚體的峰值更負(fù)，因此H2O分子與CO2分子之間的范德華作用高于H2O分子與CH4分子之間的范德華作用，這與客體分子和水合物籠子之間的相互作用趨勢一致[30]。

圖5 雙分子聚體的RDG圖Fig.5 RDG diagrams of the dimers

圖6為雙分子聚體間的IGM圖，圖中的sign(λ2)ρ=0附近有左右兩個峰，當(dāng)sign (λ2)ρ<0時，峰高表示吸引作用，當(dāng)sign (λ2)ρ>0時，峰高則表示排斥作用。從圖6可以看出，吸引作用的峰尖比排斥作用的峰尖高，表明雙分子聚體主要表現(xiàn)為吸引作用。由圖6(a)看出吸引作用的峰尖主要為藍色，即H2O- H2O聚體是以氫鍵相互作用為主。由圖6(b)和圖6(c)看出吸引作用的峰尖主要為綠色，即H2O- CH4、H2O- CO2聚體中分子間呈現(xiàn)范德華力相互作用，該結(jié)果與RDG分析一致。在分子間相互作用區(qū)域，IGM方法定義的δg函數(shù)值與相互作用強度具有正相關(guān)性。對比3種雙分子二聚體相互作用的δg值(圖6)發(fā)現(xiàn)δg, H2O- H2O>δg,H2O- CO2>δg, H2O- CH4。這與表2中相互作用能的大小具有一致性，進一步驗證了3種雙分子二聚體相互作用間的差異。

圖6 雙分子聚體的IGM圖Fig.6 IGM diagrams of the dimers

由上可知，RDG分析可表征主、客體分子之間的相互作用，IGM分析可獲得主、客體分子之間的相互作用區(qū)域，判斷相互作用大小，從而實現(xiàn)主、客體相互作用的定性及定量研究。本文研究從電子密度分析角度證實了H2O- H2O聚體是以氫鍵相互作用為主，H2O- CH4、H2O- CO2聚體中分子間是以范德華力相互作用為主。由此，主、客體雙分子聚體分子間的相互作用是影響水合物籠子穩(wěn)定性的關(guān)鍵因素。

3 結(jié)論

(1)對比不同DFT計算方法的精度，得到B3LYP- D3/6- 311+G(2d,2p)水平是CO2置換CH4水合物過程的QM研究的較優(yōu)結(jié)構(gòu)優(yōu)化參數(shù)，M062X/6- 311+G(2d,2p)水平是CO2置換CH4水合物過程的QM研究的較優(yōu)計算結(jié)合能參數(shù)。

(2)結(jié)構(gòu)優(yōu)化與靜電勢分析結(jié)果表明，在CO2置換CH4水合物過程中，H2O分子的H原子與CH4分子的C原子形成相互作用，H2O分子的O原子與CO2分子的C原子形成相互作用；主體H2O分子與客體分子相互作用的強度低于H2O- H2O相互作用，因此主體H2O分子傾向于相互作用形成籠子；H2O- CO2相互作用的強度高于H2O- CH4相互作用，因此，當(dāng)形成水合物籠子時，CO2水合物較CH4水合物更加穩(wěn)定。

(3)由靜電、色散、誘導(dǎo)和交換4種相互作用的貢獻分析得出，CO2置換CH4水合物過程中主、客體分子間的作用主要由靜電作用貢獻，色散和誘導(dǎo)作用的影響較小。置換過程中，客體分子由CH4轉(zhuǎn)變?yōu)镃O2時色散作用的影響減弱，靜電作用的影響加強。

(4)RDG分析和IGM分析從電子密度角度證實了CO2置換CH4水合物過程中主體H2O 分子之間以氫鍵相互作用為主，主、客體分子之間以范德華力相互作用為主。

由此，SPAT計算方法能較好地分析CO2置換CH4水合物過程中主、客體分子間的作用，研究結(jié)果表明靜電作用是置換的關(guān)鍵因素，提高與H2O的靜電作用是提升置換效率的有效方法。下一步我們將采用此QM方法用于天然氣水合物分解機理的微觀理論研究，與MD技術(shù)和微觀實驗結(jié)合分析CO2置換CH4水合物的動態(tài)過程，為天然氣水合物技術(shù)的發(fā)展提供理論指導(dǎo)。