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

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

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.

While it looks like you're randomly initializing the particle positions so that they are not too close to each other, you don't appear to have any guards on the particle velocities to keep the particles from approaching too close at the first update. The Lenard-Jones potential that you're using is very stiff at close distances, so if you're not careful in how you initialize the system you can have one particle fly off to infinity right away.

Instead, try initializing the system in a "cold" state: set all the velocities to zero and put the particles in a hexagonally close-packed arrangement, with the distance equal to the equilibrium distance for the potential you've chosen. You'll see some amount of oscillation about the true equilibrium configuration. This will be a much easier system to solve so that you can verify the correctness and stability of your algorithm. From there, you can slowly increase the velocities of the initial arrangement until the system reaches the liquid or gaseous state.

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.

• you don't appear to have any guards on the particle velocities to keep them from approaching too close at the next timestep Did you mean appraching too high? I thought potential should take care of velocity updation. Wouldn't changing velocities by hand change the simulation altogether? May 23 '15 at 18:16
• In the main loop, you update the positions before updating the velocities or accelerations. While you have taken precautions to keep the particles far enough apart at the start, your first update for the particle positions could throw that out the window if the initial velocities are too big. You'd have to work out what "too big" is given the parameters of your simulation. The simplest thing is to set the velocities to zero and make sure that works first. In principle, the potential should take care of the issue, provided the system is initialized in a state that isn't too unstable. May 23 '15 at 18:55
• I tried setting the initial velocities to zero. It doesn't improve the situation. May 23 '15 at 19:25
• Alright, you're still starting with 200 randomly-placed particles, which is fairly complex. Try the simplest thing possible -- 3 or 4 particles, very far apart. Then move them closer together so they actually interact. After that, add one particle at a time, with regular spacing between them. Then add randomization. At every stage, see what the maximum timestep is for your system to be stable. May 23 '15 at 20:00