結果與討論

3-1 建立 Force Field

我們以Tinker version 7.1 中 mm3 參數檔作為基礎進行參數的調整，以嘗試錯 誤方式進行參數的調整，在調整過程中是以DFT 作為目標值。首先利用 geometric combining rules 將 MEA 分子分別拆成 CH4 dimer, H2O dimer, NH3 dimer，加上 CO2

dimer，並加入 Coarse-grained 的概念，忽略 H 原子 van der Waals 作用力，目的是 當反應的系統太大時，能達到降低計算量與計算成本。

我們將作用力分成分子間和分子內作用力進行參數調整。分子內作用力的部 分是利用DFT (Density Functional Theory)進行三個小分子的幾何優化後，得到最佳 構型的振動頻率，再以Tinker 進行參數調整後計算，目標是重現分子的振動頻率。

藉由這種嘗試錯誤的方式調整分子內作用力的參數值。分子間作用力包括靜電作 用力和van der Waals 作用力。靜電力的部分是利用 DFT 計算出分子 dimer 間氫鍵 作用力後，將所得作用力作為目標值，進行Tinker 中 dipole 參數調整。van der Waals 作用力的部分則是調整Lennard Jones potential 中 ε 和 rm兩參數。CO2的部分也是 利用相同的方式進行參數的調整，Figure 3.為 Tinker 與 DFT 對比版本。

CH4…CH4

H2O…H2O

NH3…NH3

CO2…CO2

Figure 3. 利用四組分子的二聚體進行 van der Waal 參數調整。

CH4 dimer 在進行參數調整時，不需要考慮 dipole-dipole 的作用力，故 Force Field 的部分是由 van der Waals 作用力所貢獻，藉由調整 Eq.7 Lennard Jones potential 公式中σ 和 rm參數，使MM 的計算結果逼近 DFT 的計算結果。

Eq.7

ε 為勢能阱的深度，σ 是當相互作用的是能為零食兩原子的距離，r 為兩原子距離，

rm為勢能阱底部時兩原子的距離。

H2O dimer、NH3 dimer 和 CO2 dimer 的部分，因分子間具有 dipole-dipole 的作用力，

故分子間作用力為van der Waals 作用力和 dipole-dipole 作用力的總和。當所有參 數確定後即訂下我們在分子力學中所使用的力場。

3-2 利用(CO

)

建造 CO

模型

取三種不同幾何構型的(CO2)2 ，包括水平平行、垂直平行和T-shape 如 Figure 3.所示，分別用 BLYP-D2、MP2 和 MM 這三種方法進行幾何優化後計算分子間作 用力，計算結果列於Table 1.。

(CO2)2 minimum Parallel (CO2)2 T-shape (CO2)2

Figure 4. 取三種不同結構的(CO2)2進行幾何優化，分別為水平平行、垂直平行和 T-shape 構型。

首 先 用 MM 方 法 對 水 平 平 行 的 構 型 計 算 ， 得 到 分 子 間 作 用 力 為 -1.23 kcal/mol，再用 BLYP-D2 和 MP2 方法對相同構型的(CO2)2計算，得到分子的Binding Energy 分別為-1.06 和-1.27kcal/mol；垂直平行的結構是固定兩 CO2分子間的碳原 子間距，將其固定為4 Å ，一樣利用三種不同的方法計算，發現此構型分子間會有 較大的斥力，MM 方法算得其分子間作用力約為 0.16 kcal/mol，BLYP-D2 和 MP2 方法計算所得的Binding Energy 分別為 0.25 和-0.01 kcal/mol。利用三種不同的方法 對T-shape 構型計算，MM 算得到的作用力為-0.92 kcal/mol，其他兩種方法得到的 Binding Energy 分別為-0.77 和-1.07 kcal/mol。

BLYP-D2 方法中的 D 指的是 Grimme’s dispersion，表示加入了 van der Waals 作用力的考慮項，使用 BLYP 方法計算的結果當作靜電作用力的部分，接著再將

BLYP-D2 和 BLYP 方法所得作用力相減後，即可得到單獨 van der Waals 作用力。

相同的，在進行MP2 方法計算時，以 HF 作為靜電作用力的部分，將 MP2 計算結 果和HF 相減後，亦可得到單獨 van der Waals 作用力的大小。因此如 Table 1.所列，

我們將這三種方法所得到的作用力分為靜電作用力、van der Waals 作用力和分子間 總作用力和進行討論。

