1、計(jì)算材料科學(xué)中的模型、算法和多尺度關(guān)聯(lián) 系列講座 與 討論 西北工業(yè)大學(xué)科研處 西北工業(yè)大學(xué)理學(xué)院 西北工業(yè)大學(xué)材料學(xué)院 聯(lián)合舉辦 西北工業(yè)大學(xué)高性能計(jì)算中心,III、納米、亞納米尺度 一、第一原理方法及其在材料科學(xué)中的應(yīng)用 二、密度泛函方法及其應(yīng)用 三、從原子分子到納米尺度 分子力學(xué)、分子動(dòng)力學(xué)方法及其應(yīng)用,從原子分子到納米尺度 分子力學(xué)、分子動(dòng)力學(xué)方法及其應(yīng)用 Molecular Mechanics / Molecular Dynamics 一、 MM、MD理論基礎(chǔ) 二、 MM、MD計(jì)算程序 三、 MM、MD方法的應(yīng)用,一、 MM、MD基礎(chǔ)理論 原則上, 第一原理方法在理論上已經(jīng)能解決所有

2、問題 但計(jì)算量太大，計(jì)算機(jī)資源有限，原子數(shù)目較多時(shí)，如高分子、蛋白質(zhì)、原子簇以及研究表面問題、功能材料或材料的力學(xué)性能等，實(shí)際上難以完成計(jì)算 為此，發(fā)展了分子力學(xué)(Molecular Mechanics, MM)與分子動(dòng)力學(xué)(Molecular Dynamics, MD)方法 它們的應(yīng)用，又稱分子模擬(molecular simulation, molecular modeling) 或 分子設(shè)計(jì)(molecular design) MM與MD是經(jīng)典力學(xué)方法，針對(duì)的最小結(jié)構(gòu)單元不再是電子而是原子 因原子的質(zhì)量比電子大很多，量子效應(yīng)不明顯，可近似用經(jīng)典力學(xué)方法處理 20 世紀(jì) 30 年代, An

3、drews 最早提出分子力學(xué)(MM)的基本思想；40 年代以后得到發(fā)展, 并用于有機(jī)小分子研究。90年代以來得到迅猛發(fā)展和廣泛應(yīng)用,基本思想 事先構(gòu)造出簡單體系(如鏈段、官能團(tuán)等各種不同結(jié)構(gòu)的小片段)的勢(shì)能函數(shù), 簡稱 勢(shì)函數(shù) 或 力場(force field) 將勢(shì)函數(shù)建成數(shù)據(jù)庫，在形成較大分子的勢(shì)函數(shù)時(shí)，從數(shù)據(jù)庫中檢索到結(jié)構(gòu)相同的片段，組合成大體系的勢(shì)函數(shù) 利用分子勢(shì)能隨原子位置的變化有極小值的性質(zhì)，確定大分子的結(jié)構(gòu)即為分子力學(xué)(MM) 利用勢(shì)函數(shù)，建立并求解與溫度和時(shí)間有關(guān)的牛頓運(yùn)動(dòng)方程，得到一定條件下體系的結(jié)構(gòu)隨時(shí)間的演化關(guān)系即為分子動(dòng)力學(xué)(MD) 理論方法的核心是構(gòu)造勢(shì)函數(shù) 勢(shì)函數(shù)：

4、勢(shì)能與原子位置的關(guān)系。且往往是不知道的 需要通過其他方法，如量子化學(xué)方法及實(shí)驗(yàn)數(shù)據(jù)獲得,分子勢(shì)函數(shù)曲面勢(shì)能面示意圖,1、分子力場 分子片段力場的函數(shù)表達(dá)式中包含自變量和力場參數(shù) 其中自變量為分子的結(jié)構(gòu)參數(shù)，獨(dú)立參數(shù)為鍵長、鍵角和二面角，如圖,而 力場參數(shù) 一般通過與實(shí)驗(yàn)數(shù)據(jù) 和 從頭算數(shù)據(jù)進(jìn)行最小二乘法擬合來確定,勢(shì)函數(shù)形式很多，目前已被廣泛使用的力場有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING、UFF和COMPASS等 形式雖多, 但一般總表達(dá)為分子內(nèi)與分子間勢(shì)能之和: V總=V鍵合+V非鍵合 分子內(nèi)勢(shì)能(鍵合)包括鍵伸縮、鍵角彎曲和二面角扭轉(zhuǎn)勢(shì)能

