分子动力学课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,计算机分子模拟方法,第三章、分子动力学方法,本章的具体内容,1分子动力学方法根底,周期边界,最小影像,差分格式,位势截断,2L-J流体体系,3微正那么系综分子动力学,4正那么系综分子动力学,5等温等压系综分子动力学,计算机分子模拟方法,第三章、分子动力学方法,一、分子动力学方法根底,原理:计算一组分子的相空间轨道,其中每个分 子各自服从牛顿运动定律,计算机分子模拟方法,第三章、分子动力学方法,计算元胞:立方体元胞,计算机分子模拟方法,第三章、分子动力学方法,边界条件:周期性边界条件,计算机分子模拟方法,第三章、分子动力学方法,设计算元胞的限度大小为L，其体积为L3,由于引用了这样的立方体箱子， 将产生六个我们不希望出现的外表， 这些外表的存在对系统的任何一种性质都会有重大的影响模拟中碰撞这些箱子的外表的粒子被反射回到元胞内部。为了减少引入的外表效应，采用周期性边界条件，构造出一个准无穷大的体积来更精确地代表宏观系统， 即让这个小体积元胞镶嵌在一个无穷大的大块物质之中。周期性边界条件的数学表示形式为：,，,A为任意可观测量， 为任意整数。 即令根本元胞完全等同的重复无穷屡次， 当有一个粒子穿过根本MD元胞的六方体外表时， 就让这个粒子以相同的速度穿过此外表对面的外表重新进入该MD元胞内。,最小影象:,根本元胞中的一个粒子只与根本元胞中的其它N-1个粒子,或它们的最近邻影象发生作用。,条件：位势的截断距离必须小于元胞线度L的一半。,计算机分子模拟方法,第三章、分子动力学方法,L,差分格式采用有限差分法将微分方程变成有限差分方程以便数值求解 ,哈密顿表述：,牛顿表述：,(1)Verlet算法;,(2)预测-校正格式.,计算机分子模拟方法,第三章、分子动力学方法,Verlet算法:,计算机分子模拟方法,第三章、分子动力学方法,算法启动,1扰动初始位置；,2利用初始位置和速度：,原始形式的算法表述：,1规定初始位置r0,r1,2计算第n步的力Fn,3计算第n+1步的位置：,4计算第n步的速度：,5重复(2)到(4),计算机分子模拟方法,第三章、分子动力学方法,计算机分子模拟方法,第三章、分子动力学方法,计算机分子模拟方法,第三章、分子动力学方法,预测校正格式,预测:,利用预测计算得:,修正a:,计算机分子模拟方法,第三章、分子动力学方法,校正:,计算机分子模拟方法,第三章、分子动力学方法,二、L-J流体体系,计算机分子模拟方法,第三章、分子动力学方法,二无量纲化:,计算机分子模拟方法,第三章、分子动力学方法,三热力学性质:,1温度,2内能,3压力：根据哈密顿动力学：,计算机分子模拟方法,第三章、分子动力学方法,四位势截断与长程修正:,1对分布函数：,计算机分子模拟方法,第三章、分子动力学方法,(2)根据分布函数理论,内能,压力,计算机分子模拟方法,第三章、分子动力学方法,五 F3(N,E,V)程序详解,1初始化,2量纲转换,3算法,4输入输出,练习：计算文献,COMPUTER “EXPERIMENTSOF CLASSICALFLUIDS(I)：,THERMODYNAMICAL PROERTIES OF LENNEARD-JONES,MOLECULES,Verlet L.,Phys.Rev,v.159,1967,P.98-103.,中的任一状态的比照内能、压力，并画出初始和结束位型。,程序F3(N,E,V),INDEX OF CODES CCP5,F01 PERIODIC BOUNDERY CONDITIONS IN VARIOUS GEOMENTRIES,F02 5-VALUE GEAR PERDICTOR-CORRECTOR ALGORITHM,F03 LOW-STORAGE PROGRAM USING LEAPFROG VERLET ALGORITHM,F04 VELOCITY VERSION OF THE VERLET ALGORITHM,F05 QUATERNION PARAMETER PERDICTOR-CORRECTOR ALGORITHM,F06 LEAPFROG ALGORITHMS FOR ROTATIONAL MOTION,F07 CONGSTRAINT DYNAMICS A CHAIN MOLECULE,F08 SHAKE ALGORITHM FOR CONGSTRAINT DYNAMICS A CHAIN MOLECULE,F09 RATTLE ALGORITHM FOR CONGSTRAINT DYNAMICS A CHAIN MOLECULE,F10 HARD SPHERE MOLECULAR DYNAMICS PROGRAM,F11 CONSTANT-NVT MONTE CARLO FOR LENNARD JONES ATOMS,F12 CONSTANT-NVT MONTE CARLO ALGORITHM,F13 THE HEART OF A CONSTANT-MUVT MONTE CARLO PROGRAM,INDEX OF CODES,F14 ALGORITHM TO HANDLE INDICES IN CONSTAN-MUVT MONTE CARLO,F15 ROUTINE TO RANDOMLY ROTATE MOLECULES,F16 HARD DUMB-BELL MONTE CARLO PROGRAM,F17 A SIMPLE LENNARD-JONES FORCE ROUTINE,F18 ALGORITHM FOR AVOIDING THE SQUARE ROOT OPERATION,F19 THE VERLET NEIGHBOUR LIST,F20 ROUTINES TO CONSTRUCT AND USE CELL LINKED-LIST METHOD,F21 MULTIPLE TIMESTEP MOLECULAR DYNAMICS,F22 ROUTINES TO PERFORM THE EWALD SUM,F23 ROUTINES TO SET ALPHA-FCC LATTICE OF LINEAR MOLECULES,F24 INITIAL VELOCITY DISTRIBUTION,F25 ROUTINE TO CALCULATE TRANSLATIONAL ORDER PARAMETER,F26 ROUTINES TO FOLD AND UNFOLD TRAJECTORIES IN PERIODIC BOUNDARIES,INDEX OF CODES,F27 PROGRAM TO COMPUTE TIME CORRELATION FUNCTIONS,F28 CONSTRANT-NVT MOLECULAR DYNAMICS PROGRAM(EXTENDED SYSTEM METHOD),F29 CONSTRANT-NVT MOLECULAR DYNAMICS PROGRAM (CONSTRAINT METHOD),F30 CONSTRANT-NPT MOLECULAR DYNAMICS PROGRAM(EXTENDED SYSTEM METHOD),F31 CONSTRANT-NPT MOLECULAR DYNAMICS PROGRAM (CONSTRAINT METHOD),F32 CELL LINKED-LISTS METHOD IN SHEARED BOUNDARIES,F33 BROWNIAN DYNAMICS FOR A LENNARD-JONES FLUID,F34 AN EFFICIENT CLUSTERING ROUTINE,F35 THE VORONOI POLYGON ALGORITHM IN 2D AND 3D,F36 MONTE CARLO SIMULATION OF HARD LINES IN 2D,F37 ROUTINES TO CALCULATE FOURIER TRANSFORM,程序F3(N,E,V),5物理量的计算：,温度：,内能：,压力：,计算机分子模拟方法,第三章、分子动力学方法,程序详解F3(N,E,V):,(1)初始化,WRITE(*,(“ENTER RUN TITLE ),READ(*,(A) TITLE,WRITE(*,(“ENTER NUMBER OF STEPS ),READ(*,*) NSTEP,WRITE(*,(“ENTER NUMBER OF STEPS BETWEEN OUTPUT LINES ),READ(*,*) IPRINT,WRITE(*,(“ENTER CONFIGURATION FILENAME ),READ(*,(A) CNFILE,WRITE(*,(“ENTER THE FOLLOWING IN LENNARD-JONES UNITS ),WRITE(*,(“ENTER DENSITY ),READ(*, *)DENS,WRITE(*,(“ENTER POTENTIAL CUTOFF DISTANCE ),READ(*,*) RCUT,WRITE(*,(“ENTER TIME STEP ),READ(*,*)DT,CALL READCN(CNFILE),计算机分子模拟方法,第三章、分子动力学方法,周期边界,最小影象,差分格式,位势截断,计算机分子模拟方法,第三章、分子动力学方法,MD根底,最小影象,根本元胞中的一个粒子只与根本元胞中的其它N-1个粒子，或它们的近邻影象发生作用。,条件：位势的阶段距离必须小于元胞线度的一半。,计算机分子模拟方法,第三章、分子动力学方法,氩,计算机分子模拟方法,第三章、分子动力学方法,蛙跳格式,计算机分子模拟方法,第三章、分子动力学方法,蛙跳格式:,计算机分子模拟方法,第三章、分子动力学方法,蛙跳格式:,计算机分子模拟方法,第三章、分子动力学方法,蛙跳格式:,热科学中的计算机分子模拟方法,第三章、分子动力学方法,初始位型：面心立方FCC,NC=(FLOAT(N)/4.)*(1./3.),XL=1/(FLOAT(NC),Y=.5*XL,R(1)=0,0,0,R(2)=0,Y,Y,R(3)=Y,0,Y,R(4)=Y,Y,0,M=0,DO 10 I=1,NC,DO 10 J=1,NC,DO 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+4,10 CONTINUE,计算机分子模拟方法,第三章、分子动力学方法,XL,NC=4,Y,Z,X,程序F3(N,E,V):,(2)量纲转换,SIGMA=(DENS/REAL(N)*(1.0/3.0),RCUT=RCUT*SIGMA,DT=DT*SIGMA,DENS=DENS/(SIGMA*3),计算机分子模拟方法,第三章、分子动力学方法,程序F3(N,E,V):,(3)长程校正：,SR3=(SIGMA/RCUT)*3,SR9=SR3*3,SIGCUB=SIGMA*3,VLRC=(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),计算机分子模拟方法,第三章、分子动力学方法,程序F3(N,E,V):,(4)算法,算法启动：,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(OLDK),CALL MOVE(DT),CALL FORCE(DT,SIGMA,RCUT,NEWV,NEWVC,NEWW),CALL KINET(NEWK),计算机分子模拟方法,第三章、分子动力学方法,程序F3(N,E,V):,(4)算法,差分格式：,SR2=SIGSQ/RIJSQ,VIJ=4.0*(SR12-SR6),WIJ=48.0*(SR12-0.5*SR6),VELIJ=WIJ*DT/RIJSQ,DVX=VELIJ*RXIJ, ,DXI=VXI+DVX, ,VX(J)=VX(J)-DVX, ,V=V+VIJ,W=W+WIJ,CALL MOVE(DT),计算机分子模拟方法,第三章、分子动力学方法,程序F3(N,E,V):,(4)算法,最小映像：,DO 1000 I=1,N-1,DO 900 J=I+1,N,RXIJ=RXI-RX(J),RXIJ=RXIJ-ANINT(RXIJ),IF(ABS(RXIJ).LT.RCUT)THEN,RYIJ=RYI-RY(J),RYIJ=RYIJ-ANINT(RYIJ),RIJSQ=RXIJ*2+RYIJ*2,IF(RIJSQ.LT.RCUTSQ)THEN,RZIJ=RZI-RZ(J),RZIJ=RZIJ-ANINT(RZIJ),RIJSQ=RIJSQ+RZIJ*2,计算机分子模拟方法,第三章、分子动力学方法,0.6,0.4,RXIJ=RXIJ-ANINT(RXIJ),Example: -0.4=0.6-1,0.4=0.4-0,程序F3(N,E,V):,(5)物理量的计算：,动能和势能：,K=(OLDK+NEWK)/2+(OLDV-2*V+NEWV)/8,WHERE K=K(T),V=V(T),OLDK=K(T-DT/2),OLDV=V(T-DT),NEWK=K(T+DT/2),NEWV=V(T+DT),VIJ=4.0*(SR12-SR6),WIJ=48.0*(SR12-0.5*SR6),V=V+VIJ,W=W+WIJ,KN=K/REAL(N),TEMP=2.0*KN/FREE,PRES=DENS*TEMP+W,计算机分子模拟方法,第三章、分子动力学方法,具有确定的粒子数N、能量E和体积V,算法：,1规定初始位置和速度；,2对运动方程积分假设干步；,3计算势能和动能；,4假设能量不等于所需要的值，对速度进行标度；,5重复2到4，直到系综平衡。,计算机分子模拟方法,第三章、分子动力学方法,三、微正那么系综分子动力学,速度的标度:,计算机分子模拟方法,第三章、分子动力学方法,蛙跳格式的标度:,1规定初始位置、速度,2计算第n+1/2步的速度：,3计算第n+1步的位置：,4计算第n步的速度：,5计算标度因子，对速度进行标度：,重复2到5,计算机分子模拟方法,第三章、分子动力学方法,速度的标度:,SUMKX=0.0,SUMKY=0.0,SUMKZ=0.0,DO 200 I=1,N,SUMKX=SUMKX+VX(I)*VX(I)/REAL(N),SUMKY=SUMKY+VY(I)*VY(I)/REAL(N),SUMKZ=SUMKZ+VZ(I)*VZ(I)/REAL(N),200 CONTINUE,BEITA=SQRT(3.*TEMP/(SUMKX+SUMKY+SUMKZ),DO 300 I=1,N,VX(I)=VX(I)*BEITA,VY(I)=VY(I)*BEITA,VZ(I)=VZ(I)*BEITA,300 CONTINUE,计算机分子模拟方法,第三章、分子动力学方法,初始速度:,简单的选择：,计算机分子模拟方法,第三章、分子动力学方法,初始速度:,简单的选择：,计算机分子模拟方法,第三章、分子动力学方法,RTEMP=SQRT(TEMP),DO 100 I=1,N,VX(I)=RTEMP*(RANF(DUMMY)-0.5 ),100 CONTINUE,SUMKX=0.0, ,DO 200 I=1,N,SUMKX=SUMKX+VX(I)*VX(I)/REAL(N), ,200 CONTINUE,BEITAX=SQRT(TEMP/SUMKX), ,DO 300 I=1,N,VX(I)=VX(I)*BEITAX, ,300 CONTINUE,初始速度:,简单的选择：,控制总量为零：,计算机分子模拟方法,第三章、分子动力学方法,SUMX=0.0,SUMY=0.0,SUMZ=0.0,DO 200 I=1,N,SUMX=SUMX+VX(I),SUMY=SUMY+VY(I),SUMZ=SUMZ+VZ(I),200 CONTINUE,SUMX=SUMX/REAL(N),SUMY=SUMY/REAL(N),SUMZ=SUMZ/REAL(N),DO 300 I=1,N,VX(I)=VX(I)-SUMX,VY(I)=VY(I)-SUMX,VZ(I)=VZ(I)-SUMX,300 CONTINUE,初始速度:,从MAXWELL分布中抽样,高斯分布:,对于等几率随机实验(bernoulli实验)，大量的实验结果满足高斯分布：,计算机分子模拟方法,第三章、分子动力学方法,初始速度:,从MAXWELL分布中抽样,高斯分布:,MAXWELL速度分布：,由于：,单位体积的分子在每个分量上的速度分布实际上就是高斯分布,计算机分子模拟方法,第三章、分子动力学方法,初始速度:,从MAXWELL分布中抽样,高斯分布的随机数生成:,分布：,生成随机数：,计算机分子模拟方法,第三章、分子动力学方法,a1,a3,a5,a7,a9是定数,初始速度:,从MAXWELL分布中抽样,高斯分布的随机数生成:,计算机分子模拟方法,第三章、分子动力学方法,SUM=0.0,DO 10 I=1,12,SUM=SUM+RANF(DUMMY),10 CONTINUE,R=(SUM-6.0)/4.0,R2=R*R,GAUSS=(A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R,初始速度:从MAXWELL分布中抽样,从MAXWELL分布中抽样:,控制总动量为零：同前,计算机分子模拟方法,第三章、分子动力学方法,RTEMP=SQRT(TEMP),DO 100 I=1,N,VX(I)=RTEMP*GAUSS(DUMMY),VY(I)=RTEMP*GAUSS(DUMMY),VZ(I)=RTEMP*GAUSS(DUMMY),100 CONTINUE,L-J流体NVE系综分子动力学模拟程序及运行结果,real x(1:786),vh(1:786),F(1:786),real FY(1:256),FZ(1:256),Y(1:256),Z(1:256),real H,HSQ,HSQ2,TREF,real KX,KY,KZ,INTEGER CLOCK ,TIMEMX,EQUIVALENCE (FY,F(257),(FZ,F(513),(Y,X(257),(Z,X(513),open(4,file=zh4.dat),NPART=256,DEN=0.83134,SIDE=6.75284,TREF=0.722,RCOFF=2.5,H=0.064,IREP=10 !每隔10步调整一次速度,ISTOP=500 ! 50步后停止速度标度,TIMEMX=1350 !300步后退出,打印计算结果,iseed=615843,WRITE(4,*)MOLECULAR DYNAMICS SIMULATION PROGRAM,WRITE(4,*)NUMBER OF PARTICLESIS,NPART,WRITE(4,*)SIDE LENGTH OF THE BOX IS,SIDE,WRITE(4,*)CUT OFF IS,RCOFF,WRITE(4,*)TREF,TREF,WRITE(4,*)H=,H,A =SIDE/4.0,SIDEH=SIDE*0.5,HSQ=H*H,HSQ2=HSQ*0.5,NPARTM=NPART-1,RCOFFS=RCOFF*RCOFF,TSCALE=16.0/(1.0*NPART-1.0),VAVER=1.13*SQRT(TREF/24.0),N3=3*NPART,IOF1=NPART,IOF2=2*NPART,IJK=0,DO 10 LG=0,1,DO 10 I=0,3,DO 10 J=0,3,DO 10 K=0,3,IJK=IJK+1,X(IJK)=I*A+LG*A*0.5,Y(IJK)=J*A+LG*A*0.5,Z(IJK)=K*A,10 CONTINUE,DO 15 LG=1,2,DO 15 I=0,3,DO 15 J=0,3,DO 15 K=0,3,IJK=IJK+1,X(IJK)=I*A+(2-LG)*A*0.5,Y(IJK)=J*A+(LG-1)*A*0.5,Z(IJK)=K*A+A*0.5,15 CONTINUE,CALL MXWELL(VH,N3,H,TREF,iseed),DO 50 I=1,N3,F(I)=0.0,50 CONTINUE,C-,DO 200 CLOCK=1,TIMEMX,DO 210 I=1,N3,X(I)=X(I)+VH(I)+F(I),210 CONTINUE,DO 215 I=1,N3,IF(X(I).LT.0) X(I)=X(I)+SIDE,IF(X(I).GT.SIDE) X(I)=X(I)-SIDE,215 CONTINUE,DO 220 I=1,N3,VH(I)=VH(I)+F(I),220 CONTINUE,VIR=0.0,EPOT=0.0,DO 225 I=1,N3,F(I)=0.0,225CONTINUE,DO 270 I=1,NPART,XI=X(I),YI=Y(I),ZI=Z(I),DO 270 J=I+1,NPART,XX=XI-X(J),YY=YI-Y(J),ZZ=ZI-Z(J),IF(XX.LT.-SIDEH)XX=XX+SIDE,IF(XX.GT.SIDEH)XX=XX-SIDE,IF(YY.LT.-SIDEH)YY=YY+SIDE,IF(YY.GT.SIDEH)YY=YY-SIDE,IF(ZZ.LT.-SIDEH)ZZ=ZZ+SIDE,IF(ZZ.GT.SIDEH)ZZ=ZZ-SIDE,RD=XX*XX+YY*YY+ZZ*ZZ,IF(RD.GT.RCOFFS)GOTO 270,EPOT=EPOT+RD*(-6.0)-RD*(-3.0),R148=RD*(-7.0)-0.5*RD*(-4.0),VIR=VIR-RD*R148,KX=XX*R148,F(I)=F(I)+KX,F(J)=F(J)-KX,KY=YY*R148,FY(I)=FY(I)+KY,FY(J)=FY(J)-KY,KZ=ZZ*R148,FZ(I)=FZ(I)+KZ,FZ(J)=FZ(J)-KZ,270 CONTINUE,DO 275 I=1,N3,F(I)=F(I)*HSQ2,275 CONTINUE,DO 300 I=1,N3,VH(I)=VH(I)+F(I),300 CONTINUE,EKIN=0.0,DO 305 I=1,N3,EKIN=EKIN+VH(I)*VH(I),305 CONTINUE,EKIN=EKIN/HSQ,VEL=0.0,COUNT=0.0,DO 306 I=1,NPART,VX=VH(I)*VH(I),VY=VH(I+IOF1)*VH(I+IOF1),VZ=VH(I+IOF2)*VH(I+IOF2),SQ=SQRT(VX+VY+VZ),SQT=SQ/H,IF(SQT.GT.VAVER)COUNT=COUNT+1,VEL=VEL+SQ,306 CONTINUE,VEL=VEL/H,IF(CLOCK.LT.ISTOP)THEN,IF(CLOCK/IREP)*IREP.EQ.CLOCK)THEN,WRITE(4,*)VELOCITY ADJUSTMENT,TS=TSCALE*EKIN,WRITE(4,*)TEMPERATURE BEFORE SCALING IS,TS,SC=TREF/TS,SC=SQRT(SC),WRITE(4,*)SCALE FACTOR IS,SC,DO 310 I=1,N3,VH(I)=VH(I)*SC,310 CONTINUE,EKIN=TREF/TSCALE,END IF,END IF,IF(CLOCK/10)*10.EQ.CLOCK)THEN,EK=24.0*EKIN,EPOT=4.0*EPOT,ETOT=EK+EPOT,TEMP=TSCALE*EKIN,PRES=DEN*16.0*(EKIN-VIR)/NPART,VEL=VEL/NPART,RP=(COUNT/256.0)*100.0,WRITE(4,6000)CLOCK,EK,EPOT,ETOT,TEMP,PRES,VEL,RP,WRITE(6,6000)CLOCK,EK,EPOT,ETOT,TEMP,PRES,VEL,RP,END IF,200 CONTINUE,6000 FORMAT(1I6,7F15.6),close(4),END,C-,SUBROUTINE MXWELL(VH,N3,H,TREF,iseed),REAL VH(1:N3),REAL U1,U2,V1,V2,S,R,NPART=N3/3,IOF1=NPART,IOF2=2*NPART,TSCALE=16.0/(1.0*NPART-1.0),DO 10 I=1,N3,2,1 CALL GGUB(ISEED ,RAN),U1=RAN,CALL GGUB(ISEED,RAN),U2=RAN,V1=2.0*U1-1.0,V2=2.0*U2-1.0,S=V1*V1+V2*V2,IF(S.GE.1.0)GOTO 1,R=-2.0*ALOG(S)/S,VH(I)=V1*SQRT(R),VH(I+1)=V2*SQRT(R),10 CONTINUE,EKIN=0.0,SP=0.0,DO 20 I=1,NPART,SP=SP+VH(I),20 CONTINUE,SP=SP/NPART,DO 21 I=1,NPART,VH(I)=VH(I)-SP,EKIN=EKIN+VH(I)*VH(I),21 CONTINUE,WRITE(4,*)TOLAL LINEAR MOMENTUM IN X DIRECTION,SP,SP=0.0,DO 22 I=IOF1+1,IOF2,SP=SP+VH(I),22 CONTINUE,SP=SP/NPART,DO 23 I=IOF1+1,IOF2,VH(I)=VH(I)-SP,EKIN=EKIN+VH(I)*VH(I),23 CONTINUE,WRITE(4,*)TOTAL LINEAR MOMENTUM IN Y DIRECTION IS,SP,SP=0.0,DO 24 I=IOF2+1,N3,SP=SP+VH(I),24 CONTINUE,SP=SP/NPART,DO 25 I=IOF2+1,N3,VH(I)=VH(I)-SP,EKIN=EKIN+VH(I)*VH(I),25 CONTINUE,WRITE(4,*)TOTAL LINEAR MOMENTUM IN Z DIRECTION IS,SP,WRITE(4,*)VELOCITY AD JUSTIMENT,TS=TSCALE*EKIN,WRITE(4,*)TEMPERATURE BEFORE SCALING IS,TS,SC=TREF/TS,SC=SQRT(SC),WRITE(4,*)SCALE FACTOR IS,SC,SC=SC*H,DO 30 I=1,N3,VH(I)=VH(I)*SC,30 CONTINUE,END,subroutine GGUB(iseed,ran),double precision z,maxf,minf,data maxf/2147483647.d0/,data minf/0.465661287307739d-9/,z=dble(iseed),z=dmod(16807.d0*z,maxf),ran=real(z*minf),iseed=int(z),return,end,模拟结果,具有确定的粒子数N、温度T和体积V,问题的关键：温度的约束,1速度的直接标度,2热浴方法(Andersen Thermostat),3约束方法阻尼力方法(Constraint Method),4系统扩展方法(Extended System Method),计算机分子模拟方法,第三章、分子动力学方法,四、正那么系综分子动力学,1速度的直接标度,热浴方法(Andersen Thermostat),引入一个与虚拟粒子碰撞的随机力,想象系综浸没在热浴当中,1系综和热浴间的相互作用强度由随机碰撞的频率决定,2碰撞的几率等于Nu*dt,3如果一个粒子经历碰撞，它的速度将从约束温度下的Maxwell分布中随机抽取,计算机分子模拟方法,第三章、分子动力学方法,约束方法阻尼力方法：,是等动能分子动力学方法,系统的运动方程为：,引入阻尼系数 以保证将温度约束在恒定值：,根据高斯最小约束原理：,计算机分子模拟方法,第三章、分子动力学方法,(3)约束方法,系综扩展方法：,1Nose Thermostat:,设想原系统与一个耦合系统共同组成一个扩展系统；,扩展系统的拉格朗日函数为：,S为扩展坐标变量，Q为等效质量，g为扩展系统的自由度，扩展系统的哈密顿函数为：,计算机分子模拟方法,第三章、分子动力学方法,系综扩展方法：,1Nose Thermostat:,其中：,扩展系综的的运动方程为：,计算机分子模拟方法,第三章、分子动力学方法,系综扩展方法：,2Nose-Hoover Thermostat:,Hoover 通过将Nose方法中的s、 、Q整合和简化，引入了热力学阻尼系数,系综的运动方程简化为：,计算机分子模拟方法,第三章、分子动力学方法,(4),系综扩展方法,具有确定的粒子数N、压力P和温度T,问题的关键：温度和压力的约束,1温度约束同前,2压力的约束算法,3压力的系综扩展方法,计算机分子模拟方法,第三章、分子动力学方法,五、等温等压系综分子动力学,压力的扩展方法Andersen方法：,与温度的扩展方法类似(Nose方法；,扩展系统的拉格朗日函数为：,X为标度后的坐标变量,Q为等效体积,m和 为常数：,扩展系统的哈密顿函数为：,计算机分子模拟方法,第三章、分子动力学方法,压力的扩展方法(Andersen方法)：,其中：,扩展系综的的运动方程为：,计算机分子模拟方法,第三章、分子动力学方法,等温等压的约束算法:,与等温约束方法相似;,通过引入拉格朗日乘子和高斯最小约束,系统的运动方程为:,其中：,计算机分子模拟方法,第三章、分子动力学方法,定压定粒子数系综,定压定温系综,定压定温系综NPT系综,把恒压算法和速度标度结合起来。,算法 ： 同定压,对所用速度进行标度。,返回,参考文献,温度的热浴算法：,W.G.Hoover etal. Phys.Rev.leff. 1982,48:1818,D.J.Evans,J.Chem.Phys. 1983,78:3297,W.G.Hoover, A.Rev. Phys.chem. 1983， 34:103(高斯最小约束),温度的扩展算法：,S.Nose,J.Chen.Phys. 1984,51:511,WG.Hoover, Phys.Rev.A,1985,31:1695,压力的扩展算法：,H.C.Andersen,J.Chem.Phys.1980,72:2384,J.M.Haile and H.W.Graben,J.Chem.Phys,1980,73:2412,压力的约束算法：,D.J.Evans, and G.P.Morriss.Chem.Phys.1983,77:63-66,D.J.Evans, and G.P.Morriss.Chem.Phys.1984,1:297-344,