C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C  7 APR 06 - MWS - REMOVE NON-ARGUMENT ADJUSTABLE ARRAYS
C 29 MAR 06 - FLG - LOCALIZATION ROUTINES FOR ELONGATION METHOD
C
C*MODULE ELGLOC  *DECK CHKLOCL
      SUBROUTINE CHKLOCL(C,S,NOCC1,NOCC2,NVAC1,FF,
     *          FG,IFIRST,IEND,JFIRST,JEND,NUM,N,NJGK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXAO=8192)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /MOSEL / ISEL(MXAO,2),JSEL(MXAO,2)
      DIMENSION C(NUM,N),S(NUM,NUM)
      DIMENSION FF(NUM),FG(NUM)
C
C    IMXO1: WORST LOCL BAND IN A REGION OF OCC BLOCK
C    IMXO2: WORST LOCL BAND IN B REGION OF OCC BLOCK
C
C    IMXV1: WORST LOCL BAND IN A REGION OF VIRT BLOCK
C    IMXV2: WORST LOCL BAND IN B REGION OF VIRT BLOCK
C
C    A REGION: NOCC1+NVAC1
C    A + B   : N
C
      NOCC = NOCC1 + NOCC2
      NOVT = NOCC + NVAC1
C
      IF(MASWRK) THEN
      WRITE(6,*)
      WRITE(6,*)'LOCALIZATION QUALITY: <A|B>, <B|A> SHOULD BE SMALL'
      WRITE(6,*)
      WRITE(6,*)' I         <A|B>           <B|A>'
      ENDIF
C
      CALL VCLR(FF,1,NUM)
      CALL VCLR(FG,1,NUM)
C
      FLOC1 = 0.0D+00
      FLOC2 = 0.0D+00
      FLVT1 = 0.0D+00
      FLVT2 = 0.0D+00
C
      FAB = 0.0D+00
      FBA = 0.0D+00
C
      DO 70 I=1,N
         IS=ISEL(I,NJGK)
         JS=JSEL(I,NJGK)
C
         IF(IS.NE.0) THEN
C           <A|B>
            DO 40 I1=JFIRST,JEND
               DO 40 I2=JFIRST,JEND
                  FF(IS)=FF(IS)+C(I1,IS)*C(I2,IS)*S(I1,I2)
 40         CONTINUE
C
            IF(IS.LE.NOCC1) THEN
               IF(FF(IS).GT.FLOC1) THEN
                  FLOC1 = FF(IS)
                  IMXO1 = IS
               ENDIF
            ELSEIF(IS.LE.NOVT.AND.IS.GT.NOCC) THEN
               IF(FF(IS).GT.FLVT1) THEN
                  FLVT1 = FF(IS)
                  IMXV1 = IS
               ENDIF
            ENDIF
         ENDIF
C
         IF(JS.NE.0) THEN
C           <B|A>
            DO 60 I1=IFIRST,IEND
               DO 60 I2=IFIRST,IEND
                  FG(JS)=FG(JS)+C(I1,JS)*C(I2,JS)*S(I1,I2)
 60         CONTINUE
C
            IF(JS.GT.NOCC1.AND.JS.LE.NOCC) THEN
               IF(FG(JS).GT.FLOC2) THEN
                  FLOC2 = FG(JS)
                  IMXO2 = JS
               ENDIF
            ELSEIF(JS.LE.N.AND.JS.GT.NOVT) THEN
               IF(FG(JS).GT.FLVT2) THEN
                  FLVT2 = FG(JS)
                  IMXV2 = JS
               ENDIF
            ENDIF
C
         ENDIF
C
         FAB = FAB + FF(I)
         FBA = FBA + FG(I)
         IF(I.EQ.NOCC) THEN
            FABOC = FAB
            FBAOC = FBA
         END IF
C
         IF(MASWRK) WRITE(6,1)I,FF(I),FG(I)
C
 70   CONTINUE
C
      IF(MASWRK) THEN
      WRITE(6,*)
      WRITE(6,'(''TOTAL <A|B>: '',F15.10,5X,'' TOTAL <B|A>: '',
     *            F15.10)')FAB,FBA
      WRITE(6,'('' OCC  <A|B>: '',F15.10,5X,''  OCC  <B|A>: '',
     *            F15.10)')FABOC,FBAOC
      WRITE(6,'('' VAC  <A|B>: '',F15.10,5X,''  VAC  <B|A>: '',
     *            F15.10)')FAB-FABOC,FBA-FBAOC
      WRITE(6,*)
      WRITE(6,'(''WORST LMOS IN A REGION: IMXO1, FLOC1='',I3,F15.10,
     *  '', IMXV1,FLVT1='',I3,F15.10)') IMXO1,FLOC1,IMXV1,FLVT1
      WRITE(6,'(''WORST LMOS IN B REGION: IMXO2, FLOC2='',I3,F15.10,
     *  '', IMXV2,FLVT2='',I3,F15.10)') IMXO2,FLOC2,IMXV2,FLVT2
C
      WRITE(6,*)
      WRITE(6,*)'FOUR WORST LMOS LOCALIZATION QUALITY'
      WRITE(6,*)
      WRITE(6,*)' I         <A|B>           <B|A>'
C
      WRITE(6,1)IMXO1,FF(IMXO1),FG(IMXO1)
      WRITE(6,1)IMXO2,FF(IMXO2),FG(IMXO2)
      WRITE(6,1)IMXV1,FF(IMXV1),FG(IMXV1)
      WRITE(6,1)IMXV2,FF(IMXV2),FG(IMXV2)
      ENDIF
C----------------------------------------------------------------------
      RETURN
 1    FORMAT(I3,1X,2(F15.12,1X))
      END
C
C*MODULE ELGLOC  *DECK GULOCL
      SUBROUTINE GULOCL
C
C     REGIONAL ATOMIC ORBITAL LOCALIZATION SCHEME FOR ELONGATION
C
C     WRITTEN BY FENG LONG GU: APRIL 2004
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      LOGICAL LCUT
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      COMMON /ELGFLS /NEG0,NEG1,NEG2,NEG3,NEG4,NEG5,NEG6,NEG7,NEG8,NEG9
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
      COMMON /ELGINF/ NAAO,NCUTP,NCUT,NCUTT,NLATOT,NLBTOT,NLA,NLB,
     *                IMOCUT,IMOCUTT,NWA,NWM,NWMT,NHB
      COMMON /MOSEL / ISEL(MXAO,2),JSEL(MXAO,2)
      COMMON /ELGIDX/ LCUT
      COMMON /NEWPAR/ NUMNEW,NATNEW,NSHNEW,ISNEW
      COMMON /ITST/ N,N1,NOCC,NOCC1,NVAC1,NOC1X,NVC1X,NOCC2,NVAC2,IEND,
     *              NAAOX,JXXX,JWWX,NVAC
C
      PARAMETER (MXSTP=500)
      DATA UHF/8HUHF     /
C
      DIMENSION NFZATM(MXSTP)
C
      IF(MASWRK) WRITE(IW,9000)
C
         N1 = NUM
         N  = NUM
       IF(NELONG.EQ.3)  N  = NWM
C
      KELGX = 0
      JXXX = 0
      JWWX = 0
      NCUTX = 0
      NAAOX = 0
      NOC1X = 0
      NOC2X = 0
      NOC1XB = 0
      NOC2XB = 0
      NCTAOX = 0
      NCTMOX = 0
C
      IF(NELONG.EQ.3) THEN
         CALL SEQREW(NEG1)
         READ(NEG1) JXXX,JWWX,NAAOX,NOC1X,NOC2X,NOC1XB,NOC2XB,MULX,
     *              NHBX,NATX
         READ(NEG1) NATOCX,NCTMOX,NCUTX,NCTAOX,KELGX,NSHELX,NDNSLX
         READ(NEG1)
         READ(NEG1)
         READ(NEG1) (NFZATM(I),I=1,KELGX)
      ENDIF
C
      IF(NEG1.LT.0) THEN
         ISUM = NSHELX + NDNSLX + NATX + NATOCX + NCTAOX + NHBX + MULX
         WRITE(6,*) 'STRICTLY TO FOOL FTNCHEK: ISUM=',ISUM
      END IF
C
      NBTM = NAT - NATM
      IF(MASWRK)
     * WRITE(IW,'('' TOTAL NUMBER OF ATOMS: '',I4,/,
     *            ''       ATOMS IN PART A: '',I4,/,
     *            ''       ATOMS IN PART B: '',I4)')NAT,NATM,NBTM
C
      KELG = KELGX + 1
      IF(KELG.EQ.1) NCUT = 0
      IF(MASWRK) WRITE(IW,'('' NUMBER OF ELONGATION STEPS: '',I4)') KELG
C
      L1 = NUM
      L2 = (L1*(L1+1))/2
      L3 = L1*L1
C
C     ----- GET MEMORY -----
C
      LOADFM=0
      CALL VALFM(LOADFM)
      LCAL  = LOADFM + 1
      LS    = LCAL   + N*N1
      LT    = LS     + L3
      LST   = LT     + L3
      LFF   = LST    + L2
      LFG   = LFF    + L1
C
C     FOR RAO-LMO
C
      LTMP  = LFG    + L1
      LHAO  = LTMP   + L3
      LFAO  = LHAO   + L2
      LDOF  = LFAO   + L2
      LDBOF = LDOF   + L2
      LDA   = LDBOF  + L2
      LDB   = LDA    + L2
      LVEC  = LDB    + L2
      LEPS  = LVEC   + L3
      LVS   = LEPS   + L1
      LVV   = LVS    + L3
      LX    = LVV    + L3
      LY    = LX     + L3
      LES   = LY     + L3
      LCXA  = LES    + L1            ! FOR PREVIOUS A-REGION LMO
C
C       WRITE TO TAPE NEG5: B-REGION AND NEW A-REGION LMO
C
      LCWAL = LCXA   + L3
      LCXAL = LCWAL  + L3
      LAST1 = LCXAL  + L3
C
      IF(SCFTYP.EQ.UHF) THEN
       LCBT = LAST1
       LCXB = LCBT   + N*N1
       LCWBT= LCXB   + L3
       LCXBT= LCWBT  + L3
       LAST1= LCXBT  + L3
      ELSE
       LCBT = LCAL
       LCXB = LCXA
       LCWBT= LCWAL
       LCXBT= LCXAL
      ENDIF
C
C     SET SOME MEMORY FOR MOS REMOVAL DUE TO TERMINAL ATOMS
C
      LAST  = LAST1  + 5*NUM*NUM+4*NUM
C
       IF(LDOS.NE.0) THEN
        NPTS = 1001
        LDEN1 = LAST
        LDEN2 = LDEN1 + N1*NPTS
        LDEN3 = LDEN2 + NAT*NPTS
        LAST  = LDEN3 + NPTS
       ENDIF
C
      IF(LCUT) LAST  = LAST  + 2*(NUMNEW*NUMNEW+NUMNEW)
C
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C   INITIALIZATION
C
      DO 10 J=1,2
         DO 10 I=1,NUM
            ISEL(I,J) = 0
            JSEL(I,J) = 0
 10   CONTINUE
C
C      ----- READ OVERLAP FROM DAF -----
      CALL DAREAD(IDAF,IODA,XX(LST),L2,12,0)
      CALL CPYTSQ(XX(LST),XX(LS),NUM,1)
C
      CALL REGA(NATM,NBF1,NTMLB,NHB,NOCC1,NVAC1,NOCC2,NVAC2,
     *            NOCC1B,NVAC1B,NOCC2B,NVAC2B,NASPIN,NBNDAB)
      NOCC  = NOCC1  + NOCC2
      NOCCB = NOCC1B + NOCC2B
C----------------------------------------------------------------------
      IF(MASWRK) THEN
      WRITE(IW,'(40(1H-))')
      WRITE(IW,*)'         OCCUPIED MOLECULAR ORBITALS'
      WRITE(IW,'(40(1H-))')
      WRITE(IW,*)'                 ALPHA           BETA   '
      WRITE(IW,'(40(1H-))')
      WRITE(IW,'('' REGION A:'',8X,I4,11X,I4)') NOCC1,NOCC1B
      WRITE(IW,'('' REGION B:'',8X,I4,11X,I4)') NOCC2,NOCC2B
      WRITE(IW,'(40(1H-))')
      WRITE(IW,'(''TOTAL A+B:'',8X,I4,11X,I4)') NOCC,NOCCB
      WRITE(IW,'(40(1H-))')
      WRITE(IW,*)'        VACANT MOLECULAR ORBITALS'
      WRITE(IW,'(40(1H-))')
      WRITE(IW,'('' REGION A:'',8X,I4,11X,I4)') NVAC1,NVAC1B
      WRITE(IW,'('' REGION B:'',8X,I4,11X,I4)') NVAC2,NVAC2B
      WRITE(IW,'(40(1H-))')
      ENDIF
C----------------------------------------------------------------------
      NVC1X  = JXXX - NOC1X
      NVC1XB = JXXX - NOC1XB
      NVC2X  = JWWX - NOC2X
      NVC2XB = JWWX - NOC2XB
C
C     EXPLANATION OF INPUT DATA
C     NATM.... THE NO. OF ATOMS IN A-REGION
C     NTMLB....THE NO. OF TERMINAL ATOMS IN B-REGION
C     NBF1.... THE NO. OF BASIS FUNCTIONS IN A-REGION
C     NHB....  THE NO. OF BASIS FUNCTIONS IN TERMINAL ATOMS OF B-REGION
C
      NOCC1  = NOCC1  - NOC1X
      NVAC1  = NVAC1  - NVC1X
      NOCC1B = NOCC1B - NOC1XB
      NVAC1B = NVAC1B - NVC1XB