5、 分子間勢(shì)能(非鍵合)包括范德華勢(shì)和靜電勢(shì), 有的還包括H鍵： V鍵合=V鍵伸縮+V鍵角彎曲+V二面角扭轉(zhuǎn) V非鍵合=V范德華+V靜電+V氫鍵 鍵合勢(shì)函數(shù)中，一些力場還包含交叉項(xiàng)，使精度更高 交叉項(xiàng)的含義：如鍵長變化時(shí)，鍵角彎曲勢(shì)能隨鍵長的不同而不同，等,例: COMPASS-98力場(condensed-phase optimized molecular potentials for atomistic simulation studies)的表達(dá)式如下 每個(gè)k是一獨(dú)立的力場參數(shù),下標(biāo)“0”代表參考(平衡)結(jié)構(gòu)參數(shù)：,顯然是一個(gè)非諧性力場,力場參數(shù)k最小二乘法確定 基本思想 如：R-COOH

6、基團(tuán) 1) 由ab initio (構(gòu)型優(yōu)化方法)計(jì)算出 平衡結(jié)構(gòu)，得到 參考結(jié)構(gòu)參數(shù)bi0, i0, i0 2) 用偽隨機(jī)數(shù)方法將bi0, i0, i0人為改變成若干(n)組 非平衡結(jié)構(gòu)參數(shù) bi, i, il (l=1,2, , n) 3) 用各bi, i, il 結(jié)構(gòu)參數(shù)分別進(jìn)行ab initio計(jì)算，得到bi, i, il 結(jié)構(gòu)參數(shù)下對(duì)應(yīng)的能量El (l=1,2, , n) 4) 將El (l=1,2, , n)和bi, i, il代入勢(shì)能表達(dá)式 5) 用最小二乘法擬合，確定力場參數(shù)k,非鍵合勢(shì)函數(shù)中，靜電相互作用表示分子中各原子靜電荷的庫侖相互作用對(duì)勢(shì)能的貢獻(xiàn) 不同的力場，靜電相互作

7、用表達(dá)式基本相同 范德華勢(shì)也大都采用Lennard-Jones函數(shù)，但函數(shù)中的指數(shù)有所不同。如COMPASS-98的非鍵合勢(shì)函數(shù)為：,靜電相互作用：,其它力場范德華勢(shì)較多采用L-J 12-6函數(shù)：,當(dāng)然，在建立分子的勢(shì)能函數(shù)時(shí)，還有一些更細(xì)致的問題要考慮，如： 勢(shì)能展開項(xiàng)的截?cái)?周期結(jié)構(gòu)的處理 多組分混合物體系 含有離子的體系 或 金屬中的離子 等 不同的方法或程序中，分別都有更詳細(xì)的討論。參考： 德 D. 羅伯. 計(jì)算材料學(xué). 北京：化學(xué)工業(yè)出版社, 2002, 9 俞慶森,朱龍觀. 分子設(shè)計(jì)導(dǎo)論.北京：高等教育出版社,2000 楊小震. 分子模擬與高分子材料. 北京：科學(xué)出版社，2002

8、熊家炯主編. 材料設(shè)計(jì). 天津：天津大學(xué)出版社，2000 Sun H，Ren P, Fried J R. The COMPASS Force Field: Parameterization and Validation for polyphosphazenes. Computational and Theoretical Polymer Science, 1998, 8(1/2): 229 Sun H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Det

9、ails on Alkane and Benzene Compounds. J. Phys. Chem., 1998, 102: 7338,2、分子力學(xué)方法 MM是確定分子結(jié)構(gòu)的方法 利用分子勢(shì)能隨結(jié)構(gòu)的變化而變化的性質(zhì)，確定分子勢(shì)能極小時(shí)的平衡結(jié)構(gòu)(stationary point) 物理模型: 視原子為質(zhì)點(diǎn)，視化學(xué)鍵為彈簧，而彈力常數(shù)完全由數(shù)據(jù)庫中的分子力場來確定 因此是直接用勢(shì)函數(shù)研究問題，不考慮原子的動(dòng)能 不考慮動(dòng)能所對(duì)應(yīng)的結(jié)構(gòu)，相當(dāng)于體系處于 T=0 K 時(shí)的結(jié)果,由力場首先構(gòu)造并得到分子的勢(shì)函數(shù)E(x) 利用在E(x)的極小點(diǎn)處，E(x)隨各原子獨(dú)立的空間坐標(biāo)xi, i=1,2,3

