随着社会的发展，人民生活水平的提高，能源的开发和利用是我们目前生活的重中之重^［1］。天然气水合物作为一种绿色能源，引起了人们的高度重视。天然气水合物多存在于海洋沉积物、永久冻土区和大陆冰盖中。水合物中天然气的回收就是利用CO₂置换水合物中的CH₄^［2］。这是一个贮存CO₂和回收有价值的天然气^［3-4］的有效方法。

天然气水合物^［5］包含三种类型，即sI型、sII型和sH型。气体分子如CH₄、H₂S、SO₂、CO₂、Xe、C₂H₆等可以形成sI型水合物。sI型水合物由两个十二面体（5¹²）和六个十四面体(5¹²6²）组成。直径D在0.58 ～ 0.72 nm范围内的气体分子可以形成sII型水合物。sII型水合物由八个十六面体（5¹²6⁴）和十六个五边形十二面体（5¹²）组成。当气体分子非常大时，形成sH型水合物^［6］。sH型水合物由三个正十二面体（5¹²）、两个不正规十二面体（4³5⁶6³）和一个二十面体（5¹²6⁸）组成。sH型水合物在有辅助气体存在时，可以形成稳定的sH型水合物，如CH₄，占小笼。

Yezdimer等^［7］根据由分子动力学（molecular dynamics，MD）模拟得到的置换过程中吉布斯自由能的变化分析了sI型水合物中CO₂置换CH₄的热力学可行性。由于置换过程的吉布斯自由能的变化为负值，因此置换过程是自发过程。基于MD模拟确定模型参数^［8］，根据相场理论确定了在液态CO₂，如水下天然气水合物矿藏条件下，CH₄水合物转化为CO₂水合物的转化率。采用分子动力学模拟方法^［9-11］研究了sI型水合物中CO₂置换CH₄水合物的分子机理。

sH型水合物由有机大分子构成，例如2−甲基丁烷、2，3−二甲基丁烷、新己烷、甲基环己烷、金刚烷、甲基叔丁基醚（TBME），并以CH₄作为辅助气体。采用分子动力学模拟方法^［12］研究了sH型水合物的热力学性质和结构稳定性，然而，仿真结果受到MD模拟势能模型的强烈影响^［13］。Arshad采用从头算法，HF/6−31G～*基组研究了CH₄分子在十二面体(H₂O)₂₀中的稳定能^［14］，基于B3LYP/6−311++G*基组研究了CO₂在十二面体(H₂O)₂₀/十四面体(H₂O)₂₄和十六面体(H₂O)₂₈中的稳定能^［15］。Wen等^［16］基于B3LYP/6−31G基组研究了CO₂在sI型水合物的十二面体和十四面体中的稳定性。文献［17］就笼型水合物的性质和工业能源开发前景进行了讨论，强调了要特别注意笼型水合物的分子模拟的发展和进步，从而加强对笼型水合物在分子水平上的理解。

本文基于密度泛函理论，首先对比研究CH₄和CO₂在正十二面体(H₂O)₂₀小笼和不正规十二面体(H₂O)₂₀中笼中的稳定性，以及CO₂和环庚烷在正二十面体(H₂O)₃₆大笼中的稳定性。最后，计算包含不同CO₂分子个数的大笼的稳定能，从而明确CO₂置换sH型水合物中CH₄的可能性。

1 模拟方法

本文基于密度泛函理论以及自洽场（SCF）方法，对sH型水合物的各种笼子进行了几何优化，进而得到了在同一密度泛函GGA−PBE（Perdew−Burke−Ernzerhof）/TZP基组下笼子的总能量。

首先进行分子动力学模拟，模拟体系由2 × 2 × 2个水合物晶胞组成，在ME（materials explorer）平台上，等温等容(NVT)系综下，在3个正坐标方向上均使用周期性边界条件和最近映像法。各分子的初始位置根据衍射实验的结果确定，水分子和客体分子的取向随机分配。采用5阶Gear算法拟合运动方程，时间步长为0.8 fs，模拟步数为500 000步，其中前10 000步用于使体系达到平衡。各分子间Lennard−Jones作用项采用截断球近似，长程静电力应用Ewald加和算法。选定温度通过MD模拟得到sH型水合物初始结构。在ADF（amsterdam density functional）平台上，对sH型水合物进行几何结构优化计算。根据计算所得的单点能量得到在同一密度泛函GGA−PBE/TZP基组下笼子的总能量。

笼子的稳定能 $\Delta {E_{\rm GH - m}}$ 等于笼子的总能量减去水分子和客体分子的总能量，即

$ \quad\quad\Delta {E_{\rm GH \,-\, m}} = {E_{\rm GH}} - {\sum\limits_i {({E_{{{\rm H}_2}{\rm O}}})} _i} - {E_{\rm g}} $

(1)

式中：E_GH为笼子的总能量；E $_{{{\text{H}}_{\text{2}}}{\text{O}}} $ 为每个水分子的总能量；E_g为客体分子的总能量；i为水分子的个数。

范德华相互作用能ΔE_g−c等于笼子的总能量减去空笼和客体分子的总能量，即

$ \quad\quad\Delta {E_{\rm g \,-\, c}} = {E_{\rm GH}} - {E_{\rm c}} - {E_{\rm g}} $

(2)

式中，E_c为含有氢键水域的空笼的总能量。

2 结果与讨论

在几何优化过程中，只要确定了5个条件（能量变化、最大约束梯度、约束梯度有效值、最大步长、步长有效值），那么笼子和单体的最优结构就确定了。在最优结构处，笼子的键能最低。每个周期的计算时间和原子数目（客体分子和笼子）成正比。C₇H₁₄·(H₂O)₃₆的键能收敛于第550个周期，在优化结构处，键能最低。整个优化过程分为三个不同的阶段，包括有冻核方法的第一阶段，无冻核方法的第二阶段和收敛精度从4提高到5的第三阶段。可以发现，包含一个客体分子的笼子的优化结构会因初始结构和仿真过程的不同而不同。

为了验证计算结果，定量讨论了(H₂O)₃₆和C₇H₁₄·(H₂O)₃₆中水分子的结构特性。一个大笼包含36个水分子。图1为水分子的H—O—H键角、O—H $ \cdots $ O键角分布。从图1(a)中可以看出，键角分布在104.8°左右，接近实验值104.2°^［18］。文献［14］、［16］中H—O—H键角平均值均为106°。图1(b)中O—H $ \cdots $ O键角主要集中在177.6°左右。当氢键相当有效时，O—H $ \cdots $ O预期键角为180°，而文献［14］、［16］中分别为171°、176°。O—O键长分布如图2所示，峰值是0.281 nm，实验值在0.272 5～0.279 1 nm之间^［19］，而文献［14］、［16］中分别为0.291 4、0.272 nm。五边形内O—O—O键角约为107.97°（一般为108°），六边形内O—O—O键角约为119.57°（一般为120°）。采用DFT/SCF计算的sH型水合物大笼的结构特性如表1所示，显然计算结果与实验结果吻合良好。因此，所选择的计算方法正确。

2.1 结构特性

sH型水合物中，包含CH₄和CO₂的小笼、中笼的优化结构分别如图3所示。采用不同的初始结构进行结构优化，只有在达到笼子的最低键能时，才能得到相似的优化结构。由图中可知，小笼的空笼直径（平行于O=C=O方向）为0.774 nm，而包含CH₄的小笼直径增大了0.006 nm，包含CO₂的小笼直径增大了0.008 nm，这定量证明了CH₄和CO₂比小笼大些。中笼的空笼直径为0.802 nm，高度为（0.651 nm, 0.654 nm），而包含CH₄的中笼高度（0.004 nm，0.001 nm）和直径（0.006 nm）均有很小的收缩，包含CO₂的中笼，由于存在O=C=O键，高度（0.006 nm，0.020 nm）增大，直径（0.014 nm）收缩。由于包含CO₂的中笼比包含CH₄的中笼的扭曲变形大，因此CO₂分子比中笼大，而CH₄分子更适合中笼。由图3中可知，CH₄分子比CO₂分子更适合小笼和中笼。

同样，采用DFT/SCF计算包含不同CO₂分子数的大笼的几何优化过程，并与包含环庚烷分子的大笼进行对比。结果表明，大笼可以放置多达3个CO₂分子，优化结构的垂直视图和顶部视图如图4所示。由图中可知，空笼和包含1个CO₂分子的大笼较大，扭曲变形较大；而包含环庚烷分子的大笼较小，对称性较好；包含2、3个CO₂分子的大笼则居于两者之间。一般情况下，笼子越小结构越稳定。因此，空笼和包含1个CO₂分子的大笼最不稳定，包含2、3个CO₂分子的大笼较稳定，包含环庚烷分子的大笼最稳定。显然，环庚烷分子更适合大笼，并且可以形成稳定的sH型水合物^［20］。空笼和包含小分子的笼子并不意味着是较小的笼子。只有大小合适的分子才会使笼子尺寸收缩变小。对于包含2、3个CO₂分子的大笼中哪一个更稳定将参照后面的稳定能进行分析。

从文献［5］和本文计算得到的sH型水合物的实验结构特性如表2所示。由表中可以看出，计算结果与实验结果有较好的一致性。

2.2 能量特性

表3、4分别为包含CH₄和CO₂分子的小笼、中笼和sH型水合物大笼的总能量。由于采用不同的DFT/SCF方法计算获得的绝对能量不同，因此只有相对能量才有意义。由表3中可知，空的小笼比空的中笼更稳定。同时，包含CH₄分子的小笼、中笼均比包含CO₂分子的小笼、中笼稳定。结果表明：CH₄分子比CO₂分子更适合小笼、中笼；CO₂不能置换sH型水合物小笼和中笼中的CH₄。由此可知，CH₄在形成sH型水合物时易作为辅助气体占据小笼、中笼，而文献［12, 18］中已证实CH₄可作为形成sH型水合物的辅助气体。而CO₂在形成sH型水合物时不宜作为辅助气体，这一结论与文献［20］中的结果一致。

由表4中可知，包含C₇H₁₄分子的大笼比包含1、2和3个CO₂分子的大笼更稳定。这与文献［20］中的实验结果，即环庚烷可以在较低压力下形成稳定sH型水合物一致。因此，CO₂不能置换sH型水合物大笼中的环庚烷。另一方面，CO₂分子和空笼之间的范德华力是引力，与1、3个CO₂分子相比，2个CO₂分子更适合大笼，这说明CO₂可以通过每个大笼包含2个CO₂分子的形式形成sH型水合物。高压下CH₄能形成sH型水合物^［22］，这已被实验证实。

3 结　论

基于密度泛函理论，在GGA−PBE/TZP基组下对sH型水合物进行了第一性原理计算。通过分子动力学方法模拟得到sH型水合物的初始结构，进而对sH型水合物进行了几何结构优化。根据计算得到的笼子的总能量，对比研究了包含不同客体分子的笼子的稳定能和范德华相互作用能。计算结果表明：常压下CH₄可作为形成sH型水合物的辅助气体，而CO₂则不能；在sH型水合物中，CO₂既不能置换大笼中的环庚烷，也不能置换小笼、中笼里的CH₄。但是CO₂可以通过每个大笼包含2个CO₂分子的形式形成sH型水合物。