比較特別的是在靜電力的部分，Gaussian 09 和 Tinker version 7.1 這兩種計算 軟體中，對於靜電作用力的部分是採用不同的計算方法: Gaussian 09 是根據庫倫定 律計算兩電子間的靜電力和電位能，如Eq.8, Eq.9；Tinker version 7.1 則是以分子 鍵為中心，將dipole-dipole 的作用力視為靜電作用力¹⁶。

Eq.8

Eq.9

利用MM 的計算方法所得的結果，和 BLYP-D2 和 MP2 這兩種方法比較可發 現，分子間總位能的部分相差的值並不多，但若以 MM 的方法分析，CO2分子間 的作用力主要來自dipole-dipole 的作用，而利用 BLYP-D2 和 MP2 這兩種方法所得 的計算結果剛好相反，CO2分子間的Binding Energy 主要的作用力來自於 van der Waals 作用力。造成此差異的原因我們推斷，或許是因為在 Gaussian 09 和 Tinker version 7.1 這兩種計算軟體中，對於靜電作用力的部分是採用不同的計算方法，因

接著藉由計算振動頻率，判斷在動能的計算上是否能獲得一致的結果。從Table 1.可見 Tinker version 7.1 的計算結果與 Gaussian 09 的計算結果相當接近，這顯示，

利用我們所得力場以MM 方法進行計算所得的動能值與 QM 的計算結果是相當接 近的。

Table 1. 取 CO2 分子，分別用 BLYP-D2, MP2, MM 這三種方法計算 intermolecular interactions 和 binding energies。

Isomer Min Parallel T-shape CO2 monomer vibrations

1BE(BLYP-D2) -1.06 0.25 -0.77 的是dipole-dipole 作用力的部分；LJ (Lennard Jones potential)則是描述 van der Waals 作用力。左上角標示1 表示:利用 BLYP-D2 和 MP2 方法優化結構；左上角標示 2 表示:利用 MM 的方法優化結構。

3-3 (CO

)

系統的相變

利用優化兩 CO2 分子，比較三種計算方法的結果後，可確定我們所使用的力 場是足夠用來描述CO2分子系統的。接著，以13 個 CO2分子在不同的環境溫度下，

進行分子動力學模擬實驗，預測系統相變時的溫度。文獻Maillet et al. 求得 13 個 CO2分子系統的熔點為95K²²，Maillet et al. 所使用的系統中 CO2分子是剛體的無 法自由振動，而本論文中CO2分子是可以自由振動的。使用6-12 Lennard-Jones 公 式計算van der Waals 作用力，如 Eq.9，並以分子鍵中心計算 dipole-dipole 作用力。

在溫度範圍為88K 到 97K 內，每個溫度下使用 NVT 架構進行分子模擬的動 力學實驗，取達到100ns 時的能量進行計算，步長為 0.5 fs，接著利用 Eq.10 求得 熱容值，Figure 4.標示的是從 88K 到 97K 的溫度範圍內計算而得的熱容值分布。

U 為(CO2)13進行分子動力學模擬所得到的能量值，R 為理想氣體常數。

Figure 5. 在 88K 到 97K 的溫度範圍，藉由分子動力學所得能量，計算(CO2)13熱熔 值的分布圖。

Eq.10

可看見在溫度為91K 時有一明顯凸起的峰，是我們量測到固態 CO2相變的溫 度，這個值和文獻Maillet et al.所得到的熔點值 95K 差不多，藉由此結果可以確定 我們在描述CO2系統的力場是適用的。

3-4 建造 MEA 模型

取兩種不同構型的(MEA)2進行討論如Figure 6.所示，這兩種構型都是因為分 子間氫鍵作用力而形成。首先將XH···H 間角度固定成 180 度，接著分別用 MM 方 法和BLYP-D2 進行幾何優化後比較其能量。

NH···H OH···O

Figure 6. 取兩不同結構 MEA 分子，分別固定 NHN 和 OHO 的角度為 180 度，用 MM 方法與 BLYP-D2 進行幾何優化。

MEA 分子間的作用力如 Table 2.所示，利用 MM 方法和 BLYP-D2 對 OHO 構 型優化得到的能量分別為-6.90 和-6.15 kcal/mol，NHN 構型得到的能量分別是-4.62 及-4.16 kcal/mol，可發現利用這兩種方法進行優化得到的結果是相同的，OHO 構 型其分子間鍵結的能力較NHN 構型強。和計算(CO2)2分子系統一樣，將BLYP-D2 方法所得的binding energy 分成兩部分討論，一是以 BLYP 方法計算的結果當作靜

電作用力的部分，而兩方法的能量差則視為van der Waals 作用力。可得到和(CO2)2 kcal/mol，OHO 構型分別為-3.61 和-4.33 kcal/mol；van der Waals 作用力的部分，

