subroutine testforce(ameam,pos,npart,ndim,dx,cell,nvec,reall,options,rcut,nimage)

  !Teste les forces pqr derivee directe de l'energie
  !variables
  !ameam(2): parametres a du modele meam
  !eneropt: modele de calcul de l'energie
  !dx: deplacement
  !npart: nombre de particules
  !ndim: dimension de l'espace
  !pos(npart,ndim):position des particules
  !cell: cellule de simulation
  !nvec: nb de vecteurs dans la cellule
  !reall:longueur unites reeles
  !diff: erreur relative entre la force analytique et la force -de/dx
  !enerp, enerm: energie apres un deplacement de +-dx
  !input
  USE types_def
  IMPLICIT NONE
  TYPE(options_type) :: options
  integer :: ndim,npart,nvec,nimage
  double precision :: dx,reall,rcut
  double precision :: ameam(2),pos(npart,ndim),cell(nvec,ndim)
  !output

  !local variables
  double precision, allocatable :: forces(:,:)
  double precision :: enerm,enerp,diff 
  integer :: ipart,idim
  allocate(forces(npart,ndim))
  print*,'enter testforce'
  Call forc(forces,npart,ndim,pos,cell,nvec,reall,dx,options,rcut,nimage)
  do ipart=1,npart
     do idim=1,ndim
        pos(ipart,idim)=pos(ipart,idim)+dx
        call enerpot(ameam,enerp,options,pos,npart,ndim,cell,nvec,reall,rcut,nimage)
        pos(ipart,idim)=pos(ipart,idim)-2.d0*dx
        call enerpot(ameam,enerm,options,pos,npart,ndim,cell,nvec,reall,rcut,nimage)
        pos(ipart,idim)=pos(ipart,idim)+dx
        diff=forces(ipart,idim)+(enerp-enerm)/2.d0/dx
        print*,diff/enerp,forces(ipart,idim)/reall
     enddo
  enddo
  deallocate(forces)
end subroutine testforce