分子動力學模擬及相關(guān)研究

分子動力學模擬及相關(guān)研究

MolecularDynamicsSimulationandRelatedResearch報告人：P.Li23-2月-23主要內(nèi)容1.粒子的運動服從Newton運動方程2.動力學機制決定運動方程的形式3.適用范圍廣泛1.分子運動方程及其數(shù)值求解2.多體系統(tǒng)的基本概念與分子動力學方法3.分子力場1.模擬模型的設(shè)定2.給定初始條件3.趨于平衡4.宏觀物理量的計算一、分子動力學特點二、分子動力學基礎(chǔ)知識三、分子動力學模擬步驟及相關(guān)研究一.特點1.在分子動力學中，粒子的運動行為是通過經(jīng)典的Newton運動方程所描述。

分子動力學是在原子、分子水平上求解多體問題的重要的計算機模擬方法，可以預測納米尺度上的材料動力學特性。

通過求解所有粒子的運動方程，分子動力學方法可以用于模擬與原子運動路徑相關(guān)的基本過程。

分子動力學模擬方法可以看作是體系在一段時間內(nèi)的發(fā)展過程的模擬。在這樣的處理過程中我們可以看出：分子動力學方法中不存在任何隨機因素。

在分子動力學方法處理過程中，方程組的建立是通過對物理體系的微觀數(shù)學描述給出的。在這個微觀的物理體系中，每個分子都各自服從經(jīng)典的牛頓力學。每個分子運動是用理論力學上的哈密頓量或者拉格朗日量來描述，也可以直接用牛頓運動方程來描述。2.系統(tǒng)的動力學機制決定運動方程的形式這種方法可以處理與時間有關(guān)的過程，因而可以處理非平衡態(tài)問題。但是使用該方法的程序較復雜，計算量大，占內(nèi)存也多。3.適用范圍廣泛：

原則上，分子動力學方法所適用的微觀物理體系并無什么限制。這個方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點粒子體系，也可以是具有內(nèi)部結(jié)構(gòu)的體系；處理的微觀客體既可以是分子，也可以是其它的微觀粒子。

自五十年代中期開始，分子動力學方法得到了廣泛的應(yīng)用。它與蒙特卡洛方法一起已經(jīng)成為計算機模擬的重要方法。應(yīng)用分子動力學方法取得了許多重要成果，例如氣體或液體的狀態(tài)方程、相變問題、吸附問題等，以及非平衡過程的研究。二.分子動力學基礎(chǔ)知識1、分子運動方程及其數(shù)值求解

采用分子動力學方法時，必須對一組分子運動微分方程做數(shù)值求解。從計算數(shù)學的角度來看，這是個求一個初值問題的微分方程的解。實際上計算數(shù)學為了求解這種問題已經(jīng)發(fā)展了許多的算法，但是并不是所有的這些算法都可以用來解決物理問題。

以一個一維諧振子為例，來看一下如何用計算機數(shù)值計算方法求解初值問題。一維諧振子的經(jīng)典哈密頓量為：這里的哈密頓量（即能量）為守恒量。假定初始條件為x（0）P（0），則它的哈密頓方程是對時間的一階微分方程

計算在相空間中的運動軌跡(x(t),p(t))：采用有限差分法，將微分方程變?yōu)橛邢薏罘址匠?，以便在計算機上做數(shù)值求解，并得到空間坐標和動量隨時間的演化關(guān)系。

首先，取差分計算的時間步長為h，采用一階微分形式的向前差商表示，直接運用展開到h的一階泰勒展開公式：即可得到：則微分方程改寫為差分形式：我們得到解微分方程的歐拉(Euler)算法：

這是一組遞推公式。有了初始條件，就可以一步一步地使用前一時刻的坐標、動量值確定下一時刻的坐標、動量值。這個方法是一步法的典型例子。

