program verlet
      implicit none
      integer nsteps,i
      real*8 omega,tau,fc,ekin,epot,etot,fc2,f1
      real*8 x0,x1,x,v,force,potential
      external force,potential
      common/ff/omega,tau
 
      write(*,*) 'input omega, x0, v0'
      read (*,*) omega,x,v
      write(*,*) 'input tau, nsteps'
      read (*,*) tau,nsteps
 
      fc=0.5d0*tau
      fc2=fc*tau
      ekin=0.5d0*v**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
       f1=force(x)
       x=x+tau*v+fc2*f1
       v=v+fc*(f1+force(x))
       ekin=0.5d0*v**2
       epot=potential(x)
       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,x,v
       endif
      enddo
      close(1)
      stop
      end
 
      function force(x)
      implicit none
      real*8 force,x
      real*8 omega,tau
      common/ff/omega,tau
      force=-omega**2*x!-0.2*omega**2*x**2
!     force=-omega**2*8.d0*x*(x**2-1)
      return
      end
 
      function potential(x)
      implicit none
      real*8 potential,x
      real*8 omega,tau
      common/ff/omega,tau
      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