研究背景 - N-甲基乙酰胺分子上的能量耗散途徑

(d)Cr(CO)6/CHCL3。利用2 picoseconds pulse的雷射去激發溶質的CO鍵，觀察它們的振動生命周期(vibration lifetime)τ。研究發現，溶液(a)(c)(d)的τ隨熱庫溫度增高而減小，也就是溫度高能量傳遞效率會變好。然而，溶液(b)的τ會隨溫度增高上升，有別於其他溶液。為闡述此現象，讓我們觀察跟振動生命周期τ成反比的振動鬆弛常數(vibrational relaxation constant) K = ^π

ℏ ρ 〈V^{( )}〉 (n + 1) ，其中ρ 是聲子態密度，np是聲子佔據數(phonon occupation number)，<V>是和系統位能有關的非諧偶和矩陣。CCL4與CHCL3是兩個很相似的溶劑，彼此有著相近的ρ。隨著溫度升高，兩液體的液體密度會下降，造成分子間交互作用下降，使得Debye frequency的減小，因此遠離Debye frequency的聲子態密度ρ會隨溫度升高而下降 [2]。在以W(CO)6為溶質的溶液中，溶劑CCL4的ρ下降幅度並不足以蓋過np的增

加幅度，然而溶劑CHCL3的ρ下降幅度會蓋過np的增加幅度，因為CHCL3的np增加較緩和。後者是人們認為(b)的τ為何隨溫度增高上升的原因[2]。

1999 年 C.W.Rella 等人使用自由電子雷射(free electron laser，FEL)觀察血清肌紅素 (Myoglobin)的醯胺(Amid I)鍵，FEL 的脈衝半高寬大約是 1ps，樣品是結晶狀的 Myoglobin，溶劑是 D2O 與甘油(glycer(ol-d3))以 1：4 體積比混合而成。實驗指出在不同溫度下，傅立葉轉換紅外光譜儀吸收光譜(FTIR absorption spectra)的形狀寬度差不多，且醯胺鍵振動頻率位置並沒有太大的改變，暗示著醯胺鍵結對溫度相依性低。在此系統熱庫的頻率處在約 500[cm ]以下的低頻區，高於此區被稱為高頻區。理論上溫度升高，熱庫模態(bath mode)增多，振動的傳遞加快，能量 relax 的速率會增加。但既然溫度對於鍵結吸收光譜沒有影響，顯示高頻 modes 不是直接轉成熱庫模態，也就是該鍵的振動能量會先在高頻 modes 之間傳遞[3]。

2009年Shigeto等人以三階段模型來描述分子能量耗散過程經過的途徑。第一個階段是高頻振動模態(vibration mode)的能量傳至中、低頻的振動模態。第二階段為中頻振動模態能量傳至與低頻振動模態與熱庫。第三階段為低頻振動模態能量傳至熱庫。為驗證此模型，實驗上使用超快雷射脈衝，將一道紅外脈衝分別打在甘胺酸(glycine, GLY)、N-甲基乙酰胺(N-methylacetamide, NMA)、及苯甲酸離子 (Benzoate anion, BZ)樣品上，激發碳氫鏈(chain)並測量anti-Stokes Raman頻譜以觀察能量在不同振盪模態間的轉移。實驗發現GLY與BZ高頻振盪模態的能量直接耗散到熱庫的量遠高於耗散到其他低頻振盪模態量。反之，NMA的能量是間接耗散到熱庫的，其中NMA能量由CH3模態間接傳遞到熱庫模態的總生命週期為4.9 ps[4]。

2010 年 Bastida 進一步利用非線性分子模擬來驗證分子上能量的傳輸途徑。在平衡態下瞬間置換 NMA 的甲基振動模態來模擬碳氫鍵被激發，觀察計算出的遲豫時間第一階段為 0─0.14ps，第二階段為 0.14─1.84ps，第三階段為 1.84ps─10.6ps。

在模擬中考慮了三種不同熱庫，分別為 semirigid D2O、rigid D2O、及 vacuum 環境。

在實際的 D2O 分子，氫氧鍵長能微量伸縮，氫-氧-氫的鍵角能小幅擺動。在 semirigid 模型，鍵長是固定的，但鍵角依舊可以擺動。在 rigid 模型，鍵角與長度皆固定。在 vacuum 則是無熱庫的真空狀態。當能量傳遞進行到第一階段，三個組(intramolecular vibrational redistribution)與能量耗散路徑[5]。

以上能量耗散過程的模態消長資訊，在實驗上得自於時變 infrard-Raman 光譜，在

組態計算出功率譜，找出譜上特定鍵結峰值隨時間的變化，以及不同振動模態隨時間的消長，分析能量傳遞路徑和遲豫時間。在幾種不同的操作方式下，我們得到很穩定一致的遲豫時間，約為 4.5~5.0 ps，和實驗裏的遲豫時間 4.9 ps[4]及瞬時正則模模擬的總耗散時間 11.7 ps [5]很相近。「時變功率譜」分析成為描述分子內能量耗散的工具，和最近實驗的 anti-Stokes Raman 及理論的瞬時正則模分析相輔相成。