由于在實際數(shù)值計算時的大小是有限的，因而在上述算法中微分被離散化為差分形式來計算時總是有誤差的?？梢宰C明一步法的局部離散化誤差與總體誤差是相等的，都為的量級。在實際應(yīng)用中，適當?shù)剡x擇的大小是十分重要的。取得太大，得到的結(jié)果偏離也大，甚至于連能量都不守恒；取得太小，有可能結(jié)果仍然不夠好。這就要求我們改進計算方法，進一步考慮二步法。差分計算的二步法實際上泰勒展開式的一般形式為：前面敘述的歐拉算法就是取n=1。取n=2展開后可得到：這是二步法的一種,稱為Verlet方法。Verlet算法是分子動力學模擬中求解常微分方程最通用的方法2、多體系統(tǒng)的基本概念與分子動力學方法N體系統(tǒng)中，一個n體的密度函數(shù)一般可以寫為：其中是描寫系統(tǒng)的幾率函數(shù)。

通常為由系統(tǒng)中所有粒子的坐標、動量構(gòu)成的相空間中的任意一點。分子動力學元胞

分子動力學模擬方法往往用于研究大塊物質(zhì)在給定密度下的性質(zhì)，而實際計算模擬不可能在幾乎是無窮大的系統(tǒng)中進行。所以必須引進一個叫做分子動力學元胞的體積元,以維持一個恒定的密度。對氣體和液體，如果所占體積足夠大，并且系統(tǒng)處于熱平衡狀態(tài)的情況下，那么這個體積的形狀是無關(guān)緊要的。

對于晶態(tài)的系統(tǒng)，元胞的形狀是有影響的。為了計算簡便，對于氣體和液體，我們?nèi)∫粋€立方形的體積為分子動力學元胞。

設(shè)分子動力學元胞的線度大小為L，則其體積為L3。由于引進這樣的立方體箱子，將產(chǎn)生六個我們不希望出現(xiàn)的表面。模擬中碰撞這些箱的表面的粒子應(yīng)當被反射回到元胞內(nèi)部，特別是對粒子數(shù)目很少的系統(tǒng)。然而這些表面的存在對系統(tǒng)的任何一種性質(zhì)都會有重大的影響。元胞的周期性邊界條件

為了將分子動力學元胞有限立方體內(nèi)的模擬，擴展到真實大系統(tǒng)的模擬，我們通常采用周期性邊界條件。采用這種邊界條件,我們就可以消除引入元胞后的表面效應(yīng)，構(gòu)造出一個準無窮大的體積來更精確地代表宏觀系統(tǒng)。實際上，這里我們做了一個假定，即讓這個小體積元胞鑲嵌在一個無窮大的大塊物質(zhì)之中。周期性邊界條件的數(shù)學表示形式為：其中為任意的可觀測量，為任意整數(shù)。這個邊界條件就是命令基本分子動力學元胞完全等同地重復無窮多次。

該邊界條件的具體實現(xiàn)是這樣操作的：當有一個粒子穿過基本分子動力學元胞的六方體表面時，就讓這個粒子以相同的速度穿過此表面對面的表面重新進入分子動力學元胞內(nèi)。不同分子動力學元胞盒子內(nèi)粒子間的相互作用：

對于不同分子動力學元胞盒子內(nèi)粒子間的相互作用，如果相互作用是短程力，我們可以在長度rc處截斷。這里V(rc)必須要足夠小，以使截斷不會顯著地影響模擬結(jié)果。勢函數(shù)直接截斷

典型的分子動力學元胞尺度L通常選得比rc大很多。我們往往選擇元胞尺度滿足不等式條件L/2>rc，使得距離大于L/2的粒子的相互作用可以忽略，以避免有限尺寸效應(yīng)。通常L的數(shù)值應(yīng)當選得很大。

