分子動力學模擬方法

1、1957年年：基于剛球勢的分子動力學法（：基于剛球勢的分子動力學法（Alder and Wainwright） 1964年年：利用：利用Lennard-Jone勢函數(shù)法對液態(tài)勢函數(shù)法對液態(tài)氬氬性質的模擬（性質的模擬（Rahman） 1971年年：模擬具有分子：模擬具有分子團簇團簇行為的行為的水水的性質（的性質（Rahman and Stillinger） 1977年年：約束動力學方法（：約束動力學方法（Rychaert, Ciccotti & Berendsen; van Gunsteren） 1980年年：恒壓條件下的動力學方法（：恒壓條件下的動力學方法（Andersen法、法、Pa

2、rrinello-Rahman法）法） 1983年年：非平衡態(tài)動力學方法（：非平衡態(tài)動力學方法（Gillan and Dixon） 1984年年: 恒溫條件下的動力學方法恒溫條件下的動力學方法(Berendsen et al.) 1984年年：恒溫條件下的動力學方法（：恒溫條件下的動力學方法（Nos-Hoover法）法） 1985年年：第一原理分子動力學法（：第一原理分子動力學法（Car-Parrinello法）法） 1991年年：巨正則系綜巨正則系綜的分子動力學方法（的分子動力學方法（Cagin and Pettit）分子動力學簡史分子動力學簡史粒子的運動取決于經(jīng)典力學粒子的運動取決于經(jīng)典力

