分子動力學(xué)課程小結(jié)

1、分子動力學(xué)模擬課程小結(jié)一分子動力學(xué)的基本原理在分子動力學(xué)模擬中，體系原子的一系列位移是通過對牛頓運(yùn)動方程積分得到的，結(jié)果是一條運(yùn)動軌跡，它表明了系統(tǒng)內(nèi)原子的位置與速度如何隨時間而發(fā)生變化。通過解牛頓第二定律的微分方程，可以獲得原子的運(yùn)動軌跡。方程如下：這個方程描述了質(zhì)量為m：的原子i在力Fi的作用下，位置矢量為時的運(yùn)動方程。其中，F(xiàn)i可以由勢函數(shù)U的梯度給出：1系統(tǒng)的溫度則與系統(tǒng)中全部原子的總動能K通過下式相聯(lián)系：N是原子數(shù)，Nc是限制條件，kB是波爾茲曼常數(shù)。MD模擬的積分算法為了得到原子的運(yùn)動軌跡，可以采用有限差分法來求解運(yùn)動方程。有限差分法的基本思想就是將積分分成很多小步，每一小步的時間

2、固定為6to用有限差分解運(yùn)動方程有許多方法，所有的算法都假定位置與動態(tài)性質(zhì)(速度、加速度等)可以用Taylor級數(shù)展開來近似：r(f+=r(l)+Av)-f；-貂口(r)+*冷需十占撐c(/)+v(t+瘵)=vg)十西J(01-#護(hù)帆F)+按垃)+詁+西)=；+因站)十卅2斗在分子動力學(xué)模擬中，常用的有以下的幾中算法：Verlet算法運(yùn)用t時刻的位置和速度及t6t時刻的位置，計(jì)算出t+6t時刻的位置：r(t+3t)=r(t)+Stv(Z)-冷護(hù)ff(/)+.r(f-&)=F(r)-合+#用a(t)二一兩式相加并忽略高階項(xiàng)，可以得到：r(t+茂)=2r(r)-占)+護(hù)a(t)速度可以通過以下方法

3、得到：用t+6t時刻與t6t時刻的位置差除以26t：同理，半時間步t+8t時刻的速度也可以算：Verlet算法執(zhí)行簡單明了，存儲要求適度，但缺點(diǎn)是位置r(t+St)要通過小項(xiàng)與非常大的兩項(xiàng)2r(t)與r(t-6t)的差相加得到，容易造成精度損失。另外，其方程式中沒有顯示速度項(xiàng)，在沒有得到下一步的位置前速度項(xiàng)難以得到。它不是一個自啟動算法：新位置必須由t時刻與前一時刻t-8t的位置得到。在t=0時刻,只有一組位置，所以必須通過其它方法得到t-8t的位置。一般用Taylor級數(shù):TfT珂一用)=鞏0)円間譏0)Velocity-Verlet算法-+1f中+冷)二F(f)+啟譏小a(t)Leap-f

4、rog算法r(f+-F(F)+&*(f+*丹為了執(zhí)行Leap-frog算法，必須首先由t-0.58t時刻的速度與t時刻的加速度計(jì)算出速度v(t+8t),然后由方程十&)=V(f)4丄按日&+扮)+：冷7一丄帶T(f一層)計(jì)算出位置r(t+8t)。T時刻的速度可以由:得到。速度蛙跳過此t時刻的位置而得到t+0.58t時刻的速度值，而位置跳過速度值給出了t+8t時刻的位置值，為計(jì)算t+0.58t時刻的速的作準(zhǔn)備，依此類推。其缺點(diǎn)是位置與速度不同步。這意味著在位置一定時，不可能同時計(jì)算動能對總能量的貢獻(xiàn)。分子動力學(xué)計(jì)算的時間間隔時間間隔8t在積分算法中是一個非常重要的參數(shù)。為了充分利用時間，盡量選擇

