在奈米科技迅速發展的今天,體積越小、重量越輕、功能越多,以及速度越快,是現代工業產品開發的重點。在微小化的工業發展過程中,研發的思維逐漸改為由下而上的分子工程,取代傳統由上而下的研發方式。而如何藉由控制組成及分子結構,產生所需要的材料特性與元件特性,是分子工程重要的研究課題之一。
知名的理查.費曼博士(Richard P. Feynman,1965年諾貝爾物理學獎得主)在他的經典著作《費曼物理學講義》中曾提到:「我們日常生活所發生的現象都源自於原子的擾動。」(Everything that living things do can be understood in terms of the jigglings and wigglings of atoms.)也就是所有巨觀的特性都可由微觀的變化來闡釋。因此在分子工程領域中,關於分子間的交互作用、分子結構變化等微觀研究十分重要。
傳統的研究方式是以巨觀的物性量測驗證分子資訊。雖然目前已有許多儀器,如電子顯微鏡、原子力顯微鏡等,可以解析分子結構,但受限於實驗設備與條件,所能提供的分子資訊仍然有限。
分子模擬(molecular simulations)是運用牛頓和統計力學理論與電腦計算,探討系統微觀結構與能量的關係。藉由模擬系統內分子間及分子內的交互作用,把分子層級的性質(如原子內能、動能等)經統計力學運算,轉換為我們熟悉的巨觀性質(如溫度、壓力、熱容、自由能等)。分子模擬可模擬定溫定壓等真實的實驗條件,從分子微觀角度闡釋實驗結果,並與實驗驗證,甚至做出預測。本文簡介了分子模擬的類型及基本理論概念,進一步討論目前廣泛運用的進階演算法,以及分子模擬在一些領域的應用發展。
分子模擬
廣義的分子模擬包括了所有的計算化學方法,可大致分為量子化學計算、分子動態模擬(molecular dynamics simulation, MD)、蒙地卡羅(Monte Carlo, MC)模擬3種常見模擬類型。在一般文獻中,分子模擬方法是指後二種。
量子化學計算 為了描述分子的結構及對應能量,需要了解其組成原子與電子之間的交互作用。量子力學是統合了波動與粒子二元特性,以矩陣力學或波動力學描述基本粒子(如電子)的運動,最著名的就是薛丁格波動方程式。
由於量子化學計算僅需基本常數,如電子電量、普朗克常數、基本粒子質量等,因此也稱作「第一原理計算」(ab initio calculation),它的結果可用來驗證實驗結果或提出預測。薛丁格方程式可以於最簡單的氫原子系統中求得精確的解析解,但在多原子系統中,由於繁複的電子∕原子核交互作用,計算過程費時並需大量的計算資源,僅能以數值方法求得近似解。
而後Hohengerg、Kohn等人提出以系統電子密度分布計算系統基態能量,進而衍生出密度泛函理論(density functional theory, DFT)來描述電子密度分布。由於電子密度的描述大幅簡化了變數量的計算,加上描述電子密度的基底函數的準確性提升,使DFT可以在一定的精確度下,較快速地計算分子能量,是現今量子化學計算應用最多的方法之一。
儘管近年電腦運算能力大幅提升,量子化學計算仍僅能處理100~1,000個原子數目的系統,遠小於實驗系統的需求,特別是複雜的凝態系統,以及大型生物分子如蛋白質、核酸分子等。然而在材料科學、觸媒催化等奈米尺度的研究上,量子化學計算已是重要的研究工具之一。目前量子化學計算常用於計算電子軌域結構與能階、分子結構、分子光譜預測、化學反應機制與反應速率、分子熱力學性質的預測等。
分子動態模擬 分子動態模擬是以古典牛頓力學描述原子位置與速度隨時間變化的運動軌跡。各原子間的作用力如果以量子化學計算,所得到的結果雖然精確,但是計算量過於龐大,並不適用於模擬複雜分子的運動。依據波恩─歐本海默假設,電子的移動速度遠高於原子核,換句話說,如果原子核位置變動,電子迅即對應變化。在這假設下,分子的結構與能量僅需要考慮原子核的位置,電子的影響則可隱含於系統的位能描述中。
在分子模擬中,系統位能可簡化為系統的所有原子座標的數學解析函數。系統位能主要分為兩類:分子內作用力,包含化學鍵的伸縮、鍵角變化、雙面角的扭轉等;非鍵結原子間作用力,包括凡德瓦力與庫倫靜電力。
描述各種作用力的解析函數,可藉由量子化學計算、實驗數據等運用半經驗方法求得各種不同分子的參數,這些參數的集合就是分子力場。
運用分子力場可大幅減少計算系統能量的時間,並通過牛頓力學快速地計算各原子所受的淨力,進而計算系統結構隨時間的變化。因此要正確地模擬真實系統,採用準確的分子力場非常重要。目前有許多研究著重於提升分子力場的準確性,也因而衍生出各種不同的分子力場,包含常用於生物分子模擬的CHARMM及AMBER、適用於有機分子的OPLS、泛用各種分子的universal force field、適用於金屬系統的embedded atom model等。
不同的分子力場會選用不同的函數形式、不同的量子化學計算方式、不同的分子特性(如局部電荷、密度、溶解熱、氣化熱等),作為半經驗法則調整參數的依據。因此在實際運用時,使用者須考量研究課題中重要的影響因素,慎選合適的分子力場。
分子動態模擬適用於研究複雜分子系統的動態特性及熱力學特性,目前常見的應用包括:複雜分子如高分子、蛋白質等的穩定結構特性;複雜系統熱力學特性,如自由能、熱容等;系統在不同溫度、壓力及組成下,密度、結構與能量的穩定特性;動態特性如擴散係數、熱導係數與黏度,以及非平衡態的流變現象;材料的機械特性,如彎曲模數、楊氏模數等。
蒙地卡羅模擬 蒙地卡羅方法最早由梅特羅波利斯(Metropolis)等人於1953年發展出來,因此蒙地卡羅法又稱作Metropolis Monte Carlo法。稍晚於1957年,Wood等人把蒙地卡羅方法應用於液體的熱力學性質運算,近一步啟蒙了近代蒙地卡羅分子模擬的應用。蒙地卡羅方法是建構在統計熱力學的基礎上,運用隨機選取並移動系統粒子,平均計算系統在靜態平衡時的各種熱力學性質。
蒙地卡羅的隨機運算,在遍歷假設下(ergodic hypothesis),某特定物理量(如溫度、能量、壓力等)的樣本平均值等同於分子動態模擬所計算的時間平均值。而依據統計力學,這平均值也等同於實驗觀測值。在統計力學中,計算某物理量的平均值需要考量到所有可能的系統構型與狀態,以及其對應的機率:依據統計力學,某狀態出現的機率遵循波茲曼分布,和系統能量有指數關係。
實際執行蒙地卡羅模擬時,首先自系統於當下狀態中隨機選定一個系統分子,並隨機變動該分子的構形或位置產生新的系統狀態。然後重新計算系統能量,並依循波茲曼分布機率決定是否接受該變動成為狀態,或拒絕該變化而回復為原狀態。在上述簡單的流程中,複雜分子系統的能量可以用前述的分子力場來計算。
和分子動態模擬相比,蒙地卡羅模擬不需計算各原子的受力,沒有使用時間參數,也不需要對時間積分計算分子的位置變化,僅需要考慮系統能量變化,並運用隨機取樣的方式,計算系統的靜態平衡熱力學性質,因此有運算簡易且快速的優點。但由於沒有時間的描述,蒙地卡羅分子模擬無法計算系統的動力學特性。
在適當的隨機變動下,蒙地卡羅模擬可以有效地進行系統取樣,因此很適合研究系統的相態特性。目前常見的應用包括:系統平衡物性,如密度、能量、熵、焓、自由能等;系統的相平衡狀態,如氣液或固液共存的相態組成;高分子系統的自組裝結構與相變化;新熱力學性質演算法的計算與驗證。
統計力學與進階取樣計算
在分子模擬計算中,其中一個重要的物性是系統的自由能變化。在熱力學的詮釋中,自由能是描述系統由反應物變化成產物(或由起始狀態轉變為終末狀態)的自發性程度:當自由能是負值時,代表反應會自發性發生;當自由能是正值時,則不會產生自發性反應。換句話說,若一系統狀態所對應的自由能越低,代表該狀態越穩定。因此,準確的自由能計算在藥物研發、蛋白質結構、材料設計等領域可以提供非常有用的資訊。
依據統計熱力學,特定系統狀態所對應的自由能與該狀態的分布機率相關。因此實際以分子模擬進行自由能運算時,會使用進階取樣方法以更有效率、精確地統計各狀態的分布。一般而言,自由能運算可概分為兩種類型:以系統分布為基礎的進階取樣法,包括傘式取樣法(umbrella sampling)、拓展動態法等;以熱力學可逆功與系統平均力為概念的演算法,如熱力學積分法、自由能擾動法等。
系統分布為基礎的進階取樣法 如前所述,自由能與系統狀態的分布機率相關,也就是當自由能較低時,該狀態的分布機率較大。在進行分子模擬時,模擬取樣到的系統狀態大多數是自由能較低的組態;而自由能較高的組態,依照波茲曼分布,取樣的次數則以指數方式降低。因此,在有限的模擬時間限制下,會因為取樣困難產生高自由能狀態的統計誤差,降低自由能運算的準確度。
在以系統分布為基礎的諸多進階取樣演算法中,常運用的一種方法─傘式取樣法─就是為了克服高自由能狀態的取樣問題。它的概念是在系統中添加額外的偏位能,運用偏位能可以把系統限制於高自由能的狀態,進而提升統計取樣,並降低誤差。在實際運用上,可以同步進行多組模擬,每一組在不同的系統狀態加入偏位能。而後運用權重分布分析法,把各組的狀態分布統合計算出系統在不同狀態的原始分布機率,進而求得系統的自由能變化。
另一方面,和傘式取樣法中使用固定的偏位能不同,metadynamics(目前並無統一譯名,本文暫譯為拓展動態法)則採用高斯形式(Gaussian form)的累積偏位能。而每經過一段很短的時間,就會依當下系統狀態添加一個小的高斯偏位能。由於系統傾向於待在低自由能的狀態,因此在拓展動態演算法的架構下,系統中的低自由能能井會逐漸被累積的偏位能所填滿,有如以沙包填滿坑洞。最終累積的偏位能可以抵銷系統的自由能影響,達到所有狀態都是均一分布,這時系統的自由能就是所累積的偏位能的負值。
拓展動態演算法由於所需要的參數較少,不需要預先設定反應過程中的系統狀態,且可以快速拓展到多維度的自由能運算,因此近年來廣泛應用於高分子科學、生物物理、材料科學等研究領域。
熱力學可逆功為概念的演算法 在熱力學的觀念裡,自由能可以對應為系統在可逆過程中所做的功。以這概念出發,可計算兩系統狀態間的自由能差距。而在熱力學積分法(thermodynamic integration)中,可近一步採用耦合變數連接兩個系統狀態。使用耦合變數可以很簡易地產生許多的系統中間態,進而可準確計算兩狀態間的自由能變化。自由能擾動法(free energy perturbation)則是以另一種形式表述自由能代表的可逆功。若使用耦合變數,系統的總自由能變化則是所有中間態間的自由能差的總和。
使用熱力學積分法或自由能擾動法,由於計算過程費時,且無法解析反應路徑,不適用於大型分子的系統。然而可以準確計算自由能的變化,且不需要設計系統的反應途徑,因此適用於計算分子間的結合自由能,並應用在藥物設計、觸媒研究等領域。
多尺度模擬
隨著科技演進以及電腦硬體的快速發展,現今分子模擬可以描述一個約105~107個粒子的系統,而文獻中常見的模擬時間介於100奈秒至10微秒,和真實的系統尺寸與時間間距如1莫耳 = 6×1023個粒子與 1分鐘 =6×107微秒,仍有很大的差距。此外,除分子模擬之外,還有很多不同的模擬方法,例如流體力學或程序模擬。
選擇哪種方法取決於所要研究的系統大小與時間尺度,而在不同尺度下,所適用的理論敘述也應隨之調整。例如想研究化學鍵,就應該考慮量子化學計算。若是研究蛋白質摺疊,就應考慮分子模擬。若是微流體領域,則應選用流體力學運算。
結合這些不同時間及空間尺度的模擬方法,探討一系統自微觀至巨觀的性質,就是多尺度模擬法的概念。馬丁.卡普拉斯(Martin Karplus)、麥可.萊維特(Michael Levitt)和艾瑞.瓦歇爾(Arieh Warshel)發展了一貫性的多尺度模擬法,連接了量子力學、分子模擬及介觀尺度模型,催化了近代多尺度模擬研究與應用的蓬勃發展,因而於2013年一同榮獲諾貝爾化學獎。
在多尺度模擬的概念中,除了需要採用不同理論對應於不同尺度之外,更應注意的是,各尺度模擬的區分並非絕對而是相互重疊的。這意味著,在大尺度系統下如程序模擬、流體模擬所倚重的熱力、流力、反應動力學等特性,可以藉由多尺度概念串聯到小尺度系統,如分子模擬和量子模擬所測得的活性係數、自由能、黏度等微觀特性。反過來說,也可仰賴小尺度模擬系統提供準確的微觀物性資料,藉以預測大尺度系統的變化。
換句話說,運用多尺度模擬概念,在理想狀況下,可以建立「先驗的(a priori)電腦模擬實驗與程序設計」,用以預測實驗結果和程序設計的最佳化。當然,目前多尺度模擬在基本理論及實際執行上,尚有許多問題需要克服,但是隨著科技的發展,多尺度模擬已逐漸應用於高分子、材料、生醫等領域。
分子模擬的應用
隨著電腦硬體的進步,以及各種進階演算法的發展,現今分子模擬已廣泛應用於生醫研究、藥物設計、化工程序,或材料科學等領域。以下簡介以分子模擬為基礎,於3種不同領域的應用與發展。
藥物設計 小分子藥物可以產生藥效,是該藥物分子與生物體內的生物分子如蛋白質等結合後產生串聯反應(cascade reaction),進而影響生理機能。而藥物分子和生物分子的結合,通常具有立體結構的專一性,就如同鑰匙與鎖孔之間的關係,只有特定結構的分子可以與特定的生物分子結合。而分子模擬可結合分子嵌合計算,預測藥物分子的結合位置與強度,更可藉以設計藥物分子的化學修飾,調整新型藥物的生理活性。
除了分子嵌合外,定量結構活性關係(quantitative structure-activity relationships, QSAR)也廣泛運用於藥物設計。QSAR是以統計的方式,建構小分子的結構與物理化學特性的對應關係。在建構分子的資料庫後,可運用QSAR快速做出分子特性的預測。由於QSAR仰賴於結構─分子物性統計的相關性,因此不能提供物理化學機制的闡釋。在實務上,可運用QSAR篩選藥物分子,接著以分子嵌合與分子模擬預測藥物的結合位置,以及結合過程中的自由能變化。
高分子與化學工程應用 在化工程序設計時,可運用分子模擬提供相關的熱力學與物理特性,如蒸氣壓、活性係數、相態分布係數等,作為蒸餾、萃取等製程設計的參考依據。這項分子模擬─物性預測─製程設計的流程,就應用了多尺度模擬的概念。此外,分子模擬預測的結果也可用於建構分子資料庫,加速程序控制的運算時間,而實驗驗證結果也可回饋修正分子模型,提升未來模擬預測的準確性。
在高分子加工的領域,也可以運用多尺度模擬,以分子模擬計算高分子鏈的構形分布,以及高分子流體的黏彈特性,如黏度、應力、鬆弛時間等。藉由這些基本的高分子流體性質,進一步運用數值方法如有限元素法進行電腦模流分析,計算流速、壓力、溫度等,作為製程最佳化的調整參數。由國人自行研發的Moldex3D,就是廣泛應用於塑膠射出成型產業的國內外知名電腦模流軟體,應用層面含括醫療器材、汽車產業、電子產品等。
材料開發 傳統上,開發新材料多使用試誤法,整體的材料開發流程費時且費工。而分子模擬在傳統材料開發上,多用於闡釋基本物理、化學機制,並應用於結構自由能、機械模數的定量計算與預測。近年來,學術界及產業界開始建構大型的材料數位資料庫,藉以加速材料開發的流程。
建構材料資料庫,其中一種方式是運用類似於QSAR的定量結構─特性關係(quantitative structure- properties relationships, QSPR),雖然QSPR無法提供明確的物理化學意義,卻提供了統計上材料篩選的機制;藉由QSPR,可由材料資料庫中選擇合適的起始材料,而後透過前述的多尺度模擬機制,預測新材料的結構與機械特性。
國際知名廠商已實際運用材料數位資料庫於新型材料開發。其中波音公司於2014年申請「Product Chemical Profile System」專利(WO2015060960),這系統鏈結了巨觀的元件特性、基礎材料特性與微觀的材料化學組成。波音公司成功利用這系統把新材料從研發到商業化時程由原本超過12年縮短至6年,並應用於新型客機的開發。在國內,工研院也建立「高分子材料之多尺度模擬設計平台」、「混合分散技術平台」等多種材料資料庫,可藉以預測材料特性,加速材料開發的時程。
近年來人工智慧在各領域的應用引起廣泛的討論。在基礎學理上,近兩年有研究結合類神經網路與進階取樣法提升自由能計算的效率。而上述分子模擬的應用,更可結合機器學習與各種的材料及分子資料庫,提升模擬預測的準確性。隨著計算速度與準確性的提升,可以預期未來分子模擬在各種研究與應用領域中將扮演更為核心、關鍵的角色。