!-----------------------------------------------------------------!
! PURPOSE: Find unique artificial breakpoints for cyclic molecule !
!          including multi-ring compounds.                        !
! CREATED: June/July 2005, Julia Lee-TAYLOR, NCAR                 !
!-----------------------------------------------------------------!
      SUBROUTINE uniqring(nring,nca,group,bond,rank,rjg)
      IMPLICIT NONE
      include 'general.h'

! input: 
      INTEGER    nca,nring
      CHARACTER(LEN=lgr) group(mca)
! i/o:
      INTEGER    bond(mca,mca)
      INTEGER    rank(mca)  ! rank of nodes, in order supplied
! output: 
      INTEGER    rjg(mri,2) ! ring-join group pairs

! internal:
      INTEGER    ppb(mca,mca) ! product of primes of nodes in bond
      INTEGER    cprim(mca)   ! prime corresponding to node's rank
      INTEGER    tbond(mca,mca) ! working bond matrix.  Used twice:
                 ! zero 'boken' bonds; mirror for rearrangmenet
      INTEGER    i,j,ii,jj,k,maxp,n,rngflg
      INTEGER    top,sec      ! first and second nodes in new arrangement
      INTEGER    parent(mca),child(mca,3),nk(mca)
      INTEGER    iold(mca)    ! old index of group whose new index is 'k'
      INTEGER    trjg(nring,2)  ! temporary ring-join group pairs
      INTEGER    trank(mca)   ! temporary rank listing, for rearrangement
      CHARACTER(LEN=lgr) tgroup(mca) ! temporary group listing, for rearr.
      INTEGER    prim(30)    ! prime numbers
      INTEGER    rgpath(mca)
      DATA prim/  2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 
     &           31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 
     &           73, 79, 83, 89, 97,101,103,107,109,113/

!---------------------------
      !print*,'*uniqring*'

! initialise
      DO n=1,nring
        rjg(n,1)=0
        rjg(n,2)=0
      ENDDO
      DO i=1,nca
        tgroup(i)=' '
        DO j=1,nca
          tbond(i,j)=0
        ENDDO
      ENDDO

! find product of connected primes for actual bonds, setup connections 
      DO i=1,nca
        cprim(i)=prim(rank(i))
      ENDDO

      DO i=1,nca
        DO j=1,nca
          ppb(i,j) = 0
          tbond(i,j)=bond(i,j)
          IF(i.NE.j.AND.bond(i,j).GT.0) THEN
            ppb(i,j)=cprim(i)*cprim(j)
          ENDIF
        ENDDO
      ENDDO

! troubleshooting
      !DO j=1,nca
      !  write(6,'(a12,14(i2))')group(j),(tbond(i,j),i=1,nca)
      !ENDDO
      !print*,'primes bond matrix'
      !DO j=1,nca
      !  write(6,'(10(i4))') (ppb(i,j),i=1,nca)
      !ENDDO
      !print*,' '

! remove ring-joining characters from groups
      CALL rjgrm(nring,group,rjg)
      
! loop rings, finding lowest-priority on-ring connection to 'break' 
      DO n=1,nring
100     maxp = 0
        DO i=1,nca
          DO j=i+1,nca
! disallow breaking @ double bonds
            IF(bond(i,j).NE.2)THEN
              IF(ppb(i,j).GT.maxp) THEN
                ii = i
                jj = j
                maxp=ppb(ii,jj) 
              ENDIF
            ENDIF
          ENDDO
        ENDDO
        CALL findring(ii,jj,nca,tbond,rngflg,rgpath)
        ppb(ii,jj) = 0
        IF(rngflg.EQ.1) THEN ! connection is on a ring
          !print*,'breaking ring at',ii,jj
          tbond(ii,jj) = 0
          tbond(jj,ii) = 0
          rjg(n,1) = ii
          rjg(n,2) = jj