C
      IF(MASWRK) THEN
      WRITE(IW,*)'       REGION A - EXACT DESCRIPTION'
      WRITE(IW,'(40(1H-))')
      WRITE(IW,'('' FROZEN OCCUPIED: '',I4,11X,I4)') NOC1X,NOC1XB
      WRITE(IW,'('' FROZEN   VACANT: '',I4,11X,I4)') NVC1X,NVC1XB
      WRITE(IW,'('' ACTIVE OCCUPIED: '',I4,11X,I4)') NOCC1,NOCC1B
      WRITE(IW,'('' ACTIVE   VACANT: '',I4,11X,I4)') NVAC1,NVAC1B
      WRITE(IW,'(40(1H-))')
      ENDIF
C
      DO 20 I=1,NOCC1
         ISEL(I,1)=I
 20   CONTINUE
      DO 30 I=1,NVAC1
         ISEL(NOCC1+I,1)=I+NOCC1+NOCC2
 30   CONTINUE
      DO 40 I=1,NOCC2
         JSEL(I,1)=I+NOCC1
 40   CONTINUE
      DO 50 I=1,NVAC2
         JSEL(NOCC2+I,1)=I+NOCC1+NOCC2+NVAC1
 50   CONTINUE
C
      IF(SCFTYP.EQ.UHF) THEN
        DO I=1,NOCC1B
           ISEL(I,2)=I
        ENDDO
        DO I=1,NVAC1B
         ISEL(NOCC1B+I,2)=I+NOCC1B+NOCC2B
        ENDDO
        DO I=1,NOCC2B
           JSEL(I,2)=I+NOCC1B
        ENDDO
        DO I=1,NVAC2B
           JSEL(NOCC2B+I,2)=I+NOCC1B+NOCC2B+NVAC1B
        ENDDO
      ENDIF
C
C    IFIRST,IEND,JFIRST,JEND CAN BE OBTAINED FRON NBF1, NUM
C
      IFIRST = 1
      IEND = NBF1
      JFIRST = IEND + 1
      JEND = NUM
C
      CALL EGSAVE(XX(LCAL),XX(LCBT),XX(LTMP),N1,N,NOCC1,NVAC1,NOCC2,
     *   NVAC2,NOCC1B,NVAC1B,NOCC2B,NVAC2B,NOC1X,NVC1X,NOC1XB,NVC1XB)
