!================================================================================!
! This file is part of crest.
!
! Copyright (C) 2017 Stefan Grimme
!
! crest is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! crest is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with crest.  If not, see <https://www.gnu.org/licenses/>.
!================================================================================!

c determine ring and ringsize from neighbor list
c SG 2017, for QMDFF
      subroutine getring(n,nbin,a0,c,iring)
      implicit none
      integer n,nbin(20,n),a0,iring,i,nb(20,0:n)
      integer    i1,i2,i3,i4,i5,i6,i7,i8,i9,i10
      integer n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
      integer    a1,a2,a3,a4,a5,a6,a7,a8,a9,a10
      integer c(10)
      logical chk
 
      if(n.le.2) return

      nb(:,1:n)=nbin(:,:)
      do i=1,n
         if(nb(20,i).eq.1)nb(20,i)=0
      enddo

      iring =0
      c     =0

      n0=nb(20,a0)

      do i1=1,n0
         a1=nb(i1,a0)
         if(a1.eq.a0) cycle
         n1=nb(20,a1)
         do i2=1,n1
            a2=nb(i2,a1)
            if(a2.eq.a1) cycle
            n2=nb(20,a2)
            do i3=1,n2
               a3=nb(i3,a2)
               n3=nb(20,a3)
               if(a3.eq.a2) cycle
               c(1)=a1
               c(2)=a2
               c(3)=a3
               if(a3.eq.a0.and.chk(n,3,c))then
               iring=3
               goto 99
               endif
               do i4=1,n3
                  a4=nb(i4,a3)
                  n4=nb(20,a4)
                  if(a4.eq.a3) cycle
                  c(4)=a4
                  if(a4.eq.a0.and.chk(n,4,c))then
                     iring=4
                     goto 99
                  endif
                  do i5=1,n4
                     a5=nb(i5,a4)
                     n5=nb(20,a5)
                     if(a5.eq.a4) cycle
                     c(5)=a5
                     if(a5.eq.a0.and.chk(n,5,c))then
                        iring=5
                        goto 99
                     endif
                     do i6=1,n5
                        a6=nb(i6,a5)
                        n6=nb(20,a6)
                        if(a6.eq.a5) cycle
                        c(6)=a6
                        if(a6.eq.a0.and.chk(n,6,c))then
                           iring=6
                           goto 99
                        endif
                        do i7=1,n6
                           a7=nb(i7,a6)
                           n7=nb(20,a7)
                           if(a7.eq.a6) cycle
                           c(7)=a7
                           if(a7.eq.a0.and.chk(n,7,c))then
                              iring=7
                              goto 99
                           endif
                           do i8=1,n7
                              a8=nb(i8,a7)
                              n8=nb(20,a8)
                              if(a8.eq.a7) cycle
                              c(8)=a8
                              if(a8.eq.a0.and.chk(n,8,c))then
                                 iring=8
                                 goto 99
                              endif
                              do i9=1,n8
                                 a9=nb(i9,a8)
                                 n9=nb(20,a9)
                                 if(a9.eq.a8) cycle
                                 c(9)=a9
                                 if(a9.eq.a0.and.chk(n,9,c))then
                                    iring=9
                                    goto 99
                                 endif
                                 do i10=1,n9
                                    a10=nb(i10,a9)
                                    n10=nb(20,a10)
                                    if(a10.eq.a9) cycle
                                    c(10)=a10
                                    if(a10.eq.a0.and.chk(n,10,c))then
                                       iring=10
                                       goto 99
                                    endif
                                 enddo
                              enddo
                           enddo
                        enddo
                     enddo
                  enddo
               enddo
            enddo
         enddo
      enddo

      return

 99   continue
