# How can I make velocity verlet algorithm more stable?

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 :: 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

implicit none
double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
integer :: n
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
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 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

!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
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)/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 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
implicit none

integer n,i
double precision,dimension(n) :: posx,posy,posz,velx,vely,velz,ax,ay,az
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

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 + vely(i)**2 + velz(i)**2)/(2*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 + (y1-y2)**2 + (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) - (6.d0/r**8))*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) - (6.d0/r**8))*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


You can try using an adaptive timestep, where dt is adjusted at each step if the velocities get too large. However, the Verlet method is symplectic, so if the total energy is not roughly conserved then something is wrong. If two particles get too close, no reasonable timestep will make the method stable, and accurately computing the forces between particles with double-precision floats will be impossible.
You didn't ask about coding style, but all the same, I recommend organizing your code into modules and declaring the intent of every procedure argument. Fortran Best Practices has more good information on this subject.