C---------------------------------------------------------------------
      IF(IPRI.EQ.1.OR.IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(IW,'('' LCAO COEFFICIENTS BEFORE LOCALIZATION'')')
         CALL WRITN4(N,N1,XX(LCAL),N,N1)
         IF(SCFTYP.EQ.UHF) THEN
            WRITE(IW,'('' LCAO COEFFICIENTS --- BETA PART---'')')
            CALL WRITN4(N,N1,XX(LCBT),N,N1)
         ENDIF
      END IF
C
C   CALCULATE LDOS
C
      IF(LDOS.NE.0) THEN
       WRITE(6,*)'NOC1X,NVC1X,JXXX,NOC2X,NVC2X=',
     *            NOC1X,NVC1X,JXXX,NOC2X,NVC2X
       CALL DAREAD(IDAF,IODA,XX(LEPS),L1,17,0)
       IF(NELONG.LE.2) THEN
        DO I=1,L1
         XX(LES+I-1) = XX(LEPS+I-1)
        ENDDO
       ELSE
        DO I=1,NOCC-NOC1X
         XX(LES+I-1) = XX(LEPS+I+NOC1X-1)
        ENDDO
        DO I=1,L1-NOCC-NVC1X
         XX(LES+NOCC-NOC1X+I-1) = XX(LEPS+I+NOCC+NVC1X-1)
        ENDDO
       ENDIF
       CALL ELGDOS(XX(LES),XX(LCAL),XX(LS),XX(LDEN1),
     *       XX(LDEN2),XX(LDEN3),N1,N,NPTS,XX(LTMP))
      ENDIF
C---------------------------------------------------------------------
      IF(MASWRK) WRITE(IW,9001)
C
      CALL VALFM(LOADFM)
      LWRK  = LOADFM   + 1
      LIWRK = LWRK     + N1*8
      LAST2 = LIWRK    + N1
      NEED2 = LAST2 - LOADFM - 1
      CALL GETFM(NEED2)
      CALL SDGX(XX(LST),XX(LX),XX(LY),XX(LES),XX(LVV),N1,L2,
     *          XX(LWRK),XX(LIWRK))
      CALL RETFM(NEED2)
C
      NWRK = MAX(IEND,N1-IEND)
      CALL VALFM(LOADFM)
      LWRK  = LOADFM   + 1
      LIWRK = LWRK     + NWRK
      LAST2 = LIWRK    + NWRK
      NEED2 = LAST2 - LOADFM - 1
      CALL GETFM(NEED2)
C
C      RHF/ROHF, FOR OCCUPIED CMOS
C
         CALL SEQREW(NEG9)
         CALL DAO(XX(LDOF),XX(LCAL),NOCC,NOC1X,N1,N)
         CALL RAOLMO(XX(LDOF),XX(LDA),XX(LDB),XX(LCAL),XX(LT),
     *        N,N1,NOCC,NOCC1,NVAC1,NOC1X,NVC1X,NOCC2,L2,
     *        XX(LVEC),XX(LEPS),XX(LX),XX(LVS),XX(LVV),XX(LY),XX(LES),
     *        IEND,'OCC',XX(LWRK),XX(LIWRK),NWRK)
C
C      RHF/ROHF, FOR VACANT CMOS
C
         NVAC = N1 - NOCC
         CALL DAOVAC(XX(LDOF),XX(LCAL),NOCC,NOC1X,N1,N)
         CALL RAOLMO(XX(LDOF),XX(LDA),XX(LDB),XX(LCAL),XX(LT),
     *        N,N1,NVAC,NVAC1,NOCC1,NVC1X,NOC1X,NVAC2,L2,
     *        XX(LVEC),XX(LEPS),XX(LX),XX(LVS),XX(LVV),XX(LY),XX(LES),
     *        IEND,'VAC',XX(LWRK),XX(LIWRK),NWRK)
C
C      UHF, FOR BETA-OCCUPIED CMOS
C
      IF(SCFTYP.EQ.UHF) THEN
         CALL DAO(XX(LDOF),XX(LCBT),NOCCB,NOC1XB,N1,N)
         CALL RAOLMO(XX(LDOF),XX(LDA),XX(LDB),XX(LCBT),XX(LT),
     *        N,N1,NOCCB,NOCC1B,NVAC1B,NOC1XB,NVC1XB,NOCC2B,L2,
     *        XX(LVEC),XX(LEPS),XX(LX),XX(LVS),XX(LVV),XX(LY),XX(LES),
     *        IEND,'OCC',XX(LWRK),XX(LIWRK),NWRK)
C
C      UHF, FOR BETA-VACANT CMOS
C
         NVACB = N1 - NOCCB
         CALL DAOVAC(XX(LDOF),XX(LCBT),NOCCB,NOC1XB,N1,N)
         CALL RAOLMO(XX(LDOF),XX(LDA),XX(LDB),XX(LCBT),XX(LT),
     *        N,N1,NVACB,NVAC1B,NOCC1B,NVC1XB,NOC1XB,NVAC2B,L2,
     *        XX(LVEC),XX(LEPS),XX(LX),XX(LVS),XX(LVV),XX(LY),XX(LES),
     *        IEND,'VAC',XX(LWRK),XX(LIWRK),NWRK)
      ENDIF
      CALL RETFM(NEED2)
CGU====================================================================
      IF(MASWRK) WRITE(IW,*)'AFTER LOCALIZATION - ALPHA PART'
      CALL CHKLOCL(XX(LCAL),XX(LS),NOCC1,NOCC2,NVAC1,XX(LFF),
     *     XX(LFG),IFIRST,IEND,JFIRST,JEND,NUM,N,1)
C
      IF(SCFTYP.EQ.UHF) THEN
         IF(MASWRK) WRITE(IW,*)'AFTER LOCALIZATION - BETA PART'
         CALL CHKLOCL(XX(LCBT),XX(LS),NOCC1,NOCC2,NVAC1,XX(LFF),
     *        XX(LFG),IFIRST,IEND,JFIRST,JEND,NUM,N,2)
      ENDIF
C
C     PRINTING OF THE LOCALIZED ORBITAL
C
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(IW,*)
         WRITE(IW,*)'ALPHA SPIN-LCAO COEFFICIENTS AFTER LOCALIZATION'
         CALL WRITN4(N,N1,XX(LCAL),N,N1)
         IF(SCFTYP.EQ.UHF) THEN
            WRITE(IW,*)
            WRITE(IW,*)'BETA SPIN-LCAO COEFFICIENTS AFTER LOCALIZATION'
            CALL WRITN4(N,N1,XX(LCBT),N,N1)
         ENDIF
      END IF
C======================================================================
C   PREPARATION FOR NEXT ELONGATION STEP
C======================================================================
C   WRITE(6,*)'N1,N,JXXX,JWWX,NAAOX,NHBX,NCUTX,NCTAOX,KELG,KELGX,MULX'
C   WRITE(6,*) N1,N,JXXX,JWWX,NAAOX,NHBX,NCUTX,NCTAOX,KELG,KELGX,MULX
C   WRITE(6,*)'NATOCX,NCTMOX,NOCC,NOCCB,NVAC,NVACB,NBF1'
C   WRITE(6,*) NATOCX,NCTMOX,NOCC,NOCCB,NVAC,NVACB,NBF1
C   WRITE(6,*)'NOCC1,NOCC1B,NVAC1,NVAC1B,NOCC2,NOCC2B,NVAC2,NVAC2B'
C   WRITE(6,*) NOCC1,NOCC1B,NVAC1,NVAC1B,NOCC2,NOCC2B,NVAC2,NVAC2B
C   WRITE(6,*)'NOC1X,NOC1XB,NOC2X,NOC2XB,NVC1X,NVC1XB,NVC2X,NVC2XB'
C   WRITE(6,*) NOC1X,NOC1XB,NOC2X,NOC2XB,NVC1X,NVC1XB,NVC2X,NVC2XB
C
      CALL ELGPREP(XX(LHAO),XX(LFAO),XX(LS),XX(LCAL),XX(LCBT),
     *            XX(LCXA),XX(LCXB),XX(LCWAL),XX(LCWBT),
     *            XX(LCXAL),XX(LCXBT),XX(LDOF),XX(LDBOF),
     *            NFZATM,LAST1,LAST,
     *            IFIRST,IEND,JFIRST,JEND,N1,N,JXXX,NAAOX,
     *            NCUTX,KELG,KELGX,NCTMOX,
     *            NOCC,NOCCB,NBF1,
     *            NOCC1,NOCC1B,NVAC1,NVAC1B,
     *            NOCC2,NOCC2B,NVAC2,NVAC2B,
     *            NOC1X,NOC1XB,NOC2X,NOC2XB,
     *            NVC1X,NVC1XB,NVC2X,NVC2XB)
C
      CALL RETFM(NEED)
C
      RETURN
 9000 FORMAT(/5X,49(1H=)/5X,
     * '=                                               ='/5X,
     * '= ELONGATION METHOD FOR (NON-)PERIODIC POLYMERS ='/5X,
     * '=                                               ='/5X,
     * 49(1H=),/)
 9001 FORMAT(/5X,49(1H-)/5X,
     * '-                                               -'/5X,
     * '-   RAO-LMO LOCALIZATION METHOD FOR ELONGATION  -'/5X,
     * '-                                               -'/5X,
     * 49(1H-),/)
      END
C
C*MODULE ELGLOC  *DECK ELGPREP
      SUBROUTINE ELGPREP(HAO,FAO,SW,CAALPH,CABETA,CXA,CXB,CWALPH,CWBETA,
     *     CXALPH,CXBETA,DOFF,DBOFF,NFZATM,LAST1,LAST,
     *     IFIRST,IEND,JFIRST,JEND,NAAO,NAMO,JXXX,NAAOX,
     *     NCUTX,KELG,KELGX,NCTMOX,
     *     NOCC,NOCCB,NBF1,
     *     NOCC1,NOCC1B,NVAC1,NVAC1B,
     *     NOCC2,NOCC2B,NVAC2,NVAC2B,
     *     NOC1X,NOC1XB,NOC2X,NOC2XB,
     *     NVC1X,NVC1XB,NVC2X,NVC2XB)
C
C     PROGRAM TO INITIATE THE THEORETICAL SYNTHESIS OF POLYMERS
C     AT THE DENSITY FUNCTIONAL LEVEL
C          1992.12.31----
C                                        AT HEIDELBERG-DKFZ
C     RE-WRITTEN BY F. L. GU     MAY, 2003
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
      LOGICAL GOPARR,DSKWRK,MASWRK,DIRSCF,FDIFF
      LOGICAL LCUT,LTRMST
C
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      COMMON /ELGFLS /NEG0,NEG1,NEG2,NEG3,NEG4,NEG5,NEG6,NEG7,NEG8,NEG9
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
      COMMON /ELGINF/ NNAAO,NCUTP,NCUT,NCUTT,NLATOT,NLBTOT,NLA,NLB,
     *                IMOCUT,IMOCUTT,NWA,NWM,NWMT,NHB
      COMMON /ELGIDX/ LCUT
      COMMON /ELGTRM/ LTRMST,NFLTRM,NRCTRM,NPSTRM,NHTSHL
CGU   COMMON /FLDCUT/ EFCUT(3), IGFLD
C
      DIMENSION CAALPH(NAAO,NAMO),CABETA(NAAO,NAMO)
      DIMENSION SW(NUM,NUM),DOFF(NUM*(NUM+1)/2),DBOFF(NUM*(NUM+1)/2)
      DIMENSION CXA(NUM,NUM),CXB(NUM,NUM)
      DIMENSION CWALPH(NUM,NUM),CWBETA(NUM,NUM)
      DIMENSION CXALPH(NUM,NUM),CXBETA(NUM,NUM)
      DIMENSION FAO(NUM*(NUM+1)/2),HAO(NUM*(NUM+1)/2)
      PARAMETER (MXSTP=500)
      DIMENSION NFZOC(MXSTP),NFZVT(MXSTP),NFZOCB(MXSTP),NFZVTB(MXSTP)
      DIMENSION NFZATM(MXSTP)
C-----------------------------------------------------------------------
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
C-----------------------------------------------------------------------
      IF(MASWRK) THEN
      WRITE(IW,*)
      WRITE(IW,'(''IFIRST,IEND,JFIRST,JEND,NAAO,NAMO='',6(I4,2X))')
     *             IFIRST,IEND,JFIRST,JEND,NAAO,NAMO
C
      WRITE(IW,*) '     -------ALPHA ELECTRONS-------'
      WRITE(IW,'(''       NOCC1,NOCC2,NVAC1,NVAC2='',4(I4,2X))')
     *                    NOCC1,NOCC2,NVAC1,NVAC2
      IF(SCFTYP.NE.RHF) THEN
      WRITE(IW,*) '     -------BETA ELECTRONS-------'
      WRITE(IW,'(''   NOCC1B,NOCC2B,NVAC1B,NVAC2B='',4(I4,2X))')
     *                NOCC1B,NOCC2B,NVAC1B,NVAC2B
      ENDIF
C
      IF(KELG.GE.2) THEN
      WRITE(IW,*)
      WRITE(IW,*) '     -------ALPHA ELECTRONS (ON TAPE)-------'
      WRITE(IW,'(''       NOC1X,NOC2X,NVC1X,NVC2X='',4(I4,2X))')
     *                    NOC1X,NOC2X,NVC1X,NVC2X
      IF(SCFTYP.NE.RHF) THEN
      WRITE(IW,*) '     -------BETA ELECTRONS (ON TAPE)-------'
      WRITE(IW,'(''   NOC1XB,NOC2XB,NVC1XB,NVC2XB='',4(I4,2X))')
     *                NOC1XB,NOC2XB,NVC1XB,NVC2XB
      ENDIF
      ENDIF
      ENDIF
C-----------------------------------------------------------------------
      CALL VCLR(CXA,1,NUM*NUM)
      CALL VCLR(CWALPH,1,NUM*NUM)
      CALL VCLR(CXALPH,1,NUM*NUM)
C
      IF(SCFTYP.EQ.UHF) THEN
        CALL VCLR(CXB,1,NUM*NUM)
        CALL VCLR(CWBETA,1,NUM*NUM)
        CALL VCLR(CXBETA,1,NUM*NUM)
      ENDIF
C
      IF(KELG.GT.1) THEN
         CALL SEQREW(NEG5)
         CALL SEQREW(NEG1)
         READ(NEG1)
         READ(NEG1)
         READ(NEG1) (NFZOC(I),I=1,KELGX),(NFZVT(I),I=1,KELGX)
         READ(NEG1) (NFZOCB(I),I=1,KELGX),(NFZVTB(I),I=1,KELGX)
         READ(NEG1) (NFZATM(I),I=1,KELGX)
C        PREVIOUS A AND B REGION
                           READ(NEG5)((CXA(ID,JD),ID=1,NAAOX),JD=1,JXXX)
         IF(SCFTYP.EQ.UHF) READ(NEG5)((CXB(ID,JD),ID=1,NAAOX),JD=1,JXXX)
      END IF
C
      NVC1X  = JXXX - NOC1X
      NVC1XB = JXXX - NOC1XB
C-----------------------------------------------------------------------
C
C   SAVE FROZEN OCC AND VIRT
C
      NFZOC(KELG) = NOCC1
      NFZVT(KELG) = NVAC1
      NFZOCB(KELG) = NOCC1B
      NFZVTB(KELG) = NVAC1B
      IF(MASWRK) THEN
      IF(KELG.EQ.1) THEN
         WRITE(IW,*)'---- COMPUTING A1, A2, ... ----'
         WRITE(IW,'('' KELG = '',I4)')KELG
         WRITE(IW,'(''NOCC1 = '',I4,'' NVAC1 = '',I4)')NOCC1,NVAC1
         NATCT = NATCUT(1,NOCC1+NVAC1)
         WRITE(IW,'(''NATCT = '',I4)')NATCT
      ELSE
         WRITE(IW,*)'---- COMPUTING A1, A2, ... ----'
         WRITE(IW,'('' KELG = '',I4,''     KELGX = '',I4)')KELG,KELGX
         WRITE(IW,'(''NOCC1 = '',I4,''     NVAC1 = '',I4)')NOCC1,NVAC1
         WRITE(IW,'(''NOC1X = '',I4,''     NVC1X = '',I4)')NOC1X,NVC1X
         NATCT = NATCUT(NOC1X+NVC1X+1,NOCC1+NOC1X+NVAC1+NVC1X)
         WRITE(IW,'(''NATCT = '',I4)')NATCT
      ENDIF
      ENDIF
      NFZATM(KELG) = NATCT
C-----------------------------------------------------------------------
C    A+A' REGION ALPHA PART: OCCUPIED AND VIRTUAL
C-----------------------------------------------------------------------
      JXX = 0
      IF(SCFTYP.EQ.ROHF) THEN
         CALL APRTOR(CXALPH,CXA,CAALPH,NOC1X,NOCC1,NUM,NAAO,NAMO,
     *        NAAOX,NOCA,JXX,NOC1XB)
      ELSE
         CALL APRTO(CXALPH,CXA,CAALPH,NOC1X,NOCC1,NUM,NAAO,NAMO,
     *        NAAOX,NOCA,JXX)
      ENDIF
      NATOC = NOCA
      CALL APRTV(CXALPH,CXA,CAALPH,NVAC1,NUM,NAAO,NAMO,
     *     NAAOX,NOCA,JXX,JXXX,NOC1X,NOCC)
C-----------------------------------------------------------------------
C     A+A' REGION BETA PART: OCCUPIED AND VIRTUAL
C-----------------------------------------------------------------------
      IF(SCFTYP.EQ.UHF) THEN
         JXX = 0
         CALL APRTO(CXBETA,CXB,CABETA,NOC1XB,NOCC1B,NUM,NAAO,NAMO,
     *              NAAOX,NOCA_B,JXX)
         CALL APRTV(CXBETA,CXB,CABETA,NVAC1B,NUM,NAAO,NAMO,
     *              NAAOX,NOCA_B,JXX,JXXX,NOC1XB,NOCCB)
      ENDIF
C-----------------------------------------------------------------------
      IF(MASWRK) THEN
      WRITE(IW,*)
      WRITE(IW,'(''FROZEN PART OF A: '',I4,''  OCC ALPHA ELEC'')')NOC1X
      WRITE(IW,'(''ACTIVE PART OF A: '',I4,''  OCC ALPHA ELEC'')')NOCC1
      WRITE(IW,'(''FROZEN PART OF A: '',I4,''  OCC  BETA ELEC'')')NOC1XB
      WRITE(IW,'(''ACTIVE PART OF A: '',I4,''  OCC  BETA ELEC'')')NOCC1B
      WRITE(IW,*)
      WRITE(IW,'(''FROZEN PART OF A: '',I4,'' VIRT ALPHA ELEC'')')NVC1X
      WRITE(IW,'(''ACTIVE PART OF A: '',I4,'' VIRT ALPHA ELEC'')')NVAC1
      WRITE(IW,'(''FROZEN PART OF A: '',I4,'' VIRT  BETA ELEC'')')NVC1XB
      WRITE(IW,'(''ACTIVE PART OF A: '',I4,'' VIRT  BETA ELEC'')')NVAC1B
      WRITE(IW,*)
      WRITE(IW,'('' DIMENSION NBF1 = '',I4)') NBF1
      WRITE(IW,*)
      WRITE(IW,'(''NOCC1, NOCC,NVAC1,NAMO='',4I4)')NOCC1,NOCC,NVAC1,NAMO
      WRITE(IW,'(''NOCC1B,NOCCB,NVAC1B   ='',3I4)')NOCC1B,NOCCB,NVAC1B
      ENDIF
C-----------------------------------------------------------------------
C
C     B REGION OCC
C
      JWW = 0
      DO 40 J=1+NOCC1,NOCC-NOC1X
         JWW=JWW+1
         DO 40 I=1,NAAO
            CWALPH(I,JWW)=CAALPH(I,J)
 40   CONTINUE
C
C     B REGION VIRT
C
      DO 50 J=1+NOCC-NOC1X+NVAC1,NAMO
         JWW=JWW+1
         DO 50 I=1,NAAO
            CWALPH(I,JWW)=CAALPH(I,J)
 50   CONTINUE
C
      IF(SCFTYP.EQ.UHF) THEN
         JWW_B =0
         DO 60 J=1+NOCC1B,NOCCB-NOC1XB
            JWW_B=JWW_B+1
            DO 60 I=1,NAAO
               CWBETA(I,JWW_B)=CABETA(I,J)
 60      CONTINUE
C
         DO 70 J=1+NOCCB-NOC1XB+NVAC1B,NAMO
            JWW_B=JWW_B+1
            DO 70 I=1,NAAO
               CWBETA(I,JWW_B)=CABETA(I,J)
 70      CONTINUE
      ENDIF
C
      IF(SCFTYP.EQ.UHF.AND.JWW.NE.JWW_B.AND.MASWRK) THEN
         WRITE(IW,*)'STH WRONG: JWW NOT EQUAL TO JWW_B IN B REGION'
         WRITE(6,*)'JWW, JWW_B = ',JWW,JWW_B
      ENDIF
C
      IF(IPRI.EQ.1.OR.IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(IW,*)'A-REGION LMO, JXX,NOCC1,NAAO=',JXX,NOCC1,NAAO
         CALL WRITN4(JXX,NAAO,CXALPH,NUM,NUM)
         WRITE(IW,*)'BEFORE MO REMOVED,B-REGION LMO'
         CALL WRITN4(JWW,NAAO,CWALPH,NUM,NUM)
      ENDIF
C-----------------------------------------------------------------------
C    REMOVE LINEAR DEPENDENT MOS DUE TO AOS REMOVAL
C
C    1. DIAGNOLIZE MATRIX {L^+(B+M) S(B+M) L(B+M)}
C    2. REMOVE VECTORS CORRESPONDS TO THE ZERO EIGENVALUES
C-----------------------------------------------------------------------
      IF(NHB.GT.0) THEN
C
C     MEMORY ADDRESS FOR MOS REMOVAL
C
         LH    = LAST1
         LVEC  = LH     + JWW*(JWW+1)/2
         LEIG  = LVEC   + JWW*JWW
         LTMP  = LEIG   + JWW
         LBIG  = LTMP   + (NAAO-NHB)*JWW
         LJBIG = LBIG   + JWW*8
         LAST2 = LJBIG  + JWW
C
C   REMOVE AOS BELONGING TO NTMLB TERMINAL ATOMS
C
       DO 80 J=1,JWW
       DO 80 I=1,NHB
                           CWALPH(NAAO-NHB+I,J)=0.0D+00
        IF(SCFTYP.EQ.UHF)  CWBETA(NAAO-NHB+I,J)=0.0D+00
 80    CONTINUE
C
       IF(SCFTYP.EQ.UHF) THEN
          CALL MORMV(NUM,NUM,NAAO-NHB,JWW,JWW2,NOCC2,CWALPH,SW,NHB,
     *         XX(LH),XX(LVEC),XX(LEIG),XX(LTMP),XX(LBIG),XX(LJBIG))
          CALL MORMV(NUM,NUM,NAAO-NHB,JWW,JWW2,NOCC2B,CWBETA,SW,NHB,
     *         XX(LH),XX(LVEC),XX(LEIG),XX(LTMP),XX(LBIG),XX(LJBIG))
       ELSE
          CALL MORMV(NUM,NUM,NAAO-NHB,JWW,JWW2,NOCC2,CWALPH,SW,NHB,
     *         XX(LH),XX(LVEC),XX(LEIG),XX(LTMP),XX(LBIG),XX(LJBIG))
          NOCC2B = NOCC2
       ENDIF
       JWW   = JWW2
C
       IF(IPRI.EQ.1.OR.IPRI.EQ.3.AND.MASWRK) THEN
          WRITE(IW,*)'AFTER MO REMOVED FOR TERMINAL ATOMS'
          WRITE(IW,*)'B-REGION LMO, JWW,NOCC2,NAAO=',JWW,NOCC2,NAAO
          CALL WRITN4(JWW,NAAO,CWALPH,NUM,NUM)
       ENDIF
C
      ENDIF
C-----------------------------------------------------------------------
C    AO-CUT
C-----------------------------------------------------------------------
       NCTAO = 0
       NOCUT = 0
       NOCUTB = 0
       IF(SCFTYP.EQ.UHF.OR.SCFTYP.EQ.ROHF) THEN
            CALL CTRUHF(SW,CWALPH,CXALPH,CWBETA,CXBETA,NUM,
     *          NFZOC,NFZVT,NFZOCB,NFZVTB,NCTAO,NCUT,NCUTX,
     *             NAAO,JWW,JXX,KELG,NCT,NOCC1,NOC1X,NOCUT,
     *                 NOCC1B,NOC1XB,NOCUTB,NHB,DOFF,DBOFF)
       ELSE
        CALL AOMOCUT(SW,CWALPH,CXALPH,NFZOC,NFZVT,NUM,NCTAO,NCUT,
     *   NCUTX,NAAO,JWW,JXX,KELG,NCT,NOCC1,NOC1X,NOCUT,DOFF,MASWRK)
        NOCC1B = NOCC1
        NOC1XB = NOC1X
        NOCUTB = NOCUT
      ENDIF
C-----------------------------------------------------------------------
      LDNEW = LAST1
      LAST0 = LDNEW + (NUM*NUM+NUM)/2
      IF(LCUT) THEN
         CALL SEQREW(NEG4)
         READ(NEG4)
         READ(NEG4) NCUTAT,NFRGBS,I,I,I,NUMX
         NTOT1 = NUM + NFRGBS
         NTOTAL = NTOT1*(NTOT1+1)/2
C
         LDTOT = LAST0
         LDBTOT = LDTOT
           IF(SCFTYP.NE.RHF) LDBTOT = LDTOT + NTOTAL
         LHTOT = LDBTOT + NTOTAL
         LDOLD = LHTOT  + NTOTAL
         LAST2 = LDOLD  + (NUM*NUM+NUM)/2
C
          IF(LAST2.GT.LAST) THEN
           IF(MASWRK) WRITE(6,*)' MEMORY ERROR IN ELGLOC'
           CALL ABRT
          ENDIF
C
         IF(SCFTYP.EQ.RHF) THEN
            CALL ENCUT1(NTOT1,NFRGBS,NUMX,NCUTAT,XX(LDTOT),XX(LHTOT))
         ELSE
            CALL ENCUT2(NTOT1,NFRGBS,NUMX,NCUTAT,XX(LDTOT),
     *                                        XX(LDBTOT),XX(LHTOT))
         ENDIF
         IF(MASWRK)
     *     WRITE(IW,*)'CPU TIME AFTER ENCUT'
         CALL TIMIT(1)
      ENDIF
C-----------------------------------------------------------------------
      NCTMO = NCTAO + NCTMOX
      NOCC1 = NOCC1 + NOC1X - NOCUT
      NOCC1B= NOCC1B + NOC1XB - NOCUTB
C
      LL2 = NUM*(NUM+1)/2
      IF(NCTMO.NE.0) THEN
         CALL DAREAD(IDAF,IODA,FAO,LL2,16,0)
         IF(LCUT) THEN
            CALL AOCUT(NCTAO,NEG4,XX(LDTOT),NTOTAL,XX(LDBTOT))
            CALL VCLR(HAO,1,LL2)
            CALL DENRLX(NEG5,HAO,XX(LDOLD),LL2,NUMX,NFRGBS,NCTAO)
            IF(SCFTYP.NE.RHF) CALL DAREAD(IDAF,IODA,XX(LDNEW),LL2,20,0)
         ELSE
            IF(SCFTYP.NE.RHF) CALL DAREAD(IDAF,IODA,XX(LDNEW),LL2,20,0)
            CALL AOCUT(NCTAO,NEG4,FAO,LL2,XX(LDNEW))
         ENDIF
         CALL CUTSHT(FAO,LL2,NUM,NCTAO,NHB)
         CALL CUTSHT(XX(LDNEW),LL2,NUM,NCTAO,NHB)
C
      ELSE
         LL2 = NUM*(NUM+1)/2
         CALL DAREAD(IDAF,IODA,FAO,LL2,16,0)
      ENDIF
C
C  DETERMINE SHELL INFOMATION FOR TERMINAL ATOMS
C
      NSHXXX = 0
      IF(DIRSCF.AND.LCUT) THEN
         DO 110 I=1,NSHELL
            IF(KATOM(I).EQ.NAT) NSHXXX= I
 110     CONTINUE
         NSHELL = NSHXXX
      ENDIF
C
      CALL SEQREW(NEG1)
      CALL SEQREW(NEG5)
      WRITE(NEG1) JXX,JWW,NAAO,NOCC1,NOCC2,NOCC1B,NOCC2B,MUL,NHB,NAT
      WRITE(NEG1) NATOC,NCTMO,NCUT,NCTAO,KELG,NSHELL,NHTSHL,
     *     NFLTRM,NRCTRM,NPSTRM
      WRITE(NEG1) (NFZOC(I),I=1,KELG),(NFZVT(I),I=1,KELG)
      WRITE(NEG1) (NFZOCB(I),I=1,KELG),(NFZVTB(I),I=1,KELG)
      WRITE(NEG1) (NFZATM(I),I=1,KELG)
C
      CALL VALFM(LOADFM)
      LZMAT = LOADFM + 1
      LZWRK = LZMAT  + NAT*3
      LC    = LZWRK  + NAT
      LAST  = LC     + NAT*3
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL ELGZMT(NAT,NEG2,NTMLB,XX(LZMAT),XX(LZWRK))
      IF(NCUT.GT.0) THEN
         NATCT=NFZATM(NCUT)
      ELSE
         NATCT=0
      ENDIF
      CALL ELGCRT(NEG2,NAT,NTMLB,NATCT,XX(LC))
C
      CALL RETFM(NEED)
C
C     NEW A-REGION
C
                        WRITE(NEG5)((CXALPH(ID,JD),ID=1,NAAO),JD=1,JXX)
      IF(SCFTYP.EQ.UHF) WRITE(NEG5)((CXBETA(ID,JD),ID=1,NAAO),JD=1,JXX)
C
C     NEW B-REGION
C
                        WRITE(NEG5)((CWALPH(ID,JD),ID=1,NAAO),JD=1,JWW)
      IF(SCFTYP.EQ.UHF) WRITE(NEG5)((CWBETA(ID,JD),ID=1,NAAO),JD=1,JWW)
C
      IF(NCTMO.NE.0) THEN
        NNXX = (NUM-NHB-NCTAO)*(NUM-NHB-NCTAO+1)/2
        NNYY = (NUM-NCTAO)*(NUM-NCTAO+1)/2
                          WRITE(NEG5) (FAO(I),I=1,NNXX)
        IF(SCFTYP.NE.RHF) WRITE(NEG5) (XX(LDNEW+I-1),I=1,NNXX)
         IF(LCUT) CALL VADD(DOFF,1,HAO,1,DOFF,1,LL2)
                          WRITE(NEG5) (DOFF(I),I=1,NNYY)
         IF(SCFTYP.NE.RHF) THEN
            IF(LCUT) CALL VADD(DBOFF,1,XX(LDOLD),1,DBOFF,1,LL2)
                          WRITE(NEG5) (DBOFF(I),I=1,NNYY)
         ENDIF
      ELSE
         WRITE(NEG5) (FAO(I),I=1,(NUM-NHB)*(NUM+1-NHB)/2)
         IF(SCFTYP.NE.RHF) THEN
            CALL DAREAD(IDAF,IODA,FAO,LL2,20,0)
            WRITE(NEG5) (FAO(I),I=1,(NUM-NHB)*(NUM+1-NHB)/2)
         ENDIF
      ENDIF
CGU=========================================================
CGU   IF(IGFLD.EQ.0.AND.NELONG.EQ.2) THEN
CGU    CALL SEQREW(80)
CGU    WRITE(80) NOCC1,NVAC1,JXX,NOCC2,NVAC2,JWW
CGU    WRITE(80) ((CXALPH(ID,JD),ID=1,NAAO),JD=1,JXX)
CGU    WRITE(80) ((CWALPH(ID,JD),ID=1,NAAO),JD=1,JWW)
CGU   ENDIF
CGU=========================================================
CGUC
CGUC  FOR ELG-TDHF, LMOS FOR PERTURBED C'S
CGUC      SURELY THIS SHOULD BE RUNTYP
CGUC
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU       CALL CNLMO
CGU    ENDIF
C===========================================================
C
C   WRITE-OUT LMOS FOR MOLCAD PLOT
C
      LMO = 0
      IF(LMO.NE.0) THEN
         IST1 = 1
         IST2 = 1
         IF(NOCC1.GT.5) IST1 = NOCC1 - 5 + 1
         IF(NOCC2.GT.5) IST2 = NOCC2 - 5 + 1
C
      IF(MASWRK) THEN
         WRITE(6,*)
         WRITE(6,*)'LMOS OVER AOS FOR MOLCAD PLOT'
         WRITE(6,*)
         DO 115 ID=1,NAAO
            WRITE(IW,116)ID,(CXALPH(ID,JD),JD=IST1,NOCC1),
     *           (CWALPH(ID,JD),JD=IST2,NOCC2)
 115     CONTINUE
         WRITE(IW,*)
       ENDIF
      ENDIF
C      IF(MASWRK) WRITE(IW,*)'CPU TIME AFTER FINAL WRITTING'
C      CALL TIMIT(1)
C
      RETURN
 116  FORMAT(I4,2X,10(F8.4,1X))
      END
C
C*MODULE ELGLOC  *DECK EGSAVE
      SUBROUTINE EGSAVE(CA,CB,T,N1,N, NOCC1,NVAC1,NOCC2,NVAC2,
     *       NOCC1B,NVAC1B,NOCC2B,NVAC2B, NOC1X,NVC1X,NOC1XB,NVC1XB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      DATA UHF/8HUHF     /
C
      DIMENSION CA(N1,N),CB(N1,N),T(N1,N1)
C
      NOCC = NOCC1 + NOCC2
      NOCCB= NOCC1B+ NOCC2B
      L3 = N1*N1
C
C      ----- READ VEC FROM DAF -----
C
       NOCOLD=NOC1X
       CALL DAREAD(IDAF,IODA,T,L3,15,0)
       IF(N1.EQ.N) THEN
          DO 10 I=1,N1
             DO 10 J=1,N
                CA(I,J) = T(I,J)
 10       CONTINUE
       ELSE
          DO 30 I=1,N1
             DO 20 J=1,NOCC1+NOCC2
                CA(I,J) = T(I,NOCOLD+J)
                JJ = J
 20          CONTINUE
             DO 30 J=1,NVAC1+NVAC2
                CA(I,JJ+J) = T(I,NOCC+NOC1X+NVC1X+J)
 30        CONTINUE
        ENDIF
C
      IF(SCFTYP.EQ.UHF) THEN
         CALL DAREAD(IDAF,IODA,T,L3,19,0)
         IF(N1.EQ.N) THEN
            DO 40 I=1,N1
               DO 40 J=1,N
                  CB(I,J) = T(I,J)
 40         CONTINUE
         ELSE
            DO 60 I=1,N1
               DO 50 J=1,NOCC1B+NOCC2B
                  CB(I,J) = T(I,NOC1XB+J)
                  JJ = J
 50            CONTINUE
               DO 60 J=1,NVAC1B+NVAC2B
                  CB(I,JJ+J) = T(I,NOCCB+NOC1XB+NVC1XB+J)
 60         CONTINUE
         ENDIF
      ENDIF
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK REGA
      SUBROUTINE REGA(NATM,NBF1,NTMLB,NHB,NOCC1A,NVAC1A,NOCC2A,NVAC2A,
     *                NOCC1B,NVAC1B,NOCC2B,NVAC2B,NASPIN,NBNDAB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL ATOB
      PARAMETER (MXATM=2000, MXAO=8192)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /CORE/ CORE(107)
      DATA RHF/8HRHF     /
      DATA NONE/4HNONE/
      COMMON /ELGLAB/ IATM(MXAO)
      COMMON /CTSHFT/ ATOB
C-----------------------------------------------------------------------
C     HERE WE DETERMINE THE NO. OF BASIS-SET FUNCTIONS IN REGION A:  NBF
C     ADDITIONAL ASSUMPTIONS:
C      1) BASIS-SET DOES NOT INCLUDE HIGHER ANGULAR FUNCS (THE LAST F)
C         OTHERWISE ERROR
C      2) FIRST NATM ATOMS FROM THE INPUT $DATA BLOCK BELONG TO PART A
C      3) TERMINAL ATOM IS THE LAST IN THE INPUT $DATA BLOCK
C-----------------------------------------------------------------------
      NBF1 = 0
      NHB = 0
      DO 10 I=1,NUM
         IDB = IATM(I)
         IF(IDB.LE.NATM) NBF1 = NBF1 + 1
         IF((IDB+NTMLB).GT.NAT) NHB = NHB + 1
 10   CONTINUE
C---------------------------------------------------------------------
      NAELC = 0
      DO 20 I=1,NATM
         NAELC = NAELC + INT(ZAN(I))
         IF(MPCTYP.NE.NONE) THEN
            NUCZ = INT(ZAN(I))
            MPCORE = NUCZ - INT(CORE(NUCZ))
            NAELC = NAELC - MPCORE
          ENDIF
 20    CONTINUE
C
       IF(SCFTYP.EQ.RHF) THEN
C
C  ODD ELECTRON TO A OR B REGION TO FORM CLOSE-SHELL,
C  IF NAELC IS NOT EVEN
C
          NAELC = NAELC - ICH
          IF(NBNDAB.NE.0) THEN
            IF(ATOB) THEN
              NAELC = NAELC - NBNDAB
            ELSE
              NAELC = NAELC + NBNDAB
            END IF
          ENDIF
          NADTL = MOD(NAELC,2)
C
          IF(ATOB) THEN
             NAELC = NAELC - NADTL
          ELSE
             NAELC = NAELC + NADTL
          ENDIF
C
          NOCC  = NE/2
          NOCC1A = NAELC/2
          NOCC1B = NOCC1A
          NOCC2A = NOCC - NOCC1A
          NOCC2B = NOCC2A
          NVAC1A = NBF1 - NOCC1A
          NVAC1B = NVAC1A
          NVAC2A = NUM - (NOCC + NVAC1A)
          NVAC2B = NVAC2A
       ELSE
C
 101      NOCC1A = (NAELC + NASPIN)/2
          NOCC1B = NOCC1A - NASPIN
           IF((NOCC1A+NOCC1B).NE.NAELC) THEN
             NAELC = NAELC - 1
             IF(NAELC.LE.0) CALL ABRT
             GO TO 101
           ENDIF
C
          NBELC  = NE - NAELC
          NBSPIN = MUL - NASPIN - 1
          NOCC2A = (NBELC + NBSPIN)/2
          NOCC2B = NOCC2A - NBSPIN
C
         NVAC1A = NBF1 - NOCC1A
         NVAC1B = NBF1 - NOCC1B
         NVAC2A = NUM - NVAC1A - NOCC1A - NOCC2A
         NVAC2B = NUM - NVAC1B - NOCC1B - NOCC2B
C
         IF( (NOCC2A+NOCC1A-NOCC2B-NOCC1B+1).NE.MUL) CALL ABRT
      ENDIF
C
      RETURN
      END
C*MODULE ELGLOC  *DECK APRTO
      SUBROUTINE APRTO(CXALPH,CXA,CAALPH,NOC1X,NOCC1,NUM,NAAO,
     *                 NAMO,NWAX,NOCA,JXX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION CXALPH(NUM,NUM),CXA(NUM,NUM),CAALPH(NAAO,NAMO)
C
      DO 10 J=1,NOC1X
         JXX=JXX+1
         DO 10 I=1,NWAX
            CXALPH(I,JXX)=CXA(I,J)
 10   CONTINUE
C
      DO 20 J=1,NOCC1
         JXX=JXX+1
         DO 20 I=1,NAAO
            CXALPH(I,JXX)=CAALPH(I,J)
 20   CONTINUE
C
      NOCA = JXX
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK APRTV
      SUBROUTINE APRTV(CXALPH,CXA,CAALPH,NVAC1,NUM,NAAO,
     *                 NAMO,NWAX,NOCA,JXX,JXXX,NOC1X,NOCC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION CXALPH(NUM,NUM),CXA(NUM,NUM),CAALPH(NAAO,NAMO)
C
      DO 10 J=1+NOC1X,JXXX
         JXX=JXX+1
         DO 10 I=1,NWAX
            CXALPH(I,JXX)=CXA(I,J)
 10   CONTINUE
C
      DO 20 J=1,NVAC1
         JXX=JXX+1
         DO 20 I=1,NAAO
            CXALPH(I,JXX)=CAALPH(I,NOCC-NOC1X+J)
 20   CONTINUE
C
      NOCA = JXX
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK APRTOR
      SUBROUTINE APRTOR(CXALPH,CXA,CAALPH,NOC1X,NOCC1,NUM,NAAO,
     *                  NAMO,NWAX,NOCA,JXX,NOC1XB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION CXALPH(NUM,NUM),CXA(NUM,NUM),CAALPH(NAAO,NAMO)
C
      DO 10 J=1,NOC1XB
         JXX=JXX+1
         DO 10 I=1,NWAX
            CXALPH(I,JXX)=CXA(I,J)
 10   CONTINUE
C
      DO 20 J=1,NOCC1
         JXX=JXX+1
         DO 20 I=1,NAAO
            CXALPH(I,JXX)=CAALPH(I,J)
 20   CONTINUE
C
      NSOCC=NOC1X-NOC1XB
      DO 30 J=1,NSOCC
         JXX=JXX+1
         DO 30 I=1,NWAX
            CXALPH(I,JXX)=CXA(I,J+NOC1XB)
 30   CONTINUE
C
      NOCA = JXX
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK RAOLMO
      SUBROUTINE RAOLMO(D,DA,DB,C,T,N,N1,NOCC,NOCC1,NVAC1,NOC1X,
     *                  NVC1X,NOCC2,L2,VEC,EIG,X,
     *                  VS,VV,Y,ES,IEND,NAME,WRK,IWRK,NWRK)
C----------------------------------------------------------------------
C   INPUT:  D: AO-BASIS DENSITY-MAT
C           S: AO-BASIS OVERLAP-MAT
C           C: AO-BASIS CMOS
C
C  OUTPUT:  C: AO-BASIS RLMOS
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL SUCC,NEEDA
      CHARACTER*3 NAME
C
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
C
      DIMENSION D(L2),DA(L2),DB(L2)
      DIMENSION T(N1,N1),VV(N1,N1),X(N1,N1),Y(N1,N1)
      DIMENSION VEC(N1,N1),EIG(N1),VS(N1,N1),ES(N1),C(N1,N)
      DIMENSION WRK(NWRK),IWRK(NWRK)
C
      DATA NONE/4HNONE/
C----------------------------------------------------------------
      ZERO = 0.0D+00
      ONE  = 1.0D+00
      NVAC = N1 - NOCC
      NA = IEND
      NB = N1 - NA
C----------------------------------------------------------------
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)' D-MAT AT THE BEGINNING OF RAOLMO'
         CALL PRTRIL(D,N1)
      ENDIF
C
      CALL CPYTSQ(D,VS,N1,1)
C
C     FROM AO TO OAO, FOR AB INITIO (NON-ORTHOGONAL BASIS SET)
C
      IF(MPCTYP.EQ.NONE) THEN
        CALL DGEMM('T','N',N1,N1,N1,1.0D+00,X,N1,VS,N1,0.0D+00,VEC,N1)
        CALL DGEMM('N','N',N1,N1,N1,1.0D+00,VEC,N1,X,N1,0.0D+00,VS,N1)
      ENDIF
C---------------------------------------------------------------------
      IJ = 0
      IK = 0
      DO 40 I=1,N1
         DO 40 J=1,I
            IJ = IJ + 1
            IF(I.LE.NA.AND.J.LE.NA) THEN
               DA(IJ) = VS(I,J)
            ENDIF
            IF(I.GT.NA.AND.J.GT.NA) THEN
               IK = IK + 1
               DB(IK) = VS(I,J)
            ENDIF
 40   CONTINUE
      DO I=1, N1
        DO J=1,N1
          T(I,J) = 0.0D+00
        END DO
      END DO
C----------------------------------------------------------------
C
      CALL JACDG(DA,VEC,EIG,IWRK,WRK,N1,NA)
      CALL GJACORD(VEC,EIG,NA,N1)
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)'EIGENVALUES/VECTORS OF DA-MATRIX'
         CALL PREVNL(VEC,EIG,NA,NA,N1)
      ENDIF
C
C      SURELY THIS SHOULD BE RUNTYP
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU      WRITE(NEG9)((X(I,J),I=1,N1),J=1,N1)
CGU      WRITE(NEG9)((Y(I,J),I=1,N1),J=1,N1)
CGU      WRITE(NEG9)(EIG(J),J=1,NA)
CGU      WRITE(NEG9)((VEC(I,J),I=1,NA),J=1,NA)
CGU    ENDIF
C
      NDIF = NA - (NOCC1 + NVAC1 + NOC1X + NVC1X)
C
       DO 50 J=1,NOCC1
          ES(J) = EIG(J)
          DO 50 I=1,NA
             T(I,J) = VEC(I,J)
 50    CONTINUE
C
       JJ = NOCC - NOC1X
       DO 60 J=1,NVAC1
          JJ = JJ + 1
          ES(JJ) = EIG(NOCC1+NDIF+J)
          DO 60 I=1,NA
             T(I,JJ) = VEC(I,NOCC1+NDIF+J)
 60    CONTINUE
C
C  MOVE THE NOT-GOOD VAC MO TO THE END OF THE B-REGION
C
       IF(NDIF.GE.1) THEN
          DO 70 J=1,NDIF
             ES(NA-J+1) = EIG(NOCC1+J)
             DO 70 I=1,NA
                T(I,NA-J+1) = VEC(I,NOCC1+J)
 70    CONTINUE
      ENDIF
C
      NST = NOCC1+NVAC1+NDIF+1
      JJ = 0
       DO 80 J=NST,NA
          JJ = JJ + 1
          ES(N1-NDIF-JJ+1) = EIG(J)
          DO 80 I=1,NA
             T(I,N1-NDIF-JJ+1) = VEC(I,J)
 80   CONTINUE
C
      CALL JACDG(DB,VEC,EIG,IWRK,WRK,N1,NB)
      CALL GJACORD(VEC,EIG,NB,N1)
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)'EIGENVALUES/VECTORS OF DB-MATRIX'
         CALL PREVNL(VEC,EIG,NB,NB,N1)
      ENDIF
C
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU     WRITE(NEG9)(EIG(J),J=1,NB)
CGU     WRITE(NEG9)((VEC(I,J),I=1,NB),J=1,NB)
CGU    ENDIF
C
       DO 90 J=1,NB
          JJ = NOCC1 + J
          IF(J.GT.NOCC2) JJ = (NOCC1+NVAC1) + J
          ES(JJ) = EIG(J)
          DO 90 I=1,NB
             T(NA+I,JJ) = VEC(I,J)
 90    CONTINUE
C
       IF(IPRI.EQ.3.AND.MASWRK) THEN
          WRITE(6,*)'T-MAT'
          CALL WRITN4(N1,N1,T,N1,N1)
       ENDIF
C
C     CHECK TO SEE IF THERE EXISTS DEGENERACY BETWEEN OCC AND VAC
C
C     NBD: POSSIBLE BAD EIGS
C
       NBD = 0
       DO 100 I=1,NOCC-NOC1X
          IF(ES(I).LT.1.0D+00) THEN
             NBD = NBD + 1
             IF(MASWRK)
     *   WRITE(6,'(''      POSSIBLE BAD EIGS:'',I4,2X,F12.5)')I,ES(I)
          ENDIF
 100   CONTINUE
       IF(NBD.NE.0.AND.MASWRK)
     *        WRITE(6,'(''TOTAL POSSIBLE BAD EIGS:'',I4)')NBD
C
      IF(NBD.GT.0.AND.MASWRK)
     * WRITE(6,'(''   IT '',''  JO     OCCMIN'',6X,''JV      VIRMAX'')')
       DO 110 IT=1,NBD
          CALL BDCHK(ES,T,VV,VEC,N1,NOCC-NOC1X,IT)
 110   CONTINUE
C
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU     WRITE(NEG9)((T(I,J),I=1,N1),J=1,N1)
CGU    ENDIF
C
C  CALC NEW RAO-BASIS D-MAT
C
      CALL DGEMM('T','N',N1,N1,N1,ONE,T,N1,VS,N1,ZERO,VEC,N1)
      CALL DGEMM('N','N',N1,N1,N1,ONE,VEC,N1,T,N1,ZERO,VS,N1)
      CALL VCLR(VEC,1,N1*N1)
C---------------------------------------------------------------------
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)'RO-BASIS D-MAT'
         CALL WRITN4(N1,N1,VS,N1,N1)
      END IF
