分子動力學(xué)模擬中弛豫對甲烷水合物晶胞構(gòu)型的影響

廖清云，史博會，宋尚飛，王君傲，姚海元，段 旭，宮 敬，聞 峰

（1. 中國石油大學(xué)（北京） 天然氣水合物國家重點實驗室，油氣管道輸送安全國家工程實驗室，石油工程教育部重點實驗室，城市油氣輸配技術(shù)北京市重點實驗室，北京 102249；2. 中海油研究總院有限責(zé)任公司，北京 100028；3. 國家石油天然氣管網(wǎng)有限公司，北京 100007）

隨著深海油氣資源開發(fā)的快速發(fā)展，氣體水合物的開發(fā)利用及其在流動保障中的防治研究逐漸深入，其中，氣體水合物生長及分解過程的微觀機理成為了研究熱點?？紤]到常規(guī)實驗手段在時間尺度和空間尺度上的限制，國內(nèi)外已有多位學(xué)者利用分子動力學(xué)模擬進行了微觀尺度的研究，包括水合物在純水體系中的生長機理[1]、添加劑或雜質(zhì)（如無機鹽/水合物動力學(xué)抑制劑/天然果膠/蠟分子等）的存在對水合物生長機理的影響[2-5]、水合物在逐步升溫過程中的分解機理[6]、水合物在無機鹽中的分解機理[7]及不同溫度和壓力條件下抑制劑存在時水合物的分解機理[8]等。

在水合物生長及分解過程的分子模擬研究中，均需在初始模型中添加根據(jù)實驗數(shù)據(jù)[9]搭建的水合物晶胞結(jié)構(gòu)作為水合物生長種子或水合物分解前體結(jié)構(gòu)。因水合物晶胞結(jié)構(gòu)中氫原子位置具有不確定性，在分子模擬中需根據(jù)Bernal-Fowler規(guī)則、體系構(gòu)型的勢能最小原則或氫鍵水分子朝向成比例原則等，進一步確定氫原子坐標。在構(gòu)建模擬模型后，需開展能量最小化以降低體系能量并尋找體系勢能極小值以消除體系內(nèi)不合理接觸，隨后通過弛豫過程，即一段有限時長的分子動力學(xué)模擬，以進一步調(diào)整能量最小化后的模型，使系統(tǒng)從初始配置逐漸演變達到平衡狀態(tài)，才能為分子模擬計算提供更為合理的初始構(gòu)型，該過程對獲得合理的分子模擬結(jié)果至關(guān)重要。

弛豫過程以分析體系的熱力學(xué)參數(shù)（如溫度、能量和密度等）及結(jié)構(gòu)參數(shù)（如序參數(shù)、均方位移和徑向分布函數(shù)等）為主[10]。然而，不同研究者選擇的弛豫系綜及弛豫時間各有不同：Tung等[1]以具有固定溫度的正則系綜（NVT）（T= 200 K）弛豫20 ps平衡氣-液-固體系；Li等[3]以NVT系綜（T= 276 K）弛豫400 ps平衡液-固體系；Xu等[4]以NVT系綜（T= 260 K）弛豫30 ps，及具有固定溫度及壓力的系綜（NPT）（T= 260 K，p= 15 MPa）弛豫30 ps平衡液-固體系；Choudhary等[6]以NPT系綜（p= 100 MPa）弛豫5 ns平衡水合物晶胞結(jié)構(gòu)；Xu等[7]以NVT系綜在模擬溫度弛豫50 ps平衡液-固體系；Kondori等[8]以NVT系綜（T= 273 K）弛豫50 ps平衡水合物晶胞結(jié)構(gòu)；張慶東等[11]首先以NVT系綜（T= 263.15 K、T= 273.15 K、T= 283.15 K）固定氧原子位置弛豫50 ps，隨后解除氧原子約束在相同條件下弛豫1 ns平衡水合物晶胞結(jié)構(gòu)?？傊?，NVT系綜在弛豫中被廣泛應(yīng)用，且弛豫時間通常設(shè)置在皮秒量級，但關(guān)于如何確定合理的弛豫方案的研究較少。

