C***********************************************************************
      SUBROUTINE SUPPL(FW,WGT,WGTM,MAP,WGTP,IDXRES,KPAR,NRES,IRESET,
     +                 IPTR,MAXRES,ICNT,NREF,MAXREF,NRTCTL,IVALSP)
C***********************************************************************
C 10-Jun-08  EDG  Initialize routine to correct errors in beta analysis
C 25-Oct-01  CMM  Don't simultaneously delete & add a low weight structure.
C                 Don't add a reference structure more than once.
C  2-Jun-98  EDG  Fix several bugs
C 28-Jan-94  JKB  Add IVALSP=2 to restart NRT if no covalent structures and
C                    truncated matrix, else only ionic structures possible
C-----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON/NBIO/LFNIN,LFNPR,LFNAO,LFNPNA,LFNNAO,LFNPNH,LFNNHO,LFNPNB,
     +            LFNNBO,LFNPNL,LFNNLM,LFNMO,LFNDM,LFNNAB,LFNPPA,LFNARC,
     +            LFNDAF,LFNLBL,LFNDEF,LFNBRK(100)
      COMMON/NBOPT/IWDM,IW3C,IWAPOL,IWHYBS,IWPNAO,IWTNAO,IWTNAB,
     + IWTNBO,IWFOCK,IWCUBF,IPSEUD,KOPT,IPRINT,IWDETL,IWMULP,ICHOOS,
     + IWNBBP,IWMSP,IWFIXDM,IW3CHB,IWNJC,JCORE,JPRINT(100)
      COMMON/NBTHR/THRSET,PRJSET,ACCTHR,CRTSET,E2THR,ATHR,PTHR,ETHR,
     +             DTHR,DLTHR,CHSTHR,REFTHR,STTHR,PRTHR,THRNCS,THRNJC
      COMMON/NBFLAG/ROHF,UHF,CI,OPEN,COMPLX,ALPHA,BETA,MCSCF,AUHF,ORTHO
      LOGICAL ROHF,UHF,CI,OPEN,COMPLX,ALPHA,BETA,MCSCF,AUHF,ORTHO
C
      DIMENSION WGT(MAXRES,MAXREF),WGTP(MAXREF),WGTM(MAXRES*MAXREF),
     +          MAP(MAXRES*MAXREF),IDXRES(MAXRES,MAXREF),NRES(MAXREF),
     +          IRESET(MAXREF),IPTR(MAXREF),KPAR(MAXRES,MAXREF),
     +          NRTCTL(10),IDXTMP(9999)
C
      SAVE ZERO,EPS,FWMAX,KREF,IDXTMP,NTMP,NABTMP
      DATA ZERO,EPS,NTMP,IDXTMP,NABTMP/0.0D0,1.0D-5,0,9999*0,0/
C
C  Don't allow for cycles of the reference list if $NRTSTR is found
C  (NRTCTL(1)=-2):
C
      IF(NRTCTL(1).EQ.-2) RETURN
C
C  Initialize routine:
C
      IF(NRTCTL(1).EQ.-1) THEN
        FWMAX = ZERO
        KREF = 0
      END IF
C
C  Abort search for reference structures if oscillating:
C
      IVALSP = 0
      KCNT = 0
      KERR = 0
      IFORCE = 0
      DO 5 IREF = 1,MAXREF
        IF(IRESET(IREF).GE.0) KCNT = KCNT + 1
        IF(IRESET(IREF).LT.-1) KERR = KERR + 1
    5 CONTINUE
      IF(NRTCTL(1).LT.-1) THEN
        FWMAX = FW
        KREF  = KCNT
        IF(KCNT.EQ.1) FWMAX = ZERO
      ELSE
        NRTCTL(1) = 0
        IF(ABS(FW-FWMAX).LT.EPS.AND.KCNT.LE.KREF) THEN
          IF(KERR.EQ.0) RETURN
        END IF
        IF(FW.GT.FWMAX) THEN
          FWMAX = FW
          KREF  = KCNT
        END IF
      END IF
C
C  Reinitialize NTMP for beta spin analysis:
C
      IF(BETA) THEN
        IF(NABTMP.EQ.0) NTMP = 0
        NABTMP = 1
      END IF
C
C  Initialize search for additional reference structures:
C
C  IRESET(IREF):  Reference structure status flag
C        1 : new structure from last pass (accept according to wgt)
C        0 : old structure (accept according to wgt)
C       -1 : unoccupied position in reference list
C       -2 : structure with low weight (delete)
C       -3 : covalent structure from CONDNS (accept for next pass)
C       -4 : an ionic structure from CONDNS (delete)
C
      DO 10 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.0) IRESET(IREF) = 1
   10 CONTINUE