NHN 構型裡分別為-1.88 和-1.99 kcal/mol，而 OHO 構型分別為-1.4 和-1.82 kcal/mol，可發現使用這兩種方法得到的結果，無論是 dipole-dipole 作用力或是 van der Waals 作用力，都相當接近。這樣一致的結果顯示，我們使用的 MM 和 BLYP-D2 方法在進行幾何優化時是很相似的。

Table 2. 取 MEA 分子，分別用 BLYP-D2, MM 這兩種方法計算 intermolecular interactions 和 binding energies。

NHN OHO

1BE(BLYP-D2) -4.62 -6.90

Electrostatic (BLYP) -1.35 -3.50

VDW(BLYP-D2) -3.28 -3.32

potential) 則 是 描 述 van der Waals 作 用 力 。 左 上 角 標 示 1 表 示 : 利 用 BLYP-D2/aug-cc-pVTZ 方法和基組優化結構；左上角標示 2 表示:利用 MM 的方法 優化結構。

3-5 (MEA)

系統的相變

接著計算 MEA 系統發生相變時的溫度，取 27 個 MEA 分子作為反應系統，

使用NVT 架構在溫度為 280~350K 的環境下進行動力學模擬實驗。每個溫度下的 實驗進行了2.5ns 步長為 0.5fs，達平衡狀態後，取其反應位能進行 Eq.10 計算，獲 得每個溫度下系統的熱熔值，如Figure 7.所示。可清楚地看到在溫度為 305K 時有 一明顯的高峰，其為所量測到(MEA)27系統的熔點，與實驗上所得283K 相比，我 們獲得的熔點值較實驗值高出些許，但仍然在可接受的範圍內，藉此也確定了描 述MEA 系統的力場。

Figure 7. 在 280K 到 340K 的溫度範圍，藉由分子動力學所得能量，計算(MEA)27

熱熔值的分布圖。

3-6 MEA 與 CO

的鍵結

接著固定CCO2-NMEA之間的距離，一樣利用MM 和 BLYP-D2 方法，沿著 CO2

與MEA 鍵結的路徑進行計算。Figure 8.為物理吸附的相對能量和 CCO2-NMEA距離 的關係圖，紅色實線是使用MM 計算而得的曲線；藍色則是用 BLYP-D2 計算而得。

Figure 8.物理吸附的相對能量和 CCO2-NMEA距離關係圖

以 MM 方法計算時 CO2 分子與 MEA 分子的物理吸附距離為 2.9 Å ；若以 BLYP-D2 方法計算，CO2 分子與 MEA 的物理吸附距離則為 2.7Å 。由圖可見，

CCO2-NMEA距離在2.8 Å ~2.2 Å 間，兩方法計算而得的能量變化值相當一致，但當 CCO2-NMEA距離縮減至2.0Å ，MM 得到的能量和 BLYP-D2 方法所得到的能量值有 明顯的不同，ΔE 值快速上升，其值相差大約 2.5 kcal/mol。隨後，隨著 rCN距離愈 來愈接近，兩方法所得的能量變化值之間的差距，以指數的倍數增加。這表示，

我們所使用的力場足以描述在物理吸附的範圍內，CO2分子與MEA 分子的互相作

用的情形，但當兩分子愈來愈靠近接近化學吸附時，原本的力場已無法正確的描 述分子間的作用力了。

3-7 氣態 CO

溶於液態 MEA

接著我們以目前所建立的力場，進一步利用分子動力學模擬，討論氣態 CO2

與液態MEA 此兩不同相的系統反應。分別以 864 個 CO2分子組成密度為1.01g/cm³ 的氣態系統，和467 個 MEA 分子組成相同密度的液態系統，在進行幾何優化後，

取兩系統最佳結構座標配合所建立的力場進行分子動力學模擬。氣態 CO2 系統所 使用的是NVT⁶¹架構進行2.5ns，步長為 0.5fs，在常壓且環境溫度為 400K 下進行 模擬；液態MEA 系統使用相同的 NVT 架構，同樣的進行 2.5ns 模擬，步長為 0.5fs，

在常壓且環境溫度為350K 下進行模擬，並以反應的位能判定兩系統分子是否達到 平衡狀態。當皆達平衡狀態，接著將兩系統合併組合成非均相的系統。在兩系統

在常壓且環境溫度為350K 下進行模擬，並以反應的位能判定兩系統分子是否達到 平衡狀態。當皆達平衡狀態，接著將兩系統合併組合成非均相的系統。在兩系統