本文以甲烷（CH4）水合物晶胞籠型結(jié)構(gòu)的演變規(guī)律為分析參數(shù)，探討不同弛豫方案對CH4水合物晶胞構(gòu)型的影響，包括：弛豫系綜、弛豫時間、弛豫壓力與溫度等，目的是確定更為合理的弛豫方案以開展CH4水合物的分子動力學(xué)模擬，為水合物分子動力學(xué)模擬研究提供理論基礎(chǔ)。

1 模擬模型和方法

1.1 CH4水合物晶胞模型

圖1 3 × 3 × 3 sI型CH4水合物晶胞示意Fig. 1 Schematic diagram of 3 × 3 × 3 sI methane hydrate unit cell

1.2 模擬設(shè)定

基于相關(guān)文獻研究[1,3,5,16]，對于分子間相互作用，確定選用TIP4P-Ew力場[17]描述H2O分子，OPLS-AA力場[18]描述CH4分子。各原子的電荷及L-J勢函數(shù)參數(shù)如表1所示，不同原子間的L-J勢函數(shù)參數(shù)值通過Lorentz-Berthelot混合規(guī)則計算。將H2O分子視為剛體，采用SHAKE算法[19]約束H2O分子的H-O鍵長與H-O-H鍵角。范德華力采用L-J勢函數(shù)計算，截斷半徑為1 nm。短程庫倫作用截斷半徑為0.9 nm，長程庫倫相互作用通過pppm算法計算[20]，計算精度為10-4。模擬時間步長為1 fs，每0.5 ps輸出一次熱力學(xué)參數(shù)與原子位置信息。利用Nose-Hoover恒溫器及恒壓器實現(xiàn)體系的溫度壓力控制，模擬體系在x、y和z方向上均采用了周期性邊界條件。

表1 L-J勢函數(shù)參數(shù)值與原子電荷值Table 1 L-J potential function parameter value and atomic charge value

1.3 CH4水合物晶胞構(gòu)型穩(wěn)定性判定方法

本文以一種改進的籠型結(jié)構(gòu)識別方法判定CH4水合物晶胞構(gòu)型穩(wěn)定性，其具體流程如下：

（1）篩選成籠氧原子。當(dāng)相鄰氧原子間的間距小于0.35 nm，且相鄰的兩個氧原子及與其相連的一個氫原子間形成的角度小于30°，則認為這兩個氧原子之間形成氫鍵[21]，當(dāng)一個氧原子上形成的氫鍵數(shù)大于等于3時則認為其為成籠氧原子。

（2）基于氧原子篩選結(jié)果，進行五元環(huán)的識別。將成籠氧原子視為無向圖的頂點，根據(jù)各氧原子之間的連接關(guān)系來判斷是否形成五元環(huán)。

（3）512籠結(jié)構(gòu)如圖2所示，綠色實線為512半籠的中心五元環(huán)，紅色實線為半籠中與中心五元環(huán)相連的五元環(huán)，藍色實線為完整的512半籠結(jié)構(gòu)?？紤]到其為對稱結(jié)構(gòu)，可采取半籠拼接的方法進行識別[22]。根據(jù)五元環(huán)的識別結(jié)果，判斷各五元環(huán)之間的連接關(guān)系，進而確定以各五元環(huán)為中心的512半籠，再將得到的512半籠進行拼接，最終得到512全籠。

圖2 512籠結(jié)構(gòu)示意Fig. 2 Schematic diagram of 512 cage structure

水合物晶體構(gòu)型內(nèi)氫原子的位置會影響氫鍵判定中的角度計算結(jié)果，并影響所識別的環(huán)籠結(jié)構(gòu)數(shù)目。因此，以該改進的籠型識別方法分析弛豫過程中體系籠型結(jié)構(gòu)數(shù)目，若所識別的籠型結(jié)構(gòu)數(shù)目越接近理論值且波動越小，則說明該弛豫方案的設(shè)定越優(yōu)。對于3 × 3 × 3 sI型CH4水合物晶胞結(jié)構(gòu)而言，按理論分析應(yīng)有35個512籠，其連接關(guān)系如圖3所示，其中紅色小球為氧原子，藍色實線為氫鍵。

本文的成本系數(shù)(Ck)是指該方案的生命周期成本(工程造價)與所有方案成本總和之比。各方案的價值系數(shù)Vk = Fk / Ck，計算結(jié)果見表3。可見在4種防腐方案中，D方案的價值系數(shù)最高，即采用功能性玻璃鋼復(fù)合材料作為火電廠煙囪防腐工程材料是最優(yōu)解決方案。

