How to keep velocities in check in molecular dynamics simulation?

Question

I am trying to make a very simple molecular dynamics simulator with Reflective Boundary Conditions. I am assigning the initial positions in a cube randomly while making sure they are not too close to each other and keeping the initial velocities zero. I am using the Lennard-Jones potential and Velocity-Verlet algorithm. But as soon as I start the simulation, Velocities of some of the particles reach very high value and total energy also goes to very high value(not consistent with the initial total energy). Can anyone point out what are the possible reasons for such a behaviour?

Here is my code (in Fortran).

program mol_dyn_ref
  implicit none
  double precision,allocatable,dimension(:) :: posx,posy,posz,velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  integer,allocatable,dimension(:) :: seed
  double precision :: x,y,z,m,rad,eps,kin,pot,k,p,avgvx,avgvy,avgvz
  double precision :: xinit,yinit,zinit,xlen,ylen,zlen,mindist
  integer :: i,j,st,seedsize,n,iter,t
  double precision :: tot_t,dt,dist,r,Fx,Fy,Fz

  interface
     subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
       implicit none
       integer :: i,j,n
       double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
       double precision :: r,dist,Fx,Fy,Fz,eps,m
     end subroutine acc

     subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
       implicit none
       double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
       integer :: n
       double precision :: tot_t,dt,xlen,ylen,zlen,rad
     end subroutine pos_upd

     subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
       implicit none
       integer n
       double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
       double precision :: dt
     end subroutine vel_upd
  end interface



  n=50         !Number of particles
  eps=1.0d0
  m=1.0d0
  xlen=20.0d0
  ylen=20.0d0
  zlen=20.0d0
  xinit=0.0d0
  yinit=0.0d0
  zinit=0.0d0
  mindist=1.5d0
  rad=mindist/2.0d0
  tot_t=10.0d0
  dt=0.0001d0

  iter=int(tot_t/dt)

  allocate(posx(n),posy(n),posz(n),velx(n),vely(n),velz(n),ax(n),ay(n),az(n),tempx(n),tempy(n),tempz(n))

  open(100,file="pos.dat",status="replace")
  open(200,file="vel.dat",status="replace")
  open(300,file="acc.dat",status="replace")
  open(400,file="energy.dat",status="replace")

  call random_seed(size=seedsize)
  allocate(seed(seedsize))
  do i=1,seedsize
     call system_clock(st)
     seed(i)=st
  enddo
  call random_seed(put=seed)

  !Assigning initial position to first particle
10 call random_number(x)
  posx(1)=xinit+x*xlen
  call random_number(y)
  posy(1)=yinit+y*ylen
  call random_number(z)
  posz(1)=zinit+z*zlen

  if(posx(1)<rad .OR. posx(1)>xlen-rad) goto 10
  if(posy(1)<rad .OR. posy(1)>ylen-rad) goto 10
  if(posz(1)<rad .OR. posz(1)>zlen-rad) goto 10

  !Assigning initial position
  do i=2,n
20   call random_number(x)
     posx(i)=xinit+x*xlen
     call random_number(y)
     posy(i)=yinit+y*ylen
     call random_number(z)
     posz(i)=zinit+z*zlen
     if(posx(i)<rad .OR. posx(i)>xlen-rad) goto 20
     if(posy(i)<rad .OR. posy(i)>ylen-rad) goto 20
     if(posz(i)<rad .OR. posz(i)>zlen-rad) goto 20
     do j=1,i-1
        if (dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))<mindist) goto 20
     enddo
  enddo
  print*, "Position initialisation finished"

  !Assigning initial velocities
  do i=1,n
     velx(i)=0.d0
     vely(i)=0.d0
     velz(i)=0.d0
  enddo     
  print*, "Velocity initialisation finished"

  !Calculating initial acceleration
  do i=1,n
     ax(i)=0
     ay(i)=0
     az(i)=0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
     enddo
  enddo
  print*, "Acceleration initialisation finished."

  !Molecular Dynamics Simulation
  do t=1,iter
     print*, t
     do i=1,n
        write(100,*) posx(i),posy(i),posz(i)
        write(200,*) velx(i),vely(i),velz(i)
        write(300,*) ax(i),ay(i),az(i)
     enddo

     k=kin(velx,vely,velz,m,n)
     p=pot(posx,posy,posz,eps,n)
     write(400,*) t,k,p,k+p

     call pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
     call acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
     call vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)
  enddo
  call system("gnuplot --persist plot.gp")
endprogram mol_dyn_ref


