program verlet
      implicit none
      integer nsteps,i
      real*8 omega,tau,ot2,fc,ekin,epot,etot
      real*8 x0,x1,x,v0,force,potential
      external force,potential
      common/ff/omega,tau,ot2
 
      write(*,*) 'input omega, x0, v0'
      read (*,*) omega,x,v0
      write(*,*) 'input tau, nsteps'
      read (*,*) tau,nsteps
 
      ot2=omega**2*tau**2
      fc=0.5d0/tau
      x1=x-tau*v0+0.5d0*force(x)   ! x(-1) O(tau**3)
      ekin=0.5d0*v0**2
      epot=potential(x)
      etot=ekin+epot
      i=0

      open(1,file='traiettoria.dat',form='formatted')
      open(2,file='verlet.out',form='formatted')
      write(*,'(i7,3g20.10)') i,etot,ekin,epot
      write(2,'(i7,3g20.10)') i,etot,ekin,epot
      do i=1,nsteps
       x0=x
       x=2.d0*x0-x1+force(x0)
       v0=fc*(x-x1)
       x1=x0
       ekin=0.5d0*v0**2
       epot=potential(x0)
       etot=ekin+epot
       if(mod(i,10).eq.0) then
        write(*,'(i7,3g20.10)') i,etot,ekin,epot
        write(2,'(i7,3g20.10)') i,etot,ekin,epot
        write(1,'(i7,2g20.10)') i,x0,v0
       endif
      enddo
      close(1)
      stop
      end
 
      function force(x)
      implicit none
      real*8 force,x
      real*8 omega,tau,ot2
      common/ff/omega,tau,ot2
      force=-ot2*x!-0.2*ot2*x**2
!     force=-ot2*8.d0*x*(x**2-1)
      return
      end
 
      function potential(x)
      implicit none
      real*8 potential,x
      real*8 omega,tau,ot2
      common/ff/omega,tau,ot2
      potential=0.5d0*(omega*x)**2!+0.2*omega**2*x**3/3.d0
!     potential=2.d0*(1.d0+x**2*(x**2-2.d0))*omega**2
      return
      end