rings-primitive.F90

Go to the documentation of this file.

! This file is part of the 'atomes' software.
!
! 'atomes' is free software: you can redistribute it and/or modify it under the terms
! of the GNU Affero General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! 'atomes' 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 General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License along with 'atomes'.
! If not, see <https://www.gnu.org/licenses/>
!
! Copyright (C) 2022-2024 by CNRS and University of Strasbourg
!
 
INTEGER FUNCTION primitive_rings ()
 
!
! Primitive and strong ring statistics
!
#ifdef OPENMP
!$ USE OMP_LIB
#endif
USE parameters
 
IMPLICIT NONE
 
INTEGER :: rid
#ifdef OPENMP
INTEGER :: numth
LOGICAL :: doatoms=.false.
#endif
 
INTERFACE
  INTEGER FUNCTION recrings(VID)
    INTEGER, INTENT(IN) :: vid
  END FUNCTION
END INTERFACE
 
primitive_rings=0
 
! Dynamic allocation of pointers and tables used in the "RINGS" subroutine
 
rid = 3
if (calc_strings) rid = 4
 
#ifdef OPENMP
numth = omp_get_max_threads()
doatoms=.false.
! Dynamic allocation of pointers and tables used in the "RINGS" subroutine
if (ns.ge.1 .and. ns.lt.numth) then
  if (numth .ge. 2*(ns-1)) then
    doatoms=.true.
  else
    numth=ns
  endif
endif
 
if (all_atoms) doatoms=.true.
 
if (doatoms) then
  if (na.lt.numth) numth=na
#ifdef DEBUG
  write (6, *) "OpenMP on atoms, NUMTH= ",numth
#endif
  call primitive_ring_search_atoms (rid, numth)
else
#ifdef DEBUG
  write (6, *) "OpenMP on MD steps, NUMTH= ",numth
#endif
  call primitive_ring_search_steps (rid, numth)
endif
#else
call primitive_ring_search_steps (rid)
#endif
 
primitive_rings = recrings(rid)
 
END FUNCTION
 
#ifdef OPENMP
SUBROUTINE primitive_ring_search_atoms (RID, NUMTH)
 
USE parameters
!$ USE OMP_LIB
IMPLICIT NONE
INTEGER, INTENT(IN) :: numth
INTEGER, INTENT(IN) :: rid
INTEGER, DIMENSION(:), ALLOCATABLE :: tring, indte
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: savring, ordring
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: savr, ordr
INTEGER :: ri
LOGICAL, DIMENSION(2) :: fndtab
INTERFACE
  INTEGER FUNCTION rings_to_ogl_menu (IDSEARCH, NRI)
    USE parameters
    INTEGER, INTENT(IN) :: idsearch
    INTEGER, DIMENSION(TAILLR, NS), INTENT(IN) :: nri
  END FUNCTION
  INTEGER FUNCTION rings_to_ogl (STEP, IDSEARCH, NRI, RSAVED, OSAVED)
    USE parameters
    INTEGER, INTENT(IN) :: step, idsearch
    INTEGER, DIMENSION(TAILLR,NS), INTENT(IN) :: nri
    INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(IN) :: rsaved, osaved
  END FUNCTION
END INTERFACE
 
ri = 0
if(allocated(savring)) deallocate(savring)
allocate(savring(taillr,numa,taillr), stat=err)
if (err .ne. 0) then
  alc_tab="SAVRING"
  alc=.true.
  goto 001
endif
if(allocated(ordring)) deallocate(ordring)
allocate(ordring(taillr,numa,taillr), stat=err)
if (err .ne. 0) then
  alc_tab="ORDRING"
  alc=.true.
  goto 001
endif
if(allocated(cpat)) deallocate(cpat)
allocate(cpat(nna), stat=err)
if (err .ne. 0) then
  alc_tab="CPAT"
  alc=.true.
  goto 001
endif
 if(allocated(vpat)) deallocate(vpat)
allocate(vpat(nna,maxn), stat=err)
if (err .ne. 0) then
  alc_tab="VPAT"
  alc=.true.
  goto 001
endif
 
