模擬氣態二氧化碳和液態乙醇胺的分子界面

全文

(1)國立臺灣師範大學化學系碩士論文. 指導教授：蔡明剛博士 Advisor: Ming-Kang Tsai. 模擬氣態二氧化碳和液態乙醇胺的分子界面 Interfacial Simulation of CO2 Absorption by Alcoholamines. 李佳臻 Jia-Jen Li. 中華民國一百零五年七月. i.

(2) 謝誌在師大兩年的研究所生涯過得飛快，踏入新的化學領域也學習到很多從未接觸過的事物和知識，且想感謝的人真的好多! 謝謝我的指導老師蔡明剛教授，在進入實驗室後給老師帶來的問題或麻煩絕對大於我對實驗室所能提供的幫忙，但老師總是很有耐心地回答每個問題並細膩的教導，對於這從未接觸過的領域才能不過度恐懼，真的非常謝謝蔡老師!在這兩年的實驗室生活，非常謝謝黎學謙、梁哲民、詹侑得學長和學弟廖振成，從剛進入實驗室甚麼都不懂一直到編輯論文的最後階段，這一路上所給予的幫助和鼓勵，真的非常謝謝你們，特別是侑得學長和 Peter，真的謝謝你們的幫忙和鼓勵，實驗室就麻煩你們了!謝謝相處最密切的黃雨柔、黃資庭、葉相均、李鍾奇、黃怡碩同學，一起修了兩年的課程一起討論高無作業一起熬過寫論文的這段夢魘!也恭喜大家都準備踏入人生新的階段，在未來也要繼續加油。還有可愛的 Luffy，工作累時摸摸你就能放鬆好多，你的貪吃真的超級可愛，但記得不要太胖囉! 接著我要謝謝我最重要的家人，爸爸、媽媽、哥哥、妹妹、舅舅、舅媽還有三隻胖貓咪:小虎、妹妹、咪咪，謝謝你們在我念研究所的兩年時間內所付出的勞力、心力，也請你們原諒我因為忙碌時疏於照顧和關心，未來我會怒力彌補在任何方面!也要謝謝一直陪我努力的男朋友，承受我忙碌時的壞脾氣，真心感謝。要感謝的人真的太多了，也要謝謝命運的安排讓我的人生多了這兩年的精彩，謝謝所有參與的人，無論是在這兩年間或是未來的人生道路，希望大家都能保持聯繫，也祝福大家未來能一帆風順! 李佳臻謹誌中華民國 105 年 7 月 ii.

(3) 中文摘要. 隨著溫室效應加劇，近年與降低溫室氣體排放與捕捉的相關研究愈來愈盛行，其中又以二氧化碳的研究最為常見。化學家利用與化學相關方式進行二氧化碳捕捉，並以醇胺類分子作為主要的捕捉劑。本論文欲利用分子動力學模擬實驗，針對氣態 Carbon dioxide 系統與液態 Monoethanolamines 系統的氣液界面進行討論。使用 Gaussian 09 計算軟體中 BYLP-D2 和 MP2 方法，並以 aug-cc-pVTZ 作為計算基組，此結果作為量子力學的比較基準；分子力學部分則是使用 Tinker version 7.1 計算軟體，以 mm3 參數檔作為基礎進行參數的調整。利用所得力場分別計算兩個 CO2 分子和兩個 MEA 分子間的作用力，同時也使用 BYLP-D2 和 MP2 方法對相同分子系統做計算。接著利用分子動力學模擬實驗，獲得(CO2)13 系統和(MEA)27 系統達平衡狀態時的位能，進一步的計算出在不同溫度範圍下的 Cv 值，藉由 Cv 值的分布找出系統相變的溫度。所得熔點值與文獻 Maillet 團隊所得的熔點值相當接近，MEA 系統所得熔點和實驗上熔點雖有些差距，但相差值還在可接受的範圍內。接著固定單 CO2 分子和單 MEA 分子間 CN 的距離，沿著其吸附路徑計算我們所得力場的計算結果和 BLYP-D2 方法之間的能量差距，發現在物理吸附的範圍附近，以所得力場計算分子間能量和 BLYP-D2 方法計算的能量幾乎相同建立好力場後，接著，藉由分子力學的模擬實驗對氣態 CO2 系統和液態 MEA 系統的反應介面進行探討。在反應的過程中，我們確實看見氣態的 CO2 分子進入. iii.

(4) MEA 系統，且發現當表面上有較多 MEA 分子中的親水基團時，能捕捉到較多 CO2 分子的等有趣現象。. 關鍵字: Gaussian 09，Tinker version 7.1，分子動力學，Monoethanolamines。. iv.

(5) Abstract. As the global warming intensify, more and more studies that about reducing greenhouse gases emission and capture was developed recently. And the study about carbon dioxide is the most popular. Chemists utilize the chemical methods to capture carbon dioxide and regard the aminol amine as the major capturer. We utilize the molecular mechanics simulations to discuss the interface of gas carbon dioxide and liquid monoethanolamine. Using MP2 and BLYP-D2 methods from Gaussian 09 and adopted with aug-cc-pVTZ basis set. Regard the calculation results as the quantum mechanics datum. The molecular mechanics calculations were performed using Tinker version 7.1 and adopted mm3 parameters as groundwork for adjusting parameters. When the force field was decided we used it to work out the intermolecular interaction of CO2-dimer and MEA-dimer and used MP2 and BLYP-D2 method to calculated the same molecule system simultaneously. We obtained the potential energy of (CO2)13 and (MEA)27 system when these systems were achieved equilibrium during the molecular mechanics simulations. Further, we can calculate the heat capacity under different temperature and got a melting point of these systems. The solid-solid phase transition of (CO2)13 clusters was very close to the Maillet et al. that they got the melting point of (CO2)13 system. The other system of (MEA)27 was calculated by the same method. Even though there has a slight difference in melting point between experiment and the value we got, it was within the range that can accept. Then we fixed the distance between carbon and nitrogen atom from carbon dioxide and monoethanolamine, respectively. To calculate the difference of intermolecular v.

(6) interaction energy that used the BLYP-D2 method and our force field at molecular mechanics simulation. We found that the calculation result from our force field was fitting for the BLYP-D2 method within the scope of physical absorption. After setting up the force field, we started to run the molecular mechanics simulation and discuss the interface of gas carbon dioxide and liquid monoethanolamine. We saw the carbon dioxide molecule getting into the monoethanolamine system during the simulations and found interesting phenomenon was that more hydrophilic group could capture more carbon dioxide molecules on our surface of monoethanolamine system.. Key words: Gaussian 09，Tinker version 7.1，Molecular Dynamic，Monoethanolamines。. vi.

(7) 總目錄謝誌.............................................................................................................................. ii 中文摘要 ..................................................................................................................... iii Abstract ........................................................................................................................ v 圖目錄 ......................................................................................................................... ix 表目錄 .......................................................................................................................... x 第一章緒論 ................................................................................................................ 1 1-1 前言………………………………………………………………………………….1 1-2 降低二氧化碳排放 ................................................................................................ 2 1-3 研究目的. .............................................................................................................. 4 第二章計算原理與方法 ............................................................................................. 5 2-1 計算化學的理論及方法 ........................................................................................ 5 2-1-1 分子力學 (Molecular Mechanics,MM) ..................................................... 5 2-1-1-1 分子動力學(Molecular Dynamic,MD) .................................................... 6 2-1-2 量子力學 (Quantum Mechanics,QM) ....................................................... 7 2-1-2-1 第一原理方法 (ab initio methods)......................................................... 8 2-1-2-2 半經驗法 (Semi-Empirical) ................................................................... 9 2-1-2-3 密度泛函理論 (Density Functional Theory) .......................................... 9 2-1-2-4 基底函數（Basis Sets） ....................................................................... 11 2-2 計算方法…... ....................................................................................................... 15 2-2-1 幾何優化 ................................................................................................... 15 2-2-2 單點能量(Single point energy) ................................................................ 16 vii.

