1、關于分子力學和分子動力學現在學習的是第一頁，共34頁學習目標學習目標v 了解計算化學的主要方法和基本概念了解計算化學的主要方法和基本概念v 能夠讀懂相關領域的文獻，了解常見分子模擬軟件及使用能夠讀懂相關領域的文獻，了解常見分子模擬軟件及使用v 掌握基本的計算化學研究思路和方法。掌握基本的計算化學研究思路和方法。v 了解分子模擬能干什么，我們可以作那些工作，作到什么精了解分子模擬能干什么，我們可以作那些工作，作到什么精度，為在自己的研究領域開展分子模擬方面的研究打一個基度，為在自己的研究領域開展分子模擬方面的研究打一個基礎。礎?，F在學習的是第二頁，共34頁分子模擬具有理論和實驗的雙重性質分子模擬

2、不能完全取代實驗理論實驗模擬理論的正確性模擬參數的正確性模擬方法的選擇理論的更新現在學習的是第三頁，共34頁系綜系綜 系綜（ensemble）代表一大群相類似的體系的集合。對一類相同性質的體系，其微觀狀態微觀狀態（比如每個粒子的位置和速度）仍然可以大不相同。（實際上，對于一個宏觀體系，所有可能的微觀狀態數是天文數字。） 統計物理的一個基本假設（各態歷經各態歷經假設）是：對于一個處于平衡的體系，物理量的時間平均，等于對對應系綜里所有體系進行平均的結果。 體系的平衡態的物理性質可以對不同的微觀狀態求和來得到。系綜的概念是由約西亞威拉德吉布斯（J. Willard Gibbs）在1878年提出的。

3、現在學習的是第四頁，共34頁常用系綜常用系綜v微正則系綜微正則系綜 (microcanonical ensemble)v正則系綜正則系綜 (canonical ensemble) v巨正則系綜巨正則系綜 (grand canonical ensemble)v等溫等壓系綜等溫等壓系綜 (isothermal-isobaric ensemble)現在學習的是第五頁，共34頁常用系綜常用系綜微正則系綜微正則系綜 (microcanonical ensemble)：系綜里的每個體系具有相同的能量（通常每個體：系綜里的每個體系具有相同的能量（通常每個體系的粒子數和體積也是相同的）。系的粒子數和體積也是相

4、同的）。正則系綜正則系綜 (canonical ensemble)：系綜里的每個體系都可以和其他體系交換能量（每個：系綜里的每個體系都可以和其他體系交換能量（每個體系的粒子數和體積仍然是固定且相同的），但是系綜里所有體系的能量總和是固定體系的粒子數和體積仍然是固定且相同的），但是系綜里所有體系的能量總和是固定的。系綜內各體系有相同的溫度。的。系綜內各體系有相同的溫度。巨正則系綜巨正則系綜 (grand canonical ensemble)：正則系綜的推廣，每個體系都可以和：正則系綜的推廣，每個體系都可以和其他體系交換能量和粒子，但系綜內各體系的能量總和以及粒子數總和都其他體系交換能量和粒子，

5、但系綜內各體系的能量總和以及粒子數總和都是固定的。（系綜內各體系的體積相同。）系綜內各個體系有相同的溫度是固定的。（系綜內各體系的體積相同。）系綜內各個體系有相同的溫度和化學勢。和化學勢。等溫等壓系綜等溫等壓系綜 (isothermal-isobaric ensemble)：正則系綜的推廣，體系間可交換：正則系綜的推廣，體系間可交換能量和體積，但能量總和以及體積總和都是固定的。（系綜內各體系有相能量和體積，但能量總和以及體積總和都是固定的。（系綜內各體系有相同的粒子數。）正如它的名字，系綜內各個體系有相同的溫度和壓強。同的粒子數。）正如它的名字，系綜內各個體系有相同的溫度和壓強。 現在學習的是

6、第六頁，共34頁分子動力學方法工作框圖分子動力學方法工作框圖給定t 時刻的坐標和速度以及其他動力學信息，那么就可計算出 t t 時刻的坐標和速度。現在學習的是第七頁，共34頁程序構成方式 輸入指定運算條件的參數（初始溫度，粒子數，密度，時間步長） ； 體系初始化（選定初始坐標和初始速度） ； 計算作用在所有粒子上的力； 解牛頓運動方程（第 和第 步構成了模擬的核心，重復這兩步，直到體系的演化到指定的時間） ； 計算并輸出物理量的平均值，完成模擬。現在學習的是第八頁，共34頁積分算法優劣的判據積分算法優劣的判據 分子動力學中一個好的積分算法的判據主要包括： 計算速度快； 需要較小的計算機內存；