* symmetry might be broken after first bond breaking.
* therefore, call ratings to update ranks.          
          !IF (nring.GT.1) print*,"nring = ",nring,n
          IF (nring.GT.1) CALL ratings(nca,group,tbond,nring,rank)
        ELSE
          !print*,ii,jj,' not on a ring'
          GOTO 100
        ENDIF
      ENDDO ! n_rings

! find first 2 nodes in tree (linearly-expressed molecule)
      top = ii
      sec = 0
      DO j = 2,nca
        DO i = 1,j-1
          IF(tbond(i,j).GT.0)THEN
            !print*,i,j,tbond(i,j)
            IF(i.EQ.ii) sec=j
            IF(j.EQ.ii) sec=i
            IF(sec.GT.0) GOTO 200
          ENDIF
        ENDDO
      ENDDO

! rearrange into a connected skeleton (finding longest path comes later)
200   CONTINUE
      CALL findtree(tbond,top,sec,nca,parent,child)

! troubleshooting
      !DO i=1,nca
      !  print*,i,(child(i,k),k=1,3),parent(i)
      !ENDDO
      !print*,' '

      DO i=1,nca
        nk(i)=0
      ENDDO
      j=1
      iold(j)=top
      !  print*,iold(j),'old',0,'top',j,'top'
      IF(parent(iold(j)).NE.0)THEN
        !print*,'ERROR in uniqring : top has parent!'
        STOP
      ENDIF

      i=j
210   CONTINUE
      k = nk(i)+1
      IF(k.LE.3) THEN
        IF(child(iold(i),k).GT.0) THEN
!       print*,'a ',i,nk(i),j,iold(j),iold(i),k
!     &        ,parent(iold(i)),child(iold(i),k)
          nk(i) = nk(i)+1
          j=j+1
          iold(j)=child(iold(i),k)
          i=j ! step forward
          GO TO 210 ! try for next child
        ENDIF
      ENDIF
! otherwise, no further children on node
!       print*,'b ',i,nk(i),j,iold(j),iold(i),k
!     &        ,parent(iold(i)),child(iold(i),k)
      IF(parent(iold(i)).EQ.0) GO TO 212 ! finished thoroughly checking molecule
      DO ii=1,nca
        IF(parent(iold(i)).EQ.iold(ii))THEN
          i = ii
          GO TO 210 ! backtrack to parent
        ENDIF
      ENDDO
212   CONTINUE
      !print*,'finished sorting new tree'

! reassign arrays based on new ordering
      DO ii=1,nca
        IF(iold(ii).NE.0)THEN
          trank(ii) = rank(iold(ii))
          tgroup(ii) = group(iold(ii))
        ENDIF
        DO jj=1,nca
          IF(iold(jj).NE.0)THEN
            tbond(ii,jj) = bond(iold(ii),iold(jj))
          ELSE
! retain zeroes for broken bonds (e.g. for fragmentation that retains ring)
            tbond(ii,jj)=0
          ENDIF
        ENDDO
        DO n=1,nring
          DO i=1,2
            IF(rjg(n,i).EQ.iold(ii)) THEN
              trjg(n,i)=ii
            ENDIF
          ENDDO
        ENDDO
      ENDDO

! b: replace working arrays
      DO ii=1,nca
        rank(ii) = trank(ii)
        group(ii) = tgroup(ii)
        DO jj=1,nca
          bond(ii,jj) = tbond(ii,jj)
        ENDDO
      ENDDO
      DO n=1,nring
        rjg(n,1)=MIN0(trjg(n,1),trjg(n,2))
        rjg(n,2)=MAX0(trjg(n,1),trjg(n,2))
      ENDDO
     
! add ring-joining characters 
      CALL rjgadd(nring,group,rjg)

! troubleshooting
      !print*,'new BM '
      !DO j=1,nca
      !  write(6,'(a12,14(i2))')group(j),(bond(i,j),i=1,nca)
      !ENDDO
      !DO j=1,nca
      !  write(6,*)group(j),rank(j)
      !ENDDO
      !print*,''

      !print*,'*end uniqring*'

      RETURN
      END