10、,3N-6 (與內(nèi)坐標(biāo)ri、i、i等價(jià)) 的一階微分=0以及全部二階微分0的數(shù)學(xué)條件：,進(jìn)行結(jié)構(gòu)優(yōu)化 具體步驟： 1) 用各種方法構(gòu)造出一個(gè)分子的任意結(jié)構(gòu)，得到初始的結(jié)構(gòu)參數(shù)xi0, i=1,2,3,3N-6 2) 進(jìn)行坐標(biāo)變換，即根據(jù)需要將原子的直角坐標(biāo)轉(zhuǎn)變成內(nèi)坐標(biāo)ri、i、i或反之,1) 用各種方法構(gòu)造出一個(gè)分子的任意結(jié)構(gòu)，得到初始的結(jié)構(gòu)參數(shù)xi0, i=1,2,3,3N-6 2) 進(jìn)行坐標(biāo)變換，即根據(jù)需要將原子的直角坐標(biāo)轉(zhuǎn)變成內(nèi)坐標(biāo)(鍵長、鍵角、二面角)或反之 3) 建立分子體系的勢(shì)能表達(dá)式E(x) 4) 計(jì)算E(x)隨各坐標(biāo)的一階、二階導(dǎo)數(shù) 5) 計(jì)算接近數(shù)學(xué)條件的坐標(biāo)增量 6) 得到

11、新的結(jié)構(gòu)參數(shù)xi1, i=1,2,3,3N-6 重復(fù) 4)、5)、6)，直至最后兩次得到的體系勢(shì)能之差或總體勢(shì)能梯度的均方根值達(dá)到預(yù)定精度范圍為止 MM計(jì)算中，有時(shí)還要考慮體系所處的外壓條件，此時(shí)，可通過壓力因子的定義來調(diào)節(jié)原子坐標(biāo),(i=1,2,3,3N-6),3、分子動(dòng)力學(xué)方法 分子力場是分子的靜態(tài)勢(shì)函數(shù) 而實(shí)際過程通常是在一定溫度和一定壓力下發(fā)生的 為了更切實(shí)際地了解體系運(yùn)動(dòng)和演化的過程，必須考慮體系中原子的運(yùn)動(dòng)，并與溫度T和時(shí)間t建立聯(lián)系 我們知道，溫度是原子分子熱運(yùn)動(dòng)劇烈程度的量度 根據(jù)統(tǒng)計(jì)熱力學(xué)，對(duì)于n個(gè)原子的體系，體系的溫度T與各原子的運(yùn)動(dòng)速率vi的關(guān)系為：,又因體系中各原子的速

12、率為vi時(shí)，動(dòng)量pi=mivi，對(duì)應(yīng)總動(dòng)能K(p)為：,勢(shì)能由力場確定為E(x)，因此體系的Hamilton量H為：,與量子力學(xué)不同，經(jīng)典力學(xué)對(duì)Hamilton量不進(jìn)行算符化處理，也不建立和求解本征方程，而是建立并求解經(jīng)典運(yùn)動(dòng)方程：,計(jì)算過程一般為：在一定的統(tǒng)計(jì)系綜下 1) 由原子位置和連接方式，從數(shù)據(jù)庫調(diào)用力場參數(shù)并形成體系勢(shì)函數(shù) 2) 由給定溫度計(jì)算體系動(dòng)能以及總能量 3) 計(jì)算各原子的勢(shì)能梯度, 得到原子在力場中所受的力 即dp/dt=mdv/dt=ma=F,1) 由原子位置和連接方式，調(diào)用力場參數(shù)并形成體系勢(shì)函數(shù) 2) 由給定溫度計(jì)算體系動(dòng)能以及總能量 3) 計(jì)算各原子的勢(shì)能梯度, 得

13、到原子在力場中所受的力 即 dp/dt=mdv/dt=ma=F 4) 對(duì)每個(gè)原子，在一定時(shí)間間隔內(nèi)，用牛頓方程求解其運(yùn)動(dòng)行為：,5) 顯示體系能量和結(jié)構(gòu) 6) 取下一時(shí)間間隔，返回步驟 1) 不斷循環(huán)反復(fù)，可設(shè)定循環(huán)次數(shù)或強(qiáng)行終止計(jì)算 其中, 不僅時(shí)間間隔可以根據(jù)需要取不同大小(一般1 fs=10-15 s),溫度可以任意設(shè)定，而且還可以在循環(huán)過程中逐漸改變溫度，即研究體系的退火(annealing)行為,有時(shí)需要進(jìn)一步考慮外場的作用，如壓力、電場、磁場、重力場等 從原則上講，這些問題都不難解決，而且還在進(jìn)一步發(fā)展中，不再贅述 此外，即使是分子力學(xué)與分子動(dòng)力學(xué)方法，也受計(jì)算量的限制，所處理的體