第二章理論背景

2.1 速度相關函數(Velocity Auto-correlation function)

速度相關函數可用來探討原子與分子的振動及擴散性質。如果將速度向量函數寫成V t( )，速度相關函數(autocorrelation function)定義成

2.2 傅立葉變換( Fourier Transformation)

傅立葉變換是一種線性的積分變換。它可將時間域(time domain)函數轉換成方便分析的頻率域(frequency domain)函數。任何周期函數都可以用正弦函數和餘弦函數構成的無窮級數來表示，它就是所謂的傅立葉級數。假設訊號 f x( )為一連續

6 Fourier Transformation)，簡稱 FFT。最常見的快傅立葉變換演算法是 1965 年 Cooley 及 Tukey 共同提出的，此方法將一個序列長度為 N 的離散傅立葉序列拆成兩個長

2.3 功率譜(Power spectrum)

譜分析是描述振盪的基本工具。根據傅立葉定理，對於周期 T 的函數x t( )，我們

9 其中可看出 S(w)為x(t)的功率譜密度(power spectral density)。很明顯，用功率譜密度 S(w)取代振幅進行譜分析比較合適。

另外 Wiener-Khintchine theorem 告訴我們，廣義平穩(wide-sense stationary)隨機過程的功率譜密度是 2.1 節裡自相關函數的傅立葉轉換:

第三章模型建構與模擬計算

在這邊簡單的介紹 N-甲基乙酰胺(N-methylacetamide)分子，N-甲基乙酰胺分子是由十二顆原子組成的分子，圖 3.1 可以看到在 N-甲基乙酰胺分子的兩側各有一個甲基(CH3)，為了方便描述，這邊將靠近氮的甲基(CH3)命名為 CH3L，靠近氧的甲基(CH3)命名為 CH3R。在一般環境下，它會形成白色針狀結晶。熔點 30.55℃

（28℃），沸點 206℃，可溶於水、乙醇、苯、醚、氯仿。它是在合成藥劑時的中間產物，可用作溶劑也可製藥。NMA 擁有一個胜肽鍵，是一個生物相關的簡單分子，常被拿來當 Model，探討分子尺度的生物相關問題。

圖 3.1 N-甲基乙酰胺結構分子示意圖，綠色的碳(C)、藍色是氮(N)、紅色是氧(O)，

灰色是氫(H)。

3.1 系統模型建構

GROMACS 是 Groningen Machine for Chemical Simulations 的縮寫，是一套研究生物分子系統與分子動力學的集合封包，可以模擬聚合物、晶體、蛋白質等分子，並有多種的能量檔案對於分子模擬分析非常方便。GROMACS 是免費軟體可以在 www.gromacs.org 下載。

本研究的系統模型分為兩種，一種是將 N-甲基乙酰胺(N-methylacetamide)分子填滿整個盒子，如圖 3.2，此系統被稱為 NMA200，另一種是在充滿水的盒子中放一顆 NMA 分子，如圖 3.2，此系統被稱為 NMA_water。

圖 3.2：NMA200 是將兩百個 NMA 分子丟入盒子中。

圖 3.3：NMA_water 將一顆 NMA 分子丟入 2nm*2nm*2nm 充滿水的盒子中。

3.1.1 NMA200 模型建構

建構NMA200模型，需要準備系統結構資訊(.gro)、力場資訊(.itp)、環境條件 (.mdp)。.gro告訴GROMACS系統起始的分子位置、分子結構。.itp記錄分子間作用力函數、鍵結角度、鍵結長度、鍵結強度。這裡我們用的力場是OPLSAA[6]。mdp 設定模型模擬時的邊界條件、系統溫度、系統壓力、模擬時間、cutoff等等，當所有資訊導入GROMACS之後會產生拓譜檔(.top)記錄系統資訊，下面將介紹模擬的運算流程。

步驟一：創造平衡系統

grompp_d -f md1.mdp -c em.gro -p topol.top -o md1.tpr mdrun_d -nt 4 -v -deffnm md1

指令 grompp 將所有資訊導入拓譜檔(topol.top)，值得注意的是如果 GROMACS 是倍精度版本(double verson)就必須在指令後加上”_d”才能執行，在上面的指令中，

-f 導入系統環境條件(.mdp)，-c 導入系統結構(.gro)，-p 產生拓譜檔(.top)，-o 生成 GROMACS 看得懂的系統檔案(.tpr)。-nt 指定 cpu 運算核心數目，-deffnm 運算完模擬之後生成檔案的檔案名稱，模擬結束之後會產生各種紀錄系統資訊的檔案，

包括系統位置與速度檔案(.trr)、系統能量(.edr)、系統結構(.gro)、系統運算資訊(.log)、

系統狀態(.cpt)。系統能量包括動能、位能、總能、庫倫作用力、系統溫度、系統壓力、鍵能等等。系統運算資訊包括系統運算速率、系統運算效率、系統各種參數、系統各種初始值、總運算時間，直接當作文件打開就可以取得資訊。

make_ndx –f md1.gro