5、比較大的時間間隔，但是如果時間間隔太大，就會造成積分過程的不穩(wěn)定性和不精確性。時間間隔的設(shè)置同時依賴于算法和模型的情況。模型本身給時間間隔帶來的最大的限制就是最高頻率的運(yùn)動。由于算法要求在每個時間間隔內(nèi)模型的速度和加速度保持一邊，時間間隔就應(yīng)該低于振動周期的八分之一到十分之一。對大多數(shù)的有機(jī)模型來講，最高的振動頻率是鍵的伸縮振動，其振動周期的數(shù)量級為-。這樣，時間間隔就應(yīng)該是左右。如果采用受約束的或者算法，可以使用更長的時間間隔。如果研究對象是液態(tài)或者固態(tài)簡單模型，對體系內(nèi)作用模式不感興趣，也可以采用一些更長的時間間隔，比如。對離子態(tài)的材料模型，左右是合適的。時間間隔必須跟選擇的算法相匹配。比

6、如，算法的時間間隔應(yīng)該是算法的一半左算法則需要比其他算法更短的時間間隔。生物大分子的相互作用勢函數(shù)生物大分子勢一個自有度很大，結(jié)構(gòu)復(fù)雜的體系，各種物理過程極為豐富。其中的相互作用規(guī)律還難以進(jìn)行統(tǒng)一的處理。但其物理圖像是清晰的，電子的運(yùn)動和骨架有效的分開，成鍵電子是定域的，價(jià)鍵的概念依然有效。荷電集團(tuán)之間存在靜電相互作用。原子或集團(tuán)正負(fù)電荷中心不重合形成電多極矩并參與靜電相互作用。原子或集團(tuán)的電子云受靜電相互的誘導(dǎo)形成電誘導(dǎo)偶極由此產(chǎn)生色散力。靜電、范德華相互作用稱非鍵相互作用。此外有共價(jià)鍵德振動及單鍵德旋轉(zhuǎn)繞共價(jià)雙鍵德扭曲等稱為成鍵相互作用。典型德生物分子動力學(xué)的模型勢函數(shù)可以寫成一下形式：陷

7、忘矗X號叫y+另軋爐冊+e抑haiA息batd叢呼吁既址第一項(xiàng)求和是共價(jià)鍵的鍵能，第二項(xiàng)求和是鍵角能，第三項(xiàng)求和是非正常二面角能量，第四項(xiàng)求和是正常二面角能量。Kb,k9,kd,k代表式中前四項(xiàng)成鍵相互作用的力常數(shù)。式中最后一項(xiàng)式體系中的非鍵相互作用項(xiàng)，（i,j）表示對非鍵對求和，C12、C6式范德華相互作用勢中的系數(shù)。勢函數(shù)中參數(shù)qi，qj分別為第i荷第j個原子的電荷。但是，對于疏水相互作用，溶劑的極化荷對蛋白質(zhì)內(nèi)部靜電屏蔽，都需要對勢函數(shù)加以修正和補(bǔ)充。若要了解溶液中分子的動力學(xué)性質(zhì)則須將模擬的體系放到充分多的水分子做充分動力學(xué)模擬。分子動力學(xué)模擬的啟動進(jìn)行MD模擬，必須首先建立系統(tǒng)的初始

8、結(jié)構(gòu)。初始結(jié)構(gòu)可以通過實(shí)驗(yàn)數(shù)據(jù)、理論模型、或兩者的結(jié)合來獲得。除此之外，還可以給每個原子賦予初速度，它可以從一定溫度下的Maxwell-Boltzmann分布來任意選?。篗axwell-Boltzmann分布給出了質(zhì)量為mi的原子i在溫度T下沿x方向速度為v.的概率。Maxwell-Boltzmann分布勢一中Guassian分布，它可以用隨機(jī)數(shù)ix發(fā)生器得到。大多數(shù)隨機(jī)發(fā)生器產(chǎn)生的隨機(jī)數(shù)均勻分布在01之間，但是可以通過變換得到Guassian分布。均值為x和波動值為S2的Guassian分布的概率為：一種方法勢首先產(chǎn)生兩個在01之間的隨機(jī)數(shù)E1和2。運(yùn)用下列式子可產(chǎn)生兩個數(shù)下x1、x2：X,

