!====================================================================
! This is a fortran code to integrate the evolution equation of an
! harmonic oscillator using the velocity Verlet algorithm
!====================================================================
      program MDHarm
      implicit none
      include 'MD.cm'

      integer istep,nstep,ndump    

      real*8 r,v,xt,vt
      real*8 avgr, avgv
      real*8 epot, ekin, etot 
      real*8 dt, time, derpot

      open(170,file='Energy',status='unknown')
      open(171,file='Traj',status='unknown')
      open(172,file='Averages',status='unknown')
      open(273,file='ETraj',status='unknown')

      mass=1.
      omega1=1.
      omega2=7.

      dt=0.0125
      nstep=40000
      ndump=400

      avgr=0.d0
      avgv=0.d0

      call SetInCon(r,v)
        ekin=0.5*mass*v**2
        epot=0.5*mass*(omega1**2+omega2**2)*r**2
        etot=ekin+epot   
        write(170,555) 0.,epot,ekin,etot
        write(171,555) 0.,r,v 




!====================================================
! Use Velocity Verlet Algorithm
!====================================================
      call GetForce(r,derpot)

      avgr=avgr+r
      avgv=avgv+v  

      do istep=1,nstep

        v=v-0.5*dt*derpot/mass
        r=r+dt*v
        call GetForce(r,derpot)
        v=v-0.5*dt*derpot/mass

        avgr=avgr+r
        avgv=avgv+v

        if(mod(istep,ndump).eq.1) then
          ekin=0.5*mass*v**2
          epot=0.5*mass*(omega1**2+omega2**2)*r**2
          etot=ekin+epot
          write(170,555) istep*dt,epot,ekin,etot
          write(171,555) istep*dt,r,v
          write(172,555) istep*dt,avgr/real(istep),avgv/real(istep)
        end if
      enddo 
     

555   format(7E20.10)

      close(170)
      close(171)
      close(172)
      close(273) 
   

      end program MDHarm  
!==================================================================
      subroutine SetInCon(r0,v0)
!~~~~~~~~~~~~~~~~~~~~~~~~~~
! Set initial conditions
!~~~~~~~~~~~~~~~~~~~~~~~~~~
      implicit none
      real*8 r0,v0

        r0=0.
        v0=1.

      return 
      end subroutine SetInCon


!=================================================================
      subroutine GetForce(x,derp)
      implicit none
      include 'MD.cm'

      real*8 x,derp

      derp=mass*(omega1**2+omega2**2)*x

      return 
      end subroutine GetForce