c     write(*,*) 'atom ',a0,'ring found',c(1:iring)
c     call isort(iring,c)
      return

      end

      logical function chk(nn,n,c)
      implicit none
      integer n,idum(0:nn),nn,c(10),i,j
      chk=.true.
      idum=0
      do i=1,n
         idum(c(i))=idum(c(i))+1
      enddo
      j=0
      do i=1,nn
         if(idum(i).eq.1) j=j+1
      enddo
      if(j.ne.n) chk=.false.
      end


      logical function samering(n,i,j,c,s)
      implicit none
      integer n,i,j,c(10,n),s(n)
      integer m,k

      samering=.false.
      if(s(i).eq.0   ) return
      if(s(i).ne.s(j)) return

      m=0
      do k=1,s(i)
         m=m+c(k,i)-c(k,j)
      enddo
      if(m.eq.0) then
         samering=.true.
         return
      endif
      
      do k=1,s(i)
         if(c(k,i).eq.j)then
            samering=.true.
            return
         endif
      enddo
      do k=1,s(j)
         if(c(k,j).eq.i)then
            samering=.true.
            return
         endif
      enddo

      end

      subroutine neighborh(natoms,at,xyz,nb)
      implicit none  
      integer at(natoms),natoms,nb(20,natoms)
      real*8  xyz(3,natoms)

      logical da
      integer iat,i,j,k,nn
      real*8 rad(94)
      real*8 dx,dy,dz,r,rco,r2,f
      data rad/
     . 0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865,
     . 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527,
     . 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820,
     . 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730,
     . 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923,
     . 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188,
     . 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349,
     . 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216,
     . 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717,
     . 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967,
     . 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625,
     . 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657,
     . 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833,
     . 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098,
     . 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878,
     . 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790,
     . 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584,
     . 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289,
     . 3.82984466, 3.85504098, 3.88023730, 3.90543362 /

      nb=0
      nn=min(natoms,2)-1

      do i=1,natoms
      f=1.0
      k=0
 100  do iat=1,natoms
         da=.false.
         do j=1,k
            if(nb(j,i).eq.iat)da=.true.
         enddo
         if(iat.ne.i.and.(.not.da))then
            dx=xyz(1,iat)-xyz(1,i)
            dy=xyz(2,iat)-xyz(2,i)
            dz=xyz(3,iat)-xyz(3,i)
            r2=dx*dx+dy*dy+dz*dz 
            r=sqrt(r2)
            rco=rad(at(i))+rad(at(iat))
c critical step            
            if(r.lt.f*rco.and.k.lt.19)then
               k=k+1
               nb(k,i)=iat
            endif
         endif
      enddo
      if(k.lt.1.and.f.lt.1.5)then
         f=f*1.1
         goto 100
      endif
      nb(20,i)=k
      enddo

      end

      subroutine neighdist(natoms,at,xyz,nb,dist)
      implicit none  
      integer at(natoms),natoms,nb(200,natoms)
      real*8  xyz(3,natoms),dist(natoms,natoms)

      logical da
      integer iat,i,j,k,nn
      real*8 rad(94)
      real*8 dx,dy,dz,r,rco,r2,f
      data rad/
     . 0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865,
     . 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527,
     . 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820,
     . 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730,
     . 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923,
     . 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188,
     . 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349,
     . 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216,
     . 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717,
     . 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967,
     . 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625,
     . 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657,
     . 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833,
     . 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098,
     . 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878,
     . 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790,
     . 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584,
     . 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289,
     . 3.82984466, 3.85504098, 3.88023730, 3.90543362 /

      nb=0
      nn=min(natoms,2)-1

      do i=1,natoms
      f=1.0
      k=0
 100  do iat=1,natoms
         da=.false.
         do j=1,k
            if(nb(j,i).eq.iat)da=.true.
         enddo
         dx=xyz(1,iat)-xyz(1,i)
         dy=xyz(2,iat)-xyz(2,i)
         dz=xyz(3,iat)-xyz(3,i)
         r2=dx*dx+dy*dy+dz*dz 
         r=sqrt(r2)
         dist(iat,i)=r
         if(iat.ne.i.and.(.not.da))then
            rco=rad(at(i))+rad(at(iat))
c critical step            
            if(r.lt.f*rco.and.k.lt.199)then
               k=k+1
               nb(k,i)=iat
            endif
         endif
      enddo
      if(k.lt.1.and.f.lt.1.5)then
         f=f*1.1
         goto 100
      endif
      nb(200,i)=k
      enddo

      end
