分子動力學模擬

1、分子動力學模擬分子動力學模擬分子動力學就是一門結(jié)合物理,數(shù)學與化學的綜合技術(shù)。分子動力學就是一套分子模擬方法,該方法主要就是依靠牛頓力學來模擬分子體系的運動,以在由分子體系的不同狀態(tài)構(gòu)成的系統(tǒng)中抽取樣本 ,從而計算體系的構(gòu)型積分,并以構(gòu)型積分的結(jié)果為基礎(chǔ)進一步計算體系的熱力學量與其她宏觀性質(zhì)。這門技術(shù)的發(fā)展進程就是:1980 年:恒壓條件下的動力學方法(Andersen 方法、 Parrinello-Rahman 法)1983 年:非平衡態(tài)動力學方法(Gillan and Dixon)1984 年:恒溫條件下的動力學方法(能勢方法 )1985 年:第一原理分子動力學法( ? 法 )1991 年

2、:巨正則系綜的分子動力學方法(Cagin and Pettit) 、最新的巨正則系綜,即為組成系綜的系統(tǒng)與一溫度為T、化學勢為的很大的熱源、 粒子源相接觸 ,此時系統(tǒng)不僅同熱源有能量交換,而且可以同粒子源有粒子的交換,最后達到平衡 ,這種系綜稱巨正則系綜。進行分子動力學模擬的第一步就是確定起始構(gòu)型,一個能量較低的起始構(gòu)型就是進行分子模擬的基礎(chǔ) ,一般分子的其實構(gòu)型主要就是來自實驗數(shù)據(jù)或量子化學計算。在確定起始構(gòu)型之后要賦予構(gòu)成分子的各個原子速度,這一速度就是根據(jù)玻爾茲曼分布隨機生成,由于速度的分布符合玻爾茲曼統(tǒng)計,因此在這個階段 ,體系的溫度就是恒定的。另外 ,在隨機生成各個原子的運動速度之后

3、須進行調(diào)整,使得體系總體在各個方向上的動量之與為零,即保證體系沒有平動位移。由上一步確定的分子組建平衡相,在構(gòu)建平衡相的時候會對構(gòu)型、溫度等參數(shù)加以監(jiān)控。進入生產(chǎn)相之后體系中的分子與分子中的原子開始根據(jù)初始速度運動,可以想象其間會發(fā)生吸引、排斥乃至碰撞,這時就根據(jù)牛頓力學與預先給定的粒子間相互作用勢來對各個例子的運動軌跡進行計算,在這個過程中 ,體系總能量不變 , 但分子內(nèi)部勢能與動能不斷相互轉(zhuǎn)化,從而體系的溫度也不斷變化,在整個過程中 ,體系會遍歷勢能面上的各個點,計算的樣本正就是在這個過程中抽取的。用抽樣所得體系的各個狀態(tài)計算當時體系的勢能,進而計算構(gòu)型積分。作用勢的選擇與動力學計算的關(guān)系

4、極為密切,選擇不同的作用勢,體系的勢能面會有不同的形狀,動力學計算所得的分子運動與分子內(nèi)部運動的軌跡也會不同,進而影響到抽樣的結(jié)果與抽樣結(jié)果的勢能計算,在計算宏觀體積與微觀成分關(guān)系的時候主要采用剛球模型的二體勢,計算系統(tǒng)能量 ,熵等關(guān)系時早期多采用Lennard-Jones、morse 勢等雙體勢模型 ,對于金屬計算 ,主要采用 morse 勢 ,但就是由于通過實驗擬合的對勢容易導致柯西關(guān)系,與實驗不符 ,因此在后來的模擬中有人提出采用EAM 等多體勢模型 ,或者采用第一性原理計算結(jié)果通過一定的物理方法來擬合二體勢函數(shù)。 但就是對于二體勢模型,多體勢往往缺乏明確的表達式,參量很多 ,模擬收斂速

5、度很慢 ,給應用帶來很大困難,因此在一般應用中,通過第一性原理計算結(jié)果擬合勢函數(shù)的 L-J,morse 等勢模型的應用仍非常廣泛。分子動力學計算的基本思想就是賦予分子體系初始運動狀態(tài)之后,利用分子的自然運動在相空間中抽取樣本進行統(tǒng)計計算,時間步長就就是抽樣的間隔 ,因而時間步長的選取對動力學模擬非常重要。太長的時間步長會造成分子間的激烈碰撞,體系數(shù)據(jù)溢出 ;太短的時間步長會降低模擬過程搜索相空間的能力,因此一般選取的時間步長為體系各個自由度中最短運動周期的十分之一。但就是通常情況下,體系各自由度中運動周期最短的就是各個化學鍵的振動,而這種運動對計算某些宏觀性質(zhì)并不產(chǎn)生影響,因此就產(chǎn)生了屏蔽分子

