見微知著─分子模擬的應用

UUU

分子模擬是重要的科學研究手段， 可應用於生醫、材料、化工、高分子等眾多領域。 但究竟電腦是如何模擬真實的分子？ 而了解微小分子層次的資訊又如何幫助不同領域的研究呢？

在奈米科技迅速發展的今天，體積越小、重量越輕、功能越多，以及速度越快，是現代工業產品開發的重點。在微小化的工業發展過程中，研發的思維逐漸改為由下而上的分子工程，取代傳統由上而下的研發方式。而如何藉由控制組成及分子結構，產生所需要的材料特性與元件特性，是分子工程重要的研究課題之一。

知名的理查．費曼博士（Richard P. Feynman，1965年諾貝爾物理學獎得主）在他的經典著作《費曼物理學講義》中曾提到：「我們日常生活所發生的現象都源自於原子的擾動。」（Everything that living things do can be understood in terms of the jigglings and wigglings of atoms.）也就是所有巨觀的特性都可由微觀的變化來闡釋。因此在分子工程領域中，關於分子間的交互作用、分子結構變化等微觀研究十分重要。

傳統的研究方式是以巨觀的物性量測驗證分子資訊。雖然目前已有許多儀器，如電子顯微鏡、原子力顯微鏡等，可以解析分子結構，但受限於實驗設備與條件，所能提供的分子資訊仍然有限。

分子模擬（molecular simulations）是運用牛頓和統計力學理論與電腦計算，探討系統微觀結構與能量的關係。藉由模擬系統內分子間及分子內的交互作用，把分子層級的性質（如原子內能、動能等）經統計力學運算，轉換為我們熟悉的巨觀性質（如溫度、壓力、熱容、自由能等）。分子模擬可模擬定溫定壓等真實的實驗條件，從分子微觀角度闡釋實驗結果，並與實驗驗證，甚至做出預測。本文簡介了分子模擬的類型及基本理論概念，進一步討論目前廣泛運用的進階演算法，以及分子模擬在一些領域的應用發展。

UUU

分子模擬

廣義的分子模擬包括了所有的計算化學方法，可大致分為量子化學計算、分子動態模擬（molecular dynamics simulation, MD）、蒙地卡羅（Monte Carlo, MC）模擬3種常見模擬類型。在一般文獻中，分子模擬方法是指後二種。

量子化學計算 為了描述分子的結構及對應能量，需要了解其組成原子與電子之間的交互作用。量子力學是統合了波動與粒子二元特性，以矩陣力學或波動力學描述基本粒子（如電子）的運動，最著名的就是薛丁格波動方程式。

由於量子化學計算僅需基本常數，如電子電量、普朗克常數、基本粒子質量等，因此也稱作「第一原理計算」（ab initio calculation），它的結果可用來驗證實驗結果或提出預測。薛丁格方程式可以於最簡單的氫原子系統中求得精確的解析解，但在多原子系統中，由於繁複的電子∕原子核交互作用，計算過程費時並需大量的計算資源，僅能以數值方法求得近似解。

而後Hohengerg、Kohn等人提出以系統電子密度分布計算系統基態能量，進而衍生出密度泛函理論（density functional theory, DFT）來描述電子密度分布。由於電子密度的描述大幅簡化了變數量的計算，加上描述電子密度的基底函數的準確性提升，使DFT可以在一定的精確度下，較快速地計算分子能量，是現今量子化學計算應用最多的方法之一。

儘管近年電腦運算能力大幅提升，量子化學計算仍僅能處理100～1,000個原子數目的系統，遠小於實驗系統的需求，特別是複雜的凝態系統，以及大型生物分子如蛋白質、核酸分子等。然而在材料科學、觸媒催化等奈米尺度的研究上，量子化學計算已是重要的研究工具之一。目前量子化學計算常用於計算電子軌域結構與能階、分子結構、分子光譜預測、化學反應機制與反應速率、分子熱力學性質的預測等。

UUU