make_ndx 指令可以標識出系統原子或是將不同分子結合成群組，若加上-o 可指定生成檔名，預設生成檔案為 index.ndx。

步驟二：將所有的 NMA 的 CH3L 加熱至 1000K

grompp_df – md3.mdp -c md1.gro -t md1.cpt -n index.ndx -p topol.top -o md3.tpr mdrun_d -nt 4 -v -deffnm md3

在 GROMACS 裡可以將模擬分段，若想要將已計算出的結果做為起始條件，只需要將前一個系統結構(.gro)、系統狀態(.cpt)利用 grompp 指令引入 tpr 即可。

步驟三：加熱 NMA 後的 relaxation

grompp_df – md4.mdp -c md3.gro -t md3.cpt -n index.ndx -p topol.top -o md4.tpr mdrun_d -nt 4 -v -deffnm md4

步驟四：平衡態

grompp_df – md2.mdp -c md3.gro -t md1.cpt -n index.ndx -p topol.top -o md1.tpr mdrun_d -nt 4 -v -deffnm md2

指令 g_energy 可以將先前提到系統能量(.edr)內的資訊提出，方便觀察整個系統變化，唯一美中不足的是 g_energy 指令只能取得各種整個系統的平均能量，並不能對單一的局域部分做測量。

3.1.2 NMA_water 模型建構

步驟一：製造模型結構

editconf_d -f nmlgopls.gro -o box.gro -box 2.089 2.085 2.08 -center 1.0 1.0 1.0 genbox_d -cp box.gro -cs spc216.gro -o solv1.gro -p topol.top

editconf_d -f solv1.gro -o solv.gro -box 2.0 2.0 2.0 -center 1.0 1.0 1.0

指令ediconf製作結構檔，將原有的NMA分子結構引入，-box定義盒子大小，-center 定義分子位置。指令genbox可以將給定的分子填滿給定的空間，在我們的系統就是利用這個指令將水填滿整個系統空間模來擬體內環境。spc216.gro在GROMACS

裡是採用TIP3P[7]的水分子模型。

步驟二：能量最小化

grompp_d -f min.mdp -c solv.gro -p topol.top -o em.tpr -maxwarn 5 mdrun_d -nt 4 -v -deffnm em

能量最小化並不是一個力學計算，通常初始條件(粒子的位置與速度)都是亂數給定，為了讓系統有適當的位能、較完整的結構，因此需要執行能量最小化，讓系統有一個趨近於真實的起始條件。

步驟三到六則重複 NMA200 的步驟一到四。

3.2 功率譜計算(Power spectrum)

為了之後方便講述，先來定義一些符號:

T：表示 MD 模擬到某一個時刻。

dt：表示速度相關函數的運算區間範圍。為了分析取得的模擬資料，需先讀取系統記錄分子位置與速度的.trr 檔。圖 3.4(上)為在 CH3L 記錄到的速度對時間的函數。

每一個藍色的小方框代表一個時刻 T 到時刻 T+dt 的速度函數區間。利用 g_velacc 計算時刻 T 到時刻 T+dt 的速度相關函數(velocity autocorrelation function)，圖 3.4(中) 為這些藍色區間內的速度相關函數。若將圖 3.4(中)各區間個別做傅立葉轉換 (Fourier transform)，將得到不同時刻的功率譜(power spectrum)，如圖 3.4(下)。

圖 3.4：取出 CH3L 的速度函數(上圖)，將 T→T+dt 的速度函數定義為時刻 T 的速度函數。計算時刻 T 的速度相關函數(中圖)。最後將每時刻 T 的速度相關函數做傅立葉轉換得到瞬時功率譜(下圖)。

由功率譜上看出在頻率 W=3000.9982[cm^¹]的地方有一峰值(peak)，將每一時刻峰值記錄得到峰值與時間的對應關係如圖 3.5，其中峰值定義如下:

假設在時刻 T 所記錄到的峰值高度為P T( )，例如時刻 T=0[ps]的峰值記為P(0)、時刻 T=1[ps]的峰值記為P(1)，其中P T( )可為達平衡或未達平衡的峰值高度。假如系統已達平衡，平衡時的平均峰值為 ₀

[ (0) (1).... ( )] 1 ( )

P P P P N N P i

N 

  



圖 3.5：趨向平衡過程峰值隨時間的變化圖

圖 3.6：達平衡後峰值隨時間的漲落平均值即為P 。 ₀ 由圖 3.5 及 3.6 得到P T( )與P 。利用 ₀

y (T)=P(T)-P⁰ ⁰ eq(3.1) 算得圖 3.7 的紅色圓點。利用

y=A exp(-x)

 

eq(3.2) P(T)

對紅色圓點做指數(exponential)擬合(fitting)，得圖 3.7 的紅色實線。所得的 即為 CH3L 的遲豫時間(relaxation time)。

圖 3.7：W=3000.9982[cm^¹]CH3L 峰值擬合。黑點為模擬數據點，紅線為擬合曲線。

有別於峰值高度，另外一種作法為觀察峰值底下面積變化，首先找出功率譜上

在文檔中 N-甲基乙酰胺分子上的能量耗散途徑 (頁 6-0)