(8) 2-2-3 Geometric combining rules ...................................................................... 17 2-3 Tinker version 7.152 ............................................................................................. 17 2-3-1 分子動力學模擬 ........................................................................................ 17 2-3-2 本文所使用的計算方法 ............................................................................ 18 第三章結果與討論................................................................................................... 19 3-1 建立 Force Field .................................................................................................. 19 3-2 利用(CO2)2 建造 CO2 模型 ................................................................................. 22 3-3 (CO2)13 系統的相變 ............................................................................................. 24 3-4 建造 MEA 模型 .................................................................................................. 26 3-5 (MEA)27 系統的相變 ........................................................................................... 28 3-6 MEA 與 CO2 的鍵結 ........................................................................................... 29 3-7 氣態 CO2 溶於液態 MEA................................................................................... 30 第四章結論 .............................................................................................................. 35 參考文獻 .................................................................................................................... 37. viii.

(9) 圖目錄 Figure 1. Caron Capture and Storage 簡稱 CCS 技術的流程示意圖 3。 ...................... 2 Figure 2.勢能面(Potential Energy Surface)示意圖 48。 .............................................. 16 Figure 3. 利用四組分子的二聚體進行 van der Waal 參數調整。 ............................ 21 Figure 4. 取三種不同結構的(CO2)2 進行幾何優化，分別為水平平行、垂直平行和 T-shape 構型。 ......................................................................................... 22 Figure 5. 在 88K 到 97K 的溫度範圍，藉由分子動力學所得能量，計算(CO2)13 熱熔值的分布圖。 .......................................................................................... 25 Figure 6. 取兩不同結構 MEA 分子，分別固定 NHN 和 OHO 的角度為 180 度，用 MM 方法與 BLYP-D2 進行幾何優化。.................................................. 26 Figure 7. 在 280K 到 340K 的溫度範圍，藉由分子動力學所得能量，計算(MEA)27 熱熔值的分布圖。 .................................................................................. 28 Figure 8.物理吸附的相對能量和 CCO2-NMEA 距離關係圖 ......................................... 29 Figure 9.氣態(CO2)864 系統和液態(MEA)467 系統在 t=0ps 和當 t=500ps 時的模擬實驗圖。 ......................................................................................................... 30 Figure 10.重疊四張 Nitrogen atoms 的線密度分布圖。............................................ 31 Figure 11.重疊 MEA 中 Nitrogen 和 CO2 中 Carbon 的線密度分布圖。 .................. 32 Figure 12.CCO2 與 NMEA 的徑向分布函數。 ............................................................... 34. ix.

(10) 表目錄 Table 1. 取 CO2 分子，分別用 BLYP-D2, MP2, MM 這三種方法計算 intermolecular interactions 和 binding energies。 ............................................................ 24 Table 2. 取 MEA 分子，分別用 BLYP-D2, MM 這兩種方法計算 intermolecular interactions 和 binding energies。 ............................................................ 27 Table 3. 在 200~300 的界面上，在特定時間上計算 Nitrogen, Oxygen, Carbon 和 CO2 中 Carbon 原子數量。 ............................................................................. 33 Table 4. 在 700~800 的界面上，在特定時間上計算 Nitrogen, Oxygen, Carbon 和 CO2 中 Carbon 原子數量。 ............................................................................. 33. x.

(11) 第一章緒論. 1-1 前言根據環保署的環保名詞溫室效應被定義為:從太陽輻射出來的光線原本波長較小，越過大氣層可以穿透具有和玻璃一樣效應的二氧化碳、甲烷、一氧化二氮、臭氧、氟氧碳化物等氣體而到達地球表面，然而，抵達到地球表面的陽光其反射波長較長，會被二氧化碳等氣體阻擋，不易散失於大氣外，導致地球上的溫度逐年升高，造成全球暖化。全球暖化對地球造成的影響包括:極地冰原融化，造成海平面上升，淹沒較低窪的沿海路地，因此衝擊低地國以及多數國家沿海精華區。對全球的氣候變遷也造成很大的影響，不正常的暴雨、乾旱現象更加嚴重導致沙漠化現象擴大等等，對於地球上的生態系、水土資源以及全球生物的生命安全都具有很大的威脅性，因此溫室效應已成為全世界必須共同面對的問題。聯合國在 1992 年通過「氣候變化綱要合約」(Framework Convention On Climate Change,FCCC)，期望全世界共同努力抑制溫室氣提的排；1997 年氣候變化綱要公約第三次大會中通過的「京都議定書」針對六種溫室氣體進行削減，包括二氧化碳、甲烷、氧化亞氮、氟氯碳化物、全氟碳化物及六氟化硫，其中以氫氟碳化物、全氟碳化物及六氟化硫造成溫室效應的能力最強 1，但二氧化碳由於含量較多，對全球升溫的貢獻較大。故近年來，以降低二氧化碳排放作為目標，除了發展替代石化燃料的能源、降低能源的使用量及提高能源效率、支持綠色消費等等，並有更多的科學家致力於相關的研究。. 1.

(12) 1-2 降低二氧化碳排放為了降低溫室氣體的排放，近年來相關研究如雨後春筍般湧出，特別是針對在溫室氣體中佔有很大成份的二氧化碳。二氧化碳的捕捉及封存(Caron Capture and Storage, CCS)2 是國際間新發展的技術，顧名思義指的是在發火力電廠或是石化工業這些大量製造二氧化碳排放的地方進行捕捉，進一步將其加壓、封存，如 Figure 1.所示。其中二氧化碳運送技術屬於典型的氣體輸送，技術已經相當成熟，目前 CCS 技術開發主要以捕捉存為主。. Figure 1. Caron Capture and Storage 簡稱 CCS 技術的流程示意圖 3。. 化學上最常用來捕捉 CO2 分子的方法有三種: (1). 利用電子含量豐富的金屬或具有親核基的有機分子作為 Lewis base，和 CO2 分子中視為 Lewis acid 的 C 原子反 2.

(13) 應；(2). 利用電子含量缺乏的金屬或具有親電子基的有機分子作為 Lewis acid，和 CO2 分子中視為 Lewis base 的 O 原子反應；(3). C=O 雙鍵能通過形成的錯合物接受或給予電子 4。其中又以第(1). 種方法最為常見，利用 CO2 分子中的 C 作為 Lewis acid 與胺類分子反應，藉此可將 CO2 從混雜的氣體中過濾出來。文獻中有很多關於使用胺類進行 CO2 分子捕捉的研究，例如: Feng 團隊在實驗中以乙醇胺(Monoethanolamine, MEA)作為捕捉 CO2 的分子，並加入各種不同強度的酸液，觀察其對 CO2 分子的捕捉造成的影響 5；Kim 團隊用不同濃度，分別是 30%和 70%的 MEA 水溶液進行 CO2 分子捕捉，並討論有關環境溫度對 CO2 分子吸附能的影響 6；Yang 團隊在胺類分子捕捉 CO2 的反應中使用有機化合物進行催化，藉此降低反應時的壓力也達到降低反應成本的目的 7。在使用胺類進行 CO2 捕捉時，必須面對的問題是，反應所必須消耗的能量和成本很高 8，例如分離出發電廠排放出的廢氣中的 CO2 分子，必須花費約 25~40% 的能量成本，這也是造成 CCS 技術具高成本的原因。消耗過多資源且效率不高等因素，導致很多與胺類和 CO2 反應機制相關的研究漸漸發展，故有相當多實驗或理論計算針對 CO2 與胺類分子的反應機制作探討。Xie 團隊歸納出三種 CO2 吸附於 MEA 的反應機制: (1). zwitterion mechanism9 (2) single-step mechanism10 (3) carbamic acid reaction mechanism11，其中以反應機制(1).儘管在實驗上未發現此構型，但仍為大多實驗所支持。在理論計算中支持的是其他兩個反應機制，並使用 ab initio、 DFT 和 QM/MM 方法計算支持這個理論 12，認為 zwitterion 構型為 CO2 吸附於 MEA 反應中的速率決定步驟，這和實驗上未發現 zwitterion 構型的結果相符合。. 3.