圖3 3 × 3 × 3 sⅠ型CH4水合物晶胞中512籠的正交示意(a)及透視示意 (b)Fig. 3 Orthogonal schematic diagram (a) and perspective schematic diagram (b) of 512 cages in 3 × 3 × 3 sⅠmethane hydrate unit cell

2 結(jié)果與討論

2.1 弛豫系綜及弛豫時間對CH4水合物晶胞構(gòu)型的影響

為研究弛豫系綜及弛豫時間對CH4水合物晶胞構(gòu)型的影響，設(shè)定弛豫方案1-1及弛豫方案1-2。弛豫時間均為500 ps；弛豫方案1-1設(shè)定為NVT系綜，弛豫方案1-2設(shè)定為NPT系綜；以水合物-液態(tài)水-甲烷氣三相平衡溫壓條件為基礎(chǔ)，10 MPa時的水合物-液態(tài)水-甲烷氣三相平衡溫度約為281 K[1]，故暫設(shè)定弛豫溫度為250 K，弛豫壓力為10 MPa。為避免瞬時結(jié)構(gòu)波動對結(jié)果分析造成影響，將512籠識別數(shù)量在3 ps的時間窗口內(nèi)進行平均后再進行分析。

方案1-1（NVT，T= 250 K，下同）弛豫500 ps的512籠識別結(jié)果如圖4(a)所示。在500 ps弛豫過程中，512籠數(shù)量變化可分為3個階段，如表2所示。在58 ps內(nèi)快速達到近平衡狀態(tài)，在58 ps后隨著模擬時間的延長逐漸趨于更穩(wěn)定的平衡狀態(tài)。此外，體系溫度在± 5%范圍內(nèi)波動，如圖4(b)所示，表明體系溫度已達到平衡狀態(tài)。

表2 不同弛豫方案下512籠數(shù)量變化階段Table 2 Change stages of 512 cage number under different relaxation schemes

圖4 方案1-1弛豫500ps的512籠識別結(jié)果(a)及溫度變化(b)Fig. 4 512 cage recognition results (a) and temperature changes (b) of scheme 1-1 relaxation 500 ps

方案1-2（NPT，T= 250 K，p= 10 MPa，下同）弛豫500 ps的512籠識別結(jié)果如圖5(a)所示。在500 ps弛豫過程中，512籠數(shù)量在弛豫初期波動較大，其變化可分為4個階段，如表2所示。弛豫在20 ps內(nèi)快速達到近平衡狀態(tài)，在20 ps后隨著模擬時間的延長逐漸趨于512籠數(shù)量更高的近平衡狀態(tài)，在364 ps后達到更穩(wěn)定的平衡狀態(tài)。此外，體系溫度在± 5%范圍內(nèi)波動，如圖5(b)所示，表明體系溫度已達到平衡狀態(tài)。

圖5 方案1-2弛豫500 ps的512籠識別結(jié)果(a)及溫度變化(b)Fig. 5 512 cage recognition results (a) and temperature changes (b) of scheme 1-2 relaxation 500 ps

由表2可知，500 ps弛豫結(jié)束后，兩種弛豫方案下體系中識別得到的512籠數(shù)量難以穩(wěn)定達到35個。這是由于H2O分子在模擬中不斷運動，其原子存在越過模擬盒邊界的情況，進而導(dǎo)致邊界處H2O分子所形成的氫鍵數(shù)未達到成籠條件而被忽略。方案1-1和1-2弛豫后的體系結(jié)構(gòu)如圖6所示，綠色實線為氫鍵，模擬盒邊緣的各方向應(yīng)各有3個512籠，因本環(huán)籠識別代碼中的氫鍵篩選過程暫未考慮周期性邊界條件，致使方案1-1未識別出的512籠位于模擬盒的左下邊界及右中邊界，方案1-2未識別出的512籠位于模擬盒的右下邊界，故512籠數(shù)量最終呈現(xiàn)出在接近理論值附近波動變化的現(xiàn)象。

圖6 方案1-1(a)及方案1-2(b) 弛豫后的體系結(jié)構(gòu)Fig. 6 512 cages of scheme 1-1 (a) and scheme 1-2 (b) after relaxation

