The answer to this question implies that reducing the time step would make it more stable. However I have tried reducing the time step but the system is still unstable(the total energy increases to a very large value). What I understand is that if two balls are very close and the end of some step, in the next step they will have very large velocity, thus making the system unstable. How can I choose the time step so that the system remains stable? I would also appreciate if someone could look into the code I provide to figure out if there is some error in implementation of the algorithm?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  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=200          !Number of particles  eps=1.d0  m=1.d0  xlen=20.d0  ylen=20.d0  zlen=20.d0  xinit=0.d0  yinit=0.d0  zinit=0.d0  mindist=1.5d0  rad=mindist/2.d0  tot_t=10.d0  dt=0.001d0  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     call random_number(x)     velx(i)=02.1d0*(2.d0*x-1.d0)     call random_number(y)     vely(i)=02.1d0*(2.d0*y-1.d0)     call random_number(z)     velz(i)=02.1d0*(2.d0*z-1.d0)   enddo  print*, Velocity initialisation finished  !Assigning 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)/menddo  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 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 + vely(i)**2 + velz(i)**2)/(2*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 + (y1-y2)**2 + (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) - (6.d0/r**8))*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) - (6.d0/r**8))*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 can I make velocity verlet algorithm more stable?
computational physics;molecular dynamics