分子動(dòng)力學(xué)模擬入門

1、.1 多體系統(tǒng)分子動(dòng)力學(xué) Multi-body System Molecular Dynamics .2 引 言 物質(zhì)基本構(gòu)成分子、原子 在分子、原子這個(gè)微觀水平上來(lái)考察物質(zhì)：多體 世界 查清楚微觀世界，宏觀就清楚了 從微觀考慮問(wèn)題的現(xiàn)實(shí)可行性 從微觀考慮問(wèn)題的必要性 物性的觀測(cè)性參數(shù)：熱傳導(dǎo)、溫度、壓力、粘 性、. .3 微觀處理的前提 已知微觀粒子間的相互作用 假設(shè) 分子為球，惰性，分子間的作用只取決于分子間 的距離 分子動(dòng)力學(xué) (Molecular Dynamics，MD) .4 MD的應(yīng)用 領(lǐng)域：物理、化學(xué)、生物、材料等 MD方法能實(shí)時(shí)將分子的動(dòng)態(tài)行為顯示到計(jì)算機(jī)屏 幕上, 便于直觀了

2、解體系在一定條件下的演變過(guò)程 MD含溫度與時(shí)間, 因此還可得到如材料的玻璃化轉(zhuǎn) 變溫度、熱容、晶體結(jié)晶過(guò)程、輸送過(guò)程、膨脹過(guò) 程、動(dòng)態(tài)弛豫(relax)以及體系在外場(chǎng)作用下的變化 過(guò)程等 水和離子在微小硅孔中的運(yùn)動(dòng) 聚乙烯的結(jié)晶 .5 MD的基本原理 用牛頓經(jīng)典力學(xué)計(jì)算許多分子在相空間中的 軌跡 求解系統(tǒng)中的分子或原子間作用勢(shì)能和系統(tǒng)外加 約束共同作用的分子或原子的牛頓方程。 模擬系統(tǒng)隨時(shí)間推進(jìn)的微觀過(guò)程。 通過(guò)統(tǒng)計(jì)方法得到系統(tǒng)的平衡參數(shù)或輸運(yùn)性質(zhì) 計(jì)算程序較為復(fù)雜，占用較多內(nèi)存 .6 MD的主要步驟 選取要研究的系統(tǒng)及其邊界，選取系統(tǒng)內(nèi)粒 子間的作用勢(shì)能模型 設(shè)定系統(tǒng)中粒子的初始位置和初始動(dòng)

3、量 建立模擬算法，計(jì)算粒子間作用力及各粒子 的速度和位置 當(dāng)體系達(dá)到平衡后，依據(jù)相關(guān)的統(tǒng)計(jì)公式， 獲得各宏觀參數(shù)和輸運(yùn)性質(zhì) .7 分子間勢(shì)能及相互作用 N個(gè)粒子系統(tǒng)的總勢(shì)能 N ji ij NN N N rV rV rVrV rVrVrVV 1 1,- 223 11312 )( )( )()( )()()( .8 d V r V r 剛球模型 ( ) 0 rd V r rd ( )( ) d V r r 斥力力心點(diǎn)模型 d r V - ( ) ( ) rd V r d rd r Southerland模型 .9 分子間勢(shì)能及相互作用 Lennard-Jones勢(shì)能 能量尺度； 長(zhǎng)度尺度 為方便

4、，時(shí)常歸一化： 記 ； 612 4)( rr rV 612 / 1 / 1 4 )( rr rV VV/rr/ 612 11 4)( rr rV V(r) F(r) 排斥力 V(r), F(r) r 吸引力 .10 分子間勢(shì)能及相互作用 一些氣體的參數(shù) kB=1.38x1023(J/K): Boltzmann常數(shù) NeonArgonKryponXenonNitrogen (nm)0.2750.34050.3600.4100.370 /kB(K)36119.817122195 .11 分子間勢(shì)能及相互作用 相互作用 標(biāo)量形式： 直角坐標(biāo)： 至此，各粒子間相互作用已知，可進(jìn)行模擬了 )(rV rm

5、 f ( )( ) , y x f fV r xV r y mrrmrr )(rV m F .12 模擬的數(shù)學(xué)方法 Euler法和Euler-Cromer方法？ 不能用：不能保持總能量守恒 Verlet算法：速度形式 1 2 1 2 12 1 2 1 11 2 1 2 1 12 nnnn nnnn nn n nn n xxvtat vvaat vvat vvat .13 模擬的數(shù)學(xué)方法 Leap-frog算法： x的截?cái)嗾`差為 ，v的截?cái)嗾`差 1/21/2 11/2 1/2 2 nnn nnn t nnn vvta xxtv vva 4 ()Ot 2 ()Ot .14 模擬的數(shù)學(xué)方法 邊界條件