7、允許使用較長的時間步長； 表現出較好的能量守恒?，F在學習的是第九頁，共34頁分子動力學的適用范圍分子動力學方法只考慮多體系統中原子核的運動，而電子的運動不予考慮，量子效應忽略。經典近似在很寬的材料體系都較精確；但對于涉及電荷重新分布的化學反應、鍵的形成與斷裂、解離、極化以及金屬離子的化學鍵都不適用，此時需要使用量子力學方法。經典分子動力學方法（MD）也不適用于低溫，因為量子物理給出的離散能級之間的能隙比體系的熱能大，體系被限制在一個或幾個低能態中。當溫度升高或與運動相關的頻率降低（有較長的時間標度）時，離散的能級描述變得不重要，在這樣的條件下，更高的能級可以用熱激發描述?，F在學習的是第十頁，共

8、34頁 電子運動具有更高的特征頻率，必須用量子力學以及量子經典理論聯合處理。這些技術近年來取得了很大進步。在這些方法中，體系中化學反應部分用量子理論處理，而其他部分用經典模型處理。QM/MM現在學習的是第十一頁，共34頁分子力學、分子動力學方法及其應用分子力學、分子動力學方法及其應用Molecular Mechanics / Molecular Dynamics 一、 MM、MD理論基礎 二、 MM、MD計算程序 三、 MM、MD方法的應用現在學習的是第十二頁，共34頁一、一、 MM、MD基礎理論基礎理論原則上, 第一原理方法在理論上已經能解決所有問題但計算量太大，計算機資源有限，原子數目較多

9、時，如高分子、蛋白質、原子簇以及研究表面問題、功能材料或材料的力學性能等，實際上難以完成計算為此，發展了分子力學(Molecular Mechanics, MM)與分子動力學(Molecular Dynamics, MD)方法它們的應用，又稱分子模擬(molecular simulation, molecular modeling) 或 分子設計(molecular design)MM與MD是經典力學方法，針對的最小結構單元不再是電子電子而是原子 因原子的質量比電子大很多，量子效應不明顯，可近似用經典力學方法處理20 世紀 30 年代, Andrews 最早提出分子力學(MM)的基本思想；40

10、 年代以后得到發展, 并用于有機小分子研究。90年代以來得到迅猛發展和廣泛應用現在學習的是第十三頁，共34頁基本思想基本思想 事先構造出簡單體系(如鏈段、官能團等各種不同結構的小片段)的勢能函數, 簡稱 勢函數 或 力場(force field) 將勢函數建成數據庫，在形成較大分子的勢函數時，從數據庫中檢索到結構相同的片段，組合成大體系的勢函數 利用分子勢能隨原子位置的變化有極小值的性質，確定大分子的結構即為分子力學(MM) 利用勢函數，建立并求解與溫度和時間有關的牛頓運動方程，得到一定條件下體系的結構隨時間的演化關系即為分子動力學(MD) 理論方法的核心核心是構造勢函數 勢函數：勢能與原子位

11、置的關系。且往往是不知道的 需要通過其他方法，如量子化學方法及實驗數據獲得現在學習的是第十四頁，共34頁rrrE過渡狀態 E反應物產物分子勢函數曲面勢能面示意圖現在學習的是第十五頁，共34頁 1、分子力場、分子力場 分子片段力場的函數表達式中包含自變量和力場參數 其中自變量為分子的結構參數，獨立參數為鍵長、鍵角和二面角，如圖 而 力場參數 一般通過與實驗數據 和 從頭算數據進行最小二乘法擬合來確定bbb有的還使用一個非獨立參數：面外彎曲角現在學習的是第十六頁，共34頁 勢函數形式很多，目前已被廣泛使用的力場有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING

12、、UFF和COMPASS等 形式雖多, 但一般總表達為分子內與分子間勢能之和:V總總=V鍵合鍵合+V非鍵合非鍵合 分子內勢能(鍵合)包括鍵伸縮、鍵角彎曲和二面角扭轉勢能 分子間勢能(非鍵合)包括范德華勢和靜電勢, 有的還包括H鍵：V鍵合=V鍵伸縮+V鍵角彎曲+V二面角扭轉V非鍵合=V范德華+V靜電+V氫鍵 鍵合勢函數中，一些力場還包含交叉項，使精度更高 交叉項的含義：如鍵長變化時，鍵角彎曲勢能隨鍵長的不同而不同，等現在學習的是第十七頁，共34頁 例例: COMPASS-98力場(condensed-phase optimized molecular potentials for atomist

13、ic simulation studies)的表達式如下每個k是一獨立的力場參數,下標“0”代表參考(平衡)結構參數：bbbbkbbkbbkE)()()(404303202鍵伸縮：)()()(404303202kkkE鍵彎曲：)3cos(1 )2cos(1 )cos(1 3 , 032, 021 , 01kkkE二面角：2kE鍵面外彎曲：,0000)()(bbbbxkbbbbkE交叉項：,3210003cos2coscos)()(bbbkkkbbbbk, 32103cos2coscos)(bkkkbb,32103cos2coscos)(kkk,00)(cosk現在學習的是第十八頁，共34頁 力

14、場參數力場參數k最小二乘法確定最小二乘法確定 基本思想 如：R-COOH基團 1) 由ab initio (構型優化方法)計算出 平衡結構，得到 參考結構參數bi0, i0, i0 2) 用偽隨機數方法將bi0, i0, i0人為改變成若干(n)組非平衡結構參數 bi, i, il (l=1,2, , n) 3) 用各bi, i, il 結構參數分別進行ab initio計算，得到bi, i, il 結構參數下對應的能量El (l=1,2, , n) 4) 將El (l=1,2, , n)和bi, i, il代入勢能表達式 5) 用最小二乘法擬合，確定力場參數k現在學習的是第十九頁，共34頁

15、非鍵合勢函數中，靜電相互作用表示分子中各原子靜電荷的庫侖相互作用對勢能的貢獻 不同的力場，靜電相互作用表達式基本相同 范德華勢也大都采用Lennard-Jones函數，但函數中的指數有所不同。如COMPASS-98的非鍵合勢函數為：靜電相互作用：jiijjirqqEelec范德華勢：jiijijijijijvdWrrrrE,609032即Lennard-Jones 9-6函數其它力場范德華勢較多采用L-J 12-6函數：jiijijijijijvdWrrrrDE,6012002)(現在學習的是第二十頁，共34頁 當然，在建立分子的勢能函數時，還有一些更細致的問題要考慮，如： 勢能展開項的截斷

16、周期結構的處理 多組分混合物體系 含有離子的體系 或 金屬中的離子 等 不同的方法或程序中，分別都有更詳細的討論。參考：德 D. 羅伯. 計算材料學. 北京：化學工業出版社, 2002, 9 俞慶森,朱龍觀. 分子設計導論.北京：高等教育出版社,2000楊小震. 分子模擬與高分子材料. 北京：科學出版社，2002熊家炯主編. 材料設計. 天津：天津大學出版社，2000Sun H，Ren P, Fried J R. The COMPASS Force Field: Parameterization and Validation for polyphosphazenes. Computationa

17、l and Theoretical Polymer Science, 1998, 8(1/2): 229Sun H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem., 1998, 102: 7338 現在學習的是第二十一頁，共34頁 2、分子力學方法、分子力學方法 MM是確定分子結構的方法 利用分子勢能隨結構的變化而變化的性質，確定分子勢能極小時的

18、平衡結構(stationary point) 物理模型: 視原子為質點，視化學鍵為彈簧，而彈力常數完全由數據庫中的分子力場來確定 因此是直接用勢函數研究問題，不考慮原子的動能 不考慮動能所對應的結構，相當于體系處于 T=0 K 時的結果現在學習的是第二十二頁，共34頁 由力場首先構造并得到分子的勢函數E(x) 利用在E(x)的極小點處，E(x)隨各原子獨立的空間坐標xi, i=1,2,3,3N-6 (與內坐標ri、i、i等價) 的一階微分=0以及全部二階微分0的數學條件：0)(0)(22iixxExxE(i=1,2,3,3N-6) 進行結構優化 具體步驟： 1) 用各種方法構造出一個分子的任意

19、結構，得到初始的結構參數xi0, i=1,2,3,3N-6 2) 進行坐標變換，即根據需要將原子的直角坐標轉變成內坐標ri、i、i或反之現在學習的是第二十三頁，共34頁1) 用各種方法構造出一個分子的任意結構，得到初始的結構參數xi0, i=1,2,3,3N-62) 進行坐標變換，即根據需要將原子的直角坐標轉變成內坐標(鍵長、鍵角、二面角)或反之 3) 建立分子體系的勢能表達式E(x) 4) 計算E(x)隨各坐標的一階、二階導數 5) 計算接近數學條件的坐標增量 6) 得到新的結構參數xi1, i=1,2,3,3N-6 重復 4)、5)、6)，直至最后兩次得到的體系勢能之差或總體勢能梯度的均方