9、=J（-2In點(diǎn)）昌）心“4（-2In昌）sioO嘖）另一種方法勢先產(chǎn)生12個隨機(jī)數(shù)1、2、12,然后計(jì)算：這兩種方法產(chǎn)生的隨機(jī)值都服從均值為零，偏差為一個單位的正態(tài)分布。初始速度經(jīng)常被校正以滿足總動量為零，為了是總動量為零，分別計(jì)算沿三方向的動量總和，然后用每一方向的總定量除以總質(zhì)量，得到一速度值。用每個原子的速的減去此速度值，即可保持系統(tǒng)的總動量為零。在建立了系統(tǒng)的初始位形和賦予初始速度后，就具備分子動力學(xué)模擬的初步條件了。在每一步中，原子所收到的力通過對勢函數(shù)的微分可以得到。然后根據(jù)牛頓第二定律，計(jì)算加速度，再由以上提供的算法即可進(jìn)行連續(xù)的模擬計(jì)算了。分子動力學(xué)模擬的系綜采用MD模擬，必

10、須再一定的系綜下進(jìn)行，經(jīng)常用到的系綜包括微正則系綜、正則系綜、等溫等壓系綜和等溫等焓系綜。微正則系綜（NVE）是孤立的、保守的系統(tǒng)的系綜，在這種系綜中，系統(tǒng)沿著相空間中的恒定能量軌道演化。在分子動力學(xué)模擬的過程中，系統(tǒng)中的原子數(shù)（N）體積（V）、和能量（E）都保持不變。正則系綜（NVT）NVT保持不變，并且總動量為零。恒溫下，系統(tǒng)的總能量不是一個守恒量，系統(tǒng)要與外界發(fā)生能量交換。保持系統(tǒng)的溫度不變，通常運(yùn)用的方法是讓系統(tǒng)與外界的熱浴處于熱平衡狀態(tài)。由于溫度與系統(tǒng)的動能有直接的關(guān)系，通常的做法是把系統(tǒng)的動能固定在一個給定值上。等溫等壓系綜（NPT）NPT保持不變，這種系綜是我們常見的系綜，許多分

11、子動力學(xué)模擬都要在從系綜下進(jìn)行。這時，要保證系統(tǒng)的溫度恒定，還要保持它的壓力恒定。溫度恒定和以前一樣，是通過調(diào)節(jié)系統(tǒng)的速度或加一約束力來實(shí)現(xiàn)的。而對壓力進(jìn)行調(diào)節(jié)，就比較復(fù)雜。由于系統(tǒng)的壓力P與其體積V是共軛量，要調(diào)節(jié)壓力值可以通過改變系統(tǒng)的體積來實(shí)現(xiàn)。等壓等焓系綜（NPH）NPH保持不變，焓值通過H=E+P得到。在系綜下進(jìn)行模擬時要保持壓力和焓值為固定值。邊界條件正確處理邊界和邊界效應(yīng)對模擬方法時至關(guān)重要的，因?yàn)樗鼤r從模擬相對較少的原子來計(jì)算物質(zhì)的宏觀性質(zhì)的。為了減小有限尺寸的影響，經(jīng)典的方法是使用周期性邊界條件。周期性邊界條件使得可以用相對少的粒子數(shù)目來真實(shí)的模擬大塊的體系，使得粒子仿佛處在

12、一個完整的體系中，假設(shè)一個用來模擬的立方單胞，使這個單胞在各個方向上都不斷重復(fù)，看起來象有周期一樣。在二維圖像中，每一個單胞都被其他的8個單胞所包圍；在三維方向上，每一個單胞就會被26個單胞所包圍。所有的單胞中的粒子的坐標(biāo)都可以通過一個整數(shù)而得到。當(dāng)模擬的單胞中的一個粒子由于力的作用而離開這個單胞的時候，九會有另一個和它對應(yīng)的粒子運(yùn)動到這個單胞中來，這樣，模擬的整個體系的粒子數(shù)就會保持不變。在分子動力學(xué)模擬中經(jīng)常采用的邊界條件有：矩形盒子周期性邊界條件單斜盒子周期性邊界條件去頭八面體合租周期性邊界條件八.MD模擬結(jié)果分析方法MD模擬的時間能夠達(dá)到幾百皮秒或納秒，甚至更長，在運(yùn)行MD模擬時,體系