C
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU     WRITE(NEG9)((VS(I,J),I=1,N1),J=1,N1)
CGU    ENDIF
C---------------------------------------------------------------------
C   PERFORM BLOCK DIAGONALIZATION TO LOCALIZE RAO-BASIS MOS
C---------------------------------------------------------------------
       NEEDA = .FALSE.
       CALL BDJAC(VS,N1,N1,NOCC-NOC1X,EIG,VEC, NORT,SUCC,SM,NEEDA)
        IF(.NOT.SUCC) THEN
         IF(MASWRK) WRITE(6,*)'BLOCK DIAG NOT SUCCESSFUL'
         CALL ABRT
        ENDIF
C
CGU    IF(SCFTYP.EQ.TDHF) THEN
CGU     WRITE(NEG9) (EIG(I),I=1,N1)
CGU     WRITE(NEG9)((VEC(I,J),I=1,N1),J=1,N1)
CGU    ENDIF
C
      IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)'DIAGONALIZED D-MAT BY RLMO?'
         CALL WRITN4(N1,N1,VS,N1,N1)
         WRITE(6,*)'EIGENVECTOR OF RLMO-D-MAT'
         CALL WRITN4(N1,N1,VEC,N1,N1)
      ENDIF
C
      IF(MPCTYP.EQ.NONE) THEN
        CALL DGEMM('N','N',N1,N1,N1,ONE,Y,N1,T,N1,ZERO,VS,N1)
      ELSE
        CALL DCOPY(N1*N1,T,1,VS,1)
      ENDIF
