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

1、分子動力學(xué)模擬方法第1頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一1957年：基于剛球勢的分子動力學(xué)法（Alder and Wainwright） 1964年：利用Lennard-Jone勢函數(shù)法對液態(tài)氬性質(zhì)的模擬（Rahman） 1971年：模擬具有分子團(tuán)簇行為的水的性質(zhì)（Rahman and Stillinger） 1977年：約束動力學(xué)方法（Rychaert, Ciccotti & Berendsen; van Gunsteren） 1980年：恒壓條件下的動力學(xué)方法（Andersen法、Parrinello-Rahman法） 1983年：非平衡態(tài)動力學(xué)方法（Gillan a

2、nd Dixon） 1984年: 恒溫條件下的動力學(xué)方法(Berendsen et al.) 1984年：恒溫條件下的動力學(xué)方法（Nos-Hoover法） 1985年：第一原理分子動力學(xué)法（Car-Parrinello法） 1991年：巨正則系綜的分子動力學(xué)方法（Cagin and Pettit）分子動力學(xué)簡史第2頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一粒子的運(yùn)動取決于經(jīng)典力學(xué)(牛頓定律（F=ma）課程講解內(nèi)容：經(jīng)典分子動力學(xué)(Classical Molecular Dynamics)第3頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一分子動力學(xué)方法基礎(chǔ)：原理：計算一組

3、分子的相空間軌道，其中每個分子各自服從牛頓運(yùn)動定律：初始條件：第4頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一分子動力學(xué)是在原子、分子水平上求解多體問題的重要的計算機(jī)模擬方法，可以預(yù)測納米尺度上的材料動力學(xué)特性。通過求解所有粒子的運(yùn)動方程，分子動力學(xué)方法可以用于模擬與原子運(yùn)動路徑相關(guān)的基本過程。在分子動力學(xué)中，粒子的運(yùn)動行為是通過經(jīng)典的Newton運(yùn)動方程所描述。分子動力學(xué)方法是確定性方法，一旦初始構(gòu)型和速度確定了，分子隨時間所產(chǎn)生的運(yùn)動軌跡也就確定了。分子動力學(xué)方法特征：第5頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一分子動力學(xué)的算法：有限差分方法一、Verlet算法

4、粒子位置的Taylor展開式：粒子位置：粒子速度：粒子加速度：開始運(yùn)動時需要r(t-t)：+缺點(diǎn)：Verlet算法處理速度非常笨拙第6頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Verlet算法的表述： 算法啟動規(guī)定初始位置規(guī)定初始速度擾動初始位置：計算第n步的力計算第n+1步的位置：計算第n步的速度：重復(fù)至第7頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Verlet算法程序：Do 100 I = 1, N RXNEWI = 2.0 * RX(I) RXOLD(I) + DTSQ * AX(I) RYNEWI = 2.0 * RY(I) RYOLD(I) + DTSQ

5、* 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(I) RZOLD(I) = RZ(I) RX(I) = RXNEWI RY(I) = RYNEWI RZ(I) = RZNEWI100 CONTINUE第8頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一優(yōu)點(diǎn)：1、精確，誤差O(4)

6、2、每次積分只計算一次力3、時間可逆缺點(diǎn)：1、速度有較大誤差O(2)2、軌跡與速度無關(guān)，無法與熱浴耦聯(lián)Verlet算法的優(yōu)缺點(diǎn)：第9頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一二、蛙跳(Leap-frog)算法：半步算法1. 首先利用當(dāng)前時刻的加速度，計算半個時間步長后的速度：2. 計算下一步長時刻的位置：3. 計算當(dāng)前時刻的速度：t-t/2tt+t/2t+tt+3t/2t+2tvrv開始運(yùn)動時需要v(-t/2)：第10頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Leap-frog算法的表述： 算法啟動規(guī)定初始位置規(guī)定初始速度擾動初始速度：計算第n步的力計算第n+1/2

7、步的速度：計算第n+1步的位置：計算第n步的速度：重復(fù)至第11頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Leap-frog算法的優(yōu)缺點(diǎn)：優(yōu)點(diǎn)：1、提高精確度2、軌跡與速度有關(guān)，可與熱浴耦聯(lián)缺點(diǎn)：1、速度近似2、比Verlet算子多花時間第12頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一三、Velocity Verlet算法：等價于優(yōu)點(diǎn)：速度計算更加準(zhǔn)確第13頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Velocity Verlet算法的表述： 算法啟動規(guī)定初始位置規(guī)定初始速度計算第n+1步的位置：計算第n+1步的力計算第n+1步的速度：重復(fù)至第14頁，共7

8、2頁，2022年，5月20日，9點(diǎn)51分，星期一Verlet三種形式算法的比較：VerletLeap-frogVelocity Verlet第15頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一四、預(yù)測校正(Predictor-Corrector)格式算法：預(yù)測(Predictor)階段：其基本思想是Taylor展開，第16頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一根據(jù)新的原子位置rp，可以計算獲得校正后的ac(t+t),定義預(yù)測誤差:利用此預(yù)測誤差，對預(yù)測出的位置、速度、加速度等量進(jìn)行校正：校正(Corrector)階段：第17頁，共72頁，2022年，5月20日，9點(diǎn)

9、51分，星期一預(yù)測階段運(yùn)動方程的變換：定義一組矢量:第18頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一校正階段運(yùn)動方程的變換：的形式：C0， C1， C2， C3的值以及C0，取決于運(yùn)動方程的階數(shù)。第19頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一 一階運(yùn)動方程：Valuesc0c1c2c3c4c535/1211/243/813/41/65251/720111/121/31/24695/288125/2435/725/481/120第20頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一 二階運(yùn)動方程之一：Valuesc0c1c2c3c4c5301141/65/6

10、11/3519/1203/411/21/1263/20251/360111/181/61/60第21頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一 二階運(yùn)動方程之二：Valuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60第22頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一五、積分時間步長t的選擇：室溫下， t 1 fs (femtosecond 10-15s)，溫度越高，t 應(yīng)該減小 太長的時間步長會造成分子間的激烈碰撞，體系數(shù)據(jù)溢出; 太短的時間步長會降低模擬過程搜索相空間的能

11、力 第23頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一微正則系綜分子動力學(xué)(NVE MD) 它是分子動力學(xué)方法的最基本系綜 具有確定的粒子數(shù)N，能量E和體積V 算法：規(guī)定初始位置和初始速度對運(yùn)動方程積分若干步計算勢能和動能若能量不等于所需要的值，對速度進(jìn)行標(biāo)度重復(fù)至，直到系統(tǒng)平衡第24頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一微正則系綜(NVE)MD模擬算法的流程圖：給定每個分子的初始位置ri(0)和速度vi(0)計算每個分子的受力Fi和加速度ai解運(yùn)動方程并求出每個分子運(yùn)動一個時間步長后到達(dá)的位置所具有的速度統(tǒng)計系統(tǒng)的熱力學(xué)性質(zhì)及其它物理量統(tǒng)計性質(zhì)不變？ 打印結(jié)果，

12、結(jié)束YesNo移動所有分子到新的位置并具有當(dāng)前時刻的速度第25頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一微正則系綜MD模擬程序F3講解（LJ, NVE）：無因次量：第26頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一MD模擬中幾個熱力學(xué)量的計算：對于由N個單原子組成的系統(tǒng)：動能和溫度：采用對比量：第27頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一對于LJ流體：勢能：采用對比量：第28頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一內(nèi)能：內(nèi)能由勢能和動能組成：采用對比量：第29頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一壓力：采用對比量：第

13、30頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一第31頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一練習(xí)：推導(dǎo)LJ流體分子間力的表達(dá)式(fx, fy, fz及其對比量)：勢能函數(shù)形式：力：采用對比量：=x, y, z第32頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一LJ分子間的維里項：=x, y, z第33頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一采用對比量的運(yùn)動方程形式：(以蛙跳(Leap-frog)算法為例)采用對比量：第34頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一最終得到：同理得到：第35頁，共72頁，2022年，5月20日

14、，9點(diǎn)51分，星期一速度的標(biāo)度(Velocity Scaling)：根據(jù)能量均分原理，可知：標(biāo)度因子:對比量速度標(biāo)度:或第36頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一微正則系綜MD模擬程序F3講解（LJ, NVE）：初始化：READ (*,(A) TITLE ! 運(yùn)行作業(yè)題目READ (*,*) NSTEP ! 運(yùn)行步數(shù)READ (*,*) IPRINT ! 打印步數(shù)READ (*,(A) CNFILE ! 位型文件READ (*,*) DENS ! 對比密度READ (*,*) RTEMP ! 對比溫度 READ (*,*) RCUT ! 對比截斷半徑READ (*,*) D

15、T ! 對比時間步長CALL READCN ( CNFILE )第37頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一初始位型：面心立方 (face-centered cubic, FCC)：每面中心有一格點(diǎn) 第38頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一體心立方 (body-centered cubic, BCC)：簡單立方 (simple cubic, SC)：第39頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一XL初始位型：面心立方 ( FCC) (程序F23)NC=(REAL(N)/4.0)*(1.0/3.0)XL = 1.0 / REAL ( NC

16、)Y = 0.5 * XLR(1)=(0, 0, 0) R(2)=(0, Y, Y)R(3)=(Y, 0, Y) 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 CONTINUE第40頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一DO 100 I = 1, N

17、RX(I) = RX(I) - 0.5 RY(I) = RY(I) - 0.5 RZ(I) = RZ(I) - 0.5100 CONTINUE將模擬盒子的中心移到原點(diǎn)：第41頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一初始速度：簡單的選擇：V random(-0.5, 0.5)=x, y, z標(biāo)度因子:速度標(biāo)度:第42頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一FACTOR = SQRT ( RTEMP )DO 100 I = 1, N VX(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) VY(I) = FACTOR * ( RANF(DUMM

18、Y) - 0.5 ) VZ(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) CONTINUE隨機(jī)安排初始速度：第43頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一標(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

19、(RTEMP/SUMKZ)DO 300 I = 1, N VX(I) = VX(I) * BEITAX VY(I) = VY(I) * BEITAY VZ(I) = VZ(I) * BEITAZ300 CONTINUE標(biāo)度因子:第44頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一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 / REA

20、L ( N )SUMZ = SUMZ / REAL ( N ) DO 300 I = 1, N VX(I) = VX(I) - SUMX VY(I) = VY(I) - SUMY VZ(I) = VZ(I) - SUMZ300 CONTINUE控制體系的總動量為零：第45頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一從Maxwell分布中抽樣：xxdx高斯(Gauss)分布:對于等幾率隨機(jī)試驗(Bernoulli試驗)，大量的試驗結(jié)果滿足高斯分布第46頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一麥克斯韋速度分布定律:由于：=x, y, z單位體積的分子再每個分量上的速度分

21、布實際上就是高斯分布。第47頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一從Maxwell分布中抽樣：高斯(Gauss)分布的隨機(jī)數(shù)生成方法:生成隨機(jī)數(shù)：i，i=1, 2, , 12第48頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一SUM = 0.0DO 10 I = 1, 12 SUM = SUM + RANF ( DUMMY )10 CONTINUER = ( SUM - 6.0 ) / 4.0R2 = R * RGAUSS = ( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 ) * R高斯(Gauss)分布的隨機(jī)數(shù)

22、生成(程序F24)第49頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一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分布中抽樣分布中隨機(jī)安排初始速度：第50頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一微正則系綜MD模擬程序F3講解(LJ, NVE) ：量綱變換：SIGMA = ( DENS / REA

23、L ( N ) ) * ( 1.0 / 3.0 )RCUT = RCUT * SIGMADT DT * SIGMADENS = DENS / ( SIGMA * 3 )模擬盒子的邊長為1第51頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一長程校正：微正則系綜MD模擬程序F3講解(LJ, NVE) ：SR3 = ( SIGMA / RCUT ) * 3SR9 = SR3 * 3SIGCUB = SIGMA * 3VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N ) : * ( SR9 - 3.0 * SR3 )WLRC = ( 16

24、.0 / 9.0 ) * PI * DENS * SIGCUB * REAL ( N ) :* ( 2.0 * SR9 - 3.0 * SR3 )第52頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一算法：算法啟動微正則系綜MD模擬程序F3講解(LJ, NVE) ：CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )CALL MOVE ( -DT )CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W )CALL FORCE ( DT, SIGMA, RCUT, V, VC, W )CALL KINET ( OL

25、DK )CALL MOVE ( DT )CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW )CALL KINET ( NEWK )第53頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一算法：差分格式：SR2 = SIGSQ / RIJSQVIJ = 4.0 * ( SR12 - SR6 )WIJ = 24.0 * ( 2.0 * SR12 - SR6 )VELIJ = WIJ * DT / RIJSQDVX = VELIJ * RXIJ. VXI = VXI + DVX.VX(J) = VX(J) DVXV = V + VIJW = W +

