為了微觀模擬體系能夠反映宏觀實驗現(xiàn)象, 需要通過周期性邊界條件對模擬對象體系進行周期性復(fù)制, 以避免在實際中并不存在的邊緣效應(yīng)(edge effects)。原則上,對于任何分子體系的理論研究都需要求解含時薛定諤方程。但在實際工作中,更關(guān)注的是原子核的運動軌跡,這樣的軌跡可以利用波恩-奧本海默近似(Born-Oppenheimer approximation),通過求解經(jīng)典力學(xué)運動方程獲得。Alder和Wainwright曾表示計算機模擬實驗會成為聯(lián)系宏觀實驗現(xiàn)象和微觀本質(zhì)的重要橋梁,在他們首次進行分子動力學(xué)模擬實驗之后10 年,法國物理學(xué)家Verlet提出了一套牛頓運動方程的積分算法,與此同時提出的還有另一套產(chǎn)生和記錄成對近鄰原子的算法,從而大大簡化了原子間相互作用的計算。這兩種算法至今仍以一些變形的形式在實踐中被廣泛應(yīng)用[1,2]。
過去幾十年開發(fā)了多種原子級模擬方法,包括晶格靜力學(xué)、晶格動力學(xué)、蒙特卡羅和分子動力學(xué)等。其中,分子動力學(xué)特別適用于塑性變形的研究,它通過一些規(guī)定的原子間相互作用勢函數(shù)的原子相互作用系統(tǒng)的牛頓方程的解,研究變形過程中的實時行為,并包括晶格的非簡諧性、內(nèi)應(yīng)力的高度不均勻,以及系統(tǒng)的瞬態(tài)響應(yīng)等方面的影響。
分子動力學(xué)主要依靠牛頓力學(xué)來模擬分子體系的運動,以在由分子體系的不同狀態(tài)構(gòu)成的系統(tǒng)中抽取樣本,從而計算體系的構(gòu)型積分,并以構(gòu)型積分的結(jié)果為基礎(chǔ)進一步計算體系的熱力學(xué)量和其他宏觀性質(zhì)。它對原子核和電子構(gòu)成的多體系統(tǒng),求解運動方程,是一種能夠解決大量原子組成的系統(tǒng)動力學(xué)問題的計算方法,不僅可以直接模擬物質(zhì)的宏觀演變特性,得出與試驗結(jié)果相符合或相近的計算結(jié)果,還可以提供微觀結(jié)構(gòu)、粒子運動以及它們和宏觀性質(zhì)關(guān)系的明確圖像,從而為新的理論和概念的發(fā)展提供有力的技術(shù)支撐。
分子動力學(xué)的對象是一個粒子系統(tǒng),系統(tǒng)中的原子間的相互作用用勢函數(shù)來描述,因此,正確選擇勢函數(shù)的類型及其參數(shù),對于模擬的結(jié)果優(yōu)劣具有重要作用。勢能函數(shù)在大多數(shù)情況將描述分子的幾何形變最大程度地簡化為僅僅使用簡諧項和三角函數(shù)來實現(xiàn);而非鍵合原子之間的相互作用,則只采用庫侖相互作用和Lennard-Jones 勢相結(jié)合來描述。其中,對于原子間相互作用力的描述通常是經(jīng)驗或半經(jīng)驗的,這樣雖然能夠提高計算效率,但無法完全揭示電子鍵合的多體性質(zhì),尤其對于缺陷附近與自身結(jié)構(gòu)和化學(xué)性有關(guān)的復(fù)雜自洽變分函數(shù)。Daw和Baskws的EAM(Embedded-atom model)勢函數(shù)在某種程度上融合了電子鍵合的多體性質(zhì),將系統(tǒng)的總勢能表示為:
式中:Fi是原子i的嵌入能函數(shù);ρi是除第i個原子以外所有原子在i處產(chǎn)生的電子云密度之和;Φij是第i個原子與第j個原子之間的對勢作用函數(shù);rij是第i個原子與第j個原子之間的距離[1]。
勢函數(shù)的可靠性主要取決于力場參數(shù)的準(zhǔn)確性,而力場參數(shù)可以通過擬合實驗觀測數(shù)據(jù)和量子力學(xué)從頭算數(shù)據(jù)得到。目前在生物大分子體系模擬中使用最為廣泛的分子力場是CHARMM力場和AMBER力場,是早期研究生物大分子的分子力場。現(xiàn)有的力場參數(shù)仍在不斷優(yōu)化之中,并且涵蓋的分子類型也在不斷擴大。粗粒化(coarse-grained)模型在計算生物物理研究中越來越引起人們的關(guān)注,由于在該模型中定義了粗粒化粒子,對應(yīng)于全原子模型中的若干原子或原子基團甚至分子,減少了體系中的粒子數(shù),使得模擬的時間和空間尺度得以大幅度提高,但同時也將丟失原子細節(jié)信息。基于這種模型的分子動力學(xué)模擬適合于研究緩慢的生物現(xiàn)象或者依賴于大組裝體的生物現(xiàn)象。
設(shè)計一個基本力場的根本原則就是要單位時間步長(time step)內(nèi)計算能量的開銷最小化,從而實現(xiàn)模擬尺度的最大化。這一點對于全原子力場尤為重要, 即便對于所謂的粗粒化模型也是一樣。特別地,如果想要進行微秒甚至毫秒級時間尺度的模擬,這一原則極其重要。
圖1顯示了分子動力學(xué)的時間和空間尺寸的反比關(guān)系,圖中從左至右依次為:(1)水,細胞的基本成分;(2)牛胰蛋白酶抑制體,一種酶,其“呼吸”行為可以在毫秒級時間尺度上進行考察;(3)核糖體,一個復(fù)雜的生物器件,可以對基因信息進行解碼并生產(chǎn)蛋白質(zhì);(4)紫細菌光合膜片段,擁有2500 萬個原子,圖中顯示的是嵌在磷脂雙分子層上的捕光復(fù)合物及光化學(xué)反應(yīng)中心[2]。
圖1 經(jīng)典分子動力學(xué)的時間的尺度關(guān)系
隨著計算機處理器速度的快速增長以及大規(guī)模并行計算架構(gòu)的發(fā)展,大規(guī)模并行化或?qū)S玫募軜?gòu)技術(shù)與可擴展分子動力學(xué)程序的結(jié)合,計算機模擬包含從位錯到基于晶界的變形機制的整個晶粒尺寸范圍,為探索材料體系的研究前沿領(lǐng)域開辟了新的途徑。
例如,William Gon?alves等人通過選用Wolf BKS(van Beest, Kramer and van Santen)勢函數(shù)來描述原子之間的相互作用力,使用大規(guī)模原子/分子并行模擬器LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)對二氧化硅氣凝膠的彈性和強度進行了分子動力學(xué)方面的研究,他們使用velocity-Verlet算法和1.0 fs時間步長,并且在三個方向上均使用周期性邊界條件[3]。
圖2為模擬的超過7000000個原子的大體積樣品3D示意圖,以及20nm厚的樣品切片和局部放大圖(藍色為氧原子,紅色為硅原子),圖3(a)是對803nm3氣凝膠樣品進行單軸拉伸試驗,以獲得300K的應(yīng)力應(yīng)變曲線,(b-d)是典型韌性斷裂圖像,(e)抗拉強度與樣品體積的對數(shù)關(guān)系。他們分析認為,為了確保正確評估彈性等機械性能,模擬樣品的尺寸至少為孔徑的8倍,同時,表面積極高的二氧化硅氣凝膠需要相對低的應(yīng)變率以確保準(zhǔn)靜態(tài)條件。