分子動態模擬 分子動態模擬是以古典牛頓力學描述原子位置與速度隨時間變化的運動軌跡。各原子間的作用力如果以量子化學計算，所得到的結果雖然精確，但是計算量過於龐大，並不適用於模擬複雜分子的運動。依據波恩─歐本海默假設，電子的移動速度遠高於原子核，換句話說，如果原子核位置變動，電子迅即對應變化。在這假設下，分子的結構與能量僅需要考慮原子核的位置，電子的影響則可隱含於系統的位能描述中。

在分子模擬中，系統位能可簡化為系統的所有原子座標的數學解析函數。系統位能主要分為兩類：分子內作用力，包含化學鍵的伸縮、鍵角變化、雙面角的扭轉等；非鍵結原子間作用力，包括凡德瓦力與庫倫靜電力。

描述各種作用力的解析函數，可藉由量子化學計算、實驗數據等運用半經驗方法求得各種不同分子的參數，這些參數的集合就是分子力場。

運用分子力場可大幅減少計算系統能量的時間，並通過牛頓力學快速地計算各原子所受的淨力，進而計算系統結構隨時間的變化。因此要正確地模擬真實系統，採用準確的分子力場非常重要。目前有許多研究著重於提升分子力場的準確性，也因而衍生出各種不同的分子力場，包含常用於生物分子模擬的CHARMM及AMBER、適用於有機分子的OPLS、泛用各種分子的universal force field、適用於金屬系統的embedded atom model等。

不同的分子力場會選用不同的函數形式、不同的量子化學計算方式、不同的分子特性（如局部電荷、密度、溶解熱、氣化熱等），作為半經驗法則調整參數的依據。因此在實際運用時，使用者須考量研究課題中重要的影響因素，慎選合適的分子力場。

分子動態模擬適用於研究複雜分子系統的動態特性及熱力學特性，目前常見的應用包括：複雜分子如高分子、蛋白質等的穩定結構特性；複雜系統熱力學特性，如自由能、熱容等；系統在不同溫度、壓力及組成下，密度、結構與能量的穩定特性；動態特性如擴散係數、熱導係數與黏度，以及非平衡態的流變現象；材料的機械特性，如彎曲模數、楊氏模數等。

UUU

蒙地卡羅模擬 蒙地卡羅方法最早由梅特羅波利斯（Metropolis）等人於1953年發展出來，因此蒙地卡羅法又稱作Metropolis Monte Carlo法。稍晚於1957年，Wood等人把蒙地卡羅方法應用於液體的熱力學性質運算，近一步啟蒙了近代蒙地卡羅分子模擬的應用。蒙地卡羅方法是建構在統計熱力學的基礎上，運用隨機選取並移動系統粒子，平均計算系統在靜態平衡時的各種熱力學性質。

蒙地卡羅的隨機運算，在遍歷假設下（ergodic hypothesis），某特定物理量（如溫度、能量、壓力等）的樣本平均值等同於分子動態模擬所計算的時間平均值。而依據統計力學，這平均值也等同於實驗觀測值。在統計力學中，計算某物理量的平均值需要考量到所有可能的系統構型與狀態，以及其對應的機率：依據統計力學，某狀態出現的機率遵循波茲曼分布，和系統能量有指數關係。

實際執行蒙地卡羅模擬時，首先自系統於當下狀態中隨機選定一個系統分子，並隨機變動該分子的構形或位置產生新的系統狀態。然後重新計算系統能量，並依循波茲曼分布機率決定是否接受該變動成為狀態，或拒絕該變化而回復為原狀態。在上述簡單的流程中，複雜分子系統的能量可以用前述的分子力場來計算。

和分子動態模擬相比，蒙地卡羅模擬不需計算各原子的受力，沒有使用時間參數，也不需要對時間積分計算分子的位置變化，僅需要考慮系統能量變化，並運用隨機取樣的方式，計算系統的靜態平衡熱力學性質，因此有運算簡易且快速的優點。但由於沒有時間的描述，蒙地卡羅分子模擬無法計算系統的動力學特性。

