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

分子動力學(xué)模擬方法第一頁，共七十二頁，2022年，8月28日1957年：基于剛球勢的分子動力學(xué)法（AlderandWainwright）1964年：利用Lennard-Jone勢函數(shù)法對液態(tài)氬性質(zhì)的模擬（Rahman）1971年：模擬具有分子團簇行為的水的性質(zhì)（RahmanandStillinger）1977年：約束動力學(xué)方法（Rychaert,Ciccotti&Berendsen;vanGunsteren）1980年：恒壓條件下的動力學(xué)方法（Andersen法、Parrinello-Rahman法）1983年：非平衡態(tài)動力學(xué)方法（GillanandDixon）1984年:恒溫條件下的動力學(xué)方法(Berendsenetal.)1984年：恒溫條件下的動力學(xué)方法（Nosé-Hoover法）1985年：第一原理分子動力學(xué)法（→Car-Parrinello法）1991年：巨正則系綜的分子動力學(xué)方法（CaginandPettit）分子動力學(xué)簡史第二頁，共七十二頁，2022年，8月28日粒子的運動取決于經(jīng)典力學(xué)(牛頓定律（F=ma）課程講解內(nèi)容：經(jīng)典分子動力學(xué)

(ClassicalMolecularDynamics)第三頁，共七十二頁，2022年，8月28日分子動力學(xué)方法基礎(chǔ)：原理：計算一組分子的相空間軌道，其中每個分子各自服從牛頓運動定律：初始條件：第四頁，共七十二頁，2022年，8月28日分子動力學(xué)是在原子、分子水平上求解多體問題的重要的計算機模擬方法，可以預(yù)測納米尺度上的材料動力學(xué)特性。通過求解所有粒子的運動方程，分子動力學(xué)方法可以用于模擬與原子運動路徑相關(guān)的基本過程。在分子動力學(xué)中，粒子的運動行為是通過經(jīng)典的Newton運動方程所描述。分子動力學(xué)方法是確定性方法，一旦初始構(gòu)型和速度確定了，分子隨時間所產(chǎn)生的運動軌跡也就確定了。分子動力學(xué)方法特征：第五頁，共七十二頁，2022年，8月28日分子動力學(xué)的算法：有限差分方法一、Verlet算法粒子位置的Taylor展開式：粒子位置：粒子速度：粒子加速度：開始運動時需要r(t-Δt)：+缺點：Verlet算法處理速度非常笨拙第六頁，共七十二頁，2022年，8月28日Verlet算法的表述：算法啟動規(guī)定初始位置規(guī)定初始速度擾動初始位置：計算第n步的力計算第n+1步的位置：計算第n步的速度：重復(fù)④至⑥第七頁，共七十二頁，2022年，8月28日Verlet算法程序：Do100I=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))/DT2VYI=(RYNEWI–RYOLD(I))/DT2VZI=(RZNEWI–RZOLD(I))/DT2RXOLD(I)=RX(I)RYOLD(I)=RY(I)RZOLD(I)=RZ(I)RX(I)=RXNEWIRY(I)=RYNEWIRZ(I)=RZNEWI100CONTINUE第八頁，共七十二頁，2022年，8月28日優(yōu)點：1、精確，誤差O(Δ4)2、每次積分只計算一次力3、時間可逆缺點：1、速度有較大誤差O(Δ2)2、軌跡與速度無關(guān)，無法與熱浴耦聯(lián)Verlet算法的優(yōu)缺點：第九頁，共七十二頁，2022年，8月28日二、蛙跳(Leap-frog)算法：半步算法1.首先利用當(dāng)前時刻的加速度，計算半個時間步長后的速度：2.計算下一步長時刻的位置：3.計算當(dāng)前時刻的速度：t-Δt/2tt+Δt/2t+Δtt+3Δt/2t+2Δtvrv開始運動時需要v(-Δt/2)：第十頁，共七十二頁，2022年，8月28日Leap-frog算法的表述：算法啟動規(guī)定初始位置規(guī)定初始速度擾動初始速度：計算第n步的力計算第n+1/2步的速度：計算第n+1步的位置：計算第n步的速度：重復(fù)④至⑦第十一頁，共七十二頁，2022年，8月28日Leap-frog算法的優(yōu)缺點：優(yōu)點：1、提高精確度2、軌跡與速度有關(guān)，可與熱浴耦聯(lián)缺點：1、速度近似2、比Verlet算子多花時間第十二頁，共七十二頁，2022年，8月28日三、VelocityVerlet算法：等價于優(yōu)點：速度計算更加準(zhǔn)確第十三頁，共七十二頁，2022年，8月28日VelocityVerlet算法的表述：算法啟動規(guī)定初始位置規(guī)定初始速度計算第n+1步的位置：計算第n+1步的力計算第n+1步的速度：重復(fù)③至⑤第十四頁，共七十二頁，2022年，8月28日Verlet三種形式算法的比較：VerletLeap-frogVelocityVerlet第十五頁，共七十二頁，2022年，8月28日四、預(yù)測－校正(Predictor-Corrector)格式算法：預(yù)測(Predictor)階段：其基本思想是Taylor展開，第十六頁，共七十二頁，2022年，8月28日根據(jù)新的原子位置rp，可以計算獲得校正后的ac(t+t),定義預(yù)測誤差:利用此預(yù)測誤差，對預(yù)測出的位置、速度、加速度等量進行校正：校正(Corrector)階段：第十七頁，共七十二頁，2022年，8月28日預(yù)測階段運動方程的變換：定義一組矢量:第十八頁，共七十二頁，2022年，8月28日校正階段運動方程的變換：的形式：C0，C1，C2，C3的值以及C0，取決于運動方程的階數(shù)。第十九頁，共七十二頁，2022年，8月28日一階運動方程：Valuesc0c1c2c3c4c535/1211/2