26、 WIJCALL MOVE(DT)第54頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一DO 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 CONTINUEMOVE(DT):第55頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一速度的標(biāo)定（只用于平衡階段）SUMK = 0.0DO 200 I = 1, N SUMK = SUMKX + VX(I)*2+VY(I)*2 + VZ(I)*2200 CONTINUEBEITA =SQR

27、T( 3.0 * RTEMP / SUMK )DO 300 I = 1, N VX(I) = VX(I) * BEITA VY(I) = VY(I) * BEITA VZ(I) = VZ(I) * BEITA300 CONTINUE第56頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一正則系綜分子動力學(xué)(NVT MD) 具有確定的粒子數(shù)N，溫度T和體積V速度的直接標(biāo)度熱浴方法 (Andersen Thermostat)約束方法(阻尼力方法)系統(tǒng)擴(kuò)展方法(Extended Systems Method) 問題的關(guān)鍵：溫度的約束Nose-Hoover方法第57頁，共72頁，2022年，5月

28、20日，9點(diǎn)51分，星期一一、熱浴方法 (Andersen Thermostat) 引入一個與虛擬粒子碰撞的隨機(jī)力 想象系統(tǒng)浸在熱浴當(dāng)中系統(tǒng)和熱浴間的相互作用強(qiáng)度由隨機(jī)碰撞的頻率決定碰撞的幾率等于Nudt如果一個粒子經(jīng)歷碰撞，它的速度將從約束溫度下的Maxwell分布中隨機(jī)抽取 總能量和總動量均不守恒第58頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一二、約束方法 是等動能(Iso-Kinetics)分子動力學(xué)方法 系統(tǒng)的運(yùn)動方程為： 引入阻尼系數(shù)以保證將溫度約束在恒定值 根據(jù)高斯最小約束原理：第59頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一三、Nose-Hoover擴(kuò)

