分子動(dòng)力學(xué)模擬方法

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

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

3、學(xué)(牛頓定律（牛頓定律（F=ma）原理：原理：計(jì)算一組分子的相空間軌道，其中每個(gè)分子各自服從計(jì)算一組分子的相空間軌道，其中每個(gè)分子各自服從牛頓運(yùn)動(dòng)定律：牛頓運(yùn)動(dòng)定律： 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開始運(yùn)動(dòng)時(shí)需要開始運(yùn)動(dòng)時(shí)需要r(t-t)：t (0) (0) t)(ivrr+n 算法啟動(dòng)算法啟動(dòng)u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 擾動(dòng)初始位置：擾動(dòng)初始位置：u 計(jì)算第計(jì)算第n步的力步的力u 計(jì)算第計(jì)算第n+1步的位置：步的位置：u 計(jì)算第計(jì)算第n步的速度：步的速度：u 重復(fù)重復(fù)至至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. 首先利用當(dāng)前時(shí)刻的加速度，計(jì)算半個(gè)時(shí)間步長后的速度：首先利用當(dāng)前時(shí)刻的加速度，計(jì)算半個(gè)時(shí)間步長后的速度：2. 計(jì)算下一步長時(shí)刻的位置：計(jì)算下一步長時(shí)刻的位置：3. 計(jì)算當(dāng)前時(shí)刻的速度：計(jì)算當(dāng)前時(shí)刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv開始運(yùn)動(dòng)時(shí)需要開始運(yùn)動(dòng)時(shí)需要v(-t/

7、2)：t/2 (0) (0) t/2)(iavvn 算法啟動(dòng)算法啟動(dòng)u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 擾動(dòng)初始速度：擾動(dòng)初始速度：u 計(jì)算第計(jì)算第n步的力步的力u 計(jì)算第計(jì)算第n+1/2步的速度：步的速度：u 計(jì)算第計(jì)算第n+1步的位置：步的位置：u 計(jì)算第計(jì)算第n步的速度：步的速度：u 重復(fù)重復(fù)至至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等價(jià)于等價(jià)于n 算法啟動(dòng)算法啟動(dòng)u 規(guī)定初始位置規(guī)定初始位置u 規(guī)定初始速度規(guī)定初始速度u 計(jì)算第計(jì)算第n+1步的位置：步的位置：u 計(jì)算第計(jì)算第n+1步的力步的力u 計(jì)算第計(jì)算第n+1步的速度：步的速度：u 重復(fù)重復(fù)至至2iit (t) 21 t (t) (t) t)(tavrriitt)(t (t) 21 (t) t)(tiiaavviiVerletLeap-frogVelocity Verlet 預(yù)測預(yù)測(P

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

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

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

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

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

14、綜微正則系綜MD模擬程序模擬程序F3講解（講解（LJ, NVE）：）：無因次量：無因次量：/*rr 3*/*EE /*ff/*3PP /*kTT 2/12*)/(mtt 2/1*)/(vvmMD模擬中幾個(gè)熱力學(xué)量的計(jì)算：模擬中幾個(gè)熱力學(xué)量的計(jì)算：對于由對于由N個(gè)單原子組成的系統(tǒng)：個(gè)單原子組成的系統(tǒng)： 動(dòng)能和溫度：動(dò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 內(nèi)能：內(nèi)能：KEUE內(nèi)能由勢能和動(dòng)能組成：內(nèi)能由勢能和動(dòng)能組成：采用對比量：采用對比量：*KEUE 壓力：壓力：393232316 cclrcrrPlrccBPVWTkP/*rr 3*/*EE 采用對比量：采用對比量：/*3PP /*ff 111)(31NiNijijccrWWijijijijijijcrUrWrrfr)(31 )(練習(xí)：練習(xí)：推導(dǎo)推導(dǎo)LJ流體分子間力的表達(dá)式流體分子間力的表達(dá)式(fx, fy, fz及其對比量及其對比量)：

16、4)(612rrrU勢能函數(shù)形式：勢能函數(shù)形式：力：力：r)()(rUrf2612 2184 rrrrrrrUrUf/*rr 采用對比量：采用對比量：/*ff =x, y, z2*6*12* 121184 rrrrf* )(fr rWcLJ分子間的維里項(xiàng)：分子間的維里項(xiàng)：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標(biāo)度因子標(biāo)度因子:對比量對

18、比量速度標(biāo)度速度標(biāo)度:oldnewvv *old*newvv 或或微正則系綜微正則系綜MD模擬程序模擬程序F3講解（講解（LJ, NVE）：）：初始化：初始化：READ (*,(A) TITLE ! 運(yùn)行作業(yè)題目運(yùn)行作業(yè)題目READ (*,*) NSTEP ! 運(yùn)行步數(shù)運(yùn)行步數(shù)READ (*,*) IPRINT ! 打印步數(shù)打印步數(shù)READ (*,(A) CNFILE ! 位型文件位型文件READ (*,*) DENS ! 對比密度對比密度READ (*,*) RTEMP ! 對比溫度對比溫度 READ (*,*) RCUT ! 對比截?cái)喟霃綄Ρ冉財(cái)喟霃絉EAD (*,*) DT ! 對比時(shí)間

19、步長對比時(shí)間步長CALL READCN ( CNFILE ) 初始位型：初始位型：u面心立方面心立方 (face-centered cubic, FCC)：每面中心有一格點(diǎn)每面中心有一格點(diǎn) 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將模擬盒子的中心移到原點(diǎn)：將模

21、擬盒子的中心移到原點(diǎn)： 初始速度：初始速度：n 簡單的選擇：簡單的選擇：V random(-0.5, 0.5) =x, y, zmTkvB22/1*2/12*2/122/1TTvTvmTkTTreqreqrefBref標(biāo)度因子標(biāo)度因子:*old*newvv 速度標(biāo)度速度標(biāo)度:*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隨機(jī)安排初始速度：隨機(jī)安排初始速度：標(biāo)度初始速度：標(biāo)度初始速度：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標(biāo)度因子標(biāo)度因子:*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控制體系的總動(dòng)量為零：控制體系的總動(dòng)量為零：n 從從Maxwell分布中抽樣：分布中抽樣：x)(xP0 x xdx高斯高斯(Gauss)分布分布:對于等幾率隨機(jī)試驗(yàn)對于等幾率隨機(jī)試驗(yàn)(Bernoulli試驗(yàn)試驗(yàn))，大量的試驗(yàn)結(jié)果滿足高斯分布大量的試驗(yàn)結(jié)果滿足高斯分布 222)(2exp)(21xxxxxP麥克斯韋速度分布定律麥克斯韋速度分布定律: dvekTmNdvvfTkmvB22122mTkvB2由于：由于： 222)(2exp)(21xxxxxP =x,

25、y, z單位體積的分子再每個(gè)分量上的速度分布實(shí)際上就是單位體積的分子再每個(gè)分量上的速度分布實(shí)際上就是高斯分布。高斯分布。n 從從Maxwell分布中抽樣：分布中抽樣：高斯高斯(Gauss)分布的隨機(jī)數(shù)生成方法分布的隨機(jī)數(shù)生成方法: 222)(2exp)(21xxxxxP生成隨機(jī)數(shù)：生成隨機(jī)數(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)分布的隨機(jī)數(shù)生成分布的隨機(jī)數(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控制總動(dòng)量為零：同前面一樣處理?？刂瓶倓?dòng)量為零：同前面一樣處理。從從Maxwell分布中抽樣分布中隨機(jī)安排初始速度：分布中抽樣

27、分布中隨機(jī)安排初始速度：微正則系綜微正則系綜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 算法：算法啟動(dòng)算法：算法啟動(dòng)微正則系綜微正則系綜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速度的標(biāo)定（只用于平衡階段）速度的標(biāo)定（只用于平衡階段）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速度的直接標(biāo)度速度的直接標(biāo)度u熱浴方法熱浴方法 (Andersen Thermostat)u約束方法約束方法(阻尼力方法阻尼力方法)u系統(tǒng)擴(kuò)展方法系統(tǒng)擴(kuò)展方法(Extended Systems Method)u 問題的關(guān)鍵：溫度的約束問題的關(guān)鍵：溫度的約束Nose-Hoover方法方法一、熱浴方法一、熱浴方法 (Andersen Thermostat)u 引入一個(gè)與虛擬粒子碰撞的隨機(jī)力引入一個(gè)與虛擬粒子碰撞的隨機(jī)力u 想象系統(tǒng)浸在熱浴當(dāng)中想象系統(tǒng)浸在熱浴當(dāng)中u系統(tǒng)和熱浴間的相互作用強(qiáng)度由隨機(jī)碰撞的頻率決定系統(tǒng)和熱浴

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

34、展方法擴(kuò)展方法基本思想：基本思想：設(shè)想原系統(tǒng)與一個(gè)耦合系統(tǒng)共同組成一個(gè)擴(kuò)展系統(tǒng)，設(shè)想原系統(tǒng)與一個(gè)耦合系統(tǒng)共同組成一個(gè)擴(kuò)展系統(tǒng)，允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。Q：等效質(zhì)量：等效質(zhì)量S: 擴(kuò)展坐標(biāo)變量擴(kuò)展坐標(biāo)變量 ： 熱力學(xué)阻尼系數(shù)熱力學(xué)阻尼系數(shù)L： 擴(kuò)展系統(tǒng)的自由度擴(kuò)展系統(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系綜中的動(dòng)系綜中的動(dòng)量和位型，而等動(dòng)能方法只正確地描述了后者。量和位型，而等動(dòng)能方法只正確地描述了后者。擴(kuò)展系統(tǒng)的哈密頓量擴(kuò)展系統(tǒng)的哈密頓量Hamiltonian守恒：守恒：UFF, OPLS, Amber, CVFF, Compass分子模擬方法補(bǔ)充介紹：分子模擬方法補(bǔ)充介紹：kb is the spring constant of the bond.r0 is the bond lengt