分子動力學(xué)模擬的發(fā)展歷史-1
1.1.1分子動力學(xué)模擬的概念
分子模擬(molecular modeling 或 molecular simulation)是一類通過計算機 模擬來研究分子或分子體系結(jié)構(gòu)與性質(zhì)的重要研究方法,包括分子力學(xué) (molecular mechanics, MM)、Monte Carlo (MC)模擬、分子動力學(xué)(molecular dynamics, MD)模擬等。這些方法均以分子或分子體系的經(jīng)典力學(xué)模型為基 礎(chǔ),或通過優(yōu)化單個分子總能量的方法得到分子的穩(wěn)定構(gòu)型(MM);或通過反 復(fù)采樣分子體系位形空間并計算其總能量的方法,得到體系的最可幾構(gòu)型與熱力 學(xué)平衡性質(zhì)(MC);或通過數(shù)值求解分子體系經(jīng)典力學(xué)運動方程的方法得到體 系的相軌跡,并統(tǒng)計體系的結(jié)構(gòu)特征與性質(zhì)(MD)。目前,得益于分子模擬理 論、方法及計算機技術(shù)的發(fā)展,分子模擬已經(jīng)成為繼實驗與理論手段之后,從分子水平了解和認(rèn)識世界的第三種手段。
1.1.2 MD模擬的早期歷史
最早的MD模擬在1957年就已實現(xiàn)。當(dāng)時,Alder和Wainwright通過計算機模擬的方法,研究了從32個到500個剛性小球分子系統(tǒng)的運動。模擬開始時, 這些小球分子被置于有序分布的格點上,具有大小相同的速度,但速度方向隨機 分布。除相互間的完全彈性碰撞外,剛性小球分子之間沒有任何相互作用,小球 分子在碰撞間隙做勻速直線運動。在經(jīng)過一段時間的模擬,系統(tǒng)中的剛性小球分 子速度達(dá)到Maxwell-Boltzmann分布后,他們分另U根據(jù)維力定理和徑向分布函數(shù) 計算了系統(tǒng)的壓力,發(fā)現(xiàn)兩種方法得到的結(jié)果一致1959年,他們提出可以 把MD模擬方法推廣到更復(fù)雜的具有方阱勢的分子體系,模擬研究分子體系的 結(jié)構(gòu)和性質(zhì)図。
1964年,Rahman模擬研究了具有Lennard-Jones勢函數(shù)的864個Ar原子 體系,得到了與狀態(tài)方程有關(guān)的性質(zhì)、徑向分布函數(shù)、速度自相關(guān)函數(shù)、均方位移等。此后,分子模擬工作者廣泛模擬研究了具有不同勢函數(shù)參數(shù)的Len- nard-Jones模型分子體系,得到了體系的結(jié)構(gòu)及其各種熱力學(xué)性質(zhì),探討了 Lennard-Jones勢函數(shù)參數(shù)對體系結(jié)構(gòu)與性質(zhì)的影響,建立了 Lennard-Jones勢 函數(shù)參數(shù)與模型分子體系結(jié)構(gòu)及性質(zhì)之間的關(guān)系。
1.1.3分子體系運動方程的數(shù)值解
在剛性小球分子體系中,除碰撞瞬間外分子間沒有任何相互作用,分子的運 動軌跡由一系列的折線組成。因此,對剛性小球分子體系的MD模擬,算法的 核心是計算剛性小球分子間的碰撞時間及其碰撞前后的運動方向和速度的變化, 而不是直接數(shù)值求解剛性小球分子體系的運動方程。相反,對于Lennard-Jones 分子,分子間一直存在相互作用,分子的運動軌跡復(fù)雜,必須通過求解體系的運 動方程來確定,這就促進(jìn)了運動方程數(shù)值方法的發(fā)展。
單原子分子沒有內(nèi)部結(jié)構(gòu),計算量小,容易實現(xiàn)MD模擬。相反,多原子 分子具有復(fù)雜的內(nèi)部結(jié)構(gòu),運動方程更加復(fù)雜,計算量大,難以實現(xiàn)MD模擬。 因此,直到20世紀(jì)70年代,才實現(xiàn)水分子體系、正烷桂分子體系等多原 子分子體系的MD模擬。
對于多原子分子體系,如果包括化學(xué)鍵伸縮勢、鍵角彎曲勢等所有分子內(nèi)的 相互作用勢都由勢函數(shù)描述,不存在所謂的約束,則體系運動方程的數(shù)值解與單 原子分子相同,沒有本質(zhì)區(qū)別。常用的MD數(shù)值積分方法包括Verlet算法、 預(yù)測-校正算法及其由Verlet算法衍生的蛙跳法凹、速度Verlet法〔回等。但是, 如果分子中的化學(xué)鍵等被約束在固定長度或整個分子被約束成剛體,則體系的運 動方程將完全不同。
根據(jù)經(jīng)典力學(xué),沒有內(nèi)部自由度的多原子分子可以被近似為 剛體,運動狀態(tài)可以分解為質(zhì)心的平動和剛體的轉(zhuǎn)動兩種獨立運動模式。其中, 質(zhì)心的平動與質(zhì)點力學(xué)沒有任何區(qū)別,不需要另外討論。剛體的轉(zhuǎn)動通常由Euler角隨時間的演化描述,可以通過數(shù)值求解相應(yīng)的運動方程得到。但是,以 Euler角為變量描述剛體運動時,牽涉奇點問題,算法不穩(wěn)定。為此,Evans提出了以四元數(shù)描述剛體運動狀態(tài)的方法,很好地克服了奇點問題,成為常用的 MD模擬算法小?⑵。
在室溫附近,大多數(shù)原子、分子處于基態(tài),分子中的化學(xué)鍵長和鍵角均在平 衡位置附近以很小的幅度振動。因此,假設(shè)化學(xué)鍵長和鍵角在MD模擬過程中 固定不變,不會對MD模擬結(jié)果產(chǎn)生顯著影響。但是,除鍵長和鍵角以外的其 他內(nèi)部運動自由度,如二面角的旋轉(zhuǎn).仍被允許,分子不能被近似為質(zhì)點或剛 體。為了利用MD模擬研究具有固定鍵長的約束體系,可以采用SHAKE算 法⑷、RATTLE算法等。約束體系的MD模擬,是MD模擬不可缺少的重要 方法,至今仍是非常熱門的研究課題。
目前,MD模擬技術(shù)已經(jīng)成熟,不但可以模擬簡單的不具內(nèi)部自由度的單原 子分子和剛性多原子分子,也可以模擬具有內(nèi)部自由度的多原子分子。即使對于 蛋白質(zhì)、DNA這樣復(fù)雜的生物大分子的模擬,也已經(jīng)沒有任何算法上的限制, 只有計算機計算能力的限制。
1. 1.4溫度與壓力的調(diào)控與統(tǒng)計系綜的實現(xiàn)
早期MD模擬體系被限制在一定的空間之內(nèi),具有固定的體積。此外,根 據(jù)能量守恒定律,模擬體系的總能量也被固定,不允許任何波動。從統(tǒng)計力學(xué)角 度分析,這樣的MD模擬體系是微正則系綜或NVE系綜。但許多實際的化學(xué)、 生物等體系,常與一恒溫?zé)嵩礋峤佑|,具有固定溫度,但不具有固定的總能量, 這就是正則系綜或NVT系綜。此外,大多數(shù)實際或人工體系,除具有固定溫度 外,也具有固定的壓力,但具有可變的總能量和體積,這是NPT系綜。
最簡單的調(diào)控體系溫度的方法是變標(biāo)度(scaling)恒溫法。當(dāng)體系的溫度, 即總動能偏離設(shè)定值時,體系中所有原子或分子的速度被乘以稱為標(biāo)度因子的系 數(shù),使體系的總動能回歸設(shè)定動能之后,Berendsen改進(jìn)了變標(biāo)度恒溫法, 使溫度逐漸被調(diào)節(jié)到設(shè)定溫度,減小了溫度的波動幅度氏]。此外,Andersen還 提出了熱浴法,在模擬過程中讓體系中的一個或若干個原子或分子與恒溫?zé)嵩粗?的分子發(fā)生隨機碰撞,達(dá)到調(diào)整體系總體溫度并使之保持恒定的目的〔閔。但是, 由這樣的方法模擬得到的相軌跡不連續(xù),不符合實際情況。
上述溫度調(diào)控算法粗糙,沒有嚴(yán)格的理論基礎(chǔ)。具有嚴(yán)格理論基礎(chǔ)的溫度調(diào) 控算法是恒溫擴展法(extended system),通過在模擬體系廣義坐標(biāo)和廣義動量 外引入額外的自由度與熱浴耦合的方法,達(dá)到調(diào)控溫度的目的。理論上,恒溫擴 展法可以通過改變或擴展實際模擬體系的哈密頓函數(shù)或拉格朗日函數(shù)實現(xiàn),如 Nose恒溫擴展法⑵如及其經(jīng)Hoover修正的Nose-Hoover恒溫算法等購。
除了調(diào)控溫度技術(shù)外,壓力的調(diào)控也是MD模擬的基本方法。與變標(biāo)度恒 溫法相似,最直接的壓力調(diào)控算法是通過調(diào)整模擬體系的體積,達(dá)到間接調(diào)控體 系壓力的目的。但是,如果在各原子的位置坐標(biāo)直接乘以同一系數(shù),將改變體系 中分子內(nèi)原子間距或鍵長,得到錯誤的結(jié)果。因此,調(diào)控壓力的算法與變標(biāo)度恒 溫法的實現(xiàn)方法不同,需通過復(fù)雜的坐標(biāo)變換,在不改變分子內(nèi)原子間相對距離 的條件下改變分子的相對位置,實現(xiàn)改變模擬體系體積的目的,這就是標(biāo)度變換 恒壓法。同樣,壓力的調(diào)控也可通過擴展系統(tǒng)自由度,采用恒壓擴展法實現(xiàn),如 Andersen恒壓恒熔系綜〔閔等。
擴展法的提出是MD理論和方法研究的一個里程碑。目前,通過擴展模擬 體系自由度,在體系的哈密頓函數(shù)或拉格朗日函數(shù)上添加額外的項,達(dá)到調(diào)控體 系的溫度與壓力的方法,仍然是一活躍的研究領(lǐng)域。
