分子動(dòng)力學(xué)模擬及相關(guān)研究

分子動(dòng)力學(xué)模擬及相關(guān)研究

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

分子動(dòng)力學(xué)是在原子、分子水平上求解多體問題的重要的計(jì)算機(jī)模擬方法，可以預(yù)測納米尺度上的材料動(dòng)力學(xué)特性。

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

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

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

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

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

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

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

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

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

這是一組遞推公式。有了初始條件，就可以一步一步地使用前一時(shí)刻的坐標(biāo)、動(dòng)量值確定下一時(shí)刻的坐標(biāo)、動(dòng)量值。這個(gè)方法是一步法的典型例子。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

按照上面給出的運(yùn)動(dòng)方程、邊界條件和初始條件，就可以進(jìn)行分子動(dòng)力學(xué)模擬計(jì)算。

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

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

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

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

在分子動(dòng)力學(xué)模擬中，時(shí)間步長h的大小選擇是十分重要的。

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

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