43/813/41/65251/720111/121/31/24

695/288125/2435/725/481/120第二十頁，共七十二頁，2022年，8月28日二階運動方程之一：Valuesc0c1c2c3c4c5301141/65/611/3519/1203/411/21/1263/20251/360111/181/61/60第二十一頁，共七十二頁，2022年，8月28日二階運動方程之二：Valuesc0c1c2c3c4c5301141/65/611/3519/903/411/21/1263/16251/360111/181/61/60第二十二頁，共七十二頁，2022年，8月28日五、積分時間步長t的選擇：室溫下，?t≈1fs(femtosecond10-15s)，溫度越高，?t

應(yīng)該減小太長的時間步長會造成分子間的激烈碰撞，體系數(shù)據(jù)溢出;太短的時間步長會降低模擬過程搜索相空間的能力第二十三頁，共七十二頁，2022年，8月28日微正則系綜分子動力學(xué)(NVEMD)它是分子動力學(xué)方法的最基本系綜具有確定的粒子數(shù)N，能量E和體積V算法：規(guī)定初始位置和初始速度對運動方程積分若干步計算勢能和動能若能量不等于所需要的值，對速度進行標(biāo)度重復(fù)②至④，直到系統(tǒng)平衡第二十四頁，共七十二頁，2022年，8月28日微正則系綜(NVE)MD模擬算法的流程圖：給定每個分子的初始位置ri(0)和速度vi(0)計算每個分子的受力Fi和加速度ai解運動方程并求出每個分子運動一個時間步長后到達(dá)的位置所具有的速度統(tǒng)計系統(tǒng)的熱力學(xué)性質(zhì)及其它物理量統(tǒng)計性質(zhì)不變？