對于弛豫系綜，對比方案1-1和方案1-2得到的512籠識別結(jié)果，如圖7所示。由圖7(a)可知，在方案1-2下，512籠數(shù)量能更快達到較高的數(shù)值，但其在364 ps內(nèi)的數(shù)值波動較大。經(jīng)對比發(fā)現(xiàn)，在較短的模擬時長內(nèi)，方案1-1得到的結(jié)果穩(wěn)定性更好。因此，可以認為方案1-2（NPT系綜）下弛豫對籠型結(jié)構(gòu)的形成有一定的促進作用，但在方案1-1（NVT系綜）下弛豫更利于結(jié)構(gòu)的穩(wěn)定。這一現(xiàn)象與Zhang等[23]的結(jié)論一致：在NPT系綜下水合物模擬的成核速率最快，而在NVT系綜下模擬得到的水合物結(jié)晶度更高。由圖7(b)可知，方案1-1和1-2弛豫過程中體系溫度變化相近，均在平衡值附近波動。

圖7 方案1-1與方案1-2弛豫的512籠識別結(jié)果(a)及溫度變化(b)對比Fig. 7 Comparison of 512 cage recognition results (a) and temperature change (b) between scheme 1-1 and 1-2

對于弛豫時間，由圖4(a)與圖5(a)可知，弛豫時間越長，越有利于獲得穩(wěn)定的CH4水合物晶胞構(gòu)型，但弛豫時間的增加會顯著延長計算時間。適宜的弛豫時間，不僅可以避免因弛豫時間較短不能使體系達到穩(wěn)定平衡的缺點，還可以彌補弛豫時間較長而導(dǎo)致計算時長增加的缺陷。對于波動較大的弛豫方案1-2，由于512籠數(shù)量在模擬初期變化明顯，需選取適當(dāng)較長的弛豫時間以達到體系的平衡，考慮到其在100 ps后波動幅度有所減小，且此時兩種弛豫方案內(nèi)識別得到的512籠數(shù)量均可達到32個，故選用100 ps作為本文后續(xù)分析研究的弛豫時間。

基于方案1-1及方案1-2，為進一步研究弛豫溫度和壓力條件及系綜聯(lián)用對CH4水合物晶胞構(gòu)型的影響，設(shè)定多種弛豫方案，如表3所示。設(shè)定弛豫時間為100 ps，分別設(shè)置不同溫度和壓力以分析弛豫溫壓條件對弛豫過程的影響，并研究NVT及NPT系綜聯(lián)合應(yīng)用對弛豫過程的影響。

2.2 弛豫溫度對CH4水合物晶胞構(gòu)型的影響

表3中弛豫方案2-1至2-3，是在弛豫方案1-1基礎(chǔ)上，調(diào)整弛豫時間為100 ps，改變弛豫溫度設(shè)定的3個方案。方案2-1至2-3的弛豫溫度分別設(shè)置為200 K、250 K和300 K。3種不同弛豫溫度下的籠識別結(jié)果對比如圖8(a)所示。對于最低弛豫溫度200 K，在相同弛豫時間下，體系內(nèi)識別得到的512籠數(shù)量最高，且其波動幅度最小。然而，對于最高弛豫溫度300 K而言，體系內(nèi)所識別的512籠數(shù)量波動幅度最大。弛豫溫度為250 K時，512籠識別數(shù)量結(jié)果處于兩者之間?？傮w而言，在弛豫時間為75~100 ps時，弛豫方案2-1至2-3所識別的512籠數(shù)量逐漸接近。

圖8 方案2-1至2-3弛豫的512籠識別結(jié)果(a)及溫度變化(b)Fig. 8 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 2-1 to 2-3

表3 弛豫方案Table 3 Relaxation schemes

因水合物的生成條件是低溫高壓，在低溫條件下，水合物生成的驅(qū)動力更強，故較低的溫度有利于水合物籠型結(jié)構(gòu)的穩(wěn)定，更利于籠型結(jié)構(gòu)識別分析。但隨著弛豫時間達到100 ps，不同弛豫溫度所達到的水合物晶胞構(gòu)型穩(wěn)定效果逐漸接近。此模擬結(jié)果與張慶東等[11]得到的結(jié)論一致，即弛豫溫度會影響CH4水合物的穩(wěn)定性，且弛豫溫度越低，其穩(wěn)定性越高；但對于已穩(wěn)定的水合物晶胞結(jié)構(gòu)，弛豫溫度的影響將不再明顯。因此，圖8(a)的模擬結(jié)果也再次佐證了本模擬體系內(nèi)CH4水合物晶胞構(gòu)型在弛豫100 ps時已基本達到穩(wěn)定。方案2-1至2-3弛豫過程中的溫度變化如圖8(b)所示，可知各體系的溫度基本在設(shè)定值附近波動。

