C  2 FEB 06 - MWS - POLNMX: USE PRCALC TO DO THE DIPOLE INTS
C 17 JAN 06 - MWS - UPDATE ARGS TO TRFMCX
C 14 NOV 05 - DGF - PAD COMMON BLOCK SCFOPT
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 27 JUN 05 - IA  - LOCINP: REMOVE 1 JUN 05 TRAP AFTER FIXING THE BUG
C  1 JUN 05 - MWS - LOCINP: TURN OFF DYNAMIC POLARIZATION IF PARALLEL
C 30 APR 05 - DGF - SYNCHRONISE SYMBLK COMMON BLOCK
C  7 MAR 05 - IA  - LMOINP: READ IMAG FREQS FOR DYN. POLARIZABILITY
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, EDCMP
C  5 FEB 05 - MWS - PAD COMMON MCINP
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C  7 APR 03 - MWS - SYNCH FCCWFN COMMON
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 03 - MWS - SYNCH FCCWFN COMMON
C 14 JAN 03 - MWS - LOCFRZ,LMOFRZ: ENSURE PSI INVARIANT ORMAS-MCSCF LOCS
C 12 DEC 02 - MWS - ABORT LOCALIZATIONS W/O ANY ORBITALS TO LOCALIZE
C 17 AUG 02 - MWS - LMOFRZ: SYMLOC JOBS SHOULD NOT CHANGE ORBITAL PHASES
C  7 AUG 02 - PM  - LMOINP: NEW KEYWORD NODENS CONTROLS MAKING EFP
C 22 MAY 02 - MWS - PROVIDE ADDITIONAL TRFMCX ARGUMENTS
C 24 JAN 02 - JI  - DIRLMO,DIMOID,ORIPRI,ORIANAL: FIX NORBS.LT.NATOMS
C 30 NOV 01 - JI  - ORIEN: TINY FIX, LMOINP: TURN OFF ORIENT IF NOT MC
C 26 OCT 01 - JI  - IMPLEMENT ORIENTED (DIRECTED) LOCALIZED ORBITALS
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C  1 AUG 01 - JI  - TWEAKS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT AND WFNOPT
C 13 JUN 01 - PM  - LMOINP: PERMIT INPUT ABOUT OVERLAPPING REGION
C 19 NOV 00 - JHJ,HL - LMOINP: POLDCM DEFAULT MOVED, POLAPP OPTION ADDED
C 11 OCT 00 - MWS - UPDATE THE DETWFN COMMON
C  1 MAY 00 - MWS - NO SYMMETRY OPTION IN TRFOPT COMMON
C 16 FEB 00 - VK  - POLNMM,POLNMX: NEW CODE; LMOINP: NEW POLNUM KEYWORD
C 21 DEC 99 - MWS - REMOVE C1 DOWNSHIFT FROM TRANSF. BUT NOT AO INTS
C 13 MAR 99 - MWS - LOCSET: FIX CHECK RUN DAF RECORD EXISTENCE BUG
C 13 FEB 99 - MWS - LOCFRZ: CHANGE CALL TO SYMMOS
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C 12 SEP 98 - MWS - LMOPSI: CHANGES FOR DETERMINANT MCSCF
C 13 APR 98 - MWS - LMOFRZ,LOCFRZ: READ SYMMETRY LABELS FROM DAF
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 29 APR 97 - MWS - LOCSET: SKIP EIG/DEN IF RUN TYPE IS PROPERTY
C 18 MAR 97 - PMD - BOYLOC: ALLOW UP TO FIVE RESTARTS
C 14 FEB 97 - MWS - LOCENG,LOCSET MODEIFIED FOR NEW LOCANA ROUTINE
C  7 FEB 97 - MWS - LOCENG: SET DKSWRK PROPERLY FOR PARALLEL RUNS
C 17 JAN 97 - MWS - LOCENG: FIX BUG WITH INTEGRAL SKIPPING FOR MCSCF
C 18 DEC 96 - JHJ - LMOINP: IF RUNTYP IS MAKEFP, POLDCM=.T.
C  9 DEC 96 - JHJ - LMOINP,LMOLOC,BOYLOC: NPRINT=-23 OPTION
C 18 JUL 96 - MWS - LOCENG: FIX ENERGY LOC. FOR DIRECT RUNTYP=PROP JOBS
C 13 JUN 96 - MWS - CHANGES TO SUPPORT NEW CITYP INPUT
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C 14 SEP 95 - SPW - LMOINP,LMOX: CHANGES TO POLARIZABILITY DECOMPOSITION
C 26 JUL 95 - JHJ - LOCENG: IF(RDINT) WRITE BOGUS TLOC MATRIX TO DAF 73.
C 24 MAY 95 - MWS - LOCENG: WORK WITH DISTRIBUTED AO INTEGRALS
C 19 MAY 95 - SPW - BOYLOC,POPLOC: SAVE ALL VIRTUALS
C 21 APR 95 - JHJ - MLOINP: READ NEW ZDO OPTION
C  1 FEB 95 - JHJ - LMOINP: READ QUADRUPOLE DECOMPOSITION OPTION
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 14 OCT 94 - MWS - GLOBAL SUM IN LOCINT, ALL NODES MUST ENTER LOCLIZ
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  5 AUG 94 - MWS - LOCINT,LOCENG: USE DOUBLE LABEL PACKING
C 17 JUN 94 - JHJ - LMOINP: NEW INPUT FOR EDCOMP, LOCENG: READ IN $TWOEI
C                   LMOX: CALLS TO MOIDM,EDCMPM,DIPLMM
C                   LOCLIZ: PUNCH $TWOEI, WRITE TRANSF. INTEGRALS TO DAF
C                   BOYLOC,LOCLIZ,POPLOC: WRITE TRANSF. MATRIX TO DAF
C  1 JUN 94 - MWS - LOCENG: ENABLE DIRECT AND/OR PARALLEL LOCALIZING
C 17 DEC 93 - MWS - BOYLOC: CHANGE DAF RECORD NOS.
C  3 MAR 93 - MWS - LOCENG: PROP RUNS NEED TO DO 2E- INTEGRALS
C  8 SEP 92 - MWS - LMOPSI: CHANGE NUMCOR FOR LU AND LR
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  3 MAR 92 - MWS - LMOFRZ,LOCFRZ: CHANGE SYMLOC FREEZE TO USE SYMMOS
C 21 JAN 92 - MWS - UPDATE CREDITS FOR RUEDENBERG LOCALIZATION
C 11 JAN 92 - TLW - LOCINT: MAKE READS PARALLEL
C 10 JAN 92 - TLW - CHANGE REWINDS TO CALL SEQREW
C 10 JAN 92 - MWS,TLW - CHANGE OPENPK TO SEQOPN
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C  5 DEC 91 - MWS - LOCFRZ: FIX RUEDENBERG ORB FREEZING BUG
C  7 OCT 91 - MWS - CHANGE UNIX BYTE UNPACKING
C 16 AUG 91 - JAM - INCLUDE PIPEK/MEZEY POPULATION LOCALIZATION
C 24 APR 91 - MWS - COMPUTE NUMBER OF VIRTUALS IN BOYLOC, FIX CHECK RUNS
C  5 APR 91 - MWS - INCORPORATE ENERGY LOCALIZATION FROM ALIS,
C                   CLEAN UP BOYS, CHANGES TO INPUT (MOIN,SYMLOC),
C                   IMPLEMENT SINGLE JOB WAVEFUNCTION LOCALIZATION.
C 29 DEC 89 - MWS - ELIMINATE EQUIVALENCE OF IX TO /FMCOM/
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C 14 NOV 87 - STE - ADD TLOCAL TO BOYIN CALL, USE EXETYP
C  3 NOV 86 - STE - BOYLOC: USE IF/THEN LOGIC; REMOVE NAV REFERENCES
C                   USE PARAMETER STATEMENTS
C 11 JUL 86 - MWS - SAVE U IN LMORND FOR ALL MACHINES
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS,
C                   CHANGE CONVERSION FACTOR TO DEBYES
C 14 OCT 85 - STE - USE GENERIC MAX,ABS,ATAN,ATAN2,COS,SIN,SQRT
C 10 OCT 85 - STE - BOYLOC: USE DROT INSTEAD OF 900 LOOP
C  9 SEP 85 - STE - BOYIN: FIX NNAM,RNAM,JRET; DELETE 9028 FORMAT
C                   BOYLOC: FIX LOADFM. BOYLOC: FIX ZERO,NORB,K,KL,ONE
C                   LMORND: SAVE U, DFLOAT TO REAL, IDINT TO INT
C 19 FEB 85 - MWS - MOVE BOYIN FROM INPUT TO HERE
C 19 FEB 85 - STE - FIX 9030 FORMAT IN BOYLOC
C 15 NOV 84 - MWS - ALLOW VARIABLE NUMBER OF ITERS, CLEANER OUTPUT
C 15 AUG 84 - SS  - PRINT CENTROIDS OF CHARGE
C 11 FEB 83 - MWS - EXTENSIVE REWORKING, CLEANED UP, PRINT THE
C                   LOCALIZATION SUMS, TRANSFORMATION, ETC.
C 29 SEP 82 - MWS - CONVERT TO IBM
C
C     ROUTINE NAMES IMPLY: LMO*=SHARED BY MORE THAN ONE METHOD,
C         BOY*=BOYS, LOC*=RUEDENBERG, POP*=POPULATION
C
C*MODULE LOCAL   *DECK BOYLOC
      SUBROUTINE BOYLOC(IPASS,VECIN,TRAN,VECOUT,IORD,IIR,
     *                  IA,RIJ,QPIX,QPJX,NOROT,DMAO,WRK,
     *                  L1,L2,N1,N2,NOUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION LMORND
C
      DIMENSION VECIN(L1,L1),TRAN(N1,N1),VECOUT(L1,L1),DMAO(L2),
     *          WRK(L1),RIJ(N2,3),IORD(N1),IIR(N1),
     *          IA(L1),QPIX(L1),QPJX(L1),NOROT(N1,N1)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00,
     *           TENM3=1.0D-03, TENM8=1.0D-08, TENM10=1.0D-10,
     *           FAC=2.541766D+00)
C
      DATA UHF/8HUHF     /
C
C     ----- BOYS ORBITAL LOCALIZATION -----
C     THIS IMPLEMENTATION IS ORIGINALLY FROM QCPE PROGRAM 354.
C
      NORB  = N1
      NBASIS= L1
      NFROZ = NOUT
C                THIS NEXT IS DREADFUL, IT IS OVERWRITING SCF COUNT
C                I GUESS THIS WORKS BECAUSE LOCALIZATION IS DEAD LAST.
      MAXIT = MAXLOC
C
      DO 100 I=1,NBASIS
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C        BEGIN LOCALIZATION RETRYS HERE.
C
      NREDO = 0
  200 CONTINUE
      NREDO = NREDO + 1
C
C        TRANSFORM DIPOLE INTEGRALS FROM AO TO MO BASIS
C
      CALL DAREAD(IDAF,IODA,DMAO,L2,95,0)
      CALL TFTRI(RIJ(1,1),DMAO,VECIN,WRK,NORB,NBASIS,NBASIS)
      CALL DAREAD(IDAF,IODA,DMAO,L2,96,0)
      CALL TFTRI(RIJ(1,2),DMAO,VECIN,WRK,NORB,NBASIS,NBASIS)
      CALL DAREAD(IDAF,IODA,DMAO,L2,97,0)
      CALL TFTRI(RIJ(1,3),DMAO,VECIN,WRK,NORB,NBASIS,NBASIS)
C
C                    COMPUTE INITIAL LOCALIZATION SUM
C
      SUMRR = ZERO
      II=0
      DO 260 I=1,NORB
         II=II+I
         IF(NOROT(I,I).EQ.1) GO TO 260
         SUMRR = SUMRR + RIJ(II,1)**2 + RIJ(II,2)**2 + RIJ(II,3)**2
  260 CONTINUE
      SUMRR=SUMRR*FAC*FAC
      IF (MASWRK .AND. NPRINT.NE.-23) WRITE(IW,9010) SUMRR
C
C        SEED THE RANDOM FUNCTION
C
      ONEPT0 = ONE
      IF (NREDO.EQ.1) XX = LMORND(ONEPT0,VECIN,NBASIS)
      SHIFT = ATAN(ONEPT0)
C
C        INITIALIZE LOCALIZATION TRANSFORM
C
      CALL VCLR(TRAN,1,NORB*NORB)
      DO 280 I = 1,NORB
         TRAN(I,I) = ONE
  280 CONTINUE
C
C        BEGIN LOCALIZATION CYCLES
C
      ITER = 0
  300 CONTINUE
      CHANGE = ZERO
      ITER = ITER+1
      DO 320 I = 1,NORB
         IIR(I) = I
  320 CONTINUE
      NNN = NORB
      DO 330 I = 1,NORB
         XX = LMORND(CHANGE,VECIN,NBASIS)
         III = INT(XX*NNN+ONE)
         IORD(I) = IIR(III)
         IIR(III) = IIR(NNN)
         NNN = NNN-1
  330 CONTINUE
C
C        FOR EACH PAIR OF ORBITALS A 2X2 UNITARY TRANSFORMATION
C        IS PERFORMED. THE JACOBI TYPE TRANSFORMATION IS
C             PSI PRIME(I) =  COS(T)*PSI(I) + SIN(T)*PSI(J)  AND
C             PSI PRIME(J) = -SIN(T)*PSI(I) + COS(T)*PSI(J).
C        THE BOYS METHOD REQUIRES THAT T BE SUCH AS TO MAXIMIZE
C        THE SUM OF THE SQUARES OF THE ONE CENTER MOLECULAR ORBITAL
C        DIPOLE MOMENT INTEGRALS.
C
      DO 600 III = 1,NORB
         I = IORD(III)
         IF(NOROT(I,I).EQ.1) GO TO 600
C
         II = IA(I)+I
         IM  = I
         JM  = 1
         IJM = IA(IM)+JM
         RM  = ZERO
         TM  = ZERO
         SM  = ZERO
         CM  = ONE
C
C           LOOK FOR MAXIMUM ROTATION INVOLVING ORBITAL I AND ANY J
C
         DO 490 J = 1,NORB
            IF(I.EQ.J) GO TO 490
            IF(NOROT(I,J).EQ.1) GO TO 490
            IJ = IA(I) + J
            JJ = IA(J) + J
            IF(J.GT.I) IJ = IA(J) + I
            T = ZERO
            TX = ZERO
            DO 410 KK = 1,3
               T = T + FOUR*RIJ(IJ,KK)**2
     *               - RIJ(II,KK)**2 - RIJ(JJ,KK)**2
     *               + TWO*RIJ(II,KK)*RIJ(JJ,KK)
               TX = TX+RIJ(IJ,KK)*(RIJ(JJ,KK)-RIJ(II,KK))
  410       CONTINUE
            IF((ABS(T).LE.TENM10) .AND. (ABS(TX).LE.TENM10)) GO TO 490
            TX = FOUR*TX
            T = ATAN2(TX,T)
            T = T/FOUR
            SIGN = ONE
            IF (T .GT. ZERO) SIGN = -ONE
            T = T+SIGN*SHIFT
            ITIM = 0
  420       CONTINUE
            ITIM = ITIM+1
            S = SIN(T)
            C = COS(T)
            RIN = ZERO
            DO 430 KK = 1,3
               QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
               QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
               RIN = RIN+QPI*QPI+QPJ*QPJ-RIJ(II,KK)**2-RIJ(JJ,KK)**2
  430       CONTINUE
            TTEST = ABS(T)-SHIFT
            IF((ABS(T).LE.TENM8) .OR. (ABS(TTEST).LE.TENM8)) GO TO 470
            IF(RIN .GE. -TENM8) GO TO 470
            IF(ITIM .LE. 1) GO TO 460
            IF (MASWRK) THEN
               WRITE (IW,9020) I,J
               WRITE (IW,9030) T,S,C,RIN
            END IF
            RETURN
C
  460       CONTINUE
            SIGN = ONE
            IF (T .GT. ZERO) SIGN = -ONE
            T = T+SHIFT*SIGN
            GO TO 420
C
  470       CONTINUE
            IF (RIN .LE. RM) GO TO 490
            RM  = RIN
            JM  = J
            IJM = IJ
            SM  = S
            CM  = C
            TM  = T
  490    CONTINUE
C
         T   = TM
         RIN = RM
         S   = SM
         C   = CM
         J   = JM
         IJ  = IJM
         JJ  = IA(J)+J
         IF(NOROT(I,J).EQ.1) GO TO 600
C
C           ACCUMULATE THE 2X2 ROTATION
C
         CHANGE = CHANGE+T*T
         CALL DROT(NORB,TRAN(1,I),1,TRAN(1,J),1,C,S)
C
C           UPDATE THE DIPOLE INTEGRALS FOR THIS ROTATION
C
         DO 550 KK = 1,3
            QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
            QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
            QPIJ = (C*C-S*S)*RIJ(IJ,KK)+C*S*(RIJ(JJ,KK)-RIJ(II,KK))
            DO 510 K = 1,NORB
               IF(I.EQ.K  .OR.  J.EQ.K) GO TO 510
               IK = IA(I) + K
               IF(K.GT.I) IK = IA(K)+I
               JK = IA(J) + K
               IF(K.GT.J) JK = IA(K)+J
               QPIX(K) = C*RIJ(IK,KK)+S*RIJ(JK,KK)
               QPJX(K) = C*RIJ(JK,KK)-S*RIJ(IK,KK)
               RIJ(IK,KK) = QPIX(K)
               RIJ(JK,KK) = QPJX(K)
  510       CONTINUE
            RIN = RIN+QPI+QPJ-RIJ(II,KK)-RIJ(JJ,KK)
            RIJ(II,KK) = QPI
            RIJ(JJ,KK) = QPJ
            RIJ(IJ,KK) = QPIJ
  550    CONTINUE
  600 CONTINUE
C
C        TEST FOR CONVERGENCE OF THE LOCALIZATION
C
      CHANGE = SQRT(TWO*CHANGE/(NORB*(NORB-1)))
      IF((MOD(ITER,10).EQ.0) .AND. MASWRK .AND. NPRINT.NE.-23)
     *                                    WRITE(IW,9040) ITER,CHANGE
C
      IF(ITER.LT.MAXIT  .AND.  CHANGE.GT.TENM3*CVGLOC) GO TO 300
      IF(CHANGE.GT.CVGLOC) THEN
         IF(NREDO.LE.5) THEN
            IF (MASWRK) WRITE(IW,9050) NREDO
            GO TO 200
         ELSE
            IF (MASWRK) WRITE(IW,9060)
            RETURN
         END IF
      END IF
C
C        FINISHED WITH LOCALIZATION CYCLES
C
      IF (MASWRK .AND. NPRINT.NE.-23) WRITE(IW,9070) ITER
C
C        COMPUTE FINAL LOCALIZATION SUM
C
      SUMRR = ZERO
      II=0
      DO 710 I=1,NORB
         II=II+I
         IF(NOROT(I,I).EQ.1) GO TO 710
         SUMRR = SUMRR + RIJ(II,1)**2 + RIJ(II,2)**2 + RIJ(II,3)**2
  710 CONTINUE
      SUMRR=SUMRR*FAC*FAC
      IF (MASWRK .AND. NPRINT.NE.-23) WRITE(IW,9080) SUMRR
C
C        TRANSFORM TO FINAL ORBITALS, COPY VIRTUAL SPACE
C
      CALL MRARBR(VECIN,NBASIS,NBASIS,NORB,TRAN,NORB,NORB,VECOUT,NBASIS)
      NVIRT = NBASIS - NORB
      CALL DCOPY(NVIRT*NBASIS,VECIN(1,NORB+1),1,VECOUT(1,NORB+1),1)
C
C        OUTPUT LOCALIZATION TRANSFORMATION AND FINAL ORBITALS
C
      IF(IPRTLO.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9090)
         CALL PRSQ(TRAN,NORB,NORB,NORB)
      END IF
C     ----- SAVE THIS FOR EVENTUAL MP2 LOCALIZATION -----
C     WE NEED THE MATRIX WITHOUT THE CORE
C
      DO 760 I = 1,NORB
         DO 750 J = 1,NORB
            TRAN(I,J) = TRAN(I+NOUT,J+NOUT)
  750    CONTINUE
  760 CONTINUE
      NLOC = NORB - NOUT
      CALL SQZMAT(TRAN,NORB,NLOC,NORB)
      IF(IPASS.EQ.1) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,73,0)
      IF(IPASS.EQ.2) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,74,0)
C
      IF (MASWRK .AND. NPRINT.NE.-23) THEN
         WRITE(IW,9100)
         CALL PRSQL(VECOUT,NORB,NBASIS,NBASIS)
      END IF
C
      IF((IPASS.EQ.1) .AND. MASWRK) WRITE(IP,8000)
      CALL PUSQL(VECOUT,NORB,NBASIS,NBASIS)
      IF((SCFTYP.NE.UHF  .OR.  IPASS.EQ.2) .AND. MASWRK)
     *    WRITE(IP,8010)
C
C        OPTIONALLY, PRINT ORBITAL CENTROIDS
C     THIS IS JUST A FANCY NAME FOR THE ORBITAL DIPOLE MOMENT.
C
      IF(IPRTLO.EQ.0) RETURN
      IF (MASWRK) THEN
      WRITE (IW,9110)
      II = 0
      DO 850 I = 1,NORB
         II = II + I
         DO 820 K=1,NFROZ
            IF(IPASS.EQ.1) THEN
               KK=MOOUTA(K)
            ELSE
               KK=MOOUTB(K)
            END IF
            IF(KK.EQ.I) GO TO 850
  820    CONTINUE
         WRITE(IW,9120) I,(RIJ(II,IMU),IMU=1,3)
  850 CONTINUE
      END IF
      RETURN
C
 8000 FORMAT('BOYS LOCALIZED ORBITALS'/' $VEC')
 8010 FORMAT(' $END')
C                           PRINT FORMATS
 9020 FORMAT(/1X,'NO ROTATION INCREASES DIPOLE INTEGRALS',
     *           ' --- LOCALIZATION ABORTED'/
     *           10X,2HI=,I3,5X,2HJ=,I3)
 9030 FORMAT(5X,8HTHETA = ,1P,E20.10,0P/
     *          5X,12HSIN(THETA)= ,F10.7,' COS(THETA)=',F10.7/
     *          5X,'TOTAL CHANGE TO THIS POINT =',1P,E20.10)
 9040 FORMAT(10X,'BOYS ITERATION',I4,' ORBITAL CHANGE=',1P,E20.10)
 9050 FORMAT (/5X,'**** LOCALIZATION',I4,' HAS BEEN UNSUCESSFUL! ****'/
     *  5X,'RESTARTING WITH NEW RANDOM NUMBER AND ROTATION SEQUENCE')
 9060 FORMAT(/1X,'LOCALIZATION FAILED --- LOCALIZATION ABORTED')
 9070 FORMAT(10X,'LOCALIZATION CONVERGED IN',I4,' ITERATIONS')
 9010 FORMAT(/10X,'THE INITIAL LOCALIZATION SUM IS',F14.6,' DEBYE**2')
 9080 FORMAT(10X,'  THE FINAL LOCALIZATION SUM IS',F14.6,' DEBYE**2')
 9090 FORMAT(/10X,'THE LOCALIZATION TRANSFORMATION MATRIX IS')
 9100 FORMAT(/10X,'THE BOYS LOCALIZED ORBITALS ARE')
 9110 FORMAT(/4X,'CENTROIDS OF LOCALIZED ORBITALS (BOHR)'/
     *        5X,2HMO,9X,1HX,8X,1HY,8X,1HZ)
 9120 FORMAT(3X,I4,5X,3(F8.4,1X))
      END