14、系不可能太大，好的計(jì)算機(jī)可達(dá)到數(shù)十萬百萬個(gè)原子的規(guī)模 (2050 nm以下/ MD 106108個(gè)原子羅伯p110) 對(duì)于更大(50 nm、微米)尺度的問題，人們也在發(fā)展介觀尺度(mesoscale)的方法，并且已取得一些成果 分子動(dòng)力學(xué)方法也有一些本質(zhì)的缺陷 如：勢(shì)函數(shù)精度(客觀性) 勢(shì)函數(shù)形式在每次計(jì)算中都不變, 故不能模擬如分子在高溫下, 結(jié)構(gòu)發(fā)生斷裂的熱裂解過程 (國內(nèi)引入開關(guān)函數(shù)已開始做) 此時(shí), 可用QM處理化學(xué)變化或e轉(zhuǎn)移, 余用MM，稱QM/MM,二、 MM、MD計(jì)算程序 因程序結(jié)構(gòu)相對(duì)簡單，國內(nèi)外有不少, 且還在不斷產(chǎn)生 如：Chem3D, ChemOffice, Ceriu

15、s2, Material Studio, Alchemy 2000, Sybyl, Biosym, HyperChem x.x, Spartan x.x, Chemgraf, Bilder, Script, COGS, Gaussian-03; MP(Molecular Properties, 楊小震)等 輸入：體系模型建模(Builder)，選擇力場、系綜 MM: 無特殊輸入 MD: 循環(huán)次數(shù)(萬幾10萬)、溫度、時(shí)間間隔(fs), 外場,溫度變化速率，其它性質(zhì) 等 輸出： MM: 勢(shì)能-結(jié)構(gòu)曲線(數(shù)值)，動(dòng)態(tài)結(jié)構(gòu)圖, 其它性質(zhì) MD: (同上)，及 其它MD性質(zhì) 注：第一原理及MM一般只能

16、得到勢(shì)能極小點(diǎn)結(jié)構(gòu)，而MD可越過一些小的勢(shì)壘，甚至可達(dá)到最小點(diǎn),三、 MM、MD方法的應(yīng)用 領(lǐng)域：高分子、生命科學(xué)、藥物設(shè)計(jì)、催化、半導(dǎo)體 其它功能材料、結(jié)構(gòu)材料等 分子力學(xué)是用計(jì)算機(jī)在原子水平上模擬給定分子模型的結(jié)構(gòu)與性質(zhì)，進(jìn)而得到分子的各種物理性質(zhì)與化學(xué)性質(zhì)，如結(jié)構(gòu)參數(shù)、振動(dòng)頻率、構(gòu)象能量、相互作用能量、偶極矩、密度、摩爾體積、汽化焓等 分子動(dòng)力學(xué)方法能實(shí)時(shí)將分子的動(dòng)態(tài)行為顯示到計(jì)算機(jī)屏幕上, 便于直觀了解體系在一定條件下的演變過程 MD含溫度與時(shí)間, 因此還可得到如材料的玻璃化轉(zhuǎn)變溫度、熱容、晶體結(jié)晶過程、輸送過程、膨脹過程、動(dòng)態(tài)弛豫(relax)以及體系在外場作用下的變化過程等,分子模擬較早應(yīng)用于高分子問題的研究 應(yīng)用范圍主要包括: 物理性質(zhì)、結(jié)構(gòu)、構(gòu)象與彈性、晶體結(jié)構(gòu)、力學(xué)性能、玻璃態(tài)與玻璃化轉(zhuǎn)變、光譜性質(zhì)、非線性光學(xué)性質(zhì)、電性質(zhì)、共混與分子間相互作用等 應(yīng)用于生物科學(xué)和藥物設(shè)計(jì)也十分普及，如蛋白質(zhì)的多級(jí)結(jié)構(gòu)與性質(zhì)，病毒、藥物作用機(jī)理、特效藥物的大通量篩選與快速開發(fā)等 在化學(xué)領(lǐng)域，用于表面催化與催化機(jī)理、溶劑效應(yīng)、原子簇的結(jié)構(gòu)與性質(zhì)研究等 在材料科學(xué), 用于材料的優(yōu)化設(shè)計(jì)、結(jié)構(gòu)與力學(xué)性能、熱加工性能預(yù)報(bào)、界面相互作用、納米材料結(jié)構(gòu)與性能研究等,羅伯，計(jì)算材料學(xué)，p129: 已有很多研究者采用分子動(dòng)力學(xué)在原子尺度上進(jìn)行了材料