在適當的隨機變動下，蒙地卡羅模擬可以有效地進行系統取樣，因此很適合研究系統的相態特性。目前常見的應用包括：系統平衡物性，如密度、能量、熵、焓、自由能等；系統的相平衡狀態，如氣液或固液共存的相態組成；高分子系統的自組裝結構與相變化；新熱力學性質演算法的計算與驗證。

UUU

統計力學與進階取樣計算

在分子模擬計算中，其中一個重要的物性是系統的自由能變化。在熱力學的詮釋中，自由能是描述系統由反應物變化成產物（或由起始狀態轉變為終末狀態）的自發性程度：當自由能是負值時，代表反應會自發性發生；當自由能是正值時，則不會產生自發性反應。換句話說，若一系統狀態所對應的自由能越低，代表該狀態越穩定。因此，準確的自由能計算在藥物研發、蛋白質結構、材料設計等領域可以提供非常有用的資訊。

依據統計熱力學，特定系統狀態所對應的自由能與該狀態的分布機率相關。因此實際以分子模擬進行自由能運算時，會使用進階取樣方法以更有效率、精確地統計各狀態的分布。一般而言，自由能運算可概分為兩種類型：以系統分布為基礎的進階取樣法，包括傘式取樣法（umbrella sampling）、拓展動態法等；以熱力學可逆功與系統平均力為概念的演算法，如熱力學積分法、自由能擾動法等。

系統分布為基礎的進階取樣法 如前所述，自由能與系統狀態的分布機率相關，也就是當自由能較低時，該狀態的分布機率較大。在進行分子模擬時，模擬取樣到的系統狀態大多數是自由能較低的組態；而自由能較高的組態，依照波茲曼分布，取樣的次數則以指數方式降低。因此，在有限的模擬時間限制下，會因為取樣困難產生高自由能狀態的統計誤差，降低自由能運算的準確度。

在以系統分布為基礎的諸多進階取樣演算法中，常運用的一種方法─傘式取樣法─就是為了克服高自由能狀態的取樣問題。它的概念是在系統中添加額外的偏位能，運用偏位能可以把系統限制於高自由能的狀態，進而提升統計取樣，並降低誤差。在實際運用上，可以同步進行多組模擬，每一組在不同的系統狀態加入偏位能。而後運用權重分布分析法，把各組的狀態分布統合計算出系統在不同狀態的原始分布機率，進而求得系統的自由能變化。

另一方面，和傘式取樣法中使用固定的偏位能不同，metadynamics（目前並無統一譯名，本文暫譯為拓展動態法）則採用高斯形式（Gaussian form）的累積偏位能。而每經過一段很短的時間，就會依當下系統狀態添加一個小的高斯偏位能。由於系統傾向於待在低自由能的狀態，因此在拓展動態演算法的架構下，系統中的低自由能能井會逐漸被累積的偏位能所填滿，有如以沙包填滿坑洞。最終累積的偏位能可以抵銷系統的自由能影響，達到所有狀態都是均一分布，這時系統的自由能就是所累積的偏位能的負值。

拓展動態演算法由於所需要的參數較少，不需要預先設定反應過程中的系統狀態，且可以快速拓展到多維度的自由能運算，因此近年來廣泛應用於高分子科學、生物物理、材料科學等研究領域。

UUU

熱力學可逆功為概念的演算法 在熱力學的觀念裡，自由能可以對應為系統在可逆過程中所做的功。以這概念出發，可計算兩系統狀態間的自由能差距。而在熱力學積分法（thermodynamic integration）中，可近一步採用耦合變數連接兩個系統狀態。使用耦合變數可以很簡易地產生許多的系統中間態，進而可準確計算兩狀態間的自由能變化。自由能擾動法（free energy perturbation）則是以另一種形式表述自由能代表的可逆功。若使用耦合變數，系統的總自由能變化則是所有中間態間的自由能差的總和。

使用熱力學積分法或自由能擾動法，由於計算過程費時，且無法解析反應路徑，不適用於大型分子的系統。然而可以準確計算自由能的變化，且不需要設計系統的反應途徑，因此適用於計算分子間的結合自由能，並應用在藥物設計、觸媒研究等領域。