13、的速度和坐標(biāo)被保存下來。在分析時可以對熱力學(xué)參數(shù)進(jìn)行計(jì)算。熱力學(xué)參數(shù)隨時間演化的特性可以用圖形顯示，一個坐標(biāo)對應(yīng)時間，另外一個坐標(biāo)時感興趣的物理量，如能量、均方根差、原子位置的漲落，此外還能計(jì)算出平均結(jié)構(gòu)與實(shí)驗(yàn)數(shù)據(jù)相對照，這些都有助于在原子水平上直觀地理解構(gòu)象地變化。1.均方根差（RMSD）一種廣泛用于檢驗(yàn)MD模擬正確性地方法時求相對于蛋白質(zhì)晶體結(jié)果地均方根差。目前已知一部分蛋白質(zhì)地晶體結(jié)構(gòu)與在溶液中的結(jié)果存在顯著的差異，但人們普通認(rèn)為這種差異對于大多數(shù)蛋白質(zhì)時非常小的。蛋白質(zhì)模擬中常有的質(zhì)量控制時獲得既小又穩(wěn)定的蛋白質(zhì)骨架原子的RMSD值（通常v0.2nm）。最常見的方法時求出對于參考結(jié)構(gòu)的

14、旋轉(zhuǎn)和平移擬合，原子i（如Ca原子）的集合的RMSD可按下式計(jì)算：其中ri和ri0分別指擬合后原子i的位置及其參考位置。N為要計(jì)算的原子總數(shù)。第二種方法無需進(jìn)行任何擬合，而是利用了距離矩陣：r嘶叫帀客國一硝）2其中，dij是原子i和j之間的距離，而dj.0是在參考結(jié)構(gòu)中原子i和j之間的距離。這種方法幣第一種好，更常用于nmR的研究。九.分子力學(xué)與分子力場在做分子動力學(xué)模擬時，常用到分子力場參數(shù)。分子力學(xué)的經(jīng)典力學(xué)模型:分子中的化學(xué)鍵具有“自然”鍵長、鍵角，并由這些鍵長和鍵角調(diào)節(jié)構(gòu)象，給出核位置的最佳分布，即分子的平衡構(gòu)型?；诮品肿恿W(xué)計(jì)算不顯含對電子的處理。并且，量子力學(xué)從頭算的計(jì)算量隨基

15、函數(shù)數(shù)目的四次方遞增，半經(jīng)驗(yàn)方法的計(jì)算量隨基函數(shù)數(shù)目的三次方遞增，而分子力學(xué)的計(jì)算量則僅與分子中原子數(shù)目的平方成正比。因此，分子力學(xué)己成當(dāng)前研究大分子體系的分子結(jié)構(gòu)、構(gòu)象平衡與轉(zhuǎn)變的研究。分子力場的概念從量子化學(xué)中知道，對于任何一個包括原子核和電子的微粒體系，若忽略自旋軌道及其它相關(guān)效應(yīng)，則可用定態(tài)方程來描述它的運(yùn)動其中稱之為哈密頓算符是體系的總能量，甲是核坐標(biāo)和電子坐標(biāo)的函數(shù)。多電子體系的定態(tài)方程不易求解，要借助簡化近似近似，簡稱近似。這個近似的依據(jù)是原子核比電子重得多這一事實(shí)，所以電子運(yùn)動比核快得多。在很好的近似下，當(dāng)電子運(yùn)動時，我們可以把核看成是固定的。于是，電子和核運(yùn)動就可以分開來處理

16、。根據(jù)近似，體系的總波函數(shù)可以寫成其中，甲是把核坐標(biāo)作為參變量的電子波函數(shù)沖是原子核的波函數(shù)。描述電子運(yùn)動的方程為在式中，定義了一個能量，它是體系中固定原子核坐標(biāo)時給定電子狀態(tài)下的電子本征能量函數(shù)，稱為分子勢能函數(shù)或分子內(nèi)勢能函數(shù)。在分子力場中，通常我們把稱為位能面，或稱面。如果一個解析表達(dá)式能擬合這個位能面，則此解析表達(dá)式就稱為分子力場，簡稱力場。即力場是:描述分子結(jié)構(gòu)和能量之間的一種數(shù)據(jù)解析表達(dá)式，包括原子核的坐標(biāo)和一些可調(diào)節(jié)的參數(shù)。所以，一個力場的確定就是選擇解析函數(shù)形式及確定參數(shù)2）力.場能的表示形式分子力學(xué)的基本理論就是一個分子力場由分子內(nèi)相互作用和分子間相互作用兩大部分構(gòu)成，分別對

