Subroutine forc(forces,npart,ndim,pos,cell,nvec,reall,dx,options,rcut,nimage)

  !Calcule les forces
  !variables
  !eneropt: modele de calcul de l'energie
  !forces(npart,ndim): forces dur les particules
  !npart: nombre de particules
  !ndim: dimension de l'espace
  !pos(npart,ndim):position des particules
  !cell(nvec,ndim): cellule de simulation
  !nvec: nombre de vecteurs de la cellule
  !reall: cote de la maille
  !input
  USE types_def
  IMPLICIT NONE
  TYPE(options_type) :: options
  integer :: ndim,npart,nvec,nimage
  double precision :: reall,dx,rcut
  double precision :: forces(npart,ndim),pos(npart,ndim),cell(nvec,ndim)
  !local variables
  integer :: idim,ipart,jpart,lpart
  double precision, allocatable :: icord(:),jcord(:),forcesij(:),grijl(:),grjli(:),grlij(:),&
&                                  g1rijl(:),g2rijl(:),g3rijl(:),&
&                                  derdens(:),lcord(:),unitlj(:),unitli(:),unitij(:)
  double precision :: densij,densel,derembedd,distij,distil,distjl,screenij,densjl,densil,&
&                     d1gscreen,d2gscreen,d3gscreen
  allocate(icord(ndim),jcord(ndim),forcesij(ndim),grijl(ndim),grjli(ndim),grlij(ndim),&
&          g1rijl(ndim),g2rijl(ndim),g3rijl(ndim),&
&          derdens(npart),lcord(ndim),unitlj(ndim),unitli(ndim),unitij(ndim))
  !partie de paires
  if(options%forcopt=='analytique') then
     forces=0.d0
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
           if(jpart>ipart) then
              call pairforce(forcesij,ndim,icord,jcord,reall,rcut,nimage,options)
              forces(ipart,:)=forces(ipart,:)+forcesij(:)
              forces(jpart,:)=forces(jpart,:)-forcesij(:)
           endif
        enddo
     enddo
     !partie embedding
     derdens=0.d0
     do ipart=1,npart
        densel=0.d0
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
              call pairdens(densij,options%eneropt,icord,jcord,ndim,reall,pos,npart,rcut,nimage,&
&                           options)
              densel=densel+densij
        enddo
        derdens(ipart)=derembedd(densel,options)
     enddo
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
           if(jpart>ipart) then
              call densforce(forcesij,ndim,icord,jcord,reall,rcut,nimage,options)
              call screen(screenij,icord,jcord,ndim,reall,pos,npart)
              forces(ipart,:)=forces(ipart,:)+forcesij(:)*screenij*(derdens(ipart)+derdens(jpart))
              forces(jpart,:)=forces(jpart,:)-forcesij(:)*screenij*(derdens(ipart)+derdens(jpart))
           endif
        enddo
     enddo
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
           call pairdens(densij,options%eneropt,icord,jcord,ndim,reall,pos,npart,rcut,nimage,&
&                        options)
           call dist(distij,icord,jcord,ndim)
           do lpart=1,npart
              if(lpart==jpart) cycle
              if(lpart==ipart) cycle
              lcord(:)=pos(lpart,:)
              call pairdens(densjl,options%eneropt,jcord,lcord,ndim,reall,pos,npart,rcut,nimage,&
&                        options)
              call pairdens(densil,options%eneropt,icord,lcord,ndim,reall,pos,npart,rcut,nimage,&
&                        options)
              call dist(distil,icord,lcord,ndim)
              call dist(distjl,jcord,lcord,ndim)
              unitli(:)=(lcord(:)-icord(:))/distil
              unitlj(:)=(lcord(:)-jcord(:))/distjl
              grijl(:)=unitlj(:)*d3gscreen(distij,distil,distjl,reall)+&
&                      unitli(:)*d2gscreen(distij,distil,distjl,reall)
              grjli(:)=unitlj(:)*d1gscreen(distjl,distij,distil,reall)+&
&                      unitli(:)*d3gscreen(distjl,distij,distil,reall)
              grlij(:)=unitli(:)*d1gscreen(distil,distjl,distij,reall)+&
