130 likes | 149 Views
Chapter 9: Molecular-dynamics. Integrate equations of motion-- classical! Discrete form of Newton’s second law Forces from interaction potential For a simple pair potential, we get. Integrating equations of motion. This works out to give the Verlet algorithm,.
E N D
Chapter 9: Molecular-dynamics • Integrate equations of motion-- classical! • Discrete form of Newton’s second law • Forces from interaction potential • For a simple pair potential, we get
Integrating equations of motion This works out to give the Verlet algorithm,
Placing atoms on fcc lattice… non-primitive basis ! start with non-primitive fcc basis r1(1)=0.0d0 r2(1)=0.0d0 r3(1)=0.0d0 r1(2)=0.5d0 r2(2)=0.5d0 r3(2)=0.0d0 r1(3)=0.5d0 r2(3)=0.0d0 r3(3)=0.5d0 r1(4)=0.0d0 r2(4)=0.5d0 r3(4)=0.5d0
Some parameters for lattice… INTEGER, PARAMETER :: nn=10 INTEGER, PARAMETER :: Prec14=SELECTED_REAL_KIND(14) INTEGER, PARAMETER :: natoms=4*nn**3 INTEGER, PARAMETER :: nx=nn,ny=nn,nz=nn REAL(KIND=Prec14), PARAMETER :: sigma=1.0d0,epsilon=1.0d0 REAL(KIND=Prec14), PARAMETER :: L=sigma*2.0d0**(2.0d0/3.0d0)*nn
Generate the full lattice n=0 do ix=1,nx do iy=1,ny do iz=1,nz do i=1,4 n=n+1 rx(n)=(r1(i)+dble(ix-1))*sigma ry(n)=(r2(i)+dble(iy-1))*sigma rz(n)=(r3(i)+dble(iz-1))*sigma enddo enddo enddo enddo ! scale lengths by box size rx=rx/nx ry=ry/ny rz=rz/nz
Some preliminaries… ! Compute potential and force at cutoff sovr=sigma/rcut vcut=4.0d0*epsilon*(sovr**12-sovr**6) fcut=(24.0d0*epsilon/rcut)*(2.0d0*sovr**12-sovr**6)
Initialize velocities using random-number generator ! Initialize velocities by picking last position call RANDOM_SEED do n=1,natoms call RANDOM_NUMBER(ranx) call RANDOM_NUMBER(rany) call RANDOM_NUMBER(ranz) rxl(n)=rx(n)+v0*(0.5d0-ranx) ryl(n)=ry(n)+v0*(0.5d0-rany) rzl(n)=rz(n)+v0*(0.5d0-ranz) enddo Must adjust everything to give zero net momentum…
Implementing periodic boundary conditions: Computing forces and potential energy… do i=1,natoms-1 do j=i+1,natoms sxij=rx(i)-rx(j) syij=ry(i)-ry(j) szij=rz(i)-rz(j) sxij=sxij-dble(anint(sxij)) syij=syij-dble(anint(syij)) szij=szij-dble(anint(szij)) rij=L*dsqrt(sxij**2+syij**2+szij**2) if(rij.lt.rcut) then sovr=sigma/rij epot=epot+4.0d0*epsilon*(sovr**12-sovr**6)-vcut+fcut*(rij-rcut) fij=-(24.0d0*epsilon/rij)*(2.0d0*sovr**12-sovr**6)+fcut get x,y,z components of force, use Newton’s 3rd law…
Advance atoms, compute kinetic energy…. ! Verlet algorithm to move atoms, compute kinetic energy do i=1,natoms rxx=2.0d0*rx(i)-rxl(i)+fx(i)*dt**2/(amass*L) ryy=2.0d0*ry(i)-ryl(i)+fy(i)*dt**2/(amass*L) rzz=2.0d0*rz(i)-rzl(i)+fz(i)*dt**2/(amass*L) vx=(rxx-rxl(i))*L/(2.0d0*dt) vy=(ryy-ryl(i))*L/(2.0d0*dt) vz=(rzz-rzl(i))*L/(2.0d0*dt) vsq=vx**2+vy**2+vz**2 ekin=ekin+0.5d0*amass*vsq rxl(i)=rx(i) ryl(i)=ry(i) rzl(i)=rz(i) rx(i)=rxx ry(i)=ryy rz(i)=rzz enddo
Output including time and temperature ekin=ekin/natoms; epot= epot/natoms; etot=epot+ekin time=time+tunit*dt temp=(2.0d0/3.0d0)*eps*ekin write(6,100) n,time,ekin,epot,etot,temp enddo 100 format(i8,f12.2,4f18.6) eps=120 for Argon potential…. energy scale Ekin,epot, temp will fluctuate… etot should be fairly Stable if small enough time step used…
Heat capacity • Angle brackets thermal average • In MD, thermal averages done by time average • Heat capacity given by fluctuations in total energy The computed epot and temperature can be used to determine heat capacity (output in units of kB)
Radial distribution function • is the density N/ • Dirac-delta defined numerically (not infinitely sharp) • In an ideal gas, g(r) =1 • In a liquid, g(r) =1 at long ranges, short range structure • In a crystal, g( r) has sharp peaks, long-range order