6、內(nèi)部振動或其她無關(guān)運動的約束動力學 ,約束動力學可以有效地增長分子動力學模擬時間步長,提高搜索相空間的能分子動力學模擬力。以下就是做模擬的一般性步驟,具體的步驟與過程依賴于確定的系統(tǒng)或者軟件,但這不影響我們把它當做一個入門指南:1） 首先我們需要對我們所要模擬的系統(tǒng)做一個簡單的評估,三個問題就是我們必須要明確的:做什么 (what to do) 為什么做 (why to do) 怎么做 (how to do)2） 選擇合適的模擬工具,大前提就是它能夠?qū)崿F(xiàn)您所感興趣的目標,這需要您非常謹慎的查閱文獻 ,瞧瞧別人用這個工具都做了些什么,有沒有與您相關(guān)的,千萬不要做到一半才發(fā)現(xiàn)原來這個工具根本就不能

7、實現(xiàn)您所感興趣的idea,切記！考慮1:軟件的選擇 ,這通常與軟件主流使用的力場有關(guān),而軟件本身就具有一定的偏向性,比如說 ,做蛋白體系 ,Gromacs,Amber,Namd 均可 ;做 DNA,RNA 體系 ,首選肯定就是Amber; 做界面體系 ,DI_POLY比較強大 ,另外做材料體系,Lammps 會就是一個不錯的選擇?？紤]2:力場的選擇。力場就是來描述體系中最小單元間的相互作用的,就是用量化等方法計算擬合后生成的經(jīng)驗式,有人會嫌它粗糙,但就是它確確實實給我們模擬大系統(tǒng)提供了可能,只能說關(guān)注的切入點不同罷了。常見的有三類力場:全原子力場 ,聯(lián)合力場 ,粗?；?;當然還有所謂的第一

8、代,第二代 ,第三代力場的說法,這里就不一一列舉了。再次提醒注意 :必須選擇適合我們所關(guān)注體系與我們所感興趣的性質(zhì)及現(xiàn)象的力場。3） 通過實驗數(shù)據(jù)或者就是某些工具得到體系內(nèi)的每一個分子的初始結(jié)構(gòu)坐標文件,之后 ,我們需要按我們的想法把這些分子按照一定的規(guī)則或就是隨機的排列在一起,從而得到整個系統(tǒng)的初始結(jié)果,這也就是我們模擬的輸入文件。4） 結(jié)構(gòu)輸入文件得到了,我們還需要力場參數(shù)輸入文件,也就就是針對我們系統(tǒng)的力場文件,這通常就是由所選用的力場決定,比如鍵參數(shù)與非鍵參數(shù)等勢能函數(shù)的輸入?yún)?shù)。5） 體系的大小通常由您所選用的box 大小決定 ,我們必須對可行性與合理性做出評估,從而確定體系的大小,

9、這依賴于具體的體系,這里不細說了。6） 由于初始構(gòu)象可能會存在兩個原子挨的太近的情況(稱之為bad connect), 所以需要在正式模擬開始的第一步進行體系能量最小化,比較常用的能量最小化有兩種,最速下降法與共軛梯度法 ,最速下降法就是快速移除體系內(nèi)應力的好方法,但就是接近能量極小點時收分子動力學模擬斂比較慢 ,而共軛梯度法在能量極小點附近收斂相對效率高一些,所以我們一般做能量最小化都就是在最速下降法優(yōu)化完之后再用共軛梯度法優(yōu)化,這樣做能有效的保證后續(xù)模擬的進行。7） 以平衡態(tài)模擬為例,您需要設(shè)置適當?shù)哪M參數(shù),并且保證這些參數(shù)設(shè)置與力場的產(chǎn)生相一致 , 舉個簡單的例子,gromos 力場就

10、是用的范德華勢雙截斷來定范德華參數(shù)的,若您也用 gromos 力場的話也應該用雙截斷來處理范德華相互作用。常見的模擬思路就是,先在NVT 下約束住您的溶質(zhì)(劑 )做限制性模擬,這就是一個升溫的過程,當溫度達到您的設(shè)定后 ,接著做 NPT 模擬 ,此過程將調(diào)整體系的壓強進而使體系密度收斂。經(jīng)過一段時間的平衡模擬 ,在確定系統(tǒng)弛豫已經(jīng)完全消除之后 ,就可以開始取數(shù)據(jù)了 ,如何判斷體系達到平衡 ,這個問題就是比較技術(shù)性的問題 ,簡單的講可以通過以下幾種方式 ,一 ,瞧能量 (勢能 , 動能與總能 )就是否收斂;二 ,瞧系統(tǒng)的壓強,密度等等就是否收斂;三瞧系統(tǒng)的 RMSD 就是否達到您能接受的范圍,等等。8） 運行足夠長時間的模擬以確定我們所感興趣的現(xiàn)象或就是性質(zhì)能夠被觀測到,并且務(wù)必確保此現(xiàn)象的可重復性。9） 數(shù)據(jù)拿到手后,很容易通過一些可視化軟件得到軌跡動畫,但這并不能拿來發(fā)文章。真正的工作才剛剛開始分析數(shù)據(jù),您所感興趣的現(xiàn)象或性質(zhì)只就是表