3、學(牛頓定律（牛頓定律（F=ma）原理：原理：計算一組分子的相空間軌道，其中每個分子各自服從計算一組分子的相空間軌道，其中每個分子各自服從牛頓運動定律：牛頓運動定律： 11112)(21NiNijijNiiirUmpHiiiiimdtdmvrp 111111)()(NiNijijNiNijijirUrdtdijrFp)0(0itirr)0(0itidtvdr初始條件：初始條件：粒子位置的粒子位置的Taylor展開式：展開式：粒子位置粒子位置：粒子速度粒子速度： t (t) t)-(t(t)2 t)(t2iarrriii t2t)-(tt)(t (t)iiirrv粒子加速度：粒子加速度：3i2i

4、it (t) 61 t (t) 21 t (t) (t) t)(tbavrrii3i2iit (t) 61 t (t) 21 t (t) (t) t)(tbavrrii (t) (t)imiiFa開始運動時需要開始運動時需要r(t-t)：t (0) (0) t)(ivrr+n 算法啟動算法啟動u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 擾動初始位置：擾動初始位置：u 計算第計算第n步的力步的力u 計算第計算第n+1步的位置：步的位置：u 計算第計算第n步的速度：步的速度：u 重復重復至至t (0) (0) t)(ivrr t (t) t)-(t(t)2 t)(t2iarrrii

5、i t2t)-(tt)(t (t)iiirrvDo 100 I = 1, N RXNEWI = 2.0 * RX(I) RXOLD(I) + DTSQ * AX(I) RYNEWI = 2.0 * RY(I) RYOLD(I) + DTSQ * AY(I) RZNEWI = 2.0 * RZ(I) RZOLD(I) + DTSQ * AZ(I) VXI = ( RXNEWI RXOLD(I) ) / DT2 VYI = ( RYNEWI RYOLD(I) ) / DT2 VZI = ( RZNEWI RZOLD(I) ) / DT2 RXOLD(I) = RX(I) RYOLD(I) = RY

6、(I) RZOLD(I) = RZ(I) RX(I) = RXNEWI RY(I) = RYNEWI RZ(I) = RZNEWI100 CONTINUEt (t) t)21-(t t)21(tiavviit t)21(t(t) t)(tivrrii2t)21-(tt)21(t (t)iiivvv1. 首先利用當前時刻的加速度，計算半個時間步長后的速度：首先利用當前時刻的加速度，計算半個時間步長后的速度：2. 計算下一步長時刻的位置：計算下一步長時刻的位置：3. 計算當前時刻的速度：計算當前時刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv開始運動時需要開始運動時需要v(-t/

7、2)：t/2 (0) (0) t/2)(iavvn 算法啟動算法啟動u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 擾動初始速度：擾動初始速度：u 計算第計算第n步的力步的力u 計算第計算第n+1/2步的速度：步的速度：u 計算第計算第n+1步的位置：步的位置：u 計算第計算第n步的速度：步的速度：u 重復重復至至t/2 (0) (0) t/2)(iavvt (t) t)21-(t t)21(tiavviit t)21(t(t) t)(tivrrii2t)21-(tt)21(t (t)iiivvv2iit (t) 21 t (t) (t) t)(tavrriitt)(t (t) 2

8、1 (t) t)(tiiaavviit(t) 21 (t) t)21(tiavviitt)(t 21 t)21(t t)(tiavviit t)21(t (t) t)(tivrrii等價于等價于n 算法啟動算法啟動u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 計算第計算第n+1步的位置：步的位置：u 計算第計算第n+1步的力步的力u 計算第計算第n+1步的速度：步的速度：u 重復重復至至2iit (t) 21 t (t) (t) t)(tavrriitt)(t (t) 21 (t) t)(tiiaavviiVerletLeap-frogVelocity Verlet 預測預測(P

9、redictor)階段：其基本思想是階段：其基本思想是Taylor展開，展開，)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrr根據(jù)新的原子位置根據(jù)新的原子位置r rp p，可以計算獲得校正后的，可以計算獲得校正后的a ac c( (t t+ + t t),),定義預測誤差定義預測誤差: :)()()(ttttttpcaaa)()()()()()()()()()()()(3210ttctttttctttttctttttctttpcpcpcpcabbaaaavvarr利用此預測誤差，對預測

10、出的位置、速度、加速度等量進行校正：利用此預測誤差，對預測出的位置、速度、加速度等量進行校正： 校正校正(Corrector)階段：階段：u預測階段運動方程的變換：預測階段運動方程的變換：定義一組矢量定義一組矢量:333322221061d)(d21d)(dd)(d)(ttttttttttrrrrrrrr)()()()(1000310032101111)()()()(32103210ttttttttttttpppprrrrrrrr)()()()()()(21)()()()(61)(21)()()(232ttttttttttttttttttttttttppppbbbaabavvbavrru校正階

11、段運動方程的變換：校正階段運動方程的變換：rrrrrrrrr )()()()()()()()(321032103210ccccttttttttttttttttppppccccr的形式：的形式：C0， C1， C2， C3的值以及的值以及C0，取決于運動方程的階數(shù)。取決于運動方程的階數(shù)。n 一階運動方程：一階運動方程：)(rrf)()(11ttttpcrrrn 二階運動方程之一：二階運動方程之一：)(rrf )()(22ttttpcrrrn 二階運動方程之二：二階運動方程之二：),(rrr f)()(22ttttpcrrr室溫下，室溫下， t 1 fs (femtosecond 10-15s)，

12、溫，溫度越高，度越高，t 應該減小應該減小n 太長的時間步長會造成分子間的激烈碰撞，體系數(shù)據(jù)溢出太長的時間步長會造成分子間的激烈碰撞，體系數(shù)據(jù)溢出; ; n 太短的時間步長會降低模擬過程搜索相空間的能力太短的時間步長會降低模擬過程搜索相空間的能力 u 它是分子動力學方法的最基本系綜它是分子動力學方法的最基本系綜u 具有確定的粒子數(shù)具有確定的粒子數(shù)N，能量，能量E和體積和體積Vu 算法：算法：u 規(guī)定初始位置和初始速度規(guī)定初始位置和初始速度u 對運動方程積分若干步對運動方程積分若干步u 計算勢能和動能計算勢能和動能u 若能量不等于所需要的值，對速度進行標度若能量不等于所需要的值，對速度進行標度u

13、 重復重復至，直到系統(tǒng)平衡至，直到系統(tǒng)平衡微正則系綜微正則系綜(NVE)MD模擬算法的流程圖：模擬算法的流程圖：給定每個分子的初始位置給定每個分子的初始位置ri(0)和速度和速度vi(0)計算每個分子的受力計算每個分子的受力Fi和加速度和加速度ai解運動方程并求出每個分子運動一個時間步解運動方程并求出每個分子運動一個時間步長后到達的位置所具有的速度長后到達的位置所具有的速度統(tǒng)計系統(tǒng)的熱力學性質及其它物理量統(tǒng)計系統(tǒng)的熱力學性質及其它物理量統(tǒng)計性質不變？統(tǒng)計性質不變？ 打印結果，結束打印結果，結束YesNo移動所有分子到新的位置并具有當前時刻的移動所有分子到新的位置并具有當前時刻的速度速度微正則系

14、綜微正則系綜MD模擬程序模擬程序F3講解（講解（LJ, NVE）：）：無因次量：無因次量：/*rr 3*/*EE /*ff/*3PP /*kTT 2/12*)/(mtt 2/1*)/(vvmMD模擬中幾個熱力學量的計算：模擬中幾個熱力學量的計算：對于由對于由N個單原子組成的系統(tǒng)：個單原子組成的系統(tǒng)： 動能和溫度：動能和溫度：TNkmEBNiK23v2112i采用對比量：采用對比量：/*kTT 2/1*)/(vvm/*EE *12*i*23v21NTENiK21114) ru 2)(r(drNrUUUUcrNiNijijclrcc 對于對于LJ流體：流體：3933138 cclrcrrNU 勢能

15、：勢能： 4)(612rrrUc/*rr 3*/*EE 采用對比量：采用對比量： 114)(6*12*rrrUc3*9*113138 ccrrNUlrc 內能：內能：KEUE內能由勢能和動能組成：內能由勢能和動能組成：采用對比量：采用對比量：*KEUE 壓力：壓力：393232316 cclrcrrPlrccBPVWTkP/*rr 3*/*EE 采用對比量：采用對比量：/*3PP /*ff 111)(31NiNijijccrWWijijijijijijcrUrWrrfr)(31 )(練習：練習：推導推導LJ流體分子間力的表達式流體分子間力的表達式(fx, fy, fz及其對比量及其對比量)：

16、4)(612rrrU勢能函數(shù)形式：勢能函數(shù)形式：力：力：r)()(rUrf2612 2184 rrrrrrrUrUf/*rr 采用對比量：采用對比量：/*ff =x, y, z2*6*12* 121184 rrrrf* )(fr rWcLJ分子間的維里項：分子間的維里項：2*6*12* 121184 rrrrf =x, y, z 121184 )(6*12*rrrWc2* )( rrrWfct (t) t)21-(t t)21(timFvviit t)21(t(t) t)(tivrrii2t)21-(tt)21(t (t)iiivvv(t) ia/*ff2/12*)/(mtt 2/1*)/(v

17、vm采用對比量：采用對比量：/*rr 2/12*2/1*2/1*t )(t )t21-(t )t21(tmmmm*i*i*iFvv2/12/12*t )(t )t21-(t )t21(tmmm*i*i*iFvv*t )(t )t21-(t )t21(t*i*i*iFvv最終得到：最終得到：*t )t21(t )(t )t(t*i*i*ivrr2)t21-(t)t21(t (t)*i*i*ivvv同理得到：同理得到：mTkvBx2根據(jù)能量均分原理，可知：根據(jù)能量均分原理，可知：2/12*2/122/122/1333vTvmTkvNmTNkTTreqrefBrefBref標度因子標度因子:對比量對

18、比量速度標度速度標度:oldnewvv *old*newvv 或或微正則系綜微正則系綜MD模擬程序模擬程序F3講解（講解（LJ, NVE）：）：初始化：初始化：READ (*,(A) TITLE ! 運行作業(yè)題目運行作業(yè)題目READ (*,*) NSTEP ! 運行步數(shù)運行步數(shù)READ (*,*) IPRINT ! 打印步數(shù)打印步數(shù)READ (*,(A) CNFILE ! 位型文件位型文件READ (*,*) DENS ! 對比密度對比密度READ (*,*) RTEMP ! 對比溫度對比溫度 READ (*,*) RCUT ! 對比截斷半徑對比截斷半徑READ (*,*) DT ! 對比時間

19、步長對比時間步長CALL READCN ( CNFILE ) 初始位型：初始位型：u面心立方面心立方 (face-centered cubic, FCC)：每面中心有一格點每面中心有一格點 u體心立方體心立方 (body-centered cubic, BCC)：u簡單立方簡單立方 (simple cubic, SC)：XL初始位型：初始位型：面心立方面心立方 ( FCC) (程序程序F23)NC=(REAL(N)/4.0)*(1.0/3.0)XL = 1.0 / REAL ( NC )Y = 0.5 * XLR(1)=(0, 0, 0) R(2)=(0, Y, Y)R(3)=(Y, 0, Y

20、) R(4)=(Y, Y, 0)M = 0DO 10 I = 1, NCDO 10 J = 1, NCDO 10 K = 1, NC DO 11 IJ = 1, 4 RX(IJ+M)=RX(IJ) + XL*(K-1) RY(IJ+M)=RY(IJ) + XL*(J -1) RZ(IJ+M)=RZ(IJ) + XL*(I -1) 11 CONTINUE M = M + 410 CONTINUEDO 100 I = 1, N RX(I) = RX(I) - 0.5 RY(I) = RY(I) - 0.5 RZ(I) = RZ(I) - 0.5100 CONTINUE將模擬盒子的中心移到原點：將模

21、擬盒子的中心移到原點： 初始速度：初始速度：n 簡單的選擇：簡單的選擇：V random(-0.5, 0.5) =x, y, zmTkvB22/1*2/12*2/122/1TTvTvmTkTTreqreqrefBref標度因子標度因子:*old*newvv 速度標度速度標度:*2*TvFACTOR = SQRT ( RTEMP )DO 100 I = 1, N VX(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) VY(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) VZ(I) = FACTOR * ( RANF(DUMMY) - 0.5 )10

22、0 CONTINUE*2*Tv隨機安排初始速度：隨機安排初始速度：標度初始速度：標度初始速度：SUMKX = 0.0SUMKY = 0.0SUMKZ = 0.0DO 200 I = 1, N SUMKX = SUMKX + VX(I)*2 SUMKY = SUMKY + VY(I)*2 SUMKZ = SUMKZ + VZ(I)*2200 CONTINUEBEITAX =SQRT(RTEMP/SUMKX)BEITAY =SQRT(RTEMP/SUMKY) BEITAZ =SQRT(RTEMP/SUMKZ)DO 300 I = 1, N VX(I) = VX(I) * BEITAX VY(I)

23、= VY(I) * BEITAY VZ(I) = VZ(I) * BEITAZ300 CONTINUE2/1*2/12*TTvTreqreq標度因子標度因子:*old*newvv SUMX = 0.0SUMY = 0.0SUMZ = 0.0DO 200 I = 1, N SUMX = SUMX + VX(I) SUMY = SUMY + VY(I) SUMZ = SUMZ + VZ(I) CONTINUESUMX = SUMX / REAL ( N )SUMY = SUMY / REAL ( N )SUMZ = SUMZ / REAL ( N ) DO 300 I = 1, N VX(I) =

24、 VX(I) - SUMX VY(I) = VY(I) - SUMY VZ(I) = VZ(I) - SUMZ300 CONTINUE控制體系的總動量為零：控制體系的總動量為零：n 從從Maxwell分布中抽樣：分布中抽樣：x)(xP0 x xdx高斯高斯(Gauss)分布分布:對于等幾率隨機試驗對于等幾率隨機試驗(Bernoulli試驗試驗)，大量的試驗結果滿足高斯分布大量的試驗結果滿足高斯分布 222)(2exp)(21xxxxxP麥克斯韋速度分布定律麥克斯韋速度分布定律: dvekTmNdvvfTkmvB22122mTkvB2由于：由于： 222)(2exp)(21xxxxxP =x,

25、y, z單位體積的分子再每個分量上的速度分布實際上就是單位體積的分子再每個分量上的速度分布實際上就是高斯分布。高斯分布。n 從從Maxwell分布中抽樣：分布中抽樣：高斯高斯(Gauss)分布的隨機數(shù)生成方法分布的隨機數(shù)生成方法: 222)(2exp)(21xxxxxP生成隨機數(shù)：生成隨機數(shù)： i，i=1, 2, , 1246121iiR)R a R ) a R ) a )R a Ra ( 123252729SUM = 0.0DO 10 I = 1, 12 SUM = SUM + RANF ( DUMMY )10 CONTINUER = ( SUM - 6.0 ) / 4.0R2 = R *

26、RGAUSS = ( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 ) * R高斯高斯(Gauss)分布的隨機數(shù)生成分布的隨機數(shù)生成(程序程序F24)FACTOR = SQRT ( RTEMP )DO 100 I = 1, N VX(I) = FACTOR * GAUSS(DUMMY) VY(I) = FACTOR * GAUSS(DUMMY) VZ(I) = FACTOR * GAUSS(DUMMY) CONTINUE控制總動量為零：同前面一樣處理?？刂瓶倓恿繛榱悖和懊嬉粯犹幚?。從從Maxwell分布中抽樣分布中隨機安排初始速度：分布中抽樣

27、分布中隨機安排初始速度：微正則系綜微正則系綜MD模擬程序模擬程序F3講解講解(LJ, NVE) ： 量綱變換：量綱變換：SIGMA = ( DENS / REAL ( N ) ) * ( 1.0 / 3.0 )RCUT = RCUT * SIGMADT DT * SIGMADENS = DENS / ( SIGMA * 3 )333NVN*模擬盒子的邊長為模擬盒子的邊長為1VN3*L2/1*)/(mtt 長程校正：長程校正：微正則系綜微正則系綜MD模擬程序模擬程序F3講解講解(LJ, NVE) ：SR3 = ( SIGMA / RCUT ) * 3SR9 = SR3 * 3SIGCUB = S

28、IGMA * 3VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N ) : * ( SR9 - 3.0 * SR3 )WLRC = ( 16.0 / 9.0 ) * PI * DENS * SIGCUB * REAL ( N ) :* ( 2.0 * SR9 - 3.0 * SR3 )393*398cclrcrrNU393*32916 cclrclrcrrNVPW 算法：算法啟動算法：算法啟動微正則系綜微正則系綜MD模擬程序模擬程序F3講解講解(LJ, NVE) ：CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NE

29、WVC, NEWW )CALL MOVE ( -DT )CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W )CALL FORCE ( DT, SIGMA, RCUT, V, VC, W )CALL KINET ( OLDK )CALL MOVE ( DT )CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )CALL KINET ( NEWK ) 算法：差分格式：算法：差分格式：SR2 = SIGSQ / RIJSQVIJ = 4.0 * ( SR12 - SR6 )WIJ = 24.0 * ( 2.0 * SR12

30、- SR6 )VELIJ = WIJ * DT / RIJSQDVX = VELIJ * RXIJ. VXI = VXI + DVX.VX(J) = VX(J) DVXV = V + VIJW = W + WIJCALL MOVE(DT) 114)(6*12*crrrU 121184 )(6*12*rrrWc2* )( rrrWfc*t )(t )t21-(t )t21(t*i*i*iFvvDO 1000 I = 1, N RX(I) = RX(I) + VX(I) * DT RY(I) = RY(I) + VY(I) * DT RZ(I) = RZ(I) + VZ(I) * DT1000 C

31、ONTINUEMOVE(DT):*t )t21(t )(t )t(t*i*i*ivrr速度的標定（只用于平衡階段）速度的標定（只用于平衡階段）SUMK = 0.0DO 200 I = 1, N SUMK = SUMKX + VX(I)*2+VY(I)*2 + VZ(I)*2200 CONTINUEBEITA =SQRT( 3.0 * RTEMP / SUMK )DO 300 I = 1, N VX(I) = VX(I) * BEITA VY(I) = VY(I) * BEITA VZ(I) = VZ(I) * BEITA300 CONTINUE2/1*2/12*3TTvTreqreq*old*

32、newvv u 具有確定的粒子數(shù)具有確定的粒子數(shù)N，溫度，溫度T和體積和體積Vu速度的直接標度速度的直接標度u熱浴方法熱浴方法 (Andersen Thermostat)u約束方法約束方法(阻尼力方法阻尼力方法)u系統(tǒng)擴展方法系統(tǒng)擴展方法(Extended Systems Method)u 問題的關鍵：溫度的約束問題的關鍵：溫度的約束Nose-Hoover方法方法一、熱浴方法一、熱浴方法 (Andersen Thermostat)u 引入一個與虛擬粒子碰撞的隨機力引入一個與虛擬粒子碰撞的隨機力u 想象系統(tǒng)浸在熱浴當中想象系統(tǒng)浸在熱浴當中u系統(tǒng)和熱浴間的相互作用強度由隨機碰撞的頻率決定系統(tǒng)和熱浴

33、間的相互作用強度由隨機碰撞的頻率決定u碰撞的幾率等于碰撞的幾率等于Nudtu如果一個粒子經(jīng)歷碰撞，它的速度將從約束溫度下的如果一個粒子經(jīng)歷碰撞，它的速度將從約束溫度下的Maxwell分布中隨機抽取分布中隨機抽取u 總能量和總動量均不守恒總能量和總動量均不守恒二、約束方法二、約束方法u 是等動能是等動能(Iso-Kinetics)分子動力學方法分子動力學方法u 系統(tǒng)的運動方程為：系統(tǒng)的運動方程為：m/pr pprfp),(u 引入阻尼系數(shù)引入阻尼系數(shù) 以保證將溫度約束在恒定值以保證將溫度約束在恒定值u 根據(jù)高斯最小約束原理：根據(jù)高斯最小約束原理：iiiii2pfp三、三、Nose-Hoover擴

34、展方法擴展方法基本思想：基本思想：設想原系統(tǒng)與一個耦合系統(tǒng)共同組成一個擴展系統(tǒng)，設想原系統(tǒng)與一個耦合系統(tǒng)共同組成一個擴展系統(tǒng)，允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。Q：等效質量：等效質量S: 擴展坐標變量擴展坐標變量 ： 熱力學阻尼系數(shù)熱力學阻尼系數(shù)L： 擴展系統(tǒng)的自由度擴展系統(tǒng)的自由度u Predictor-corrector algorithm is straightforwardu Verlet algorithm is feasible, but tricky to implement積分方案：積分方案：update of depends on pupdate of p depends on Nos-Hoover方法正確地描述了方法正確地描述了NVT系綜中的動系綜中的動量和位型，而等動能方法只正確地描述了后者。量和位型，而等動能方法只正確地描述了后者。擴展系統(tǒng)的哈密頓量擴展系統(tǒng)的哈密頓量Hamiltonian守恒：守恒：UFF, OPLS, Amber, CVFF, Compass分子模擬方法補充介紹：分子模擬方法補充介紹：kb is the spring constant of the bond.r0 is the bond lengt