C
      IF(NAME.EQ.'OCC') THEN
        CALL DGEMM('N','N',N1,NOCC-NOC1X,N1,ONE,VS,N1,VEC,N1,ZERO,C,N1)
        CALL DAO(DA,C,NOCC,NOC1X,N1,N)
      ELSE
        CALL DGEMM('N','N',N1,NOCC-NOC1X,N1,ONE,VS,N1,VEC,N1,ZERO,T,N1)
          DO J=1,NOCC-NOC1X
          DO I=1,N1
           C(I,NVAC-NVC1X+J) = T(I,J)
          ENDDO
          ENDDO
        CALL DAOVAC(DA,C,NVAC,NVC1X,N1,N)
      ENDIF
C
       IF(IPRI.EQ.3.AND.MASWRK) THEN
          WRITE(6,*)' D-MAT CALC FROM AO-BASIS RLMOS'
          CALL PRTRIL(DA,N1)
       ENDIF
C
       IJ = 0
       DIFF = 0.0D+00
       DMAX = 0.0D+00
       DO 190 I=1,N1
       DO 190 J=1,I
       IJ = IJ + 1
          DIFPRT=ABS(DA(IJ) - D(IJ))
          DIFF = DIFF + DIFPRT
          IF(DIFPRT.GT.DMAX) THEN
            DMAX=DIFPRT
            IX = I
            JX = J
          ENDIF
 190   CONTINUE