C
C  Find the resonance structure of highest weight:
C
      WGTMAX = ZERO
      DO 15 IREF = 1,NREF
        IF(IRESET(IREF).GE.0) THEN
          DO 14 IRES = 1,NRES(IREF)
            WGTMAX = MAX(WGTMAX,WGT(IRES,IREF)*WGTP(IREF))
   14     CONTINUE
        END IF
   15 CONTINUE
C
C  Count the number of structures in MAP (WGTM) having weight larger
C  than THRESH:
C
      JCNT = 0
      THRESH = ABS(REFTHR) * WGTMAX
      DO 20 I = 1,ICNT
        IF(WGTM(I).LT.THRESH) GOTO 30
        JCNT = I
   20 CONTINUE
   30 CONTINUE
C
C  Accept all reference structures with weight greater than THRESH:
C
      DO 50 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.1) THEN
          DO 40 I = 1,JCNT
            IF(IDXRES(1,IREF).EQ.ABS(MAP(I))) THEN
              IF(WGTM(I).GT.THRESH) THEN
                IRESET(IREF) = 0
                MAP(I) = -ABS(MAP(I))
                GOTO 50
              END IF
            END IF
   40     CONTINUE
          IFORCE = 1
        END IF
   50 CONTINUE
      N = 0
      DO 60 IREF = 1,MAXREF
        IF(IRESET(IREF).NE.0.AND.IRESET(IREF).NE.-1) THEN
          N = N + 1
          IPTR(N) = IREF
        END IF
   60 CONTINUE
      IF(N.GT.0) THEN
        IF(JPRINT(57).NE.0) THEN
          WRITE(LFNPR,910)
          WRITE(LFNPR,920) (IPTR(I),I=1,N)
        END IF
      END IF
C
C  CMM  If a structure is deleted, don't add it back.
C
      DO 62 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.-2) THEN
          DO 61 I = 1,JCNT
            IF(IDXRES(1,IREF).EQ.MAP(I)) MAP(I) = -ABS(MAP(I))
   61     CONTINUE
        END IF
   62 CONTINUE
C
C  Add covalent partners of the ionic structures deleted in CONDNS:
C  CMM Don't add in a structure for the second time.
C
      DO 65 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.-3) THEN
          IDLFLG = 0
          DO 63 ITMP = 1,NTMP
            IF(IDXRES(1,IREF).EQ.IDXTMP(ITMP)) IDLFLG = 1
   63     CONTINUE
          IF(IDLFLG.EQ.0) THEN
            NTMP = NTMP + 1
            IDXTMP(NTMP) = IDXRES(1,IREF)
            IRESET(IREF) = 2
            DO 64 I = 1,JCNT
              IF(MAP(I).EQ.IDXRES(1,IREF)) MAP(I) = -ABS(MAP(I))
   64       CONTINUE
          ELSE
            IRESET(IREF) = -4
          END IF
        END IF
   65 CONTINUE
C
C  CMM If a structure has been added in the past, don't add it again.
C
          DO 67 ITMP = 1,NTMP
            DO 66 I=1,JCNT
              IF(MAP(I).EQ.IDXTMP(ITMP)) MAP(I) = -ABS(MAP(I))
   66       CONTINUE
   67     CONTINUE
C
C  Add secondary structures of high weight to the reference list
C  (neglect secondary structures from KEKULE/APPEND):
C
      MORE = 0
      DO 80 I = 1,JCNT
        IF(MAP(I).GT.0) THEN
          DO 70 IREF = 1,MAXREF
            IF(IRESET(IREF).NE.-2.AND.
     +         IRESET(IREF).NE. 0.AND.
     +         IRESET(IREF).NE. 1.AND.
     +         IRESET(IREF).NE. 2) THEN
              DO 69 JREF = 1,MAXREF
                IF(IRESET(JREF).EQ.-2.OR.
     +             IRESET(JREF).EQ. 0.OR.
     +             IRESET(JREF).EQ. 1) THEN
                  DO 68 JRES = 1,NRES(JREF)
                    IF(KPAR(JRES,JREF).EQ.0) THEN
                      IF(IDXRES(JRES,JREF).EQ.MAP(I)) THEN
                        IF(WGTM(I).GT.THRESH) THEN
                          NTMP = NTMP + 1
                          IDXTMP(NTMP) = MAP(I)
                          IDXRES(1,IREF) = MAP(I)
                          IRESET(IREF) = 2
                          GOTO 80
                        END IF
                      END IF
                    END IF
   68             CONTINUE
                END IF
   69         CONTINUE
              GOTO 80
            END IF
   70     CONTINUE
          MORE = MORE + 1
        END IF
   80 CONTINUE
      IF(MORE.GT.0) THEN
        WRITE(LFNPR,900) MAXREF,MAXREF+MORE
        NRTCTL(1) = 0
        RETURN
      END IF
C
      DO 85 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.-2) IRESET(IREF) = -1
        IF(IRESET(IREF).EQ. 1) IRESET(IREF) = -1
   85 CONTINUE