(14) 1-3 研究目的就目前文獻中的研究而言，對於 CO2 與 MEA 的反應機制在實驗和理論計算上雖然支持不同的理論，且未能得知確切的反應機制與過程，但擁有的共同目標是降低反應的成本以達更好的捕捉效率，甚至更進一步的發展出新的捕捉 CO2 分子的物質。本篇論文即針對文獻中較少研究的主題，氣態 CO2 與液態 MEA 的反應界面進行探討，對兩分子在反應時，系統表面的分子狀態能有更深入的了解，希望藉此能對新分子的設計有所幫助。同時利用分子力學作為實驗出發點，並以重現量子力學所得實驗結果作為目標。. 4.

(15) 第二章計算原理與方法. 2-1 計算化學的理論及方法計算化學的理論背景主要可分為兩大類: 1.分子力學(Molecular Mechanics) 2. 量子力學 (Quantum Mechanics) ，量子力學中又包含 : 第一原理方法 (ab initio methods)、半經驗法(Semi-Empirical)、密度泛函理論(Density Functional Theory) 13 。. 2-1-1 分子力學 (Molecular Mechanics,MM) 1930 年 Andrews 提出:分子內的化學鍵都具有自然的鍵角與鍵長值，分子內部會自動調整其幾何形狀，使鍵角與鍵長盡可能接近自然值，且分子內部還有非鍵結作用(van der Waals)，亦使其為最小狀態，如此能讓原子處於最佳的分布位置。雖然在早期就已經建立分子力學的思想和方法，但直到 50 年代以後，隨著電腦的發展，利用分子力學分析分子結構和性質的研究才日增月益。分子力學方法是利用古典力學來描述原子和化學鍵之間的關係計算出能量，利用一系列的數學函數計算描繪出分子的特徵，此函數被稱為位能函數(potential energy functions14)。其與量子力學在計算能量上較大的不同處在於，分子力學中不考慮電荷分布等細節，是利用位能函數的組合計算分子的總位能，即為力場(force field15)。在描述分子在位能面上的運動狀態時，其總能可由鍵長、鍵角、二面角、靜電力、凡德瓦作用力等相關能量的總和表示:. E = Ebond + Eangle + Etor + Eelectrostatic + Evdw 5. Eq.1.

(16) E 為力場(force field)計算的能量，Ebond 代表鍵的伸縮(Bond Stretching)，Eangle 代表鍵的彎曲(Bond Bending)，Etor 代表扭曲力(torsion)，Eelectrostatic 代表靜電作用力，Evdw 則為凡得瓦力(van der Waals force)或氫鍵(Hydrogen Bonding)。分子力場可以根據量子力學的 Born-Oppenheimer Approximation16 而得，此近似法考慮到原子中的原子核質量遠大於電子，在相同的交互作用下，電子的移動速度較原子核的速度快很多，這極大的速度差異導致電子每一刻彷彿是繞著靜止的原子核運動，故可將原子核的運動與電子的運動視為兩獨立系統，因此一個分子的能量可近似為構成分子的各個原子空間座標的函數，亦即分子的能量會隨著構型的改變而變化，描述這種分子能量和結構之間的關係就是分子力場函數。. 2-1-1-1 分子動力學(Molecular Dynamic,MD) 在古典力學分子計算中，分子是由一顆顆原子和彈簧所組成的彈簧組，其運動模式可用牛頓運動定律來表示兩原子間之位能關係。一分子在定溫定壓下經過一段時間後，各個原子的位移可透過牛頓運動定律由位能函數和合適的力場得到。藉由計算短且不連續的時間間隔使物體移動，經由依次計算的速率可得知下一步的加速度，初始速度則是隨意值，亦可利用原子間的作用力而得。藉由分子動力學能得知更多有關於分子的資訊，其可用於研究分子結構的穩定性和任何與分子運動相關的性質，例如本篇論文所進行的模擬實驗，就是以兩不同相的系統進行分子動力學模擬，進一步觀察其分子界面中分子的各種狀態與性質。. 6.

(17) 2-1-2 量子力學 (Quantum Mechanics,QM) 十七世紀牛頓的古典力學可用來解釋物體的運動行為，但隨著世代的發展，以古典力學理論已無法成功的解釋微觀粒子的運動，顯然需要一套全新適用於微觀粒子系統的理論。1990 年，Max Plank 研究出可以解釋黑體輻射的理論. 17. 。理. 論當中 Max Plank 認為黑體輻射能量的吸收和放射是因為微小的粒子被限制在不連續的值上，也就是量子化的概念。其量子理論和傳統的古典力學的理論最大的差異就在能量的描述，古典力學認為能量是連續的而量子力學認為能量是不連續的。1905 年，Albert Einstein 提出光量子的觀念且解釋了光電效應的原理 18，認為光具有粒子的行為；1923 年，De Broglie 則提出了二重性的概念，認為物質皆含有波和粒子這兩種性質，就是所謂的波動-粒子二相性(wave-particle duality)19。藉此可說，是 Albert Einstein 的光子說和 De Broglie 的物質波理論推動了量子力學的發展。量子力學是描述微觀世界的基本理論，能有效的解釋原子結構、原子光譜的規律性、化學元素的性質、光的吸收與輻射等。而今，量子化學結合電腦的使用已成為科學家不可或缺的工具，利用各種方程式計算出分子的鍵長、鍵角、偶極矩、電量分布、電子密度、軌域能量、分子間反應的能障以及分子間作用力 20；也可以藉此尋找過渡態的性質及能量，推斷反應所需的能量和反應常數等，還能用來分析未知的系統，例如化合物、表面等，或是分析化合物的 UV、IR、NMR 圖譜 20。雖然在計算上所使用的方程式相當複雜，且計算量龐大，但隨著科技日新月異的發展，利用電腦解決以上問題也變得容易許多。因此量子力學和電腦的結合運用，在研究未知的反應系統或是開創新的研究方向上都具有很大的幫助。. 7.

(18) 2-1-2-1 第一原理方法 (ab initio methods) 在量子力學中，了解一分子的能量與性質可利用薛丁格方程式(Schrödinger equation21)解出，但對於多電子的原子或分子，因其用薛丁格方程式解會過於複雜，以至於無法求得其解，而第一原理方法是使用數值求得與時間無關之薛丁格方程式的近似解。計算過程中完全不做任何假設和帶入任何實驗上的參數 22，所有計算都建立在量子力學上，此即是為甚麼這種計算方法被稱為 ab initio 的原因。理論上用第一原理方法求得的計算結果較準確，但在計算成本的消耗上卻非常多，所以只適用於計算較小分子，其計算方程式為：. Hψ = Eψ. Eq.2. E 為能量值，H 為 Hamiltonian operator，而 ψ 為波函數，主要的目的是描述電子在原子核外的位置。最常用的 ab initio 計算為 Hartree-Fock23 (簡稱 HF)，HF 方程式的基本概念為多電子體系波函數，是由系統分子軌道波函數為基礎所組成的 Slater 行列式(Slater determinant)，而系統分子軌道波函數則是由系統中所有原子軌道波函數經線性組合而成，若只改變原子波函數的係數，就能讓系統的能量達到最低點，那麼這一最低能量便是系統電子總能量的近似值，故在這一點上獲得的多電子體系波函數即為系統波函數的近似。此方法可用來計算分子的幾何結構、能量、振動頻率、IR 光譜、UV 光譜、NMR 光譜、游離能和電子親合力等等. 12. ，但其最大的限制是在. 計算時，並沒有考慮到電子相關性(Electron Correlation)，為此有許多方法來協助. 8.