!Updating Acceleration
subroutine acc(n,posx,posy,posz,m,eps,tempx,tempy,tempz)
  implicit none
  integer :: i,j,n
  double precision,dimension(n) :: posx,posy,posz,ax,ay,az,tempx,tempy,tempz
  double precision :: r,dist,Fx,Fy,Fz,eps,m
  do i=1,n
     tempx(i)=ax(i)
     tempy(i)=ay(i)
     tempz(i)=az(i)
     ax(i)=0.d0
     ay(i)=0.d0
     az(i)=0.d0
     do j=1,n
        if (j==i) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        ax(i)=ax(i)+ Fx(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        ay(i)=ay(i)+ Fy(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
        az(i)=az(i)+ Fz(r,posx(i)-posx(j),posy(i)-posy(j),posz(i)-posz(j),eps)/m
     enddo
  enddo
end subroutine acc

!Updating Position
subroutine pos_upd(n,posx,posy,posz,velx,vely,velz,ax,ay,az,dt,xlen,ylen,zlen,rad)
  implicit none

  integer n,i
  double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
  double precision :: dt,xlen,ylen,zlen,rad
  do i=1,n
     posx(i)=posx(i) + velx(i)*dt + (ax(i)*(dt**2))/2
     posy(i)=posy(i) + vely(i)*dt + (ay(i)*(dt**2))/2
     posz(i)=posz(i) + velz(i)*dt + (az(i)*(dt**2))/2

     if(posx(i)<rad .OR. posx(i)>xlen-rad) velx(i)=-velx(i)
     if(posy(i)<rad .OR. posy(i)>ylen-rad) vely(i)=-vely(i)
     if(posz(i)<rad .OR. posz(i)>zlen-rad) velz(i)=-velz(i)
  enddo
end subroutine pos_upd

!Updating Velocity
subroutine vel_upd(n,velx,vely,velz,ax,ay,az,dt,tempx,tempy,tempz)

  implicit none
  integer n,i
  double precision,dimension(n) :: velx,vely,velz,ax,ay,az,tempx,tempy,tempz
  double precision :: dt
  do i=1,n
     velx(i)=velx(i) + 0.5d0*(ax(i)+tempx(i))*dt
     vely(i)=vely(i) + 0.5d0*(ay(i)+tempy(i))*dt
     velz(i)=velz(i) + 0.5d0*(az(i)+tempz(i))*dt
  enddo
end subroutine vel_upd

function pot(posx,posy,posz,eps,n)
  implicit none
  double precision pot,r,dist,eps
  integer i,j,n
  double precision, dimension(n) :: posx,posy,posz
  pot=0.d0
  do i=1,n
     do j=1,n
        if(i==j) cycle
        r=dist(posx(i),posy(i),posz(i),posx(j),posy(j),posz(j))
        pot= pot + (4.d0*eps*((1.d0/r)**12 - (1.d0/r)**6))                 !r is relative distance. x,y,z are components of r.
     enddo
  enddo
end function pot

function kin(velx,vely,velz,m,n)
  implicit none
  double precision :: kin,m
  integer :: i,n
  double precision, dimension(n) :: velx,vely,velz
  kin=0.d0
  do i=1,n
     kin = kin + ((velx(i)**2.d0 + vely(i)**2.d0 + velz(i)**2.d0)/(2.d0*m))
  enddo
end function kin

function dist(x1,y1,z1,x2,y2,z2)
  implicit none
  double precision :: dist,x1,y1,z1,x2,y2,z2
  dist = sqrt((x1-x2)**2.d0 + (y1-y2)**2.d0 + (z1-z2)**2)
end function dist

function Fx(r,x,y,z,eps)
  implicit none
  double precision :: Fx,r,x,y,z,eps
  Fx = 4.d0*eps*((12.d0/r**14.d0) - (6.d0/r**8.d0))*x
end function Fx

function Fy(r,x,y,z,eps)
  implicit none
  double precision :: Fy,r,x,y,z,eps
  Fy = 4.d0*eps*((12.d0/r**14.d0) - (6.d0/r**8.d0))*y
end function Fy

function Fz(r,x,y,z,eps)
  implicit none
  double precision :: Fz,r,x,y,z,eps
  Fz = 4.d0*eps*((12.d0/r**14) - (6.d0/r**8))*z
end function Fz

No, the kinetic energy reaches very high value and so does the total energy. — Yogesh Yadav, May 24 '15 at 15:01
What happens if you set the LJ interactions to zero (that is, simulate an ideal gas)? Is the energy conserved in that case? If it is, it would suggest that there's an implementation error in your forces / potential. If the kinetic energy is not conserved, it could be that your implementation of Verlet is buggy. Also, what happens if you halve the time step length? Does energy conservation improve? — Juan M. Bello-Rivas, May 24 '15 at 15:11
Also, are the particles causing the spikes in kinetic energy close to the boundary? Elastic collisions with the boundary can lead to particles being too close to each other in a single time step. You would see spikes in potential energy too in that case. One way around this would be to endow the boundary with some sort of Weeks-Chandler-Andersen potential instead of perfectly elastic collisions. — Juan M. Bello-Rivas, May 24 '15 at 15:32
@Juan I tried to simulate the ideal gas by setting initial acceleration to zero, preventing any updation of acceleration and giving random initial velocities. The simulation runs fine in this case and the kinetic and total energy are conserved. — Yogesh Yadav, May 24 '15 at 16:31
Then my diagnostic is that it's either elastic collisions or some bug related to the forces and the corresponding updates. One way to rule out the latter cause of problems would be to re-enable your LJ potential and replace the reflecting boundaries by periodic boundaries, then see if that keeps your total energy conserved. — Juan M. Bello-Rivas, May 24 '15 at 16:45
@Juan One more thing I noticed that the total energy is conserved a lot better if I reduce the number of particles from 50 to very small (5). Does this give any additional hint? — Yogesh Yadav, May 24 '15 at 17:11

How to keep velocities in check in molecular dynamics simulation?

0 Answers0