越低的弛豫溫度越有利于CH4水合物晶胞構(gòu)型的穩(wěn)定，而過高的弛豫溫度會減弱CH4水合物晶胞構(gòu)型的穩(wěn)定性。對于TIP4P-Ew水分子模型，冰點約為245 K[24]，若弛豫體系中含有液態(tài)水，需設(shè)置高于冰點的模擬溫度以避免冰的產(chǎn)生。對于本文所選參數(shù)，三相平衡溫度在10 MPa時為281 K，弛豫溫度需低于所選力場模型的相平衡溫度以避免水合物結(jié)構(gòu)失穩(wěn)。對于本文研究的體系，250 K的弛豫溫度較合理，為此，表3中的弛豫方案3-1、3-2、4-1、4-2和4-3均選取250 K作為弛豫溫度，以開展后續(xù)研究。

2.3 弛豫壓力對CH4水合物晶胞構(gòu)型的影響

表3中弛豫方案3-1和3-2，是在弛豫方案1-2基礎(chǔ)上，調(diào)整弛豫時間為100 ps，改變弛豫壓力設(shè)定的兩個方案。方案3-1與3-2的弛豫壓力分別設(shè)置為10 MPa與30 MPa。方案3-1和方案3-2在100 ps內(nèi)的512籠識別結(jié)果對比如圖9(a)所示。數(shù)據(jù)表明，增大壓力后512籠數(shù)量的變化更加劇烈。在100 ps模擬末，30 MPa弛豫壓力下的512籠識別數(shù)量仍有顯著波動。這是因為分子動力學(xué)模擬通過調(diào)整模擬盒尺寸并隨之縮放各原子位置以實現(xiàn)壓力的控制，在高壓下體系的體積調(diào)整較大，故其結(jié)構(gòu)的波動變化也越劇烈，具體表現(xiàn)為在高壓下弛豫時，所識別的512籠數(shù)量穩(wěn)定性差。方案3-1和方案3-2下體系溫度的變化如圖9(b)所示，數(shù)據(jù)表明體系溫度基本在設(shè)定值附近波動。因此，對于NPT系綜下弛豫壓力的確定，一方面，要避免高壓所引起的體系結(jié)構(gòu)變化；另一方面，應(yīng)滿足后續(xù)正式模擬壓力需求。為此，針對本文研究體系，為保證結(jié)構(gòu)穩(wěn)定，選取10 MPa作為NPT系綜弛豫壓力。

圖9 方案3-1和3-2弛豫的512籠識別結(jié)果(a)及溫度變化(b)Fig. 9 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 3-1 to 3-2

2.4 雙系綜弛豫方案對CH4水合物晶胞構(gòu)型的影響

經(jīng)前續(xù)模擬結(jié)果分析，在適宜的弛豫時間和弛豫溫度下，NVT系綜下的弛豫能得到穩(wěn)定的結(jié)構(gòu)，但水合物結(jié)構(gòu)達到一定數(shù)量所需的時間更長；在適宜的弛豫溫度和壓力下，NPT系綜下可更快得到一定數(shù)量的水合物結(jié)構(gòu)，但水合物結(jié)構(gòu)的波動變化較大；兩個系綜下的弛豫均可平衡體系溫度。因此，考慮將兩個系綜相結(jié)合，借助NVT系綜易于得到穩(wěn)定結(jié)構(gòu)的優(yōu)勢，加以NPT系綜對籠型結(jié)構(gòu)形成的促進作用，得到表3中的弛豫方案4-1，總弛豫時長為200 ps，首先在NVT下弛豫100 ps以穩(wěn)定初始結(jié)構(gòu)，再在NPT下弛豫100 ps以進一步達到接近理論結(jié)構(gòu)的穩(wěn)定狀態(tài)。