(19) Hartree-Fock，常見的有 MPn (Moller-Plesset perturbation)24,25，例如 MP2、MP3、 MP4 等等。. 2-1-2-2 半經驗法 (Semi-Empirical) 半經驗法是一種電子結構的計算方法，和第一原理法的 Hartree-Fock 計算類似，但在計算過程中所使用的參數大多都是使用實驗的數據，其將部分波函數積分用參數代替，藉此可減少極大的計算量，可適用於大分子系統中，也是唯一實際使用量子力學方法於大系統分子上的計算法。但也因為如此，半經驗演算法的適用範圍會隨著所採用的實驗數據而有不同的應用範圍，故使用時必須根據研究系統的具體情況進行選擇，且因所使用的是實驗數據也導致精準度較第一原理方法下降一些。. 2-1-2-3 密度泛函理論 (Density Functional Theory) 古典電子結構理論方法中，如 Hartree-Fock 方法，是以 N 組單電子波函數構成多電子波函數由 Slater 行列式(Slater determinant)來表示，Hartree-Fock 方法可計算出全部電子波函數，但這樣在數學上的計算會非常複雜，因此，1927 年 Thomas 和 Fermi 以電子密度為變數，利用統計力學方法研究原子中電子的分布，藉此以電子密度取代電子波函數，計算總電子能量和整體電荷密度分佈，即以電子密度表示能量的理論，稱為 Thomas-Fermi Model26。但實際上此理論模型有些粗糙，更嚴重的是用在化學方面，甚至出現分子不會形成鍵結的計算結果，加上並沒有考慮到波函數的特性，所以準確度不高。. 9.

(20) 1964 年，Kohn 和 Hohenberg 證明了兩個非常重要的定理，第一定理告訴我們處於非簡併基態的多電子體系，其電子密度確定了該體系的所有物理性質，換句話說就是利用電荷密度定義基態能量 27，系統的總能量是電荷密度的泛函，所以解出一個系統的電荷密度分佈就可以知道其總能量。根據 Hohenberg-Kohn 的理論 28. ，DFT 的能量泛函可寫成下列兩項：. E[  (r )]   Vext (r )  (r )dr  F[  (r )]. Eq.3. 第一項是電子與外勢場(external potential)的相互作用，第二項是電子的動能和電子間相互作用的貢獻總和。第二理論相當於波函數量子力學中的變分原理，可證明用平均場觀念轉換成單粒子的平均勢能，將多體作用全部納入，最後可導出 Kohn-Sham 方程式 29:. F[  (r )]  EKE [  (r )]  EH [  (r )]  EXC [  (r )]. Eq.4. EKE [  (r )] 為電子的動能 (kinetic energy) 、 EH [  (r )] 為電子間的庫倫能量. (electron-electron Coulombic energy) 、 E XC [  (r ) ] 為電子相關性能量 (Electron Correlation Energy)和電子交換能量(Electron Exchange Energy)25,30； EKE [  (r )] 定義為在相同密度下系統中非相互作用的電子其動能， EH [  (r )] 是 Hatree electrostatic energy，交互關聯能項 E XC [  (r ) ]和整體系統的電子密度分布相關，對於處理多電子系統的問題而言是相當複雜，求出準確的交互關聯能項的方法目前尚未求得，. 10.

(21) 故發展出不同的泛函來解決這個問題。先前提到的第一原理法之 Hartree-fock 計算時，只考慮平均的電子密度，忽略了電子的瞬間交互作用，而 DFT 方法包含了電子交互作用的效應，雖然較花時間成本也較高，但卻有更準確的結果。本篇使用的 functional 為 BLYP-D2 (BLYP exchange–correlation functional 31,32 corrected with Grimme’s D2 dispersion33)；D 指的是 Grimme’s dispersion，由於密度泛函無法從波動電荷之間的動態相互關係正確描述 van der Waals interactions，所以在原本的 Kohn-Sham 方程式裡加入 semi-empirical dispersion potential34。. 2-1-2-4 基底函數（Basis Sets）在第一原理方法(ab initio)或是密度泛函理論(DFT)的計算中，由兩個因素來決定計算成果：一為電子相關性(Electron correlation)，另一個為基底函數(Basis Sets)。電子相關性的方面是決定要用什麼方法或是 functional 來處理波函數，稱為 operator；基底函數則是用來進行波函數的描述。基底函數為一個數學函數，每個原子有不同的基底函數，利用這樣的基底函數來描述原子的軌域 (Atomic Orbitals)，再藉由線性組合(Linear Combination)組成分子軌域(Molecular Orbitals)，像這樣的過程稱為原子軌道線性組合(Linear Combination of Atomic Orbitals， LCAO)35。常見的基底函數有兩種，一種是 STO (Slater Type Orbitals)36,37,38，另一種是 GTO (Gaussian Type Orbitals). 36,37,38. ，目前最常使用的是 GTO 類型。STO 和. GTO 分別是利用 Slater Function 和 Gaussian Function 去描述軌域:. Slater function:.   a exp(br ) 11. Eq.5.

(22) Gaussian function:.   a exp(br 2 ). Eq.6. 基底函數的發展可分為最小基底函數(Minimal Basis Set)、分裂價殼層基底函數 (Split-Valence Basis Set)39 、極化函數 (Polarization Function) 、擴散函數 (Diffuse Function)、相關組成基底函數(Correlation-Consistent Basis Set) 和 Effective Core Potentials (ECP)。最小基底函數最常見的是 STO-1G 和 STO-3G，由 S. F. Boys. 33. 發展出來，原. 理為利用 STO orbital 的形狀，分別用一個和三個 GTO 的線性組合而成，由於基組比較小，對軌域的描述程度有限；後來由 J. A. Pople 發展出分裂價殼層基底函數 39，典型形式如下：k-nmG，k 指的是中心軌域(Core Orbitals)使用 k 個高斯函數(Gaussian Function)40 的線性組合而成的基底函數，而 nm 分別是將價電子軌域(Valence Orbital) 分為兩個部分：n 為內殼層(Inner Valence Shell)、m 為外殼層(Outer Valence Shell)，代表由 n 個內殼層和 m 個外殼層的高斯函數來組成一個基底函數，常見的有 3-21G、6-31G 和 6-311G 等。以碳的 3-21G 基底函數為例子，碳的殼層可分為:. C1s: 1. (core orbitals). C2s', 2px', 2py', 2pz':  2 , 3 ,  4 , 5 C2s", 2px", 2py", 2pz":  6 ,  7 , 8 , 9. 12. (inner valence shell) (outer valence shell).