C
       DIFF = DIFF/(N1*N1)
C
      IF(MASWRK) THEN
       WRITE(6,'(''MAX DIFFERENCE IN D-MATS: '',E12.5,'' ,FOR I,J='',
     *          I4,1X,I4,/,''      AVERAGE DIFFERENCE: '',E12.5,/)')
     *           DMAX,IX,JX,DIFF
      ENDIF
C
       IF(DIFF.GT.1.0D-04) THEN
       IF(MASWRK) WRITE(6,*)'LOCALIZATION IS NOT UNITARY TRANSFORMATION'
          CALL ABRT
       ENDIF
C-------------------------------------------------------------------
      RETURN
      END
C
C*MODULE ELGLOC  *DECK GJACORD
      SUBROUTINE GJACORD(VEC,EIG,N,LDVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VEC(LDVEC,N),EIG(N)
C
C     ---- SORT EIGENDATA INTO DECENDING ORDER -----
C
      DO 290 I = 1, N
         JJ = I
         DO 270 J = I, N
            IF (EIG(J) .GT. EIG(JJ)) JJ = J
 270     CONTINUE
         IF (JJ .EQ. I) GO TO 290
         T = EIG(JJ)
         EIG(JJ) = EIG(I)
         EIG(I) = T
         DO 280 J = 1, N
            T = VEC(J,JJ)
            VEC(J,JJ) = VEC(J,I)
            VEC(J,I) = T
 280     CONTINUE
 290  CONTINUE
      RETURN
      END
C
C*MODULE ELGLOC  *DECK BDCHK
      SUBROUTINE BDCHK(ES,T,VV,VEC,N1,NOCC,IT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION EIG(2),VS(2,2),DA(3),JBIG(2),BIG(2)
      DIMENSION T(N1,N1),ES(N1),VV(N1,N1),VEC(N1,N1)
C
C  CHECK IF OCCUPANCIES ARE IN ORDER:
C
      NVIRST = NOCC + 1
      VIRMAX = 0.0D+00
      DO 60 J = NVIRST,N1
         IF(ES(J).LT.VIRMAX) GO TO 60
        VIRMAX = ES(J)
        JV = J
 60   CONTINUE
      OCCMIN = 100.0D+00
      DO 70 I = 1,NOCC
         IF(ES(I).GT.OCCMIN) GO TO 70
         OCCMIN = ES(I)
         JO = I
 70   CONTINUE
      X = OCCMIN - VIRMAX
C
      IF(X.GT.1.0D-03) THEN
         IF(MASWRK) WRITE(6,'(''NO BAD FOUND FOR IT='',I4)')IT
         RETURN
      ENDIF
C
C   IF OCCMIN < VIRMAX, SIMPLY EXCHANGE JO AND JV TWO T'S COLUMNS.
C
      IF(ABS(X).GT.1.0D-05) THEN
      IF(MASWRK) WRITE(6,'(I4,2X,I4,F12.5,4X,I4,F12.5,
     *  ''  EXCHANGE THESE TWO BAD T COLUMNS'')')IT,JO,OCCMIN,JV,VIRMAX
         DO 80 I=1,N1
            TX = T(I,JO)
            T(I,JO) = T(I,JV)
            T(I,JV) = TX
 80      CONTINUE
         TX = ES(JO)
         ES(JO) = ES(JV)
         ES(JV) = TX
         RETURN
      ENDIF
C
      DO 90 I=1,N1
         VV(I,1) = T(I,JO)
         VV(I,2) = T(I,JV)
 90   CONTINUE
C
C   CALC D-MAT IN RAO BASIS: D=T^+ D T USING THESE TWO DEGENERATE VECS
C
      DO 100 I=1,2
         DO 100 J=1,2
            VS(I,J) = 0.0D+00
            DO 100 K=1,N1
               DO 100 L=1,N1
                  VS(I,J) = VS(I,J) + VV(K,I)*VEC(K,L)*VV(L,J)
 100  CONTINUE
C
      DA(1) = VS(1,1)
      DA(2) = VS(2,1)
      DA(3) = VS(2,2)
C
      CALL JACDG(DA,VS,EIG,JBIG,BIG,2,2)
      IF(EIG(1).GT.EIG(2)) THEN
         ES(JO) = EIG(1)
         ES(JV) = EIG(2)
      ELSE
         ES(JO) = EIG(2)
         ES(JV) = EIG(1)
      ENDIF
C
C     NEW T(AO,RAO)
C
      DO 110 I=1,N1
         DO 110 J=1,2
            VEC(I,J) = 0.0D+00
            DO 110 K=1,2
               VEC(I,J) = VEC(I,J) + VV(I,K)*VS(K,J)
 110  CONTINUE
C
      DO 120 I=1,N1
         IF(EIG(1).GT.EIG(2)) THEN
            T(I,JO) = VEC(I,1)
            T(I,JV) = VEC(I,2)
         ELSE
            T(I,JO) = VEC(I,2)
            T(I,JV) = VEC(I,1)
         ENDIF
 120  CONTINUE
C
      IF(MASWRK) THEN
      WRITE(6,*)'T(AO,RAO) AND EIGS AFTER REMOVING DEGENERACY'
      CALL PREVNL(T,ES,N1,N1,N1)
      END IF
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK MORMV
      SUBROUTINE MORMV(NUM,MUM,NWA,NWMM,NWM,NOCC,CWALPH,SW,NHB,
     *          H,VEC,EIG,TMP,WRK,IWRK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION CWALPH(NUM,MUM),SW(NUM,NUM)
      DIMENSION H(NWMM*(NWMM+1)/2),EIG(NWMM)
      DIMENSION VEC(NWMM,NWMM),TMP(NWA,NWMM)
      DIMENSION WRK(NWMM,8),IWRK(NWMM)
C----------------------------------------------------------------------
C
      CALL DGEMM('N','N',NWA,NWMM,NWA,1.0D+00,SW,NUM,CWALPH,NUM,
     *           0.0D+00,TMP,NWA)
      CALL DGEMM('T','N',NWMM,NWMM,NWA,1.0D+00,CWALPH,NUM,TMP,NWA,
     *           0.0D+00,VEC,NWMM)
      CALL LINEAR(VEC,H,NWMM,NWMM)
C
      DO I=1,NWMM
        VEC(I,I) = 1.0D+00
      DO J=1,NWMM
        IF(I.NE.J) VEC(I,J) = 0.0D+00
      ENDDO
      ENDDO
C
      NB1 = NWMM
      NB2 = (NB1*NB1+NB1)/2
      NMIN = 1
      NMAX = NB1
C
      CALL JACDIA(H,VEC,NB1,NB2,NWMM,NMIN,NMAX,WRK,IWRK)
C
C      CALL GLDIAG(NWMM,NWMM,NWMM,H,WRK,EIG,VEC,IERR,IWRK)
C
      CALL DGEMM('N','N',NWA,NWMM,NWMM,1.0D+00,CWALPH,NUM,VEC,NWMM,
     *           0.0D+00,TMP,NWA)
C
C       WRITE(6,*)'H REMOVAL EIGEN VECTOR AND VALUES,NWMM=',NWMM
      DO I=1,NWMM
        EIG(I) = H((I*I+I)/2)
      ENDDO
C
      II = 0
      IMV1 = 0
      IMV2 = 0
      NOCCN = 0
      NVIT = 0
      DO 50 I=1,NWMM
         IF(EIG(I).LE.1.0D-01) THEN
            IF(I.LE.NOCC) IMV1 = IMV1 + 1
            IF(I.GT.NOCC) IMV2 = IMV2 + 1
         ELSE
            II = II + 1
            IF(I.LE.NOCC) NOCCN = NOCCN + 1
            IF(I.GT.NOCC) NVIT = NVIT + 1
            DO 40 K=1,NWA
               CWALPH(K,II) = TMP(K,I)
 40         CONTINUE
         ENDIF
 50   CONTINUE
C
      NWM = II
      NOCC = NOCCN
C
C  NORMALIZE CW
C
      DO I=1,NWM
         FCOR = 0.0D+00
         DO K=1,NWA
            DO L=1,NWA
               FCOR=FCOR + CWALPH(K,I)*SW(K,L)*CWALPH(L,I)
            ENDDO
         ENDDO
C
         DO K=1,NWA
            CWALPH(K,I) = CWALPH(K,I)/SQRT(FCOR)
         ENDDO
      ENDDO
C
      IF(MASWRK) THEN
      WRITE(6,'(''MOS REMOVED IN OCC:'',I4, ''  IN VIR:'',I4)')IMV1,IMV2
      WRITE(6,'(''REMOVED TERMINAL, NOCC='',I4,''  NWM='',I4)')NOCC,NWM
      END IF
      NMVT = IMV1 + IMV2
      IF(NMVT.NE.NHB) THEN
         IF(MASWRK)
     *     WRITE(6,*)'REMOVED ORBITALS',NMVT,' NOT EQUAL TO NHB',NHB
         STOP
      ENDIF
C---------------------------------------------------------------
      RETURN
      END
C
C*MODULE ELGLOC  *DECK AOMOCUT
      SUBROUTINE AOMOCUT(SW,CWALPH,CXALPH,NFZOC,NFZVT,NUM,NCTAO,NCUT,
     *  NCUTX,NAAO,JWW,JXX,NCYCLE,NCT,NOCA,NOCXA,NOCUT,DOFF,MASWRK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL MASWRK
      DIMENSION CWALPH(NUM,NUM),CXALPH(NUM,NUM),SW(NUM,NUM)
      DIMENSION DOFF(NUM*(NUM+1)/2)
      PARAMETER (MXSTP=500)
      PARAMETER (CTHRSH=1.0D-05)
      DIMENSION NFZOC(MXSTP),NFZVT(MXSTP)
C
      IF(MASWRK) THEN
       WRITE(6,'(''PREPARING FOR AO-CUT,NCUTX='',I4)')NCUTX
       WRITE(6,'('' NSTEP   NFZOC  NFZVT  '')')
        DO I=1,NCYCLE
         WRITE(6,'(3I6,2X)')I,NFZOC(I),NFZVT(I)
        ENDDO
      ENDIF
C
      IF(NCYCLE.LE.1) RETURN
C
      NCUT = NCUTX
      NCTAO = NFZOC(NCUT+1) + NFZVT(NCUT+1)
      FCUT = 0.0D+00
      DO 20 I=1,JWW
         DO 20 K=1,NCTAO
            DO 20 L=1,NCTAO
               FCUT = FCUT + ABS(CWALPH(K,I)*SW(K,L)*CWALPH(L,I))
 20   CONTINUE
       IF(MASWRK) WRITE(6,'(''FCUT,CTHRSH='',2(E10.3,2X))')FCUT,CTHRSH
C
      IF(FCUT.LT.CTHRSH.AND.NCT.NE.0) THEN
      IF(MASWRK)WRITE(6,*)'AO-CUT CAN BE PERFORMED IN THE NEXT ELG STEP'
C
C     PERFORM AO-CUT IN B-REGION
C
         NCUT = NCUT + 1
         NAAO  = NAAO  - NCTAO
         NOCUT = NFZOC(NCUT)
         NVCUT = NFZVT(NCUT)
         NOCRST = NOCA + NOCXA - NOCUT
         DO 30 I=1,JWW
            DO 30 J=1,NAAO
               CWALPH(J,I) = CWALPH(NCTAO+J,I)
 30      CONTINUE
CGU
C    CALCULATE CUT-OFF DENSITY MATRIX
CGU
         II = 0
         DO 40 I=1,NUM
            DO 40 J=1,I
               II = II + 1
               DOFF(II) = 0.0D+00
               DO 40 K=1,NOCUT
                  DOFF(II) = DOFF(II) + CXALPH(I,K)*CXALPH(J,K)*2.0D+00
 40      CONTINUE
C
         II = 0
         JJ = 0
         DO 50 I=1,NUM
            DO 50 J=1,I
               II = II  + 1
               IF(J.GT.NCTAO.AND.I.GT.NCTAO) THEN
                  JJ = JJ + 1
                  DOFF(JJ) = DOFF(II)
               ENDIF
 50      CONTINUE
C
C   PERFORM MO-CUT IN A-REGION
C
         JXX = JXX - NCTAO
         DO 70 J=1,NAAO
            DO 60 I=1,NOCRST
               CXALPH(J,I) = CXALPH(NCTAO+J,NOCUT+I)
 60         CONTINUE
            DO 70 I=NOCRST+1,JXX
               CXALPH(J,I) = CXALPH(NCTAO+J,NVCUT+NOCUT+I)
 70      CONTINUE
C
      ELSE
       IF(MASWRK)  WRITE(6,*)'NO AO-CUT FOR NEXT ELONGATION STEP'
         NCTAO = 0
      ENDIF
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK AOCUT
      SUBROUTINE AOCUT(NCTAO,NEG4,D,NNN,DB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      PARAMETER (MXSH=5000, MXGTOT=20000,MXATM=2000, MXAO=8192)
      LOGICAL LCUT
      COMMON /ACUT  / NCTATM,NNNCTAO,NCTSHL,NCTELC,NCTGAS
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /CORE /  CORE(107)
      COMMON /ELGIDX/ LCUT
      DIMENSION D(NNN),DB(*)
      DIMENSION INFSHL(MXSH,5)
      DATA NONE/4HNONE/
      COMMON /ELGLAB/ IATM(MXAO)
      COMMON /MASSES/ ZMASS(MXATM)
      DATA RHF/8HRHF     /
CGU
      NNNCTAO = NCTAO
C====================================================================
C INFORMATION CONCERNING CURRENTLY CUT-OFF FRAGMENT
C====================================================================
      DO 10 I=1,NCTAO
         NCTATM = IATM(I)
 10   CONTINUE
C
      DO 20 I=1,NSHELL
         IF(KATOM(I).EQ.NCTATM) NCTSHL = I
 20   CONTINUE
C
      NCTELC = 0
      DO 30 I=1,NCTATM
         NCTELC = NCTELC + INT(ZAN(I))
         IF(MPCTYP.NE.NONE) THEN
            NUCZ = INT(ZAN(I))
            MPCORE = NUCZ - INT(CORE(NUCZ))
            NCTELC = NCTELC - MPCORE
         ENDIF
 30   CONTINUE
      NCTELC = NCTELC - ICH
      NADTL = -MOD(NCTELC,2)
      NCTELC = NCTELC + NADTL
C
      NCTGAS = KSTART(NCTSHL) + KNG(NCTSHL) - 1
C
      IF(MASWRK) THEN
      WRITE(6,'(60(1H=))')
      WRITE(6,'(8X,''INFORMATION ABOUT CURRENTLY CUT-OFF FRAGMENT'')')
      WRITE(6,'(60(1H-))')
      WRITE(6,'(5X,''ATOMS'',9X,''SHELLS'',7X,''ELECTRONS'',
     *                                             6X,''GAUSSIANS'')')
      WRITE(6,'(60(1H-))')
      WRITE(6,'(4(3X,I8,4X))') NCTATM,NCTSHL,NCTELC,NCTGAS
      WRITE(6,'(60(1H=))')
      END IF
C
C---------------------------------------------------------------------
      IF(LCUT) THEN
C---------------------------------------------------------------------
         CALL SEQREW(NEG4)
         READ(NEG4)
         READ(NEG4) NCTATX,NCTAOX,NCTSHX,NCTELX,NCTGSX,NUMX
C
         IF(NEG4.LT.0) THEN
            WRITE(6,*) 'STRICTLY TO FOOL FTNCHEK: NUMX=',NUMX
         END IF
C
         NFRGAT = NCTATX + NCTATM
         NFRGSH = NCTSHL + NCTSHX
         NFRGGS = NCTGAS + NCTGSX
         NFRGBS = NCTAOX + NCTAO
         NFRGEL = NCTELC + NCTELX
         NTOT = NCTAOX + NUM
         NRECT = NFRGBS*(NFRGBS+1) + NFRGBS*(NTOT-NFRGBS)
         L2 = NTOT*(NTOT+1)/2
         LL = MXGTOT
C---------------------------------------------------------------------
         IF(MASWRK) THEN
         WRITE(6,'(60(1H=))')
         WRITE(6,'(8X,''          TOTAL CUT-OFF INFORMATION'')')
         WRITE(6,'(60(1H-))')
         WRITE(6,'(5X,''ATOMS'',9X,''SHELLS'',7X,''ELECTRONS'',
     *                                             6X,''GAUSSIANS'')')
         WRITE(6,'(60(1H-))')
         WRITE(6,'(4(3X,I8,4X))') NFRGAT,NFRGSH,NFRGEL,NFRGGS
         WRITE(6,'(60(1H=))')
         END IF
C---------------------------------------------------------------------
         LOADFM=0
         CALL VALFM(LOADFM)
         LWRK =  LOADFM  +   1
         LL1  =  LWRK    +   L2
         LL2  =  LL1     +   NFRGBS
         LL3  =  LL2     +   NFRGAT
         LL4  =  LL3     +   NFRGAT
         LL5  =  LL4     +   NFRGAT
         LL6  =  LL5     +   NFRGAT
         LL7  =  LL6     +   NFRGGS
         LL8  =  LL7     +   NFRGGS
         LL9  =  LL8     +   NFRGGS
         LL10 =  LL9     +   NFRGGS
         LL11 =  LL10    +   NFRGGS
         LAST =  LL11    +   NFRGGS
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
         READ(NEG4) (XX(LL1+I-1),I=1,NCTATX)
         READ(NEG4) (XX(LL2+I-1),XX(LL3+I-1),XX(LL4+I-1),XX(LL5+I-1),
     *        I=1,NCTATX)
         READ(NEG4) (XX(LL6+I-1),XX(LL7+I-1),XX(LL8+I-1),XX(LL9+I-1),
     *        XX(LL10+I-1),XX(LL11+I-1),I=1,NCTGSX)
         READ(NEG4) ((INFSHL(I,J),J=1,5),I=1,NCTSHX)
         DO 40 I = NCTATX+1,NFRGAT
            XX(LL1 + I - 1) = ZMASS(I-NCTATX)
            XX(LL2 + I - 1) = ZAN(I-NCTATX)
            XX(LL3 + I - 1) = C(1,I-NCTATX)
            XX(LL4 + I - 1) = C(2,I-NCTATX)
            XX(LL5 + I - 1) = C(3,I-NCTATX)
 40      CONTINUE
         DO 50 I = NCTGSX+1,NFRGGS
            XX(LL6 + I - 1)  = EX(I-NCTGSX)
            XX(LL7 + I - 1)  = CS(I-NCTGSX)
            XX(LL8 + I - 1)  = CP(I-NCTGSX)
            XX(LL9 + I - 1)  = CD(I-NCTGSX)
            XX(LL10 + I - 1) = CF(I-NCTGSX)
            XX(LL11 + I - 1) = CG(I-NCTGSX)
 50      CONTINUE
         DO 60 I = NCTSHX+1,NFRGSH
            INFSHL(I,1)  =  KATOM(I-NCTSHX) + NCTATX
            INFSHL(I,2)  =  KTYPE(I-NCTSHX)
            INFSHL(I,3)  =    KNG(I-NCTSHX)
            INFSHL(I,4)  =   KMIN(I-NCTSHX)
            INFSHL(I,5)  =   KMAX(I-NCTSHX)
 60      CONTINUE
         CALL SEQREW(NEG4)
         WRITE(NEG4) NCTATM,NCTAO,NCTSHL,NCTELC,NCTGAS,NUM
         WRITE(NEG4) NFRGAT,NFRGBS,NFRGSH,NFRGEL,NFRGGS,NTOT
         WRITE(NEG4) (XX(LL1+I-1),I=1,NFRGAT)
         WRITE(NEG4) (XX(LL2+I-1),XX(LL3+I-1),XX(LL4+I-1),
     *        XX(LL5+I-1),I=1,NFRGAT)
         WRITE(NEG4) (XX(LL6+I-1),XX(LL7+I-1),XX(LL8+I-1),
     *        XX(LL9+I-1),XX(LL10+I-1),XX(LL11+I-1),
     *        I=1,NFRGGS)
         WRITE(NEG4) ((INFSHL(I,J),J=1,5),I=1,NFRGSH)
C
C CUT OFF PART OF D MATRIX
C
         II = 1
         JJ = 0
         DO 80 I = 1,NTOT
            DO 80 J = 1,I
               IF(J.GT.NFRGBS) GOTO 70
               IF(JJ.GT.(NRECT-1)) GOTO 70
               XX(LWRK+JJ) = D(II)
               JJ = JJ + 1
 70            II = II + 1
 80      CONTINUE
         WRITE(NEG4) (XX(LWRK+I-1),I=1,NRECT)
C
         IF(SCFTYP.NE.RHF) THEN
            II = 1
            JJ = 0
            DO 100 I = 1,NTOT
               DO 100 J = 1,I
                  IF(J.GT.NFRGBS) GOTO 90
                  IF(JJ.GT.(NRECT-1)) GOTO 90
                  XX(LWRK+JJ) = DB(II)
                  JJ = JJ + 1
 90               II = II + 1
 100        CONTINUE
            WRITE(NEG4) (XX(LWRK+I-1),I=1,NRECT)
         ENDIF
C
         CALL RETFM(NEED)
C---------------------------------------------------------------------
      ELSE
C---------------------------------------------------------------------
         CALL SEQREW(NEG4)
         WRITE(NEG4) NCTATM,NCTAO,NCTSHL,NCTELC,NCTGAS,NUM
         WRITE(NEG4) NCTATM,NCTAO,NCTSHL,NCTELC,NCTGAS,NUM
         WRITE(NEG4) (ZMASS(J),J=1,NCTATM)
         WRITE(NEG4) (ZAN(J),(C(I,J),I=1,3),J=1,NCTATM)
         WRITE(NEG4) (EX(I),CS(I),CP(I),CD(I),CF(I),CG(I),I=1,NCTGAS)
         WRITE(NEG4) (KATOM(I),KTYPE(I),KNG(I),KMIN(I),KMAX(I),
     *        I=1,NCTSHL)
         L2  = NUM*(NUM+1)/2
         LL  = (NUM-NCTAO)*NCTAO+NCTAO*(NCTAO+1)/2
         LOADFM=0
         CALL VALFM(LOADFM)
         LWRK = LOADFM + 1
         LDO  = LWRK   + L2
         LAST = LDO    + LL
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
C CUT OFF PART OF D MATRIX
C
         II = 1
         JJ = 0
         DO 120 I = 1,NUM
            DO 120 J = 1,I
               IF(J.GT.NCTAO) GOTO 110
               IF(JJ.GT.(LL-1)) GOTO 110
               XX(LDO+JJ) = D(II)
               JJ = JJ + 1
 110           II = II + 1
 120     CONTINUE
         WRITE(NEG4) (XX(LDO+I),I=0,LL-1)
C
         IF(SCFTYP.NE.RHF) THEN
            II = 1
            JJ = 0
            DO 140 I = 1,NUM
               DO 140 J = 1,I
                  IF(J.GT.NCTAO) GOTO 130
                  IF(JJ.GT.(LL-1)) GOTO 130
                  XX(LDO+JJ) = DB(II)
                  JJ = JJ + 1
 130              II = II + 1
 140        CONTINUE
            WRITE(NEG4) (XX(LDO+I),I=0,LL-1)
         ENDIF
C
         CALL RETFM(NEED)
C---------------------------------------------------------------------
      ENDIF
      RETURN
      END
C
C*MODULE ELGLOC  *DECK CUTSHT
      SUBROUTINE CUTSHT(D,NN,NUM,NFRGBS,NHB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(NN)
C
      JJ = 0
      II = 0
      DO 10 I = 1,NUM - NHB
         DO 10 J = 1,I
            II = II + 1
            IF(I.GT.NFRGBS.AND.J.GT.NFRGBS) THEN
               JJ = JJ + 1
               D(JJ) = D (II)
            ENDIF
 10   CONTINUE
C
      DO 20 I = JJ+1,NN
         D(I) = 0.0D+00
 20   CONTINUE
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK CTRUHF
      SUBROUTINE CTRUHF(SW,CWALPH,CXALPH,CWBETA,CXBETA,NUM,
     *          NFZOC,NFZVT,NFZOCB,NFZVTB,NCTAO,NCUT,NCUTX,
     *            NAAO,JWW,JXX,NCYCLE,NCT,NOCA,NOCXA,NOCUT,
     *                    NOCB,NOCXB,NOCUTB,NHB,DAOFF,DBOFF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION CWALPH(NUM,NUM),CXALPH(NUM,NUM),SW(NUM,NUM)
      DIMENSION CWBETA(NUM,NUM),CXBETA(NUM,NUM)
      DIMENSION DAOFF(NUM*(NUM+1)/2),DBOFF(NUM*(NUM+1)/2)
      PARAMETER (MXSTP=500)
      PARAMETER(CTHRSH=1.0D-05)
      DIMENSION NFZOC(MXSTP),NFZVT(MXSTP),NFZOCB(MXSTP),NFZVTB(MXSTP)
C
      IF(MASWRK) THEN
       WRITE(6,'(''PREPARING FOR AO-CUT,NCUTX='',I4)')NCUTX
       WRITE(6,'('' NSTEP   NFZOC  NFZVT   '',''  NFZOCB   NFZVTB'')')
         DO I=1,NCYCLE
          WRITE(6,'(I4,3X,4(I6,2X))')
     *              I,NFZOC(I),NFZVT(I),NFZOCB(I),NFZVTB(I)
         ENDDO
      ENDIF
C
      IF(NCYCLE.LE.1) RETURN
C
      NCUT = NCUTX
      NCTAO = NFZOC(NCUT+1) + NFZVT(NCUT+1)
      FCUTA = 0.0D+00
      DO 20 I=1,JWW
         DO 20 K=1,NCTAO
            DO 20 L=1,NCTAO
               FCUTA = FCUTA + ABS(CWALPH(K,I)*SW(K,L)*CWALPH(L,I))
 20   CONTINUE
C
      NCTAOB = NFZOCB(NCUT+1) + NFZVTB(NCUT+1)
      IF(NCTAOB.NE.NCTAO) CALL ABRT
      FCUTB = 0.0D+00
      DO 30 I=1,JWW
         DO 30 K=1,NCTAO
            DO 30 L=1,NCTAO
               FCUTB = FCUTB + ABS(CWBETA(K,I)*SW(K,L)*CWBETA(L,I))
 30   CONTINUE
       IF(MASWRK) WRITE(6,'(''FCUTA,FCUTB,CTHRSH='',3(E10.3,2X))')
     *                        FCUTA,FCUTB,CTHRSH
C
      IF(FCUTA.LT.CTHRSH.AND.FCUTB.LT.CTHRSH.AND.NCT.NE.0) THEN
         IF(MASWRK)
     *      WRITE(6,*)'AO-CUT CAN BE PERFORMED IN THE NEXT ELG STEP'
C
C  PERFORM AO-CUT IN B-REGION
C
         NCUT = NCUT + 1
         NAAO  = NAAO  - NCTAO
         NOCUT = NFZOC(NCUT)
         NVCUT = NFZVT(NCUT)
         NOCRST = NOCA + NOCXA - NOCUT
         NOCUTB = NFZOCB(NCUT)
         NVCUTB = NFZVTB(NCUT)
         NOCRSTB = NOCB + NOCXB - NOCUTB
         IF(MASWRK) THEN
         WRITE(6,*)
         WRITE(6,'(''     NOCUT,  NVCUT = '',I4,2X,I4)')  NOCUT,NVCUT
         WRITE(6,'(''    NOCUTB, NVCUTB = '',I4,2X,I4)') NOCUTB,NVCUTB
         WRITE(6,'(''    NOCRST,NOCRSTB = '',I4,2X,I4)')NOCRST,NOCRSTB
         WRITE(6,*)
         WRITE(6,'(''  NCTAO,NCUT,NCUTX ='',3(I4,2X))')NCTAO,NCUT,NCUTX
         WRITE(6,'(''  NOCA,NOCXA,NOCUT ='',3(I4,2X))')NOCA,NOCXA,NOCUT
         WRITE(6,'('' NOCB,NOCXB,NOCUTB ='',3(I4,2X))')NOCB,NOCXB,NOCUTB
         WRITE(6,'('' NAAO,JWW,JXX,NCYCLE,NCT,NHB ='',6(I4,2X))')
     *                NAAO,JWW,JXX,NCYCLE,NCT,NHB
         END IF
C
         DO 40 I=1,JWW
            DO 40 J=1,NAAO
               CWALPH(J,I) = CWALPH(NCTAO+J,I)
               CWBETA(J,I) = CWBETA(NCTAO+J,I)
 40      CONTINUE
CGU
C    CALCULATE CUT-OFF DENSITY MATRIX
CGU
         II = 0
         DO 50 I=1,NUM
            DO 50 J=1,I
               II = II + 1
               DAOFF(II) = 0.0D+00
               DO 50 K=1,NOCUT
                  DAOFF(II) = DAOFF(II) + CXALPH(I,K)*CXALPH(J,K)
 50      CONTINUE
C
         II = 0
         DO 60 I=1,NUM
            DO 60 J=1,I
               II = II + 1
               DBOFF(II) = 0.0D+00
               DO 60 K=1,NOCUTB
                  DBOFF(II) = DBOFF(II) + CXBETA(I,K)*CXBETA(J,K)
 60      CONTINUE
C
         II = 0
         JJ = 0
         DO 70 I=1,NUM
            DO 70 J=1,I
               II = II  + 1
               IF(J.GT.NCTAO.AND.I.GT.NCTAO) THEN
                  JJ = JJ + 1
                  DAOFF(JJ)  = DAOFF(II)
                  DBOFF(JJ)  = DBOFF(II)
               ENDIF
 70      CONTINUE
         DO 80 J=JJ+1,NUM*(NUM+1)/2
            DAOFF(J)=0.0D+00
            DBOFF(J)=0.0D+00
 80      CONTINUE
C
C   PERFORM MO-CUT IN A-REGION
C
         JXX = JXX - NCTAO
         DO 120 J=1,NAAO
            DO 90 I=1,NOCRST
               CXALPH(J,I) = CXALPH(NCTAO+J,NOCUT+I)
 90         CONTINUE
            DO 100 I=NOCRST+1,JXX
               CXALPH(J,I) = CXALPH(NCTAO+J,NOCUT+NVCUT+I)
 100        CONTINUE
            DO 110 I=1,NOCRSTB
               CXBETA(J,I) = CXBETA(NCTAO+J,NOCUTB+I)
 110        CONTINUE
            DO 120 I=NOCRSTB+1,JXX
               CXBETA(J,I) = CXBETA(NCTAO+J,NOCUTB+NVCUTB+I)
 120     CONTINUE
C
      ELSE
         IF(MASWRK) WRITE(6,*)'NO AO-CUT FOR NEXT ELONGATION STEP'
         NCTAO = 0
      ENDIF
C
      RETURN
      END
C
C*MODULE ELGLOC  *DECK DAO
      SUBROUTINE DAO(D,C,NOCC,NOC1X,N1,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N1*(N1+1)/2),C(N1,N)
C
      OCC = 2.0D+00
      CALL VCLR(D,1,N1*(N1+1)/2)
      IJ = 0
      DO 10 I=1,N1
         DO 10 J=1,I
            IJ = IJ + 1
              D(IJ) = 0.0D+00
            DO K=1,NOCC-NOC1X
               D(IJ) = D(IJ) + C(I,K)*C(J,K)*OCC
            ENDDO
 10   CONTINUE
      RETURN
      END
C*MODULE ELGLOC  *DECK DAOVAC
      SUBROUTINE DAOVAC(D,C,NOCC,NOC1X,N1,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N1*(N1+1)/2),C(N1,N)
C
      OCC = 2.0D+00
      CALL VCLR(D,1,N1*(N1+1)/2)
      IJ = 0
      DO 10 I=1,N1
         DO 10 J=1,I
            IJ = IJ + 1
               D(IJ) = 0.0D+00
            DO K=1+NOCC-NOC1X,N
               D(IJ) = D(IJ) + C(I,K)*C(J,K)*OCC
            ENDDO
 10   CONTINUE
      RETURN
      END
C
C*MODULE ELGLOC  *DECK BDJAC
      SUBROUTINE BDJAC(A, N, NP, M, D, V, NROT, SUCC, SM, NEEDA)
C
C     ROUTINE FOR BLOCK DIAGONALIZATION BASED ON JACOBI SCHEME
C     A - SYMMETRIC MATRIX OF LOGICAL SIZE N X N AND
C       PHYSICAL SIZE NP X NP
C       ON OUTPUT UPPER TRIANGULAR PART OF A CONTAINS TRANSFORMED
C                       MATRIX (SEE ALSO NEEDA)
C     M - SIZE OF THE FIRST BLOCK (M X M)
C     D - RETURNS PSEUDOEIGENVALUES
C     V - THE MATRIX OF THE SAME SIZES AS A, ON OUTPUT CONTAINS
C       PSEUDOEIGENVECTORS IN COLUMNS
C     NROT - NUMBER OF ROTATIONS MADE
C     SUCC - .TRUE. IF COMPLETED SUCCESSFULLY, .FALSE. IN ANOTHER CASE
C     SM - SUM OF ABSOLUTE VALUES OF OFF-BLOCK-DIAGONAL ELEMENTS
C     NEEDA - IF .FALSE. ALSO DIAGONAL PART OF A WILL CONTAIN ON OUTPUT
C        THE TRANSFORMED ELEMENTS,
C       IF .TRUE. ORIGINAL ONES WILL BE SAVED
C
      IMPLICIT NONE
      INTEGER N, M, NP, NROT, MXAO, MAXROT
      LOGICAL SUCC, NEEDA
      DOUBLE PRECISION SM
      DOUBLE PRECISION A(NP, NP), D(NP), V(NP, NP)
C
      INTEGER I, IP, JP, J, II
      DOUBLE PRECISION C, G, H, S, T, TAU, THETA, TRESH, ZERO, ONE
      PARAMETER (MXAO=8192,MAXROT=100)
      PARAMETER(ZERO = 0.0D+00, ONE = 1.0D+00)
      DOUBLE PRECISION B(MXAO), Z(MXAO)
C
      SUCC = .FALSE.
      DO IP = 1, N
         DO JP = 1, N
            V(IP, JP) = ZERO
         END DO
         V(IP, IP) = ONE
      END DO
      DO IP = 1, N
         B(IP) = A(IP, IP)
         D(IP) = B(IP)
         Z(IP) = ZERO
      END DO
      NROT = 0
      DO I = 1, MAXROT
         SM = ZERO
         DO IP = 1, M
            DO JP = M + 1, N
               SM = SM + ABS(A(IP, JP))
            END DO
         END DO
C
         IF(SM.LT.1.0D-10) THEN
            IF (.NOT.NEEDA) THEN
               DO II = 1, N
                  A(II, II) = D(II)
               END DO
            END IF
            SUCC = .TRUE.
            RETURN
         END IF
         IF(I.LT.4) THEN
            TRESH = 2.0D-01 * SM / N ** 2
         ELSE
            TRESH = ZERO
         END IF
         DO IP = 1, M
            DO JP = M + 1, N
               G = 1.D2 * ABS(A(IP, JP))
               IF((I.GT.4).AND.((ABS(D(IP)) + G).EQ.ABS(D(IP))).AND.
     N              ((ABS(D(JP)) + G).EQ.ABS(D(JP)))) THEN
                  A(IP,JP) = ZERO
               ELSE IF(ABS(A(IP, JP)).GT.TRESH) THEN
                  H = D(JP) - D(IP)
                  IF((ABS(H) + G).EQ.ABS(H)) THEN
                     T = A(IP, JP) / H
                  ELSE
                     THETA = 5.0D-01 * H / A(IP, JP)
                     T = ONE / (ABS(THETA) + SQRT(ONE + THETA ** 2))
                     IF(THETA.LT.ZERO) T = -T
                  END IF
                  C = ONE / SQRT(ONE + T ** 2)
                  S = T * C
                  TAU = S / (ONE + C)
                  H = T * A(IP, JP)
                  Z(IP) = Z(IP) - H
                  Z(JP) = Z(JP) + H
                  D(IP) = D(IP) - H
                  D(JP) = D(JP) + H
                  A(IP, JP) = ZERO
                  DO J = 1, IP - 1
                     G = A(J, IP)
                     H = A(J, JP)
                     A(J, IP) = G - S * (H + G * TAU)
                     A(J, JP) = H + S * (G - H * TAU)
                  END DO
                  DO J = IP + 1, JP - 1
                     G = A(IP, J)
                     H = A(J, JP)
                     A(IP, J) = G - S * (H + G * TAU)
                     A(J, JP) = H + S * (G - H * TAU)
                  END DO
                  DO J = JP + 1, N
                     G = A(IP, J)
                     H = A(JP, J)
                     A(IP, J) = G - S * (H + G * TAU)
                     A(JP, J) = H + S * (G - H * TAU)
                  END DO
                  DO J = 1,N
                     G = V(J, IP)
                     H = V(J, JP)
                     V(J, IP) = G - S * (H + G * TAU)
                     V(J, JP) = H + S * (G - H * TAU)
                  END DO
                  NROT = NROT + 1
               END IF
            END DO
         END DO
         DO IP = 1, N
            B(IP) = B(IP) + Z(IP)
            D(IP) = B(IP)
            Z(IP) = ZERO
         END DO
      END DO
      IF (.NOT.NEEDA) THEN
         DO II = 1, N
            A(II, II) = D(II)
         END DO
      END IF
      RETURN
      END
C
C*MODULE ELGLOC  *DECK SDGX
      SUBROUTINE SDGX(ST,X,Y,EIG,VV,N1,L2,WRK,IWRK)
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DATA NONE/4HNONE/
C
      DIMENSION ST(L2),X(N1,N1),Y(N1,N1),EIG(N1),VV(N1,N1),
     *          WRK(N1,8),IWRK(N1)
C----------------------------------------------------------------
      IF(MPCTYP.EQ.NONE) THEN
C
      CALL GLDIAG(N1,N1,N1,ST,WRK,EIG,VV,IERR,IWRK)
C
        IF(IPRI.EQ.3.AND.MASWRK) THEN
         WRITE(6,*)'EIGENVALUES/VECTORS OF WHOLE AO-BASIS S-MATRIX'
         CALL PREVNL(VV,EIG,N1,N1,N1)
        ENDIF
C
         DO 10 I=1,N1
            IF(EIG(I).LE.1.0D-08) THEN
               IF(MASWRK) WRITE(6,500) I,EIG(I)
               CALL ABRT
            END IF
 10      CONTINUE
C
C     X = V * S( 1/2) * V(DAGGER)
C     Y = V * S(-1/2) * V(DAGGER)
C
         DO 30 I=1,N1
            DO 30 J=1,N1
               X(I,J) = 0.0D+00
               Y(I,J) = 0.0D+00
               DO 20 K=1,N1
                  X(I,J) = X(I,J) + VV(I,K)*SQRT(EIG(K))*VV(J,K)
                  Y(I,J) = Y(I,J) + VV(I,K)*VV(J,K)/SQRT(EIG(K))
 20            CONTINUE
 30      CONTINUE
C
         IF(IPRI.EQ.3.AND.MASWRK) THEN
            WRITE(6,*)'X-MAT IN RAOLMO'
            CALL WRITN4(N1,N1,X,N1,N1)
         ENDIF
C
      ENDIF
C
 500  FORMAT('EIGENVALUE OF S-MAT IS NEGATIVE: I,S(I,I)=',
     *     I4,2X,F10.6)
      RETURN
      END