方案4-1模擬過程中512籠識別數(shù)量的變化如圖10(a)所示。其亦可大致分為3個階段，如表4所示。弛豫在58 ps內(nèi)快速達到近平衡狀態(tài)，在58 ps后隨著模擬時間的延長逐漸趨于更穩(wěn)定的平衡狀態(tài)。方案4-1弛豫過程中體系的溫度變化如圖10(b)所示，弛豫過程體系溫度基本在設(shè)定值附近波動。因此，在弛豫方案4-1下，首先進行NVT弛豫再進行NPT弛豫可有效地平衡體系。

圖10 方案4-1弛豫的512籠識別結(jié)果(a)及溫度變化(b)Fig. 10 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 4-1

表4 方案4-1 512籠數(shù)量變化階段Table 4 Change stages of 512 cage number of relaxation scheme 4-1

將方案4-1、4-2及4-3進行對比，NVT弛豫后再進行NPT弛豫，512籠數(shù)量在NPT弛豫初期會出現(xiàn)一定波動，但在167 ps即達到33個512籠，且數(shù)值基本穩(wěn)定，這一階段快于如圖11(a)所示的單獨NVT弛豫，且其在100~200 ps弛豫時間內(nèi)波動范圍明顯小于如圖11(b)所示的單獨NPT弛豫。在200 ps的弛豫時間下，單獨進行NPT弛豫識別得到的512籠平均值更大，單獨進行NVT弛豫識別得到的512籠平均數(shù)量最少。以上結(jié)果進一步表明，NPT弛豫有利于籠型結(jié)構(gòu)的形成，但其結(jié)構(gòu)穩(wěn)定性較差，而NVT弛豫可優(yōu)化穩(wěn)定水合物結(jié)構(gòu)，兩者結(jié)合后有利于水合物弛豫平衡。

圖11 方案4-3 (a) 及方案4-2 (b) 與方案4-1弛豫的512籠識別結(jié)果對比Fig. 11 Comparison of 512 cage recognition results of relaxation between scheme 4-1 and scheme 4-3 (a) and scheme 4-2 (b)

此外，在NPT弛豫過程中，由于體系體積處于波動變化，故以分段密度變化考察體系平衡狀態(tài)，對所選方案的弛豫結(jié)果進行驗證。對于NPT弛豫過程，從體系基本穩(wěn)定（512籠達到31個）后開始每隔10 ps取一次密度的平均值，計算每段平均值的標準差，當(dāng)標準差小于0.001 g/cm3時[25]，認為體系達到平衡狀態(tài)。對于單獨NPT弛豫，從模擬20 ps后（表2）開始取密度平均值；對于NVT + NPT弛豫，從NVT弛豫結(jié)束后，即100 ps開始取密度平均值。密度計算結(jié)果如表5所示，單獨NPT弛豫在200 ps內(nèi)分段密度標準差仍大于0.001 g/cm3，到500 ps時分段密度達到平衡狀態(tài)。而首先進行NVT弛豫，再進行NPT弛豫，可在200 ps時達到密度的平衡狀態(tài)。因此，與水合物籠型結(jié)構(gòu)數(shù)量的變化類似，從分段密度的角度衡量，將NVT系綜及NPT系綜結(jié)合進行弛豫，水合物結(jié)構(gòu)更易達到平衡狀態(tài)。

2.5 雙系綜弛豫后水合物分解溫度的模擬分析

經(jīng)不同弛豫方案的模擬比較，擬推薦采用雙系綜弛豫方案4-1，即首先進行100 ps的NVT弛豫，再進行100 ps的NPT弛豫，故對雙系綜弛豫后處于平衡態(tài)的水合物構(gòu)型進行分解模擬以驗證弛豫方案的適用性。采用空隙誘導(dǎo)法進行水合物分解模擬，Choudhary等[6]的模擬表明，水合物分解溫度在空隙率為10%左右時處于平衡階段，此時的分解溫度可確定為該力場下的水合物分解溫度?；陔p系綜弛豫后的水合物構(gòu)型，隨機移除H2O分子及CH4分子共計115個，獲得空隙率約為8%的水合物結(jié)構(gòu)，如圖12所示，橙線為氫鍵，綠球為CH4分子中的碳原子，藍線為隨機移除H2O分子后破壞的氫鍵，藍球為隨機移除的CH4分子。