C*MODULE LOCAL   *DECK LMOFRZ
      SUBROUTINE LMOFRZ(NOUT,IPASS,MCORE,MDOC,NUMLOC,NOROT,LABMO,L1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXAO=8192)
C
      LOGICAL IDOC,IACT,JDOC,JACT,GOPARR,DSKWRK,MASWRK,ABEL,
     *        CANONC,FCORE,FORS,NOCI,EKT,LINSER,FDIRCT,QCORR,
     *        ORIENT,ORMFUL
C
      DIMENSION NOROT(NUMLOC,NUMLOC),LABMO(NUMLOC)
C
      PARAMETER (MXNORO=250)
C
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIMFCC(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NORTMC(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /ORNTMO/ ORIENT,ORMFUL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMQMT/ IRPLAB(14),IRPNUM(14),IRPDIM(14),IRPDEG(14)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RMC,ORMAS/8HMCSCF   ,8HORMAS   /
C
      DO 120 I=1,NUMLOC
         DO 110 J=1,NUMLOC
            NOROT(I,J) = 0
  110    CONTINUE
  120 CONTINUE
C
C     BLOCK OUT ANY ROTATIONS INVOLVING THE FROZEN ORBITALS
C
      DO 170 K=1,NOUT
         IF(IPASS.EQ.1) THEN
            KK = MOOUTA(K)
         ELSE
            KK = MOOUTB(K)
         END IF
         DO 160 I=1,NUMLOC
            NOROT(I,KK) = 1
            NOROT(KK,I) = 1
  160    CONTINUE
  170 CONTINUE
C
C     LOCALIZE AS MUCH AS POSSIBLE, LEAVING WAVEFUNCTION INVARIANT.
C
      NDBOCC = MCORE+MDOC
      DO 220 I=1,NUMLOC
         IDOC  = I.LE.NDBOCC
         IACT  = I.GT.NDBOCC
         DO 210 J=1,NUMLOC
            JDOC  = J.LE.NDBOCC
            JACT  = J.GT.NDBOCC
            IF(IDOC .AND. JACT) NOROT(I,J)=1
            IF(IACT .AND. JDOC) NOROT(I,J)=1
  210    CONTINUE
  220 CONTINUE
C
C         ORMAS-MCSCF IS NOT INVARIANT TO ROTATIONS THAT MIX
C         ORBITALS BETWEEN THE MULTIPLE ACTIVE SPACES, SO UNLESS
C         THE USER SPECIFICALLY WANTS TO, FREEZE THESE ROTATIONS.
C         STORAGE -LABMO- USED TO DEFINE ORBITAL SPACE COUNTER
C
      IF(.NOT.(SCFTYP.EQ.RMC  .AND.  CISTEP.EQ.ORMAS)) GO TO 500
      IF(ORMFUL) GO TO 500
C
      DO I=1,NDBOCC
         LABMO(I) = 1
      ENDDO
      DO ISPACE=1,NSPACE
         DO I=MSTA(ISPACE),NUMLOC
            LABMO(I) = ISPACE+1
         ENDDO
      ENDDO
C
      DO I=1,NUMLOC
         ISPACE = LABMO(I)
         DO J=1,I
            IF(ISPACE.NE.LABMO(J)) THEN
               NOROT(I,J)=1
               NOROT(J,I)=1
            END IF
         ENDDO
      ENDDO
C
C        OPTIONALLY, RESTRICT ROTATIONS SO THAT ORBITALS WITH
C        DIFFERENT SYMMETRIES ARE NOT MIXED.
C
  500 CONTINUE
      IF(ISYMLO.EQ.0) RETURN
C
      IF(ILABMO.EQ.0) THEN
         IF (MASWRK) WRITE(IW,9010)
         RETURN
      ELSE
         IF(IPASS.EQ.1) CALL DAREAD(IDAF,IODA,LABMO,L1,255,1)
         IF(IPASS.EQ.2) CALL DAREAD(IDAF,IODA,LABMO,L1,256,1)
      END IF
C
      IDEG=0
      DO 390 I=1,NUMLOC
         ISYM=LABMO(I)
         IDEG=IDEG+1
         DO 310 IRP=1,NIRRED
            IF(ISYM.EQ.IRPLAB(IRP)) IRREP = IRP
  310    CONTINUE
         IDIM = IRPDIM(IRREP)
         IF(IDIM.EQ.1) IDEG=0
C
         JDEG = 0
         DO 380 J=1,I
            JSYM = LABMO(J)
            JDEG = JDEG + 1
            DO 330 IRP=1,NIRRED
               IF(JSYM.EQ.IRPLAB(IRP)) IRREP = IRP
  330       CONTINUE
            JDIM = IRPDIM(IRREP)
            IF(JDIM.EQ.1) JDEG=0
            IF(ISYM.NE.JSYM  .OR.  IDEG.NE.JDEG) THEN
               NOROT(I,J)=1
               NOROT(J,I)=1
            END IF
            IF(JDEG.EQ.JDIM) JDEG=0
  380    CONTINUE
         IF(IDEG.EQ.IDIM) IDEG=0
  390 CONTINUE
C
      RETURN
 9010 FORMAT(//1X,55(1H*)/
     *       1X,'YOUR GROUP DOES NOT SUPPORT ORBITAL SYMMETRY LABELS.'/
     *       1X,'-SYMLOC- OPTION IS DEACTIVATED.  MAXIMUM LOCALIZATION'/
     *       1X,'CONSISTENT WITH WAVEFUNCTION INVARIANCE WILL BE DONE.'/
     *       1X,55(1H*)//)
      END
C*MODULE LOCAL   *DECK LMOINP
      SUBROUTINE LMOINP
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION LOCAL,MAKEFP,METHOD
C
      LOGICAL PRTLOC,FCORE,SYMLOC,GOPARR,DSKWRK,MASWRK
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,POLDCM,POLANG
      LOGICAL POLNUM,POLAPP,KMIDPT,POLDYN,ORIENT,ORMFUL
      LOGICAL CANONC,FCORMC,FORS,NOCI,EKT,LINSER
C
      PARAMETER (MXATM=2000, NMO=500, MXAO=8192)
      PARAMETER (MXMID=50, MXNORO=250)
      PARAMETER (ONE=1.0D+00)
C
      DIMENSION MOINA(50),MOINB(50),TYPES(3)
C
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      COMMON /EDCMP / ZIJ(NMO),ZMO(5,NMO),OCCUP(NMO),DPFREQ(50),
     *                MOIDNO(5,NMO),IJMO(2,NMO),MOIJ(NMO),NMOIJ(NMO),
     *                NMOAT(NMO),NDPFREQ,IPROT(5),NPROT,
     *                MOIDON,EDCOMP,DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,
     *                POLANG,POLAPP,KMIDPT,POLDYN
      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 /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORMC,FORS,NOCI,EKT,LINSER
      COMMON /NOPOL / NOPATM(MXATM),KPOINT(MXMID),NODENS(MXATM)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /ORNTMO/ ORIENT,ORMFUL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /POLNMP/ POLNUM
      COMMON /PRPOPT/ ILOCAL
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
C     ----- SET UP NAMELIST SIMULATION -----
C
      PARAMETER (NNAM=37)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA LOCAL/8HLOCAL   /
      DATA QNAM/8HNOUTA   ,8HNOUTB   ,8HMOOUTA  ,8HMOOUTB  ,
     *          8HNINA    ,8HNINB    ,8HMOINA   ,8HMOINB   ,
     *          8HCVGLOC  ,8HMAXLOC  ,8HPRTLOC  ,8HFCORE   ,
     *          8HSYMLOC  ,8HIJMO    ,8HZIJ     ,8HMOIDON  ,
     *          8HEDCOMP  ,8HDIPDCM  ,8HMOIJ    ,8HNMOIJ   ,
     *          8HIPROT   ,8HDEPRNT  ,8HQADDCM  ,8HZDO     ,
     *          8HPOLDCM  ,8HPOLANG  ,8HPOLNUM  ,8HPOLAPP  ,
     *          8HNOPATM  ,8HKPOINT  ,8HKMIDPT  ,8HORIENT  ,
     *          8HNODENS  ,8HORMFUL  ,8HPOLDYN  ,8HDPFREQ  ,
     *          8HNDPFRQ  /
      DATA KQNAM/1,1,-1,-1,1,1,501,501,3,1,0,0,0,10001,5003,0,0,0,
     *           5001,5001,51,0,0,0,0,0,0,0,1001,101,0,0,5001,0,0,
     *           503,1/
C
      DATA TYPES/8HBOYS    ,8HRUEDNBRG,8HPOP     /
      DATA GUGA,GVB,UHF/8HGUGA    ,8HGVB     ,8HUHF     /
      DATA RHF,ROHF,RMC/8HRHF     ,8HROHF    ,8HMCSCF   /
      DATA ORMAS/8HORMAS   /
      DATA MAKEFP/8HMAKEFP  /
C
C     ----- READ IN LOCALIZATION INPUT PARAMETERS -----
C
      IF(ILOCAL.EQ.0 .AND. NBUFMO.EQ.0) RETURN
C
      IF(CITYP.EQ.GUGA  .OR.  SCFTYP.EQ.GVB) THEN
         IF (MASWRK) WRITE(IW,9000) SCFTYP
         CALL ABRT
      END IF
C
      KQNAM(3)=10*MXAO + 1
      KQNAM(4)=10*MXAO + 1
      NOUTA = 0
      NOUTB = 0
      DO 100 I = 1,NUM
         MOOUTA(I) = 0
         MOOUTB(I) = 0
  100 CONTINUE
      NINA = 0
      NINB = 0
      DO 110 I=1,50
         MOINA(I) = 0
         MOINB(I) = 0
  110 CONTINUE
      CVGLOC= 1.0D-06
      MAXLOC= 250
      PRTLOC= .FALSE.
      FCORE = .TRUE.
      SYMLOC= .FALSE.
      DO 120 I = 1,NMO
         IJMO(1,I) = 0
         IJMO(2,I) = 0
         ZIJ(I) = -ONE
         MOIJ(I) = 0
         NMOIJ(I) = -1
  120 CONTINUE
      MOIDON = .FALSE.
      EDCOMP = .FALSE.
      DIPDCM = .FALSE.
      QADDCM = .FALSE.
      DO 130 I = 1,5
         IPROT(I) = 0
  130 CONTINUE
      NPROT=0
      DEPRNT = .FALSE.
      ZDO = .FALSE.
C
      POLDCM=.FALSE.
      POLDYN=.FALSE.
      POLANG=.TRUE.
      POLNUM=.FALSE.
      POLAPP=.FALSE.
      IF (RUNTYP.EQ.MAKEFP) POLDCM = .TRUE.
      IF (RUNTYP.EQ.MAKEFP) POLDYN = .TRUE.
C
      KMIDPT=.TRUE.
      ORIENT=.FALSE.
      ORMFUL=.FALSE.
C
      DO I = 1,MXATM
         NOPATM(I) = 0
      ENDDO
      DO I = 1,MXMID
         KPOINT(I) = 0
      ENDDO
      DO I = 1,MXATM
         NODENS(I) = 0
      ENDDO
C
      IF (RUNTYP.EQ.MAKEFP) THEN
         NDPFREQ=12
         DPFREQ(1) = 0.00279164292546976127D+00
         DPFREQ(2) = 0.015106644820649056973D+00
         DPFREQ(3) = 0.03900169130663191517D+00
         DPFREQ(4) = 0.077996102406373782231D+00
         DPFREQ(5) = 0.13865052123704475618D+00
         DPFREQ(6) = 0.233222703362246034002D+00
         DPFREQ(7) = 0.38589725057860362524D+00
         DPFREQ(8) = 0.649114040084500804964D+00
         DPFREQ(9) = 1.15390381343780158400D+00
         DPFREQ(10)= 2.307592234716658147562D+00
         DPFREQ(11)= 5.95764321386442380950D+00
         DPFREQ(12)= 32.239080141259587698936D+00
      ELSE
         NDPFREQ=1
         DPFREQ(1) = 0.0D+00
      END IF
C
C     ----- READ IN $LOCAL PARAMETERS -----
C
      JRET=0
      CALL NAMEIO(IR,JRET,LOCAL,NNAM,QNAM,KQNAM,
     *            NOUTA,NOUTB,MOOUTA,MOOUTB,NINA,NINB,MOINA,MOINB,
     *            CVGLOC,MAXLOC,PRTLOC,FCORE,SYMLOC,IJMO,ZIJ,MOIDON,
     *            EDCOMP,DIPDCM,MOIJ,NMOIJ,IPROT,DEPRNT,QADDCM,
     *            ZDO,POLDCM,POLANG,POLNUM,POLAPP,NOPATM,KPOINT,KMIDPT,
     *            ORIENT,NODENS,ORMFUL,POLDYN,DPFREQ,NDPFREQ,
     *   0,0,  0,0,0,0,0,
     *   0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
C
      IF (JRET .EQ. 2) THEN
         IF (MASWRK) WRITE (IW,9010)
         CALL ABRT
      END IF
C
C     ----- THE ENERGY DECOMPOSITION IS ONLY IMPLEMENTED FOR RHF/ROHF
C
      IF (EDCOMP .AND. (SCFTYP.NE.RHF .AND. SCFTYP.NE.ROHF)) THEN
         IF (MASWRK) WRITE(IW,9015)
         CALL ABRT
      END IF
C
C           ORIENTATION APPLIES ONLY TO MCSCF ACTIVE ORBITALS
C           AND NEEDS THE INTEGRALS OF A RUEDENBERG LOCALIZATION
C
      IF(SCFTYP.NE.RMC  .OR.  ILOCAL.NE.2) ORIENT=.FALSE.
C
      NERR=0
      IF(SCFTYP.NE.UHF) THEN
         NOUTB=0
         NINB=0
      END IF
C
      IPRTLO=0
      IF(PRTLOC .OR. EDCOMP) IPRTLO=1
      ISYMLO=0
      IF(SYMLOC) ISYMLO=1
      IFCORE=0
      IF (EDCOMP) FCORE = .FALSE.
      IF(FCORE)  IFCORE=1
      DO 200 I = 1,5
         IF (IPROT(I) .NE. 0) NPROT = NPROT + 1
  200 CONTINUE
C
C        GENERATE NOUTA/NOUTB DATA FOR FROZEN CORES ONLY
C
      IF(FCORE) THEN
         IF(NINA .GT.0) THEN
            IF (MASWRK) WRITE(IW,9060) 'NINA '
            NERR=NERR+1
         END IF
         IF(NOUTA.GT.0) THEN
            IF (MASWRK) WRITE(IW,9060) 'NOUTA'
            NERR=NERR+1
         END IF
         CALL LMOPSI(1,MCORE,MDOC,MACT,NUMLOC)
         NOUTA = MCORE
         DO 210 I=1,MCORE
            MOOUTA(I) = I
  210    CONTINUE
         IF(SCFTYP.EQ.UHF) THEN
            IF(NINB .GT.0) THEN
               IF (MASWRK) WRITE(IW,9060) 'NINB '
               NERR=NERR+1
            END IF
            IF(NOUTB.GT.0) THEN
               IF (MASWRK) WRITE(IW,9060) 'NOUTB'
               NERR=NERR+1
            END IF
            CALL LMOPSI(2,MCORE,MDOC,MACT,NUMLOC)
            NOUTB = MCORE
            DO 220 I=1,MCORE
               MOOUTB(I) = I
  220       CONTINUE
         END IF
      END IF
      IF(NERR.GT.0) CALL ABRT
C
C        CONVERT NINA/NINB TO EQUIVALENT NOUTA/NOUTB DATA
C
      IF(NINA.GT.0) THEN
         IF(NOUTA.GT.0) THEN
            IF (MASWRK) WRITE(IW,9050) 'NINA','NOUTA'
            NERR=NERR+1
         END IF
         CALL LMOPSI(1,MCORE,MDOC,MACT,NUMLOC)
         NOUTA = NUMLOC-NINA
         MO = 0
         DO 320 I=1,NUMLOC
            DO 310 J=1,NINA
               IF(I.EQ.MOINA(J)) GO TO 320
  310       CONTINUE
            MO = MO+1
            MOOUTA(MO)=I
  320    CONTINUE
      END IF
      IF(NINB.GT.0) THEN
         IF(NOUTB.GT.0) THEN
            IF (MASWRK) WRITE(IW,9050) 'NINB','NOUTB'
            NERR=NERR+1
         END IF
         CALL LMOPSI(2,MCORE,MDOC,MACT,NUMLOC)
         NOUTB = NUMLOC-NINB
         MO = 0
         DO 360 I=1,NUMLOC
            DO 350 J=1,NINB
               IF(I.EQ.MOINB(J)) GO TO 360
  350       CONTINUE
            MO = MO+1
            MOOUTB(MO)=I
  360    CONTINUE
      END IF
      IF(NERR.GT.0) CALL ABRT
C
C        AT THIS POINT ALL FREEZING DATA, EXCEPT SYMLOC,
C        IS INCORPORATED INTO NOUTA/NOUTB FORM.
C
      IF (MASWRK .AND. NPRINT.NE.-23) THEN
         WRITE(IW,9020) TYPES(ILOCAL),NOUTA,PRTLOC,NOUTB,MAXLOC,CVGLOC,
     *                  FCORE,SYMLOC,MOIDON,EDCOMP,DIPDCM,NPROT,QADDCM,
     *                  DEPRNT,POLNUM,ORIENT,POLDYN
         IF(SCFTYP.EQ.RMC  .AND.  CISTEP.EQ.ORMAS) WRITE(IW,9025) ORMFUL
         IF(NOUTA.NE.0) WRITE (IW,9030) (MOOUTA(I),I = 1,NOUTA)
         IF(NOUTB.NE.0) WRITE (IW,9040) (MOOUTB(I),I = 1,NOUTB)
         IF(NPROT.NE.0) WRITE (IW,9045) (IPROT(I),I=1,NPROT)
      END IF
      RETURN
C
 9000 FORMAT(/1X,'ERROR! LOCALIZATION DOES NOT WORK WITH SCFTYP=',A8)
 9010 FORMAT(/1X,'ERROR! YOUR $LOCAL INPUT IS INCORRECT')
 9015 FORMAT(/1X,'ERROR! EDCOMP ONLY WORKS FOR SCFTYP=RHF/ROHF')
 9020 FORMAT(/10X,29(1H-)/
     *       10X,'LOCALIZATION INPUT PARAMETERS'/10X,29(1H-)/
     *       10X,' LOCAL=',A8,5X,' NOUTA=',I8/
     *       10X,'PRTLOC=',L8,5X,' NOUTB=',I8/
     *       10X,'MAXLOC=',I8,5X,'CVGLOC=',1P,E8.2,0P/
     *       10X,' FCORE=',L8,5X,'SYMLOC=',L8/
     *       10X,'MOIDON=',L8,5X,'EDCOMP=',L8/
     *       10X,'DIPDCM=',L8,5X,' NPROT=',I8/
     *       10X,'QADDCM=',L8,5X,'DEPRNT=',L8/
     *       10X,'POLNUM=',L8,5X,'ORIENT=',L8/
     *       10X,'POLDYN=',L8)
 9025 FORMAT(10X,'ORMAS-MCSCF FULL LOCALIZATION FLAG ORMFUL=',L1)
 9030 FORMAT(10X,'MOOUTA =',2X,12I5,30(/20X,12I5))
 9040 FORMAT(10X,'MOOUTB =',2X,12I5,50(/20X,12I5))
 9045 FORMAT(10X,'IPROT  =',2X,5I5)
 9050 FORMAT(1X,'ERROR! ',A4,' AND ',A5,' ARE MUTUALLY EXCLUSIVE')
 9060 FORMAT(1X,'ERROR!  FCORE IS INCOMPATIBLE WITH',A5)
      END
C*MODULE LOCAL   *DECK LMOLOC
      SUBROUTINE LMOLOC(ILOCAL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION SPIN(2)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /FMCOM / X(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 /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA UHF/8HUHF     /, SPIN/8H ALPHA  ,8H* BETA  /
      DATA CHECK/8HCHECK   /
C
C     ----- MEMORY DRIVER FOR BOYS/POPULATION LOCALIZATION -----
C
      IF (NPRINT.NE.-23) THEN
         IF((ILOCAL.EQ.1) .AND. MASWRK) WRITE(IW,9000)
         IF((ILOCAL.EQ.3) .AND. MASWRK) WRITE(IW,9100)
      END IF
C
      IPASS=0
  200 CONTINUE
      IPASS=IPASS+1
      IF(IPASS.GT.2) GO TO 700
      IF((SCFTYP.EQ.UHF) .AND. MASWRK) WRITE(IW,9010) SPIN(IPASS)
C
      CALL LMOPSI(IPASS,MCORE,MDOC,MACT,NUMLOC)
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      N1 = NUMLOC
      N2 = (N1*N1+N1)/2
      N3 = N1*N1
C
C     PARTITION FAST MEMORY
C
      CALL VALFM(LOADFM)
      LVIN  = 1     + LOADFM
      LTRAN = LVIN  + L1*L1
      LVOUT = LTRAN + N1*N1
      LTRI  = LVOUT + L1*L1
      LWRK  = LTRI  + L2
      LIORD = LWRK  + L1
      LIIR  = LIORD + N1
      LIA   = LIIR  + L1
      LRIJ  = LIA   + L1
      IF(ILOCAL.EQ.1) THEN
         LAST  = LRIJ  + N2*3
      ELSE
         LAST  = LRIJ  + N2*NAT
      END IF
      LQPIX = LAST
      LQPJX = LQPIX + L1
      LNOROT= LQPJX + L1
      LAST  = LNOROT+ N3
      NEED = LAST - LOADFM -1
      IF (MASWRK .AND. NPRINT.NE.-23) WRITE(IW,9020) NEED
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 600
C
C        READ IN THE ORBITALS TO BE LOCALIZED
C
      IF(IPASS.EQ.1) THEN
         NOUT=NOUTA
         CALL DAREAD(IDAF,IODA,X(LVIN),L3,15,0)
      ELSE
         NOUT=NOUTB
         CALL DAREAD(IDAF,IODA,X(LVIN),L3,19,0)
      END IF
C
      CALL LMOFRZ(NOUT,IPASS,MCORE,MDOC,N1,X(LNOROT),X(LIIR),L1)
C
C        CARRY OUT THE LOCALIZATION
C
      IF (MASWRK .AND. NPRINT.NE.-23) WRITE(IW,9030) NUMLOC,NOUT
      IF (NUMLOC-NOUT .LE.0) THEN
         IF(MASWRK) WRITE(IW,9040)
         CALL ABRT
      END IF
C
      IF(ILOCAL.EQ.1) THEN
         CALL BOYLOC(IPASS,X(LVIN),X(LTRAN),X(LVOUT),
     *               X(LIORD),X(LIIR),X(LIA),X(LRIJ),
     *               X(LQPIX),X(LQPJX),X(LNOROT),X(LTRI),X(LWRK),
     *               L1,L2,N1,N2,NOUT)
      END IF
      IF(ILOCAL.EQ.3) THEN
         LSAO=LTRI
         LMAP=LIA
         CALL POPLOC(IPASS,X(LVIN),X(LTRAN),X(LVOUT),
     *               X(LIORD),X(LIIR),X(LMAP),X(LRIJ),
     *               X(LQPIX),X(LQPJX),X(LNOROT),X(LSAO),
     *               L1,L2,N1,N2,NAT,NOUT)
      END IF
C
      IF(IPASS.EQ.1) THEN
         CALL DAWRIT(IDAF,IODA,X(LVOUT),L3,71,0)
      ELSE
         CALL DAWRIT(IDAF,IODA,X(LVOUT),L3,72,0)
      END IF
C
  600 CONTINUE
      CALL RETFM(NEED)
      IF(SCFTYP.EQ.UHF) GO TO 200
C
  700 CONTINUE
C
      IF (MASWRK .AND. NPRINT.NE.-23)
     * WRITE(IW, FMT='('' ...... END OF ORBITAL LOCALIZATION ......'')')
      CALL TIMIT(1)
      RETURN
 9000 FORMAT(/10X,25(1H-)/10X,'BOYS ORBITAL LOCALIZATION'/10X,25(1H-)/
     *        10X,'S.F.BOYS, IN "QUANTUM THEORY OF ATOMS, MOLECULES,',
     *            ' AND SOLIDS"'/
     *        10X,'P.O.LOWDIN, ED. ACADEMIC PRESS, NY, 1966 PP253-266.')
 9100 FORMAT(/10X,53(1H-)/,10X,
     * 'MOLECULAR ORBITALS LOCALIZED BY THE POPULATION METHOD'/
     *  10X,53(1H-)/10X,'REFERENCE: ',
     * 'J. PIPEK AND P. G. MEZEY, J. CHEM. PHYS. 90, 4916 (1989)')
 9010 FORMAT(/10X,'*****',A6,' ORBITAL LOCALIZATION *****')
 9020 FORMAT(/10X,'THIS LOCALIZATION REQUIRES',I10,' WORDS OF MEMORY.')
 9030 FORMAT(10X,'THIS LOCALIZATION HAS',I4,' ORBITALS, OF WHICH',
     *        I4,' ARE FROZEN.')
 9040 FORMAT(/1X,'*** ERROR ***'/
     *       1X,'THIS LOCALIZATION HAS NO ORBITALS TO LOCALIZE.'/
     *       1X,'PLEASE RECONSIDER YOUR FCORE/MOIN/MOOUT INPUT.')
      END
C*MODULE LOCAL   *DECK LMOPSI
      SUBROUTINE LMOPSI(IPASS,MCORE,MDOC,MACT,MOCC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      LOGICAL CANONC,FCORE,FORS,NOCI,EKT,LINSER
C
      PARAMETER (MXATM=2000, MXRT=100, MXNORO=250)
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACTDT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
      DATA RMC/8HMCSCF   /
      DATA ALDET,ORMAS,GUGA,GENCI
     *    /8HALDET   ,8HORMAS   ,8HGUGA    ,8HGENCI   /
C
C     ----- SET UP DEFAULTS FOR THE WAVEFUNCTION -----
C     -MCORE- IS THE NUMBER OF CHEMICAL CORE (NONVALENCE) ORBITALS,
C     -MDOC- IS THE NUMBER OF DOUBLY OCCUPIED VALENCE ORBITALS,
C     -MACT- IS THE NUMBER OF PARTIALLY OCCUPIED ORBITALS.
C     FOR UHF, -MCORE- IS STILL THE CHEMICAL CORE, BUT -MDOC- IS THE
C     NUMBER OF SINGLY OCCUPIED VALENCE ORBITALS, OF ALPHA OR BETA
C     SPIN DEPENDING ON THE VALUE FOR -IPASS-.
C
      MCORE= 0
      MDOC = 0
      MACT = 0
C
C     FIRST COUNT THE CHEMICAL CORE ORBITALS.
C
      MCORE = NUMCOR()
C
C     ASSIGN VALENCE ORBITALS ACCORDING TO THE SCF TYPE.
C
      IF(SCFTYP.EQ.RHF) THEN
         MDOC = NA - MCORE
         MACT = 0
      END IF
      IF(SCFTYP.EQ.UHF) THEN
         IF(IPASS.EQ.1) THEN
            MDOC = NA - MCORE
         ELSE
            MDOC = NB - MCORE
         END IF
         MACT = 0
      END IF
      IF(SCFTYP.EQ.ROHF) THEN
         MDOC = NB - MCORE
         MACT = NA - NB
      END IF
      IF(SCFTYP.EQ.RMC) THEN
         IF(CISTEP.EQ.ALDET   .OR.  CISTEP.EQ.ORMAS
     *                        .OR.  CISTEP.EQ.GENCI) THEN
            MDOC = NCORSV - MCORE
            MACT = NACTDT
         END IF
         IF(CISTEP.EQ.GUGA) THEN
            MDOC = NMCC - MCORE
            MACT = NDOC + NALP + NAOS + NBOS + NVAL
         END IF
      END IF
C
      MOCC = MCORE + MDOC + MACT
      RETURN
      END
C*MODULE LOCAL   *DECK LMORND
      DOUBLE PRECISION FUNCTION LMORND(XX,D,L1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NBASIS,NQMT,NE,NAE,NBE,
     *                ZAN(MXATM),COORD(3,MXATM),IAN(MXATM)
C
      DIMENSION D(L1,L1),U(1)
C
      SAVE U
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      PI = ACOS(-ONE)
      IF (XX .EQ. ZERO) GO TO 100
         N = ABS(NAE-NBASIS)+1
         M = N+5
         XY = D(N,M)*ATAN(ONE)
         U(1) = (PI+XY)**5
         IU1 = INT(U(1))
         XY = IU1
         U(1) = U(1)-XY
         LMORND = U(1)
         RETURN
C
  100 CONTINUE
      U(1) = (PI+U(1))**5
      IU1 = INT(U(1))
      XY = IU1
      U(1) = U(1)-XY
      LMORND = U(1)
      RETURN
      END
C*MODULE LOCAL   *DECK LMOX
      SUBROUTINE LMOX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (NMO=500, MXMOL=10)
C
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,POLDCM,POLANG,
     *        POLNUM,POLAPP,KMIDPT,POLDYN,ORIENT,ORMFUL
C
      COMMON /EDCMP / ZIJ(NMO),ZMO(5,NMO),OCCUP(NMO),DPFREQ(50),
     *                MOIDNO(5,NMO),IJMO(2,NMO),MOIJ(NMO),NMOIJ(NMO),
     *                NMOAT(NMO),NDPFREQ,IPROT(5),NPROT,
     *                MOIDON,EDCOMP,DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,
     *                POLANG,POLAPP,KMIDPT,POLDYN
      COMMON /MOINFO/ NMORO,NMOL,MONO(MXMOL),NSHLLS(MXMOL),NBAS(MXMOL),
     *                NEM(MXMOL),ICHM(MXMOL),NLOC(MXMOL),IEQM(MXMOL),
     *                IBSSE,ICTPPL,ICTPLX,IMOROX,IRVS,IDUM
      COMMON /ORNTMO/ ORIENT,ORMFUL
      COMMON /POLNMP/ POLNUM
      COMMON /PRPOPT/ ILOCAL
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA AMOROK /8HMOROKUMA/
      DATA RHF/8HRHF     /
C
      IF (ILOCAL.EQ.0) RETURN
      IF (ILOCAL.EQ.1) CALL LMOLOC(ILOCAL)
      IF (ILOCAL.EQ.2) CALL LOCENG
      IF (ILOCAL.EQ.3) CALL LMOLOC(ILOCAL)
C
      IF (RUNTYP.EQ.AMOROK .AND. NMORO.LE.2) RETURN
C
C     ----- ASSIGN ORBITALS TO ATOM(S) -----
C
      IF (MOIDON)      CALL MOIDM
C
C     ----- ORBITAL ENERGY DECOMPOSITION -----
C
      IF (EDCOMP)      CALL EDCMPM
C
C     ----- MOMENTS LMO-DECOMPOSITION -----
C
      IF (DIPDCM)      CALL DIPLMM
      IF (QADDCM)      CALL QUADMM
C
C     ----- ALPHA POLARIZABILITY LMO-DECOMPOSITION -----
C
      IF(POLDYN) POLNUM=.FALSE.
      IF(POLDCM.OR.POLDYN) THEN
         IF(SCFTYP.EQ.RHF) THEN
            CALL MOIDM
            CALL DIPLMM
         END IF
         IF(POLNUM) THEN
            CALL POLNMM
         ELSE
            CALL LPOLX
         END IF
      END IF
C
C     ----- ORIENTATION OF ORBITALS -----
C     ONLY IF EDMISTON-RUEDENBERG LOCALIZATION PERFORMED
C     AS THIS IS THE ONLY METHOD THAT SAVES THE DENSITY.
C
      IF (ORIENT.AND.ILOCAL.EQ.2) CALL DIRLMO
C
      RETURN
      END
C*MODULE LOCAL   *DECK LOCANA
      SUBROUTINE LOCANA(TLOC,TRI,TRILOC,WRK,M1,M2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION TLOC(M1,M1),TRI(M2),TRILOC(M2),WRK(M1)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RHF,UHF,AMCSCF/8HRHF     ,8HUHF     ,8HMCSCF   /
      DATA ROHF,GVB/8HROHF    ,8HGVB     /
C
C     ----- ANALYSIS ROUTINE FOR RUEDENBERG LMOS -----
C     FOR RHF AND UHF, ROTATE THE FOCK MATRIX TO THE LOCALIZED ORBITAL
C     BASIS, WHILE FOR MCSCF ROTATE THE DENSITY MATRIX.  NOTE THAT
C     ONLY THE LOCALIZED ORBITAL SPACE IS PRESENT IN THE -TRI- MATRIX
C     HOLDING THE FOCK/DENSITY, AND IN THE TRANSFORMATION MATRIX -TLOC-.
C
      IF(SCFTYP.EQ.ROHF) RETURN
      IF(SCFTYP.EQ.GVB)  RETURN
      IF(.NOT.MASWRK) RETURN
C
      IF(SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF) WRITE(IW,9010)
      IF(SCFTYP.EQ.AMCSCF) WRITE(IW,9020)
      CALL TFTRI(TRILOC,TRI,TLOC,WRK,M1,M1,M1)
      CALL PRTRI(TRILOC,M1)
C
C          WRITE OUT THE LOCALIZED DENSITY TO DAF RECORD 285
C
      CALL DAWRIT(IDAF,IODA,TRILOC,M2,285,0)
      RETURN
C
 9010 FORMAT(1X,'FOCK OPERATOR FOR THE LOCALIZED ORBITALS IS')
 9020 FORMAT(1X,'DENSITY MATRIX FOR THE LOCALIZED ORBITALS IS')
      END
C*MODULE LOCAL   *DECK LOCCHK
      SUBROUTINE LOCCHK(IW,TWOEI,IA,M1,M2,M4)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION TWOEI(M4),IA(M2)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- COMPUTE SUBSUMS OF INTEGRALS DURING LOCALIZATION -----
C     M1    NUMBER OF ORBITALS BEING LOCALIZED
C     TWOEI ARRAY OF INTEGRALS STORED WITH TOTAL SYMMETRY
C
C     D =  LOCALIZATION SUM   GOES TO A MAXIMUM
C     C =  COULOMB SUM        INVARIANT
C     X =  EXCHANGE SUM       INVARIANT
C
      C=ZERO
      X=ZERO
      D=ZERO
      IT=0
      DO 110 I=1,M1
         II=IT + I
         IIT=IA(II)
         IIII=IIT + II
         D=D + TWOEI(IIII)
         JJ=0
         DO 100 J=1,I
            JJ=JJ+J
            IF(I.EQ.J) GO TO 100
               IJ=IT + J
               IJIJ=IA(IJ) + IJ
               X=X + TWOEI(IJIJ)
               IIJJ=IIT+JJ
               C=C + TWOEI(IIJJ)
  100    CONTINUE
         IT=IT+I
  110 CONTINUE
      C=D + C + C
      X=D + X + X
      IF (MASWRK) WRITE(IW,9000) D,C,X
      RETURN
C
 9000 FORMAT(5X,'DIAGONAL SUM D=',F20.10/
     *       5X,' COULOMB SUM C=',F20.10/
     *       5X,'EXCHANGE SUM X=',F20.10)
      END
C*MODULE LOCAL   *DECK LOCENG
      SUBROUTINE LOCENG
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL PACK2E,PK,PANDK,BLOCK,
     *        GOPARR,DSKWRK,MASWRK,SVDSKW,SVGPAR,
     *        ABEL,ABELPT,DIRSCF,FDIFF,DIRTRF,SVDIRT,
     *        RDINT,DOEXCH
C
      DIMENSION SPIN(2)
C
      PARAMETER (MXATM=2000, MXSH=5000, MXAO=8192)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA UHF/8HUHF     /, SPIN/6H ALPHA,6H* BETA/
      DATA RMC/8HMCSCF   /
      DATA CHECK/8HCHECK   /, HSSIAN,PROP/8HHESSIAN ,8HPROP    /
C
C     ----- MAIN DRIVER FOR ENERGY LOCALIZATION -----
C     THE CODE ASSOCIATED WITH ENERGY LOCALIZATION IS TAKEN
C     PROGRAM 'LOCL' FROM THE 'ALIS' (AMES LAB, IOWA STATE)
C     PROGRAM SYSTEM, IN A 1984 VERSION.
C
C     MANY PEOPLE FROM THE RUEDENBERG GROUP PARTICIPATED
C     IN THE DEVELOPMENT OF THIS PROGRAM.  THE MAIN PORTION
C     OF THIS PROGRAM (THE PRESENT ROUTINE LOCLIZ) WAS WRITTEN
C     BY RICH RAFFENETTI BETWEEN 1966-1968.  OF COURSE, MANY
C     OTHERS IN THE RUEDENBERG GROUP, SUCH AS CLYDE EDMISTON,
C     KEN SUNDBERG, STEVE ELBERT, ... ALSO CONTRIBUTED.
C
C     THE ENERGY LOCALIZATION CODE WAS ADAPTED TO GAMESS
C     BY MIKE SCHMIDT AT NORTH DAKOTA STATE UNIVERSITY
C     IN MARCH OF 1991, BY PERMISSION OF KLAUS RUEDENBERG.
C
C     TWO PAPERS GERMANE TO THIS METHOD AND ITS IMPLEMENTATION ARE
C     C.EDMISTON, K.RUEDENBERG, REV.MOD.PHYS. 1963, 35, 457-465
C     R.C.RAFFENETTI, K.RUEDENBERG, C.L.JANSSEN, H.F.SCHAEFER,
C     THEORET.CHIM.ACTA 1993, 86, 149-165.
C
      IF (MASWRK) WRITE(IW,9000)
      SVDIRT = DIRTRF
C
C     ----- CHECK TO SEE IF THE TRANSFORMED TWO-ELECTRON INTEGRALS
C           ARE GIVEN -----
C
      RDINT = .FALSE.
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $TWOEI ',IEOF)
      IF (IEOF.EQ.0) RDINT = .TRUE.
      IF (RDINT) GO TO 100
C
C     ----- PREPARE CORRECT INTEGRAL LIST ---
C     THIS MUST BE -J- FORMAT
C     IN SOME CASES, THIS FILE MAY BE ALREADY THERE.
C
      IF(PK) THEN
         IF(MASWRK) WRITE(IW,*) 'LOCENG: SOMEHOW THIS IS A -PK- FILE'
         CALL ABRT
      END IF
      IF(SCFTYP.EQ.RMC  .AND.  RUNTYP.NE.PROP) GO TO 100
      IF(RUNTYP.EQ.HSSIAN) GO TO 100
C
C     IF DUPLICATED INTEGRAL FILES DESIRED FOR TRANSF., GENERATE THEM.
C     IF POINT GROUP IS NON-ABELIAN, GENERATE A C1 INTEGRAL FILE.
C     IF THIS IS A PROPERTIES RUN, NO 2E- INTEGRALS HAVE BEEN CREATED.
C
      ABEL = ABELPT()
      NTSAVE = NT
      SVDSKW = DSKWRK
      SVGPAR = GOPARR
      DSKWRK = .TRUE.
      IF((GOPARR.AND.ITRFAO.EQ.1)
     *   .OR.  .NOT.ABEL
     *   .OR.   RUNTYP.EQ.PROP) THEN
         IF(ITRFAO.EQ.1) GOPARR = .FALSE.
         IF(.NOT.ABEL) NT = 1
         IF(.NOT.DIRSCF) CALL JANDK
      END IF
      GOPARR = SVGPAR
      NT = NTSAVE
C
C     LOOP OVER SEPARATE ALPHA AND (MAYBE) BETA LOCALIZATIONS
C
  100 CONTINUE
      DSKWRK = .TRUE.
      IPASS = 0
  110 CONTINUE
      IPASS = IPASS+1
      IF(IPASS.GT.2) GO TO 800
      IF((SCFTYP.EQ.UHF) .AND. MASWRK) WRITE(IW,9010) SPIN(IPASS)
C
      CALL LMOPSI(IPASS,MCORE,MDOC,MACT,NUMLOC)
      IF(IPASS.EQ.1) THEN
         NOUT=NOUTA
      ELSE
         NOUT=NOUTB
      END IF
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      M1 = NUMLOC - NOUT
      M2 = (M1*M1+M1)/2
      M3 = M1*M1
      M4 = (M2*M2+M2)/2
      NUMLC2 = (NUMLOC*NUMLOC+NUMLOC)/2
C
C     ----- SELECT THE ORBITALS WHICH ARE TO BE LOCALIZED -----
C     CANONICAL ORBITALS AT -LCLMO- ARE REORDERED INTO -LCMO-
C
      CALL VALFM(LOADFM)
      LCLMO  = LOADFM + 1
      LCMO   = LCLMO  + L3
      LRHOC  = LCMO   + L3
      LEIGC  = LRHOC  + NUMLC2
      LTRIC  = LEIGC  + L1
      LMAP   = LTRIC  + M2
      LAST   = LMAP   + M1
      NEED1 = LAST - LOADFM -1
      CALL GETFM(NEED1)
C
      CALL LOCSET(IPASS,X(LCLMO),X(LCMO),X(LRHOC),X(LEIGC),X(LTRIC),
     *            X(LMAP),L1,M1,M2,NUMLOC,NUMLC2)
C
C     ----- INTEGRAL TRANSFORMATION -----
C     THIS IS PERFORMED ONLY FOR THE ORBITALS TO BE LOCALIZED.
C
      IF(DIRSCF) DIRTRF=.TRUE.
      IPRT=0
      ICORE = -NOUT
      DOEXCH = .FALSE.
      IF(RUNTYP.EQ.PROP  .AND.  DIRTRF) DOEXCH=.TRUE.
      IF(.NOT. RDINT) THEN
         CALL TRFMCX(IPRT,ICORE,NUMLOC,NUMLOC,.FALSE.,DOEXCH,
     *               .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,
     *               .FALSE.,.FALSE.,.TRUE.)
      END IF
C
C     ----- ALLOCATE MEMORY FOR THE LOCALIZATION -----
C
      IF (MASWRK) THEN
         WRITE(IW,9000)
         IF(SCFTYP.EQ.UHF) WRITE(IW,9010) SPIN(IPASS)
      END IF
      CALL VALFM(LOADFM)
      LIA    = LOADFM + 1
      LXX    = LIA    + M2
      LIX    = LXX    + NINTMX
      LTWOEI = LIX    + NINTMX
      LCINT  = LTWOEI + M4
      LXINT  = LCINT  + M2
      LTLOC  = LXINT  + M2
      LNOROT = LTLOC  + L3
      LAST   = LNOROT + M2
      IF(ISYMLO.EQ.1) THEN
         LS     = LAST
         LQ     = LS     + L2
         LWRK   = LQ     + L3
         LLABMO = LWRK   + L1
         LAST   = LLABMO + M1
      ELSE
         LS     = LAST
         LQ     = LAST
         LWRK   = LAST
         LLABMO = LAST
      END IF
C
      NEED2 = LAST - LOADFM - 1
      NEED = NEED1+NEED2
      IF (MASWRK) WRITE(IW,9020) NEED
      CALL GETFM(NEED2)
      IF(EXETYP.EQ.CHECK) GO TO 600
C
C     ----- READ ALL TRANSFORMED INTEGRALS INTO MEMORY ----
C
      IF (RDINT) THEN
         CALL DAWRIT(IDAF,IODA,X(LTLOC),M3,73,0)
         CALL TWEIIN(X(LCINT),X(LXINT),M2)
         CALL DAREAD(IDAF,IODA,X(LCLMO),L3,15,0)
         GO TO 400
      ELSE
         CALL LOCINT(X(LTWOEI),IW,IJKT,NINTMX,X(LXX),X(LIX),
     *               X(LIA),M2,M4)
      END IF
C
C     ----- CARRY OUT THE LOCALIZATION -----
C
      IF (MASWRK) WRITE(IW,9030) NUMLOC,NOUT
      IF (NUMLOC-NOUT .LE.0) THEN
         IF(MASWRK) WRITE(IW,9040)
         CALL ABRT
      END IF
C
      CALL DAREAD(IDAF,IODA,X(LCMO),L3,15,0)
      CALL LOCFRZ(IPASS,MCORE,MDOC,X(LNOROT),X(LS),X(LQ),
     *            X(LWRK),X(LCMO+L1*NOUT),X(LLABMO),L0,L1,L2,L3,M1,M2)
      CALL LOCLIZ(X(LCMO+L1*NOUT),X(LTLOC),X(LCLMO),
     *            X(LCINT),X(LXINT),X(LTWOEI),X(LIA),X(LNOROT),
     *            IW,L1,M1,M2,M4,CVGLOC,IPRTLO,IPASS)
C
C     ----- ANALYSIS OF FOCK OR DENSITY MATRIX -----
C     STORAGE AT -LCINT- AND -LXINT- USED AS SCRATCH HERE
C
      IF(RUNTYP.NE.PROP) CALL LOCANA(X(LTLOC),X(LTRIC),X(LCINT),
     *                               X(LXINT),M1,M2)
C
C     ----- OUTPUT THE FINAL ORBITALS -----
C
  400 CONTINUE
      CALL LOCFIN(IPASS,X(LCMO),X(LCLMO),L1,M1,NUMLOC)
  600 CONTINUE
      CALL RETFM(NEED2)
      CALL RETFM(NEED1)
      IF(SCFTYP.EQ.UHF) GO TO 110
C
  800 CONTINUE
      DIRTRF = SVDIRT
      DSKWRK = SVDSKW
      IF (MASWRK) WRITE(IW,9050)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/10X,39(1H-)/
     *       10X,'EDMISTON-RUEDENBERG ENERGY LOCALIZATION'/10X,39(1H-)/
     *       10X,'C.EDMISTON AND K.RUEDENBERG,',
     *           ' REV.MOD.PHYS. 35, 457-465 (1963)')
 9010 FORMAT(/1X,'*****',A6,' ORBITAL LOCALIZATION *****')
 9020 FORMAT(/1X,'THIS LOCALIZATION REQUIRES',I10,' WORDS OF MEMORY.')
 9030 FORMAT(1X,'THIS LOCALIZATION HAS',I4,' ORBITALS, OF WHICH',
     *        I4,' ARE FROZEN.')
 9040 FORMAT(/1X,'*** ERROR ***'/
     *       1X,'THIS LOCALIZATION HAS NO ORBITALS TO LOCALIZE.'/
     *       1X,'PLEASE RECONSIDER YOUR FCORE/MOIN/MOOUT INPUT.')
 9050 FORMAT(1X,'... DONE WITH ENERGY LOCALIZATION ...')
      END
C*MODULE LOCAL   *DECK LOCFIN
      SUBROUTINE LOCFIN(IPASS,CMO,CLMO,L1,M1,NUMLOC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXAO=8192)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION CMO(L1,L1),CLMO(L1,M1)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA UHF/8HUHF     /
C
C     ----- FINAL ENERGY LOCALIZATION OUTPUT -----
C     ON EXIT, IPASS=1, 15 = CANONICAL ALPHA MO-S
C                       71 = LOCALIZED ALPHA MO-S
C                       19 = CANONICAL BETA MO-S
C                       72 = NOTHING
C     ON EXIT, IPASS=2, 15 = CANONICAL ALPHA MO-S
C                       71 = LOCALIZED ALPHA MO-S
C                       19 = CANONICAL BETA MO-S
C                       72 = LOCALIZED BETA MO-S
C
      IF(IPASS.EQ.1) THEN
         NOUT=NOUTA
         NDAFC=15
         NDAFL=71
      ELSE
         NOUT=NOUTB
         CALL DAREAD(IDAF,IODA,CMO,L1*L1,19,0)
         CALL DAWRIT(IDAF,IODA,CMO,L1*L1,15,0)
         NDAFC=19
         NDAFL=72
      END IF
C
C     RECOVER THE ORIGINAL CANONICAL ORBITALS.
C
      CALL DAREAD(IDAF,IODA,CMO,L1*L1,NDAFL,0)
      CALL DAWRIT(IDAF,IODA,CMO,L1*L1,NDAFC,0)
C
C     COPY THE MOS WHICH HAVE BEEN LOCALIZED TO THEIR ORIGINAL SPOT.
C
      NMO = 0
      DO 120 J=1,NUMLOC
         DO 110 I=1,NOUT
            IF(IPASS.EQ.1) THEN
               IF(MOOUTA(I).EQ.J) GO TO 120
            ELSE
               IF(MOOUTB(I).EQ.J) GO TO 120
            END IF
  110    CONTINUE
         NMO=NMO+1
         CALL DCOPY(L1,CLMO(1,NMO),1,CMO(1,J),1)
  120 CONTINUE
C
C     OUTPUT THE LOCALIZED ORBITALS
C
      CALL DAWRIT(IDAF,IODA,CMO,L1*L1,NDAFL,0)
C
      IF (MASWRK) WRITE(IW,9010)
      CALL PRSQL(CMO,NUMLOC,L1,L1)
C
      IF((IPASS.EQ.1) .AND. MASWRK) WRITE(IP,9020)
      CALL PUSQL(CMO,NUMLOC,L1,L1)
      IF((IPASS.EQ.2  .OR.  SCFTYP.NE.UHF) .AND. MASWRK)
     *   WRITE(IP,9030)
      RETURN
C
 9010 FORMAT(/1X,'EDMISTON-RUEDENBERG ENERGY LOCALIZED ORBITALS')
 9020 FORMAT('ENERGY LOCALIZED ORBITALS'/' $VEC')
 9030 FORMAT(' $END')
      END
C*MODULE LOCAL   *DECK LOCFRZ
      SUBROUTINE LOCFRZ(IPASS,MCORE,MDOC,NOROT,S,Q,
     *                  WRK,CMO,LABMO,L0,L1,L2,L3,M1,M2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,ABEL,CANONC,FCORE,FORS,NOCI,
     *        EKT,LINSER,FDIRCT,QCORR,ORIENT,ORMFUL
C
      DIMENSION NOROT(M2),S(L2),Q(L3),WRK(L1),CMO(L1,M1),LABMO(L1)
C
      PARAMETER (MXNORO=250)
C
      COMMON /FCCWFN/ NSPACE,MSTA(51),MNUM(51),MINI(51),MAXI(51),
     *                IAMI(51),IAMA(51),IBMI(51),IBMA(51),IDIMFCC(51),
     *                LBST(51),NREF0,FDIRCT,QCORR,C0SQ
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NORTMC(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /ORNTMO/ ORIENT,ORMFUL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMBLK/ NIRRED,NSALC,NSALC2,NSALC3,NSAFMO
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMQMT/ IRPLAB(14),IRPNUM(14),IRPDIM(14),IRPDEG(14)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RMC,ORMAS/8HMCSCF   ,8HORMAS   /
C
C     ----- SET UP DEFAULTS FOR FREEZING ORBITALS -----
C
      DO 100 I=1,M2
         NOROT(I) = 0
  100 CONTINUE
C
C     SET UP AN INVARIANT LOCALIZATION, BY RESTRICTING ANY
C     ROTATIONS MIXING DOUBLY AND PARTIALLY OCCUPIED ORBITALS.
C     LEAVE OUT ANY ROW AND COLUMN INVOLVING A FROZEN ORBITAL.
C
      NUMLOC = M1
      NDBOCC = MCORE + MDOC
      IJ=0
      DO 280 I=1,NUMLOC+NOUTA
         IF(IPASS.EQ.1) THEN
            DO KK=1,NOUTA
               IF(I.EQ.MOOUTA(KK)) GO TO 280
            ENDDO
         ELSE
            DO KK=1,NOUTB
               IF(I.EQ.MOOUTB(KK)) GO TO 280
            ENDDO
         END IF
         DO 270 J=1,I
            IF(IPASS.EQ.1) THEN
               DO 240 KK=1,NOUTA
                  IF(J.EQ.MOOUTA(KK)) GO TO 270
  240          CONTINUE
            ELSE
               DO 250 KK=1,NOUTB
                  IF(J.EQ.MOOUTB(KK)) GO TO 270
  250          CONTINUE
            END IF
            IJ = IJ+1
            IF(I.GT.NDBOCC  .AND.  J.LE.NDBOCC) NOROT(IJ)=1
  270    CONTINUE
  280 CONTINUE
C
C         ORMAS-MCSCF IS NOT INVARIANT TO ROTATIONS THAT MIX
C         ORBITALS BETWEEN THE MULTIPLE ACTIVE SPACES, SO UNLESS
C         THE USER SPECIFICALLY WANTS TO, FREEZE THESE ROTATIONS.
C         STORAGE -LABMO- USED TO DEFINE ORBITAL SPACE COUNTER
C
      IF(.NOT.(SCFTYP.EQ.RMC  .AND.  CISTEP.EQ.ORMAS)) GO TO 500
      IF(ORMFUL) GO TO 500
C
      DO I=1,MDOC
         LABMO(I) = 1
      ENDDO
      DO 330 ISPACE=1,NSPACE
         DO 320 I=MSTA(ISPACE),M1+NOUTA
            IF(IPASS.EQ.1) THEN
               DO KK=1,NOUTA
                  IF(I.EQ.MOOUTA(KK)) GO TO 320
               ENDDO
            ELSE
               DO KK=1,NOUTB
                  IF(I.EQ.MOOUTB(KK)) GO TO 320
               ENDDO
            END IF
            LABMO(I-MCORE) = ISPACE+1
  320    CONTINUE
  330 CONTINUE
C
      IJ=0
      DO I=1,M1
         ISPACE = LABMO(I)
         DO J=1,I
            IJ=IJ+1
            IF(ISPACE.NE.LABMO(J)) NOROT(IJ)=1
         ENDDO
      ENDDO
C
C        OPTIONALLY, RESTRICT ROTATIONS SO THAT ORBITALS WITH
C        DIFFERENT SYMMETRIES ARE NOT MIXED.
C
  500 CONTINUE
      IF(ISYMLO.EQ.0) RETURN
C                     ******
C
C        -SYMMOS- WILL CHANGE ORBITAL PHASES, MAKING THE MOS
C        INCONSISTENT WITH THOSE USED ALREADY IN THE INTEGRAL
C        TRANSFORMATION, AND WITH DATA LIKE DENSITY MATRICES
C        READ FROM DISK FOR EXTENDED ANALYSIS.  SO, WE PRESERVE
C        THE ORIGINAL ORBITAL PHASES BY DUMPING THEM TO DISK,
C        AND THEN RELOADING THEM.  ALL WE WANT IS -LABMO-
C
                     NDAFC=15
      IF(IPASS.EQ.2) NDAFC=19
C
      IF(ILABMO.EQ.0) THEN
         IF (MASWRK) WRITE(IW,9010)
         RETURN
      ELSE
         CALL DAREAD(IDAF,IODA,S,L2,12,0)
         CALL DAREAD(IDAF,IODA,Q,L3,45,0)
         CALL DAWRIT(IDAF,IODA,CMO,L3,NDAFC,0)
         CALL SYMMOS(LABMO,Q,S,CMO,WRK,L0,L1,NUMLOC,L1)
         CALL DAREAD(IDAF,IODA,CMO,L3,NDAFC,0)
      END IF
C
      IJ=0
      IDEG=0
      DO 590 I=1,NUMLOC
         ISYM=LABMO(I)
         IDEG=IDEG+1
         DO 510 IRP=1,NIRRED
            IF(ISYM.EQ.IRPLAB(IRP)) IRREP = IRP
  510    CONTINUE
         IDIM = IRPDIM(IRREP)
         IF(IDIM.EQ.1) IDEG=0
C
         JDEG = 0
         DO 580 J=1,I
            IJ=IJ+1
            JSYM = LABMO(J)
            JDEG = JDEG + 1
            DO 530 IRP=1,NIRRED
               IF(JSYM.EQ.IRPLAB(IRP)) IRREP = IRP
  530       CONTINUE
            JDIM = IRPDIM(IRREP)
            IF(JDIM.EQ.1) JDEG=0
            IF(ISYM.NE.JSYM  .OR.  IDEG.NE.JDEG) NOROT(IJ)=1
            IF(JDEG.EQ.JDIM) JDEG=0
  580    CONTINUE
         IF(IDEG.EQ.IDIM) IDEG=0
  590 CONTINUE
      RETURN
C
 9010 FORMAT(//1X,55(1H*)/
     *       1X,'YOUR GROUP DOES NOT SUPPORT ORBITAL SYMMETRY LABELS.'/
     *       1X,'-SYMLOC- OPTION IS DEACTIVATED.  MAXIMUM LOCALIZATION'/
     *       1X,'CONSISTENT WITH WAVEFUNCTION INVARIANCE WILL BE DONE.'/
     *       1X,55(1H*)//)
      END
C*MODULE LOCAL   *DECK LOCINT
      SUBROUTINE LOCINT(TWOEI,IW,NFTI,NINTMX,XX,IX,IA,M2,M4)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION IA(M2),TWOEI(M4),XX(NINTMX),IX(*)
C
      COMMON /PCKLAB/ LABSIZ
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- READ ALL MO INTEGRALS -----
C     THE ARRAY -TWOEI- SHOULD BE IN CONVENTIONAL TRIANGULAR
C     CANONICAL ORDER, I.GE.J, K.GE.L, AND IJ.GE.KL.
C
      CALL VCLR(TWOEI,1,M4)
C
      II = 0
      DO 100 I=1,M2
         IA(I) = II
         II = II + I
  100 CONTINUE
C
C     TRANSFORMED INTEGRALS ASSUMED IN REVERSE CANONICAL ORDER
C
      CALL SEQREW(NFTI)
      IF (MASWRK) READ(NFTI)
  200 CONTINUE
      CALL PREAD(NFTI,XX,IX,NX,NINTMX)
      IF(NX.EQ.0) GO TO 400
      MX=ABS(NX)
      IF(MX.GT.NINTMX) GO TO 9000
C
C     ----- UNPACK LABEL AND PROCESS INTEGRAL -----
C           FIRST UNPACK 16-BIT LABEL FOR LARGER BASIS SETS
C           IF NOT LARGE ENOUGH, USE 8-BIT LABEL
C
      DO 300 M=1,MX
C
         VAL   = XX(M)
C
                       NPACK = M
                       IF (LABSIZ .EQ. 2) THEN
*I32                     LABEL1 = IX( 2*NPACK - 1 )
*I32                     LABEL2 = IX( 2*NPACK     )
*I32                     IPACK = ISHFT( LABEL1, -16 )
*I32                     JPACK = IAND( LABEL1, 65535 )
*I32                     KPACK = ISHFT( LABEL2, -16 )
*I32                     LPACK = IAND( LABEL2, 65535 )
*I64                     LABEL = IX(NPACK)
*I64                     IPACK = ISHFT( LABEL, -48 )
*I64                     JPACK = IAND( ISHFT( LABEL, -32 ), 65535 )
*I64                     KPACK = IAND( ISHFT( LABEL, -16 ), 65535 )
*I64                     LPACK = IAND( LABEL, 65535 )
                       ELSE IF (LABSIZ .EQ. 1) THEN
*I32                     LABEL = IX(NPACK)
*I32                     IPACK = ISHFT( LABEL, -24 )
*I32                     JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I32                     KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I32                     LPACK = IAND( LABEL, 255 )
*I64                     IF ( MOD(NPACK,2) .EQ. 0 ) THEN
*I64                       LABEL = IX( NPACK/2 )
*I64                       IPACK = IAND( ISHFT( LABEL, -24 ), 255 )
*I64                       JPACK = IAND( ISHFT( LABEL, -16 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL,  -8 ), 255 )
*I64                       LPACK = IAND( LABEL, 255 )
*I64                     ELSE
*I64                       LABEL = IX( (NPACK/2)+1 )
*I64                       IPACK = ISHFT( LABEL, -56 )
*I64                       JPACK = IAND( ISHFT( LABEL, -48 ), 255 )
*I64                       KPACK = IAND( ISHFT( LABEL, -40 ), 255 )
*I64                       LPACK = IAND( ISHFT( LABEL, -32 ), 255 )
*I64                     END IF
                       END IF
                       K = IPACK
                       L = JPACK
                       I = KPACK
                       J = LPACK
C
         IJ = IA(I) + J
         KL = IA(K) + L
         IF(I.EQ.K  .AND.  J.LE.L) THEN
            IJ = IA(I) + L
            KL = IA(K) + J
         END IF
         IJKL = IA(IJ) + KL
         TWOEI(IJKL) = VAL
  300 CONTINUE
C
      IF(NX.GT.0) GO TO 200
  400 CONTINUE
      IF(GOPARR) CALL DDI_GSUMF(1300,TWOEI,M4)
      RETURN
C
 9000 CONTINUE
      IF (MASWRK) WRITE(IW,*)
     *   'LOCINT: ERROR READING TRANSFORMED INTEGRAL FILE'
      CALL ABRT
      STOP
      END
C*MODULE LOCAL   *ECK LOCLIZ
      SUBROUTINE LOCLIZ(CMO,TLOC,CLMO,CINT,XINT,TWOEI,IA,NOROT,
     *                  IW,L1,M1,M2,M4,CVGLOC,IPRTLO,IPASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION CMO(L1,M1),TLOC(M1,M1),CLMO(L1,M1),IA(M2),
     *          CINT(M2),XINT(M2),TWOEI(M4),NOROT(M2)
C
      COMMON /IOFILE/ IR,IWDUM,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00, TENTH=0.1D+00,
     *           TEN6=1.0D+06)
C
C        ----- CONVERGENCE DRIVER FOR ENERGY LOCALIZATION -----
C
      II = 0
      DO 100 I=1,M2
         IA(I) = II
         II = II + I
  100 CONTINUE
C
      CALL VCLR(TLOC,1,M1*M1)
      DO 130 J=1,M1
         TLOC(J,J) = ONE
  130 CONTINUE
C
C        PRINT INITIAL COULOMB AND EXCHANGE INTEGRALS.
C        COMPUTE INITIAL LOCALIZATION, COULOUMB, AND EXCHANGE SUMS.
C
      IF(IPRTLO.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9050) 'INITIAL'
         CALL LOCPCX(IW,TWOEI,CINT,XINT,IA,M1,M2,M4)
      END IF
      IF (MASWRK) WRITE(IW,9060) 'INITIAL'
      CALL LOCCHK(IW,TWOEI,IA,M1,M2,M4)
C
C        CARRY OUT LOCALIZATION ROTATIONS AT 6 CONVERGENCE LEVELS
C        ONLY THE MASTER NODE CARRIES OUT THE JACOBI INTEGRAL UPDATES,
C        SO LOCALIZING IS A SEQUENTIAL BOTTLENECK.
C
      IF(MASWRK) THEN
         WRITE(IW,9040)
         NTROT=0
         LOCOPT=0
         CVG = CVGLOC*TEN6
         IF(CVGLOC.LT.0) GO TO 301
         DO 300 LOCYCL=1,6
            CVG=TENTH*CVG
            CALL LOCROT(IW,TWOEI,TLOC,IA,NOROT,CVG,LOCOPT,NROT,M1,M2,M4)
            NTROT=NTROT+NROT
  300    CONTINUE
  301    CONTINUE
         WRITE(IW,9000) NTROT,CVG
      END IF
C
C        PRINT FINAL COULOMB AND EXCHANGE INTEGRALS.
C        COMPUTE FINAL LOCALIZATION, COULOUMB, AND EXCHANGE SUMS.
C
      IF(IPRTLO.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9050) '  FINAL'
         CALL LOCPCX(IW,TWOEI,CINT,XINT,IA,M1,M2,M4)
         IF (MASWRK) THEN
            WRITE(IP,*) '$TWOEI'
            WRITE(IP,FMT='(1X,1P,5E15.8)') (CINT(J),J=1,M2)
            WRITE(IP,FMT='(1X,1P,5E15.8)') (XINT(J),J=1,M2)
            WRITE(IP,*) '$END'
         END IF
         CALL DAWRIT(IDAF,IODA,CINT,M2,59,0)
         CALL DAWRIT(IDAF,IODA,XINT,M2,60,0)
      END IF
      IF (MASWRK) WRITE(IW,9060) '  FINAL'
      CALL LOCCHK(IW,TWOEI,IA,M1,M2,M4)
C
C        TRANSFORM INPUT ORBITALS TO LOCALIZED ORBITALS.
C
      CALL MRARBR(CMO,L1,L1,M1,TLOC,M1,M1,CLMO,L1)
C
C        PRINT LOCALIZATION TRANSFORMATION
C
      IF(IPRTLO.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9020)
         CALL PRSQ(TLOC,M1,M1,M1)
      END IF
C
C     ----- SAVE FOR MP2 LOCALIZATIONS -----
C
      IF(IPASS.EQ.1) CALL DAWRIT(IDAF,IODA,TLOC,M1*M1,73,0)
      IF(IPASS.EQ.2) CALL DAWRIT(IDAF,IODA,TLOC,M1*M1,74,0)
      RETURN
C
 9000 FORMAT(1X,I10,' TOTAL PAIR ROTATIONS REQUIRED TO CONVERGE TO',
     *       1P,E12.2)
 9020 FORMAT(/1X,'LOCALIZATION TRANSFORMATION')
 9040 FORMAT(/1X,'BEGINNING 2X2 JACOBI ROTATIONS...')
 9050 FORMAT(/1X,A7,' COULOMB AND EXCHANGE INTEGRALS')
 9060 FORMAT(/1X,A7,' LOCALIZATION SUMS ARE')
      END
C*MODULE LOCAL   *DECK LOCPCX
      SUBROUTINE LOCPCX(IW,TWOEI,CINT,XINT,IA,M1,M2,M4)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION TWOEI(M4),CINT(M2),XINT(M2),IA(M2)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        ----- PRINT COULOMB AND EXCHANGE INTEGRALS -----
C        ON ENTRY, ALL TWO ELECTRON MO INTEGRALS ARE IN -TWOEI-.
C        COPY COULOMB AND EXCHANGE INTEGRALS INTO -CINT- AND
C        -XINT-, AND PRINT THESE FOR THE USER'S CONVENIENCE.
C
      CALL VCLR(CINT,1,M2)
      CALL VCLR(XINT,1,M2)
C
      II=0
      IJ=0
      DO 120 I=1,M1
         II=II+I
         IIA=IA(II)
         JJ=0
         DO 110 J=1,I
            IJ=IJ+1
            JJ=JJ+J
            IJIJ=IA(IJ)+IJ
            IIJJ=IIA+JJ
            CINT(IJ)=TWOEI(IIJJ)
            XINT(IJ)=TWOEI(IJIJ)
  110    CONTINUE
  120 CONTINUE
C
      IF (MASWRK) WRITE(IW,9000) ' COULOMB'
      CALL PRTRI(CINT,M1)
      IF (MASWRK) WRITE(IW,9000) 'EXCHANGE'
      CALL PRTRI(XINT,M1)
      RETURN
C
 9000 FORMAT(/1X,A8,' INTEGRAL MATRIX')
      END
C*MODULE LOCAL   *DECK LOCROT
      SUBROUTINE LOCROT(IW,TWOEI,TLOC,IA,NOROT,CVG,LOCOPT,NROT,
     *                  M1,M2,M4)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION TLOC(M1,M1),IA(M2),TWOEI(M4),NOROT(M2)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (PIOVR4=7.853981633974D-01, ZERO=0.0D+00,
     *           TENTH=0.1D+00, QUARTR=0.25D+00, HALF=0.5D+00,
     *           ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
C
C     ----- ENERGY LOCALIZATION VIA A JACOBI TECHNIQUE -----
C     TWO BY TWO ORBITAL ROTATIONS, USING INTEGRAL UPDATING.
C     C.EDMISTON AND K.RUEDENBERG, REV.MOD.PHYS. 35, 457 (1963)
C
C  TWOEI    ARRAY OF MINIMUM NUMBER OF MOLECULAR INTEGRALS STORED
C           WITH ORDINARY SYMMETRY; I.GE.J, K.GE.L, (IJ).GE.(KL)
C  TLOC     TRANSFORMATION MATRIX (INITIALIZE TO UNITY BEFORE CALLING)
C  M1       NUMBER OF MOS TO BE LOCALIZED
C  CVG      CONVERGENCE CRITERION ON THE TRANSFORMATION MATRIX.
C           1.0D-06 IMPLIES BETTER THAN SIX FIGURES IN THE TRANS.
C  LOCOPT   NONZERO TO DELOCALIZE ORBITALS, OTHERWISE LOCALIZE
C  NROT     RETURNS THE NUMBER OF PAIR ROTATIONS ACTUALLY CARRIED OUT
C
C           LL > LM > MM
C           LLLL > LLLM > LLMM > LMLM > LMMM > MMMM
C
      NROT=0
      NTRI=IA(M1) + M1
      DIV=TENTH
      THR=ONE
C
  130 CONTINUE
      THR=THR*DIV
  140 CONTINUE
      IND = 0
C
      DO 370 L=2,M1
         LMO=L - 1
         LT=IA(L)
         LL=LT + L
         LLT=IA(LL)
         LLLL=LLT + LL
         DO 360 M=1,LMO
            LM=LT + M
            IF(NOROT(LM).EQ.1) GO TO 360
            LMT=IA(LM)
            CLLLL=TWOEI(LLLL)
            LLLM=LLT + LM
            CLLLM=TWOEI(LLLM)
            MT=IA(M)
            MM=MT + M
            LMMM=LMT + MM
            CLMMM=TWOEI(LMMM)
            BLM=CLLLM - CLMMM
            LMLM=LMT + LM
            CLMLM=TWOEI(LMLM)
            LLMM=LLT + MM
            CLLMM=TWOEI(LLMM)
            MMT=IA(MM)
            MMMM=MMT + MM
            CMMMM=TWOEI(MMMM)
            ALM=QUARTR*(CLLLL + CMMMM - CLLMM - CLLMM) - CLMLM
            ALPHA=QUARTR*ATAN2(BLM,ALM)
            IF (LOCOPT.NE.0 .AND. ALM.GT.ZERO) ALPHA=ALPHA+PIOVR4
            IF(ABS(ALPHA).LT.THR)GO TO 360
               NROT=NROT + 1
C
               IND=1
C
               SINA=SIN(ALPHA)
               COSA=COS(ALPHA)
               SINA2=SINA*SINA
               COSA2=COSA*COSA
               SCA=SINA*COSA
               COS2A=COSA2 - SINA2
               SIN2A=SCA + SCA
               COS2A2=COS2A*COS2A
               SIN2A2=SIN2A*SIN2A
               SC2A=SIN2A*COS2A
               SIN4A=SC2A + SC2A
               COSA4=COSA2*COSA2
               SINA4=SINA2*SINA2
C
               DO 230 IJ=1,NTRI
                  IJT=IA(IJ)
                  IF(IJ.GT.LL)GO TO 170
                     IF(IJ.EQ.LL)GO TO 230
                        LLIJ=LLT + IJ
                        GO TO 180
C
  170                LLIJ=IJT + LL
  180                IF(IJ.GT.LM)GO TO 190
                        IF(IJ.EQ.LM)GO TO 230
                           LMIJ=LMT + IJ
                           GO TO 200
C
  190                   LMIJ=IJT + LM
  200                   IF(IJ.GT.MM)GO TO 210
                           IF(IJ.EQ.MM)GO TO 230
                              MMIJ=MMT + IJ
                              GO TO 220
C
  210                      MMIJ=IJT + MM
  220                      CONTINUE
                           CLLIJ=TWOEI(LLIJ)
                           CLMIJ=TWOEI(LMIJ)
                           CMMIJ=TWOEI(MMIJ)
                           TWOEI(LMIJ)=COS2A*CLMIJ-SCA*(CLLIJ-CMMIJ)
                           T=COSA2*CMMIJ+SINA2*CLLIJ- SIN2A*CLMIJ
                           TWOEI(MMIJ)=T
                           TWOEI(LLIJ)=CLLIJ + CMMIJ - T
  230          CONTINUE
C
          TWOEI(LLLL)=COSA4*CLLLL + SINA4*CMMMM + SIN2A2*(HALF*
     *         CLLMM + CLMLM) + TWO*SIN2A*(COSA2*CLLLM + SINA2*CLMMM)
          TWOEI(MMMM)=COSA4*CMMMM + SINA4*CLLLL + SIN2A2*(HALF*
     *         CLLMM + CLMLM) - TWO*SIN2A*(SINA2*CLLLM + COSA2*CLMMM)
          TWOEI(LMLM)=COS2A2*CLMLM - SC2A*(CLLLM - CLMMM)
     *         + QUARTR*SIN2A2*(CLLLL + CMMMM - CLLMM - CLLMM)
          TWOEI(LLMM)=(COSA4 + SINA4)*CLLMM + QUARTR*SIN2A2*
     *         (CLLLL+ CMMMM-FOUR*CLMLM)-SC2A*(CLLLM-CLMMM)
          TWOEI(LLLM)=(COS2A*COSA2 - HALF*SIN2A2)*CLLLM
     *         + (COS2A*SINA2 + HALF*SIN2A2)*CLMMM
     *         + QUARTR*SIN4A*(CLMLM + CLMLM + CLLMM)
     *         - SCA*(COSA2*CLLLL - SINA2*CMMMM)
          TWOEI(LMMM)=(COS2A*COSA2 - HALF*SIN2A2)*CLMMM
     *         + (COS2A*SINA2 + HALF*SIN2A2)*CLLLM
     *         - QUARTR*SIN4A*(CLLMM + CLMLM + CLMLM)
     *         + SCA*(COSA2*CMMMM - SINA2*CLLLL)
C
          DO 340 K=1,M1
             IF(K.EQ.L.OR.K.EQ.M)GO TO 340
                KT=IA(K)
                IF(K.GT.L)GO TO 240
                   KL=LT + K
                   GO TO 250
C
  240           KL=  KT + L
  250           IF(K.GT.M)GO TO 260
                   KM=MT + K
                   GO TO 270
C
  260           KM=KT + M
  270           CONTINUE
                KLT=IA(KL)
                KMT=IA(KM)
                DO 330 IJ=1,NTRI
                   IJT=IA(IJ)
                   IF(KL.GT.IJ)GO TO 280
                      KLIJ=IJT + KL
                      GO TO 290
C
  280              KLIJ=KLT + IJ
  290              IF(KM.GT.IJ)GO TO 300
                      IF(KM.EQ.IJ)GO TO 330
                         KMIJ=IJT + KM
                         GO TO 310
C
  300                 KMIJ=KMT + IJ
  310                 CONTINUE
                      CKLIJ=TWOEI(KLIJ)
                      CKMIJ=TWOEI(KMIJ)
                      IF(KL.EQ.IJ)GO TO 320
                         TWOEI(KLIJ)=COSA*CKLIJ + SINA*CKMIJ
                         TWOEI(KMIJ)=COSA*CKMIJ - SINA*CKLIJ
                         GO TO 330
C
  320                 KMKM=KMT + KM
                      KLKL=KLIJ
                      KLKM=KMIJ
                      CKMKM=TWOEI(KMKM)
                      CKLKL=CKLIJ
                      CKLKM=CKMIJ
                      T=COSA2*CKLKL + SINA2*CKMKM + SIN2A*CKLKM
                      TWOEI(KLKL)=T
                      TWOEI(KMKM)=CKLKL + CKMKM - T
                      TWOEI(KLKM)=COS2A*CKLKM - SCA*(CKLKL - CKMKM)
  330           CONTINUE
  340     CONTINUE
C
            DO 350 K=1,M1
               T1 =  COSA*TLOC(K,L) + SINA*TLOC(K,M)
               T2 = -SINA*TLOC(K,L) + COSA*TLOC(K,M)
               TLOC(K,L) = T1
               TLOC(K,M) = T2
  350       CONTINUE
  360    CONTINUE
  370 CONTINUE
C
      IF(IND.EQ.1) GO TO 140
      IF(THR.GT.CVG) GO TO 130
C
      IF (MASWRK) WRITE(IW,9000) NROT, CVG
      RETURN
C
 9000 FORMAT(1X,I10,' PAIR ROTATIONS WHEN CONVERGED TO',1P,E12.2)
      END
C*MODULE LOCAL   *DECK LOCSET
      SUBROUTINE LOCSET(IPASS,VCANON,VLOC,RHOC,EIGC,TRIC,MAP,
     *                  L1,M1,M2,NUMLOC,NUMLC2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION VCANON(L1,L1),VLOC(L1,L1),RHOC(NUMLC2),EIGC(L1),
     *          TRIC(M2),MAP(M1)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RHF,UHF,AMCSCF/8HRHF     ,8HUHF     ,8HMCSCF   /
      DATA CHECK,PROP/8HCHECK   ,8HPROP    /
C
C     ----- REORDER THE ORBITALS FOR ENERGY LOCALIZATION -----
C     ALL FROZEN ORBITALS MUST BE IN FRONT OF ALL ORBITALS
C     WHICH ARE TO BE LOCALIZED, FOR THE INTEGRAL TRANSFORMATION.
C     THE TRANSFORMATION WORKS WITH WHATEVER ORBITALS ARE STORED
C     ON DAF RECORD 15.
C
      IF(IPASS.EQ.1) THEN
         NOUT=NOUTA
         NDAFC=15
         NDAFE=17
         NDAFL=71
      ELSE
         NOUT=NOUTB
         NDAFC=19
         NDAFE=21
         NDAFL=72
      END IF
C
C     LOAD ORBITALS, AND MAYBE EIGENVALUES OR DENSITY MATRIX
C
      CALL DAREAD(IDAF,IODA,VCANON,L1*L1,NDAFC,0)
      IF(RUNTYP.EQ.PROP  .OR.  EXETYP.EQ.CHECK) THEN
         IF(SCFTYP.EQ.RHF)    CALL VCLR(EIGC,1,L1)
         IF(SCFTYP.EQ.UHF)    CALL VCLR(EIGC,1,L1)
         IF(SCFTYP.EQ.AMCSCF) CALL VCLR(RHOC,1,NUMLC2)
      ELSE
         IF(SCFTYP.EQ.RHF)    CALL DAREAD(IDAF,IODA,EIGC,L1,NDAFE,0)
         IF(SCFTYP.EQ.UHF)    CALL DAREAD(IDAF,IODA,EIGC,L1,NDAFE,0)
         IF(SCFTYP.EQ.AMCSCF) CALL DAREAD(IDAF,IODA,RHOC,NUMLC2,68,0)
      END IF
C
C     A COPY OF THE ORIGINAL ORDER ORBITALS IS KEPT.
C
      CALL DAWRIT(IDAF,IODA,VCANON,L1*L1,NDAFL,0)
      CALL VCLR(VLOC,1,L1*L1)
C
C     PUT THE ORBITALS WHICH ARE TO BE FROZEN FIRST...
C
      NMO=0
      DO 120 J=1,NUMLOC
         DO 110 I=1,NOUT
            IF((IPASS.EQ.1  .AND.  MOOUTA(I).EQ.J)  .OR.
     *         (IPASS.EQ.2  .AND.  MOOUTB(I).EQ.J)) THEN
                  NMO = NMO+1
                  CALL DCOPY(L1,VCANON(1,J),1,VLOC(1,NMO),1)
                  GO TO 120
            END IF
  110    CONTINUE
  120 CONTINUE
C
C     ...THEN PUT IN THE ORBITALS TO BE LOCALIZED...
C
      LL=0
      DO 220 J=1,NUMLOC
         DO 210 I=1,NOUT
            IF(IPASS.EQ.1) THEN
               IF(MOOUTA(I).EQ.J) GO TO 220
            ELSE
               IF(MOOUTB(I).EQ.J) GO TO 220
            END IF
  210    CONTINUE
         NMO = NMO+1
         CALL DCOPY(L1,VCANON(1,J),1,VLOC(1,NMO),1)
         LL = LL+1
         MAP(LL) = J
  220 CONTINUE
C
      IF(NMO.NE.NUMLOC) THEN
         IF (MASWRK) WRITE(IW,9000) NMO,NUMLOC
         CALL ABRT
      END IF
C
C     ...AND FINISH UP WITH THE VIRTUAL SPACE.
C
      LVIRT = (L1-NUMLOC)*L1
      CALL DCOPY(LVIRT,VCANON(1,NMO+1),1,VLOC(1,NMO+1),1)
C
C        RHF/UHF WILL LATER TRANSFORM CANONICAL FOCK MATRIX TO LMO BASIS
C        STORE THE EIGENVALUES DOWN THE DIAGONAL OF A SYMMETRIC MATRIX
C
      IF(SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF) THEN
         CALL VCLR(TRIC,1,M2)
         II = 0
         DO 310 I=1,M1
            II = II+I
            TRIC(II) = EIGC(MAP(I))
  310    CONTINUE
      END IF
C
C        MCSCF WILL LATER TRANSFORM DENSITY MATRIX TO LMO BASIS
C        PICK OUT THE DENSITY ELEMENTS FOR THE ORBITALS BEING LOCALIZED
C
      IF(SCFTYP.EQ.AMCSCF) THEN
         IJ=0
         DO 330 I=1,M1
            IMAP = MAP(I)
            DO 320 J=1,I
               IJ=IJ+1
               JMAP = MAP(J)
               IJMAP = (IMAP*IMAP-IMAP)/2 + JMAP
               TRIC(IJ) = RHOC(IJMAP)
  320       CONTINUE
  330    CONTINUE
      END IF
C
C     THE INTEGRAL TRANSFORMATION WORKS WITH WHATEVER ORBITALS
C     HAPPEN TO BE STORED IN RECORD 15.  IF THIS IS A UHF RUN,
C     WE MUST GRAB THE ALPHA ORBITALS, AND SAVE THEM FOR
C     LATER RESTORATION.  ON EXIT, WE WILL HAVE
C        IPASS=1      15 = ALPHA'S READY TO BE LOCALIZED
C                     71 = CANONICAL ALPHA MO-S
C                     19 = CANONICAL BETA MO-S
C                     72 = NOTHING
C        IPASS=2      15 = BETA'S READY TO BE LOCALIZED
C                     71 = LOCALIZED ALPHA MO-S
C                     19 = CANONICAL ALPHA MO-S
C                     72 = CANONICAL BETA MO-S
C
      IF(IPASS.EQ.2) THEN
         CALL DAREAD(IDAF,IODA,VCANON,L1*L1,15,0)
         CALL DAWRIT(IDAF,IODA,VCANON,L1*L1,19,0)
      END IF
      CALL DAWRIT(IDAF,IODA,VLOC,L1*L1,15,0)
      RETURN
C
 9000 FORMAT(1X,'PROBLEM GETTING ORBITALS TO BE LOCALIZED'/
     *       1X,'NMO=',I5,' BUT NUMLOC=',I5)
      END
C*MODULE LOCAL   *DECK POPLOC
      SUBROUTINE POPLOC(IPASS,VECIN,TRAN,VECOUT,IORD,IIR,
     *                  MAP,RIJ,QPIX,QPJX,NOROT,SAO,
     *                  L1,L2,N1,N2,NATOMS,NOUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LMORND
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION VECIN(L1,N1),TRAN(N1,N1),SAO(L2),NOROT(N1,N1),
     *          VECOUT(L1,N1),RIJ(N2,NATOMS),
     *          IIR(L1),MAP(L1),IORD(L1),QPIX(L1),QPJX(L1)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /ATLIM / LIMLOW(MXATM),LIMSUP(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00,
     *           TENM3=1.0D-03, TENM8=1.0D-08, TENM10=1.0D-10)
C
      DATA UHF/8HUHF     /
C
      INDX(I,J) = ((MAX(I,J)*(MAX(I,J)-1))/2+MIN(I,J))
C
C     POPLOC CALCULATES THE LOCALIZED MOLECULAR ORBITALS BY THE
C     POPULATION METHOD OF J. PIPEK AND P. G. MEZEY,
C     J. CHEM. PHYS. 90, 4916 (1989),
C     THIS IMPLEMENTATION IS DERIVED FROM QCPE PROGRAM 354,
C     AND WAS WRITTEN BY JOHN MONTGOMERY, JULY 1991
C
      NORB  = N1
      NBASIS= L1
      NFROZ = NOUT
      MAXIT = MAXLOC
C
      CALL AOLIM
      DO 20 K=1,NBASIS
      DO 20 I=1,NATOMS
      IF(LIMLOW(I).LE.K.AND.LIMSUP(I).GE.K) MAP(K)=I
   20 CONTINUE
C
      NREDO = 0
  110 CONTINUE
      NREDO=NREDO+1
C
C          CONSTRUCT INITIAL ATOMIC POPULATIONS
C
      CALL VCLR(RIJ,1,N2*NATOMS)
      CALL DAREAD(IDAF,IODA,SAO,L2,12,0)
      IJ = 0
      DO 280 I = 1,NORB
         DO 280 J = 1,I
            IJ = IJ+1
            DO 260 K = 1,NBASIS
               KK = MAP(K)
               DO 260 L = 1,NBASIS
                  LL = MAP(L)
                  KL = INDX(K,L)
                  SUM = VECIN(K,I)*VECIN(L,J)*SAO(KL)/TWO
                  RIJ(IJ,KK) = RIJ(IJ,KK) + SUM
                  RIJ(IJ,LL) = RIJ(IJ,LL) + SUM
  260       CONTINUE
  280 CONTINUE
C
C          COMPUTE INITIAL LOCALIZATION SUM
C
      SUMRR = ZERO
      DO 320 I=1,NORB
         IF(NOROT(I,I).EQ.1) GO TO 320
         II = INDX(I,I)
         DO 310 K=1,NATOMS
            SUMRR = SUMRR + RIJ(II,K)**2
  310    CONTINUE
  320 CONTINUE
      IF (MASWRK) WRITE(IW,9010) SUMRR
      IF(IPRTLO.EQ.1) CALL POPQAT(RIJ,VECOUT,N2,NORB,NATOMS)
C
C          SEED THE RANDOM FUNCTION, INITIALIZE
C          THE LOCALIZATION TRANSFORMATION, ETC.
C
      ONEPT0 = ONE
      IF (NREDO.EQ.1) XX = LMORND(ONEPT0,VECIN,L1)
C
      CALL VCLR(TRAN,1,NORB*NORB)
      DO 340 I = 1,NORB
         TRAN(I,I) = ONE
  340 CONTINUE
      ITER = 0
      SHIFT = ATAN(ONEPT0)
C
C          BEGIN LOCALIZATION CYCLES
C
  360 CONTINUE
      CHANGE = ZERO
      ITER = ITER+1
      DO 380 I = 1,NORB
         IIR(I) = I
  380 CONTINUE
      NNN = NORB
      DO 400 I = 1,NORB
         XX = LMORND(CHANGE,VECIN,L1)
         III = INT(XX*NNN+ONE)
         IORD(I) = IIR(III)
         IIR(III) = IIR(NNN)
         NNN = NNN-1
  400 CONTINUE
C
C        FOR EACH PAIR OF ORBITALS A TWO DIMENSIONAL UNITARY
C        TRANSFORMATION IS PERFORMED. THE TRANSFORMATION IS
C
C           PSI'(I) =  COS(T)*PSI(I) + SIN(T)*PSI(J)  AND
C           PSI'(J) = -SIN(T)*PSI(I) + COS(T)*PSI(J).
C
C        LOCALIZATION REQUIRES THAT T BE SUCH AS TO MAXIMIZE
C        THE SUM OF THE SQUARES OF THE ATOMIC POPULATIONS.
C
      DO 920 III = 1,NORB
      I = IORD(III)
      IF(NOROT(I,I).EQ.1) GO TO 920
      II = INDX(I,I)
      JM = 1
      RM = ZERO
      TM = ZERO
      SM = ZERO
      CM = ONE
      DO 580 J = 1,NORB
      IF(I.EQ.J) GO TO 580
      IF(NOROT(I,J).EQ.1) GO TO 580
      IJ = INDX(I,J)
      JJ = INDX(J,J)
      T = ZERO
      TX = ZERO
      DO 480 KK = 1,NATOMS
         T = T + FOUR*RIJ(IJ,KK)**2 - RIJ(II,KK)**2 - RIJ(JJ,KK)**2
     *         + TWO*RIJ(II,KK)*RIJ(JJ,KK)
         TX = TX + RIJ(IJ,KK)*(RIJ(JJ,KK) - RIJ(II,KK))
  480 CONTINUE
      IF ((ABS(T) .LE. TENM10) .AND. (ABS(TX) .LE. TENM10)) GO TO 580
      TX = FOUR*TX
      T = ATAN2(TX,T)/FOUR
      SIGN = ONE
      IF (T .GT. ZERO) SIGN = -ONE
      T = T+SIGN*SHIFT
      ITIM = 0
  500 ITIM = ITIM+1
      S = SIN(T)
      C = COS(T)
      RIN = ZERO
      DO 520 KK = 1,NATOMS
         QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
         QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
         RIN = RIN+QPI*QPI+QPJ*QPJ-RIJ(II,KK)**2-RIJ(JJ,KK)**2
  520 CONTINUE
      TTEST = ABS(T)-SHIFT
      IF ((ABS(T) .LE. TENM8) .OR. (ABS(TTEST) .LE. TENM8)) GO TO 560
      IF (RIN .GE. -TENM8) GO TO 560
      IF (ITIM .LE. 1) GO TO 540
      IF (MASWRK) THEN
         WRITE (IW,9020) I,J
         WRITE (IW,9030) T,S,C,RIN
      END IF
      RETURN
C
  540 SIGN = ONE
      IF (T .GT. ZERO) SIGN = -ONE
      T = T+SHIFT*SIGN
      GO TO 500
C
  560 IF (RIN .LE. RM) GO TO 580
      RM = RIN
      TM = T
      SM = S
      CM = C
      JM = J
  580 CONTINUE
C
      RIN = RM
      T = TM
      S = SM
      C = CM
      J = JM
      IJ = INDX(I,J)
      JJ = INDX(J,J)
      IF(NOROT(I,J).EQ.1) GO TO 920
C
C        ACCUMULATE THE 2X2 ROTATION
C
      CHANGE = CHANGE+T*T
      CALL DROT(NORB,TRAN(1,I),1,TRAN(1,J),1,C,S)
C
C        UPDATE THE ATOMIC POPULATIONS
C
      DO 880 KK = 1,NATOMS
         QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
         QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
         QPIJ = (C*C-S*S)*RIJ(IJ,KK)+C*S*(RIJ(JJ,KK)-RIJ(II,KK))
         DO 720 K = 1,NORB
            IF (I.EQ.K.OR.J.EQ.K) GO TO 720
            IK = INDX(I,K)
            JK = INDX(J,K)
            QPIX(K) = C*RIJ(IK,KK)+S*RIJ(JK,KK)
            QPJX(K) = C*RIJ(JK,KK)-S*RIJ(IK,KK)
            RIJ(IK,KK) = QPIX(K)
            RIJ(JK,KK) = QPJX(K)
  720    CONTINUE
         RIN = RIN+QPI+QPJ-RIJ(II,KK)-RIJ(JJ,KK)
         RIJ(II,KK) = QPI
         RIJ(JJ,KK) = QPJ
         RIJ(IJ,KK) = QPIJ
  880 CONTINUE
  920 CONTINUE
C
C          TEST FOR CONVERGENCE OF LOCALIZATION PROCEDURE
C
      CHANGE = SQRT(TWO*CHANGE/(NORB*(NORB-1)))
      IF((NPRINT.EQ.-7) .AND. MASWRK) WRITE(IW,9050) ITER,CHANGE
      IF((MOD(ITER,10).EQ.0) .AND. MASWRK) WRITE(IW,9050) ITER,CHANGE
      IF(ITER.LT.MAXIT  .AND.  CHANGE.GT.TENM3*CVGLOC) GO TO 360
      IF(CHANGE.LE.CVGLOC) GO TO 1000
         IF((NREDO.LE.2) .AND. MASWRK) WRITE(IW,9060)
         IF(NREDO.LE.2) GO TO 110
            IF (MASWRK) WRITE(IW,9070)
            RETURN
C
C          FINISHED WITH LOCALIZATION CYCLES
C
 1000 CONTINUE
      IF (MASWRK) WRITE (IW,9080) ITER
C
C        TRANSFORM TO FINAL ORBITALS, COPY VIRTUAL SPACE
C
      CALL MRARBR(VECIN,NBASIS,NBASIS,NORB,TRAN,NORB,NORB,VECOUT,NBASIS)
      NVIRT = NBASIS - NORB
      CALL DCOPY(NVIRT*NBASIS,VECIN(1,NORB+1),1,VECOUT(1,NORB+1),1)
C
C          PRINT/PUNCH FINAL LOCALIZED ORBITALS
C
      IF (MASWRK) WRITE(IW,9090)
      CALL PRSQL(VECOUT,NORB,L1,L1)
C
      IF(IPRTLO.EQ.1) THEN
         IF (MASWRK) WRITE(IW,9120)
         CALL PRSQ(TRAN,NORB,NORB,N1)
      END IF
C
C     ----- SAVE THIS FOR EVENTUAL MP2 LOCALIZATION -----
C     WE NEED THE MATRIX WITHOUT THE CORE
C
      DO 1110 I = 1,NORB
         DO 1100 J = 1,NORB
            TRAN(I,J) = TRAN(I+NOUT,J+NOUT)
 1100    CONTINUE
 1110 CONTINUE
      NLOC = NORB - NOUT
      CALL SQZMAT(TRAN,NORB,NLOC,NORB)
      IF(IPASS.EQ.1) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,73,0)
      IF(IPASS.EQ.2) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,74,0)
C
      IF (MASWRK) WRITE (IW,9130) NORB,NFROZ,NBASIS
C
      IF((IPASS.EQ.1) .AND. MASWRK) WRITE(IP,8000)
      CALL PUSQL(VECOUT,NORB,L1,L1)
      IF((SCFTYP.NE.UHF  .OR.  IPASS.EQ.2) .AND. MASWRK)
     *   WRITE(IP,8010)
C
C          CONSTRUCT FINAL ATOMIC POPULATIONS
C
      CALL VCLR(RIJ,1,N2*NATOMS)
      CALL DAREAD(IDAF,IODA,SAO,L2,12,0)
      IJ = 0
      DO 1280 I = 1,NORB
         DO 1280 J = 1,I
            IJ = IJ+1
            DO 1260 K = 1,NBASIS
               KK = MAP(K)
               DO 1260 L = 1,NBASIS
                  LL = MAP(L)
                  KL = INDX(K,L)
                  SUM = VECOUT(K,I)*VECOUT(L,J)*SAO(KL)/TWO
                  RIJ(IJ,KK) = RIJ(IJ,KK) + SUM
                  RIJ(IJ,LL) = RIJ(IJ,LL) + SUM
 1260       CONTINUE
 1280 CONTINUE
C
C               COMPUTE FINAL LOCALIZATION SUM (SKIP FROZEN MO-S)
C
      SUMRR = ZERO
      DO 1380 I=1,NORB
         IF(NOROT(I,I).EQ.1) GO TO 1380
         II = INDX(I,I)
         DO 1370 K = 1,NATOMS
            SUMRR = SUMRR + RIJ(II,K)**2
 1370    CONTINUE
 1380 CONTINUE
      IF (MASWRK) WRITE(IW,9140) SUMRR
      IF(IPRTLO.EQ.1) CALL POPQAT(RIJ,VECIN,N2,NORB,NATOMS)
      RETURN
C
 8000 FORMAT('PIPEK-MEZEY POPULATION LOCALIZED ORBITALS'/' $VEC')
 8010 FORMAT(' $END')
C
 9010 FORMAT(/10X,'THE INITIAL LOCALIZATION SUM IS',F14.6)
 9020 FORMAT(1H0,'NO ROTATION INCREASES ATOMIC POPULATIONS',
     *           ' --- LOCALIZATION ABORTED'/
     *           10X,2HI=,I3,5X,2HJ=,I3)
 9030 FORMAT(5X,8HTHETA = ,G20.10/5X,12HSIN(THETA)= ,F10.7,
     +     15H   COS(THETA)= ,F10.7/5X,29HTOTAL CHANGE TO THIS POINT = ,
     +     G20.10)
 9050 FORMAT(10X,'ITERATION',I4,'   ORBITAL CHANGE=',G20.10)
 9060 FORMAT (/10X,
     +     114HLOCALIZATION HAS BEEN UNSUCESSFUL. PROGRAM WILL RESTART W
     +ITH NEW RANDOM NUMBER AND ROTATION SEQUENCE FOR ORBITALS.)
 9070 FORMAT(1H0,'LOCALIZATION FAILS --- LOCALIZATION ABORTED')
 9080 FORMAT(10X,'LOCALIZATION CONVERGED IN',I4,' ITERATIONS')
 9090 FORMAT(/10X,'THE PIPEK-MEZEY POPULATION LOCALIZED ORBITALS ARE')
 9120 FORMAT(//10X,'THE TRANSFORMATION FROM INITIAL MO-S (ROWS)',
     *             ' TO LOCALIZED MO-S (COLUMNS) IS')
 9130 FORMAT(//10X,'THIS LOCALIZATION HAD',I3,' M.O. S (',
     *       I3,' FROZEN) AND',I3,' BASIS FUNCTIONS')
 9140 FORMAT(10X,'THE FINAL LOCALIZATION SUM IS',F16.6/)
      END
C*MODULE LOCAL   *DECK POPQAT
      SUBROUTINE POPQAT(RIJ,POP,N2,NORB,NATOMS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION RIJ(N2,1),POP(NATOMS,NORB)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00)
      DATA UHF/8HUHF     /
C
C     COMPUTE MULLIKEN GROSS ATOMIC POPULATIONS
C
      INDX(I,J) = ((MAX(I,J)*(MAX(I,J)-1))/2+MIN(I,J))
      OCC = TWO
      IF (SCFTYP.EQ.UHF) OCC = ONE
      DO 10 I=1,NORB
      II = INDX(I,I)
      DO 10 J=1,NATOMS
   10 POP(J,I) = OCC*RIJ(II,J)
      IF (MASWRK) WRITE(IW,1000)
 1000 FORMAT(/'          MULLIKEN GROSS ATOMIC POPULATIONS')
      CALL PRSQ(POP,NORB,NATOMS,NATOMS)
      IF (MASWRK) WRITE(IW,1001)
 1001 FORMAT(/)
      RETURN
      END
C*MODULE LOCAL   *DECK POLNMM
      SUBROUTINE POLNMM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      LNA = NA - NOUTA
C
C     PARTITION FAST MEMORY
C
      CALL VALFM(LOADFM)
      LVA  = 1     + LOADFM
      LVX  = LVA  + L3
      LARR  = LVX  + L3
      LWRK  = LARR + L2
      LWRK2 = LWRK + L1
      LWRK3 = LWRK2+ L1
      LIWRK = LWRK3+ L3
      LCCH0 = LIWRK+ L1
      LCCHX = LCCH0+ 3*NA
      LPOL  = LCCHX+ 3*NA
      LNMOA0= LPOL + 9*LNA
      LMOID0= LNMOA0+ NA
      LXDIP = LMOID0+ 5*NA
      LYDIP = LXDIP + L2
      LZDIP = LYDIP + L2
      LWINT = LZDIP + L2
      LAST  = LWINT + 784*3
C
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL POLNMX(X(LVA),X(LVX),X(LARR),X(LCCH0),X(LCCHX),X(LWRK),
     *            X(LWRK2),X(LWRK3),X(LIWRK),X(LPOL),X(LNMOA0),
     *            X(LMOID0),X(LXDIP),X(LYDIP),X(LZDIP),X(LWINT),
     *            L0,L1,L2,LNA,NA)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE LOCAL   *DECK POLNMX
      SUBROUTINE POLNMX(VA,VX,ARRAY,CCHG,CCHX,WRK,WRK2,WRK3,IWRK,DLPOL,
     *                  NMOA0,MOID0,XDIP,YDIP,ZDIP,WINT,
     *                  L0,L1,L2,LNA,NORB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL EFLDL,MFRZ,MFRZBK
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,POLDCM,POLANG
      LOGICAL GOPARR,DSKWRK,MASWRK,POLAPP,KMIDPT,POLDYN
C
      PARAMETER (MXATM=2000, NMO=500, MXAO=8192)
      PARAMETER (ZERO=0.0D+00,THREE=3.0D+00)
      PARAMETER (UNIT=0.52917724924D+00)
      PARAMETER (SMFLD=1.0D-04)
C
      DIMENSION VA(L1,L1), VX(L1,L1), ARRAY(L2)
      DIMENSION CCHG(3,LNA), CCHX(3,LNA), WRK(L1), WRK2(L1), WRK3(L1,L1)
      DIMENSION IWRK(L1),DLPOL(9,LNA),NMOA0(NORB),MOID0(5,NORB),
     *          XDIP(L2),YDIP(L2),ZDIP(L2),WINT(784*3)
C
      COMMON /EDCMP / ZIJ(NMO),ZMO(5,NMO),OCCUP(NMO),DPFREQ(50),
     *                MOIDNO(5,NMO),IJMO(2,NMO),MOIJ(NMO),NMOIJ(NMO),
     *                NMOAT(NMO),NDPFREQ,IPROT(5),NPROT,
     *                MOIDON,EDCOMP,DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,
     *                POLANG,POLAPP,KMIDPT,POLDYN
      COMMON /EFLDC / EVEC(3),EFLDL
      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 /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PRPOPT/ ILOCAL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DIPMX,DIPMY,DIPMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      DATA RHF,RMC/8HRHF     ,8HMCSCF   /
      DATA RNONE/8HNONE    /
      DATA ELMOM/8HELMOM   /
C
C----THE DRIVER FOR CALCULATING LMO NUMERIC POLARIZABILITIES----
C    D.R.GARMER, W.J.STEVENS J.PHYS.CHEM. 93, 8263-8270(1989)
C
C   THE ZERO FIELD DIPOLES (IN AO BASIS) WERE CALCULATED BEFORE,
C   WE JUST READ THEM IN AND THEN TRANSFORM TO MO BASIS. THIS IS ALMOST
C   THE SAME AS IN THE PAULIX ROUTINE.
C
C GET LMOS
C
      CALL DAREAD(IDAF,IODA,VA,L1*NA,71,0)
C    SAVE LMOS ON DAF RECORD 280.
C    READ THE FOCK MATRIX, ALPHA AND BETA DENSITY MATRICES AND SAVE THEM
C    FOR THE LATER USE IN STONE ANALYSIS.
C
      CALL DAWRIT(IDAF,IODA,VA,L1*L1,280,0)
      CALL DAREAD(IDAF,IODA,ARRAY,L2,14,0)
      CALL DAWRIT(IDAF,IODA,ARRAY,L2,283,0)
      CALL DAREAD(IDAF,IODA,ARRAY,L2,16,0)
      CALL DAWRIT(IDAF,IODA,ARRAY,L2,281,0)
      IF(SCFTYP.EQ.RHF .OR. SCFTYP.EQ.RMC .OR. CITYP.NE.RNONE)GOTO 120
        CALL DAREAD(IDAF,IODA,ARRAY,L2,20,0)
        CALL DAWRIT(IDAF,IODA,ARRAY,L2,282,0)
 120  CONTINUE
C
      IF (NOUTA.GT.0) THEN
         DO I = 1,LNA
            CALL DCOPY(L1,VA(1,NOUTA+I),1,VA(1,I),1)
         END DO
      END IF
C
C  COPY ZERO FIELD MOIDNO AND NMOAT INTO AUXILIARY ARRAYS
C  TO BE USED SOLELY FOR THE PRINTOUT
C
      CALL ICOPY(NORB,NMOAT,1,NMOA0,1)
      CALL ICOPY(5*NORB,MOIDNO,1,MOID0,1)
C
C   CLEAN UP JUNK FROM THE "VIRTUAL SPACE" OF VA MATRIX
C   AND FILL IT WITH IDENTICAL UNIT VECTORS SO THAT
C   LATER THE OVERLAP CAN BE MORE SAFELY CALCULATED.
C
      CALL VCLR(VA(1,LNA+1),1,L1*(L1-LNA))
      DO I=LNA+1,L1
        VA(1,I)=1.0D+00
      END DO
C
C   CALCULATE CENTROIDS OF CHARGE AND SAVE THE ZERO ELECTRIC FIELD
C   CENTROIDS ON DAF RECORDS 252-254  WRITTEN LATER
C   BACK ON RECORDS 95-97, TO BE USED WITH THE STONE ANALYSIS
C   DURING THE MAKEFP RUN.
C
      DO I = 1,3
         CALL DAREAD(IDAF,IODA,ARRAY,L2,94+I,0)
         CALL DAWRIT(IDAF,IODA,ARRAY,L2,251+I,0)
         DO J = 1,LNA
            CALL TFTRI(CCHG(I,J),ARRAY,VA(1,J),WRK2,1,L1,L1)
         END DO
      END DO
C      WRITE(IW,*)'CCHG:'
C      DO J=1,LNA
C          WRITE(IW,*)J,(CCHG(I,J),I=1,3)
C      END DO
C
C   NOW WE JUST HAVE TO APPLY FIELD IN 3 DIRECTIONS
C
      NPRBAK=NPRINT
      NPUNBK=NPUNCH
      IF(MASWRK) WRITE(IW,9000)
      EFLDL=.TRUE.
      EVEC(1)=ZERO
      EVEC(2)=ZERO
      EVEC(3)=ZERO
C
      NPUNCH=0
      DO 100 IEF=1,3
      EVEC(IEF)=SMFLD
      IF(MASWRK) WRITE(IW,9010) (EVEC(I),I=1,3)
      NPRBAK=NPRINT
      NPRINT=-5
C
      CALL ENERGX
C
C        MAKE SURE DIPOLE INTEGRALS EXIST (EVEN IF NPRINT.EQ.-5).
C        ORIGINAL VERSION DIDN'T ATTEMPT TO CONTROL CENTER OF MOMENT
C        EVALUATION, BUT WHY NOT USE CENTER OF MASS?
C        NOTE, X,Y,Z DIPOLES MUST BE BACK TO BACK TO BACK STORAGE.
C
      CALL CALCOM(XP,YP,ZP)
      CALL PRCALC(ELMOM,XDIP,WINT,3,L2,.FALSE.)
      CALL DAWRIT(IDAF,IODA,XDIP,L2,95,0)
      CALL DAWRIT(IDAF,IODA,YDIP,L2,96,0)
      CALL DAWRIT(IDAF,IODA,ZDIP,L2,97,0)
C
C   LMOLOC USES UNDOCUMENTED NPRINT=-23 FOR MINIMAL PRINTOUT
C
      NPRINT=-23
      IF(SCFTYP.EQ.RHF) THEN
         CALL MOIDM
         CALL DIPLMM
      END IF
      IF (ILOCAL.EQ.1) CALL LMOLOC(ILOCAL)
      IF (ILOCAL.EQ.2) CALL LOCENG
      IF (ILOCAL.EQ.3) CALL LMOLOC(ILOCAL)
C
C GET LMOS
C
      CALL DAREAD(IDAF,IODA,VX,L1*NA,71,0)
      IF (NOUTA.GT.0) THEN
         DO I = 1,LNA
            CALL DCOPY(L1,VX(1,NOUTA+I),1,VX(1,I),1)
         END DO
      END IF
      CALL VCLR(VX(1,LNA+1),1,L1*(L1-LNA))
      DO I=LNA+1,L1
         VX(1,I)=1.0D+00
      END DO
C
C   CALCULATE CENTROIDS OF CHARGE
C
      DO I = 1,3
         CALL DAREAD(IDAF,IODA,ARRAY,L2,94+I,0)
         DO J = 1,LNA
            CALL TFTRI(CCHX(I,J),ARRAY,VX(1,J),WRK2,1,L1,L1)
CCC            CCHX(I,J) = CCHX(I,J) +XYZP(I)
         END DO
      END DO
C
C   KEEP THE ORDER OF MOS IN VX THE SAME AS IN VA -
C   IN CASE THE ORBITAL ORDER CHANGES UNDER THE EXTERNAL FIELD
C   WE HAVE TO SWITCH OFF MO FREEZING IN OVLSEL
C
      MFRZBK=.FALSE.
      IF(MFRZ) THEN
          MFRZBK=MFRZ
          MFRZ=.FALSE.
      END IF
      CALL DAREAD(IDAF,IODA,ARRAY,L2,12,0)
      CALL OVLSEL(VX,WRK2,VA,ARRAY,WRK3,WRK,IWRK,WRK,L0,L1,L2)
      IF(MFRZBK) MFRZ=MFRZBK
C
C   CENTROIDS NEED TO BE REORDERED TOO
C
      CALL ICOPY(L1,WRK,1,IWRK,1)
      CALL REORDR(CCHX,IWRK,LNA,3)
C      WRITE(IW,*)'REORDERED CCHX:'
C     DO J=1,LNA
C         WRITE(IW,*)J,(CCHX(I,J),I=1,3)
C     END DO
      IEFPR=(IEF-1)*3
      DO I=1,LNA
        DLPOL(IEFPR+1,I)=2*(CCHG(1,I)-CCHX(1,I))/SMFLD
        DLPOL(IEFPR+2,I)=2*(CCHG(2,I)-CCHX(2,I))/SMFLD
        DLPOL(IEFPR+3,I)=2*(CCHG(3,I)-CCHX(3,I))/SMFLD
      END DO
      EVEC(IEF)=ZERO
 100  CONTINUE
      EFLDL=.FALSE.
      CALL DAWRIT(IDAF,IODA,DLPOL,LNA*9,251,0)
C
C ----- PRINT OUT POLARIZABILITIES----
C
      IF (MASWRK) THEN
         IF(POLANG) THEN
            WRITE(IW,9020)
         ELSE
            WRITE(IW,9025)
         END IF
      END IF
      DO I=1,LNA
         DO 520 J=1,9
            IF(POLANG) THEN
               DLPOL(J,I)=DLPOL(J,I)*UNIT**THREE
            END IF
  520    CONTINUE
C
C
       ILO=I+NOUTA
       IF (MASWRK) THEN
            IF(NMOA0(ILO).EQ.1) WRITE(IW,9003) I, MOID0(1,ILO)
            IF(NMOA0(ILO).GE.2)
     *               WRITE(IW,9005) I, MOID0(1,ILO),MOID0(2,ILO)
            WRITE(IW,9050)
            WRITE(IW,9060) (DLPOL(J,I),J=1,3)
            WRITE(IW,9080) (DLPOL(J,I),J=4,6)
            WRITE(IW,9100) (DLPOL(J,I),J=7,9)
            ALPHAM=(DLPOL(1,I)+DLPOL(5,I)+DLPOL(9,I))/THREE
            WRITE(IW,9110)ALPHAM
        END IF
      END DO
      NPRINT=NPRBAK
      NPUNCH=NPUNBK
C
C   FILL THE DAF RECORDS 71, 14, 16, 20, 95-7 WITH THE ARRAYS SAVED
C   BEFORE, THEY WILL BE USED IN THE PAULIX ROUTINE
C
      CALL DAREAD(IDAF,IODA,VA,L1*L1,280,0)
      CALL DAWRIT(IDAF,IODA,VA,L1*L1,71,0)
      CALL DAREAD(IDAF,IODA,ARRAY,L2,283,0)
      CALL DAWRIT(IDAF,IODA,ARRAY,L2,14,0)
      CALL DAREAD(IDAF,IODA,ARRAY,L2,281,0)
      CALL DAWRIT(IDAF,IODA,ARRAY,L2,16,0)
      IF(SCFTYP.EQ.RHF .OR. SCFTYP.EQ.RMC .OR. CITYP.NE.RNONE) GOTO 130
         CALL DAREAD(IDAF,IODA,ARRAY,L2,282,0)
         CALL DAWRIT(IDAF,IODA,ARRAY,L2,20,0)
 130  CONTINUE
      DO I = 1,3
         CALL DAREAD(IDAF,IODA,ARRAY,L2,251+I,0)
         CALL DAWRIT(IDAF,IODA,ARRAY,L2,94+I,0)
      END DO
C
C
      RETURN
 9000 FORMAT(/1X,'*** NUMERIC COMPUTATION OF LOCALIZED ORBITAL ',
     *       'POLARIZABILITIES ***'/5X,
     *'D.R.GARMER AND W.J.STEVENS J.PHYS.CHEM. 93, 8263-8270(1989)')
 9003 FORMAT(//1X,' LMO NUMBER ',I3,' ALPHA POLARIZABILITY TENSOR ',
     *        'FOR CORE OR LONE PAIR ON ATOM',I3)
 9005 FORMAT(//1X,' LMO NUMBER ',I3,' ALPHA POLARIZABILITY TENSOR ',
     *        'FOR BOND BETWEEN ATOM',I3,' AND ATOM',I3)
 9010 FORMAT(/1X,'APPLIED ELECTRIC FIELD: X=',F7.4,' Y=',F7.4,
     *       ' Z=',F7.4/)
 9020 FORMAT(/10X,35(1H-)/
     *        10X,' LOCALIZED ALPHA POLARIZABILITIES '/
     *        10X,'       IN ANGSTROMS**3            '/
     *        10X,35(1H-)/)
 9025 FORMAT(//10X,35(1H-)/
     *        10X,' LOCALIZED ALPHA POLARIZABILITIES '/
     *        10X,'       IN ATOMIC UNITS            '/
     *        10X,35(1H-)/)
 9050 FORMAT(22X,'UX',13X,'UY',13X,'UZ')
 9060 FORMAT(10X,     ' UX ',3F15.9)
 9080 FORMAT(10X,     ' UY ',3F15.9)
 9100 FORMAT(10X,     ' UZ ',3F15.9)
 9110 FORMAT(/10X,     ' MEAN ALPHA POLARIZABILITY = ',3F15.9)
      END
C
C*MODULE LOCAL   *DECK DIRLMO
      SUBROUTINE DIRLMO
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA UHF/8HUHF     /, CHECK/8HCHECK   /
C
C   SUBROUTINE TO ORIENT A SET OF LOCALIZED ORBITALS:
C
C   1/  THE LOCALIZED ORBITALS ARE ANALYSED TO DETERMINE WHICH ARE
C       QUASI-ATOMIC.  THIS IS DONE BY DETERMINING THE MULLIKEN
C       POPULATIONS AND IS LIKE THE METHOD IN MOIDM WHICH WAS WRITTEN
C       BY JAN JENSEN.
C   2/  THESE QASI-AO'S ARE THEN DIRECTED BY SIMPLIFICATION OF THE
C       DENSITY MATRIX USING THE FOURTH POWER METHOD. ONLY ROTATIONS
C       BETWEEN ORBITALS ON THE SAME ATOM ARE CONSIDERED. ONLY THE
C       ALPHA ORBITALS ARE DIRECTED.
C   3/  ORIENTED DENSITY IS PRINTED OUT, FOR ALL ORBITALS NOT FROZEN
C       IN THE LOCALIZATION.
C   4/  BONDING ANALYSIS IS PERFORMED AND PRINTED OUT.
C   5/  ALL OCCUPIED ORBITALS ARE PRINTED OUT.
C   6/  ORBITALS AND DENSITY WRITTEN OUT TO .F10 AND .DAT.
C
C   REFERENCE FOR THE METHOD WILL COME SOON,
C   JOE IVANIC AND KLAUS RUEDENBERG.
C   CODE WRITTEN BY J. IVANIC.
C
C   CVGLOC USED AS LOCALIZATION CRITERION IS USED AS CRITERION IN
C          ORIENTATION.
C   MAXLOC USED FOR MAXIMUM NUMBER OF CYCLES IN LOCALIZATION METHODS
C          WHICH ARE NOT EDMISTON-RUEDENBERG.  IT IS USED HERE FOR
C          MAXIMUM NUMBER OF SWEEPS IN ORIENTATION.  DEFAULT VALUE
C          FOR MAXLOC IN LOCALIZATION METHODS IS 250.  IF DEFAULT
C          VALUE USED THEN DEFAULT VALUE FOR ORIENTATION IS 100,000
C          OTHERWISE USE READ IN VALUE.
C
C   READ IN:
C       71: LOCALIZED ORBITALS
C      285: LOCALIZED DENSITY
C   WRITTEN OUT:
C      286: ORIENTED ORBITALS
C      287: ORIENTED DENSITY
C
      NDAFL =  71
      NDADL = 285
      NDAGL = 286
      NDAHL = 287
C
      IF (MASWRK) WRITE(IW,9000)
      IF (MAXLOC.EQ.250) THEN
         MITER=100000
      ELSE
         MITER=MAXLOC
      ENDIF
      PRITOL = 0.02D+00
C
C  DETERMINE SIZES OF ORBITAL SPACES.
C
      IPASS = 1
      CALL LMOPSI(IPASS,MCORE,MDOC,MACT,NUMLOC)
      NOUT = NOUTA
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      M1 = NUMLOC - NOUT
      M2 = (M1*M1+M1)/2
      M3 = M1*M1
C
C  MEMORY REQUIREMENTS.
C
      CALL VALFM(LOADFM)
      LCLMO  = LOADFM + 1
      LRLMO  = LCLMO  + L3
      LSTRI  = LRLMO  + L3
      LTRAN  = LSTRI  + L2
      LDEN   = LTRAN  + M3
      LIWHI  = LDEN   + M2
      LINAT  = LIWHI  + M1 + NAT
      LIATB  = LINAT  + M1 + NAT
      LATMU  = LIATB  + (M1*NAT)/NWDVAR + 1
      LMAPT  = LATMU  + NAT
      LIATM  = LMAPT  + M1
      LAST   = LIATM  + (M1*NAT)/NWDVAR + 1
      NEED1 = LAST - LOADFM -1
      IF (MASWRK) WRITE(IW,9010) NEED1
      CALL GETFM(NEED1)
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C  PRINT OUT NUMBER OF OCCUPIED AND CORE ORBITALS.
C
      IF (MASWRK) WRITE(IW,9020) NUMLOC,MCORE
      IF (MASWRK.AND.SCFTYP.EQ.UHF) WRITE(IW,9025)
C
C  READ LOCALIZED ORBITALS AND LOCALIZED DENSITY.
C
      CALL DAREAD(IDAF,IODA,X(LCLMO),L3,NDAFL,0)
      CALL DAREAD(IDAF,IODA,X(LDEN),M2,NDADL,0)
C
C  PRINT OUT THE LIST OF FROZEN ORBITALS. THESE WILL NOT
C  BE INCLUDED IN THE ORIENTATION ANALYSIS.
C
      IF (MASWRK) THEN
         WRITE(IW,9030) NOUTA,(MOOUTA(I),I=1,NOUTA)
         WRITE(IW,9040)
      ENDIF
C
C  REORDER ORBITALS SUCH THAT FROZEN ONES OCCUR FIRST.
C
      CALL GETMAP(NUMLOC,X(LCLMO),X(LRLMO),L1)
      NATS = NAT
C
C  IDENTIFY WHICH ORBITALS ARE QUASI-ATOMIC AND REORDER SO THAT
C  THESE OCCUR FIRST IN THE UNFROZEN SET.
C
      CALL DIMOID(X(LDEN),X(LRLMO),X(LCLMO),X(LSTRI),X(LATMU),X(LIATM),
     *    X(LIWHI),X(LMAPT),X(LINAT),X(LIATB),L1,L2,M1,M2,NATS,NOSI,
     *    NCAT,NSWE)
C
C  NOW TO PRINT OUT THE CONVERGENCE CRITERION FOR THE ANGLE, AND
C  THEN PERFORM THE ORIENTATION PROCEDURE.
C
      IF (MASWRK) WRITE(IW,9050) M1-NOSI,CVGLOC,NSWE
      CALL ORIEN(IW,NCAT,M1,X(LINAT),X(LDEN),X(LTRAN),CVGLOC,MITER)
C
C  NOW TO PRINT OUT THE DENSITY MATRIX FOR NON-FROZEN ORBITALS.
C
      IF (MASWRK) WRITE(IW,9060) PRITOL
      CALL ORIPRI(IW,NCAT,NATS,M1,X(LINAT),X(LDEN),X(LIWHI),X(LIATB),
     *  NOSI,PRITOL)
C
C  PRINT OUT THE BONDING ANALYSIS.
C
      IF (MASWRK) WRITE(IW,9065)
      CALL ORIANAL(IW,NCAT,NATS,M1,X(LINAT),X(LDEN),X(LIATB),NOSI)
C
C  MAKE NEW ORBITAL BASIS COEFFICIENTS.
C
      CALL MRARBR(X(LRLMO+NOUTA*L1),L1,L1,M1,X(LTRAN),M1,M1,
     *            X(LCLMO+NOUTA*L1),L1)
C
C  PRINT OUT FINAL ORIENTED ORBITALS TO GAMESS OUTPUT FILE..
C
      IF (MASWRK) THEN
         WRITE(IW,9070)
         CALL PRSQL(X(LCLMO),NUMLOC,L1,L1)
      ENDIF
C
C  PRINT OUT ORBITALS TO .DAT FILE.
C
      IF(MASWRK) THEN
         WRITE(IP,9080)
         CALL PUSQL(X(LCLMO),NUMLOC,L1,L1)
         WRITE(IP,9085)
      ENDIF
C
C  WRITE OUT ORIENTED ORBITALS TO NDAGL = 286
C
      CALL DAWRIT(IDAF,IODA,X(LCLMO),L1*L1,NDAGL,0)
C
C  WRITE OUT ORIENTED DENSITY TO NDAHL = 287
C
      CALL DAWRIT(IDAF,IODA,X(LDEN),M2,NDAHL,0)
C
C  FINISH UP.
C
  800 CONTINUE
      CALL RETFM(NEED1)
      IF (MASWRK) WRITE(IW,9075)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/10X,50(1H-)/
     * 10X,'IVANIC-RUEDENBERG QUASI-ATOMIC ORBITAL ORIENTATION'/
     * 10X,50(1H-))
 9010 FORMAT(/1X,'THIS ORIENTATION REQUIRES',I10,' WORDS OF MEMORY.')
 9020 FORMAT(/1X,'THIS WAVEFUNCTION HAS',I4,' OCCUPIED ORBITALS',
     *       ' OF WHICH',I4,' ARE CORE.')
 9025 FORMAT(1X,'ONLY ALPHA LOCALIZED ORBITALS WILL BE ORIENTED')
 9030 FORMAT(1X,'THERE WERE',I3,' ORBITALS FROZEN IN THE ',
     * 'LOCALIZATION.',/1X,'THESE ARE: ',100I3)
 9040 FORMAT(1X,'THEY WILL BE MOVED TO THE BEGINNING OF THE ',
     * 'ORBITAL LIST'/1X,'AND OMMITTED FROM THE ORIENTATION ANALYSIS.')
 9050 FORMAT(/1X,'***********************************************'
     *       /1X,'    ORIENTING FIRST ',I3,' NON-FROZEN ORBITALS    '
     *       /1X,'***********************************************'
     *  //1X,'ROTATION ANGLES WILL BE CONVERGED TO ',E9.2,
     *  /1X,'NUMBER OF ANGLES IN A SWEEP = ',I6)
 9060 FORMAT(
     *   /1X,'----------------------------------------------------'
     *   /1X,'DENSITY MATRIX FOR QUASI-ATOMIC AND BONDING ORBITALS'
     *   /1X,'ONLY ELEMENTS BIGGER THAN ',F4.2,' ARE PRINTED'
     *   /1X,'----------------------------------------------------')
 9065 FORMAT(
     *   /1X,'----------------'
     *   /1X,'BONDING ANALYSIS'
     *   /1X,'----------------')
 9070 FORMAT(//10X,'ORIENTED LOCALIZED ORBITALS')
 9075 FORMAT(1X,'... DONE WITH ORIENTATION ...')
 9080 FORMAT('ORIENTED LOCALIZED ORBITALS'/' $VEC')
 9085 FORMAT(' $END')
C
      END
C
C*MODULE LOCAL   *DECK GETMAP
      SUBROUTINE GETMAP(NUMLOC,CLMO,RLMO,L1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION CLMO(L1,L1),RLMO(L1,L1)
C
      PARAMETER (MXAO=8192)
C
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
C
C  SUBROUTINE TO DETERMINE WHICH ALPHA ORBITALS HAVE BEEN LOCALIZED,
C  AND REORDER THE ORBITALS SO THAT FROZEN ONES OCCUR FIRST, AND
C  VIRTUAL LAST.  THIS HAS BEEN RIPPED OUT OF SUBROUTINE LOCSET.
C
      NOUT = NOUTA
      CALL VCLR(RLMO,1,L1*L1)
C
C     PUT THE ORBITALS WHICH ARE TO BE FROZEN FIRST...
C
      NMO=0
      DO 120 J=1,NUMLOC
         DO 110 I=1,NOUT
            IF(MOOUTA(I).EQ.J) THEN
                  NMO = NMO+1
                  CALL DCOPY(L1,CLMO(1,J),1,RLMO(1,NMO),1)
                  GO TO 120
            END IF
  110    CONTINUE
  120 CONTINUE
C
C     ...THEN PUT IN THE ORBITALS TO BE LOCALIZED...
C
      LL=0
      DO 220 J=1,NUMLOC
         DO 210 I=1,NOUT
            IF(MOOUTA(I).EQ.J) GO TO 220
  210    CONTINUE
         NMO = NMO+1
         CALL DCOPY(L1,CLMO(1,J),1,RLMO(1,NMO),1)
         LL = LL+1
  220 CONTINUE
C
C     ...AND FINISH UP WITH THE VIRTUAL SPACE.
C
      LVIRT = (L1-NUMLOC)*L1
      CALL DCOPY(LVIRT,CLMO(1,NMO+1),1,RLMO(1,NMO+1),1)
C
      RETURN
      END
C
C*MODULE LOCAL   *DECK DIMOID
      SUBROUTINE DIMOID(DEN,RLMO,SSQU,STRI,ATMU,IATM,IWHI,MAPT,INAT,
     *   IATB,L1,L2,M1,M2,NATS,NOSI,NCAT,NSWE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION RLMO(L1,L1),SSQU(L1,L1),STRI(L2),ATMU(NATS),DEN(M2)
      DIMENSION IATM(NATS,M1),IWHI(M1+NATS),MAPT(M1),INAT(M1+NATS)
      DIMENSION IATB(NATS,M1)
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000, MXAO=8192)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      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 /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
C
C   SUBROUTINE TO DETERMINE ATOMIC MULLIKEN POPULATIONS FOR
C   EACH NON-FROZEN LOCALIZED MO.
C   LMO'S ARE THEN SORTED AS EITHER QUASI-ATOMIC OR BOND MO'S.
C   QUASI-ATOMIC ORBITALS ARE PLACED FIRST AND BOND MO'S PLACED
C   LAST IN THE NON-FROZEN LMO LIST.
C
      TOLR = 0.30D+00
      NUMLOC = M1+NOUTA
C
      CALL DAREAD(IDAF,IODA,STRI,L2,12,0)
      CALL CPYTSQ(STRI,SSQU,L1,1)
C
      DO 920 II=1,M1
         INAT(II) = 0
  920 CONTINUE
C
      IF (MASWRK) WRITE(IW,8000)
      NST = 1
 2000 CONTINUE
      NEND = MIN(NAT,NST+8)
C
      IF (MASWRK) THEN
         WRITE(IW,9000) (I,A(I),I=NST,NEND)
         WRITE(IW,9005)
      ENDIF
C
      DO 900 IO = NOUTA+1,NUMLOC
         IZ = IO - NOUTA
         DO 895 II=NST,NEND
            ATMU(II) = 0.0D+00
            IATM(II,IZ) = 0
  895    CONTINUE
         IFUNC = 0
         DO 890 ISHELL = 1,NSHELL
            IAT = KATOM(ISHELL)
            IST = KMIN(ISHELL)
            IEN = KMAX(ISHELL)
            DO 880 INO = IST,IEN
               IFUNC = IFUNC + 1
               IF (IAT.LT.NST.OR.IAT.GT.NEND) GOTO 880
               ZINT  = 0.0D+00
               DO 870 II = 1,L1
                  ZINT = ZINT + RLMO(II,IO)*SSQU(II,IFUNC)
  870          CONTINUE
               ATMU(IAT) = ATMU(IAT) + RLMO(IFUNC,IO)*ZINT
  880       CONTINUE
  890    CONTINUE
         IF (MASWRK) WRITE(IW,9010) IZ,(ATMU(II),II=NST,NEND)
         DO 790 IA=NST,NEND
            IF (ATMU(IA).GT.TOLR) THEN
               INAT(IZ) = INAT(IZ) + 1
               IND = INAT(IZ)
               IATM(IND,IZ) = IA
            ENDIF
  790    CONTINUE
  900 CONTINUE
C
      NST = NST + 9
      IF (NEND.LT.NAT) GOTO 2000
C
      IF (MASWRK) THEN
         WRITE(IW,9015) TOLR
         WRITE(IW,9020)
      ENDIF
C
      NOSI = 0
      DO 700 II=1,M1
         NO=0
         DO 720 JJ=1,NAT
            IF (IATM(JJ,II).EQ.0) GOTO 740
            NO = NO + 1
  720    CONTINUE
  740    CONTINUE
         IF (NO.GT.1.OR.NO.EQ.0) THEN
            NOSI = NOSI + 1
            IWHI(NOSI) = II
         ENDIF
        IF (MASWRK)
     *     WRITE(IW,9030) II,(IATM(J,II),A(IATM(J,II)),J=1,NO)
  700 CONTINUE
C
      IF (MASWRK) THEN
         WRITE(IW,9035) NOSI
         IF (NOSI.GT.0) THEN
            WRITE(IW,9040) (IWHI(I),I=1,NOSI)
            WRITE(IW,9045)
         ELSE
            WRITE(IW,9050)
         ENDIF
      ENDIF
C
      CALL DCOPY(L1*L1,RLMO,1,SSQU,1)
      CALL DCOPY(M2,DEN,1,STRI,1)
C
      IP2 = NOUTA
      IS2 = M1+NOUTA-NOSI
      DO 695 II=1,NAT
         INAT(II) = 0
  695 CONTINUE
C
      DO 690 IAT=1,NAT
         DO 680 IORB=1,M1
            IP1 = IORB + NOUTA
            IF (IATM(1,IORB).NE.IAT) GOTO 680
            IF (IATM(2,IORB).NE.0) GOTO 680
            INAT(IAT) = INAT(IAT) + 1
            IP2 = IP2 + 1
            CALL DCOPY(L1,SSQU(1,IP1),1,RLMO(1,IP2),1)
            CALL ICOPY(NAT,IATM(1,IORB),1,IATB(1,IP2-NOUTA),1)
            MAPT(IORB) = IP2-NOUTA
  680    CONTINUE
         DO 670 IORB=1,NOSI
            IS1 = IWHI(IORB) + NOUTA
            IF (IAT.EQ.NAT.AND.IATM(1,IWHI(IORB)).EQ.0) GOTO 675
            IF (IATM(1,IWHI(IORB)).NE.IAT) GOTO 670
  675       CONTINUE
            IS2 = IS2 + 1
            CALL DCOPY(L1,SSQU(1,IS1),1,RLMO(1,IS2),1)
            CALL ICOPY(NAT,IATM(1,IS1-NOUTA),1,IATB(1,IS2-NOUTA),1)
            MAPT(IWHI(IORB)) = IS2-NOUTA
  670    CONTINUE
  690 CONTINUE
C
      IJ=0
      DO 650 I=1,M1
         IMAP = MAPT(I)
         DO 640 J=1,I
            IJ=IJ+1
            JMAP = MAPT(J)
            IN1 = MAX(IMAP,JMAP)
            IN2 = MIN(IMAP,JMAP)
            IJMAP = (IN1*IN1-IN1)/2 + IN2
            DEN(IJMAP) = STRI(IJ)
  640    CONTINUE
  650 CONTINUE
C
      DO 630 II=1,NAT
         IWHI(II) = INAT(II)
  630 CONTINUE
C
      NSWE = 0
      NCAT = 0
      LASP = 1
      NLAST = 0
      DO 620 II=1,NAT
         NSWE = NSWE + (IWHI(II)*(IWHI(II)-1))/2
         IF (IWHI(II).EQ.0) GOTO 620
         NCAT = NCAT + 1
         INAT(NCAT) = LASP + NLAST
         LASP = INAT(NCAT)
         NLAST = IWHI(II)
         IWHI(NCAT) = II
  620 CONTINUE
C
      DO 610 II=1,NOSI
         NCAT = NCAT + 1
         INAT(NCAT) = LASP + NLAST
         LASP = INAT(NCAT)
         NLAST = 1
         IWHI(NCAT) = 0
  610 CONTINUE
C
      RETURN
C
 8000 FORMAT(/1X,'** MULLIKEN ATOMIC POPULATIONS FOR EACH NON-FROZEN ',
     *       'LOCALIZED ORBITAL **')
 9000 FORMAT(/3X,'ATOM',2X,100(I2,1X,A4))
 9005 FORMAT(1X,'LMO')
 9010 FORMAT(1X,I3,3X,100F7.3)
 9015 FORMAT(/1X,'** ATOMIC POPULATIONS GREATER THAN ',F4.2,
     *   ' ARE CONSIDERED MAJOR **')
 9020 FORMAT(/2X,'LMO',3X,'MAJOR CONTRIBUTIONS FROM ATOM(S)')
 9030 FORMAT(2X,I3,2X,100(I2,1X,A2,2X))
 9035 FORMAT(/1X,'NO OF LMOS INVOLVING MORE THAN ONE ATOM =',I3)
 9040 FORMAT(1X,'THESE ARE LMOS :',100I3)
 9045 FORMAT(
     *1X,'THESE WILL BE MOVED TO THE END OF THE NON-FROZEN LMO LIST'/
     *1X,'AND REMOVED FROM THE ORIENTATION PROCEDURE.'/
     *1X,'REMAINING ORBITALS WILL BE REORDERED ACCORDING TO WHICH ATOM'/
     *1X,'THEY BELONG TO, I.E. LMOS ON ATOM 1 COME FIRST, ETC..')
 9050 FORMAT(1X,
     * 'LMOS WILL BE REORDERED ACCORDING TO WHICH ATOM THEY BELONG TO'/
     * 1X,'I.E. LMOS ON ATOM 1 COME FIRST, ETC..')
C
      END
C
C*MODULE LOCAL   *DECK ORIEN
C     --------------------------------------------------------------
      SUBROUTINE ORIEN(IW,NATOT,NTOTORB,NATORB,P,T,CRIT,MXITER)
C     ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION NATORB(NATOT),P(NTOTORB*(NTOTORB+1)/2)
      DIMENSION T(NTOTORB,NTOTORB)
C
C    DIRECTS ATOMIC ORBITALS SUCH THAT THEY ARE BOND-ORIENTED.
C    TAKES A DENSITY MATRIX, EXPRESSED IN TERMS OF ATOMIC ORBITALS,
C    AND REFINES IT AS FOLLOWS:
C    ORBITALS ON EACH ATOM ARE ROTATED SUCH THAT THE OFF-DIAGONAL
C    INTER-ATOMIC BLOCKS CONTAIN 1) AS FEW LARGE ELEMENTS AS POSSIBLE,
C    AND 2) AS MANY NEAR-ZERO ELEMENTS AS POSSIBLE.
C
C    *  IW    = UNIT WHERE ANY OUTPUT IS WRITTEN TO.
C    *  NATOT = NUMBER OF CATEGORIES OF ORBITALS.
C    *  NTOTORB = TOTAL NUMBER OF ORBITALS.
C    *  NATORB(I) = WHERE THE VALENCE ATOMIC ORBS OF CATEGORY I BEGIN.
C    *  P IS DENSITY MATRIX, STORED IN LOWER TRIANGLE FORM.
C       ORBITALS SHOULD OCCUR AT THE BEGINNING.
C    *  T IS THE RETURNED TRANSFORMATION MATRIX FOR THE ORBITALS.
C    *  CRIT = ROTATION EXECUTION THRESHOLD.  IE. FOR A PAIR OF
C       ORBITALS TO BE UPDATED, THE ABSOLUTE VALUE OF THE ANGLE MUST
C       BE LARGER THAN CRIT (I LIKE THE VALUE OF 1.0D-06)
C    *  MXITER = MAXIMUM NUMBER OF SWEEPS.  THIS CODE FLIES ALONG
C       SO MAKE IT BIG, E.G. 100,000 OR MAGNITUDES BIGGER.
C
      ITER = 0
      CALL INITT(T,NTOTORB)
C
C    LOOP OVER ALL ATOMS
C
  100 CONTINUE
      ITER = ITER + 1
      NUMROT = 0
      DO 9000 IATOM=1,NATOT
C
C     LOOP OVER ALL PAIRS OF AOS ON ATOM IATOM.
C
         ILAST = NTOTORB
         IF (IATOM.NE.NATOT) ILAST=NATORB(IATOM+1)-1
C
         DO 8000 IAOI=NATORB(IATOM),ILAST
            DO 7000 IAOJ = IAOI+1,ILAST
C
               R2 = 0.0D+00
               R3 = 0.0D+00
C
               DO 6000 INOTA=1,NATOT
                  IF (INOTA.EQ.IATOM) GOTO 6000
                  ILAST2 = NTOTORB
                  IF (INOTA.NE.NATOT) ILAST2=NATORB(INOTA+1)-1
C
                  DO 5000 IK=NATORB(INOTA),ILAST2
                     IMAI=MAX(IK,IAOI)
                     IMII=MIN(IK,IAOI)
                     IMAJ=MAX(IK,IAOJ)
                     IMIJ=MIN(IK,IAOJ)
                     IKI=(IMAI*(IMAI-1))/2 + IMII
                     IKJ=(IMAJ*(IMAJ-1))/2 + IMIJ
                     PIKI=P(IKI)
                     PIKJ=P(IKJ)
C
       R2 = R2 + (PIKI**4)-6*(PIKI*PIKI*PIKJ*PIKJ)+(PIKJ**4)
       R3 = R3 + ((PIKI**3)*PIKJ) - ((PIKJ**3)*PIKI)
C
 5000             CONTINUE
 6000          CONTINUE
C
               R2 = (R2/4.0D+00)
               Q = SQRT(R2*R2 + R3*R3)
               IF (Q.LT.1.0D-08) GO TO 7000
               THETA = (ATAN2((R3/Q),(R2/Q)))/4.0D+00
C
C    CHECK IF THE ROTATION IS WORTH EXECUTING.
C
               IF (ABS(THETA).LT.CRIT) GOTO 7000
C
C    NOW HAVE THE ANGLE, THETA, HAVE TO ROTATE DENSITY MATRIX
C    AND UPDATE THE TRANSFORMATION T.
C
               NUMROT = NUMROT + 1
               A = COS(THETA)
               B = -SIN(THETA)
               CALL ROT1INT(NTOTORB,IAOI,IAOJ,A,B,P)
               CALL ADJTR(IAOI,IAOJ,NTOTORB,A,B,T)
C
 7000       CONTINUE
 8000    CONTINUE
 9000 CONTINUE
C
      IF (ITER.GT.MXITER) THEN
         WRITE(IW,*) '*************************************************'
         WRITE(IW,*) 'ORBITALS NOT CONVERGED AFTER ',MXITER,' ITERS.'
         WRITE(IW,*) '*************************************************'
         CALL ABRT
      ENDIF
C
C  IF THERE WERE ANY ROTATIONS, DO ANOTHER SWEEP.
C
      IF (NUMROT.GT.0) GOTO 100
C
      WRITE(IW,'(/A,I7)')
     * ' -- NUMBER OF SWEEPS REQUIRED FOR CONVERGENCE ='
     *     ,ITER
C
C  NOW ORDER THE ATOMIC ORBITALS ACCORDING TO THEIR OCCUPATION NUMBER.
C
      DO 9100 IATOM=1,NATOT
         ILAST = NTOTORB
         IF (IATOM.NE.NATOT) ILAST=NATORB(IATOM+1)-1
         IF (ILAST-NATORB(IATOM).EQ.0) GOTO 9100
         DO 9150 IO1=NATORB(IATOM),ILAST
            IMAX = IO1
            INDX = LTPOS(IO1,IO1)
            PMAX = P(INDX)
            DO 9160 IO2=IO1+1,ILAST
               IND = LTPOS(IO2,IO2)
               IF (P(IND).GT.PMAX) THEN
                  PMAX = P(IND)
                  IMAX = IO2
               ENDIF
 9160       CONTINUE
            IF (IMAX.EQ.IO1) GOTO 9150
            DO 9170 II=1,NTOTORB
               TEMP = T(II,IO1)
               T(II,IO1) = T(II,IMAX)
               T(II,IMAX) = TEMP
               IF (II.EQ.IO1.OR.II.EQ.IMAX) GOTO 9170
               IND1 = LTPOS(IO1,II)
               IND2 = LTPOS(IMAX,II)
               TEMP = P(IND1)
               P(IND1) = P(IND2)
               P(IND2) = TEMP
 9170       CONTINUE
            IND1 = LTPOS(IO1,IO1)
            IND2 = LTPOS(IMAX,IMAX)
            TEMP = P(IND1)
            P(IND1) = P(IND2)
            P(IND2) = TEMP
 9150    CONTINUE
 9100 CONTINUE
C
C  NOW TRY TO MAKE AS MANY OFF-DIAGONAL ELEMENTS AS POSSIBLE POSITIVE.
C
C
      DO 9200 IATOM=1,NATOT
         ILAST = NTOTORB
         IF (IATOM.NE.NATOT) ILAST = NATORB(IATOM+1)-1
         DO 9210 IO1=NATORB(IATOM),ILAST
            IMAX = 1
            ZMAX = 0.0D+00
            DO 9230 IO2=IO1+1,NTOTORB
               DO 9240 II=NATORB(IATOM),ILAST
                  IF (IO2.EQ.II) GOTO 9230
 9240          CONTINUE
               IND1 = LTPOS(IO1,IO2)
               IF (ABS(P(IND1)).GT.ZMAX) THEN
                  IOZ = IO2
                  IMAX = IND1
                  ZMAX = ABS(P(IND1))
               ENDIF
 9230       CONTINUE
            IF (P(IMAX).LT.0.0D+00) THEN
               DO 9250 IO2=1,NTOTORB
                  IPOS = LTPOS(IO2,IOZ)
                  T(IO2,IOZ) = -T(IO2,IOZ)
                  P(IPOS) = -P(IPOS)
 9250          CONTINUE
               IPOS = LTPOS(IOZ,IOZ)
               P(IPOS) = ABS(P(IPOS))
            ENDIF
 9210    CONTINUE
 9200 CONTINUE
      RETURN
      END
C
C*MODULE LOCAL   *DECK INITT
C     ---------------------
      SUBROUTINE INITT(T,N)
C     ---------------------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION T(N,N)
C
      DO 13 I=1,N
         DO 113 J=1,N
            T(I,J) = 0.0D+00
  113    CONTINUE
   13 CONTINUE
      DO 14 I=1,N
         T(I,I) = 1.0D+00
   14 CONTINUE
C
      RETURN
      END
C
C*MODULE LOCAL   *DECK ROT1INT
C     --------------------------------------------------------
      SUBROUTINE ROT1INT(N,I,J,A,B,CL)
C     --------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,T-Z)
      IMPLICIT INTEGER(I-S)
      DIMENSION CL(1)
C
C     APPLICATION OF JACOBI TRANSFORMATION TO A CANONICAL
C     ARRAY OF 1-ELECTRON INTEGRALS.
C     IE FOR <P|Q> WHERE P.GE.Q
C
C     N       ORBITAL DIMENSION OF THE INTEGRAL MATRIX
C     I,J     ORBITALS TO BE ROTATED, I<J
C     A,B     TOP LEFT AND RIGHT ELEMENTS, RESPECTIVELY, OF THE JACOBI
C             TRANSFORMATION MATRIX
C     CL      ARRAY CONTAINING INTEGRALS.
C
C     CONSTANTS REQUIRED FOR FULL TRANSFORMATION.
C
      A2 = A*A
      AB = A*B
      B2 = B*B
      T1 = 2.0D+00*AB
      T2 = A2 - B2
C
C     END OF CONSTANTS
C
C     GROUP A/  BOTH ORBITALS IN INTEGRAL ARE EITHER I OR J.
C
      II = I*(I+1)/2
      JJ = J*(J+1)/2
      JI = JJ - J + I
C
C    UP TO HERE OKAY
C
      C1 = CL(II)
      C2 = CL(JI)
      C3 = CL(JJ)
C
      CL(II) = A2*C1 - T1*C2 + B2*C3
      CL(JI) = AB*C1 + T2*C2 - AB*C3
      CL(JJ) = B2*C1 + T1*C2 + A2*C3
C
C    GROUP B/  ONE ORBITAL IN INTEGRAL IS EITHER I OR J.
C
      DO 13 K=1,N
         IF (K.EQ.I.OR.K.EQ.J) GOTO 13
         IF (K.LT.I) THEN
            I1 = II - I + K
            I2 = JJ - J + K
         ELSE IF(K.GT.I.AND.K.LT.J) THEN
            I1 = K*(K-1)/2 + I
            I2 = JJ - J + K
         ELSE
            KX = K*(K-1)/2
            I1 = KX + I
            I2 = KX + J
         ENDIF
C
         C1 = CL(I1)
         C2 = CL(I2)
         CL(I1) = A*C1 - B*C2
         CL(I2) = B*C1 + A*C2
C
   13 CONTINUE
C
      RETURN
      END
C
C
C*MODULE LOCAL   *DECK ORIPRI
C     --------------------------------------------------------------
      SUBROUTINE ORIPRI(IW,NATOT,NATS,NTOTORB,NATORB,P,IWRK2,IATB,
     *    NOSI,CRIT)
C     ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION NATORB(NATOT),P(NTOTORB*(NTOTORB+1)/2)
      DIMENSION IWRK2(NATOT),IATB(NATS,NTOTORB)
      CHARACTER*3 MOL
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
C
C     PRINTS 1-PARTICLE DENSITY MATRIX.
C
C    *  IW IS UNIT NUMBER FOR WRITING OUTPUT.
C    *  NATOT = NUMBER OF CATEGORIES OF ORBITALS.
C    *  NATS  = NUMBER OF ACTUAL ATOMS.
C    *  NTOTORB = TOTAL NUMBER OF ORBITALS.
C    *  NATORB(I) IS WHERE THE ORBITALS OF CATEGORY I BEGIN.
C    *  P IS DENSITY MATRIX, STORED IN LOWER TRIANGLE FORM.
C    *  IWRK2 IS A SCRATCH INTEGER ARRAY.
C    *  IATB(I,J) CONTAINS INFORMATION ABOUT WHICH ATOMS ARE INVOLVED
C       IN ORBITAL J.
C    *  NOSI IS THE NUMBER OF ORBITALS THAT ARE NOT ATOMIC IN NATURE.
C    *  CRIT = ELEMENTS WITH ABS. VALUE GREATER THAN CRIT ARE PRINTED.
C
      NOAO = NTOTORB - NOSI
C
      DO 41 II=1,NATOT
         IEND = NTOTORB
         IF (II.LT.NATOT) IEND = NATORB(II+1)-1
         IWRK2(II) = IEND-NATORB(II)+1
   41 CONTINUE
C
      NST = 1
C
   10 CONTINUE
      ICOL = 8
      WRITE(IW,8050)
C
      DO 13 II=NST,NATOT
         IF (NATORB(II).GT.NOAO.
     *   AND.NATORB(II-1).LE.NOAO.AND.(II-NST).GT.0) GOTO 14
         NY1 = IWRK2(II)
         ICOL = ICOL + 6*(NY1-1)+8
         IF (ICOL.GT.85) GOTO 14
         ICOL = ICOL + 1
   13 CONTINUE
C
   14 IAT = II-1
C
      DO 113 IQ = NST,IAT
         IF (IWRK2(IQ).LT.1) GOTO 113
         IF (NATORB(IQ).LE.NOAO) THEN
            WRITE(IW,9000) 'ATOM',IATB(1,NATORB(IQ))
         ELSE
            WRITE(IW,9000) 'BOND',NATORB(IQ)-NOAO
         ENDIF
         IF (IWRK2(IQ).EQ.1) THEN
            WRITE(IW,9003)
            GOTO 113
         ENDIF
         DO 32 IO = 2,IWRK2(IQ)
            WRITE(IW,9006)
   32    CONTINUE
         WRITE(IW,9003)
  113 CONTINUE
C
      WRITE(IW,9020)
      DO 313 IQ = NST,IAT
         IF (IWRK2(IQ).LT.1) GOTO 313
         IF (IWRK2(IQ).EQ.1) THEN
            WRITE(IW,9030) NATORB(IQ)
            GOTO 313
         ENDIF
         DO 332 IO=1,IWRK2(IQ)
            WRITE(IW,9032) NATORB(IQ)+IO-1
  332    CONTINUE
         WRITE(IW,9004)
  313 CONTINUE
C
      IF (NATORB(IAT).LE.NOAO) THEN
         WRITE(IW,8060)
      ELSE
         WRITE(IW,8061)
      ENDIF
      DO 413 IQ = NST, IAT
         IF (IWRK2(IQ).LT.1) GOTO 413
         IF (IWRK2(IQ).EQ.1) THEN
            IF (NATORB(IQ).LE.NOAO) THEN
               WRITE(IW,9037) A(IATB(1,NATORB(IQ))),
     *     IATB(1,NATORB(IQ))
            ELSE
               IF (IATB(1,NATORB(IQ)).EQ.0) WRITE(IW,7039)
               IF (IATB(1,NATORB(IQ)).NE.0.AND.
     *             NATS.GT.2.AND.
     *             IATB(3,NATORB(IQ)).EQ.0)
     *             WRITE(IW,7037)
     *  IATB(1,NATORB(IQ)),IATB(2,NATORB(IQ))
               IF (IATB(1,NATORB(IQ)).NE.0.AND.
     *             NATS.LE.2)
     *             WRITE(IW,7037)
     *  IATB(1,NATORB(IQ)),IATB(2,NATORB(IQ))
               IF (IATB(1,NATORB(IQ)).NE.0.AND.
     *             NATS.GT.2.AND.
     *             IATB(3,NATORB(IQ)).NE.0)
     *             WRITE(IW,7038)
     *  IATB(1,NATORB(IQ)),IATB(2,NATORB(IQ))
            ENDIF
            GOTO 413
         ENDIF
         DO 432 IO=1,IWRK2(IQ)
            WRITE(IW,9039) A(IATB(1,NATORB(IQ))),IO
  432    CONTINUE
         WRITE(IW,9005)
  413 CONTINUE
C
      WRITE(IW,7062)
      DO 513 IQ = NST, IAT
         IF (IWRK2(IQ).LT.1) GOTO 513
         IF (IWRK2(IQ).EQ.1) THEN
            WRITE(IW,7040)
            GOTO 513
         ENDIF
         DO 532 IO=1,IWRK2(IQ)
            WRITE(IW,7042)
  532    CONTINUE
         WRITE(IW,7044)
  513 CONTINUE
C
      DO 613 IQA = 1,NATOT
         IF (IWRK2(IQA).LT.1) GOTO 613
         DO 676 IQO = 1,IWRK2(IQA)
            ICT = NATORB(IQA)+IQO-1
            IF (NATORB(IQA).LE.NOAO) THEN
            WRITE(IW,9052) ICT,A(IATB(1,NATORB(IQA))),
     *      IATB(1,NATORB(IQA))
            ELSE
            WRITE(IW,9052) ICT,'BD',ICT-NOAO
            ENDIF
C
            DO 713 IQ = NST,IAT
               IF (IWRK2(IQ).LT.1) GOTO 713
               IF (IWRK2(IQ).EQ.1) THEN
                  IJT = NATORB(IQ)
                  IPOS = LTPOS(ICT,IJT)
                  IF (ABS(P(IPOS)).LT.CRIT) THEN
                     MOL = ' | '
                     WRITE(IW,9053) MOL
                     GOTO 713
                  ENDIF
                  IF (ABS(P(IPOS)).LT.0.3D+00) THEN
                     MOL = ' | '
                  ELSE
                     MOL = '#| '
                  ENDIF
                  IF (ICT.EQ.IJT) MOL = '*| '
                  WRITE(IW,9054) P(IPOS),MOL
                  GOTO 713
               ENDIF
               DO 732 IO = 1,IWRK2(IQ)
                  IJT = NATORB(IQ)+IO-1
                  IPOS = LTPOS(ICT,IJT)
                  IF (ABS(P(IPOS)).LT.CRIT) THEN
                     MOL = ' '
                     WRITE(IW,9056) MOL
                     GOTO 732
                  ENDIF
                  IF (ABS(P(IPOS)).LT.0.3D+00) THEN
                     MOL = ' '
                  ELSE
                     MOL = '#'
                  ENDIF
                  IF (ICT.EQ.IJT) MOL = '*'
                  WRITE(IW,9055) P(IPOS),MOL
  732          CONTINUE
               WRITE(IW,9057)
  713       CONTINUE
  676    CONTINUE
C
C
      WRITE(IW,8062)
      DO 913 IQ = NST, IAT
         IF (IWRK2(IQ).LT.1) GOTO 913
         IF (IWRK2(IQ).EQ.1) THEN
            WRITE(IW,9040)
            GOTO 913
         ENDIF
         DO 932 IO=1,IWRK2(IQ)
            WRITE(IW,9042)
  932    CONTINUE
         WRITE(IW,9044)
  913 CONTINUE
  613 CONTINUE
      WRITE(IW,'(/)')
C
      NST = II
      IF (IAT.LT.NATOT) GOTO 10
C
C
C
 8050 FORMAT(/1X,'       |',$)
 8060 FORMAT(/1X,'       | ',$)
 8061 FORMAT(/1X,' ATOMS | ',$)
 8062 FORMAT(/1X,'-------|',$)
 7062 FORMAT(/1X,'=======|',$)
 9000 FORMAT(A4,I3,$)
 9003 FORMAT('|',$)
 9004 FORMAT(' |',$)
 9005 FORMAT('| ',$)
 9006 FORMAT('      ',$)
 9020 FORMAT(/1X,' OLMO  |',$)
 9030 FORMAT(2X,I2,'   |',$)
 9032 FORMAT(2X,I2,'  ',$)
 9037 FORMAT(1X,A2,I2,' | ',$)
 9039 FORMAT(1X,A2,I1,'  ',$)
 7037 FORMAT(I2,'-',I2,' | ',$)
 7038 FORMAT(I2,'-',I2,'-| ',$)
 7039 FORMAT(1X,'  ?  | ',$)
 9040 FORMAT('-------|',$)
 7040 FORMAT('=======|',$)
 9042 FORMAT('------',$)
 7042 FORMAT('======',$)
 9044 FORMAT('-|',$)
 7044 FORMAT('=|',$)
 9052 FORMAT(/1X,I2,1X,A2,I2,'| ',$)
 9053 FORMAT('     ',A3,$)
 9054 FORMAT(1F5.2,A3,$)
 9055 FORMAT(1F5.2,A1,$)
 9056 FORMAT('     ',A1,$)
 9057 FORMAT('| ',$)
C
      RETURN
      END
C
C*MODULE LOCAL   *DECK ORIANAL
C     --------------------------------------------------------------
      SUBROUTINE ORIANAL(IW,NATOT,NATS,NTOTORB,NATORB,P,IATB,
     *                   NOSI)
C     ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION NATORB(NATOT),P(NTOTORB*(NTOTORB+1)/2)
      DIMENSION IATB(NATS,NTOTORB)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
C
C     PERFORMS ANALYSIS ON DENSITY MATRIX.
C
C    *  IW IS UNIT NAME FOR WRITING OUTPUT
C    *  NATOT = NUMBER OF CATEGORIES OF ORBITALS.
C    *  NATS  = ACTUAL NUMBER OF ATOMS.
C    *  NTOTORB = TOTAL NUMBER OF ORBITALS.
C    *  NATORB(I) IS WHERE THE ORBITALS OF CATEGORY I BEGIN.
C    *  P IS DENSITY MATRIX, STORED IN LOWER TRIANGLE FORM.
C    *  IATB(I,J) CONTAINS INFORMATION ABOUT WHICH ATOMS ARE INVOLVED
C       IN ORBITAL J.
C    *  NOSI IS THE NUMBER OF ORBITALS THAT ARE NOT ATOMIC IN NATURE.
C
      NOAO = NTOTORB - NOSI
C
      IF (NOSI.EQ.0) THEN
         WRITE(IW,9080)
         WRITE(IW,*)
         DO 67 IU = 1,NATOT
            WRITE(IW,9081) IATB(1,NATORB(IU)),
     *    A(IATB(1,NATORB(IU)))
   67    CONTINUE
         WRITE(IW,*)
C
         DO 89 II=1,NATOT
               IST = NATORB(II)
            IEND = NTOTORB
            IF (II.LT.NATOT) IEND = NATORB(II+1)-1
            TPOP = 0.0D+00
            DO 92 JJ=IST,IEND
               IND = LTPOS(JJ,JJ)
               TPOP = TPOP + P(IND)
   92       CONTINUE
            WRITE(IW,9085) TPOP
   89    CONTINUE
         WRITE(IW,*)
      ENDIF
C
      WRITE(IW,*)
      WRITE(IW,*) 'BONDING CHARACTERISTICS OF ORIENTED LOCALIZED MOS'
      WRITE(IW,*)
      WRITE(IW,*) 'COLUMNS OLMO,COMP: POPULATIONS OF THE OLMO ',
     *    ' AND ITS COMPLEMENTARY MO.'
      WRITE(IW,*) 'COLUMN MXBO      : MAXIMUM TOTAL BOND ORDER ',
     *      'COMPATIBLE WITH THE ORBITAL '
      WRITE(IW,*) '                   OCCUPANCIES FOR ANTISYMMETRIC ',
     *     ' WAVEFUNCTIONS.'
      WRITE(IW,*) 'COLUMN ACBO      : ACTUAL TOTAL BOND ORDER OF ',
     *            'THE OLMO.'
      WRITE(IW,*) 'COLUMN %         : COLUMN ACBO AS A PERCENTAGE OF ',
     *    'COLUMN MXBO.'
      WRITE(IW,*)
     * 'COLUMN NOBDS     : NUMBER OF OTHER INTERACTING OLMOS ',
     * 'X/Y IMPLIES '
      WRITE(IW,'(2A/A)') '            ',
     *   '        X PIJ ELEMENTS >= 0.5 AND Y PIJ ELEMENTS >= 0.2',
     *           '                    BUT < 0.5.'
      WRITE(IW,'(2A/2A/2A)')
     *   ' COLUMNS T11,T21  : TRANSFORMATION MATRIX ',
     *      ' ELEMENTS FOR FORMING THE',
     *'                    BONDING(+) AND ANTIBONDING(-) ',
     *      'COMBINATIONS OF THE OLMO ',
     *         '                    AND ITS COMPLIMENTARY',
     * ' ORBITAL (T22 = T11, T12 = -T21).'
      WRITE(IW,'(2A/A)')
     *  ' COLUMNS N+,N-    : OCCUPATION NUMBERS OF THE ',
     *     'BONDING(+) AND ANTIBONDING(-)',
     *   '                    MOS GIVEN BY COLUMNS OF THE MATRIX T.'
      WRITE(IW,*)
      WRITE(IW,*)
      WRITE(IW,9060)
      WRITE(IW,9061)
      WRITE(IW,9062)
      WRITE(IW,9064)
C
      DO 678 II=1,NATOT
         IST  = NATORB(II)
         IEND = NTOTORB
         IF (II.LT.NATOT) IEND = NATORB(II+1)-1
         IL = 0
         DO 689 IO = IST,IEND
            IL = IL + 1
            IF (NATORB(II).LE.NOAO) THEN
            WRITE(IW,9052) IO,A(IATB(1,IO))
            ELSE
            WRITE(IW,9052) IO,'BD',IL
            ENDIF
            PAN = 0.0D+00
            P11 = 0.0D+00
            IND = LTPOS(IO,IO)
            PAA = P(IND)
            MS = 0
            MW = 0
            DO 786 IY = 1,NATOT
               IF (IY.EQ.II) GOTO 786
               IST1 = NATORB(IY)
               IEND1 = NTOTORB
               IF (IY.LT.NATOT) IEND1 = NATORB(IY+1)-1
               DO 777 IT = IST1,IEND1
                  IND = LTPOS(IT,IO)
                  IF (ABS(P(IND)).GT.0.5D+00) THEN
                     MS = MS + 1
                  ELSE IF (ABS(P(IND)).GT.0.2D+00) THEN
                     MW = MW + 1
                  ENDIF
                  PAN = PAN + P(IND)*P(IND)
                  DO 789 IZ = 1,NATOT
                     IF (IZ.EQ.II) GOTO 789
                     IST2 = NATORB(IZ)
                     IEND2 = NTOTORB
                     IF (IZ.LT.NATOT) IEND2 = NATORB(IZ+1)-1
                     DO 877 IQ = IST2,IEND2
                     INDQ = LTPOS(IQ,IO)
                     INDT = LTPOS(IT,IQ)
                     P11 = P11 + P(IND)*P(INDQ)*P(INDT)
  877                CONTINUE
  789             CONTINUE
  777          CONTINUE
  786       CONTINUE
            IF (ABS(PAN).LT.1.0D-08) THEN
               P11 = 0.0D+00
               PAN = 0.0D+00
               DEL = 0.0D+00
               S2G = 0.0D+00
            ELSE
               P11 = P11/PAN
               PAN = SQRT(PAN)
               DEL = (PAA - P11)/2.0D+00
               S2G = PAN/SQRT(DEL*DEL + PAN*PAN)
            ENDIF
            TOT = P11 + PAA
            PBA = TOT/2.0D+00
            BMX =
     *  SQRT(MIN((PBA*PBA - DEL*DEL),((2-PBA)**2.0D+00 - DEL*DEL)))
            GAM = ASIN(S2G)/2.0D+00
            CG = COS(GAM)
            SG = SIN(GAM)
            BON = PBA + SQRT((DEL*DEL + PAN*PAN))
            ANT = PBA - SQRT((DEL*DEL + PAN*PAN))
            P3 = (PAN/BMX)*100
            WRITE(IW,9070) PAA,P11,BMX,PAN,P3,MS,MW,CG,SG,BON,ANT
  689    CONTINUE
         WRITE(IW,9068)
  678 CONTINUE
C
 9052 FORMAT(/I2,1X,A2,' | ',$)
 9060 FORMAT(6X,'|',2X,'POPULATIONS','  |       BOND ORDERS',10X,
     *       '| BONDING/ANTIBONDING MOS',$)
 9061 FORMAT(/' OLMO ','|',15X,'|',28X,'|')
 9062 FORMAT(6X,'|  OLMO   COMP  |  MXBO   ACBO    %    NOBDS | T11',
     *        '   T21    N+    N-')
 9064 FORMAT('======|===============|============================|',
     *       '========================',$)
 9068 FORMAT(/'------|---------------|----------------------------|',
     *       '------------------------',$)
 9070 FORMAT(1X,F4.2,3X,F4.2,1X,' |  ',
     *    F4.2,3X,F4.2,1X,'  ',F4.1,' ',
     *       2X,I1,'/',I1,2X,'| ',F4.2,2X,F4.2,2X,F4.2,2X,F4.2,$)
 9080 FORMAT(/' QUASIATOMIC VALENCE POPULATIONS')
 9081 FORMAT(' AT',I2,'(',A2,')',$)
 9085 FORMAT(2X,F5.3,'  ',$)
C
      RETURN
      END
C
C*MODULE LOCAL   *DECK LTPOS
C     ---------------------------
      INTEGER FUNCTION LTPOS(I,J)
C     ---------------------------
      IMPLICIT INTEGER(A-Z)
      IMA=MAX(I,J)
      IMI=MIN(I,J)
      LTPOS=(IMA*(IMA-1))/2 + IMI
      RETURN
      END