&                      unitlj(:)*d2gscreen(distil,distjl,distij,reall)
              forces(lpart,:)=forces(lpart,:)+derdens(ipart)*densij*grijl(:)+&
&                                             derdens(jpart)*densjl*grjli(:)+&
&                                             derdens(lpart)*densil*grlij(:)
          enddo         
       enddo
    enddo     
  else if (options%forcopt=='rapide') then
     forces=0.d0
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
           if(jpart>ipart) then
              call pairforce(forcesij,ndim,icord,jcord,reall,rcut,nimage,options)
              forces(ipart,:)=forces(ipart,:)+forcesij(:)
              forces(jpart,:)=forces(jpart,:)-forcesij(:)
           endif
        enddo
     enddo
     !partie embedding
     derdens=0.d0
     do ipart=1,npart
        densel=0.d0
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
              call pairdens(densij,options%eneropt,icord,jcord,ndim,reall,pos,npart,rcut,nimage,&
&                        options)
              densel=densel+densij
        enddo
        derdens(ipart)=derembedd(densel,options)
     enddo
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           jcord(:)=pos(jpart,:)
           if(jpart>ipart) then
              call densforce(forcesij,ndim,icord,jcord,reall,rcut,nimage,options)
              screenij=1.0d0
              if(options%eneropt=='cai')call screen(screenij,icord,jcord,ndim,reall,pos,npart)
              forces(ipart,:)=forces(ipart,:)+forcesij(:)*screenij*(derdens(ipart)+derdens(jpart))
              forces(jpart,:)=forces(jpart,:)-forcesij(:)*screenij*(derdens(ipart)+derdens(jpart))
           endif
        enddo
     enddo
     if(options%eneropt=='cai') then
     do ipart=1,npart
        icord(:)=pos(ipart,:)
        do jpart=1,npart
           if(ipart==jpart) cycle
           jcord(:)=pos(jpart,:)
           call pairdens(densij,options%eneropt,icord,jcord,ndim,reall,pos,npart,rcut,nimage,&
&                        options)
           call dist(distij,icord,jcord,ndim)
           unitij(:)=(icord(:)-jcord(:))/distij
              do lpart=1,npart
              if(lpart==jpart) cycle
              if(lpart==ipart) cycle
              lcord(:)=pos(lpart,:)
!              call dist(distil,icord,lcord,ndim)
!              call dist(distjl,jcord,lcord,ndim)
              distil=0.d0
              do idim=1,ndim;distil=distil+(lcord(idim)-icord(idim))**2;enddo
                 distil=sqrt(distil)
              distjl=0.d0
              do idim=1,ndim;distjl=distjl+(lcord(idim)-jcord(idim))**2;enddo
                 distjl=sqrt(distjl)
              unitli(:)=(lcord(:)-icord(:))/distil
              unitlj(:)=(lcord(:)-jcord(:))/distjl
              g1rijl(:)= unitij(:)*d1gscreen(distij,distil,distjl,reall)-&
&                        unitli(:)*d2gscreen(distij,distil,distjl,reall)
              g2rijl(:)=-unitij(:)*d1gscreen(distij,distil,distjl,reall)-&
&                        unitlj(:)*d3gscreen(distij,distil,distjl,reall)
              g3rijl(:)= unitlj(:)*d3gscreen(distij,distil,distjl,reall)+&
&                        unitli(:)*d2gscreen(distij,distil,distjl,reall)
              forces(ipart,:)=forces(ipart,:)+derdens(ipart)*densij*g1rijl(:)
              forces(jpart,:)=forces(jpart,:)+derdens(ipart)*densij*g2rijl(:)
              forces(lpart,:)=forces(lpart,:)+derdens(ipart)*densij*g3rijl(:)
          enddo         
       enddo
    enddo 
    endif    
  endif
  return
end subroutine forc