(23) 1 = a1exp(-α1r2)+a2exp(-α2r2)+a3exp(-α3r2).  2 , 3 ,  4 , 5 = b1exp(-β1r2)+b2exp(-β2r2)  6 ,  7 , 8 , 9 = c1exp(-γ1r2). 碳的 1s 軌域為中心軌域，1 由三個高斯函數的線性組合組成基底函數(因此標記成 "3")；2s'-2p'內殼層軌域，  2 - 5 分別由兩個高斯函數的線性組合組成基底函數； 2s"-2p"外殼層軌域，  6 - 9 分別由一個高斯函數的線性組合組成基底函數(標記成 "21")。如果分子結構中有第二週期之外的元素(氖以後的元素)，以 3-21G 的基底函數無法正確描述其分子軌域，電子雲的分布情形可能改變，為了要克服此問題，計算化學家利用極化函數，在價殼層中增加較高角動量的基底函數來加強彈性使其可以改變軌域形狀。極化函數是在基底函數後面加上一個星號"*"，代表非氫原子加入極化函數，第一個星號"*"為加入六個 d 形式(d-type)的基底函數，3-21G*加入 d 函數，使電子分布比較容易被極化。原子氫到氖 3-21G 和 3-21G*的基底函數是沒有差別的。若後面加兩個星號"*"，代表氫和氦原子也加入三個 p 形式(p-type)的基底函數來描述電子分布，標記成 3-21G**，也可以寫成 3-21G(d,p)。陰離子(Anions)分子有較多的電子密度分布，在高能量激發態的時候，電子離原子核較遠，使高斯函數(   a exp( br 2 ) )中的 r 變大，自然數(Exponents)值變小。因此較難以描述離原子核較遠的電子分布情形，所以應用擴散函數來修正這樣的問題。擴散函數的表示方法是在基底函數中加入一個加號"+"，一般來說和極化函數一起使用碳的 3-21G*和 3-21+G*差別在於 3-21+G*多 2s+, 2px+, 2py+和 2pz+四 13.

(24) 個基底函數。若加入兩個加號"+"，代表氫和氦原子加入擴散函數來描述電子分布的情形。通常擴散函數也可以用來敘述分子間較長距離的交互作用，像是凡得瓦力(van der Waals force)。 T. H. Dunning 團隊發展出另一種基底函數，稱為 Correlation-Consistent (cc) basis set41,42,43 ，可表示成 aug-cc-pVnZ ， aug 代表擴散函數； cc-p 為 correlation-consistent polarized 表示極化函數；V 代表的是 Valence；n 為分裂殼層數代表的為 D(double)、T(triple)、Q(quadruple)、5(quintuple)、6(sextuple)；Z 則為 Zeta。. Atoms. cc-pVDZ. cc-pVTZ. cc-Pvqz. H. 2s,1p. 3s,2p,1d. 4s,3p,2d,1f. He. 2s,1p. 3s,2p,1d. 4s,3p,2d,1f. Li-Be. 3s,2p,1d. 4s,3p,2d,1f. 5s,4p,3d,2f,1g. B-Ne. 3s,2p,1d. 4s,3p,2d,1f. 5s,4p,3d,2f,1g. Na-Ar. 4s,3p,1d. 5s,4p,2d,1f. 6s,5p,3d,2f,1g. Atoms. cc-pVDZ. Ca. 5s,4p,2d. Sc-Zn. 6s,5p,3d,1f. Ga-Kr. cc-pVTZ. cc-Pvqz. 6s,5p,3d,1f. 7s,6p,4d,2f,1g. 7s,6p,4d,2f,1g. 5s,4p,2d. 6s,5p,3d,1f. 以使用 aug-cc-pVDZ 的碳為例：. 14. 8s,7p,5d,3f,2g,1h 7s,6p,4d,2f,1g.

