!     Neighbour list for all atoms
!     There are several schemes for building the Neighbour list.
!     The one used here is deiscribed in Frenkel & Smit pp. 363-367
!     It is NOT memory efficient but easy to understand and better
!     suited for parallel algorithms.
!
!     In the method described in Frenkel & Smit pp. 363-363 they fill
!     all the neighbors of both atoms. For pair potentials, half the
!     list is sufficient.  The array NNNG, however, is less than half
!     full in this case.  For some many-body potentials, such as S-W
!     for Si, and in Monte Carlo calculations, the full neighbor list
!     is needed.  In this code logical variable FullList is used to
!     choose between halfthe list and full list calculation.
!
!     Box filtering is used.  We evaluate dz first and see if it is
!     larger than RList.  Then dx and dy are checked and only then r^2
!     is checked against RList2.  Box filtering is more efficient when
!     the system size is large or one or two dimensions are larger
!     than the others (here we assume that z dimension is the largest.
!
!     Other schemes:
!     Method with a one-dimensional array.  It is more memory efficient
!     but probably less transparent AND, it cannot be made parallel.
!
!     Linked cell method, Frenkel & Smit pp. 368-371, is better/faster
!     for really large systems.
!
!     MSE 6270, Leonid Zhigilei
!**********************************************************************

      SUBROUTINE NbList()
      INCLUDE 'common.h'
      real*8 wcl_st, wcl_end, cpu_st, cpu_end

      !SMP$ SCHEDULE (affinity)
!     call cpu_time(cpu_st)
!     cpu_st=mclock()
!     wcl_st=rtc()

      nc=0
      NNG(1:NAN)= 0  ! contains the number of neighbors of particle i

      !SMP$ ASSERT(NODEPS)
      !SMP$ PARALLEL DO &
      !SMP$&   PRIVATE(I,I3,IP1,J,J3,K,DX,DY,DZ,RR), &
      !SMP$&   SHARED(RList,RList2,LIDX,LIDY,LIDZ,XL,YL,ZL, &
      !SMP$&     XLHALF,YLHALF,ZLHALF,R,X,NNNG,NNG,NAN), &
      !SMP$&   REDUCTION(NC)
      outer_loop: DO I=1,NAN
        ktypei=KTYPE(I)
        IP1=I+1
        inner_loop: DO J=IP1,NAN
          IJ = ijindex(KTYPEI, KTYPE(J))

          DZ=XD(3,I)-XD(3,J)
          IF(LIDZ.EQ.1) THEN
            IF(DZ.GT.ZLHALF) DZ=DZ-ZL
            IF(DZ.LT.-ZLHALF) DZ=DZ+ZL
          ENDIF

          IF (ABS(DZ).LT.RList(IJ)) THEN
            DX=XD(1,I)-XD(1,J)
            IF(LIDX.EQ.1) THEN
              IF(DX.GT.XLHALF) DX=DX-XL
              IF(DX.LT.-XLHALF) DX=DX+XL
            ENDIF

            IF (ABS(DX).LT.RList(IJ)) THEN
              DY=XD(2,I)-XD(2,J)
              IF(LIDY.EQ.1) THEN
                IF(DY.GT.YLHALF) DY=DY-YL
                IF(DY.LT.-YLHALF) DY=DY+YL
              ENDIF

              IF(ABS(DY).LT.RList(IJ)) THEN
                RR=DX*DX+DY*DY+DZ*DZ
                IF (RR.LT.RList2(IJ)) THEN
!                 We've found a neighbor close enough to be on the list
                  NNG(I)=NNG(I)+1  ! increase the count of neighbors
                  NNNG(NNG(I),I)=J      ! add the neighbor of atom I

!                 Many-body potentials (e.g. SW,EAM) need the complete list
                  full:  IF (FullList) THEN
                    NNG(J)=NNG(J)+1
                    NNNG(NNG(J),J)=I
                  END IF full

                END IF
              END IF
            END IF
          END IF
        END DO inner_loop
        IF (NNG(I).gt.MAXNNB) THEN
          print *, ' TOO MANY NEIGHBOURS: ',NNG(I),' > ', maxnnb
          STOP
        ENDIF
        nc=nc+nng(i)
      END DO outer_loop
!     print 55, NC
!     WRITE(17,55) NC
! 55  FORMAT(5X,'  ***** Neighbour list ----> NC=',I8)
!     call cpu_time(cpu_end)
!     cpu_end=mclock()
!     wcl_end=rtc()
!     write(17,*) ' Wall clock used in sos3: ',wcl_end-wcl_st,' sec'
!     write(17,*) ' CPU time  used in sos3: ',cpu_end-cpu_st,' sec'
      RETURN
      END