水合物分解模擬過程相關(guān)參數(shù)設(shè)置與弛豫過程一致，如1.2節(jié)所述。模擬壓力為10 MPa，模擬溫度由隨機移除分子后的平衡溫度185 K勻速升溫至350 K，升溫速率為0.165 K/ps，總模擬時長為10 ns。體系的溫度及勢能變化如圖13所示，勢能突變位置即為相變點，故該模型在10 MPa下的分解溫度約為292 K。對比10 MPa下水合物實際分解溫度286.2 K，模擬分解溫度高約5.8 K。這一偏差是由于水合物分解溫度模擬研究依賴于模擬體系大小，3 × 3 × 3 sI型CH4水合物晶胞模擬系統(tǒng)較小，而較小體系模擬所得到的水合物分解溫度較高[6]。此外，Choudhary等[6]使用3 × 3 × 3 sI型CH4水合物晶胞進行分解模擬時，所得分解溫度偏差范圍為5~35 K，故本模擬結(jié)果在偏差范圍內(nèi)且偏差較小。對于水合物分解過程，某一水合物籠坍塌后，水合物開始分解，如圖14(a)所示（6.5 ns時），在隨后的1 ns內(nèi)完全分解，由于分解后CH4分子過飽和，形成氣液兩相體系，如圖14(b)所示（7.3 ns時），其中橙線為氫鍵及H2O分子，綠球為CH4分子中的碳原子。以上模擬結(jié)果表明，雙系綜弛豫方案所得的穩(wěn)定水合物構(gòu)型可用于后續(xù)產(chǎn)生相的模擬分析，該弛豫方案可作為產(chǎn)生相模擬的基礎(chǔ)。

圖14 水合物分解模擬中6.5 ns (a)及7.3 ns (b)時體系構(gòu)型Fig. 14 System configuration at 6.5 ns (a) and 7.3 ns (b)during the hydrate decomposition simulation

3 結(jié)論

在進行產(chǎn)生相分子模擬前，需要對體系進行弛豫以合理化初始結(jié)構(gòu)。對于不同體系與模擬需求，弛豫過程的設(shè)置不盡相同。

（1）需選取適宜的弛豫時間。一定時間弛豫后，構(gòu)型變化并不明顯，過長弛豫時間會顯著增加計算時長。

（2）在NVT系綜下，弛豫有利于提高體系的結(jié)晶度，更低的弛豫溫度可更有效優(yōu)化并穩(wěn)定構(gòu)型，但溫度的設(shè)置需綜合考慮所選力場模型的相平衡條件及冰點。在體系基本穩(wěn)定后（合理的弛豫時間下），溫度影響不再顯著。故應(yīng)結(jié)合后續(xù)模擬需求考慮，選擇適中的弛豫溫度。

（3）NPT系綜下的弛豫有利于水合物結(jié)構(gòu)的形成，越高的弛豫壓力可在越短弛豫時間內(nèi)得到可識別的籠型結(jié)構(gòu)，但其穩(wěn)定性越差。弛豫壓力應(yīng)根據(jù)后續(xù)模擬壓力需求，結(jié)合相平衡條件及選定的弛豫溫度而確定。

（4）不同弛豫系綜所產(chǎn)生的弛豫效果不同。在穩(wěn)定溫度的基礎(chǔ)上，NVT弛豫更利于結(jié)構(gòu)的穩(wěn)定，而NPT系綜對籠型結(jié)構(gòu)的形成具有促進作用。在體系弛豫時可優(yōu)先考慮NVT弛豫以平衡溫度與穩(wěn)定結(jié)構(gòu)，再結(jié)合NPT弛豫以促進水合物結(jié)構(gòu)達到更接近理論數(shù)量的穩(wěn)定狀態(tài)。

因此，對于本文3 × 3 × 3 sI型CH4水合物晶胞，250 K的弛豫溫度及10 MPa的弛豫壓力可以得到較好的弛豫效果。在此溫度和壓力條件下，首先進行100 ps NVT弛豫，再進行100 ps NPT弛豫，即共200 ps的弛豫處理，所得到的平衡結(jié)構(gòu)可基本擬合水合物分解溫度及分解過程，具備開展產(chǎn)生相分子模擬的適用性。

符號說明

p―弛豫壓力，MPa；T―弛豫溫度，K；q―原子電荷值，e；σ―L-J勢函數(shù)參數(shù)，nm；ε―L-J勢函數(shù)參數(shù)，kJ/mol。