6、 模擬能力限制，不能模擬大量分子，只能模擬有 限空間中的有限個(gè)分子：有限空間邊界 固體（剛性）邊界條件 不僅僅有分子間的相互作 用，還引入了壁面的作用 分子量大時(shí)，壁面作用可 忽略不計(jì) 2 3 3 611a aaN 和壁面作用分子數(shù)壁面積 總分子數(shù)體積 3 aN 總分子數(shù) 壁面積和壁面作用分子數(shù) .15 模擬的數(shù)學(xué)方法 取 ， 前比值為0.20.01。取前值，模 擬粗糙；取后值，模擬計(jì)算量太大 處理方法：使用周期性邊界條件 周期性邊界條件 2 10N 6 10 a .16 模擬的數(shù)學(xué)方法 兩個(gè)不同粒子在x或y方向上的最大分離距離為a/2 最小像約定：兩粒子分離距離最大分離距離， 相互作用力可以

7、忽略，而加入其中像粒子之一相 互作用力來(lái)考慮 ? 2 a 2/a a .17 模擬的數(shù)學(xué)方法 考慮的粒子總數(shù)不變 初始條件 隨機(jī)初始條件給法之一 要求 大?。?條 件 一： 規(guī) 則 給 法 條 件 二： 隨 機(jī) 給 法 max |Vv ) 1random2( max Vv .18 模擬的數(shù)學(xué)方法 random：隨機(jī)數(shù)產(chǎn)生函數(shù)，產(chǎn)生(0，1)之間的隨 機(jī)數(shù)。 方向(按球坐標(biāo)給法)： 分量 ： x y z )5 . 0random(2arccos( )5 . 0random(2arccos)0.5random(sign 180180;1800 cossinvvx sinsinvvy cosvvz

8、.19 模擬 微觀量 溫度 根據(jù)統(tǒng)計(jì)熱力學(xué)，平衡態(tài)下經(jīng)典系統(tǒng)的能量中的 每一個(gè)二次項(xiàng)具有平均值kBT/2,即 注意：上式在系統(tǒng)質(zhì)心速度為0時(shí)適用 i iiB vmTNk d 2 2 1 2 空間維數(shù) 粒子個(gè)數(shù) ：取時(shí)間平均 .20 模擬 問(wèn)題：如何給定系統(tǒng)的初始條件，得到所需要的平 衡態(tài)溫度Teq？ 解決方法之一：速度標(biāo)定法 任給初始條件，模擬到平衡，得到系統(tǒng)平衡態(tài)溫 度T。一般TTeq。令 用速度 再模擬直到平衡，若所得溫度仍不等于Teq，再進(jìn) 行上述過(guò)程 TTf eq / 1/d ii fvv .21 給定初始條件：xi，vi 計(jì)算到平衡態(tài) |T - Teq| f = Teq/ T vi

9、= vi f 1/d 計(jì)算結(jié)束 Y N .22 其它方法：Gaussian熱浴法（約束溫度調(diào)節(jié)方法） 其基本原理在運(yùn)動(dòng)方程中加入“摩擦力”項(xiàng)， 并將其與粒子速度聯(lián)系起來(lái)。 平衡態(tài)時(shí)，系統(tǒng)溫度不變，因此dEk/dt=0 iii fmamv 0 ii i v a 2 ii i i i f v mv .23 宏觀性質(zhì)的統(tǒng)計(jì) 系統(tǒng)的勢(shì)能 系統(tǒng)的內(nèi)能 系統(tǒng)的總能 E = Ep+Ek 系統(tǒng)的溫度 1 ( ) pij ij N EV r 2 2 k i i p E m 2 1 ii i B Tmv dNk .24 模擬 熱容 定義熱容 計(jì)算系統(tǒng)在溫度T和T+T時(shí)的總能ET、ET +T， V v T E C

10、E：系統(tǒng)總能 T EE C TTT v .25 模擬 壓強(qiáng) 對(duì)壁面的壓強(qiáng) t時(shí)間里作用在單位面積壁上的壓力 t i v 時(shí)刻，速度為 tt i v 時(shí)刻，速度為 i ii i t vv m dA p 1 剛 性 壁 .26 模擬 粒子速度分布 選速度間隔v，模擬nt個(gè)時(shí)間步，記錄在每個(gè)速 度間隔中的粒子數(shù)，最后歸一化。 v v )(vN vvvv2 .27 模擬 氣、液狀態(tài)方程 維里定理(Virial Theorem) i iiB Fr d TNkpV 1 壓強(qiáng) 體積 粒子i的位置矢量 粒子i所受到的其 它粒子的合相互 作用力 溫度的模 擬 可得此項(xiàng) 在溫度的模擬基 礎(chǔ)上再模擬此項(xiàng) .28 模

11、擬 例：用此可確定高密度氣體和液體狀態(tài)方程(van der Waals方程)中的系數(shù) 理想氣體狀態(tài)方程在高密度情況下不可用 確定系數(shù)a和b 2 1 a b Tkp B 氣體密度 .29 CASE(1) - Couette Flow Size of domain is:12.51x7.22x16.71 1,1,1.5,0.6, 0.02,1.4 llwwwlllwwwl wall FT 612 612 4() ijij ijij E rr 612 712 4( 612) ijij ijij dE F drrr If is less than 0, then the two species are

12、 immiscible. i, j represent different species .30 CASE(1) - Couette Flow z time .31 CASE(1) - Couette Flow z .32 CASE(2) - Contact Angle Simulation Mass: m1=1, m2=8, m3=0.8 L=25.05, W=6.56 H=10.29 T=1.2 110.2 111 0.211 ij 126 ( ) ijij LJijij Vr rr i, j=1,2,3, 1 red fluid, 2 - wall, 3 green fluid .33

13、 V .34 CASE(3)- Rayleigh-Taylor instability generation(重力場(chǎng)) Dzwinel, W., Alda, W., Pogoda, M., and Yuen, D.A., 2000, Turbulent mixing in the microscale: a 2D molecular dynamics simulation, Physica D, Vol. 137, pp. 157-171. Fig. 1. The snapshots of MD simulations of the RT instability for two particl

14、e systems: (A) closed; (B) open (the heavy fluid is coloured in light grey while the lighter one is dark grey, the black part of the figure is empty, i.e., it does not contain particles). Fig. 2. The effect of converging geometry obtained by MD simulation of one million particles in the microscale.

15、0.5 m .35 CASE (4)Typical translocation event A 1.4V bias applied to membrane. 20 base-pair fragment of double stranded DNA placed in front of a nanopore. End of DNA nearest to the pore is pulled into the pore by its charged backbone (a,b) System reaches a meta-stable state (c) and translocation hal

16、ts. Base-pairs start to split. Some freed nucleotides adhere to pore surface. Voltage increased momentarily to drive system out of metastable state. DNA exits pore. One of the bases holds firmly to the pore surface. After 50ns, most of DNA has left pore. Nine of twenty base pairs are split. .36 Bubb

17、le nucleation on solid surfaces Maruyama, S, and Kimura, T., 2000, A Molecular Dynamics Simulation of a Bubble Nucleation on Solid Surface, Heat and Technology, Vol. 18, pp. 69-73. Molecular transport in droplets Maruyama, S., Matsumoto, S., and Ogita, A., 1994, Surface Phenomena of Molecular Clusters by Molecular Dynamics Method, Thermal Science and