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 particle10 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,n20   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 Accelerationsubroutine 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  enddoend subroutine acc!Updating Positionsubroutine 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)  enddoend subroutine pos_upd!Updating Velocitysubroutine 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  enddoend subroutine vel_updfunction 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  enddoend function potfunction 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))  enddoend function kinfunction 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 distfunction 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))*xend function Fxfunction 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))*yend function Fyfunction 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))*zend function Fz
How to keep velocities in check in molecular dynamics simulation?
computational physics;time integration;molecular dynamics;integration