do i=1, ns
 
  savring(:,:,:)=0
  ordring(:,:,:)=0
  call setup_cpat_vpat_ring (nna, i, contj, voisj, cpat, vpat)
  ! OpenMP on atoms only
  !$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
  !$OMP& PRIVATE(FNDTAB, MAXAT, MINAT, SAUT, PATH, PATHOUT, &
  !$OMP& h, j, k, l, m, n, o, p, INDTE, APNA, RES_LIST, &
  !$OMP& ERR, TRING, SAVR, ORDR, PRINGORD, NPRING, MATDIST, QUEUE, QUERNG) &
  !$OMP& SHARED(NUMTH, i, RID, CALC_STRINGS, NS, NA, NNA, NNP, TLT, NSP, LOT, TAILLR, CONTJ, VOISJ, &
  !$OMP& NUMA, MAXPNA, MINPNA, ABAB, NO_HOMO, TBR, ALC, ALC_TAB, SAVRING, ORDRING, CPAT, VPAT, &
  !$OMP& NCELLS, THE_BOX, FULLPOS, PBC, MAXN, NRING, INDRING, PNA, ri)
  if(allocated(matdist)) deallocate(matdist)
  allocate(matdist(nna), stat=err)
  if (err .ne. 0) then
    alc_tab="MATDIST"
    alc=.true.
    goto 002
  endif
  if(allocated(queue)) deallocate(queue)
  allocate(queue(nna), stat=err)
  if (err .ne. 0) then
    alc_tab="QUEUE"
    alc=.true.
    goto 002
  endif
  if (allocated(pringord)) deallocate(pringord)
  allocate(pringord(numa*10,taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="PRINGORD"
    alc=.false.
    goto 002
  endif
  if (allocated(npring)) deallocate(npring)
  allocate(npring(nna), stat=err)
  if (err .ne. 0) then
    alc_tab="NPRING"
    alc=.false.
    goto 002
  endif
  if(allocated(res_list)) deallocate(res_list)
  allocate(res_list(taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="RES_LIST"
    alc=.true.
    goto 002
  endif
  if(allocated(savr)) deallocate(savr)
  allocate(savr(taillr,numa,taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="SAVR"
    alc=.true.
    goto 002
  endif
  if(allocated(ordr)) deallocate(ordr)
  allocate(ordr(taillr,numa,taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="ORDR"
    alc=.true.
    goto 002
  endif
  if(allocated(indte)) deallocate(indte)
  allocate(indte(numa), stat=err)
  if (err .ne. 0) then
    alc_tab="INDTE"
    goto 002
  endif
  if (allocated(tring)) deallocate(tring)
  allocate(tring(taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="TRING"
    alc=.true.
    goto 002
  endif
  if(allocated(apna)) deallocate(apna)
  allocate(apna(taillr), stat=err)
  if (err .ne. 0) then
    alc_tab="APNA"
    alc=.true.
    goto 002
  endif
  savr(:,:,:)=0
  ordr(:,:,:)=0
  tring(:)=0
  !$OMP DO SCHEDULE(STATIC,NA/NUMTH)
  do j=nnp+1, nnp+na ! atoms-loop
 
    if (tbr .or. alc) goto 003
    if (tlt .eq. nsp+1 .or. lot(j-nnp) .eq. tlt) then
 
      apna(:)=0
      minat=taillr
      maxat=1
      saut=.true.
 
      call dijkstra (j, cpat, vpat, queue, matdist)
      if (tbr .or. alc) goto 003
 
      do k=1, taillr/2 + mod(taillr,2) ! ring-sizes-loop
 
        indte(:)=0
        res_list(:)=0
        path=0
        do l=1, nna
          if (matdist(l) .eq. k) then
            call spath_rec (path,l,k,k,matdist,cpat,vpat,npring,pringord)
          endif
        enddo
        h = path*(path-1)/2
        if (allocated(querng)) deallocate(querng)
        allocate(querng(h,2), stat=err)
        if (err .ne. 0) then
          alc_tab="QUERNG"
          alc=.true.
          goto 003
        endif
 
        l=0
        do m=1, path-1
          do n=m+1, path
            l=l+1
            querng(l,1)=m
            querng(l,2)=n
!           Paths which share atoms ... so with overlap ... are deleted
            pathout=.false.
            do o=1, k-1
              do p=1, k-1
                if (pringord(m,o) .eq. pringord(n,p)) then
                  pathout=.true.
                  exit
                endif
              enddo
              if (pathout) exit
            enddo
            if (pathout) then
              querng(l,1)=0
              querng(l,2)=0
              l=l-1
            endif
          enddo
        enddo
        fndtab(:)=.false.
        call prim_ring (fndtab, j, l, k, h, cpat, vpat, querng, pringord, matdist, &
                        savr, ordr, tring, indte, res_list)
        if (tbr .or. alc) goto 003
        m = 2*k
        if (fndtab(1)) then
          if (apna(m).eq.0) then
            apna(m)=1
            minat=min(minat,m)
            maxat=max(maxat,m)
            saut=.false.
          endif
        endif
        if (fndtab(2)) then
          m = m + 1
          if (apna(m).eq.0) then
            apna(m)=1
            minat=min(minat,m)
            maxat=max(maxat,m)
            saut=.false.
          endif
        endif
        if (allocated(querng)) deallocate(querng)
 
      enddo ! end ring-sizes-loop
 
      if (.not. saut) then
        do k=3, taillr
          do l=3, taillr
            if (apna(k).eq.1 .and. apna(l).eq.1) then
              !$OMP ATOMIC
              pna(k,l,i)=pna(k,l,i)+1
            endif
          enddo
        enddo
        !$OMP ATOMIC
        maxpna(maxat,i)=maxpna(maxat,i)+1
        !$OMP ATOMIC
        minpna(minat,i)=minpna(minat,i)+1
      endif
 
    endif
 
    003 continue
 
  enddo ! end atoms-loop
  !$OMP END DO NOWAIT
 
  if (tbr .or. alc) goto 002
 
  !$OMP CRITICAL
  do k=3, taillr
    if (tring(k).gt.0) then
      if (nring(k,i).gt.0) then
        o = 0
        do l=1, tring(k)
          do m=1, nring(k,i)
            saut=.true.
            do n=1, k
              if (savring(k,m,n) .ne. savr(k,l,n)) then
                saut=.false.
                exit
              endif
            enddo
            if (saut) exit
          enddo
          if (.not.saut) then
            o = o + 1
            if (nring(k,i)+o .gt. numa) then
              tbr=.true.
              goto 004
            endif
            do m=1, k
              savring(k,nring(k,i)+o,m) = savr(k,l,m)
              ordring(k,nring(k,i)+o,m) = ordr(k,l,m)
            enddo
          endif
        enddo
        nring(k,i)=nring(k,i)+o
      else
        do l=1, tring(k)
          do m=1, k
            savring(k,l,m) = savr(k,l,m)
            ordring(k,l,m) = ordr(k,l,m)
          enddo
        enddo
        nring(k,i) = tring(k)
      endif
    endif
  enddo
 
  004 continue
  !$OMP END CRITICAL
 
  002 continue
 
  if (allocated(indte)) deallocate (indte)
  if (allocated(apna)) deallocate (apna)
  if (allocated(tring)) deallocate (tring)
  if (allocated(savr)) deallocate (savr)
  if (allocated(ordr)) deallocate (ordr)
  if (allocated(matdist)) deallocate(matdist)
  if (allocated(queue)) deallocate(queue)
  if (allocated(pringord)) deallocate(pringord)
 
  !$OMP END PARALLEL
 
  if (alc .or. tbr) goto 001
  ri = ri + rings_to_ogl(i, rid, nring, savring, ordring)
 
  do k=3, taillr
    if (nring(k,i) .ne. 0) then
!      do j=1, NRING(k,i)
!!  The algorithm implies that you may have a problem for
!!  the biggest size of ring therefore lets put this size out of the proof checking.
!!  Furthermore we check results only if all atoms are used to initiate the search.
!        if (mod(INDRING(k,j,i), k).ne.0 .and. k.lt.TAILLR .and. LTLT.eq.0) then
!          write (6, 003) i, k, j, INDRING(k,j,i)
!          write (6, 007)
!          do l=1, k
!            write (6, '(i4,2x)', advance='no') RINGORD(k,j,l,i)
!          enddo
!          write (6, *)
!          write (6, 004)
!          if (ABAB) then
!            write (6, 008)
!          else
!            write (6, 005)
!          endif
!          write (6, *)
!        endif
!      enddo
    endif
  enddo
 
enddo
 
001 continue
 
if (alc) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Subroutine: PRIMITIVE_RING_SEARCH_ATOMS"//char(0), "Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
endif
 
if (allocated(cpat)) deallocate (cpat)
if (allocated(vpat)) deallocate (vpat)
if (allocated(savring)) deallocate (savring)
if (allocated(ordring)) deallocate (ordring)
 
if (ri .eq. ns) ri = rings_to_ogl_menu(rid, nring)
 
END SUBROUTINE
#endif
 
#ifdef OPENMP
SUBROUTINE primitive_ring_search_steps (RID, NUMTH)
 
USE parameters
!$ USE OMP_LIB
IMPLICIT NONE
INTEGER, INTENT(IN) :: numth
#else
SUBROUTINE primitive_ring_search_steps (RID)
 
USE parameters
IMPLICIT NONE
#endif
INTEGER, INTENT(IN) :: RID
INTEGER, DIMENSION(:), ALLOCATABLE :: TRING, INDTE
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: SAVRING, ORDRING
INTEGER :: ri
LOGICAL, DIMENSION(2) :: FNDTAB
INTERFACE
  INTEGER FUNCTION rings_to_ogl_menu (IDSEARCH, NRI)
    USE parameters
    INTEGER, INTENT(IN) :: IDSEARCH
    INTEGER, DIMENSION(TAILLR, NS), INTENT(IN) :: NRI
  END FUNCTION
  INTEGER FUNCTION rings_to_ogl (STEP, IDSEARCH, NRI, RSAVED, OSAVED)
    USE parameters
    INTEGER, INTENT(IN) :: STEP, IDSEARCH
    INTEGER, DIMENSION(TAILLR,NS), INTENT(IN) :: NRI
    INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(IN) :: RSAVED, OSAVED
  END FUNCTION
END INTERFACE
 
ri = 0
#ifdef OPENMP
! OpenMP on steps only
!$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
!$OMP& PRIVATE(FNDTAB, MAXAT, MINAT, SAUT, PATH, PATHOUT, &
!$OMP& h, i, j, k, l, m, n, o, p, INDTE, APNA, RES_LIST, &
!$OMP& ERR, TRING, SAVRING, ORDRING, CPAT, VPAT, &
!$OMP& PRINGORD, NPRING, MATDIST, QUEUE, QUERNG) &
!$OMP& SHARED(NUMTH, RID, CALC_STRINGS, NS, NA, NNA, NNP, TLT, NSP, LOT, TAILLR, CONTJ, VOISJ, &
!$OMP& NUMA, MAXPNA, MINPNA, ABAB, NO_HOMO, TBR, ALC, ALC_TAB, &
!$OMP& NCELLS, THE_BOX, FULLPOS, PBC, MAXN, NRING, INDRING, PNA, ri)
#endif
if(allocated(matdist)) deallocate(matdist)
allocate(matdist(nna), stat=err)
if (err .ne. 0) then
  alc_tab="MATDIST"
  alc=.true.
  goto 001
endif
if(allocated(queue)) deallocate(queue)
allocate(queue(nna), stat=err)
if (err .ne. 0) then
  alc_tab="QUEUE"
  alc=.true.
  goto 001
endif
if (allocated(pringord)) deallocate(pringord)
allocate(pringord(numa*10,taillr), stat=err)
if (err .ne. 0) then
  alc_tab="PRINGORD"
  alc=.false.
  goto 001
endif
if (allocated(npring)) deallocate(npring)
allocate(npring(nna), stat=err)
if (err .ne. 0) then
  alc_tab="NPRING"
  alc=.false.
  goto 001
endif
if(allocated(cpat)) deallocate(cpat)
allocate(cpat(nna), stat=err)
if (err .ne. 0) then
  alc_tab="CPAT"
  alc=.true.
  goto 001
endif
if(allocated(vpat)) deallocate(vpat)
allocate(vpat(nna,maxn), stat=err)
if (err .ne. 0) then
  alc_tab="VPAT"
  alc=.true.
  goto 001
endif
if(allocated(res_list)) deallocate(res_list)
allocate(res_list(taillr), stat=err)
if (err .ne. 0) then
  alc_tab="RES_LIST"
  alc=.true.
  goto 001
endif
if(allocated(savring)) deallocate(savring)
allocate(savring(taillr,numa,taillr), stat=err)
if (err .ne. 0) then
  alc_tab="SAVRING"
  alc=.true.
  goto 001
endif
if(allocated(ordring)) deallocate(ordring)
allocate(ordring(taillr,numa,taillr), stat=err)
if (err .ne. 0) then
  alc_tab="ORDRING"
  alc=.true.
  goto 001
endif
if(allocated(indte)) deallocate(indte)
allocate(indte(numa), stat=err)
if (err .ne. 0) then
  alc_tab="INDTE"
  goto 001
endif
if (allocated(tring)) deallocate(tring)
allocate(tring(taillr), stat=err)
if (err .ne. 0) then
  alc_tab="TRING"
  alc=.true.
  goto 001
endif
if(allocated(apna)) deallocate(apna)
allocate(apna(taillr), stat=err)
if (err .ne. 0) then
  alc_tab="APNA"
  alc=.true.
  goto 001
endif
 
!$OMP DO SCHEDULE(STATIC,NS/NUMTH)
do i=1, ns
 
  if (tbr .or. alc) goto 002
  savring(:,:,:)=0
  ordring(:,:,:)=0
  tring(:)=0
  call setup_cpat_vpat_ring (nna, i, contj, voisj, cpat, vpat)
 
  do j=nnp+1, nnp+na ! atoms-loop
 
    if (tlt .eq. nsp+1 .or. lot(j-nnp) .eq. tlt) then
 
      apna(:)=0
      minat=taillr
      maxat=1
      saut=.true.
 
      call dijkstra(j, cpat, vpat, queue, matdist)
      if (tbr .or. alc) goto 002
 
      do k=1, taillr/2 + mod(taillr,2) ! ring-sizes-loop
 
        indte(:)=0
        res_list(:)=0
        path=0
        do l=1, nna
          if (matdist(l) .eq. k) then
            call spath_rec (path,l,k,k,matdist,cpat,vpat,npring,pringord)
          endif
        enddo
        h = path*(path-1)/2
        if (allocated(querng)) deallocate(querng)
        allocate(querng(h,2), stat=err)
        if (err .ne. 0) then
          alc_tab="QUERNG"
          alc=.true.
          goto 002
        endif
 
        l=0
        do m=1, path-1
          do n=m+1, path
            l=l+1
            querng(l,1)=m
            querng(l,2)=n
!           Paths which share atoms ... so with overlap ... are deleted
            pathout=.false.
            do o=1, k-1
              do p=1, k-1
                if (pringord(m,o) .eq. pringord(n,p)) then
                  pathout=.true.
                  exit
                endif
              enddo
              if (pathout) exit
            enddo
            if (pathout) then
              querng(l,1)=0
              querng(l,2)=0
              l=l-1
            endif
          enddo
        enddo
        fndtab(:)=.false.
        call prim_ring (fndtab, j, l, k, h, cpat, vpat, querng, pringord, matdist, &
                        savring, ordring, tring, indte, res_list)
        if (tbr .or. alc) goto 002
        m = 2*k
        if (fndtab(1)) then
          if (apna(m).eq.0) then
            apna(m)=1
            minat=min(minat,m)
            maxat=max(maxat,m)
            saut=.false.
          endif
        endif
        if (fndtab(2)) then
          m = m + 1
          if (apna(m).eq.0) then
            apna(m)=1
            minat=min(minat,m)
            maxat=max(maxat,m)
            saut=.false.
          endif
        endif
        if (allocated(querng)) deallocate(querng)
 
      enddo ! end ring-sizes-loop
 
      if (.not. saut) then
        do k=3, taillr
          do m=3, taillr
            if (apna(k).eq.1 .and. apna(m).eq.1) then
              pna(k,m,i)=pna(k,m,i)+1
            endif
          enddo
        enddo
        maxpna(maxat,i)=maxpna(maxat,i)+1
        minpna(minat,i)=minpna(minat,i)+1
      endif
 
    endif
 
  enddo ! end atoms-loop
 
  do j=3, taillr
    nring(j,i) = tring(j)
  enddo
 
  k = rings_to_ogl(i, rid, nring, savring, ordring)
  !$OMP ATOMIC
  ri = ri + k
 
  do k=3, taillr
    if (nring(k,i) .ne. 0) then
!      do j=1, NRING(k,i)
!!  The algorithm implies that you may have a problem for
!!  the biggest size of ring therefore lets put this size out of the proof checking.
!!  Furthermore we check results only if all atoms are used to initiate the search.
!        if (mod(INDRING(k,j,i), k).ne.0 .and. k.lt.TAILLR .and. LTLT.eq.0) then
!          write (6, 003) i, k, j, INDRING(k,j,i)
!          write (6, 007)
!          do l=1, k
!            write (6, '(i4,2x)', advance='no') RINGORD(k,j,l,i)
!          enddo
!          write (6, *)
!          write (6, 004)
!          if (ABAB) then
!            write (6, 008)
!          else
!            write (6, 005)
!          endif
!          write (6, *)
!        endif
!      enddo
    endif
  enddo
 
  002 continue
 
enddo
#ifdef OPENMP
!$OMP END DO NOWAIT
#endif
 
001 continue
 
if (alc) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Subroutine: PRIMITIVE_RING_SEARCH_STEPS"//char(0), "Table: "//alc_tab(1:len_trim(alc_tab))//char(0))
endif
 
if (allocated(tring)) deallocate (tring)
if (allocated(savring)) deallocate (savring)
if (allocated(ordring)) deallocate (ordring)
if (allocated(cpat)) deallocate (cpat)
if (allocated(vpat)) deallocate (vpat)
if (allocated(indte)) deallocate (indte)
if (allocated(apna)) deallocate (apna)
if(allocated(matdist)) deallocate(matdist)
if(allocated(queue)) deallocate(queue)
if (allocated(pringord)) deallocate(pringord)
 
#ifdef OPENMP
!$OMP END PARALLEL
#endif
 
if (ri .eq. ns) ri = rings_to_ogl_menu(rid, nring)
 
END SUBROUTINE
 
SUBROUTINE dijkstra(NODE, CPT, VPT, QUE, MATDIS)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: NODE
INTEGER, DIMENSION(NNA), INTENT(INOUT) :: QUE, MATDIS
INTEGER, DIMENSION(NNA), INTENT(IN) :: CPT
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: VPT
INTEGER :: QBEGIN, QEND, QID, AT1, AT2, DAT1
 
que(:)=0
matdis(:)=nna+2
que(1)=node
matdis(node)=0
qbegin=0
qend=1
 
! Trouver les plus cours chemin reliant un atome aux autres
! Find the shortest paths between atoms
 
do while (qbegin < qend)
 
  qbegin=qbegin+1
  at1=que(qbegin)
  dat1=matdis(at1)+1
 
  do qid=1, cpt(at1)
 
    at2=vpt(at1,qid)
    if (matdis(at2) .gt. dat1) then
 
      matdis(at2) = dat1
      if (dat1 < nna) then
 
        qend=qend+1
        que(qend)= at2
 
      endif
 
    endif
 
  enddo
 
enddo
 
END SUBROUTINE
 
RECURSIVE SUBROUTINE spath_rec (PTH, NODE, LENGTH, LNGTH, MATDIS, CPT, VPT, NPRI, PORDR)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(INOUT) :: pth
INTEGER, INTENT(IN) :: node, length, lngth
INTEGER, DIMENSION(NNA), INTENT(IN) :: cpt, matdis
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: vpt
INTEGER, DIMENSION(NNA), INTENT(INOUT) :: npri
INTEGER, DIMENSION(NUMA*10,TAILLR), INTENT(INOUT) :: pordr
INTEGER :: distnn, idv, vdi, idt
 
npri(length) = node
do idv=1, cpt(node)
 
! Boucle sur tous les voisins du noeud "NODE"
! Loop on all the neighbors of node "NODE"
 
  vdi=vpt(node,idv)                  ! Selection du voisin - neighbor selection
  distnn=matdis(vdi)                 ! Distance to starting atom = VDI in the node distance table
 
  if (distnn .eq. 0) then
 
    pth=pth+1
    do idt=1, lngth
      pordr(pth,idt)=npri(idt)
    enddo
 
  elseif (distnn .eq. length-1) then
 
    call spath_rec (pth, vdi, distnn, lngth, matdis, cpt, vpt, npri, pordr)
 
  endif
 
enddo
 
END SUBROUTINE
 
INTEGER FUNCTION real_atom_id (IND, NATS)
 
IMPLICIT NONE
INTEGER, INTENT(IN) :: ind, nats
 
real_atom_id=  ind - (ind/nats)*nats
if (real_atom_id .eq. 0) real_atom_id=nats
 
END FUNCTION
 
SUBROUTINE prim_ring (FNDTAB, NODE, PTH, LGTH, NPT, CPT, VPT, QRNG, PORD, MATDIS, &
                      RSAVED, OSAVED, TRIN, INDP, RESLP)
 
USE parameters
 
IMPLICIT NONE
 
LOGICAL, DIMENSION(2), INTENT(INOUT) :: FNDTAB
INTEGER, INTENT(IN) :: NODE, PTH, LGTH, NPT
INTEGER, DIMENSION(NNA), INTENT(IN) :: CPT, MATDIS
INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: TRIN
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: VPT
INTEGER, DIMENSION(NPT,2), INTENT(IN) :: QRNG
INTEGER, DIMENSION(NUMA*10,TAILLR), INTENT(INOUT) :: PORD
INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDP
INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESLP
INTEGER :: RM, RN, PR, PTH1, PTH2, IRS, IRX
INTEGER :: ATB, ATA, ATC, ATD, ATE, MAXD, MIND
INTEGER :: PROBE
INTEGER, DIMENSION(TAILLR) :: TOPRIM, PRIMTO
LOGICAL:: GOAL, TOSAVE
LOGICAL, DIMENSION(NNA) :: CHK
INTERFACE
  INTEGER FUNCTION real_atom_id (IND, NATS)
    INTEGER, INTENT(IN) :: IND, NATS
  END FUNCTION
END INTERFACE
 
goal=.false.
 
do pr=1, pth
 
! Premier chemin i->j - First path i->j
  pth1= qrng(pr,1)
! Second chemin i->k - second path i->k
  pth2= qrng(pr,2)
 
  probe=0
  ata= pord(pth1,lgth)
  atb= pord(pth2,lgth)
  do rm=1, cpt(ata)
    ate=vpt(ata,rm)
    if (ate .eq. atb) probe=1
  enddo
 
  if (ata .eq. atb .or. probe .eq. 1) then
 
    do rm=1, lgth-1
    do rn=rm, lgth-1+probe
 
      atc= pord(pth1,rm)
      atd= pord(pth2,rn)
      maxd=matdis(atc)+matdis(atd)
      mind=2*lgth+probe-maxd
      chk(:)=.false.
 
      call pair_search (goal, atc, atd, 1, maxd, mind, cpt, vpt, chk)
      if (goal) then
        goal= .false.
        goto 001
      endif
 
    enddo
    enddo
 
    do rn=1, lgth-1
    do rm=rn, lgth-1+probe
 
      atc= pord(pth1,rm)
      atd= pord(pth2,rn)
      maxd=matdis(atc)+matdis(atd)
      mind=2*lgth+probe-maxd
      chk(:)=.false.
 
      call pair_search (goal, atc, atd, 1, maxd, mind, cpt, vpt, chk)
      if (goal) then
        goal= .false.
        goto 001
      endif
 
    enddo
    enddo
 
    if ((probe.eq.0 .and. 2*lgth.le.taillr) .or. (probe.eq.1 .and. 2*lgth.lt.taillr)) then
 
      toprim(:)=0
      primto(:)=0
      do irx=1, lgth
        toprim(irx)=pord(pth1,irx)
      enddo
      if (probe .eq. 1) toprim(lgth+1)=pord(pth2,lgth)
      irs=lgth+probe
      do irx=lgth-1, 1, -1
        irs=irs+1
        toprim(irs)=pord(pth2,irx)
      enddo
      toprim(irs+1)=node
 
!      To find atom id in the primitive unit cell
      do irx=1, 2*lgth+probe
        primto(irx)=real_atom_id(toprim(irx),na)
      enddo
      tosave=.true.
      if (abab) then
        if (probe .eq. 0) then
          call testabab (tosave,lgth,primto)
        else
          tosave=.false.
        endif
      endif
      if (no_homo) then
        call testhomo(tosave,lgth,primto)
      endif
      if (tosave .and. lgth*2+probe.le.taillr) then
 
        if (calc_strings) then
           call strong_rings (fndtab, lgth, probe, toprim, primto, rsaved, osaved, trin, indp, reslp, cpt, vpt)
        else
          call save_dijkstra_ring (primto, lgth*2+probe, rsaved, osaved, trin, indp, reslp)
          fndtab(1+probe)=.true.
        endif
        if (tbr .or. alc) goto 002
 
      endif
 
    endif
 
  endif
 
001 continue
 
enddo
 
002 continue
 
END
 
RECURSIVE SUBROUTINE pair_search (GOAL, AT1, AT2, LG, MAXM, MINM, CPT, VPT, CHK)
 
USE parameters
 
IMPLICIT NONE
 
LOGICAL, INTENT(INOUT) :: goal
INTEGER, INTENT(IN) :: at1, at2, lg, maxm, minm
INTEGER, DIMENSION(NNA), INTENT(IN) :: cpt
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: vpt
LOGICAL, DIMENSION(NNA), INTENT(INOUT) :: chk
INTEGER :: psc, at3
 
chk(at1)=.true.
 
if (at1 .eq. at2) then
 
  goal=.true.
 
else
 
  do psc=1, cpt(at1)
 
    at3=vpt(at1,psc)
 
    if (.not.chk(at3)) then
 
      if (at3.eq.at2) then
 
        goal=.true.
        goto 001
 
      elseif (lg.lt.maxm-1 .and. lg.lt.minm-1) then
 
        call pair_search (goal, at3, at2, lg+1, maxm, minm, cpt, vpt, chk)
 
        if (goal) then
           goto 001
        endif
 
      endif
 
    endif
 
  enddo
 
endif
 
001 continue
 
chk(at1)=.false.
 
END SUBROUTINE
 
SUBROUTINE strong_rings (FNDTAB, RLGTH, RPROBE, TOPRIM, PRIMTO, ASRING, OSRING, TRNG, INDT, RESL, CPT, VPT)
 
USE parameters
 
IMPLICIT NONE
 
LOGICAL, DIMENSION(2), INTENT(INOUT) :: FNDTAB
INTEGER, INTENT(IN) :: RLGTH, RPROBE
INTEGER, DIMENSION(TAILLR), INTENT(IN) :: TOPRIM, PRIMTO
INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: ASRING, OSRING
INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDT
INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: TRNG, RESL
INTEGER, DIMENSION(NNA), INTENT(IN) :: CPT
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: VPT
INTEGER :: RM, RN, STLGT
LOGICAL, DIMENSION(NNA) :: CHKS
 
LOGICAL :: DVSTR
LOGICAL :: SGOAL=.false.
LOGICAL :: FGOAL=.false.
 
stlgt=2*rlgth+rprobe
 
chks(:)=.false.
 
do rm=1, stlgt-1
 
  do rn=rm+1, stlgt
 
    dvstr=.false.
    call search_strong_rings (dvstr, rm, rn, stlgt, stlgt, sgoal, fgoal, toprim, cpt, vpt, chks)
    if (tbr .or. alc) goto 002
    if (fgoal) goto 001
    if (sgoal) goto 002
 
  enddo
 
enddo
 
001 continue
 
call save_dijkstra_ring (primto, rlgth*2+rprobe, asring, osring, trng, indt, resl)
fndtab(1+rprobe)=.true.
 
if (tbr .or. alc) goto 002
 
002 continue
 
END SUBROUTINE
 
RECURSIVE SUBROUTINE search_strong_rings (DLOW, IDX, IDY, DLGTR, LGTR, SRGOAL, FRGOAL, TOTER, CPT, VPT, CHKS)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: dlgtr, lgtr, idx, idy
LOGICAL, INTENT(INOUT) :: dlow, srgoal, frgoal
INTEGER, DIMENSION(DLGTR), INTENT(IN) :: toter
INTEGER, DIMENSION(NNA), INTENT(IN) :: cpt
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: vpt
LOGICAL, DIMENSION(NNA), INTENT(INOUT) :: chks
INTEGER :: dxy, dyx, dtest, dmin
INTEGER :: rg, rf, rh, rj, rl, ro, rp, rq, rr
INTEGER :: nbpath,  idz
INTEGER, DIMENSION(:), ALLOCATABLE :: strst, newter
INTEGER, DIMENSION(:,:), ALLOCATABLE :: totpath
LOGICAL, DIMENSION(:), ALLOCATABLE :: cutring
LOGICAL, DIMENSION(NNA) :: chk
TYPE(ring), TARGET :: st_ring
TYPE(ring), POINTER :: sring, tmp_ri
 
INTERFACE
  RECURSIVE SUBROUTINE path_search (AT0,AT1,AT2,IDT1,IDT2,DMAX,DMED,DARING,NPATH,DLPATH, &
                                    STPATH,ACRING,TPATH,CUTPATH,CPT,VPT,CHK,CHKS,SRING)
    USE parameters
    INTEGER, INTENT(IN) :: at0, at1, at2, idt1, idt2, dmax, dmed, daring
    INTEGER, INTENT(INOUT) :: npath
    LOGICAL, INTENT(INOUT) :: dlpath
    INTEGER, DIMENSION(DMAX*DMAX), INTENT(INOUT) :: stpath
    INTEGER, DIMENSION(DMAX*DMAX,DMAX), INTENT(INOUT) :: tpath
    INTEGER, DIMENSION(DMED), INTENT(IN) :: acring
    INTEGER, DIMENSION(NNA), INTENT(IN) :: cpt
    INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: vpt
    LOGICAL, DIMENSION(DMAX*DMAX), INTENT(INOUT) :: cutpath
    LOGICAL, DIMENSION(NNA), INTENT(INOUT) :: chk, chks
    TYPE (ring), POINTER, INTENT(INOUT) :: sring
  END
  SUBROUTINE shortcut_ring (IDA, IDB, IDC, DLT, TABAB, SRING)
    USE parameters
    INTEGER, INTENT(IN) :: ida, idb, idc, dlt
    INTEGER, DIMENSION(DLT), INTENT(IN) :: tabab
    TYPE(ring), POINTER, INTENT(INOUT) :: sring
  END
  SUBROUTINE creat_ring (RING_INIT, ELEM_CR)
    USE parameters
    TYPE (ring), INTENT(INOUT) :: ring_init
    INTEGER, INTENT(IN) :: elem_cr
  END SUBROUTINE
  SUBROUTINE do_ring (THE_RING, ELEM_DO)
    USE parameters
    TYPE (ring), POINTER, INTENT(INOUT) :: the_ring
    INTEGER, INTENT(IN) :: elem_do
  END SUBROUTINE
END INTERFACE
 
do rg=1, dlgtr
  chks(toter(rg))=.true.
enddo
 
dxy=abs(idy-idx)
dyx=dlgtr-dxy
dmin=min(dxy,dyx)
dtest=lgtr-1-dmin
 
!   First we find find a ring starting from an atom i and enclosing
!   all tested atoms, then this ring will be tested and so on
!   if the size of the enclosing ring became smaller than
!   DLGTR then the initial ring is not strong
 
if (dtest.ge.2) then
 
  do rj=2, dtest
 
    if (allocated(totpath)) deallocate(totpath)
    allocate(totpath(rj*rj,rj), stat=err)
    if (err .ne. 0) then
      alc_tab="TOTPATH"
      alc=.true.
      goto 001
    endif
    if (allocated(cutring)) deallocate(cutring)
    allocate(cutring(rj*rj), stat=err)
    if (err .ne. 0) then
      alc_tab="CUTRING"
      alc=.true.
      goto 001
    endif
    if (allocated(strst)) deallocate(strst)
    allocate(strst(rj*rj), stat=err)
    if (err .ne. 0) then
      alc_tab="STRST"
      alc=.true.
      goto 001
    endif
 
    nbpath=0
    strst(:)=0
    cutring(:)=.false.
    totpath(:,:)=0
    chk(:)=.false.
 
    rr=1
    call creat_ring (st_ring, toter(idx))
    sring => st_ring
! On recherche tous les chemins de taille RJ entre les atomes TOTER(IDX) et TOTER(IDY)
! On stock ces NBPATH chemins dans le tableau TOTPATH
 
! We are looking for all path of size RJ between atoms TOTER(IDX) and TOTER(IDY)
! Then these NBPATH are saved in the tab TOTPATH
 
    call path_search (toter(idx),toter(idx),toter(idy),idx,idy,rj,dlgtr,lgtr,nbpath,dlow, &
                      strst,toter,totpath,cutring,cpt,vpt,chk,chks,sring)
    if (alc) goto 001
    if (nbpath .ne. 0) then
 
      do rh=1, nbpath
 
        idz=0
        if (cutring(rh)) then
          do rg=1, dlgtr
            if (toter(rg) .eq. strst(rh)) then
              idz=rg
              exit
            endif
          enddo
        endif
 
        if(idz.eq.0)then
 
          if (.not.cutring(rh)) then
            call creat_ring (st_ring, toter(idx))
            sring => st_ring
            do rg=1, rj
              rr=rr+1
              call do_ring (sring, totpath(rh,rg))
              if (alc) goto 001
            enddo
            if (dxy .le. dyx) then
              do rg=idy+1, dlgtr
                call do_ring (sring, toter(rg))
                if (alc) goto 001
              enddo
              do rg=1, idx
                call do_ring (sring, toter(rg))
                if (alc) goto 001
              enddo
            else
              do rg=idy-1, idx, -1
                call do_ring (sring, toter(rg))
                if (alc) goto 001
              enddo
            endif
          else
            call creat_ring (st_ring, toter(1))
            sring => st_ring
            do rg=2, dlgtr
              call do_ring (sring, toter(rg))
              if (alc) goto 001
            enddo
            chks(totpath(rh,1))=.true.
          endif
 
        else
 
          call creat_ring (st_ring, totpath(rh,1))
          sring => st_ring
          call shortcut_ring (idx, idy, idz, dlgtr, toter, sring)
          if (alc) goto 001
 
        endif
 
!    A new path as been found, so a new rings as been built according
!    to the new list of atoms.
!    The size of this new ring has to be tested
 
        if (sring%SPEC+1 .lt. lgtr) then
 
          tmp_ri => sring
          do rf=1, sring%SPEC
            tmp_ri => tmp_ri%PAST
            deallocate (tmp_ri%NEXT)
          enddo
          srgoal=.true.
          goto 001
 
        else
 
! If all the atoms are CHK and ring is bigger than the initial ring
! the initial ring is strong .. we have to check this
          rf=0
          do rg=1, nna
            if (chks(rg)) rf=rf+1
          enddo
          if (rf .eq. nna) then
 
            tmp_ri => sring
            do rg=1, sring%SPEC
              tmp_ri => tmp_ri%PAST
              deallocate (tmp_ri%NEXT)
            enddo
            frgoal=.true.
            goto 001
 
          else
 
!       Else the new ring has to be tested
            rr = sring%SPEC+1
            if (allocated(newter)) deallocate(newter)
            allocate(newter(rr), stat=err)
            if (err .ne. 0) then
              alc_tab="NEWTER"
              alc=.true.
              goto 001
            endif
            srgoal=.false.
            tmp_ri => sring
            do rg=1, rr
              newter(rr-rg+1)=tmp_ri%ATOM
              if (rg .lt. rr) then
                tmp_ri => tmp_ri%PAST
                deallocate (tmp_ri%NEXT)
              endif
            enddo
            dlow=.false.
            sring => st_ring
            if (idz.eq.0) then
              if (dmin .eq. 1) then
                ro=2
                rp=rj
                dlow=.true.
              else
                ro=1
                rp=rr-1
              endif
            else
              ro=1
              rp=rr-1
            endif
 
            do rg=ro, rp
 
              if (idz.eq.0) then
                if (dmin .eq. 1) then
                  dlow=.true.
                  rq=rj+2
                else
                  rq=rg+1
                endif
              else
                 rq=rg+1
              endif
 
              do rf=rq, rr
                call search_strong_rings (dlow, rg, rf, rr, lgtr, srgoal, frgoal, newter, cpt, vpt, chks)
                if (tbr .or. alc) goto 001
                if (srgoal .or. frgoal) goto 001
                do rl=1, dlgtr
                  chks(toter(rl))=.true.
                enddo
                if (idz.eq.0 .and. cutring(rh)) chks(totpath(rh,1))=.true.
              enddo
 
            enddo
 
            if (allocated(newter)) deallocate(newter)
 
          endif
 
        endif
 
        if (idz.eq.0 .and. cutring(rh)) chks(totpath(rh,1))=.false.
 
      enddo
 
    else
 
      srgoal=.false.
 
    endif
 
    if (allocated(totpath)) deallocate(totpath)
    if (allocated(cutring)) deallocate(cutring)
    if (allocated(strst)) deallocate(strst)
 
  enddo
 
endif
 
001 continue
 
if (allocated(totpath)) deallocate(totpath)
if (allocated(cutring)) deallocate(cutring)
if (allocated(strst)) deallocate(strst)
if (allocated(newter)) deallocate(newter)
 
do rg=1, dlgtr
  chks(toter(rg))=.false.
enddo
 
END SUBROUTINE
 
RECURSIVE SUBROUTINE path_search (AT0,AT1,AT2,IDT1,IDT2,DMAX,DMED,DARING,NPATH,DLPATH, &
                                  STPATH,ACRING,TPATH,CUTPATH,CPT,VPT,CHK,CHKS,SRING)
USE parameters
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: at0, at1, at2, idt1, idt2, dmax, dmed, daring
INTEGER, INTENT(INOUT) :: npath
LOGICAL, INTENT(INOUT) :: dlpath
INTEGER, DIMENSION(DMAX*DMAX), INTENT(INOUT) :: stpath
INTEGER, DIMENSION(DMAX*DMAX,DMAX), INTENT(INOUT) :: tpath
INTEGER, DIMENSION(DMED), INTENT(IN) :: acring
INTEGER, DIMENSION(NNA), INTENT(IN) :: cpt
INTEGER, DIMENSION(NNA,MAXN), INTENT(IN) :: vpt
LOGICAL, DIMENSION(DMAX*DMAX), INTENT(INOUT) :: cutpath
LOGICAL, DIMENSION(NNA), INTENT(INOUT) :: chk, chks
TYPE (ring), POINTER, INTENT(INOUT) :: sring
 
INTEGER :: at3, at4, at5, at6, at7
INTEGER :: duv, duw, dvw
LOGICAL :: val1, val2
LOGICAL :: touch
TYPE (ring), POINTER :: tmpr
 
INTERFACE
  SUBROUTINE do_ring (THE_RING, ELEM_DO)
    USE parameters
    TYPE (ring), POINTER, INTENT(INOUT) :: the_ring
    INTEGER, INTENT(IN) :: elem_do
  END SUBROUTINE
END INTERFACE
 
chk(at1)=.true.
 
if (at1.eq.at2 .and. sring%SPEC.eq.dmax) then
 
 
  npath=npath+1
  tmpr => sring
  do at3=1, dmax
    tpath(npath,dmax-at3+1)=tmpr%ATOM
    if (at3 .lt. dmax) tmpr => tmpr%PAST
  enddo
 
else
 
  do at4=1, cpt(at1)
 
    at3=vpt(at1,at4)
    if (.not.chk(at3)) then
 
      call do_ring (sring, at3)
      if (alc) goto 001
      if (at3 .eq. at2 .and. sring%SPEC.eq.dmax) then
 
        npath=npath+1
        tmpr => sring
        do at5=1, dmax
          tpath(npath,dmax-at5+1)=tmpr%ATOM
          if (at5 .lt. dmax) tmpr => tmpr%PAST
        enddo
 
      elseif (sring%SPEC.lt.dmax .and. .not.chks(at3)) then
 
        touch=.false.
        at7=0
        val1=.false.
        val2=.false.
        do at5=1, cpt(at3)
          at6=vpt(at3,at5)
          if (at6 .eq. at0) val1=.true.
          if (at6 .eq. at2) val2=.true.
          if (chks(at6)) then
            at7=at7+1
            if (at6.ne.at0 .and. at6.ne.at2) stpath(npath+1)=at6
          endif
        enddo
        at6=0
        do at5=1, dmed
          if (acring(at5) .eq. stpath(npath+1)) then
            at6=at5
            exit
          endif
        enddo
        if (at7 .eq. 0) then
          touch=.false.
        elseif (at7.eq.1 .and. val1 .and. .not.val2) then
          touch=.false.
        elseif (at7.eq.1 .and. val2 .and. .not.val1) then
          touch=.false.
        elseif (at7.eq.1 .and. .not.val1 .and. .not.val2) then
          touch=.true.
        elseif (at7.eq.2 .and. val1 .and. val2 .and. dmax.eq.2) then
          touch=.false.
        elseif (at7.gt.2 .and. val1 .and. val2 .and. dmax.eq.2) then
          if (dlpath) then
            touch=.true.
          else
            if (at6 .ne. 0) then
              if (at6 .lt. idt1) then
                duv=idt1-at6
                dvw=idt2-idt1
                duw=dmed-(idt2-at6)
              elseif(at6 .lt. idt2) then
                duv=at6-idt1
                dvw=dmed-(at2-at1)
                duw=idt2-at6
              elseif(at6 .gt. idt2) then
                duv=dmed-(at6-idt1)
                dvw=idt2-idt1
                duw=at6-idt2
              endif
              if (duv.lt.daring-1 .and. duw.lt.daring-1) then
                cutpath(npath+1)=.true.
                touch=.false.
              elseif (dvw.lt.daring-1 .and. duw.lt.daring-1) then
                cutpath(npath+1)=.true.
                touch=.false.
              elseif (duv.lt.daring-1 .and. dvw.lt.daring-1) then
                cutpath(npath+1)=.true.
                touch=.false.
              else
                touch=.true.
              endif
            else
               cutpath(npath+1)=.true.
               touch=.false.
            endif
          endif
        else
          touch=.true.
        endif
        if (.not.touch) then
          call path_search (at0,at3,at2,idt1,idt2,dmax,dmed,daring,npath,dlpath, &
                            stpath,acring,tpath,cutpath,cpt,vpt,chk,chks,sring)
          if (alc) goto 001
        endif
      endif
      sring => sring%PAST
      deallocate (sring%NEXT)
 
    endif
 
  enddo
 
endif
 
chk(at1)=.false.
001 continue
 
END SUBROUTINE
 
SUBROUTINE shortcut_ring (IDA, IDB, IDC, DLT, TABAB, SRING)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, INTENT(IN) :: IDA, IDB, IDC, DLT
INTEGER, DIMENSION(DLT), INTENT(IN) :: TABAB
TYPE(ring), POINTER, INTENT(INOUT) :: SRING
INTEGER :: IDD
INTEGER :: LXY, LXZ, LYZ
 
INTERFACE
  SUBROUTINE do_ring (THE_RING, ELEM_DO)
    USE parameters
    TYPE (RING), POINTER, INTENT(INOUT) :: THE_RING
    INTEGER, INTENT(IN) :: ELEM_DO
  END SUBROUTINE
END INTERFACE
 
if (ida .gt. idc) then
! IDB > IDA > IDC
  lxz=ida-idc
  lxy=idb-ida
  lyz=dlt-(idb-idc)
  if (lxy .gt. lyz) then
    if (lxy .gt. lxz) then
      do idd=ida, idb
        call do_ring (sring, tabab(idd))
        if (alc) goto 001
      enddo
    else
      do idd=idc, ida
        call do_ring (sring, tabab(idd))
        if (alc) goto 001
      enddo
    endif
  else
    if (lxz .gt. lyz) then
      do idd=idc, ida
        call do_ring (sring, tabab(idd))
        if (alc) goto 001
      enddo
    else
      do idd=idb, dlt
        call do_ring (sring, tabab(idd))
        if (alc) goto 001
      enddo
      do idd=1, idc
        call do_ring (sring, tabab(idd))
        if (alc) goto 001
      enddo
    endif
  endif
else
! IDB > IDA .and. IDC .gt. IDA
  if (idc .gt. idb) then
! IDB > IDA .and. IDC .gt. IDA .and. IDC .gt. IDB
! =>  IDC > IDB > IDA
    lxy=idb-ida
    lyz=idc-idb
    lxz=dlt-(idc-ida)
    if (lxy .gt. lyz) then
      if (lxy .gt. lxz) then
        do idd=ida, idb
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      else
        do idd=idc, dlt
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
        do idd=1, ida
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      endif
    else
      if (lyz .ge. lxz) then
        do idd=idb, idc
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      else
        do idd=ida, idb
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      endif
    endif
  else
! IDB > IDA .and. IDC .gt. IDA .and. IDB .gt. IDC
! =>  IDB > IDC > IDA
    lxy=dlt-(idb-ida)
    lyz=idb-idc
    lxz=idc-ida
    if (lxy .gt. lyz) then
      if (lxy .gt. lxz) then
        do idd=idb, dlt
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
        do idd=1, ida
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      else
        do idd=ida, idc
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      endif
    else
      if (lyz .gt. lxz) then
        do idd=idc, idb
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      else
        do idd=ida, idc
          call do_ring (sring, tabab(idd))
          if (alc) goto 001
        enddo
      endif
    endif
  endif
endif
 
001 continue
 
END SUBROUTINE
 
SUBROUTINE testabab(VALTEST, LGTEST, PRIMT)
 
USE parameters
 
IMPLICIT NONE
 
LOGICAL, INTENT(INOUT) :: VALTEST
INTEGER, INTENT(IN) :: LGTEST
INTEGER, DIMENSION(TAILLR), INTENT(IN) :: PRIMT
INTEGER :: LTA, LTB, LTC, LTD
 
lta=lot(primt(1))
ltb=lot(primt(2))
 
valtest=.true.
if (lta .ne. ltb) then
  do ltd=3, 2*lgtest
    ltc=lot(primt(ltd))
    if (mod(ltd,2).eq.0 .and. ltc.ne.ltb) then
      valtest=.false.
      exit
    elseif (mod(ltd,2).ne.0 .and. ltc.ne.lta) then
      valtest=.false.
      exit
    endif
  enddo
else
  valtest=.false.
endif
 
END SUBROUTINE
 
SUBROUTINE testhomo(VALTEST, LGTEST, PRIMT)
 
USE parameters
 
IMPLICIT NONE
 
LOGICAL, INTENT(INOUT) :: VALTEST
INTEGER, INTENT(IN) :: LGTEST
INTEGER, DIMENSION(TAILLR), INTENT(IN) :: PRIMT
INTEGER :: LTA, LTB
 
valtest=.true.
do lta=1, 2*lgtest
  ltb=lta+1
  if (lta .eq. 2*lgtest) ltb=1
  if (lot(primt(lta)) .eq. lot(primt(ltb))) then
    valtest=.false.
    exit
  endif
enddo
 
END SUBROUTINE
 
SUBROUTINE save_dijkstra_ring (TAB, TLES, RSAVED, OSAVED, TRING, INDT, RESL)
 
USE parameters
 
IMPLICIT NONE
 
INTEGER, DIMENSION(TAILLR), INTENT(IN) :: TAB
INTEGER, INTENT(IN) :: TLES
INTEGER, DIMENSION(TAILLR,NUMA,TAILLR), INTENT(INOUT) :: RSAVED, OSAVED
INTEGER, DIMENSION(NUMA), INTENT(INOUT) :: INDT
INTEGER, DIMENSION(TAILLR), INTENT(INOUT) :: RESL, TRING
INTEGER :: idx, idy
LOGICAL :: NEWRING
INTEGER, DIMENSION(TLES) :: TOTRI, TOSAV
 
! A ring has been found, we need to check if it has already been found or not
 
do idx=1, tles
 
! Creation of two tab first for the numerical sorting of the list
! The other without sorting for output
  totri(idx)= tab(idx)
  tosav(idx)= tab(idx)
 
enddo
 
call tri(totri, tles)
 
if (tring(tles) .ne. 0) then
 
  do idx=1, tring(tles)
!   do-loop on existing rings to check if the new on has already been found
 
    newring=.false.
 
    do idy=1, tles
      if (totri(idy) .ne. rsaved(tles,idx,idy)) then
        newring=.true.
        exit
      endif
    enddo
 
    if (.not.newring) then
 
      ! Already been found n-times, increment of the counter
      indt(idx)=indt(idx)+1
      exit
 
    endif
 
  enddo
 
else
 
  newring=.true.
 
endif
 
if (newring) then
 
! A new ring has been found
! We save the data
  resl(tles)=resl(tles)+1
  tring(tles)=tring(tles)+1
  if (tring(tles) .gt. numa) then
    tbr=.true.
    goto 001
  endif
  indt(tring(tles))=indt(tring(tles))+1
  do idx=1, tles
    rsaved(tles,tring(tles),idx)=totri(idx)
    osaved(tles,tring(tles),idx)=tosav(idx)
  enddo
 
endif
 
001 continue
 
END SUBROUTINE