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