打印結(jié)果，結(jié)束YesNo移動所有分子到新的位置并具有當(dāng)前時刻的速度第二十五頁，共七十二頁，2022年，8月28日微正則系綜MD模擬程序F3講解（LJ,NVE）：無因次量：第二十六頁，共七十二頁，2022年，8月28日MD模擬中幾個熱力學(xué)量的計算：對于由N個單原子組成的系統(tǒng)：動能和溫度：采用對比量：第二十七頁，共七十二頁，2022年，8月28日對于LJ流體：勢能：采用對比量：第二十八頁，共七十二頁，2022年，8月28日內(nèi)能：內(nèi)能由勢能和動能組成：采用對比量：第二十九頁，共七十二頁，2022年，8月28日壓力：采用對比量：第三十頁，共七十二頁，2022年，8月28日第三十一頁，共七十二頁，2022年，8月28日練習(xí)：推導(dǎo)LJ流體分子間力的表達(dá)式(fx,fy,fz及其對比量)：勢能函數(shù)形式：力：采用對比量：=x,y,z第三十二頁，共七十二頁，2022年，8月28日LJ分子間的維里項：=x,y,z第三十三頁，共七十二頁，2022年，8月28日采用對比量的運動方程形式：

(以蛙跳(Leap-frog)算法為例)采用對比量：第三十四頁，共七十二頁，2022年，8月28日最終得到：同理得到：第三十五頁，共七十二頁，2022年，8月28日速度的標(biāo)度(VelocityScaling)：根據(jù)能量均分原理，可知：標(biāo)度因子:對比量速度標(biāo)度:或第三十六頁，共七十二頁，2022年，8月28日微正則系綜MD模擬程序F3講解（LJ,NVE）：初始化：READ(*,‘(A)’)TITLE!運行作業(yè)題目READ(*,*)NSTEP!運行步數(shù)READ(*,*)IPRINT!打印步數(shù)READ(*,‘(A)’)CNFILE!位型文件READ(*,*)DENS!對比密度READ(*,*)RTEMP

!對比溫度

READ(*,*)RCUT!對比截斷半徑READ(*,*)DT!對比時間步長CALLREADCN(CNFILE)第三十七頁，共七十二頁，2022年，8月28日初始位型：面心立方(face-centeredcubic,FCC)：每面中心有一格點

第三十八頁，共七十二頁，2022年，8月28日體心立方(body-centeredcubic,BCC)：簡單立方(simplecubic,SC)：第三十九頁，共七十二頁，2022年，8月28日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)R(4)=(Y,Y,0)M=0DO10I=1,NCDO10J=1,NCDO10K=1,NCDO11IJ=1,4RX(IJ+M)=RX(IJ)+XL*(K-1)RY(IJ+M)=RY(IJ)+XL*(J-1)RZ(IJ+M)=RZ(IJ)+XL*(I-1)11CONTINUEM=M+410CONTINUE第四十頁，共七十二頁，2022年，8月28日DO100I=1,NRX(I)=RX(I)-0.5RY(I)=RY(I)-0.5RZ(I)=RZ(I)-0.5100CONTINUE將模擬盒子的中心移到原點：第四十一頁，共七十二頁，2022年，8月28日初始速度：簡單的選擇：V＝random(-0.5,0.5)=x,y,z標(biāo)度因子:速度標(biāo)度:第四十二頁，共七十二頁，2022年，8月28日FACTOR=SQRT(RTEMP)DO100I=1,NVX(I)=FACTOR

*(RANF(DUMMY)-0.5)VY(I)=FACTOR

*(RANF(DUMMY)-0.5)VZ(I)=FACTOR

*(RANF(DUMMY)-0.5)CONTINUE隨機安排初始速度：第四十三頁，共七十二頁，2022年，8月28日標(biāo)度初始速度：SUMKX=0.0SUMKY=0.0SUMKZ=0.0DO200I=1,NSUMKX=SUMKX+VX(I)**2SUMKY=SUMKY+VY(I)**2SUMKZ=SUMKZ+VZ(I)**2200CONTINUEBEITAX=SQRT(RTEMP/SUMKX)BEITAY=SQRT(RTEMP/SUMKY)

