!-------------------------------------------------------------!
! Main program
Program twoparticles
! Velocity verlet is used to solve the initial value problem
! to find the classical motion of two particles interacting
! through an Lennard-Jones potential

Implicit None
Real*8 :: x(2,3),v(2,3),xo(2,3),vo(2,3),f(2,3),fo(2,3),dt,pot,Energy ! double precision
Integer :: i,j,k,l,nit

! Input data
Write(*,*)" Input the number of timesteps"
Read(*,*) nit
Write(*,*) nit
dt=0.004			! timestep

! Initialize positions, forces and velocities
xo=0.0				! set all initial positions to zero					
xo(1,1) = -1.0			! set x-coordinate of first particle
xo(2,1) = 1.0			! set x-coordinate of second particle 
Call Force(2,xo,fo)		! call force subroutine to find initial forces
vo=0.0				! set all initial velocities to zero

Do i=1,nit
	
	! velocity Verlet updates
	x = xo + dt*vo + dt**2*fo/2.0	! update positions
	Call Force(2,x,f)		! update forces
	v = vo + dt*(f+fo)/2.0		! update velocities

	! reset forces, velocities, positions
	fo=f ; vo=v ; xo=x		
	
	call Output(i,2,dt,x,v)
	
	End Do

End Program twoparticles

!---------------------------------------------------------------!
! Force subroutine
Subroutine Force(n,xs,fs)
Implicit None
Real*8 :: xs(2,3),fs(2,3),r2 ! double precision
Integer :: i,n

r2 = (xs(1,1)-xs(2,1))**2. + (xs(1,2)-xs(2,2))**2. + (xs(1,3)-xs(2,3))**2.
fs(1,:) =  -24.0*(2./r2**7. - 1./r2**4.)*(xs(2,:)-xs(1,:))
fs(2,:)=  - fs(1,:)
! write (*,*) fs

End Subroutine

!----------------------------------------------------------------!
! Potential energy subroutine
Subroutine Potential(n,xs,pots)
Implicit None
Real*8 :: xs(2,3),r2,pots ! double precision
Integer :: i,n

r2 = (xs(1,1)-xs(2,1))**2. + (xs(1,2)-xs(2,2))**2. + (xs(1,3)-xs(2,3))**2.
pots =  4.0*(1./r2**6. - 1./r2**3.)

End Subroutine

!-----------------------------------------------------------------!
! Output subroutine
Subroutine Output(i,n,dts,xs,vs)
Implicit None
Real*8 :: xs(2,3),vs(2,3),pot,dts ! double precision
Integer :: i,n
Open(10,file="mdout1.dat")
Open(11,file="mdout2.dat")
Open(12,file="mdout3.dat")
Open(13,file="mdout4.dat")

Write(10,*) dts*float(i), xs(1,1)
Write(11,*) dts*float(i), xs(2,1)
Write(12,*) dts*float(i), vs(1,1)
Call Potential(2,xs,pot)		      ! Call potential energy subroutine
Write(13,*) dts*float(i), pot + vs(1,1)**2.0  ! Total energy - should be constant

End Subroutine