29、展方法基本思想：設(shè)想原系統(tǒng)與一個耦合系統(tǒng)共同組成一個擴(kuò)展系統(tǒng)，允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。Q：等效質(zhì)量S: 擴(kuò)展坐標(biāo)變量： 熱力學(xué)阻尼系數(shù)L： 擴(kuò)展系統(tǒng)的自由度第60頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一 Predictor-corrector algorithm is straightforward Verlet algorithm is feasible, but tricky to implement積分方案：update of x depends on pupdate of p depends on x第61頁，共72頁，2022年，5月20日，9點(diǎn)51分，星

30、期一Nos-Hoover方法正確地描述了NVT系綜中的動量和位型，而等動能方法只正確地描述了后者。擴(kuò)展系統(tǒng)的哈密頓量Hamiltonian守恒：第62頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一1. 復(fù)雜分子體系的勢能函數(shù)形式：Potential Energy = Stretching Energy + Bending Energy + Torsion Energy + Non-Bonded Interaction Energy這些方程與描述原子或鍵各種不同行為的參數(shù)就構(gòu)成了力場， force-field.UFF, OPLS, Amber, CVFF, Compass分子模擬方法補(bǔ)

31、充介紹：第63頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Bond Stretching Energykb is the spring constant of the bond.r0 is the bond length at equilibrium.Unique kb and r0 assigned for each bond pair, i.e. C-C, O-H第64頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Bending Energyk is the spring constant of the bend.0 is the bond length at equ

32、ilibrium.Unique parameters for angle bending are assigned to each bonded triplet of atoms based on their types (e.g. C-C-C, C-O-C, C-C-H, etc.)第65頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一The “Hookeian” potentialkb and k broaden or steepen the slope of the parabola. The larger the value of k, the more energy is required to deform an angle (or bond) from its equilibrium value. 第66頁，共72頁，2022年，5月20日，9點(diǎn)51分，星期一Torsion EnergyA controls the amplitude of the curven controls its periodicity shifts the entire curve along the rotation angle axis (). The parameters are determined from curve fit