subroutine SANN

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Fortran implementation of the SANN algorithm                        !
!                                                                         !
!     van Meel, Filion, Valeriani and Frenkel November (2011)             !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     


!     declare all variable used in the subroutine      
      implicit none
      
!     npart = total number of particles in the system
      integer npart
!     m = tentative number of neighbours
      integer i,j,k,m
!     countneighbors = number of neighbours of particle i
      integer countneighbors(1000)
!     neighbor = list of neighbours of particles i
      integer neighbor(1000,100)
!     sortneighbor = sorted neighbours 
      integer sortneighbor(1000,100)
!     selectedneighbors = list of selected neighbours    
      integer selectedneighbors(1000,100)
!     Nb = final number of neighbours of particle i
      integer Nb(1000)
!     edge of the simulation box
      double precision box
!     distance = list of distances between each 
!     neighbour of particle i and particle i 
      double precision  distance(1000,100)
!     distancesorted = sorted distances
      double precision distancesorted(1000,100)       
!     R(m) as in Eq.3 in the manuscript
      double precision rm,rm1
!     x,y,z component of every particle i
      double precision x(1000),y(1000),z(1000)
!     distance between particle i and particle j
      double precision dr
!     cutoff distance to identify all potential neighbours
      double precision rcutoff
      
      
!     Step 1:
!     first we identify the particles within a cutoff radius rcutoff
      do i=1,npart
!     loop over all particles different from i
         do j =1,npart
            if (j.ne.i)then
!     compute x,y,z component of the distance between particle i and j
               dx = x(j) - x(i)
               dy = y(j) - y(i)
               dz = z(j) - z(i)
!     applying periodic boundary conditions
               dx=dx-nint(dx/box)*box
               dy=dy-nint(dy/box)*box
               dz=dz-nint(dz/box)*box
!     compute distance dr  between particle i and j
               dr =  sqrt(dx*dx+dy*dy+dz*dz)
!     identify neighbours that are within a cutoff (rcutoff)
               if(dr.lt.rcutoff)then
!     j is a neighbour of i
                  countneighbors(i) = countneighbors(i) + 1
!     build a list of neighbours
                  neighbor(i,countneighbors(i))= j 
!     create a list with the distance between i and j 
                  distance(i,countneighbors(i))=dr
               endif
            endif
         enddo
      enddo	 
  
!     Step 2:
!     for every particle i sort all (countneighbors) 
!     neighbours (neighbor) according to their 
!     distances (distance) and create  a new list of 
!     particle i's (sortneighbor)
!     and a new sorted list of distances (distancesorted)
      do i=1,npart
         call sort(i,countneighbors,distance,neighbor,
     &        sortneighbor,distancesorted)
      enddo

      do i=1,npart
!     Step 3: 
!     start with 3 neighbours
         m = 3
!     Step 4: 
!     compute R(m) as in Eq.3 
         rm = 0
         do k=1,m
            rm = rm + distancesorted(i,k)
         enddo
         rm = rm/(m-2)
!     compute r(m+1)
         do j = 1,countneighbors(i)      
            rm1 = 0
            do k=1,m
               rm1 = rm1 + distancesorted(i,k)
            enddo
            rm1 = rm1/(m-2)
!     Step 5:  
!     if rm > rm1     
            if(rm.ge.rm1)then     
               rm = rm1
!     increase m
               m = m+1
            else
!     Step 6:
!     if rm < rm1, m is the final number of neighbours
               exit
            endif
         enddo
!     the final number of neighbours is m = Nb(i) 
!     and the neighbours are  selectedneighbors
         Nb(i) = m
         do j=1,Nb(i)
            selectedneighbors(i,j) = sortneighbor(i,j)
         enddo
      enddo
      
      
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!