1 / 13

Molecular Dynamics Simulation: Integrating Equations with the Verlet Algorithm for Lennard-Jones Potential

This chapter delves into molecular dynamics simulation, covering the integration of equations of motion using the classical approach, focusing on the Verlet algorithm for noble gas atoms. Starting with a non-primitive FCC lattice basis, the text explains the placement of atoms and parameters for lattice generation. It details force calculations, implementing periodic boundary conditions, computing forces, and potential energy. The Verlet algorithm is utilized to move atoms and compute kinetic energy. The output includes time, temperature, and energy components. Additionally, heat capacity and the radial distribution function are discussed, highlighting the thermal average and the density-based function g(r) in various states of matter.

hazelgreen
Download Presentation

Molecular Dynamics Simulation: Integrating Equations with the Verlet Algorithm for Lennard-Jones Potential

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 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

  2. Integrating equations of motion This works out to give the Verlet algorithm,

  3. Lennard-Jones potential for noble gas

  4. 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

  5. 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

  6. 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

  7. 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)

  8. 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…

  9. 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…

  10. 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

  11. 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…

  12. 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)

  13. 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

More Related