在考慮粒子間的相互作用時，通常采用最小像力約定。最小像力約定是在由無窮重復的分子動力學基本元胞中，每一個粒子只同它所在的基本元胞內(nèi)的另外N-1個中（設(shè)在此元胞內(nèi)有N個粒子）的每個粒子或其最鄰近的影像粒子發(fā)生相互作用。分子力學，又叫力場方法（forcefieldmethod），是基于經(jīng)典牛頓力學方程的一種計算分子的平衡結(jié)構(gòu)和能量的方法。與量子力學不同，它求解的是Newton方程，而不是薛定諤方程。基本假設(shè)：⑴.TheBorn-OppenheimerApproximation對勢能面的經(jīng)驗性擬合。量子力學中的薛定諤方程為:(非相對論和無時間依賴的情況下）

體系的哈密頓算符與原子核（R）和電子（r）位置相關(guān)的波函數(shù)分子力學----簡介3.分子力場

基于Born-Oppenheimer近似，其物理模型可描述為：原子核的質(zhì)量是電子質(zhì)量的103~105倍，電子速度遠遠大于原子核的運動速度，每當核的分布形式發(fā)生微小變更，電子立刻調(diào)整其運動狀態(tài)以適應(yīng)新的核場。

這意味著，在任一確定的核分布形式下，電子都有相應(yīng)的運動狀態(tài)；同時核間的相對運動可視為所有電子運動的平均結(jié)果。所以電子的波函數(shù)只依賴于原子核的位置，而不是他們的動能。于是這個近似認為，電子的運動與原子核的運動可以分開處理，可以將上式分解為電子運動的波函數(shù)核運動的波函數(shù)①②電子運動方程：核運動方程：基本假設(shè)方程①中的能量Ee(勢能面)僅僅是原子核坐標有關(guān)。相應(yīng)的，方程②所表示的為在核勢能面E(R)上的核運動方程。直接求解方程①，采用的是從頭算或者是半經(jīng)驗，這樣的量化計算都是把電子的波函數(shù)和能量處理成原子核坐標的函數(shù)。由于量子化學求解電子波函數(shù)和勢能面耗時巨大，常常將勢能面進行經(jīng)驗性的擬合，成為力場，由此構(gòu)成分子力學的基礎(chǔ)。將方程②用牛頓運動方程代替，勢能面采用力場擬合，就構(gòu)成了分子動力學的基礎(chǔ)?；炯僭O(shè)①②電子運動方程：核運動方程：⑵.簡單作用模型對體系相互作用的貢獻來自諸如鍵伸縮、鍵角的開合、單鍵的旋轉(zhuǎn)等等。即使使用類似Hooke定律這樣的簡單函數(shù)，也能令力場運轉(zhuǎn)良好。⑶.力場的可移植性（關(guān)鍵屬性）僅在少數(shù)情況下通過測試的一套函數(shù)，可以用來解決更廣范圍內(nèi)的問題。進一步講，從小分子得來的數(shù)據(jù)可以用來研究類似高分子的大分子。簡單分子力場分子力場是分子力學的核心。分子力學的基本理論就是一個分子力場由分子內(nèi)相互作用和分子間相互作用兩大部分構(gòu)成，即力場的勢能包括成鍵和非鍵相互作用，所有的勢能的總和即為分子的構(gòu)象能?；炯僭O(shè)力場分類力場分類規(guī)則定義：把鍵長、鍵角、二面角分成不同類型，同種類型的立場參數(shù)相同分類規(guī)則總是對原子進行。理論上說，如果一個力場有N個原子類型，那么就有N(N-1)/2種鍵，N(N-1)(N-2)/2種鍵角，N(N-1)(N-2)(N-3)/2種二面角舉例：乙烷分子1個C-C鍵，6個C-H鍵6個C-C-H鍵角，6個H-C-H鍵角9個H-C-C-H二面角傳統(tǒng)力場AMBER力場

Kollmangroup,1984最初僅為蛋白質(zhì)和核酸體系提供相應(yīng)的原子類型和力場參數(shù)1990，發(fā)展了適用于多糖模擬的力場參數(shù)（Homan1990）1995,2000，加入了適用于有機小分子的原子類型和參數(shù)傳統(tǒng)力場CHARMm力場（ChemistryatHarvardMacromolecularmechanics）Karplusgroup，1983適用于各種分子性質(zhì)的計算和模擬對于從孤立的小分子到溶劑化的大生物體系的多種模擬體系都可以給出較好的結(jié)果但不適合于有機金屬配合物傳統(tǒng)力場CVFF力場（ConsistentValenceForceField）Dauber-Osguthorpegroup,1988適用范圍包括有機小分子和蛋白質(zhì)體系擴展后可用于某些無機體系的模擬，如硅酸鹽、鋁硅酸鹽、磷鋁化合物主要用于預測分子的結(jié)構(gòu)和結(jié)合自由能傳統(tǒng)力場MMX力場Allingergroup,1989包括MM2和MM3主要針對有機小分子函數(shù)形式比較復雜，包含交叉項也可用于生物大分子體系，但是速度會比較慢傳統(tǒng)力場第二代力場CFF力場CFF91主要適用于模擬有機小分子、蛋白質(zhì)以及小分子-蛋白質(zhì)之間的相互作用CFF95除了適用于蛋白質(zhì)和有機小分子體系，還可用于有機高分子體系的模擬，如聚碳酸酯及多糖PCFF在CFF91的基礎(chǔ)上，還適用于聚碳酸酯、三聚氰胺甲醛樹脂、多糖、核酸、分子篩等其他無機和有機材料體系的模擬。同時還提供20種金屬離子的參數(shù)第二代力場COMPASS力場（Condensed-phaseOptimizedMolecularPotentialsforAtomisticSimulationStudies用于原子水平模擬研究的凝聚態(tài)優(yōu)化的分子力場）Sungroup,1994第一個能同時預測氣態(tài)和凝聚態(tài)性質(zhì)的從頭力場適用于常見的有機小分子、無機小分子和高分子最新的力場也適合于金屬、金屬氧化物及金屬鹵化物的模擬可較為準確的預測晶體的多種性質(zhì)第二代力場MMFF94力場（MerkMolecularForceField）Hagler,1996目前最為準確的力場之一定義了非常完備的原子類型既適用于有機小分子，也適用于大分子體系，如蛋白質(zhì)第二代力場通用力場基于規(guī)則的力場（Rule-basedForceFields）較好的可移植性，適用范圍廣精度差ESFF力場（ExtensiveSystematicForceField）用于有機分子、無機分子及有機金屬化合物的結(jié)構(gòu)預測UFF力場（UniversalForceField）覆蓋了周期表中所有元素，應(yīng)用最廣泛的通用力場Dreiding力場適用于有機分子、生物大分子及位于主族的無機分子通用力場分子力場的選擇蛋白質(zhì)分子的模擬：首選AMBER力場和CHARMm力場，也可用CFF力場、CVFF力場和MMFF94力場核酸分子的模擬：采用AMBER力場、CHARMm力場、MMFF94力場或用戶自定義的力場小分子-蛋白質(zhì)復合物體系的模擬：首選CHARMm力場和MMFF94力場，也可用CVFF力場和CFF力場高分子的模擬：首選COMPASS力場，也可用PCFF力場和CFF95力場分子力場的選擇三.分子動力學模擬的基本步驟在計算機上對分子系統(tǒng)的分子動力學模擬的實際步驟可以劃分為四步：首先是設(shè)定模擬所采用的模型；第二，給定初始條件；第三，趨于平衡的計算過程；最后是宏觀物理量的計算。1．模擬模型的設(shè)定設(shè)定模型是分子動力學模擬的第一步工作。例如在一個分子系統(tǒng)中，假定兩個分子間的相互作用勢為硬球勢，實際上，更常用的是Lennard-Jones型勢。它的勢函數(shù)表示為:其中ε?是位勢的最小值，r=σ時位勢為零.第i個粒子與元胞內(nèi)其它N-1個粒子或其最鄰近的影像粒子的相互作用力在x方向的分量為:Y和Z分量的表達式與上式相同模型確定后,根據(jù)經(jīng)典物理學的規(guī)律我們就可以知道在系綜模擬中的守恒量。例如對在微正則系綜的模擬中能量、動量和角動量均為守恒量。為了計算方便，取分子動力學元胞為立方形，元胞的線度大小為L。例如，Verlet方法需要兩組坐標來啟動計算：一組是零時刻的坐標，另一組是前進一個時間步長時的坐標，或者是一組零時刻的速度值。2．給定初始條件分子動力學模擬的過程進入對系統(tǒng)微分方程組做數(shù)值求解時，需要知道粒子的初始位置和速度的數(shù)值。不同的算法要求不同的初始條件。

但是，一般來說系統(tǒng)的初始條件都是不可能知道的。表面上看這是一個難題。實際上，精確選擇待求系統(tǒng)的初始條件是沒有什么意義的，因為模擬時間足夠長時，系統(tǒng)就會忘掉初始條件。但是初始條件的合理選擇將可以加快系統(tǒng)趨于平衡。常用的初始條件可以選擇為：（1）令初始位置在差分劃分網(wǎng)格的格子上，初始速度則從玻爾茲曼分布隨機抽樣得到。（2）令初始位置隨機地偏離差分劃分網(wǎng)格的格子，初始速度為零。3．趨于平衡

按照上面給出的運動方程、邊界條件和初始條件，就可以進行分子動力學模擬計算。

但是，這樣計算出的系統(tǒng)不會具有所要求的系統(tǒng)能量，并且這個狀態(tài)本身也還不是一個平衡態(tài)。

為了使系統(tǒng)達到平衡，模擬中需要一個趨衡過程。

在這個過程中，我們增加或從系統(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。

然后，再對運動方程中的時間向前積分若干步，使系統(tǒng)持續(xù)給出確定能量值。我們稱這時系統(tǒng)已經(jīng)達到平衡態(tài)。這段達到平衡所需的時間稱為弛豫時間。

在分子動力學模擬中，時間步長h的大小選擇是十分重要的。

為了減小誤差，步長必須取得小一些；但是取得太小，系統(tǒng)模擬的弛豫時間就很長。這里需要積累一定的模擬經(jīng)驗，選擇適當?shù)臅r間步長h。4．宏觀物理量的計算

實際計算宏觀物理量往往是在分子動力學模擬的最后階段進行的。它是沿著相空間軌跡求平均來計算得到的。例如對于一個宏觀物理量A，它的測量值為平均值如果已知初始位置和動量為和(上標N表示系綜N個粒子的對應(yīng)坐標和動量參數(shù)），選擇某種分子動力學算法求解具有初值問題的運動方程，便得到相空間軌跡:對軌跡平均的宏觀物理量A的表示為:如果宏觀物理量為動能，它的平均為:分子動力學相關(guān)軟件CHARMm:由哈佛大學開發(fā)的分子力學計算工具?？梢詾橛脩籼峁┨幚砀鞣N小分子、大分子(包括蛋白質(zhì)、核酸和糖)的經(jīng)驗化能量計算,模擬過程提供了有關(guān)分子結(jié)構(gòu)、相互作用、能量等信息。Amber：用于分子（特別是生物大分子）的動力學模擬以及小分子與大分子相互作用的研究工具ABINIT：主程序使用贗勢和平面波，用密度泛函理論計算總能量，電荷密度，分子和周期性固體的電子結(jié)構(gòu)，進行幾何優(yōu)化和分子動力學模擬，用TDDFT（對分子）或GW近似（多體微擾理論）計算激發(fā)態(tài)。此外還提供了大量的工具程序。程序的基組庫包括了元素周期表1-109號所有元素。ABINIT適于固體物理，材料科學，化學和材料工程的研究，包括固體，分子，材料的表面，以及界面，如導體、半導