BEITAZ

=SQRT(RTEMP/SUMKZ)DO300I=1,NVX(I)=VX(I)*BEITAXVY(I)=VY(I)*BEITAYVZ(I)=VZ(I)*BEITAZ300CONTINUE標(biāo)度因子:第四十四頁，共七十二頁，2022年，8月28日SUMX=0.0SUMY=0.0SUMZ=0.0DO200I=1,NSUMX=SUMX+VX(I)SUMY=SUMY+VY(I)SUMZ=SUMZ+VZ(I)CONTINUESUMX=SUMX/REAL(N)SUMY=SUMY/REAL(N)SUMZ=SUMZ/REAL(N)DO300I=1,NVX(I)=VX(I)-SUMXVY(I)=VY(I)-SUMYVZ(I)=VZ(I)-SUMZ300CONTINUE控制體系的總動量為零：第四十五頁，共七十二頁，2022年，8月28日從Maxwell分布中抽樣：xxdx高斯(Gauss)分布:對于等幾率隨機試驗(Bernoulli試驗)，大量的試驗結(jié)果滿足高斯分布第四十六頁，共七十二頁，2022年，8月28日麥克斯韋速度分布定律:由于：=x,y,z單位體積的分子再每個分量上的速度分布實際上就是高斯分布。第四十七頁，共七十二頁，2022年，8月28日從Maxwell分布中抽樣：高斯(Gauss)分布的隨機數(shù)生成方法:生成隨機數(shù)：i，i=1,2,…,12第四十八頁，共七十二頁，2022年，8月28日SUM=0.0DO10I=1,12SUM=SUM+RANF(DUMMY)10CONTINUER=(SUM-6.0)/4.0R2=R*RGAUSS=((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R高斯(Gauss)分布的隨機數(shù)生成(程序F24)第四十九頁，共七十二頁，2022年，8月28日FACTOR=SQRT(RTEMP)DO100I=1,NVX(I)=FACTOR

*GAUSS(DUMMY)VY(I)=FACTOR

*GAUSS(DUMMY)VZ(I)=FACTOR

*GAUSS(DUMMY)CONTINUE控制總動量為零：同前面一樣處理。從Maxwell分布中抽樣分布中隨機安排初始速度：第五十頁，共七十二頁，2022年，8月28日微正則系綜MD模擬程序F3講解(LJ,NVE)

：量綱變換：SIGMA=(DENS/REAL(N))**(1.0/3.0)RCUT=RCUT*SIGMADT＝DT*

SIGMADENS=DENS/(SIGMA**3)模擬盒子的邊長為1第五十一頁，共七十二頁，2022年，8月28日長程校正：微正則系綜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.0/9.0)*PI*DENS*SIGCUB*REAL(N): *(2.0*SR9-3.0*SR3)第五十二頁，共七十二頁，2022年，8月28日算法：算法啟動微正則系綜MD模擬程序F3講解(LJ,NVE)

：CALLFORCE(-DT,SIGMA,RCUT,NEWV,NEWVC,NEWW)CALLMOVE(-DT)CALLFORCE(-DT,SIGMA,RCUT,V,VC,W)CALLFORCE(DT,SIGMA,RCUT,V,VC,W)CALLKINET(OLDK)CALLMOVE(DT)CALLFORCE(DT,SIGMA,RCUT,NEWV,NEWVC,NEWW)CALLKINET(NEWK)第五十三頁，共七十二頁，2022年，8月28日算法：差分格式：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)–DVX………V=V+VIJW=W+WIJCALLMOVE(DT)第五十四頁，共七十二頁，2022年，8月28日DO1000I=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):第五十五頁，共七十二頁，2022年，8月28日速度的標(biāo)定（只用于平衡階段）SUMK=0.0DO200I=1,NSUMK=SUMKX+VX(I)**2+VY(I)**2+VZ(I)**2200CONTINUEBEITA=SQRT(3.0*RTEMP/SUMK)DO300I=1,NVX(I)=VX(I)*BEITAVY(I)=VY(I)*BEITAVZ(I)=VZ(I)*BEITA300CONTINUE第五十六頁，共七十二頁，2022年，8月28日正則系綜分子動力學(xué)(NVTMD)具有確定的粒子數(shù)N，溫度T和體積V速度的直接標(biāo)度熱浴方法(AndersenThermostat)約束方法(阻尼力方法)系統(tǒng)擴展方法(ExtendedSystemsMethod)問題的關(guān)鍵：溫度的約束Nose-Hoover方法第五十七頁，共七十二頁，2022年，8月28日一、熱浴方法(AndersenThermostat)引入一個與虛擬粒子碰撞的隨機力想象系統(tǒng)浸在熱浴當(dāng)中系統(tǒng)和熱浴間的相互作用強度由隨機碰撞的頻率決定碰撞的幾率等于Nu×dt如果一個粒子經(jīng)歷碰撞，它的速度將從約束溫度下的Maxwell分布中隨機抽取總能量和總動量均不守恒第五十八頁，共七十二頁，2022年，8月28日二、約束方法是等動能(Iso-Kinetics)分子動力學(xué)方法系統(tǒng)的運動方程為：引入阻尼系數(shù)以保證將溫度約束在恒定值根據(jù)高斯最小約束原理：第五十九頁，共七十二頁，2022年，8月28日三、Nose-Hoover擴展方法基本思想：設(shè)想原系統(tǒng)與一個耦合系統(tǒng)共同組成一個擴展系統(tǒng)，允許熱流在原系統(tǒng)和耦合系統(tǒng)之間交換。Q：等效質(zhì)量S:擴展坐標(biāo)變量：熱力學(xué)阻尼系數(shù)L：擴展系統(tǒng)的自由度第六十頁，共七十二頁，2022年，8月28日Predictor-correctoralgorithmisstraightforwardVerletalgorithmisfeasible,buttrickytoimplement積分方案：updateofxdependsonpupdateofpdependsonx第六十一頁，共七十二頁，2022年，8月28日Nosé-Hoover方法正確地描述了NVT系綜中的動量和位型，而等動能方法只正確地描述了后者。擴展系統(tǒng)的哈密頓量Hamiltonian守恒：第六十二頁，共七十二頁，2022年，8月28日1.復(fù)雜分子體系的勢能函數(shù)形式：PotentialEnergy=StretchingEnergy+BendingEnergy+TorsionEnergy+ Non-BondedInteractionEnergy這些方程與描述原子或鍵各種不同行為的參數(shù)就構(gòu)成了力場，

force-field.UFF,OPLS,Amber,CVFF,Compass分子模擬方法補充介紹：第六十三頁，共七十二頁，2022年，8月28日BondStretchingEnergykbisthespringconstantofthebond.r0isthebondlengthatequilibrium.Uniquekbandr0assignedforeachbondpair,i.e.C-C,O-H第六十四頁，共七十二頁，2022年，8月28日BendingEnergykisthespringconstantofthebend.0isthebondlengthatequilibrium.Uniqueparametersforanglebendingareassignedtoeachbondedtripletofatomsbasedontheirtypes(e.g.C-C-C,C-O-C,C-C-H,etc.)第六十五頁，共七十二頁，2022年，8月28日The“Hookeian”potentialkbandkbroadenorsteepentheslopeoftheparabola.Thelargerthevalueofk,themoreenergyisrequiredtodeformanangle(orbond)fromitsequilibriumvalue.第六十六頁，共七十二頁，2022年，8月28日TorsionEnergyAcontrolstheamplitudeofthecurvencontrolsitsperiodicityshiftstheentirecurvealongtherotationangleaxis().Theparametersaredeterminedfromcurvefitting.Uniqueparametersfortorsionalrotationareassignedtoeachbondedquartetofatomsbasedontheirtypes(e.