C
C  Now search the secondary lists for structures of high weight:
C
      DO 140 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.2) THEN
C
C  First locate the secondary structure of highest weight:
C
          KRES = 0
          KREF = 0
          WGTMAX = ZERO
          DO 100 JREF = 1,MAXREF
            IF(IRESET(JREF).EQ.0) THEN
              DO 90 JRES = 2,NRES(JREF)
                IF(KPAR(JRES,JREF).EQ.0) THEN
                  IF(IDXRES(JRES,JREF).EQ.IDXRES(1,IREF)) THEN
                    WGTRES = WGTP(JREF) * WGT(JRES,JREF)
                    IF(WGTRES.GE.THRESH.AND.WGTRES.GT.WGTMAX) THEN
                      WGTMAX = WGTRES
                      KRES = JRES
                      KREF = JREF
                    END IF
                  END IF
                END IF
   90         CONTINUE
            END IF
  100     CONTINUE
C
C  Insert the KRESth structure of the KREFth reference manifold into
C  the IREFth location of the reference list:
C
          IRESET(IREF) = 1
          NRES(IREF)   = 1
          KPAR(1,IREF) = KREF*10000 + KRES
          IF(KPAR(1,IREF).EQ.0) GOTO 140
C
C  If the KRESth structure of KREF came from KEKULE, mark it so that
C  it won't be used again:
C
          IF(KPAR(KRES,KREF).EQ.0) KPAR(KRES,KREF) = 1
C
C  And search the KREFth manifold for other structures of identical
C  TOPO matrix with high weight.  Also add these to the reference list:
C
          DO 130 JRES = 2,NRES(KREF)
            IF(JRES.NE.KRES) THEN
              IF(KPAR(JRES,KREF).EQ.0) THEN
                IF(IDXRES(JRES,KREF).EQ.IDXRES(1,IREF)) THEN
                  WGTRES = WGTP(KREF) * WGT(JRES,KREF)
                  IF(WGTRES.GT.THRESH) THEN
                    DO 110 JREF = 1,MAXREF
                      IF(IRESET(JREF).LE.-1) THEN
                        IDXRES(1,JREF) = IDXRES(1,IREF)
                        NRES(JREF)   = 1
                        IRESET(JREF) = 1
                        KPAR(1,JREF) = KREF*10000 + JRES
                        GOTO 120
                      END IF
  110               CONTINUE
                    CALL NBHALT('Increase NRTMEM in the $NBO keylist.')
  120               CONTINUE
                  END IF
                END IF
              END IF
            END IF
  130     CONTINUE
        END IF
  140 CONTINUE
C
      N = 0
      DO 150 IREF = 1,MAXREF
        IF(IRESET(IREF).EQ.1) THEN
          NRTCTL(1) = 1
          N = N + 1
          IPTR(N) = IREF
        END IF
  150 CONTINUE
      IF(IFORCE.EQ.1) NRTCTL(1) = 1
      IF(N.GT.0) THEN
        IF(JPRINT(57).NE.0) THEN
          WRITE(LFNPR,930)
          WRITE(LFNPR,920) (IPTR(I),I=1,N)
        END IF
      END IF
C
C  Reset:
C
      DO 160 IREF = 1,MAXREF
        IRESET(IREF) = MAX(IRESET(IREF),-1)
  160 CONTINUE
      NREF = 0
      DO 170 I = 1,ICNT
        MAP(I) = ABS(MAP(I))
  170 CONTINUE
      DO 180 IREF = MAXREF,1,-1
        IF(NREF.EQ.0.AND.IRESET(IREF).GE.0) THEN
          NREF = IREF
          RETURN
        END IF
  180 CONTINUE
C
C If SUPPL got this far without RETURNing, only ionic structures were
C    found and no corresponding covalent structure was generated.  Try
C    restarting with full AO matrix if not done already, else quit.
C
      IF(JPRINT(77).EQ.0) THEN
        WRITE(LFNPR,940)
        IVALSP = 2
      ELSE
        WRITE(LFNPR,950)
        IVALSP = 3
      END IF
      RETURN
C
  900 FORMAT(/1X,'The NRT program is currently configured to handle ',
     + I2,' reference structures',/1X,'which is insufficient for the',
     + ' present calculation.  Increase this number',/1X,'to least ',
     + I2,' using the NRTMEM keyword of the $NBO keylist.')
  910 FORMAT(/1X,'Deleting the following reference structures:')
  920 FORMAT(1X,19I4)
  930 FORMAT(/1X,'Adding the following reference structures:')
  940 FORMAT(/1X,'No covalent structures available within minimal ',
     +   'valence space (SR SUPPL).',/1X,'Reallocate scratch vector',
     +   ' and restart using full AO density matrix.')
  950 FORMAT(/1X,'No covalent reference structures found; ',
     +   'abandon NRT analysis.')
      END