subroutine pairforce(force,ndim,icord,jcord,reall,rcut,nimage,options)
  !Calcule les forces de paire entre icord et jcord (-deij/dri)
  !variables
  !force(ndim): forces sur les particules
  !icord,jcord: les coordonnees de la paire d'atomes
  !ndim: dimension de l'espace
  !reall: cote de la maille
  USE types_def
  IMPLICIT NONE
  TYPE(options_type) :: options
  !input
  integer :: ndim,nimage
  double precision :: force(ndim),icord(ndim),jcord(ndim)
  double precision :: reall,rcut
  !local variables
  double precision :: dist,vprim
  integer :: idim,iimage,ishift(ndim),div
  force(:)=0.0d0
  do iimage=0,(2*nimage+1)**3-1
     ishift(1)=mod(iimage,2*nimage+1); div=(iimage-ishift(1))/(2*nimage+1)
     ishift(2)=mod(div,2*nimage+1); div=(div-ishift(2))/(2*nimage+1)
     ishift(3)=div
     ishift(:)=ishift(:)-nimage
     dist=0.0d0
!     ishift(:)=0
     do idim=1,ndim
        dist=dist+(icord(idim)-jcord(idim)-dfloat(ishift(idim)))**2
     enddo
     dist=sqrt(dist)
     force(:)=force(:)-(icord(:)-jcord(:)-dfloat(ishift(:)))*vprim(dist,reall,rcut,options)/dist
!     return
 enddo
  return
end subroutine pairforce



subroutine densforce(force,ndim,icord,jcord,reall,rcut,nimage,options)
  !Calcule les forces de paire entre icord et jcord (-dfij/dri)
  !dues a la derivee de la densite
  !variables
  !force(ndim): forces sur les particules
  !icord,jcord: les coordonnees de la paire d'atomes
  !ndim: dimension de l'espace
  !reall: cote de la maille
  USE types_def
  IMPLICIT NONE
  TYPE(options_type) :: options
  !input
  integer :: ndim, nimage
  double precision :: force(ndim),icord(ndim),jcord(ndim)
  double precision :: reall,rcut
  !local variables
  double precision :: dist,dereldens
  integer :: idim,iimage,ishift(ndim),div
  force(:)=0.0d0
  do iimage=0,(2*nimage+1)**3-1
     ishift(1)=mod(iimage,2*nimage+1); div=(iimage-ishift(1))/(2*nimage+1)
     ishift(2)=mod(div,2*nimage+1); div=(div-ishift(2))/(2*nimage+1)
     ishift(3)=div
     ishift(:)=ishift(:)-nimage
     dist=0.0d0
     do idim=1,ndim
        dist=dist+(icord(idim)-jcord(idim)-dfloat(ishift(idim)))**2
     enddo
     dist=sqrt(dist)
     force(:)=force(:)-(icord(:)-jcord(:)-dfloat(ishift(:)))*dereldens(dist,reall,rcut,&
&                       options)/dist
  enddo
  return
end subroutine densforce



subroutine dist(distij,icord,jcord,ndim)
  !Calcule la distance entre  icord et jcord
  !variables
  !distij: distance entre les particules
  !icord,jcord: les coordonnees de la paire d'atomes
  !ndim: dimension de l'espace
  USE types_def
  IMPLICIT NONE
  !input
  integer :: ndim
  double precision :: icord(ndim),jcord(ndim)
  double precision :: distij
  !local variables
  integer :: idim
  distij=0.0d0
  do idim=1,ndim
     distij=distij+(icord(idim)-jcord(idim))**2
  enddo
  distij=sqrt(distij)
  return
end subroutine dist

subroutine unit(unitij,icord,jcord,ndim)
  !Calcule le vecteur unitaire entre  icord-jcord
  !variables
  !unitij: distance entre les particules
  !icord,jcord: les coordonnees de la paire d'atomes
  !ndim: dimension de l'espace
  USE types_def
  IMPLICIT NONE
  TYPE(options_type) :: options
  !input
  integer :: ndim
  double precision :: icord(ndim),jcord(ndim),unitij(ndim)
  double precision :: distij
  !local variables
  call dist(distij,icord,jcord,ndim)
  unitij(:)=(icord(:)-jcord(:))/distij
  return
end subroutine unit