20、根值達到預定精度范圍為止 MM計算中，有時還要考慮體系所處的外壓條件，此時，可通過壓力因子的定義來調節原子坐標0)( , 0)(22iixxExxE(i=1,2,3,3N-6)現在學習的是第二十四頁，共34頁 3、分子動力學方法、分子動力學方法 分子力場是分子的靜態勢函數 而實際過程通常是在一定溫度和一定壓力下發生的 為了更切實際地了解體系運動和演化的過程，必須考慮體系中原子的運動，并與溫度T和時間t建立聯系 我們知道，溫度是原子分子熱運動劇烈程度的量度 根據統計熱力學，對于n個原子的體系，體系的溫度T與各原子的運動速率vi的關系為：niiivmnTk12B31 又因體系中各原子的速率為vi時

21、，動量pi=mivi，對應總動能K(p)為：niziyixiipppmpK12,2,2,)(21)(現在學習的是第二十五頁，共34頁 勢能由力場確定為E(x)，因此體系的Hamilton量H為：)()(),(xEpKpxH 與量子力學不同，經典力學對Hamilton量不進行算符化處理，也不建立和求解本征方程，而是建立并求解經典運動方程：ppxHtxxpxHtp),(dd ),(dd 計算過程一般為：在一定的統計系綜下 1) 由原子位置和連接方式，從數據庫調用力場參數并形成體系勢函數 2) 由給定溫度計算體系動能以及總能量 3) 計算各原子的勢能梯度, 得到原子在力場中所受的力 即dp/dt=m

22、dv/dt=ma=F現在學習的是第二十六頁，共34頁 1) 由原子位置和連接方式，調用力場參數并形成體系勢函數2) 由給定溫度計算體系動能以及總能量3) 計算各原子的勢能梯度, 得到原子在力場中所受的力 即dp/dt=mdv/dt=ma=F 4) 對每個原子，在一定時間間隔內，用牛頓方程求解其運動行為：ppxHtxxpxHtp),(dd ,),(dd 22ddddtxmtvmamFiiiiiii 5) 顯示體系能量和結構 6) 取下一時間間隔，返回步驟 1) 不斷循環反復，可設定循環次數或強行終止計算 其中, 不僅時間間隔可以根據需要取不同大小(一般1 fs=10-15 s),溫度可以任意設定

23、，而且還可以在循環過程中逐漸改變溫度，即研究體系的退火(annealing)行為現在學習的是第二十七頁，共34頁 有時需要進一步考慮外場的作用，如壓力、電場、磁場、重力場等 從原則上講，這些問題都不難解決，而且還在進一步發展中，不再贅述此外，即使是分子力學與分子動力學方法，也受計算量的限制，所處理的體系不可能太大，好的計算機可達到數十萬百萬個原子的規模 (2050 nm以下/ MD 106108個原子羅伯p110) 對于更大(50 nm、微米)尺度的問題，人們也在發展介觀尺度(mesoscale)的方法，并且已取得一些成果 分子動力學方法也有一些本質的缺陷 如：勢函數精度(客觀性) 勢函數形式

24、在每次計算中都不變, 故不能模擬如分子在高溫下, 結構發生斷裂的熱裂解過程 (國內引入開關函數已開始做) 此時, 可用QM處理化學變化或e轉移, 余用MM，稱QM/MM現在學習的是第二十八頁，共34頁二、二、 MM、MD計算程序計算程序 因程序結構相對簡單，國內外有不少, 且還在不斷產生如：Chem3D, ChemOffice, Cerius2, Material Studio, Alchemy 2000, Sybyl, Biosym, HyperChem x.x, Spartan x.x, Chemgraf, Bilder, Script, COGS, Gaussian-03; MP(Mol

25、ecular Properties, 楊小震)等 輸入：輸入：體系模型建模(Builder)，選擇力場、系綜 MM: 無特殊輸入 MD: 循環次數(萬幾10萬)、溫度、時間間隔(fs), 外場,溫度變化速率，其它性質 等 輸出：輸出： MM: 勢能-結構曲線(數值)，動態結構圖, 其它性質 MD: (同上)，及 其它MD性質注：注：第一原理及MM一般只能得到勢能極小點結構，而MD可越過一些小的勢壘，甚至可達到最小點現在學習的是第二十九頁，共34頁三、三、 MM、MD方法的應用方法的應用 領域：高分子、生命科學、藥物設計、催化、半導體 其它功能材料、結構材料等 分子力學是用計算機在原子水平上模擬給定分子模型的結構與性質，進而得到分子的各種物理性質與化學性質，如結構參數、振動頻率、構象能量、相互作用能量、偶極矩、密度、摩爾體積、汽化焓等 分子動力學方法能實時將分子的動態行為顯示到計算機屏幕上, 便于直觀了解體系在一定條件下的演變過程 MD含溫度與時間, 因此還可得到如材料的玻璃化轉變溫度、熱容、晶體結晶過程、輸送過程、膨脹過程、動態弛豫(relax)以及體系在外場作用下的變化過程等現在學習的是第三十頁，共