(25) Orbital: [1s, 2s, 2p] cc-pVDZ: [1s, 2s, 2s', 2p, 2p', d']. aug-cc-pVDZ 的碳有 3 個 s 軌域(L = 0)、2 個 p 軌域(L = 1)、1 個 d 軌域(L = 2)，共用 14 個函數。這類的基底函數比較大，著重在高階 correlation energy25,30 的計算，以及外插至 complete basis set (CBS) limit 的收斂情形，常用來做比較精準的電子相關性的計算。對於重原子分子的計算，因為過多的電子在計算電子與電子之間的作用力時太複雜，會導致產生龐大的計算量，需要縮減計算量，因此重原子會使用 Effective Core Potentials (ECP)基底函數 44 來做計算。ECP 將重原子的核電子能量給近似掉，例如:Sc-Zn、Y-Cd 和 La-Hg 分別使用[Ne]、[Ar]和[Kr]的核。一般最常見的 ECP 為 Los Alamos National Laboratory (LANL)，像是 LANL2DZ45,46 和 LANL0847,44。. 2-2 計算方法 2-2-1 幾何優化自然界中分子主要是以最穩定的形式存在，即為最低能量形式，所以我們在進行任何計算前都要先對分子進行結構的優化後，才能會進行下一步的研究。幾何優化就是根據化學鍵的計算模型找到空間中原子的排列，使每個原子與原子間的力接近於零，其產生的能量變化稱為勢能面(Potential Energy Surface，PES) 48，以座標和其能量的形式作圖進行幾何結構優化時，結構會依照位能曲面的趨勢調整至穩定點為止，目的是求極小值，表示此時的幾何結構是處於穩定狀態。. 15.

(26) Figure 2.勢能面(Potential Energy Surface)示意圖 48。. 如 Figure 2 所示，在勢能面上有幾個重要的點：最大值(global maximum)是勢能面上的最高點，該點能量是整個勢能面最大，最小值(global minimum)是勢能面上的最低點，該點能量是整個勢能面最小，也是最穩定的點，而極大值 (local maximum)是區域的高點，極小值(local minimum)是區域的低點，鞍點(saddle point) 則是勢能面上兩個高點連線和兩個低點連線的交叉點，表示是兩個低點的過渡狀態(transition state)。優化的過程中會找到最小值、極小值或鞍點，其能量對於分子結構的一階導數都為零，但鞍點上的幾何結構並非為所求的結構，因此通常會加入頻率分析 (Frequency)的指令一起計算，計算能量的二階導數來確定是否為鞍點，若是得到正頻率，即為極小值，若得到負頻率，則該點為鞍點。優化過程若不給以限制將無止盡的計算下去，因此要設收斂標準已終止計算，當計算出兩個能量值的差落在程式所設的標準內時，分子優化收斂而結束優化過程。. 2-2-2 單點能量(Single point energy) 單點能量的計算目的是為了得到一個分子的能量和基本資訊。根據所給的分子座標進行計算，可以得到分子結構的能量和相關訊息，做為下一步計算上可使 16.

(27) 用的資訊。因單點能量的計算可以在不同的 functional 採用不同基組進行計算，所以可由低階計算得到的結果進行高階的計算。. 2-2-3 Geometric combining rules Combining rules 是計算化學和分子動力學中常會使用到的方法，通常用來代表兩不雷同且未鍵結的原子，其之間的 van der Waals 作用力，且在模擬實驗中選擇不同的 Combining rules 會影響實驗的結果。常見的 Combining rules 有很多，例如: Lorentz-Berthelot rules49、Waldman-Hagler rules50、Fender-Halsey51。. 2-3 Tinker version 7.152 Tinker version 7.1 是以分子力學和分子動力學作為背景的分子模擬軟體，具有些特殊計算功能，可運用於大分子聚合物例如:蛋白質，或是小分子模擬。軟體中有多種參數檔可隨著不同系統選擇適合的參數進行計算，也可自行更改所需參數後進行實驗模擬。. 2-3-1 分子動力學模擬在進行分子動力學模擬之前，必須將所得座標優化至最佳狀態即為幾何優化，獲得優化座標後才可進一步進行模擬實驗。在 Tinker version 7.1 中會利用 Minimize52 指令進行結構的優化，輸入初步獲得的座標與配合的參數檔後即可執行。軟體將構型優化至最佳並同時製造出新的座標檔，利用此新的座標檔進行分子動力學模擬。. 17.

(28) 在軟體中執行 Dynamic52 指令進行模擬，同時需要輸入適用於系統的參數檔，接著告知軟體欲進行的步數、步長、系統環境等資訊後，即可開始動力學的模擬。本論文欲針對非均相的系統進行分子動力學模擬，兩系統分別為氣態 CO2 和液態 MEA 系統，但在進行非均相系統反應之前，兩個別系統內的分子必須為均勻分布的狀態，即達動力學平衡態，在此我們是以反應的位能作為判定依據。. 2-3-2 本文所使用的計算方法用 tinker version 7.152 進行分子動力學的實驗 53,54,55，以 mm3 參數檔作為基礎 56,57,58,59. ，針對特定的參數進行調整，目標是接近 BLYP-D2 的計算結果。系統中分. 子的參數包括 :Bond Stretching 、 Bond Bending 、 inter-atomic Lennard-Jones interactions、Bond-centered dipole-dipole interactions。Lennard-Jones interactions 的部分，是利用嘗試錯誤的方式和 geometric combining rules 分別求得 radius 和 epsilon 值，van der Waals 的 cutoff 值設成 10.0 Å ，而論文中電子結構計算的部分是利用 Gaussian 09. 60. 而得，計算方法為 MP2 和 BLYP-D2 並使用 aug-cc-pVTZ 作. 為計算基組 29,30,31。. 18.

(29) 第三章結果與討論. 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 對比版本。. 19.

(30) CH4…CH4. H2O…H2O. NH3…NH3. 20.

(31) 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 作用力的總和。當所有參數確定後即訂下我們在分子力學中所使用的力場。. 21.

(32) 3-2 利用(CO2)2 建造 CO2 模型取三種不同幾何構型的(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 方法計算的結果當作靜電作用力的部分，接著再將 22.

(33) 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 的作用力視為靜電作用力 16。. Eq.8. Eq.9 r 為兩電荷間的距離，ke 為庫倫常數。. 利用 MM 的計算方法所得的結果，和 BLYP-D2 和 MP2 這兩種方法比較可發現，分子間總位能的部分相差的值並不多，但若以 MM 的方法分析，CO2 分子間的作用力主要來自 dipole-dipole 的作用，而利用 BLYP-D2 和 MP2 這兩種方法所得的計算結果剛好相反，CO2 分子間的 Binding Energy 主要的作用力來自於 van der Waals 作用力。造成此差異的原因我們推斷，或許是因為在 Gaussian 09 和 Tinker version 7.1 這兩種計算軟體中，對於靜電作用力的部分是採用不同的計算方法，因此在計算結果上有很大的差異。 23.

(34) 接著藉由計算振動頻率，判斷在動能的計算上是否能獲得一致的結果。從 Table 1.可見 Tinker version 7.1 的計算結果與 Gaussian 09 的計算結果相當接近，這顯示，利用我們所得力場以 MM 方法進行計算所得的動能值與 QM 的計算結果是相當接近的。. Table 1. 取 CO2 分子，分別用 BLYP-D2, MP2, MM 這三種方法計算 intermolecular interactions 和 binding energies。 Isomer BE(BLYP-D2). Min -1.06. Parallel 0.25. T-shape -0.77. CO2 monomer vibrations. Electrostatic (BLYP). 0.29. 0.76. 0.09. 636, 1302, 2310. VDW(BLYP-D2). -1.36. -0.50. -0.86. (cm-1). BE(MP2). -1.27. -0.01. -1.07. Electrostatic (HF). -0.05. 0.83. -0.16. 668, 1310, 2374. VDW(MP2). -1.22. -0.84. -0.91. (cm-1). Eint(MM). -1.23. 0.16. -0.92. µµ(MM). -1.00. 0.32. -0.70. 629, 1240, 2374. LJ (MM). -0.23. -0.15. -0.22. (cm-1). 1. 1. 2. BE 為分子的 binding energies；Eint 為分子間相互作用能；在分子力學中 µµ 描述的是 dipole-dipole 作用力的部分；LJ (Lennard Jones potential)則是描述 van der Waals 作用力。左上角標示 1 表示:利用 BLYP-D2 和 MP2 方法優化結構；左上角標示 2 表示:利用 MM 的方法優化結構。. 3-3 (CO2)13 系統的相變利用優化兩 CO2 分子，比較三種計算方法的結果後，可確定我們所使用的力場是足夠用來描述 CO2 分子系統的。接著，以 13 個 CO2 分子在不同的環境溫度下， 24.

(35) 進行分子動力學模擬實驗，預測系統相變時的溫度。文獻 Maillet et al. 求得 13 個 CO2 分子系統的熔點為 95K22，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 的溫度範圍內計算而得的熱容值分布。. Eq.10. U 為(CO2)13 進行分子動力學模擬所得到的能量值，R 為理想氣體常數。. Figure 5. 在 88K 到 97K 的溫度範圍，藉由分子動力學所得能量，計算(CO2)13 熱熔值的分布圖。. 25.

(36) 可看見在溫度為 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 方法計算的結果當作靜 26.

(37) 電作用力的部分，而兩方法的能量差則視為 van der Waals 作用力。可得到和(CO2)2 分子系統一樣的結果，以 MM 的方法分析，MEA 分子間的作用力主要來自 dipole-dipole 的作用力，而利用 BLYP-D2 方法所得的計算結果剛好相反，反而是 van der Waals 作用力才是 MEA 分子間 Binding Energy 主要的作用力來來源。 Table 2 中，左上方標示 1、2 分別是用 BLYP-D2 和 MM 方法進行幾何優化所得到的分子間的作用力，NHN 構型中 dipole-dipole 作用力分別是-1.85 和-2.17 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。. BE(BLYP-D2). NHN -4.62. OHO -6.90. Electrostatic (BLYP). -1.35. -3.50. VDW(BLYP-D2). -3.28. -3.32. -1.85. -3.61. LJ(MM). -1.88. -1.40. Total(MM). -4.16. -6.15. 2. -2.17. -4.33. -1.99. -1.82. 1. 1. µµ(MM). 1 2. µµ(MM). 2. LJ(MM). BE 為分子的 binding energies；Total 為分子間相互作用能，是 µµ 和 LJ 兩作用力的總和；在分子力學中描述 µµ 的是 dipole-dipole 作用力的部分；LJ (Lennard Jones 27.

(38) potential) 則是描述 van der Waals 作用力。左上角標示 1 表示 : 利用 BLYP-D2/aug-cc-pVTZ 方法和基組優化結構；左上角標示 2 表示:利用 MM 的方法優化結構。. 3-5 (MEA)27 系統的相變接著計算 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 熱熔值的分布圖。 28.

(39) 3-6 MEA 與 CO2 的鍵結接著固定 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 分子的互相作. 29.

(40) 用的情形，但當兩分子愈來愈靠近接近化學吸附時，原本的力場已無法正確的描述分子間的作用力了。. 3-7 氣態 CO2 溶於液態 MEA 接著我們以目前所建立的力場，進一步利用分子動力學模擬，討論氣態 CO2 與液態 MEA 此兩不同相的系統反應。分別以 864 個 CO2 分子組成密度為 1.01g/cm3 的氣態系統，和 467 個 MEA 分子組成相同密度的液態系統，在進行幾何優化後，取兩系統最佳結構座標配合所建立的力場進行分子動力學模擬。氣態 CO2 系統所使用的是 NVT61 架構進行 2.5ns，步長為 0.5fs，在常壓且環境溫度為 400K 下進行模擬；液態 MEA 系統使用相同的 NVT 架構，同樣的進行 2.5ns 模擬，步長為 0.5fs，在常壓且環境溫度為 350K 下進行模擬，並以反應的位能判定兩系統分子是否達到平衡狀態。當皆達平衡狀態，接著將兩系統合併組合成非均相的系統。在兩系統中間與前後留下約 1angstrom 的距離，形成邊長為 25x25x180 的長方體盒子後，使用 NVT 架構進行 2.5ns，步長為 0.5fs，在常溫常壓下開始兩非均相系統的分子動力學模擬。Figure 9.左方為初始示意圖，右方為當反應達 500ps 時，有部分 CO2 分子已滲入 MEA 系統中。. Figure 9. t=0ps 和當 t=500ps 時氣態(CO2)864 系統和液態(MEA)467 系統的反應界面。 30.

(41) 接著我們針對兩系統反應的界面做討論。在討論界面分子之前，為了確保系統在計算過程中，是否保有原系統的品質，我們將分子動力學模擬所進行的步長平均成四等分，分別標示為圖中的 snapshot1~snapshot4，再以每等分的原子座標進行分析。將邊長為 180 angstrom 的反應箱子平均分成一千等分，如圖 Figure 10.的 x 軸所示，而其每一等分相當於 25x25x0.18 的反應箱子，接著以特定原子作為積分目標後進行線積分，獲得其分布的趨勢線。圖 Figure 10.所標示，可看見相同範圍的步長，其 Nitrogen atoms 分布線之間的差異不大，表示 CO2 在滲入 MEA 系統的過程中，MEA 系統保有原本的液體的狀態且箱子沒有太大的型變。確認完系統品質後，接著討論界面上的分子。. Figure 10.重疊四張 Nitrogen atoms 的線密度分布圖。. 31.

(42) 首先，分別求得 MEA 中的 N 和 CO2 中的 C 原子的分布趨勢線，再將兩線重疊可知道兩分子作用的界面範圍如圖 Figure 11.，因為週期性重複的關係，故會有兩分子界面，而這兩界面大約在圖中 x 軸為 200~300 和 700~800 之間. Figure 11.重疊 MEA 中 Nitrogen 和 CO2 中 Carbon 的線密度分布圖。. 接著分別計算在 0ps、20ps、40ps、100ps、150ps、200ps、250ps、300ps、350ps、 400ps、450ps、500ps 這些特定時間的界面上，MEA 分子內 N、O、C 原子的個數和 CO2 中 Carbon 原子數量，如 Table3、4.所示。我們發現在 200~300 範圍的界面上，N 原子數量隨著反應時間增加而下降，我們推測是因為 MEA 中的 NH2 基團形成氫鍵的最用力較 OH 基團強，故隨著反映進行，在界面上的 NH2 基團會迫使 MEA 分子轉向，導致在表面上 N 原子數量減少的現象。在 700~800 的範圍界面上也有相同現象發生，當反應從 100ps~150ps，N 原子數量降為零，直到反應進行到 350ps~400ps 時，N 原子數量才慢慢上升，但速度相對較緩慢。 32.

(43) Table 3. 在 200~300 的界面上，在特定時間上計算 Nitrogen, Oxygen, Carbon 和 CO2 中 Carbon 原子數量。. (200-300) Nitrogen 0ps 10 20ps 9 40ps 10 100ps 10 150ps 13 200ps 11 250ps 15 300ps 16 350ps 14 400ps 7 450ps 3 500ps 1. Oxygen 11 16 18 19 15 15 17 14 19 14 11 19. (Ni+Ox) 21 25 28 29 28 26 32 30 33 21 14 20. Carbon(MEA) Carbon(CO2) 21 36 21 27 25 28 22 25 28 38 26 28 25 39 24 31 27 49 18 44 12 38 10 28. Table 4. 在 700~800 的界面上，在特定時間上計算 Nitrogen, Oxygen, Carbon 和 CO2 中 Carbon 原子數量。. (700-800) 0ps 20ps 40ps 100ps 150ps 200ps 250ps 300ps 350ps 400ps 450ps 500ps. Nitrogen 13 17 13 13 0 0 0 0 1 5 7 8. Oxygen 11 13 10 7 9 9 10 13 7 11 11 10. (Ni+Ox) 24 30 23 20 9 9 10 13 8 16 18 18. 33. Carbon(MEA) Carbon(CO2) 28 42 30 21 28 23 20 20 16 25 16 20 16 49 14 22 16 36 21 48 22 52 20 17.

(44) 在 200~300 和 700~800 範圍的界面上，另一個有趣現象是捕捉到 CO2 分子的數量和界面上 N 原子數量與原子數量和呈現正比關係，而 MEA 中 NH2 與 OH 這兩官能基和 C 原子數相比是屬於親水性的基團，但是否界面上有較多親水性基團就能捕捉到較多 CO2 分子，這部分還沒辦法下定論，只是發現有這樣的趨勢。分析完反應系統界面後，接著利用分子動力學模擬的各個分子軌跡，計算 CO2 分子中的 C 原子與 MEA 分子中 N 原子的徑向分布函數(Radial distribution functions)，用來分析其兩分子間吸附的關係，如圖 Figure 12。. Figure 12.CCO2 與 NMEA 的徑向分布函數。圖中清楚可見，在 CCO2-NMEA 距離小於 3Å 處有一明顯的峰值，此峰顯示我們的反應系統中 CO2 分子與 MEA 分子為物理吸附的關係，但峰值卻不大，造成此現象是因為我們的動力學模擬時間不夠充裕，吸附於 MEA 分子的 CO2 分子數量不夠多所導致。隨著 CCO2-NMEA 距離增加徑向分布函數的曲線慢慢逼近 1，因當兩分子距離太遠，其彼此的作用力趨近於零，原子便以隨機的形式出現。. 34.

(45) 第四章結論. 我們選擇以分子動力學模擬實驗，以氣態 CO2 與液態 MEA 的反應界面作為探討目標。在進行模擬實驗前，以 mm3 參數檔作為基礎，主要以嘗試錯誤的方式進行參數的調整，同時與 ab initio 和 DFT 的計算結果做比較，或得到一組誤差值最小即最佳的參數後，完成分子力學中力場的建立，也成功地計算出和文獻值與實驗值非常接近的 CO2 系統與 MEA 系統的熔點，分別為 91K 和 305K，文獻中的 CO2 系統熔點為 95 K，而實驗上所測得 MEA 熔點為 283K。接著固定單 CO2 分子和單分子 MEA 間 CN 距離，利用 BLYP-D2 方法和分子力學計算方法，沿著吸附路計算兩分子間的作用力。以 DFT 計算得到兩分子物理吸附的距離為 2.9Å ，若以分子力學計算方法所得的物理吸附距離為 2.7Å ，且在物理吸附範圍兩方法的計算結果幾乎一致，但當兩分子距離逐漸靠近，約為 2.2 Å 時，兩計算方法所得的誤差值快速上升。建立好系統力場後，開始進行分子力學模擬實驗，實驗過程中確實清楚看見氣態 CO2 分子進入液態 MEA 系統，也確定反應過程中，系統並沒有太大的型變或破壞。表面上有一有趣現象，MEA 分子中屬於親水性基團的粒子數和與捕捉到 CO2 量呈現正比關係。另一有趣現象為，表面上 N 原子數量會隨著反應時間增長而降低，我們推測，NH2 基較 OH 基容易形成氫鍵，在與 CO2 分子反應的界面上，NH2 基團會較偏愛傾向於 MEA 系統本身，導致表面上的 N 原子數隨著反應時間的增加而下降。接著計算 CO2 分子中的 C 原子與 MEA 分子中 N 原子的徑向分布函數，可確認在我們的系統中 CO2 分子與 MEA 分子式以物理吸附的方式進行吸附。. 35.

(46) 我們在分子力學中建立的力場雖然在某程度上確實獲得和量子力學中 MP2, DFT 方法非常接近的計算結果，但仍然不夠完整，必須再考慮分子間更多更細的作用力或是克服軟體上的限制，以達到更好的計算結果。探討氣液界面上分子的性質或分布情況所得的資訊，希望在設計新分子進行 CO2 分子捕捉的研究上，會有些許的幫助。. 36.

(47) 參考文獻. (1). O'Neill, B. C.; Oppenheimer, M. Science 2002, 296, 1971.. (2). D'Alessandro, D. M.; Smit, B.; Long, J. R. Angewandte Chemie 2010, 49, 6058.. (3). Pires, J. C. M.; Martins, F. G.; Alvim-Ferraz, M. C. M.; Simões, M. Chemical Engineering Research and Design 2011, 89, 1446.. (4). Rochelle, G. T. Science 2009, 325, 1652.. (5). Feng, B.; Du, M.; Dennis, T. J.; Anthony, K.; Perumal, M. J. Energy & Fuels 2010, 24, 213.. (6). Kim, I.; Svendsen, H. F. Industrial & Engineering Chemistry Research 2007, 46, 5803.. (7). Yang, Z.-Z.; He, L.-N.; Gao, J.; Liu, A.-H.; Yu, B. Energy & Environmental Science 2012, 5, 6602.. (8). Haszeldine, R. S. Science 2009, 325, 1647.. (9). Caplow, M. Journal of the American Chemical Society 1968, 90, 6795.. (10). da Silva, E. F.; Svendsen, H. F. Industrial & Engineering Chemistry Research 2004, 43, 3413.. (11). Arstad, B.; Blom, R.; Swang, O. The Journal of Physical Chemistry A 2007, 111, 1222.. (12). Xie, H.-B.; Zhou, Y.; Zhang, Y.; Johnson, J. K. The Journal of Physical Chemistry A 2010, 114, 11844.. (13). Lewars, E. G. Computational Chemistry: Introduction to the Theory and Applications of Molecular and Quantum Mechanics 2nd ed., 2011.. (14). Halgren, T. A. Current Opinion in Structural Biology 1995, 5, 205.. (15). Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. Journal of the American Chemical Society 1992, 114, 10024.. (16). Woolley, R. G.; Sutcliffe, B. T. Chemical Physics Letters 1977, 45, 393.. (17). Cramer, C. J. Essentials of Computational Chemistry: Theories and Models 2nd Edition., 2004.. (18). Adawi, I. Physical Review 1964, 134, A788.. (19). In Particle Metaphysics: A Critical Account of Subatomic Reality; Falkenburg, 37.

(48) B., Ed.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2007, p 265. (20). I. N. Levine 2013.. (21). Feit, M. D.; Fleck, J. A.; Steiger, A. Journal of Computational Physics 1982, 47, 412.. (22). Young, D. C. 2001.. (23). Baerends, E. J.; Ellis, D. E.; Ros, P. Chemical Physics 1973, 2, 41.. (24). Møller, C.; Plesset, M. S. Physical Review 1934, 46, 618.. (25). Binkley, J. S.; Pople, J. A. International Journal of Quantum Chemistry 1975, 9, 229.. (26). Lieb, E. H.; Simon, B. Advances in Mathematics 1977, 23, 22.. (27). Hohenberg, P.; Kohn, W. Physical Review 1964, 136, B864.. (28). Parr, R. G.; Yang, W. Annual Review of Physical Chemistry 1995, 46, 701.. (29). Castro, A.; Marques, M. A. L.; Rubio, A. The Journal of Chemical Physics 2004, 121, 3425.. (30). Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. The Journal of Physical Chemistry A 2007, 111, 10439.. (31). Becke, A. D. Physical Review A 1988, 38, 3098.. (32). Lee, C.; Yang, W.; Parr, R. G. Physical Review B 1988, 37, 785.. (33). Grimme, S. Journal of Computational Chemistry 2006, 27, 1787.. (34). Foster, M. E.; Sohlberg, K. Physical chemistry chemical physics : PCCP 2010, 12, 307.. (35). Pipek, J.; Mezey, P. G. The Journal of Chemical Physics 1989, 90, 4916.. (36). Gill, P. M. W. In Advances in Quantum Chemistry; John, R. S., Michael, C. Z., Eds.; Academic Press: 1994; Vol. Volume 25, p 141.. (37). Boys, S. F. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 1950, 200, 542.. (38). Gill, P. M. W.; Pople, J. A. International Journal of Quantum Chemistry 1991, 40, 753.. (39). Ditchfield, R.; Hehre, W. J.; Pople, J. A. The Journal of Chemical Physics 1971, 54, 724.. (40). Guo, H. IEEE Signal Processing Magazine 2011, 28, 134.. (41). Dunning, T. H. The Journal of Chemical Physics 1989, 90, 1007.. (42). Woon, D. E.; Dunning, T. H. The Journal of Chemical Physics 1993, 98, 1358. 38.

(49) (43). Woon, D. E.; Dunning, T. H. The Journal of Chemical Physics 1995, 103, 4572.. (44). Roy, L. E.; Hay, P. J.; Martin, R. L. Journal of Chemical Theory and Computation 2008, 4, 1029.. (45). Hay, P. J.; Wadt, W. R. The Journal of Chemical Physics 1985, 82, 270.. (46). Wadt, W. R.; Hay, P. J. The Journal of Chemical Physics 1985, 82, 284.. (47). Hay, P. J.; Wadt, W. R. The Journal of Chemical Physics 1985, 82, 299.. (48). James B. Foresman, A. F. Exploring Chemistry With Electronic Structure Methods Second Edition, 1996.. (49). Delhommelle, J.; MilliÉ , P. Molecular Physics 2001, 99, 619.. (50). Peng, Z.; Ewig, C. S.; Hwang, M.-J.; Waldman, M.; Hagler, A. T. The Journal of Physical Chemistry A 1997, 101, 7243.. (51). Schnabel, T.; Vrabec, J.; Hasse, H. Journal of Molecular Liquids 2007, 135, 170.. (52). Ponder, J. W.. (53). Kundrot, C. E.; Ponder, J. W.; Richards, F. M. Journal of Computational. (accessed Feb. 2015).. Chemistry 1991, 12, 402. (54). Dudek, M. J.; Ponder, J. W. Journal of Computational Chemistry 1995, 16, 791.. (55). Ren, P.; Ponder, J. W. The Journal of Physical Chemistry B 2003, 107, 5933.. (56). Allinger, N. L.; Yuh, Y. H.; Lii, J. H. Journal of the American Chemical Society 1989, 111, 8551.. (57). Lii, J. H.; Allinger, N. L. Journal of the American Chemical Society 1989, 111, 8566.. (58). Lii, J.-H.; Allinger, N. L. Journal of Physical Organic Chemistry 1994, 7, 591.. (59). Lii, J. H.; Allinger, N. L. Journal of the American Chemical Society 1989, 111, 8576.. (60). Frisch, M. J. T., G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; J. A. Montgomery, J.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; 39.

(50) Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. (61). Fraternali, F. Biopolymers 1990, 30, 1083.. 40. 2009..

(51)