17、應(yīng)于鍵伸縮、角彎曲、扭轉(zhuǎn)運(yùn)動、偶合相互作用以及相互作用和靜電相互作用或氫鍵能。力場的參數(shù)和力常數(shù)可從一個分子遷移至相似結(jié)構(gòu)的另一分子。一個力場中諸項(xiàng)的能量表達(dá)式基本上是相似的。其中，是各種成鍵相互作用能。成鍵相互作用能由鍵伸縮、角彎曲、二面角扭轉(zhuǎn)和鍵角面外彎曲能等子項(xiàng)組成。即是其中是鍵、角之間的偶合所相應(yīng)的能量在力場中稱為交叉項(xiàng)。這是因?yàn)榭紤]到分子內(nèi)鍵長及鍵角的變化受到相鄰鍵、角的影響。包括鍵一鍵偶合、鍵一角偶合等。其中是非鍵相互作用，由范德華能、靜電效應(yīng)及氫鍵能部分組成。具體到各種分子結(jié)構(gòu)和分子勢能函數(shù)的計(jì)算，由于研究的對象和考慮的因素不同而其解析表達(dá)式有所不同3）常.用的分子力場介紹第一代

18、經(jīng)典的分子力場。主要是，及F其特點(diǎn)是函數(shù)形式簡單;應(yīng)用的范圍比較特定（大部分適合于生物分子）;優(yōu)化力場參數(shù)的方法比較多，結(jié)果也比較好能合理地預(yù)測分子結(jié)構(gòu)，構(gòu)象性質(zhì)，凝聚態(tài)性質(zhì)。B第二代分子力場。及。這一類力場的共同特點(diǎn)是函數(shù)形式較復(fù)雜，附加項(xiàng)多、適用范圍寬，特別是大而復(fù)雜的分子模型的多能量極小和勢壘、優(yōu)化得到的力場參數(shù)比較合理;能比較好地預(yù)測分子結(jié)構(gòu)，振動頻率，構(gòu)象性質(zhì)。針對周期表中的所有元素的力場。所有的力學(xué)參數(shù)基于元素，雜化和化合的規(guī)則而產(chǎn)生，并為許多結(jié)構(gòu)形式所證實(shí)。這類力場有，F(xiàn)F。特點(diǎn)是:函數(shù)形式簡單;適用的范圍廣;能比較合理地預(yù)測分子結(jié)構(gòu);根照某些規(guī)則使力場參數(shù)化。力場是針對整個周期

19、表的分子力學(xué)和動態(tài)模擬的力場，力場參數(shù)是依靠元素、元素的雜化及化合性而計(jì)算出來的，與電荷平衡計(jì)算法相結(jié)合。用于特定目的的分子力場。如用于優(yōu)化玻璃質(zhì)的用于聚合體的，沸石結(jié)構(gòu)的吸附作用的等。這類的特點(diǎn)自然是只局限于某一特定的范圍內(nèi)，不能在其他場合中使用。十練習(xí)：真實(shí)水溶劑環(huán)境中的腺嘌吟一胸腺嘧啶十聚體的分子動力學(xué)模擬用.建立NUCDA5A命令創(chuàng)建1AA雙螺旋結(jié)構(gòu)A3文件AAAAANUC2DT5TTTTTTTTT3END$ABDNA運(yùn)行命令創(chuàng)建文件$AMBEROME/exe/nucgen-O-ia-dna.nucgen.in-oa-dna.nucgen.out-d$AMBERHOME/dat/lea

20、p/parm/nucgen.dat-pa-dna.nucgen.pdb這樣產(chǎn)生個文件。創(chuàng)建輸入文件下一步是為了給創(chuàng)建拓?fù)湮募妥鴺?biāo)文件。在用創(chuàng)建之前，我們需要給a-dna.nucgen.pdb文件中的第10、11個殘基之間加TER標(biāo)記。修改后的pdb文件為a-dna.pdb運(yùn)行：$AMBERHOME/exe/xleap-s-f$AMBERHOME/dat/leap/cmd/leaprc.ff99在xleap窗口創(chuàng)建一個工作單元adna打開a-dna.pdb文件：adna=loadpdba-dna.pdb給我們的系統(tǒng)加抗衡粒子（中性體系）：addionsadnaNa+0現(xiàn)在產(chǎn)生中性體系（非溶劑化

21、體系）的拓?fù)浜妥鴺?biāo)文件：saveamberparmadnaa-dna_charges_only.prmtopa-dna_charges_only.inpcrd下一步用一個截面八面體的水盒子使系統(tǒng)溶劑化：solvateoctadnaWATBOX2168.0#amber7.0不支持TIP3PBOX保存溶劑化后的拓?fù)浜妥鴺?biāo)文件：saveamberparmadnaa-dna_wat.prmtopa-dna_wat.inpcrd創(chuàng)）對系統(tǒng)進(jìn)行能量優(yōu)化：在對系統(tǒng)做之前要對系統(tǒng)做一次能量優(yōu)化，為了除去加氫和溶劑化時傳生的不良觸點(diǎn)。第一步：優(yōu)化離子和水的位置。優(yōu)化的前50步0用最陡下降法后50步0用共軛梯度a

22、-dna_min1.inA-DNA10-mer:initialminimisationsolvent+ions&cntrlimin=1,maxcyc=1000,TOC o 1-5 h zncyc=500,ntb=1,ntr=1,cut=10HoldtheDNAfixed500.0RES120ENDENDsander-O-ia-dna_min1.in-oa-dna_min1.out-pa-dna_wat.prmtop-ca-dna_wat.inpcrd-ra-dna_min1.rst-refa-dna_wat.inpcrd第二步：全局優(yōu)化。前100步0采用最陡下降法，后150步0采用共軛梯度法，

23、這時不需要約束力。輸入文件如下：a-dna_min2.inA-DNA10-mer:initialminimisationsolvent+ions&cntrlTOC o 1-5 h zimin=1,maxcyc=2500,ncyc=1000,ntb=1,ntr=0,cut=10sander-O-ia-dna_min2.in-oa-dna_min2.out-pa-dna_wat.prmtop-ca-dna_min1.rst-ra-dna_min2.rst初始平衡這一步開始運(yùn)行分子動力學(xué)。我們首先運(yùn)行的恒加熱平衡。采用的方的弱限制力，采用恒周期邊界條件做始溫度是，最終溫度是0采用溫度控制，在原子上采

24、用限制，步長為。adna_md1.inA-DNA10-mer:20psMDwithresonDNA&cntrlimin=0,irest=0,ntx=1,ntb=1,cut=10,ntr=1,ntc=2,ntf=2,tempi=0.0,temp0=300.0,ntt=3,gamma_ln=1.0,nstlim=10000,dt=0.002,ntpr=250,ntwx=250,ntwr=10000/KeepDNAfixedwithweakrestraints10.0RES120ENDENDsander-O-ia-dna_md1.in-oa-dna_md1.out-pa-dna_wat.prmtop

25、-ca-dna_min2.rst-ra-dna_md1.rst-refa-dna_min2.rst-xa-dna_md1.mdcrd進(jìn)一步的平衡這一步將運(yùn)行更長的的，對不加以限制，采用Langevin恒溫法使溫度保持在，采用恒壓周期邊界條件，對原子采用限制，步長為fadna_md_1800ps.inA-DNA10-mer:1800psofMD&cntrlimin=0,irest=1,ntx=7,ntb=2,pres0=1.0,ntp=1,taup=2.0,cut=10,ntr=0,ntc=2,ntf=2,tempi=300.0,temp0=300.0,ntt=3,gamma_ln=1.0,ns

26、tlim=900000,dt=0.002,ntpr=250,ntwx=250,ntwr=10000/$AMBERHOME/exe/sander-O-ia-dna_md_1800ps.in-oa-dna_md2.out-pa-dna_wat.prmtop-ca-dna_md1.rst-ra-dna_md2.rst-xa-dna_md2.mdcrd結(jié)果分析：_從輸出文件a-dna_md1.outa-dna_md2.out中可以看書系統(tǒng)的很多特性，如:系統(tǒng)的勢能、動能、總能量、溫度、壓力、體積、密度、均方根差。先用process-mdout.perl腳本對a-dna_md2.out進(jìn)行數(shù)據(jù)處理，產(chǎn)生一系列summary文件，再用繪圖程序xmgrace-5.99.0對上面產(chǎn)生的summary文件