C 22 FEB 06 - HL  - SEMIIT: FIX THE DIIS DYN. POLARIZABILITY SOLVER
C 17 JAN 06 - MWS - UPDATE ARGS TO TRFMCX
C 14 NOV 05 - LVS - ADD GMRES, BCGSTAB, DIIS SOLVERS FOR DYN. POL. COMP.
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  - POLRHF: ENSURE ALL NODES DO DYN.POLAR. I/O CORRECTLY
C  1 JUN 05 - MWS - POLRHF: RESTORE GOPARR AFTER SOLVCPDYN
C  7 MAR 05 - IA  - ADD COMPUTATION OF THE DYNAMIC POLARIZABILITY
C 13 FEB 05 - MWS - PAD COMMON BLOCK EDCMP
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 22 MAY 02 - MWS - PROVIDE ADDITIONAL TRFMCX ARGUMENTS
C 26 MAR 02 - MWS - POLRHF: RETURN MEMORY WHEN DONE
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - MWS - PAD OUT EDCMP COMMON
C 19 NOV 00 - RMM,HL - POLAP,POLAPM: CALC POLARIZ. W/APPROX TO CPHF
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C  1 MAY 00 - MWS - POLRHF: UPDATE ARGS TO CPWXYO FOR SPHER. HARM.
C 16 FEB 00 - VK  - LPOLX: MCORE SHOULD NOT BE ZERO IF NOUTA.GT.0
C                   LPOLX: ABORT IF MO-S ARE FROZEN AND SCFTYP=RHF
C 21 DEC 99 - MWS - REMOVE C1 DOWNSHIFT FROM TRANSF. BUT NOT AO INTS
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C 27 OCT 98 - MAF - LPOLX: ALLOW FOR USE OF SPHERICAL HARMONICS
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C  2 SEP 97 - MWS - POLRHF,POLROHF: WATCH FOR 32 BIT INTEGER OVERFLOWS
C 19 JUL 97 - PND - LPOLX: RESET VALUE FOR DSKWRK ON EXIT
C 18 DEC 96 - MWS - POLRHF,POLROHF: FIX BUGS DATING FROM FEB 16
C 12 DEC 96 - JHJ - TRANSFORMED DIPOLE INTS WRITTEN TO DAF 252-4
C 13 JUN 96 - MWS - CORRECT CPHF CALLS TO PASS DIAGONAL HESSIAN STORAGE
C  6 MAR 96 - SPW - LAPOL,LOAPOL,LPOLX: CHANGES FOR LOCALIZED CORE
C 14 SEP 95 - SPW - NEW MODULE FOR LOCALIZED ORBITAL POLARIZABILITIES
C
C*MODULE LOCPOL    *DECK LAPOL
      SUBROUTINE LAPOL(U,UL,HF,HFL,TRAN,DLPOL,VCOC,IA,NXYZ,
     *                  NUM2,NOCC,NVIR,MCORE,NLOC,NSTEP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,MOIDON,EDCOMP,DIPDCM,QADDCM,
     *        DEPRNT,ZDO,POLDCM,POLANG,POLAPP,KMIDPT,POLDYN
C
      PARAMETER (MXATM=2000, NMO=500, MXAO=8192)
C
      DIMENSION U(NOCC,NVIR,NXYZ),
     *          UL(NVIR,3),HF(NUM2,3),HFL(NVIR,3),TRAN(NLOC,NLOC),
     *          IA(NUM),CAPOL(9),APOLANG(9),DLPOL(9,NLOC),VCOC(3*NA)
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 /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 /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, FOUR=4.0D+00, TWO=2.0D+00)
      PARAMETER (THREE=3.0D+00, UNIT=0.52917724924D+00)
C
      CALL DAREAD(IDAF,IODA,HF(1,1),NUM2,252,0)
      CALL DAREAD(IDAF,IODA,HF(1,2),NUM2,253,0)
      CALL DAREAD(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C     ---- CALCULATE ALPHA POLARIZABILITY TENSOR FOR THE CORE ----
C
      IF (IFCORE.EQ.1) THEN
      DO 400 K=1,3
         DO 410 L=1,3
         KL=(K-1)*3 + L
         CAPOL(KL)=ZERO
            DO 420 IVIR=1,NVIR
               DO 430 IOCC=1,MCORE
                  IJ=IA(IVIR+NOCC)+IOCC
                     IF(NSTEP.EQ.0) THEN
                     CAPOL(KL)=CAPOL(KL)-(FOUR*U(IOCC,IVIR,K)*HF(IJ,L))
                     ELSE
                     CAPOL(KL)=CAPOL(KL)-(TWO*U(IOCC,IVIR,K)*HF(IJ,L))
                     END IF
  430          CONTINUE
  420       CONTINUE
  410    CONTINUE
  400 CONTINUE
C
C     ----- CONVERT UNITS FROM BOHRS TO ANGTROMS ----
C
      DO 435 I=1,9
         IF(POLANG) THEN
            APOLANG(I)=CAPOL(I)*UNIT**THREE
         ELSE
            APOLANG(I)=CAPOL(I)
         END IF
  435 CONTINUE
C
C     ----- PRINT CORE ALPHA POLARIZABILITY MATRIX -----
C
      IF (MASWRK) THEN
         IF(POLANG) THEN
            WRITE(IW,9020)
         ELSE
            WRITE(IW,9025)
         END IF
         WRITE(IW,9040)
         WRITE(IW,9060) (APOLANG(I),I=1,3)
         WRITE(IW,9080) (APOLANG(I),I=4,6)
         WRITE(IW,9100) (APOLANG(I),I=7,9)
      END IF
C
C     ---- CALCULATE THE MEAN CORE  POLARIZABILITY ----
C
      ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
      IF(MASWRK) WRITE(IW,9110) ALPHAM
C
      ELSE
         CALL VCLR(CAPOL,1,9)
      END IF
C
C     -- CALCULATE ALPHA POLARIZABILITY TENSOR FOR LOCALIZED ORBITALS --
C
      CALL DAREAD(IDAF,IODA,TRAN,NLOC*NLOC,73,0)
C
      DO 440 LOC=1,NLOC
         CALL VCLR(HFL,1,3*NVIR)
         CALL VCLR(UL,1,3*NVIR)
         DO 450 IVIR=1,NVIR
            IL=0
            DO 460 IOCC=MCORE+1,NOCC
               IL=IL+1
               IJ=IA(IVIR+NOCC)+IOCC
               HFL(IVIR,1)=HFL(IVIR,1)+(HF(IJ,1)*TRAN(IL,LOC))
               HFL(IVIR,2)=HFL(IVIR,2)+(HF(IJ,2)*TRAN(IL,LOC))
               HFL(IVIR,3)=HFL(IVIR,3)+(HF(IJ,3)*TRAN(IL,LOC))
               UL(IVIR,1)=UL(IVIR,1)+(U(IOCC,IVIR,1)*TRAN(IL,LOC))
               UL(IVIR,2)=UL(IVIR,2)+(U(IOCC,IVIR,2)*TRAN(IL,LOC))
               UL(IVIR,3)=UL(IVIR,3)+(U(IOCC,IVIR,3)*TRAN(IL,LOC))
  460       CONTINUE
  450    CONTINUE
C
         DO 470 K=1,3
            DO 480 L=1,3
               KL=(K-1)*3 + L
               DLPOL(KL,LOC)=ZERO
               DO 490 IVIR=1,NVIR
                  IF(NSTEP.EQ.0) THEN
                  DUM=(FOUR*UL(IVIR,K)*HFL(IVIR,L))
                  ELSE
                  DUM=(TWO*UL(IVIR,K)*HFL(IVIR,L))
                  END IF
                  DLPOL(KL,LOC)=DLPOL(KL,LOC)-DUM
                  CAPOL(KL)=CAPOL(KL)-DUM
  490          CONTINUE
  480       CONTINUE
  470    CONTINUE
  440 CONTINUE
C
C
C     ---- WRITE ALPHA POLARIZABILITY TENSORS TO DAF ----
C
      IF(NSTEP.EQ.0) CALL DAWRIT(IDAF,IODA,DLPOL,NLOC*9,251,0)
C
      DO 510 I=1,NLOC
         DO 520 J=1,9
            IF(POLANG) THEN
               APOLANG(J)=DLPOL(J,I)*UNIT**THREE
            ELSE
               APOLANG(J)=DLPOL(J,I)
            END IF
  520    CONTINUE
C
         ILO=I+MCORE
         IF (MASWRK) THEN
            IF(NMOAT(ILO).EQ.1) WRITE(IW,9000) MOIDNO(1,ILO)
            IF(NMOAT(ILO).EQ.2) WRITE(IW,9005)
     *                          MOIDNO(1,ILO),MOIDNO(2,ILO)
            WRITE(IW,9050)
            WRITE(IW,9060) (APOLANG(K),K=1,3)
            WRITE(IW,9080) (APOLANG(K),K=4,6)
            WRITE(IW,9100) (APOLANG(K),K=7,9)
            ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
            WRITE(IW,9110) ALPHAM
         END IF
  510 CONTINUE
C
C     ----- PRINT TOTAL ALPHA POLARIZABILITY MATRIX -----
C
      DO 550 I=1,9
         IF(POLANG) THEN
            APOLANG(I)=CAPOL(I)*UNIT**THREE
         ELSE
            APOLANG(I)=CAPOL(I)
         END IF
  550 CONTINUE
C
      IF (MASWRK) THEN
         WRITE(IW,9055)
         WRITE(IW,9060) (APOLANG(I),I=1,3)
         WRITE(IW,9080) (APOLANG(I),I=4,6)
         WRITE(IW,9100) (APOLANG(I),I=7,9)
         ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
         WRITE(IW,9110) ALPHAM
      END IF
C
C     ---- ROTATE BONDS AND LONE PAIRS ONTO Z AXIS ----
C
      IF(NSTEP.NE.0) GO TO 601
      IF(MASWRK) WRITE(IW,9010)
      DO 600 IDLOC=1,NLOC
C
         CALL ROTPOL(DLPOL,VCOC,NLOC,IDLOC,MCORE)
C
         DO 610 J=1,9
            IF(POLANG) THEN
               APOLANG(J)=DLPOL(J,IDLOC)*UNIT**THREE
            ELSE
               APOLANG(J)=DLPOL(J,IDLOC)
            END IF
  610    CONTINUE
C
         IF (MASWRK) THEN
            ILO=IDLOC+MCORE
            IF(NMOAT(ILO).EQ.1) WRITE(IW,9000) MOIDNO(1,ILO)
            IF(NMOAT(ILO).EQ.2) WRITE(IW,9005)
     *                       MOIDNO(1,ILO),MOIDNO(2,ILO)
            WRITE(IW,9050)
            WRITE(IW,9060) (APOLANG(I),I=1,3)
            WRITE(IW,9080) (APOLANG(I),I=4,6)
            WRITE(IW,9100) (APOLANG(I),I=7,9)
         END IF
C
         ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
         IF(MASWRK) WRITE(IW,9110) ALPHAM
  600 CONTINUE
  601 RETURN
C
 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-)/)
 9000 FORMAT(//1X,' LMO ALPHA POLARIZABILITY TENSOR ',
     *        'FOR CORE OR LONE PAIR ON ATOM',I3)
 9005 FORMAT(//1X,' LMO ALPHA POLARIZABILITY TENSOR ',
     *        'FOR BOND BETWEEN ATOM',I3,' AND ATOM',I3)
 9010 FORMAT(//10X,47(1H-)/
     *        10X,' LOCALIZED ALPHA POLARIZABILITY TENSORS FOR '/
     *        10X,' LONE PAIRS AND BONDS PROJECTED ONTO THE Z-AXIS '/
     *        10X,47(1H-))
 9040 FORMAT(10X,' CORE ALPHA POLARIZABILITY TENSOR '/
     *       22X,'UX',13X,'UY',13X,'UZ')
 9050 FORMAT(22X,'UX',13X,'UY',13X,'UZ')
 9055 FORMAT(//15X,35(1H-)/
     *        15X,' TOTAL ALPHA POLARIZABILITY TENSOR '/
     *        15X,35(1H-)/
     *        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*MODULE LOCPOL    *DECK LOAPOL
      SUBROUTINE LOAPOL(U,UL,HF,HFL,TRAN,DLPOL,VCOC,IA,NXYZ,
     *                  NUM2,MCORE,NLOC,LOCS,NROT,IROT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,MOIDON,EDCOMP,DIPDCM,QADDCM,
     *        DEPRNT,ZDO,POLDCM,POLANG,POLAPP,KMIDPT,POLDYN
C
      PARAMETER (MXATM=2000, NMO=500, MXAO=8192)
C
      DIMENSION U(NROT,NXYZ),
     *          UL(LOCS,3),HF(NUM2,3),HFL(LOCS,3),TRAN(NLOC,NLOC),
     *          IA(NUM),CAPOL(9),APOLANG(9),DLPOL(9,NLOC),VCOC(3*NA),
     *          IROT(NUM,NUM)
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 /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /GSSORD/ IORDER(MXAO),JORDER(MXAO),NORDER
      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 /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, FOUR=4.0D+00)
      PARAMETER (THREE=3.0D+00, UNIT=0.52917724924D+00)
C
      L1=NUM
C
C     ----- SET UP ROTATION INDEX ARRAY -----
C
      MROT=0
      DO 220 I=1,L1
         ISH = NCONF(I)
         DO 210 J=1,I
            IF(NCONF(J).EQ.ISH) THEN
               IROT(I,J) = 0
               IROT(J,I) = 0
            ELSE
               MROT=MROT+1
               IROT(I,J) = MROT
               IROT(J,I) = MROT
            END IF
  210    CONTINUE
  220 CONTINUE
C
      CALL DAREAD(IDAF,IODA,HF(1,1),NUM2,252,0)
      CALL DAREAD(IDAF,IODA,HF(1,2),NUM2,253,0)
      CALL DAREAD(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C     ---- CALCULATE ALPHA POLARIZABILITY TENSOR FOR THE CORE ----
C
      IF (IFCORE.EQ.1) THEN
      JMAX=MCORE
      DO 400 K=1,3
         DO 410 L=1,3
         KL=(K-1)*3 + L
         CAPOL(KL)=ZERO
            DO 420 I=1,NUM
               ISH = NCONF(I)
               DO 430 J=1,I
                  IF (J .GT. JMAX) GO TO 420
                  JSH = NCONF(J)
                  IF(ISH.EQ.JSH) GO TO 430
                  FIJ  = F(ISH) - F(JSH)
                  IJ = IA(I) + J
                  IJROT = IROT(I,J)
                  CAPOL(KL)=CAPOL(KL)+(FOUR*FIJ*U(IJROT,K)*HF(IJ,L))
  430          CONTINUE
  420       CONTINUE
  410    CONTINUE
  400 CONTINUE
C
C     ----- CONVERT UNITS FROM BOHRS TO ANGTROMS ----
C
      DO 435 I=1,9
         IF(POLANG) THEN
            APOLANG(I)=CAPOL(I)*UNIT**THREE
         ELSE
            APOLANG(I)=CAPOL(I)
         END IF
  435 CONTINUE
C
C     ----- PRINT CORE ALPHA POLARIZABILITY MATRIX -----
C
      IF (MASWRK) THEN
         IF(POLANG) THEN
            WRITE(IW,9020)
         ELSE
            WRITE(IW,9025)
         END IF
         WRITE(IW,9040)
         WRITE(IW,9060) (APOLANG(I),I=1,3)
         WRITE(IW,9080) (APOLANG(I),I=4,6)
         WRITE(IW,9100) (APOLANG(I),I=7,9)
      END IF
C
C     ---- CALCULATE THE MEAN CORE  POLARIZABILITY ----
C
      ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
      IF(MASWRK) WRITE(IW,9110) ALPHAM
C
      ELSE
         CALL VCLR(CAPOL,1,9)
      END IF
C
C     -- CALCULATE ALPHA POLARIZABILITY TENSOR FOR LOCALIZED ORBITALS --
C
      CALL DAREAD(IDAF,IODA,TRAN,NLOC*NLOC,73,0)
C
      CALL VCLR(DLPOL,1,9*NLOC)
      DO 440 LOC=1,NLOC
         CALL VCLR(HFL,1,3*LOCS)
         CALL VCLR(UL,1,3*LOCS)
         DO 450 I=1,LOCS
            IL=0
            JSH = NCONF(MCORE+I)
            DO 460 IOCC=MCORE+1,NA
               ISH = NCONF(IOCC)
               IL=IL+1
               IF(IOCC.GT.I+MCORE) GO TO 460
               IJ=IA(I+MCORE)+IOCC
               IJROT = IROT(I+MCORE,IOCC)
               IF (IJROT.EQ.0) GO TO 460
               HFL(I,1)=HFL(I,1)+(HF(IJ,1)*TRAN(IL,LOC))
               HFL(I,2)=HFL(I,2)+(HF(IJ,2)*TRAN(IL,LOC))
               HFL(I,3)=HFL(I,3)+(HF(IJ,3)*TRAN(IL,LOC))
               UL(I,1)=UL(I,1)+(U(IJROT,1)*TRAN(IL,LOC))
               UL(I,2)=UL(I,2)+(U(IJROT,2)*TRAN(IL,LOC))
               UL(I,3)=UL(I,3)+(U(IJROT,3)*TRAN(IL,LOC))
  460       CONTINUE
  450    CONTINUE
C
         JSH = NCONF(MCORE+LOC)
         DO 470 K=1,3
            DO 480 L=1,3
               KL=(K-1)*3 + L
               DO 490 I=1,LOCS
                  ISH = NCONF(I+MCORE)
                  FIJ  = F(ISH) - F(JSH)
                  IF(LOC.GT.I) GO TO 490
                  IF(ISH.EQ.JSH) GO TO 490
                  DUM=(FOUR*FIJ*UL(I,K)*HFL(I,L))
                  DLPOL(KL,LOC)=DLPOL(KL,LOC)+DUM
                  CAPOL(KL)=CAPOL(KL)+DUM
  490          CONTINUE
  480       CONTINUE
  470    CONTINUE
  440 CONTINUE
C
C     ---- WRITE ALPHA POLARIZABILITY TENSORS TO DAF ----
      CALL DAWRIT(IDAF,IODA,DLPOL,NLOC*9,251,0)
C
C     ---- REORDER ----
      CALL MOIDM
      CALL DIPLMM
      CALL ICOPY(L1,IORDER,1,JORDER,1)
      ICNT = 1
      DO 505 I = 1,NA
         J = IORDER(I) - MCORE
         IF (J .GT. 0) THEN
            JORDER(ICNT) = J
            ICNT = ICNT + 1
         END IF
  505 CONTINUE
      CALL REORDR(DLPOL,JORDER,NLOC,9)
C
      DO 510 I=1,NLOC
         DO 520 J=1,9
            IF(POLANG) THEN
               APOLANG(J)=DLPOL(J,I)*UNIT**THREE
            ELSE
               APOLANG(J)=DLPOL(J,I)
            END IF
  520    CONTINUE
C
         ILO=I+MCORE
         IF (MASWRK) THEN
            IF(NMOAT(ILO).EQ.1) WRITE(IW,9000) MOIDNO(1,ILO)
            IF(NMOAT(ILO).EQ.2) WRITE(IW,9005)
     *                          MOIDNO(1,ILO),MOIDNO(2,ILO)
            WRITE(IW,9050)
            WRITE(IW,9060) (APOLANG(K),K=1,3)
            WRITE(IW,9080) (APOLANG(K),K=4,6)
            WRITE(IW,9100) (APOLANG(K),K=7,9)
            ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
            WRITE(IW,9110) ALPHAM
         END IF
  510 CONTINUE
C
C     ----- PRINT TOTAL ALPHA POLARIZABILITY MATRIX -----
C
      DO 550 I=1,9
         IF(POLANG) THEN
            APOLANG(I)=CAPOL(I)*UNIT**THREE
         ELSE
            APOLANG(I)=CAPOL(I)
         END IF
  550 CONTINUE
C
      IF (MASWRK) THEN
         WRITE(IW,9055)
         WRITE(IW,9060) (APOLANG(I),I=1,3)
         WRITE(IW,9080) (APOLANG(I),I=4,6)
         WRITE(IW,9100) (APOLANG(I),I=7,9)
         ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
         WRITE(IW,9110) ALPHAM
      END IF
C
C     ---- ROTATE BONDS AND LONE PAIRS ONTO Z AXIS ----
C
      IF(MASWRK) WRITE(IW,9010)
      DO 600 IDLOC=1,NLOC
         CALL ROTPOL(DLPOL,VCOC,NLOC,IDLOC,MCORE)
         DO 610 J=1,9
            IF(POLANG) THEN
               APOLANG(J)=DLPOL(J,IDLOC)*UNIT**THREE
            ELSE
               APOLANG(J)=DLPOL(J,IDLOC)
            END IF
  610    CONTINUE
C
         IF (MASWRK) THEN
            ILO=IDLOC+MCORE
            IF(NMOAT(ILO).EQ.1) WRITE(IW,9000) MOIDNO(1,ILO)
            IF(NMOAT(ILO).EQ.2) WRITE(IW,9005)
     *                          MOIDNO(1,ILO),MOIDNO(2,ILO)
            WRITE(IW,9050)
            WRITE(IW,9060) (APOLANG(I),I=1,3)
            WRITE(IW,9080) (APOLANG(I),I=4,6)
            WRITE(IW,9100) (APOLANG(I),I=7,9)
         END IF
         ALPHAM=(APOLANG(1)+APOLANG(5)+APOLANG(9))/THREE
         IF(MASWRK) WRITE(IW,9110) ALPHAM
  600 CONTINUE
C
      RETURN
C
 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-)/)
 9000 FORMAT(//14X,' LMO ALPHA POLARIZABILITY TENSOR '/
     *        10X,'FOR CORE/LONE PAIR/UNPAIRED ELECTRON ON ATOM',I3)
 9005 FORMAT(//1X,' LMO ALPHA POLARIZABILITY TENSOR ',
     *        'FOR BOND BETWEEN ATOM',I3,' AND ATOM',I3)
 9010 FORMAT(//10X,47(1H-)/
     *        10X,' LOCALIZED ALPHA POLARIZABILITY TENSORS FOR '/
     *        10X,' LONE PAIRS AND BONDS PROJECTED ONTO THE Z-AXIS '/
     *        10X,47(1H-))
 9040 FORMAT(10X,' CORE ALPHA POLARIZABILITY TENSOR '/
     *       22X,'UX',13X,'UY',13X,'UZ')
 9050 FORMAT(22X,'UX',13X,'UY',13X,'UZ')
 9055 FORMAT(//15X,35(1H-)/
     *        15X,' TOTAL ALPHA POLARIZABILITY TENSOR '/
     *        15X,35(1H-)/
     *        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*MODULE LOCPOL    *DECK LPOLX
      SUBROUTINE LPOLX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL PPGVB,TCSCF,ROGVB,SOME,DBG,GOPARR,DSKWRK,MASWRK,
     *        MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,POLDCM,POLANG,
     *        POLAPP,KMIDPT,POLDYN,
     *        FIRST,SECND,CPHF,BOTH,MFIRST,MSECND,MCPHF,
     *        PK,PANDK,BLOCK,DIRSCF,DIRTRF,FDIFF,SVGPAR,SVDSKW,
     *        ABEL,ABELPT,MFRZ
C
      PARAMETER (MXATM=2000, MXSH=5000, NMO=500, MXAO=8192)
C
      DIMENSION NOCP(3)
C
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /HSSPAR/ FIRST,SECND,CPHF,BOTH,MFIRST,MSECND,MCPHF
      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,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      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 CPWORD/8HCPHF    /
      DATA RHF/8HRHF     /, ROHF/8HROHF    /, GVB/8HGVB     /,
     *     UHF/8HUHF     /, RMC/8HMCSCF   /, GUGA/8HGUGA    /
C
C     ----- MAIN DRIVER FOR LOCALIZED ORBITAL POLARIZABILITIES -----
C     CODE WRITTEN BY SIMON WEBB AT ISU, SEPTEBER 1995.
C     THE THREE ELECTRIC FIELD RESPONSES ARE BASICALLY WHAT IS NEEDED.
C
C               ADDITION OF NEW DYNAMIC POLARIZABILITY
C           ON THE END OF STATIC POLARIZABILITY CALCULATIONS
C
      IF(POLDYN) POLAPP = .FALSE.
      IF(POLAPP) GO TO 100
      CALL CPINP
C
      FIRST=.TRUE.
      SECND=.TRUE.
      CPHF=.TRUE.
      BOTH=.TRUE.
      MFIRST=.TRUE.
      MSECND=.TRUE.
      MCPHF=.TRUE.
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,*) 'LPOLX: SOMEHOW THIS IS A -PK- FILE'
         CALL ABRT
      END IF
C
C     IF DUPLICATED INTEGRAL FILES DESIRED FOR TRANSF., GENERATE THEM.
C     IF POINT GROUP IS NON-ABELIAN, GENERATE A C1 INTEGRAL FILE.
C
      ABEL = ABELPT()
      NTSAVE = NT
      SVDSKW = DSKWRK
      SVGPAR = GOPARR
      IF(GOPARR) DSKWRK = .TRUE.
      IF(CPHF  .AND.
     *   ((GOPARR.AND.ITRFAO.EQ.1)
     *    .OR.  .NOT.ABEL)) 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     ----- INTEGRAL TRANSFORMATION -----
C
      NOCC = 0
      IF(SCFTYP.EQ.RHF)  NOCC=NA
      IF(SCFTYP.EQ.UHF)  NOCC=NA
      IF(SCFTYP.EQ.ROHF) NOCC=NA
      IF(SCFTYP.EQ.GVB)  NOCC=NORB
      IF(SCFTYP.EQ.RMC)  NOCC=NMCC+NDOC+NAOS+NBOS+NALP+NVAL
      IF(DIRSCF) DIRTRF=.TRUE.
      IF(CPHF) CALL TRFMCX(0,0,NOCC,NQMT,.FALSE.,.FALSE.,
     *                     .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,
     *                     .FALSE.,.FALSE.,.TRUE.)
      IF(DIRSCF) DIRTRF=.FALSE.
C
C     ----- MAIN DRIVER FOR CPHF CONTRIBUTION TO HESSIAN -----
C
  100 CONTINUE
      DBG=EXETYP.EQ.CPWORD .AND. MASWRK
      SOME = NPRINT.NE.-5 .AND. MASWRK
      NFT16  = 16
      NFTH2  = 67
      NFTH21 = 68
C
C     ----- DEFINE UNIQUE -XYZ- FOR WHICH -CPHF- IS NEEDED -----
C
      NXYZ = 3
      NUNIQ= 3
      NOCP(1)=0
      NOCP(2)=0
      NOCP(3)=0
C
      CALL LMOPSI(0,MCORE,MDOC,MACT,NUMLOC)
      IF (IFCORE.EQ.0) MCORE=0
C
C     THE POLARIZATION CODE DID NOT CORRECTLY TAKE INTO ACCOUNT CASE
C     OF NOUTA BEING NON-ZERO! FIXING IT...
C
      IF(NOUTA.GT.0 .AND. IFCORE.EQ.0) THEN
        IFCORE=1
        MCORE=NOUTA
        IF(MASWRK) THEN
           WRITE(IW,*)
           WRITE(IW,*)' ADJUSTING MCORE IN LPOLX'
         END IF
        IF(MOOUTA(NOUTA).NE.NOUTA) THEN
          IF(MASWRK) WRITE(IW,9100)
          CALL ABRT
        END IF
      END IF
C
      NLOC = NUMLOC-MCORE
C
C     ----- DO THE APPROPRIATE CPHF PROBLEM -----
C
      ROGVB = (SCFTYP.EQ.GVB  .AND.  NPAIR.EQ.0)  .OR.  SCFTYP.EQ.ROHF
      TCSCF = (SCFTYP.EQ.GVB  .AND.  NPAIR.EQ.1)
      PPGVB = (SCFTYP.EQ.GVB  .AND.  NPAIR.GT.1)
C
      IF(SCFTYP.EQ.RHF .AND. MFRZ .AND. MASWRK .AND. .NOT.POLAPP) THEN
        WRITE(IW,*)
        WRITE(IW,*)'            ***************************************'
        WRITE(IW,*)'            *   *ERROR*  THE POLARIZABILITIES     *'
        WRITE(IW,*)'            *   WILL BE INCORRECT BECAUSE OF THE  *'
        WRITE(IW,*)'            *   NON-DIAGONAL FOCK MATRIX DUE TO   *'
        WRITE(IW,*)'            *   FREEZING MO-S DURING HARTREE-FOCK.*'
        WRITE(IW,*)'            *   FIX: REPEAT RUN WITH RUNTYP=ROHF. *'
        WRITE(IW,*)'            ***************************************'
        WRITE(IW,*)
        CALL ABRT
      END IF
C
      IF(SCFTYP.EQ.RHF.AND..NOT.POLAPP) THEN
                        CALL POLRHF(SOME,DBG,NOCP,NXYZ,NUNIQ,
     *                       NFT16,MCORE,NLOC,NFTH2,NFTH21)
      END IF
      IF(POLAPP)        CALL POLAPM(SOME,NXYZ,MCORE,NLOC)
      IF(SCFTYP.EQ.UHF) CALL ABRT
      IF(ROGVB)         CALL POLROHF(TCSCF,SOME,DBG,NOCP,NXYZ,NUNIQ,
     *                               NFT16,MCORE,NLOC)
      IF(TCSCF)         CALL ABRT
      IF(PPGVB)         CALL ABRT
      IF(SCFTYP.EQ.RMC) CALL ABRT
      IF(CITYP.EQ.GUGA) CALL ABRT
C
      DSKWRK = SVDSKW
      IF (MASWRK) WRITE(IW,9000)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/1X,'...... DONE WITH CPHF CONTRIBUTIONS ......')
 9100 FORMAT(//1X,'ERROR IN LPOLX:'/
     *   1X,'YOUR MOOUTA ARRAY SHOULD CONSIST OF CONSECUTIVE NUMBERS'/
     *   1X,'STARTING WITH 1 (I.E. IMITATE CORE),'/
     *   1X,'OTHERWISE THE POLARIZABILITIES WILL BE INCORRECT.'/
     *   1X,'IF YOU REALLY HAVE TO FREEZE OUT OF LOCALIZATION',
     *   1X,'SOME ARBITRARY ORBITALS,'/
     *   1X,'USE OF $MOFRZ MIGHT HELP.')
      END
C*MODULE LOCPOL    *DECK OPOLB
      SUBROUTINE OPOLB(HF,THF,V,WRK,L1,L2,L3,NXYZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION HF(L2),THF(L2,NXYZ),V(L3),WRK(L1)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
C     READ MO VECTORS.
C
      CALL DAREAD(IDAF,IODA,V,L3,15,0)
C
C     ---- GET DIPOLE INTEGRALS AND TRANSFORM TO MO BASIS   ----
C     ---- WRITE TRANSFORMED INTEGRALS TO DAF FOR LATER USE ----
C
      CALL DAREAD(IDAF,IODA,HF,L2,95,0)
      CALL TFTRI(THF(1,1),HF,V,WRK,L1,L1,L1)
      CALL DAREAD(IDAF,IODA,HF,L2,96,0)
      CALL TFTRI(THF(1,2),HF,V,WRK,L1,L1,L1)
      CALL DAREAD(IDAF,IODA,HF,L2,97,0)
      CALL TFTRI(THF(1,3),HF,V,WRK,L1,L1,L1)
C
      CALL DAWRIT(IDAF,IODA,THF(1,1),L2,252,0)
      CALL DAWRIT(IDAF,IODA,THF(1,2),L2,253,0)
      CALL DAWRIT(IDAF,IODA,THF(1,3),L2,254,0)
C
      RETURN
      END
C*MODULE LOCPOL    *DECK POLB
      SUBROUTINE POLB(B,C,T,D,HF,IA,NOCC,
     *                    NVIR,NORB,NXYZ,NUM,NUM2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION B(NOCC,NVIR,NXYZ),C(NUM,NUM),T(NUM),D(NUM2),
     *          HF(NUM2,3),IA(*)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
      PARAMETER (FOUR=4.0D+00)
C
C     ---- GET DIPOLE INTEGRALS AND TRANSFORM TO MO BASIS   ----
C     ---- WRITE TRANSFORMED INTEGRALS TO DAF FOR LATER USE ----
C     ---- THEN FORM THE MATRIX -B-                         ----
C
      CALL DAREAD(IDAF,IODA,D,NUM2,95,0)
      CALL TFTRI(HF(1,1),D,C,T,NORB,NUM,NUM)
      CALL DAREAD(IDAF,IODA,D,NUM2,96,0)
      CALL TFTRI(HF(1,2),D,C,T,NORB,NUM,NUM)
      CALL DAREAD(IDAF,IODA,D,NUM2,97,0)
      CALL TFTRI(HF(1,3),D,C,T,NORB,NUM,NUM)
C
      CALL DAWRIT(IDAF,IODA,HF(1,1),NUM2,252,0)
      CALL DAWRIT(IDAF,IODA,HF(1,2),NUM2,253,0)
      CALL DAWRIT(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C     ---- NOTE THAT MULTIPLICATION BY FOUR NEEDED TO RECONCILE ----
C     ---- SCHAEFERS AND KORMORNIKIS FORMULATIONS               ----
C
         DO 650 IVIR=1,NVIR
            DO 660 IOCC=1,NOCC
               IJ=IA(IVIR+NOCC)+IOCC
               B(IOCC,IVIR,1)=FOUR*HF(IJ,1)
               B(IOCC,IVIR,2)=FOUR*HF(IJ,2)
               B(IOCC,IVIR,3)=FOUR*HF(IJ,3)
  660       CONTINUE
  650    CONTINUE
C
      RETURN
      END
C*MODULE LOCPOL    *DECK POLRHF
      SUBROUTINE POLRHF(SOME,DBG,NOCP,NXYZ,NUNIQ,NFT16,MCORE,NLOC,
     *                  NFTH2,NFTH21)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG,SOME,CHKMEM,PACK2E,INCORE,GOPARR,DSKWRK,MASWRK
      LOGICAL SVDSKW,SVGPAR,POLAR
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,
     *        POLDCM,POLANG,POLAPP,KMIDPT,POLDYN
C
      PARAMETER (MXATM=2000, MXAO=8192, NMO=500)
C
      DIMENSION NOCP(3)
C
      COMMON /CHFINF/ CHFSLV,CPTYPE,NWDCHF,POLAR
      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 /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C     ----- CLOSED SHELL CPHF CALCULATION -----
C
      IF(SOME) WRITE(IW,9000)
      NXYZ = 3
      NORB = NQMT
      NOCC = NA
      NVIR = NORB - NOCC
      NROT = NOCC * NVIR
      IF(SOME) WRITE(IW,9010) NROT
C
      NUM2 =(NUM*NUM+NUM)/2
      NUM3 = NUM*NUM
C
      IF(NVIR .EQ. 0) RETURN
C
C     ----- BEGIN CPHF TERMS -----
C     FIRST BRANCH ALL ABOUT TO CHECK MEMORY NEEDS,
C     THEN COME BACK TO -90- AND ACTUALLY BEGIN TO DO THE WORK.
C
      TIM1 = 0.0D+00
      CHKMEM = .TRUE.
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      IF(NWDCHF.GT.0) NGOTMX=MIN(NGOTMX,NWDCHF)
      INCORE=.FALSE.
      NDB    = 0
      NDA    = 0
      NDAI = 0
      GO TO 125
C
   90 CONTINUE
      CHKMEM = .FALSE.
      IF(SOME) THEN
         WRITE(IW,9020) NDB,NDA,NGOTMX
         IF(.NOT.INCORE) WRITE(IW,9030) NDAI
      END IF
      NDMAX = MAX(NDB,NDA)
      IF(NDMAX.GT.NGOTMX) THEN
         IF (MASWRK) WRITE(IW,9040) NDMAX,NGOTMX
         CALL ABRT
         STOP
      END IF
C
      IF(EXETYP.EQ.CHECK) GO TO 900
C
C     ----- FORM -B- MATRIX -----
C
      CALL TSECND(TIM1)
  125 CONTINUE
      LB   = 1+LOADFM
      LC   = LB   + NROT *NXYZ
      LT   = LC   + NUM3
      LD   = LT   + NUM
      LHF  = LD   + NUM2
      LAST = LHF  + NUM2*3
C
      NEED = LAST - LOADFM - 1
      NDB   = NEED
      IF(CHKMEM) GO TO 135
      CALL GETFM(NEED)
C
      CALL DAREAD(IDAF,IODA,X(LC),NUM3,15,0)
C
      CALL POLB(X(LB),X(LC),X(LT),X(LD),X(LHF),IA,NOCC,NVIR,NORB,NXYZ,
     *          NUM,NUM2)
C
      IF(DBG) THEN
         WRITE(IW,*) '-B- MATRIX'
         CALL PRSQ(X(LB),NXYZ,NROT,NROT)
      END IF
C
      CALL TSECND(TIM2)
      TIMFCK=TIM2-TIM1
      IF(SOME) WRITE(IW,9988) ' -B-  ',TIMFCK
C
      CALL RETFM(NEED)
C
C     ----- CALCULATE -A- MATRIX AND SOLVE CPHF EQNS -----
C     NOTE THAT -B- AND -U- MUST BE THE FIRST ARRAYS ALLOCATED,
C     AS ASSUMED BY THE CODE BELOW THE SOLUTION STEP.
C
      CALL TSECND(TIM1)
  135 CONTINUE
      LB   = 1+LOADFM
      LAST = LB   + NROT*NXYZ
C                         OVERLAP PART OF -A- FORMATION STORAGE...
      LE   = LAST
      LXX  = LE   + NUM
      LIX  = LXX  + NINTMX
      LAST1= LIX  + NINTMX
C                         ...WITH PART OF THE CPHF SOLUTION STORAGE...
      LU   = LAST
      LAST2= LU   + NROT*NXYZ
C
      LAST = MAX(LAST1,LAST2)
      NNEED= LAST-LB
      NFREE= NGOTMX-NNEED
C
C         DECIDE IF PROBLEM FITS IN AVAILABLE MEMORY,
C         WATCH OUT FOR 32 BIT SYSTEMS
C
      IF(NWDVAR.EQ.1  .OR.  NROT.LE.46300) THEN
         INMEM  = NROT*NROT + 2 + 5*NROT
         INCORE = INMEM.LE.NFREE
         NDAI = NNEED + INMEM
      ELSE
         INMEM = 0
         INCORE = .FALSE.
         NDAI = 0
      END IF
C
      IF(POLDYN) THEN
         INCORE=.FALSE.
         INMEM = 0
         NDAI = 0
      END IF
C                         ...FOLLOWED BY MORE SOLVER STORAGE,
C                         WITH -A- ITSELF ALLOCATED LAST.
C
C        IF -A- FITS IN MEMORY, WE SOLVE FOR ONE RESPONSE
C        AT A TIME TO MINIMIZE WORKING STORAGE -LWRK-.
C        IF NOT, WE MINIMIZE DISK I/O ON THE EXTERNAL FILE OF
C        -A- BY WORKING ON ALL RESPONSES SIMULTANEOUSLY.
C
C        BY USING THE NEXT FOUR LINES, SMALL PROBLEMS CAN BE FORCED
C        INTO AN OUT OF CORE, 2 SLICE STORAGE OF THE -A- MATRIX:
C
C        IF(INCORE) THEN
C           INCORE=.FALSE.
C           NFREE = NUM*NUM +2*NXYZ +4*NXYZ*NROT +NROT +(NROT*NROT+10)/2
C        END IF
C
      MAXEQU=128
C
C IMPORTANT FOR STATIC POLARIZABILITY
C
      NSOLV = NXYZ
      IF(NROT.LE.MAXEQU .AND. INCORE) NSOLV=1
C
      IF(INCORE) THEN
         NWIDTH = NROT
         IF(SOME  .AND.  CHKMEM) WRITE(IW,9050)
         LIROT= LAST
         LBNRM= LAST
         LRESD= LBNRM + NSOLV
         LWRK1= LRESD + NSOLV
         LWRK2= LWRK1 + NROT*NSOLV
         LWRK3= LWRK2 + NROT*NSOLV
         LWRK4= LWRK3 + NROT*NSOLV
         LWII = LWRK4 + NROT*NSOLV
         LA   = LWII  + NROT
         LAST = LA    + NROT*NROT
      ELSE
         NWIDTH = (NFREE - NUM*NUM - 2*NXYZ- 4*NROT*NXYZ - NROT)/NROT
C
C TAKE IN ACCOUNT FACT THAT WE NEED 3 TIMES MORE MEMORY
C
         IF(POLDYN) NWIDTH =  NWIDTH/3
         IF(NWIDTH.GT.NROT) NWIDTH = NROT
C
C   DEBUG
C        INCORE = .FALSE.
C        NWIDTH = 1000
C
         IF(NWIDTH.GT.0) THEN
            NSLICE = (NROT-1)/NWIDTH + 1
         ELSE
            NSLICE=-1
         END IF
         IF(NWIDTH.LE.0  .OR.  (NSLICE.GT.NAT .AND. NAT.GT.5)) THEN
            IF (MASWRK) THEN
               NEEDX = NUM*NUM + 2*NXYZ + 4*NROT*NXYZ + NROT
               NNEED = NNEED + NEEDX
               NFREE = NFREE - NEEDX
               NEEDW = (NROT-1)/(NAT-1)
               NEED  = NNEED + NROT*NEEDW
               WRITE(IW,9070) NWIDTH,NSLICE,NAT,NFREE,NNEED,NROT,NEED
            END IF
            CALL ABRT
            STOP
         END IF
         IF(SOME  .AND.  CHKMEM) WRITE(IW,9060) NSLICE,NWIDTH
         LIROT= LAST
         LBNRM= LIROT + NUM*NUM
         LRESD= LBNRM + NSOLV
         LWRK1= LRESD + NSOLV
         LWRK2= LWRK1 + NROT*NSOLV
         LWRK3= LWRK2 + NROT*NSOLV
         LWRK4= LWRK3 + NROT*NSOLV
         LWII = LWRK4 + NROT*NSOLV
         LA   = LWII  + NROT
         LAST = LA    + NROT*NWIDTH
      END IF
C
      NEED = LAST - LOADFM - 1
      NDA   = NEED
      IF(CHKMEM) GO TO 90
      CALL GETFM(NEED)
C
C        GET RHF ORBITAL ENERGIES
C
      CALL DAREAD(IDAF,IODA,X(LE),NUM,17,0)
C
      IF(INCORE) THEN
         NWIDTH = NROT
         CALL CPWXYI(X(LA),X(LE),X(LXX),X(LIX),
     *               NUM,NOCC,NVIR,IJKT,NINTMX)
         IF(DBG) THEN
            WRITE(IW,9998)
            CALL PRSQ(X(LA),NROT,NROT,NROT)
         END IF
      ELSE
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQOPN(NFT16,'WORK16','UNKNOWN',.FALSE.,'UNFORMATTED')
         DSKWRK = SVDSKW
         CALL CPWXYO(X(LA),X(LE),X(LXX),X(LIX),X(LIROT),X(LWII),
     *               NROT,NWIDTH,NUM,NOCC,NVIR,IJKT,NINTMX,NFT16,DBG)
      END IF
C
      CALL TSECND(TIM2)
      TIMA=TIM2-TIM1
      IF(SOME) THEN
         WRITE(IW,9988) ' -A-  ',TIMA
         CALL FLSHBF(IW)
      END IF
      TIM1 = TIM2
C
C        SOLVE FOR SYMMETRY UNIQUE CPHF RESPONSE VECTORS
C          THIS IS DONE SEQUENTIALLY ON THE MASTER NODE AT PRESENT
C
      NITER=0
      NPASS=0
      SVGPAR = GOPARR
      SVDSKW = DSKWRK
      GOPARR = .FALSE.
      DSKWRK = .FALSE.
      IF(MASWRK) THEN
      CALL SOLVCP(X(LA),X(LB),X(LU),NROT,NXYZ,NOCP,X(LWII),
     *            X(LWRK1),X(LWRK2),X(LWRK3),X(LWRK4),X(LBNRM),X(LRESD),
     *            NITER,NPASS,NWIDTH,NSOLV,INCORE,NFT16,MAXEQU)
      END IF
      GOPARR = SVGPAR
      DSKWRK = SVDSKW
C
      IF(DBG) THEN
         WRITE(IW,*) 'CPHF RESPONSE VECTORS -U-'
         CALL PRSQ(X(LU),NXYZ,NROT,NROT)
      END IF
      CALL RETFM(NEED)
      CALL TSECND(TIM2)
      TIMYA=TIM2-TIM1
      IF(SOME) WRITE(IW,9988) ' -U-   ',TIMYA
      IF(SOME) THEN
         IF(NITER.GT.0) WRITE(IW,9080) NUNIQ,NITER
         IF(.NOT.INCORE) WRITE(IW,9090) NPASS
      END IF
C
C     ---- MEMORY ALLOCATION FOR CALCULATION OF POLARIZABILITY ----
C
      LB   = 1+LOADFM
      LU   = LB   + NROT*NXYZ
      LH   = LU   + NROT*NXYZ
      LUL  = LH   + NUM2*3
      LHL  = LUL  + NVIR*3
      LTR  = LHL  + NVIR*3
      LPOL = LTR  + NLOC*NLOC
      LVCOC= LPOL + NLOC*9
      LAST = LVCOC+ NA*3
      NEED = LAST - LB
      CALL GETFM(NEED)
C
C     ----- CALCULATE LOCALIZED ALPHA POLARIZABILITIES ----
C
      NSTEP=0
      CALL LAPOL(X(LU),X(LUL),X(LH),X(LHL),X(LTR),X(LPOL),X(LVCOC),
     *           IA,NXYZ,NUM2,NOCC,NVIR,MCORE,NLOC,NSTEP)
      CALL RETFM(NEED)
      IF(POLDYN) THEN
         IF(MASWRK) WRITE(IW,*)
     *       '.... DONE WITH STATIC POLARIZABILITY ....'
         CALL TIMIT(1)
      END IF
C
C       START OF IMAGINARY FREQUENCY DEPENDENT (DYNAMIC) POLARIZABILITY
C
      IF(POLDYN) THEN
         LH2   = 1+LOADFM
         LE    = LH2 + NROT*NWIDTH
         LXX   = LE   + NUM
         LIX   = LXX  + NINTMX
         LIROT = LIX  + NINTMX
         LH2II = LIROT + NUM*NUM
         LAST  = LH2II  + NROT
         NEED  = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
C        GET RHF ORBITAL ENERGIES
C
         CALL DAREAD(IDAF,IODA,X(LE),NUM,17,0)
C
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQOPN(NFTH2,'DCPHFH2','UNKNOWN',.FALSE.,'UNFORMATTED')
         DSKWRK = SVDSKW
         CALL CPH2O(X(LH2),X(LE),X(LXX),X(LIX),X(LIROT),X(LH2II),
     *               NROT,NWIDTH,NUM,NOCC,NVIR,IJKT,NINTMX,NFTH2,DBG)
C
         CALL RETFM(NEED)
C
C ALLOCATE MEMORY FOR  H2*HF AND CALL POLDB
C
      LDB   = 1+LOADFM
      LH2   = LDB   + NROT*3
      LHF   = LH2   + NROT*NWIDTH
      LAST  = LHF   + NUM2*3
      NEED = LAST - LDB
      CALL GETFM(NEED)
C
      SVGPAR = GOPARR
      SVDSKW = DSKWRK
      GOPARR = .FALSE.
      DSKWRK = .FALSE.
      IF(MASWRK) THEN
      CALL POLDB(X(LDB),X(LH2),X(LHF),IA,NOCC,
     *                    NVIR,NUM2,
     *                    NROT,NWIDTH,NFTH2)
      END IF
      GOPARR = SVGPAR
      DSKWRK = SVDSKW
C
      IF(DBG) THEN
         WRITE(IW,*) 'DYNAMIC  -DB-'
         CALL PRSQ(X(LDB),NXYZ,NROT,NROT)
      END IF
C
      CALL RETFM(NEED)
C
C   MEMORY ALLOCATION FOR SOLVCPDYN
      LDB   = 1+LOADFM
      LDPOL = LDB + NROT*3
      LH    = LDPOL + 6
      LZ    = LH    + NUM2*3
      LZL   = LZ    + NROT*3
      LHL   = LZL  + NVIR*3
      LTR   = LHL  + NVIR*3
      LPOL  = LTR  + NLOC*NLOC
      LVCOC = LPOL + NLOC*9*NDPFREQ
      LAST  = LVCOC+ NA*3
      NEED = LAST - LDB
      CALL GETFM(NEED)
C
      NITER=0
      NPASS=0
      SVGPAR = GOPARR
      SVDSKW = DSKWRK
      GOPARR = .FALSE.
      DSKWRK = .FALSE.
C
      IF(MASWRK) THEN
        CALL SOLVCPDYN(X(LH),X(LDB),X(LZ),NROT,NXYZ,NOCP,
     *               X(LDPOL),NWIDTH,NSOLV,
     *               NUM2,NOCC,NVIR,NFT16,NFTH2,NFTH21,
     *               X(LZL),X(LHL),X(LTR),X(LPOL),X(LVCOC),NLOC,MCORE)
      END IF
      DSKWRK = SVDSKW
      GOPARR = SVGPAR
      CALL DAWRIT(IDAF,IODA,X(LPOL),NLOC*9*NDPFREQ,401,0)
C
      CALL RETFM(NEED)
      END IF
  900 CONTINUE
      RETURN
C
 9998 FORMAT(' -A- MATRIX ')
 9988 FORMAT(1X,'       TIME FOR',A8,'= ',F12.3)
 9000 FORMAT(/5X,43(1H-)/
     *        5X,'COUPLED-PERTURBED CLOSED SHELL HARTREE-FOCK'/
     *        5X,43(1H-))
 9010 FORMAT(1X,'THE CPHF HAS',I8,' INDEPENDENT ORBITAL ROTATIONS.')
 9020 FORMAT(9X,'-B- WILL USE',I10,' WORDS,'/
     *        1X,'-A- AND -U- WILL USE',I10,' WORDS,'/
     *       10X,'  THERE ARE',I10,' WORDS AVAILABLE.')
 9030 FORMAT(1X,'IN MEMORY FORMATION OF -A- WOULD REQUIRE',I10,
     *          ' WORDS')
 9040 FORMAT(1X,'NOT ENOUGH MEMORY TO DO CPHF, NEED=',I10,' GOT=',I10)
 9050 FORMAT(1X,'CHOOSING IN MEMORY CPHF ALGORITHM')
 9060 FORMAT(1X,'CHOOSING OUT OF MEMORY CPHF ALGORITHM,'/
     *   1X,'-A- WILL BE FORMED IN',I5,' SLICES OF',I7,' COLUMNS.')
 9070 FORMAT(/1X,'*** NOT ENOUGH MEMORY TO RUN THE OUT OF MEMORY',
     *          ' CPHF CODE ***'/
     *       1X,'NWIDTH=',I10,'  NSLICE=',I10,'   NAT=',I10/
     *       1X,' NFREE=',I10,'   NNEED=',I10,'  NROT=',I10//
     *       1X,'INCREASE MEMORY IN $SYSTEM TO AT LEAST',I10/
     *       1X,'LARGER AMOUNTS OF MEMORY ARE EVEN BETTER!')
 9080 FORMAT(1X,'SOLVING FOR ALL',I4,' UNIQUE RESPONSES TOOK',
     *          I5,' ITERATIONS.')
 9090 FORMAT(1X,'CONJUGATE GRADIENT SOLVER MADE',I4,
     *          ' READS OF THE EXTERNALLY STORED -A-')
      END
C*MODULE LOCPOL    *DECK POLROHF
      SUBROUTINE POLROHF(TCSCF,SOME,OUT,NOCP,NXYZ,NUNIQ,
     *                   NFT16,MCORE,NLOC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL TCSCF,SOME,OUT,PACK2E,CHKMEM,INCORE,GOPARR,DSKWRK,MASWRK,
     *        SVGPAR,SVDSKW,POLAR
C
      DIMENSION NOCP(3)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /CHFINF/ CHFSLV,CPTYPE,NWDCHF,POLAR
      COMMON /FMCOM / X(1)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /IJPAIR/ IA(MXAO)
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C          NOTE USE OF INCORE SOLVER
      DATA MAXEQU/128/
C
C     ----- RESTRICTED OPEN SHELL CPHF CODE -----
C     BASED ON CODE WRITTEN BY MIKE SCHMIDT AT NORTH DAKOTA STATE U.,
C     IN FITS AND SPURTS, OVER THE PERIOD MAY-NOVEMBER, 1989.
C            THIS IS CODED FROM
C     Y.OSAMURA, Y.YAMAGUCHI, P.SAXE, D.J.FOX, M.A.VINCENT,
C     H.F.SCHAEFER,   J.MOL.STRUCT. 103,183-186(1983)
C            BUT SEE ALSO
C     P.SAXE, Y.YAMAGUCHI, H.F.SCHAEFER
C        J.CHEM.PHYS. 77, 5647-5654(1982)
C     Y.OSAMURA, Y.YAMAGUCHI, P.SAXE, M.A.VINCENT, J.F.GAW,
C     H.F.SCHAEFER    CHEM.PHYS. 72, 131-9(1982)
C            WITH FORMULAE SUMMARIZED IN
C     Y.YAMAGUCHI, M.FRISCH, J.GAW, H.F.SCHAEFER, J.S.BINKLEY
C        J.CHEM.PHYS. 84, 2262-2278(1986)
C
C     ----- OPEN SHELL TWO CONFIGURATION CODE -----
C     WRITTEN BY JOHN MONTGOMERY AND MIKE SCHMIDT BETWEEN
C     JUNE 22 AND JULY 22, 1990.
C            THIS IS CODED FROM
C     M.DURAN, Y.YAMAGUCHI, H.F.SCHAEFER III
C     J.PHYS.CHEM. 92, 3070-3075(1988).
C            BUT SEE ALSO
C     M.DURAN, Y.YAMAGUCHI, R.B.REMINGTON, H.F.SCHAEFER
C     CHEM.PHYS. 122, 201-231(1988).
C
      IF(SOME  .AND.  .NOT.TCSCF) WRITE(IW,9000)
      IF(SOME  .AND.       TCSCF) WRITE(IW,9010)
      TIM1 = 0.0D+00
      DO 100 I=1,NUM
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     ----- INITIALIZE -NCONF- FOR A VIRTUAL SHELL -----
C
      DO 150 I=NORB+1,NUM
         NCONF(I) = NHAM+1
  150 CONTINUE
C
      IF(OUT) THEN
         WRITE(IW,8150) NCO,NSETO,NOPEN,NPAIR,NORB,NHAM
         WRITE(IW,8160) (NCONF(I),I=1,NORB)
         WRITE(IW,8170) '    F=',(F(I),I=1,NHAM)
         NTRI = (NHAM*NHAM+NHAM)/2
         WRITE(IW,8170) 'ALPHA=',(ALPHA(I),I=1,NTRI)
         WRITE(IW,8170) ' BETA=',( BETA(I),I=1,NTRI)
      END IF
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      LOCS=L1-MCORE
C
C     ----- DETERMINE THE NUMBER OF INTER-SHELL ROTATIONS -----
C
      NROT=0
      DO 210 I=1,L0
         NCI = NCONF(I)
         DO 200 J=1,I
            IF(NCONF(J).NE.NCI) NROT=NROT+1
  200    CONTINUE
  210 CONTINUE
      IF(SOME) WRITE(IW,9020) NROT
      IF(TCSCF) THEN
         NROT = NROT+2
         IF (MASWRK) WRITE(IW,9030)
      END IF
C
C     FIRST, ESTABLISH MEMORY NEEDS.
C     THEN, COME BACK TO 300 AND BEGIN TO DO THE ACTUAL WORK.
C
      CHKMEM=.TRUE.
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      IF(NWDCHF.GT.0) NGOTMX=MIN(NGOTMX,NWDCHF)
      IF (MASWRK) WRITE(IW,9040) NGOTMX
      NEED1 = 0
      NEED2 = 0
      GO TO 310
C
  300 CONTINUE
      CHKMEM=.FALSE.
      NEED = MAX(NEED1,NEED2)
      IF(NEED.GT.NGOTMX) THEN
         IF (MASWRK) THEN
            WRITE(IW,*)
     *         '**** INSUFFICIENT MEMORY FOR CP-ROHF *****'
            WRITE(IW,*) 'NEED',NEED,' WORDS.'
            WRITE(IW,*) ' GOT',NGOTMX,' WORDS.'
         END IF
         CALL ABRT
         STOP
      END IF
      IF(EXETYP.EQ.CHECK) GO TO 800
      TIM1=0.0D+00
      CALL TSECND(TIM1)
C
C     ----- TRANSFORM DIPOLE INTEGRALS TO MO BASIS -----
C
  310 CONTINUE
      LHF    = 1      + LOADFM
      LTHF   = LHF    + L2
      LV     = LTHF   + L2*NXYZ
      LWRK   = LV     + L3
      LAST   = LWRK   + L1
      NEED1  = LAST - LOADFM -1
      IF(CHKMEM) GO TO 400
      CALL GETFM(NEED1)
C
      CALL OPOLB(X(LHF),X(LTHF),X(LV),X(LWRK),L1,L2,L3,NXYZ)
C
      CALL RETFM(NEED1)
      CALL TSECND(TIM2)
      THYME = TIM2-TIM1
      IF (MASWRK) WRITE(IW,9050)
     *   'TRANSFORMING DIPOLE INTEGRALS TO MO BASIS',THYME
      TIM1 = TIM2
C
C     ----- MEMORY FOR THE CPHF -----
C     NOTE THAT -SA- FROM THE PREVIOUS STEP MUST BE THE 1ST ARRAY
C
  400 CONTINUE
      LXX    = 1      + LOADFM
      LIXX   = LXX    + NINTMX
      LEPS   = LIXX   + NINTMX
      LZETA  = LEPS   + L3
      LALP   = LZETA  + L2*(NHAM+1)
      LBET   = LALP   + L2
      LIROT  = LBET   + L2
      LHMO   = LIROT  + L3/NWDVAR + 1
      LW     = LHMO   + L2
      LDDM   = LW     + L1*NORB*NXYZ
      LDDMA  = LDDM   + NAT*9
      LSK    = LDDMA  + NAT*9
      LTHF   = LSK    + L1
      LDLPOL = LTHF   + L2*3
      LAST   = LDLPOL + NLOC*9
C
      IF(TCSCF) THEN
         LEPSIJ = LAST
         LAST   = LEPSIJ + 3*L3
      ELSE
         LEPSIJ = LAST
      END IF
C
      LB     = LAST
      LU     = LB     + NROT*NXYZ
      LAII   = LU     + NROT*NXYZ
      LUL    = LAII   + NROT
      LHL    = LUL    + LOCS*3
      LTR    = LHL    + LOCS*3
      LPOL   = LTR    + NLOC*NLOC
      LVCOC  = LPOL   + NLOC*9
      IROT   = LVCOC  + NA*3
      LAST   = IROT   + L1*L1
C
C        THE VECTOR VERSION NEEDS 4 EXTRA ARRAYS FOR TFTRI LATER.
C
      NEEDX  = LAST -LOADFM -1 + (3*L3+L1)
C
C        DECIDE UPON IN CORE OR OUT OF CORE METHOD
C
      NSOLV = NXYZ
      IF(NROT.LE.MAXEQU) NSOLV=1
C
      DMY = 2.0D+00**31 - 2*NSOLV - 4*NROT*NSOLV - 1000
      N32MAX = INT(SQRT(DMY))
      IF(NWDVAR.EQ.1  .OR.  NROT.LE.N32MAX) THEN
         NEEDIN = NEEDX + NROT*NROT + 2*NSOLV + 4*NROT*NSOLV
         INCORE = NEEDIN.LT.NGOTMX
      ELSE
         NEEDIN = 0
         INCORE = .FALSE.
      END IF
C
C        TO RUN SMALL CASES OUT OF CORE (IN 2 SLICES) FOR
C        TESTING PURPOSES, UNCOMMENT THE NEXT FEW LINES.
C            IF(INCORE) THEN
C               INCORE = .FALSE.
C               NGOTMX = NEEDX +2*NXYZ +4*NROT*NXYZ +(NROT*NROT+1)/2
C               IF(TCSCF) NGOTMX = NGOTMX + 2*NROT
C            END IF
C
      IF(INCORE) THEN
         NWIDTH = NROT
         NSLICE = 1
         LACI   = LAST
         LWRK1  = LAST
         LWRK2  = LWRK1 + NROT*NSOLV
         LWRK3  = LWRK2 + NROT*NSOLV
         LWRK4  = LWRK3 + NROT*NSOLV
         LBNRM  = LWRK4 + NROT*NSOLV
         LRESD  = LBNRM + NSOLV
         LA     = LRESD + NSOLV
         LAST   = LA    + NROT*NROT
         IF(CHKMEM .AND. MASWRK) WRITE(IW,9060)
      ELSE
         LACI   = LAST
         IF(TCSCF) LAST = LACI  + 2*NROT
         LWRK1  = LAST
         LWRK2  = LWRK1 + NROT*NSOLV
         LWRK3  = LWRK2 + NROT*NSOLV
         LWRK4  = LWRK3 + NROT*NSOLV
         LBNRM  = LWRK4 + NROT*NSOLV
         LRESD  = LBNRM + NSOLV
         LAST   = LRESD + NSOLV
         NEEDX  = LAST -LOADFM -1 + (3*L3+L1)
         NFREE  = NGOTMX - NEEDX
         NWIDTH = NFREE/NROT
         IF(NWIDTH.GT.0) THEN
            NSLICE = (NROT-1)/NWIDTH + 1
         ELSE
            NSLICE = -1
         END IF
         LA     = LAST
         LAST   = LA    + NROT*NWIDTH
         IF(NWIDTH.LE.0  .OR.  (NSLICE.GT.NAT .AND. NAT.GT.5)) THEN
            IF (MASWRK) WRITE(IW,9070)
     *         NWIDTH,NSLICE,NAT,NFREE,NEEDX,NROT
            CALL ABRT
            STOP
         END IF
         IF(CHKMEM .AND. MASWRK) WRITE(IW,9080) NSLICE,NWIDTH
      END IF
C
      NEED2  = LAST - LOADFM -1
      IF(SOME  .AND.  CHKMEM) THEN
         WRITE(IW,9090) 'DS/DA TRANSFORM',NEED1
         WRITE(IW,9090) ' CPHF EQUATIONS',NEED2
         IF(.NOT.INCORE) WRITE(IW,9100) NEEDIN
      END IF
      IF(CHKMEM) GO TO 300
      CALL GETFM(NEED2)
C
C     ----- READ BARE NUCLEUS HAMILTONIAN, PUT IN MO BASIS -----
C     TEMP ARRAYS: LALP=HAO, LEPS=V, LBET=WRK
C
      CALL DAREAD(IDAF,IODA,X(LALP),L2,11,0)
      CALL DAREAD(IDAF,IODA,X(LEPS),L3,15,0)
      CALL TFTRI(X(LHMO),X(LALP),X(LEPS),X(LBET),L0,L1,L1)
C
C     ----- SET UP THE OPEN SHELL/TCSCF CPHF EQUATIONS -----
C
      CALL ROCPA(TCSCF,X(LXX),X(LIXX),NINTMX,IA,X(LALP),X(LBET),
     *            X(LEPS),X(LZETA),X(LHMO),X(LIROT),X(LA),X(LB),
     *            X(LW),X(LEPSIJ),X(LACI),L0,L1,L2,L3,NROT,NXYZ,
     *            NHAM,NORB,OUT,INCORE,X(LTHF))
      IF(.NOT.INCORE) THEN
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQOPN(NFT16,'WORK16','UNKNOWN',.FALSE.,'UNFORMATTED')
         DSKWRK = SVDSKW
         CALL ROCPOM(TCSCF,X(LXX),X(LIXX),NINTMX,IA,X(LALP),X(LBET),
     *               X(LEPS),X(LZETA),X(LIROT),X(LA),X(LACI),X(LAII),
     *               L0,L1,L2,NROT,NHAM,NWIDTH,NFT16,OUT)
      END IF
C
      CALL TSECND(TIM2)
      THYME = TIM2-TIM1
      IF (MASWRK) WRITE(IW,9050)
     *   ' SETTING UP THE CPHF EQUATIONS',THYME
      CALL FLSHBF(IW)
      TIM1 = TIM2
C
C     ----- SOLVE THE OPEN SHELL/TCSCF CPHF EQUATIONS -----
C
      NITER=0
      NPASS=0
      SVGPAR = GOPARR
      SVDSKW = DSKWRK
      GOPARR = .FALSE.
      DSKWRK = .FALSE.
      IF(MASWRK) THEN
      CALL SOLVCP(X(LA),X(LB),X(LU),NROT,NXYZ,NOCP,X(LAII),
     *            X(LWRK1),X(LWRK2),X(LWRK3),X(LWRK4),
     *            X(LBNRM),X(LRESD),NITER,NPASS,NWIDTH,NSOLV,
     *            INCORE,NFT16,MAXEQU)
      END IF
      GOPARR = SVGPAR
      DSKWRK = SVDSKW
      IF(OUT) THEN
         WRITE(IW,*) 'CPHF RESPONSE VECTORS -U-'
         CALL PRSQ(X(LU),NXYZ,NROT,NROT)
      END IF
C
      CALL TSECND(TIM2)
      THYME = TIM2-TIM1
      IF (MASWRK) THEN
         WRITE(IW,9050) 'SOLVING FOR THE CPHF RESPONSES',THYME
         IF(NITER.NE.0) WRITE(IW,9110) NUNIQ,NITER
         IF(.NOT.INCORE) WRITE(IW,9120) NPASS
      END IF
      TIM1 = TIM2
C
      CALL LOAPOL(X(LU),X(LUL),X(LTHF),X(LHL),X(LTR),X(LPOL),X(LVCOC),
     *           IA,NXYZ,L2,MCORE,NLOC,LOCS,NROT,X(IROT))
C
      CALL RETFM(NEED2)
  800 CONTINUE
      RETURN
C
 8150 FORMAT(1X,'NCO,NSETO,NOPEN,NPAIR,NORB,NHAM=',6I4)
 8170 FORMAT(1X,A6,10F10.5/(7X,10F10.5))
 8160 FORMAT(1X,'NCONF=',10I5/(7X,10I5))
C
 9000 FORMAT(/5X,41(1H-)/
     *        5X,'COUPLED-PERTURBED OPEN SHELL HARTREE-FOCK'/5X,41(1H-))
 9010 FORMAT(/5X,48(1H-)/
     *        5X,'COUPLED-PERTURBED TWO CONFIGURATION HARTREE-FOCK'/
     *        5X,48(1H-))
 9020 FORMAT(1X,'THE CPHF HAS',I8,' INDEPENDENT ORBITAL ROTATIONS.')
 9030 FORMAT(1X,'THE CPHF HAS 2 ADDITIONAL CSF ROTATIONS.')
 9040 FORMAT(1X,'THE CPHF CODE HAS',I12,' WORDS OF MEMORY AVAILABLE.')
 9050 FORMAT(5X,A30,' TOOK',F12.1,' SECONDS')
 9060 FORMAT(1X,'CHOOSING IN MEMORY CPHF ALGORITHM')
 9070 FORMAT(/1X,'*** NOT ENOUGH MEMORY TO RUN THE OUT OF MEMORY',
     *          ' CPHF CODE ***'/
     *       1X,'NWIDTH=',I10,'  NSLICE=',I10,'   NAT=',I10/
     *       1X,' NFREE=',I10,'   NEEDX=',I10,'  NROT=',I10//
     *       1X,'INCREASE THE VALUE OF MEMORY IN $CONTRL SO THAT'/
     *       1X,'NFREE=MEMORY-NEEDX AT LEAST EXCEEDS NROT*NROT/NAT'//
     *       1X,'LARGER AMOUNTS OF MEMORY ARE EVEN BETTER!')
 9080 FORMAT(1X,'CHOOSING OUT OF MEMORY CPHF ALGORITHM'/
     *       1X,'THE -A- MATRIX WILL BE FORMED IN',I5,
     *          ' SLICES OF',I7,' COLUMNS.')
 9090 FORMAT(1X,A15,' STEP WILL REQUIRE',I12,' WORDS OF MEMORY.')
 9100 FORMAT(1X,'(THE IN CORE ALGORITHM WOULD NEED',I12,' WORDS)')
 9110 FORMAT(1X,'SOLVING FOR ALL',I5,' UNIQUE RESPONSES TOOK',
     *          I5,' ITERATIONS.')
 9120 FORMAT(1X,'CONJUGATE GRADIENT SOLVER MADE',I4,
     *          ' READS OF THE EXTERNALLY STORED -A- MATRIX')
      END
C*MODULE LOCPOL    *DECK ROCPA
      SUBROUTINE ROCPA(TCSCF,XX,IX,NINTMX,IA,ALP,BET,EPS,ZETA,
     *                  HMO,IROT,A,B,W,EPSIJ,
     *                  ACI,L0,L1,L2,L3,NROT,NXYZ,NHAM,NORB,
     *                  OUT,INCORE,HF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL INCORE,OUT,TCSCF,IOCC,JOCC,KOCC,
     *        ICOP,JCOP,KCOP,LCOP,IPR1,JPR1,KPR1,LPR1,
     *        IPR2,JPR2,KPR2,LPR2,GOPARR,MASWRK,DSKWRK
C
      DIMENSION XX(NINTMX),IX(NINTMX),IA(L1),ALP(L2),BET(L2),
     *          EPS(L1,L1),ZETA(L2,NHAM),HMO(L2),IROT(L0,L0),
     *          A(NROT,NROT),B(NROT,NXYZ),
     *          W(L1,NORB,NXYZ),
     *          EPSIJ(L1,L1,3),
     *          ACI(NROT,2),HF(L2,3)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /FUNCT / ETOT,EG(3,MXATM)
      COMMON /GVBCI / CIHAM(91),CILOW(12),JPAIR,KONE,KTWO,KCORB(2,12)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORBB,NCONF(MXAO),NHAMM
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
C
C     ----- SET UP GENERAL OPEN SHELL CPHF PROBLEM -----
C     IF -INCORE- IS TRUE, SET UP THE ENTIRE ORBITAL HESSIAN -A-
C     IN MEMORY, -ACI- IS NOT USED IN THIS CASE.
C     OTHERWISE, -A- IS NOT COMPUTED, ALTHOUGH THE POSSIBLE TCSCF
C     CI TERMS ARE COMPUTED IN -ACI-.
C
C     ALL THE OTHER ARRAYS SUCH AS -EPS- OR -ZETA- AND -W- AND SO ON
C     ARE ALWAYS FORMED IN MEMORY BY THIS ROUTINE.
C
      VALHF = ZERO
      VAL4  = ZERO
      AFACT = ZERO
      BFACT = ZERO
      ICOP = .FALSE.
      JCOP = .FALSE.
      KCOP = .FALSE.
      LCOP = .FALSE.
      IPR1 = .FALSE.
      JPR1 = .FALSE.
      KPR1 = .FALSE.
      LPR1 = .FALSE.
      IPR2 = .FALSE.
      JPR2 = .FALSE.
      KPR2 = .FALSE.
      LPR2 = .FALSE.
      NCLOP = NCO+NOPEN
C
C     FIRST FILL ALPHA AND BETA COUPLING COEFFICIENT MATRICES
C
      CALL VCLR(ALP,1,L2)
      CALL VCLR(BET,1,L2)
      IJ=0
      DO 110 I=1,NORB
         DO 100 J=1,I
            IJ = IJ + 1
            ISH = NCONF(I)
            JSH = NCONF(J)
            IJSH = IA(ISH) + JSH
            ALP(IJ) = ALPHA(IJSH)
            BET(IJ) =  BETA(IJSH)
  100    CONTINUE
  110 CONTINUE
C
C     ----- SET UP ROTATION INDEX ARRAY -----
C
      MROT=0
      DO 220 I=1,L0
         ISH = NCONF(I)
         DO 210 J=1,I
            JSH = NCONF(J)
            IF(JSH.EQ.ISH) THEN
               IROT(I,J) = 0
               IROT(J,I) = 0
            ELSE
               MROT=MROT+1
               IROT(I,J) = MROT
               IROT(J,I) = MROT
            END IF
  210    CONTINUE
  220 CONTINUE
      IF(TCSCF) MROT = MROT+2
      IF(MROT.NE.NROT) CALL ABRT
C
C     ZERO OUT CPHF -A- AND -B- MATRICES, AND THE -W- MATRIX
C     ZERO OUT LAGRANGIAN -EPS- AND GENERALIZED LAGRANGIANS -ZETA-
C     FOR TCSCF, ZERO OUT BARE LAGRANGIANS -EPSIJ-
C
      IF(INCORE) THEN
         CALL VCLR(A,1,NROT*NROT)
      ELSE
         IF(TCSCF) CALL VCLR(ACI,1,NROT*2)
      END IF
      CALL VCLR(B,1,NROT*NXYZ)
      CALL VCLR(EPS,1,L3)
      CALL VCLR(ZETA,1,L2*NHAM)
      IF(TCSCF) CALL VCLR(EPSIJ,1,3*L3)
C
C     ----- READ IN 2 ELECTRON MO INTEGRALS -----
C
      CALL SEQREW(IJKT)
      IF (MASWRK) READ(IJKT)
C
  300 CONTINUE
      CALL PREAD(IJKT,XX,IX,NX,NINTMX)
      IF(NX.EQ.0) GO TO 900
      MX=ABS(NX)
      IF(MX.GT.NINTMX) CALL ABRT
C
C     ----- UNPACK LABEL AND PROCESS INTEGRAL -----
C     TRANSFORMED INTEGRALS ASSUMED IN REVERSE CANONICAL ORDER
C
      DO 800 M=1,MX
      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
      IF(J.EQ.I           ) VAL=VAL*HALF
      IF(L.EQ.K           ) VAL=VAL*HALF
      IF(K.EQ.I.AND.L.EQ.J) VAL=VAL*HALF
      VAL2 = VAL+VAL
C
      IJ = IA(I)+J
      IK = IA(I)+K
      IL = IA(I)+L
      KL = IA(K)+L
      JK = IA(MAX(J,K))+MIN(J,K)
      JL = IA(MAX(J,L))+MIN(J,L)
C
      ISH = NCONF(I)
      JSH = NCONF(J)
      KSH = NCONF(K)
      LSH = NCONF(L)
C
      IOCC = I.LE.NORB
      JOCC = J.LE.NORB
      KOCC = K.LE.NORB
C
      IF(TCSCF) THEN
         VALHF = VAL*HALF
         VAL4 = VAL2 + VAL2
C
         ICOP = I.LE.NCLOP
         JCOP = J.LE.NCLOP
         KCOP = K.LE.NCLOP
         LCOP = L.LE.NCLOP
C
         IPR1 = I.EQ.NORB-1
         JPR1 = J.EQ.NORB-1
         KPR1 = K.EQ.NORB-1
         LPR1 = L.EQ.NORB-1
C
         IPR2 = I.EQ.NORB
         JPR2 = J.EQ.NORB
         KPR2 = K.EQ.NORB
         LPR2 = L.EQ.NORB
      END IF
C
C                  COULOMB CONTRIBUTIONS TO -EPS- AND -ZETA-
C                                             (I,I//K,L)
C
      IF(IOCC  .AND.  I.EQ.J) THEN
         IF(TCSCF) THEN
            IF(ICOP) THEN
               IF(KCOP) THEN
                  EPSIJ(K,L,1) = EPSIJ(K,L,1) + ALP(IK)*VAL2
                  EPSIJ(K,L,3) = EPSIJ(K,L,3) + ALP(IK)*VAL2
               ELSE IF(KPR1) THEN
                  EPSIJ(K,L,1) = EPSIJ(K,L,1) + F(NCONF(I))*VAL4
               ELSE IF(KPR2) THEN
                  EPSIJ(K,L,3) = EPSIJ(K,L,3) + F(NCONF(I))*VAL4
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,K,1) = EPSIJ(L,K,1) + ALP(IL)*VAL2
                  EPSIJ(L,K,3) = EPSIJ(L,K,3) + ALP(IL)*VAL2
               ELSE IF(LPR1) THEN
                  EPSIJ(L,K,1) = EPSIJ(L,K,1) + F(NCONF(I))*VAL4
               ELSE IF(LPR2) THEN
                  EPSIJ(L,K,3) = EPSIJ(L,K,3) + F(NCONF(I))*VAL4
               END IF
            ELSE IF(IPR1) THEN
               IF(KCOP) THEN
                  EPSIJ(K,L,1) = EPSIJ(K,L,1) + F(NCONF(K))*VAL4
               ELSE IF(KPR1) THEN
                  EPSIJ(K,L,1) = EPSIJ(K,L,1) + VAL2
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,K,1) = EPSIJ(L,K,1) + F(NCONF(L))*VAL4
               ELSE IF(LPR1) THEN
                  EPSIJ(L,K,1) = EPSIJ(L,K,1) + VAL2
               END IF
            ELSE IF(IPR2) THEN
               IF(KCOP) THEN
                  EPSIJ(K,L,3) = EPSIJ(K,L,3) + F(NCONF(K))*VAL4
               ELSE IF(KPR2) THEN
                  EPSIJ(K,L,3) = EPSIJ(K,L,3) + VAL2
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,K,3) = EPSIJ(L,K,3) + F(NCONF(L))*VAL4
               ELSE IF(LPR2) THEN
                  EPSIJ(L,K,3) = EPSIJ(L,K,3) + VAL2
               END IF
            END IF
         ELSE
            EPS(K,L)=EPS(K,L) + ALP(IK)*VAL2
            EPS(L,K)=EPS(L,K) + ALP(IL)*VAL2
         END IF
         DO 410 MSH=1,NHAM
            IM = IA(MAX(ISH,MSH))+MIN(ISH,MSH)
            ZETA(KL,MSH)=ZETA(KL,MSH) + ALPHA(IM)*VAL2
  410    CONTINUE
      END IF
C                                             (I,J//K,K)
      IF(KOCC  .AND.  K.EQ.L) THEN
         IF(TCSCF) THEN
            IF(KCOP) THEN
               IF(ICOP) THEN
                  EPSIJ(I,J,1) = EPSIJ(I,J,1) + ALP(IK)*VAL2
                  EPSIJ(I,J,3) = EPSIJ(I,J,3) + ALP(IK)*VAL2
               ELSE IF(IPR1) THEN
                  EPSIJ(I,J,1) = EPSIJ(I,J,1) + F(NCONF(K))*VAL4
               ELSE IF(IPR2) THEN
                  EPSIJ(I,J,3) = EPSIJ(I,J,3) + F(NCONF(K))*VAL4
               END IF
               IF(JCOP) THEN
                  EPSIJ(J,I,1) = EPSIJ(J,I,1) + ALP(JK)*VAL2
                  EPSIJ(J,I,3) = EPSIJ(J,I,3) + ALP(JK)*VAL2
               ELSE IF(JPR1) THEN
                  EPSIJ(J,I,1) = EPSIJ(J,I,1) + F(NCONF(K))*VAL4
               ELSE IF(JPR2) THEN
                  EPSIJ(J,I,3) = EPSIJ(J,I,3) + F(NCONF(K))*VAL4
               END IF
            ELSE IF(KPR1) THEN
               IF(ICOP) THEN
                  EPSIJ(I,J,1) = EPSIJ(I,J,1) + F(NCONF(I))*VAL4
               ELSE IF(IPR1) THEN
                  EPSIJ(I,J,1) = EPSIJ(I,J,1) + VAL2
               END IF
               IF(JCOP) THEN
                  EPSIJ(J,I,1) = EPSIJ(J,I,1) + F(NCONF(J))*VAL4
               ELSE IF(JPR1) THEN
                  EPSIJ(J,I,1) = EPSIJ(J,I,1) + VAL2
               END IF
            ELSE IF(KPR2) THEN
               IF(ICOP) THEN
                  EPSIJ(I,J,3) = EPSIJ(I,J,3) + F(NCONF(I))*VAL4
               ELSE IF(IPR2) THEN
                  EPSIJ(I,J,3) = EPSIJ(I,J,3) + VAL2
               END IF
               IF(JCOP) THEN
                  EPSIJ(J,I,3) = EPSIJ(J,I,3) + F(NCONF(J))*VAL4
               ELSE IF(JPR2) THEN
                  EPSIJ(J,I,3) = EPSIJ(J,I,3) + VAL2
               END IF
            END IF
         ELSE
            EPS(I,J)=EPS(I,J) + ALP(IK)*VAL2
            EPS(J,I)=EPS(J,I) + ALP(JK)*VAL2
         END IF
         DO 420 MSH=1,NHAM
            KM = IA(MAX(MSH,KSH))+MIN(MSH,KSH)
            ZETA(IJ,MSH)=ZETA(IJ,MSH) + ALPHA(KM)*VAL2
  420    CONTINUE
      END IF
C
C                  EXCHANGE CONTRIBUTIONS TO -EPS- AND -ZETA-
C                                             (I,J//I,L)
      IF(IOCC  .AND.  I.EQ.K) THEN
         IF(TCSCF) THEN
            IF(ICOP) THEN
               IF(JCOP) THEN
                  EPSIJ(J,L,1) = EPSIJ(J,L,1) + BET(IJ)*VAL
                  EPSIJ(J,L,3) = EPSIJ(J,L,3) + BET(IJ)*VAL
               ELSE IF(JPR1) THEN
                  EPSIJ(J,L,1) = EPSIJ(J,L,1) - F(NCONF(I))*VAL
               ELSE IF(JPR2) THEN
                  EPSIJ(J,L,3) = EPSIJ(J,L,3) - F(NCONF(I))*VAL
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,J,1) = EPSIJ(L,J,1) + BET(IL)*VAL
                  EPSIJ(L,J,3) = EPSIJ(L,J,3) + BET(IL)*VAL
               ELSE IF(LPR1) THEN
                  EPSIJ(L,J,1) = EPSIJ(L,J,1) - F(NCONF(I))*VAL
               ELSE IF(LPR2) THEN
                  EPSIJ(L,J,3) = EPSIJ(L,J,3) - F(NCONF(I))*VAL
               END IF
            ELSE IF(IPR1) THEN
               IF(JCOP) THEN
                  EPSIJ(J,L,1) = EPSIJ(J,L,1) - F(NCONF(J))*VAL
               ELSE IF(JPR2) THEN
                  EPSIJ(J,L,2) = EPSIJ(J,L,2) + VALHF
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,J,1) = EPSIJ(L,J,1) - F(NCONF(L))*VAL
               ELSE IF(LPR2) THEN
                  EPSIJ(L,J,2) = EPSIJ(L,J,2) + VALHF
               END IF
            ELSE IF(IPR2) THEN
               IF(JCOP) THEN
                  EPSIJ(J,L,3) = EPSIJ(J,L,3) - F(NCONF(J))*VAL
               ELSE IF(JPR1) THEN
                  EPSIJ(J,L,2) = EPSIJ(J,L,2) + VALHF
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,J,3) = EPSIJ(L,J,3) - F(NCONF(L))*VAL
               ELSE IF(LPR1) THEN
                  EPSIJ(L,J,2) = EPSIJ(L,J,2) + VALHF
               END IF
            END IF
         ELSE
            EPS(J,L)=EPS(J,L) + BET(IJ)*VAL
            EPS(L,J)=EPS(L,J) + BET(IL)*VAL
         END IF
         DO 510 MSH=1,NHAM
            IM = IA(MAX(ISH,MSH))+MIN(ISH,MSH)
            ZETA(JL,MSH)=ZETA(JL,MSH) + BETA(IM)*VAL
  510    CONTINUE
      END IF
C                                             (I,J//K,I)
      IF(IOCC  .AND.  I.EQ.L) THEN
         IF(TCSCF) THEN
            IF(ICOP) THEN
               IF(JCOP) THEN
                  EPSIJ(J,K,1) = EPSIJ(J,K,1) + BET(IJ)*VAL
                  EPSIJ(J,K,3) = EPSIJ(J,K,3) + BET(IJ)*VAL
               ELSE IF(JPR1) THEN
                  EPSIJ(J,K,1) = EPSIJ(J,K,1) - F(NCONF(I))*VAL
               ELSE IF(JPR2) THEN
                  EPSIJ(J,K,3) = EPSIJ(J,K,3) - F(NCONF(I))*VAL
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,J,1) = EPSIJ(K,J,1) + BET(IK)*VAL
                  EPSIJ(K,J,3) = EPSIJ(K,J,3) + BET(IK)*VAL
               ELSE IF(KPR1) THEN
                  EPSIJ(K,J,1) = EPSIJ(K,J,1) - F(NCONF(I))*VAL
               ELSE IF(KPR2) THEN
                  EPSIJ(K,J,3) = EPSIJ(K,J,3) - F(NCONF(I))*VAL
               END IF
            ELSE IF(IPR1) THEN
               IF(JCOP) THEN
                  EPSIJ(J,K,1) = EPSIJ(J,K,1) - F(NCONF(J))*VAL
               ELSE IF(JPR2) THEN
                  EPSIJ(J,K,2) = EPSIJ(J,K,2) + VALHF
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,J,1) = EPSIJ(K,J,1) - F(NCONF(K))*VAL
               ELSE IF(KPR2) THEN
                  EPSIJ(K,J,2) = EPSIJ(K,J,2) + VALHF
               END IF
            ELSE IF(IPR2) THEN
               IF(JCOP) THEN
                  EPSIJ(J,K,3) = EPSIJ(J,K,3) - F(NCONF(J))*VAL
               ELSE IF(JPR1) THEN
                  EPSIJ(J,K,2) = EPSIJ(J,K,2) + VALHF
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,J,3) = EPSIJ(K,J,3) - F(NCONF(K))*VAL
               ELSE IF(KPR1) THEN
                  EPSIJ(K,J,2) = EPSIJ(K,J,2) + VALHF
               END IF
            END IF
         ELSE
            EPS(J,K)=EPS(J,K) + BET(IJ)*VAL
            EPS(K,J)=EPS(K,J) + BET(IK)*VAL
         END IF
         DO 520 MSH=1,NHAM
            IM = IA(MAX(ISH,MSH))+MIN(ISH,MSH)
            ZETA(JK,MSH)=ZETA(JK,MSH) + BETA(IM)*VAL
  520    CONTINUE
      END IF
C                                             (I,J//J,L)
      IF(JOCC  .AND.  J.EQ.K) THEN
         IF(TCSCF) THEN
            IF(JCOP) THEN
               IF(ICOP) THEN
                  EPSIJ(I,L,1) = EPSIJ(I,L,1) + BET(IJ)*VAL
                  EPSIJ(I,L,3) = EPSIJ(I,L,3) + BET(IJ)*VAL
               ELSE IF(IPR1) THEN
                  EPSIJ(I,L,1) = EPSIJ(I,L,1) - F(NCONF(J))*VAL
               ELSE IF(IPR2) THEN
                  EPSIJ(I,L,3) = EPSIJ(I,L,3) - F(NCONF(J))*VAL
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,I,1) = EPSIJ(L,I,1) + BET(JL)*VAL
                  EPSIJ(L,I,3) = EPSIJ(L,I,3) + BET(JL)*VAL
               ELSE IF(LPR1) THEN
                  EPSIJ(L,I,1) = EPSIJ(L,I,1) - F(NCONF(J))*VAL
               ELSE IF(LPR2) THEN
                  EPSIJ(L,I,3) = EPSIJ(L,I,3) - F(NCONF(J))*VAL
               END IF
            ELSE IF(JPR1) THEN
               IF(ICOP) THEN
                  EPSIJ(I,L,1) = EPSIJ(I,L,1) - F(NCONF(I))*VAL
               ELSE IF(IPR2) THEN
                  EPSIJ(I,L,2) = EPSIJ(I,L,2) + VALHF
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,I,1) = EPSIJ(L,I,1) - F(NCONF(L))*VAL
               ELSE IF(LPR2) THEN
                  EPSIJ(L,I,2) = EPSIJ(L,I,2) + VALHF
               END IF
            ELSE IF(JPR2) THEN
               IF(ICOP) THEN
                  EPSIJ(I,L,3) = EPSIJ(I,L,3) - F(NCONF(I))*VAL
               ELSE IF(IPR1) THEN
                  EPSIJ(I,L,2) = EPSIJ(I,L,2) + VALHF
               END IF
               IF(LCOP) THEN
                  EPSIJ(L,I,3) = EPSIJ(L,I,3) - F(NCONF(L))*VAL
               ELSE IF(LPR1) THEN
                  EPSIJ(L,I,2) = EPSIJ(L,I,2) + VALHF
               END IF
            END IF
         ELSE
            EPS(I,L)=EPS(I,L) + BET(IJ)*VAL
            EPS(L,I)=EPS(L,I) + BET(JL)*VAL
         END IF
         DO 530 MSH=1,NHAM
            JM = IA(MAX(JSH,MSH))+MIN(JSH,MSH)
            ZETA(IL,MSH)=ZETA(IL,MSH) + BETA(JM)*VAL
  530    CONTINUE
      END IF
C                                             (I,J//K,J)
      IF(JOCC  .AND.  J.EQ.L) THEN
         IF(TCSCF) THEN
            IF(JCOP) THEN
               IF(ICOP) THEN
                  EPSIJ(I,K,1) = EPSIJ(I,K,1) + BET(IJ)*VAL
                  EPSIJ(I,K,3) = EPSIJ(I,K,3) + BET(IJ)*VAL
               ELSE IF(IPR1) THEN
                  EPSIJ(I,K,1) = EPSIJ(I,K,1) - F(NCONF(J))*VAL
               ELSE IF(IPR2) THEN
                  EPSIJ(I,K,3) = EPSIJ(I,K,3) - F(NCONF(J))*VAL
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,I,1) = EPSIJ(K,I,1) + BET(JK)*VAL
                  EPSIJ(K,I,3) = EPSIJ(K,I,3) + BET(JK)*VAL
               ELSE IF(KPR1) THEN
                  EPSIJ(K,I,1) = EPSIJ(K,I,1) - F(NCONF(J))*VAL
               ELSE IF(KPR2) THEN
                  EPSIJ(K,I,3) = EPSIJ(K,I,3) - F(NCONF(J))*VAL
               END IF
            ELSE IF(JPR1) THEN
               IF(ICOP) THEN
                  EPSIJ(I,K,1) = EPSIJ(I,K,1) - F(NCONF(I))*VAL
               ELSE IF(IPR2) THEN
                  EPSIJ(I,K,2) = EPSIJ(I,K,2) + VALHF
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,I,1) = EPSIJ(K,I,1) - F(NCONF(K))*VAL
               ELSE IF(KPR2) THEN
                  EPSIJ(K,I,2) = EPSIJ(K,I,2) + VALHF
               END IF
            ELSE IF(JPR2) THEN
               IF(ICOP) THEN
                  EPSIJ(I,K,3) = EPSIJ(I,K,3) - F(NCONF(I))*VAL
               ELSE IF(IPR1) THEN
                  EPSIJ(I,K,2) = EPSIJ(I,K,2) + VALHF
               END IF
               IF(KCOP) THEN
                  EPSIJ(K,I,3) = EPSIJ(K,I,3) - F(NCONF(K))*VAL
               ELSE IF(KPR1) THEN
                  EPSIJ(K,I,2) = EPSIJ(K,I,2) + VALHF
               END IF
            END IF
         ELSE
            EPS(I,K)=EPS(I,K) + BET(IJ)*VAL
            EPS(K,I)=EPS(K,I) + BET(JK)*VAL
         END IF
         DO 540 MSH=1,NHAM
            JM = IA(MAX(JSH,MSH))+MIN(JSH,MSH)
            ZETA(IK,MSH)=ZETA(IK,MSH) + BETA(JM)*VAL
  540    CONTINUE
      END IF
C
C                  INTEGRAL CONTRIBUTIONS TO -A- MATRIX
C                                              (I,J//K,L)
      IF(.NOT.INCORE) GO TO 800
      IF(ISH.NE.JSH  .AND.  KSH.NE.LSH) THEN
         IJROT = IROT(I,J)
         KLROT = IROT(K,L)
         IF(IJROT.LT.KLROT) THEN
            NN    = IJROT
            IJROT = KLROT
            KLROT = NN
         END IF
         AFACT = ALP(IK)-ALP(JK)-ALP(IL)+ALP(JL)
         A(IJROT,KLROT) = A(IJROT,KLROT) + AFACT*VAL2
      END IF
C                                         (I,K//J,L) AND (I,L//J,K)
      IF(ISH.NE.KSH  .AND.  JSH.NE.LSH) THEN
         IKROT=IROT(I,K)
         JLROT=IROT(J,L)
         IF(IKROT.LT.JLROT) THEN
            NN    = IKROT
            IKROT = JLROT
            JLROT = NN
         END IF
         IF(J.GE.L) BFACT =  BET(IJ) -BET(JK) -BET(IL) +BET(KL)
         IF(J.LT.L) BFACT =  BET(IL) -BET(KL) -BET(IJ) +BET(JK)
         A(IKROT,JLROT) = A(IKROT,JLROT) + BFACT*VAL
      END IF
C                                         (I,L//K,J) AND (I,K//L,J)
      IF(ISH.NE.LSH  .AND.  JSH.NE.KSH) THEN
         ILROT=IROT(I,L)
         JKROT=IROT(J,K)
         IF(ILROT.LT.JKROT) THEN
            NN    = ILROT
            ILROT = JKROT
            JKROT = NN
         END IF
         IF(K.GE.J) BFACT =  BET(IK) -BET(KL) -BET(IJ) +BET(JL)
         IF(K.LT.J) BFACT =  BET(IJ) -BET(JL) -BET(IK) +BET(KL)
         A(ILROT,JKROT) = A(ILROT,JKROT) + BFACT*VAL
      END IF
C
C                  AND GO BACK FOR ANOTHER INTEGRAL BUFFER.
C
  800 CONTINUE
      IF(NX.GT.0) GO TO 300
C
C                  DONE WITH 2 ELECTRON INTEGRALS
C
  900 CONTINUE
      CALL SEQREW(IJKT)
C
C     IF RUNNING IN PARALLEL, WE MUST SUM ALL INTEGRAL CONTRIBUTIONS
C
      IF(GOPARR) THEN
         IF(TCSCF)  CALL DDI_GSUMF(2111,EPSIJ,3*L3)
                    CALL DDI_GSUMF(2112,EPS  ,L3)
                    CALL DDI_GSUMF(2113,ZETA ,L2*NHAM)
                    CALL DDI_GSUMF(2114,B    ,NROT*NXYZ)
         IF(INCORE) CALL GSUMCP(2116,A,NROT*NROT)
      END IF
C
C     DOUBLE THE DIAGONALS, SYMMETRIZE, AND SIGN CHANGE.
C
      II = 0
      DO 920 I=1,L1
         II = II + I
         DO 910 MSH=1,NHAM
            ZETA(II,MSH) = ZETA(II,MSH)+ZETA(II,MSH)
  910    CONTINUE
  920 CONTINUE
C
      IF(.NOT.INCORE) GO TO 980
      DO 950 IJROT=1,NROT
         A(IJROT,IJROT) = A(IJROT,IJROT)+A(IJROT,IJROT)
         DO 940 KLROT=1,IJROT
            A(KLROT,IJROT) = A(IJROT,KLROT)
  940    CONTINUE
  950 CONTINUE
C
  980 CONTINUE
      CALL DSCAL(L1*NORB*NXYZ,-ONE,W,1)
C
C        FOR TCSCF, ADD THE ONE ELECTRON TERM TO -EPSIJ-,
C                   AND FORM THE TOTAL LAGRANGIAN -EPS-.
C
      IF(TCSCF) THEN
         M = NORB-1
         N = NORB
         DO 1010 I=1,NORB
            FI = F(NCONF(I))
            IF(I.EQ.M  .OR.  I.EQ.N) FI = ONE
            DO 1000 J=1,L1
               IJ = IA(MAX(I,J)) + MIN(I,J)
               IF(I.NE.N) EPSIJ(I,J,1) = FI*HMO(IJ) + EPSIJ(I,J,1)
               IF(I.NE.M) EPSIJ(I,J,3) = FI*HMO(IJ) + EPSIJ(I,J,3)
 1000       CONTINUE
 1010    CONTINUE
C
         DO 1030 I=1,L1
            DO 1020 J=1,L1
               EPS(I,J) = CICOEF(1,1)*CICOEF(1,1)*EPSIJ(I,J,1) +
     *              TWO * CICOEF(1,1)*CICOEF(2,1)*EPSIJ(I,J,2) +
     *                    CICOEF(2,1)*CICOEF(2,1)*EPSIJ(I,J,3)
 1020       CONTINUE
 1030    CONTINUE
C
         IF(OUT) THEN
            WRITE(IW,*) 'CICOEF=', CICOEF(1,1),CICOEF(2,1)
            WRITE(IW,*) 'EPS(1,1) MATRIX'
            CALL PRSQ(EPSIJ(1,1,1),L1,L1,L1)
            WRITE(IW,*) 'EPS(1,2) MATRIX'
            CALL PRSQ(EPSIJ(1,1,2),L1,L1,L1)
            WRITE(IW,*) 'EPS(2,2) MATRIX'
            CALL PRSQ(EPSIJ(1,1,3),L1,L1,L1)
         END IF
C
C        FOR NON-TCSCF, SIMPLY ADD THE ONE ELECTRON TERM TO -EPS-
C
      ELSE
         IJ=0
         DO 1060 I=1,L1
            FI = ZERO
            IF(I.LE.NORB) FI = F(NCONF(I))
            DO 1050 J=1,L1
               IJ = IA(MAX(I,J)) + MIN(I,J)
               EPS(I,J) = FI*HMO(IJ) + EPS(I,J)
 1050       CONTINUE
 1060    CONTINUE
      END IF
C
C        ADD THE ONE ELECTRON CONTRIBUTIONS TO -ZETA-
C
      DO 1080 MSH=1,NHAM
         CALL DAXPY(L2,F(MSH),HMO,1,ZETA(1,MSH),1)
 1080 CONTINUE
C
C        THE VIRTUAL SHELL -ZETA- IS IDENTICALLY ZERO
C
      CALL VCLR(ZETA(1,NHAM+1),1,L2)
C
      IF(OUT) THEN
         WRITE(IW,*) 'LAGRANGIAN -EPS- IN MO BASIS IS'
         CALL PRSQ(EPS,L1,L1,L1)
         DO 1090 MSH=1,NHAM
            WRITE(IW,*) 'GEN.LAG. -ZETA- FOR SHELL',MSH
            CALL PRTRI(ZETA(1,MSH),L1)
 1090    CONTINUE
      END IF
C
      DSKWRK = .FALSE.
C
C     ----- LAGRANGIAN CONTRIBUTIONS TO -A- MATRIX -----
C
      IF(.NOT.INCORE) GO TO 1300
      DO 1250 I=1,L0
         ISH = NCONF(I)
         DO 1240 J=1,I
            JSH = NCONF(J)
            IF(ISH.EQ.JSH) GO TO 1240
            IJROT = IROT(I,J)
            DO 1230 K=1,L0
               KSH = NCONF(K)
               DO 1220 L=1,K
                  LSH = NCONF(L)
                  IF(KSH.EQ.LSH) GO TO 1220
                  KLROT = IROT(K,L)
C
                  IF(J.EQ.K) THEN
                     IL = IA(MAX(I,L)) + MIN(I,L)
                     A(IJROT,KLROT) = A(IJROT,KLROT) + EPS(I,L)
     *                                               - ZETA(IL,JSH)
                  END IF
C
                  IF(I.EQ.K) THEN
                     JL = IA(MAX(J,L)) + MIN(J,L)
                     A(IJROT,KLROT) = A(IJROT,KLROT) - EPS(J,L)
     *                                               + ZETA(JL,ISH)
                  END IF
C
                  IF(J.EQ.L) THEN
                     IK = IA(MAX(I,K)) + MIN(I,K)
                     A(IJROT,KLROT) = A(IJROT,KLROT) - EPS(I,K)
     *                                               + ZETA(IK,JSH)
                  END IF
C
                  IF(I.EQ.L) THEN
                     JK = IA(MAX(J,K)) + MIN(J,K)
                     A(IJROT,KLROT) = A(IJROT,KLROT) + EPS(J,K)
     *                                               - ZETA(JK,ISH)
                  END IF
 1220          CONTINUE
 1230       CONTINUE
 1240    CONTINUE
 1250 CONTINUE
C
 1300 CONTINUE
C
C        ----- ADD TCSCF CI DERIVATIVE TERMS -----
C        THE A11 AND B1-A BLOCKS ARE COMPLETE AT THIS POINT.
C
      IF(.NOT.TCSCF) GO TO 2000
      MROT = NROT-1
C
C                                 A22 BLOCK
C
      IF(INCORE) THEN
         A(MROT,MROT) = HALF*(CIHAM(1) -ETOT + CICOEF(1,1)*CICOEF(1,1))
         A(MROT,NROT) = HALF*(CIHAM(2)       + CICOEF(1,1)*CICOEF(2,1))
         A(NROT,MROT) = A(MROT,NROT)
         A(NROT,NROT) = HALF*(CIHAM(3) -ETOT + CICOEF(2,1)*CICOEF(2,1))
      ELSE
         ACI(MROT,1) = HALF*(CIHAM(1) -ETOT + CICOEF(1,1)*CICOEF(1,1))
         ACI(MROT,2) = HALF*(CIHAM(2)       + CICOEF(1,1)*CICOEF(2,1))
         ACI(NROT,1) = ACI(MROT,2)
         ACI(NROT,2) = HALF*(CIHAM(3) -ETOT + CICOEF(2,1)*CICOEF(2,1))
      END IF
C
C                                 A21 BLOCK
C
      DO 1520 I=1,L0
         DO 1510 J=1,I
            IJROT=IROT(I,J)
            IF(IJROT.EQ.0) GO TO 1510
C
            DUM = TWO*(CICOEF(1,1)*(EPSIJ(J,I,1)-EPSIJ(I,J,1))
     *               + CICOEF(2,1)*(EPSIJ(J,I,2)-EPSIJ(I,J,2)))
            IF(INCORE) THEN
               A(IJROT,MROT) = DUM
               A(MROT,IJROT) = DUM
            ELSE
               ACI(IJROT,1) = DUM
            END IF
C
            DUM = TWO*(CICOEF(1,1)*(EPSIJ(J,I,2)-EPSIJ(I,J,2))
     *               + CICOEF(2,1)*(EPSIJ(J,I,3)-EPSIJ(I,J,3)))
            IF(INCORE) THEN
               A(IJROT,NROT) = DUM
               A(NROT,IJROT) = DUM
            ELSE
               ACI(IJROT,2) = DUM
            END IF
 1510    CONTINUE
 1520 CONTINUE
C
C
 2000 CONTINUE
C
C      ---- IF CALCULATING ALPHA POLARIZABILITIES NEED TO SET UP ----
C      ---- 3 EXTRA COLUMNS OF -B- MATRIX                         ----
C
         CALL DAREAD(IDAF,IODA,HF(1,1),L2,252,0)
         CALL DAREAD(IDAF,IODA,HF(1,2),L2,253,0)
         CALL DAREAD(IDAF,IODA,HF(1,3),L2,254,0)
C
         DO 2100 I=1,L0
            ISH = NCONF(I)
            DO 2110 J=1,I
               JSH = NCONF(J)
               IJROT = IROT(I,J)
               IF(IJROT.EQ.0) GO TO 2100
               IF(ISH.EQ.JSH) GO TO 2110
               FIJ  = F(ISH) - F(JSH)
               IJ = IA(I) + J
                   B(IJROT,1)=FIJ*HF(IJ,1)
                   B(IJROT,2)=FIJ*HF(IJ,2)
                   B(IJROT,3)=FIJ*HF(IJ,3)
 2110       CONTINUE
 2100    CONTINUE
C
C     ---- ADD POSSIBLE CI CONTRIBUTIONS TO -B- MATRIX ----
C
         IF(TCSCF) THEN
            NCOPO=NORBB-2*NPAIR
            DUMX = ZERO
            DUMY = ZERO
            DUMZ = ZERO
            DO 2150 I=1,NCOPO
               II = IA(I) + I
               ISH = NCONF(I)
               DUMX = DUMX + F(ISH)*HF(II,1)
               DUMY = DUMY + F(ISH)*HF(II,2)
               DUMZ = DUMZ + F(ISH)*HF(II,3)
 2150       CONTINUE
C
            HMNX=ZERO
            HMNY=ZERO
            HMNZ=ZERO
            DO 2155 I=NORBB-1,NORBB
               ISH = NCONF(I)
               II = IA(I) + I
               HMNX=HMNX+F(ISH)*HF(II,1)
               HMNY=HMNY+F(ISH)*HF(II,2)
               HMNZ=HMNZ+F(ISH)*HF(II,3)
 2155       CONTINUE
C
            DO 2160 ICI=1,2
               IJROT=NROT-2+ICI
               NP=NCOPO+ICI
               ISH = NCONF(NP)
               II = IA(NP) + NP
               DUMAX = TWO*(DUMX + HF(II,1))
               DUMAY = TWO*(DUMY + HF(II,2))
               DUMAZ=  TWO*(DUMZ + HF(II,3))
               B(IJROT,1)=CICOEF(ICI,1)*(DUMX+HMNX)
     *                         -HALF*CICOEF(ICI,1)*DUMAX
               B(IJROT,2)=CICOEF(ICI,1)*(DUMY+HMNY)
     *                         -HALF*CICOEF(ICI,1)*DUMAY
               B(IJROT,3)=CICOEF(ICI,1)*(DUMZ+HMNZ)
     *                         -HALF*CICOEF(ICI,1)*DUMAZ
 2160       CONTINUE
         END IF
C
C     ----- REVERSE THE SIGN OF -B- FROM THAT IN THE PAPERS -----
C
      CALL DSCAL(NROT*NXYZ,-ONE,B,1)
C
      IF(OUT) THEN
         IF(INCORE) THEN
            WRITE(IW,*) 'THE CPHF -A- MATRIX IS'
            CALL PRSQ(A,NROT,NROT,NROT)
         END IF
         IF(.NOT.INCORE  .AND.  TCSCF) THEN
            WRITE(IW,*) 'THE CPHF -ACI- MATRIX IS'
            CALL PRSQ(ACI,2,NROT,NROT)
         END IF
         WRITE(IW,*) 'THE CPHF -B- MATRIX IS'
         CALL PRSQ(B,NXYZ,NROT,NROT)
      END IF
      RETURN
      END
C*MODULE LOCPOL    *DECK ROTPOL
      SUBROUTINE ROTPOL(DLPOL,VCOC,NLOC,IDLOC,MCORE)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,
     *        POLDCM,POLANG,POLAPP,KMIDPT,POLDYN
C
      DIMENSION DLPOL(9,NLOC),VCOC(NA*3)
      DIMENSION T(3,3),Q1(9),Q2(9),WRK(3)
C
      PARAMETER (TM6=1.0D-06)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (MXATM=2000, NMO=500)
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 /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      CALL DAREAD(IDAF,IODA,VCOC,NA*3,250,0)
C
      IJ=((IDLOC-1)*3) + (3*MCORE)
      IK=IDLOC+MCORE
      IF (NMOAT(IK).EQ.1) THEN
         DX=VCOC(IJ+1)-C(1,MOIDNO(1,IK))
         DY=VCOC(IJ+2)-C(2,MOIDNO(1,IK))
         DZ=VCOC(IJ+3)-C(3,MOIDNO(1,IK))
      ELSE
         DX=C(1,MOIDNO(1,IK))-C(1,MOIDNO(2,IK))
         DY=C(2,MOIDNO(1,IK))-C(2,MOIDNO(2,IK))
         DZ=C(3,MOIDNO(1,IK))-C(3,MOIDNO(2,IK))
      END IF
C
      DTOT = SQRT(DX**2+DY**2+DZ**2)
      DXY = SQRT(DX**2+DY**2)
C
      CALL VCLR(T,1,9)
C
      IF (DXY/TWO .LE.TM6) THEN
         T(1,1) = ONE
         T(2,2) = ONE
         T(3,3) = ONE
      ELSE
         T(1,1) =  DY/DXY
         T(2,1) = -DX/DXY
         T(3,1) =  ZERO
         T(1,2) =  DX*DZ/(DTOT*DXY)
         T(2,2) =  DY*DZ/(DTOT*DXY)
         T(3,2) = -DXY  /DTOT
         T(1,3) =  DX/DTOT
         T(2,3) =  DY/DTOT
         T(3,3) =  DZ/DTOT
      END IF
C
      Q1(1)=DLPOL(1,IDLOC)
      Q1(2)=DLPOL(2,IDLOC)
      Q1(3)=DLPOL(3,IDLOC)
      Q1(4)=DLPOL(4,IDLOC)
      Q1(5)=DLPOL(5,IDLOC)
      Q1(6)=DLPOL(6,IDLOC)
      Q1(7)=DLPOL(7,IDLOC)
      Q1(8)=DLPOL(8,IDLOC)
      Q1(9)=DLPOL(9,IDLOC)
C
      CALL TFSQU(Q2,Q1,T,WRK,3,3)
C
      DLPOL(1,IDLOC)=Q2(1)
      DLPOL(2,IDLOC)=Q2(2)
      DLPOL(3,IDLOC)=Q2(3)
      DLPOL(4,IDLOC)=Q2(4)
      DLPOL(5,IDLOC)=Q2(5)
      DLPOL(6,IDLOC)=Q2(6)
      DLPOL(7,IDLOC)=Q2(7)
      DLPOL(8,IDLOC)=Q2(8)
      DLPOL(9,IDLOC)=Q2(9)
C
      RETURN
      END
C*MODULE LOCPOL    *DECK POLAPM
      SUBROUTINE POLAPM(SOME,NXYZ,MCORE,NLOC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      NXYZ = 3
      NORB = NQMT
      NOCC = NA
      NVIR = NORB - NOCC
      NROT = NOCC * NVIR
C
      NUM2 =(NUM*NUM+NUM)/2
      NUM3 = NUM*NUM
C
      IF(NVIR .EQ. 0) RETURN
C
C     ---- MEMORY ALLOCATION FOR CALCULATION OF U AND ALPHA
C          POLARIZABILITY MATRICES ----
C
      CALL VALFM(LOADFM)
C
      LB   = 1+LOADFM
      LU   = LB   + NROT*NXYZ
      LH   = LU   + NROT*NXYZ
      LUL  = LH   + NUM2*3
      LHL  = LUL  + NVIR*3
      LE   = LHL  + NVIR*3
      LHP  = LE   + NUM
      LTR  = LHP  + NUM*NUM
      LPOL = LTR  + NLOC*NLOC
      LVCOC= LPOL + NLOC*9
      LC   = LVCOC+ NA*3
      LT   = LC   + NUM3
      LD   = LT   + NUM
      LAST = LD   + NUM2
      NEED = LAST - LOADFM - 1
C
      CALL GETFM(NEED)
C
      CALL POLAP(X(LU),X(LH),X(LUL),X(LHL),X(LE),X(LHP),X(LTR),
     *           X(LPOL),X(LVCOC),X(LC),X(LT),X(LD),SOME,
     *           NXYZ,MCORE,NLOC,NORB,NUM2,NROT,NOCC,NVIR)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE LOCPOL    *DECK POLAP
      SUBROUTINE POLAP(U,HF,UL,HFL,E,HFP,TRAN,DLPOL,VCOC,
     *                 CX,T,D,SOME,NXYZ,MCORE,NLOC,
     *                 NORB,NUM2,NROT,NOCC,NVIR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL SOME
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION U(NOCC,NVIR,NXYZ),UL(NVIR,3),HF(NUM2,3),HFL(NVIR,3),
     *          E(NUM),HFP(NUM,NUM),CX(NUM,NUM),T(NUM),D(NUM2)
C
      COMMON /IJPAIR/ IA(MXAO)
      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)
C
C  ---CALCULATE U CHEAPER, USE U=MUIJ/EI-EJ
C
      IF(SOME) WRITE(IW,9000)
      NXYZ = 3
      NORB = NQMT
      NOCC = NA
      NVIR = NORB - NOCC
      NROT = NOCC * NVIR
C
      NUM2 =(NUM*NUM+NUM)/2
      NUM3 = NUM*NUM
C
      IF(NVIR .EQ. 0) RETURN
C
C  TRANSFORM DIPOLE INTEGRALS TO MO BASIS AND WRITE FOR LATER USE
C
C    READ MO VECTORS
C
      CALL DAREAD(IDAF,IODA,CX,NUM3,15,0)
C
      CALL DAREAD(IDAF,IODA,D,NUM2,95,0)
      CALL TFTRI(HF(1,1),D,CX,T,NORB,NUM,NUM)
      CALL DAREAD(IDAF,IODA,D,NUM2,96,0)
      CALL TFTRI(HF(1,2),D,CX,T,NORB,NUM,NUM)
      CALL DAREAD(IDAF,IODA,D,NUM2,97,0)
      CALL TFTRI(HF(1,3),D,CX,T,NORB,NUM,NUM)
C
      CALL DAWRIT(IDAF,IODA,HF(1,1),NUM2,252,0)
      CALL DAWRIT(IDAF,IODA,HF(1,2),NUM2,253,0)
      CALL DAWRIT(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C   GET THE DIPOLE TENSORS, PUT INTO HF
C
C      CALL DAREAD(IDAF,IODA,HF(1,1),NUM2,252,0)
C      CALL DAREAD(IDAF,IODA,HF(1,2),NUM2,253,0)
C      CALL DAREAD(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C   GET RHF ORBITAL ENERGIES, PUT INTO E
C
      CALL DAREAD(IDAF,IODA,E,NUM,17,0)
C
C   CALCULATE U
C
      DO 300 IXYZ=1,NXYZ
         CALL CPYTSQ(HF(1,IXYZ),HFP,NUM,1)
         IL=0
         DO 320 IVIR=NA+1,NUM
            IL=IL+1
            DO 340 IOCC=1,NOCC
               U(IOCC,IL,IXYZ)=HFP(IOCC,IVIR)/(E(IOCC)-E(IVIR))
  340       CONTINUE
  320    CONTINUE
  300 CONTINUE
C
C     ----- CALCULATE LOCALIZED ALPHA POLARIZABILITIES ----
C
      NSTEP=0
      CALL LAPOL(U,UL,HF,HFL,TRAN,DLPOL,VCOC,
     *           IA,NXYZ,NUM2,NOCC,NVIR,MCORE,NLOC,NSTEP)
C
      RETURN
C
 9000 FORMAT(/5X,43(1H-)/
     *        5X,' CALCULATING U USING TRUNCATED PERTURBATION '/
     *        5X,43(1H-))
      END
C
C*MODULE LOCPOL    *DECK POLDB
      SUBROUTINE POLDB(DB,H2,HF,IA,NOCC,
     *                    NVIR,NUM2,
     *                    NROT,NWIDTH,NFTH2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DB(NROT,3),HF(NUM2,3),IA(*),H2(NROT,NWIDTH)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, EIGHT=8.0D+00)
C
C     ----- READ IN MO DIPOLE INTEGARLS -------
C
      CALL DAREAD(IDAF,IODA,HF(1,1),NUM2,252,0)
      CALL DAREAD(IDAF,IODA,HF(1,2),NUM2,253,0)
      CALL DAREAD(IDAF,IODA,HF(1,3),NUM2,254,0)
C
C     --- READ IN H2 MATRIX IN STRIPES AND FORM DYNAMIC B (DB=H2*HF) ---
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
C
      CALL SEQREW(NFTH2)
C
      KOL2 = 0
  100 CONTINUE
      KOL0 = KOL2
      KOL1 = KOL2 + 1
      KOL2 = KOL2 + NWIDTH
      IF(KOL2.GT.NROT) KOL2=NROT
            NCOL = KOL2-KOL1+1
C
            LEN = NROT*NCOL
            CALL SQREAD(NFTH2,H2,LEN)
C
           DO 150 K=1,NCOL
C
               SUM1 = ZERO
               SUM2 = ZERO
               SUM3 = ZERO
C
               IROT = 0
               DO 140 IVIR=1,NVIR
                   DO 130   IOCC=1,NOCC
                   IROT=IROT+1
                   IJ=IA(IVIR+NOCC)+IOCC
                   SUM1 = SUM1 + H2(IROT,K)*HF(IJ,1)
                   SUM2 = SUM2 + H2(IROT,K)*HF(IJ,2)
                   SUM3 = SUM3 + H2(IROT,K)*HF(IJ,3)
C
  130              CONTINUE
  140          CONTINUE
C
                   DB(K+KOL0,1) = EIGHT*SUM1
                   DB(K+KOL0,2) = EIGHT*SUM2
                   DB(K+KOL0,3) = EIGHT*SUM3
  150       CONTINUE
C
         IF(KOL2.LT.NROT) GO TO 100
C
      CALL SEQREW(NFTH2)
      DSKWRK = SVDSKW
      RETURN
      END
C
C*MODULE LOCPOL    *DECK POLH21FULL
      SUBROUTINE POLH21FULL(H1,H2,H21,
     *                    NROT,NWIDTH,NFT16,NFTH2,NFTH21)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW
C
      DIMENSION H2(NROT,NWIDTH), H1(NROT,NWIDTH),
     *          H21(NROT,NWIDTH)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C
C
C     ----- READ IN H1 MATRIX (OR A OR WXY) IN STRIPES   --------
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
C
      CALL SEQREW(NFTH21)
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
C
      KOL2 = 0
  100 CONTINUE
      KOL0 = KOL2
      KOL2 = KOL2 + NWIDTH
      IF(KOL2.GT.NROT) KOL2=NROT
      NCOL = KOL2 - KOL0
            LEN = NROT*NCOL
            CALL SQREAD(NFT16,H1,LEN)
            CALL SEQREW(NFTH2)
C
      LKOL2 = 0
  101 CONTINUE
      LKOL0 = LKOL2
      LKOL2 = LKOL2 + NWIDTH
      IF(LKOL2.GT.NROT) LKOL2=NROT
            LNCOL = LKOL2 - LKOL0
C
            LLEN = NROT*LNCOL
            CALL SQREAD(NFTH2,H2,LLEN)
      DO 160 K=1,NCOL
           DO 150 L=1,LNCOL
               SUM = 0.0D+00
C
C
               DO 140 IROT=1,NROT
C
                   SUM = SUM + H2(IROT,L)*H1(IROT,K)
C
  140          CONTINUE
                 H21(L+LKOL0,K) = SUM
  150       CONTINUE
  160    CONTINUE
C
C
         IF(LKOL2.LT.NROT) GO TO 101
C
      CALL SQWRIT(NFTH21,H21,NROT*NCOL)
C
         IF(KOL2.LT.NROT) GO TO 100
C
      CALL SEQREW(NFTH21)
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
      DSKWRK = SVDSKW
      RETURN
C
C 9050 FORMAT(1X,'SLICE OF -H21- N POLH21 FOR ROTATIONS',I5,' TO',I5)
C 9051 FORMAT(1X,'SLICE OF -H1- IN POLH21 FOR ROTATIONS',I5,' TO',I5)
      END
C
C
C*MODULE LOCPOL    *DECK POLH21
      SUBROUTINE POLH21(H1,H2,DIAG,
     *                  NROT,NWIDTH,NFT16,NFTH2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW
C
      DIMENSION H2(NROT,NWIDTH), H1(NROT,NWIDTH),
     *          DIAG(NROT)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C
C
C     ----- CALCULATE DIAGONAL OF H2*H1    --------
C     ASSUMES THAT H2 AND H1 ARE SYMMETRIC MATRICES
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
C
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
C
      KOL2 = 0
  100 CONTINUE
      KOL0 = KOL2
      KOL2 = KOL2 + NWIDTH
      IF(KOL2.GT.NROT) KOL2=NROT
      NCOL = KOL2 - KOL0
      LEN = NROT*NCOL
      CALL SQREAD(NFT16,H1,LEN)
      CALL SQREAD(NFTH2,H2,LEN)
C
      DO 160 K=1,NCOL
         DIAG(KOL0+K) = DDOT(NROT,H2(1,K),1,H1(1,K),1)
  160 CONTINUE
      IF(KOL2.LT.NROT) GO TO 100
C
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
      DSKWRK = SVDSKW
      RETURN
C
      END
C
C*MODULE LOCPOL    *DECK SOLVCPDYN
      SUBROUTINE SOLVCPDYN(H,DB,ZA,NROT,NXYZ,NOCP,
     *                  DPOL,NWIDTH,NSOLV,
     *                  NUM2,NOCC,NVIR,
     *                  NFT16,NFTH2,NFTH21,
     *                  ZAL,HL,TRAN,POL,VCOC,NLOC,MCORE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C
      LOGICAL MOIDON,EDCOMP,DIPDCM,QADDCM,DEPRNT,ZDO,
     *        POLDCM,POLANG,POLAPP,KMIDPT,POLDYN,POLAR
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW,SVGPAR
C
      PARAMETER (NMO=500, MXAO=8192, TWO=2.0D+00)
      PARAMETER (ZERO=0.0D+00, FOUR=4.0D+00,TOL=1.0D-10)
C
      COMMON /CHFINF/ CHFSLV,CPTYPE,NWDCHF,POLAR
      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 /IJPAIR/ IA(MXAO)
      COMMON /FMCOM/ X(1)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C
      DIMENSION DPOL(6),DB(NROT,NXYZ),ZA(NROT,NXYZ),
     *          NOCP(3),H(NUM2,3),
     *          POL(9,NLOC,NDPFREQ)
C
      DATA EXACT/8HGAUSS   /,  GMRES/8HGMRES   /
      DATA BCGST/8HBCGST   /, DIIS/8HDODIIS  /
C
      IF((CHFSLV .NE. EXACT) .AND. (CHFSLV .NE. BCGST) .AND.
     *   (CHFSLV .NE. GMRES) .AND. (CHFSLV .NE. DIIS)) THEN
         CHFSLV = GMRES
      END IF
C
C
      DO 420 IXYZ=1,NXYZ
         IF(NOCP(IXYZ).EQ.0) THEN
            DBMAX=ZERO
            DO 410 IROT=1,NROT
               DBI = -DB(IROT,IXYZ)
               IF(ABS(DBI).GT.DBMAX) DBMAX=ABS(DBI)
               DB(IROT,IXYZ) = DBI
                ZA(IROT,IXYZ) = DBI
  410       CONTINUE
            IF(DBMAX.LT.TOL) NOCP(IXYZ)=2
         END IF
         IF(NOCP(IXYZ).GT.0) CALL VCLR(ZA(1,IXYZ),1,NROT)
  420 CONTINUE
C
      IF(MASWRK) WRITE(IW,9040) CHFSLV,NROT
      IF(MASWRK) WRITE(IW,9020) NDPFREQ
      CALL FLSHBF(IW)
C
C
      CALL VALFM(LOADFM)
      LAST = LOADFM
C
C     ------------CALCULATION OF---
C     -----DYNAMIC POLARIZABILITIES---------------
C
C   GAUSS-JORDAN ELIMINATION
      IF(CHFSLV .EQ. EXACT) THEN
C
C     ------ ALLOCATE MEMORY FOR H2*H1
C
         IF(MASWRK) WRITE(IW,*) 'MEMORY ALLOCATION FOR H2*H1'
         CALL FLSHBF(IW)
         LH1    = 1+LOADFM
         LH2    = LH1    + NROT*NWIDTH
         LH21   = LH2    + NROT*NWIDTH
         LAST   = LH21   + NROT*NWIDTH
         NEED = LAST - LOADFM -1
         CALL GETFM(NEED)
C
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQOPN(NFTH21,'DCPHF21','UNKNOWN',.FALSE.,'UNFORMATTED')
         DSKWRK = SVDSKW
C
         SVGPAR = GOPARR
         SVDSKW = DSKWRK
         GOPARR = .FALSE.
         DSKWRK = .FALSE.
         IF(MASWRK) THEN
            CALL POLH21FULL(X(LH1),X(LH2),X(LH21),
     *                    NROT,NWIDTH,NFT16,NFTH2,NFTH21)
         END IF
         GOPARR = SVGPAR
         DSKWRK = SVDSKW
         CALL RETFM(NEED)
C
C    ----- PREPARING MEMORY FOR EXACT SOLVER
         NSOLV=1
         LH21IM  = 1 + LOADFM
         LWRK1 = LH21IM + NROT*NROT
         LWRK2 = LWRK1 + NROT*NSOLV
         LAST  = LWRK2 + NROT*NSOLV
         NEEDEX = LAST - LOADFM
         CALL GETFM(NEEDEX)
C
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQREW(NFTH21)
         LEN21 = NROT*NROT
         CALL SQREAD(NFTH21,X(LH21IM),LEN21)
         CALL SEQREW(NFTH21)
C
C
      ELSE
C
C     BCGSTAB
         IF(CHFSLV .EQ. BCGST) THEN
C
           IF(MASWRK) WRITE(IW,*)'ALLOCATING MEMORY FOR BCGSTAB SOLVER'
            CALL FLSHBF(IW)
C    MEMORY ALLOCATION FOR BCGSTAB SOLVER
            LD    = LOADFM + 1
            LP    = LD   + NROT
            LV    = LP  + NROT
            LS    = LV   + NROT
            LT    = LS  + NROT
            LR0   = LT  + NROT
            LR    = LR0 + NROT
            LRR   = LR  + NROT
            LAST  = LRR + NROT
            NEEDBCG  = LAST - LOADFM - 1
            CALL GETFM(NEEDBCG)
            IF(MASWRK) WRITE(IW,*)'ALLOCATED'
            CALL FLSHBF(IW)
C
C
         ELSE IF (CHFSLV .EQ. GMRES) THEN
C
            IF(MASWRK) WRITE(IW,*) 'ALLOCATING MEMORY FOR GMRES SOLVER'
            CALL FLSHBF(IW)
C
C    MEMORY ALLOCATION FOR GMRES SOLVER
            LD    = LOADFM + 1
            LV    = LD + NROT
            LWRK1 = LV  + NROT
            LVV   = LWRK1  + NROT
            LAST    = LVV  + NROT*51
            NEEDGM  = LAST - LOADFM - 1
            CALL GETFM(NEEDGM)
C
C   DIIS SOLVER
         ELSE IF (CHFSLV .EQ. DIIS) THEN
C
            IF(MASWRK) WRITE(IW,*) 'ALLOCATING MEMORY FOR DIIS SOLVER'
            CALL FLSHBF(IW)
            MXDII = 10
C
            LQOUT  = LOADFM +1
            LQIN   = LQOUT   + NROT
            LM1Z    = LQIN    + NROT
            LDIMAT = LM1Z    + NROT
            LQREP  = LDIMAT  + (MXDII+1)*(MXDII+1)
            LTMP1  = LQREP   + 2*NROT*MXDII
            LTMP2  = LTMP1   + MXDII + 1
            LIPVT  = LTMP2   + (MXDII+1)*(MXDII+1)
            LAST   = LIPVT   + MXDII + 1
            NEEDDIIS   = LAST    - LOADFM - 1
            CALL GETFM(NEEDDIIS)
C
         END IF
C
C  ----  ALLOACTE MEMORY FOR H1 AND H2
C  ---- AND COMPUTE THE DIAGONAL ELEMENTS OF H21, STORE IN LDIAG
C
         LH1   = 1+ LAST
         LH2   = LH1   + NROT*NWIDTH
         LDIAG = LH2 + NROT*NWIDTH
         LASTD  = LDIAG   + NROT
         NEEDDG  = LASTD - LAST - 1
         CALL GETFM(NEEDDG)
C
         IF(MASWRK) WRITE(IW,*) 'CALCULATING H2*H1 DIAGONAL'
         CALL FLSHBF(IW)
         IF(MASWRK) WRITE(IW,*) 'WIDTH=',NWIDTH,'NROT=',NROT
         CALL FLSHBF(IW)
         SVDSKW = DSKWRK
         SVGPAR = GOPARR
         GOPARR = .FALSE.
         DSKWRK = .FALSE.
         IF(MASWRK) THEN
            CALL POLH21(X(LH1),X(LH2),X(LDIAG),
     *               NROT,NWIDTH,NFT16,NFTH2)
         END IF
         GOPARR = SVGPAR
         DSKWRK = SVDSKW
C        CALL RETFM(NEEDDG)
C
      END IF
C
C
C    ------- LOOP 110 OVER FREQUENCIES
      DO 110 I=1,NDPFREQ
      IF(MASWRK) WRITE(IW,9030) DPFREQ(I)
         CALL FLSHBF(IW)
         FREQI=DPFREQ(I)
         FREQI2=FREQI*FREQI
C
C
C      ------ BCGSTAB
         IF(CHFSLV .EQ. BCGST) THEN
C
            CALL FLSHBF(IW)
C            TOLREL = 1.0D-06
C            TOLABS = 1.0D-06
            NITER  = 100
C
C           IF(MASWRK) WRITE(IW,*)'CALLING BCGSTAB SOLVER'
C           CALL FLSHBF(IW)
            DO 440 IXYZ=1,NXYZ
               IF(NOCP(IXYZ).GT.0) GO TO 440
               CALL BCGSTB(X(LH1),X(LH2),ZA(1,IXYZ),DB(1,IXYZ),
     *               X(LDIAG),X(LD),NROT,NITER,
     *               NWIDTH,X(LP),X(LV),X(LS),X(LT),X(LR0),X(LR),
     *               X(LRR),
     *               NFT16,NFTH2,.TRUE.,FREQI)
 440        CONTINUE
C            CALL RETFM(NEEDBCG)
C
C     --------- GMRES
         ELSE IF (CHFSLV .EQ. GMRES) THEN
C
            MAXITS = 100
            IM = 45
C
C           IF(MASWRK) WRITE(IW,*) 'CALLING GMRES SOLVER'
C           CALL FLSHBF(IW)
C
            DO 450 IXYZ=1,NXYZ
               IF(NOCP(IXYZ).GT.0) GO TO 450
               CALL PGMRES(X(LH1),X(LH2),DB(1,IXYZ),ZA(1,IXYZ),NROT,IM,
     *                 X(LV),X(LVV),X(LWRK1),
     *                 MAXITS,X(LDIAG),X(LD),NWIDTH, NFT16, NFTH2,
     *                 IERR,.TRUE., FREQI)
  450       CONTINUE
C         CALL RETFM(NEEDGM)
C
C    --------- DIIS
         ELSE IF (CHFSLV .EQ. DIIS) THEN
C
            MXDII = 10
            FREQI2=FOUR*DPFREQ(I)*DPFREQ(I)
            DO 460 IXYZ=1,NXYZ
               IF(NOCP(IXYZ).GT.0) GO TO 460
               CALL SEMIIT(X(LH1),X(LH2),FREQI2,NWIDTH,NFT16,NFTH2,
     *               ZA(1,IXYZ),NROT,X(LQOUT),X(LQIN),
     *               X(LDIAG),X(LM1Z),X(LDIMAT),X(LQREP),
     *               X(LTMP1),X(LTMP2),X(LIPVT),MXDII)
               DO J = 1, NROT
                  ZA(J,IXYZ) = X(LQOUT + J -1)
               END DO
  460       CONTINUE
C        CALL RETFM(NEEDDIIS)
C
C     ---------- EXACT
         ELSE IF (CHFSLV .EQ. EXACT) THEN
C
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
         CALL SEQREW(NFTH21)
         LEN21 = NROT*NROT
         CALL SQREAD(NFTH21,X(LH21IM),LEN21)
         CALL SEQREW(NFTH21)
C
         CALL ADDFR(X(LH21IM),NROT,FREQI,.TRUE.)
C
C        IF(DBG) THEN
C            WRITE(IW,*) 'INITIAL ZA VECTOR'
C            CALL PRSQ(ZA,NXYZ,NROT,NROT)
C        END IF
C
            NODCMP=1
            DO 475 IXYZ=1,NXYZ
               IF(NOCP(IXYZ).GT.0) GO TO 475
               CALL LINEQU(X(LH21IM),NROT,ZA(1,IXYZ),NROT,
     *                    X(LWRK1),X(LWRK2),DETWXY,IERR,NODCMP)
               NODCMP=0
  475       CONTINUE
C            CALL ADDFR(X(LH21IM),NROT,FREQI,.FALSE.)
C
         END IF
C
C   ----------- ADD POLARIZABILITY
C
C     ---- GET TRANSFORMED DIPOLE INTEGRALS AND CALCULATE RESPONSE ----
C
C
         CALL DAREAD(IDAF,IODA,H(1,1),NUM2,252,0)
         CALL DAREAD(IDAF,IODA,H(1,2),NUM2,253,0)
         CALL DAREAD(IDAF,IODA,H(1,3),NUM2,254,0)
C
C     ---- CALCULATE ALPHA POLARIZABILITY TENSOR ----
C
         DO 402 K=1,3
            DO 412 L=1,K
               KL=IA(K)+L
               DPOL(KL)=ZERO
               IROT=0
               DO 422 IVIR=1,NVIR
                  DO 432 IOCC=1,NOCC
                     IROT=IROT + 1
                     IJ=IA(IVIR+NOCC)+IOCC
                     DPOL(KL)=DPOL(KL)-(TWO*ZA(IROT,K)*H(IJ,L))
  432             CONTINUE
  422          CONTINUE
  412       CONTINUE
  402    CONTINUE
C
C
C    -----  CALCULATE LMO CONTRIBUTION TO DYNAMIC POLARIZABILITY
C
         NSTEP=I
         CALL LAPOL(ZA,ZAL,H,HL,TRAN,POL(1,1,I),VCOC,
     *              IA,NXYZ,NUM2,NOCC,NVIR,MCORE,NLOC,NSTEP)
C
C
         DO 600 IXYZ=1,NXYZ
            IF(NOCP(IXYZ).EQ.2) NOCP(IXYZ)=0
  600    CONTINUE
C
C  ----- RESTORE OLD VALUE FOR DB  ---------
C       (INITIAL ZA HAS DB MATRIX)
C
         IF (CHFSLV .EQ. EXACT .OR. CHFSLV .EQ. DIIS) THEN
            DO 421 IXYZ=1,NXYZ
               DO 411 IROT=1,NROT
                  ZA(IROT,IXYZ) = DB(IROT,IXYZ)
  411          CONTINUE
  421       CONTINUE
         END IF
C
  110 CONTINUE
C
C
C      ------ RETURNING MRMORY -----
C
      IF (CHFSLV .EQ. EXACT) THEN
         CALL RETFM(NEEDEX)
      ELSE
         CALL RETFM(NEEDDG)
         IF (CHFSLV .EQ. BCGST) THEN
            CALL RETFM(NEEDBCG)
         ELSE IF (CHFSLV .EQ. GMRES) THEN
            CALL RETFM(NEEDGM)
         ELSE IF (CHFSLV .EQ. DIIS) THEN
            CALL RETFM(NEEDDIIS)
         END IF
      END IF
C
      RETURN
C
 9020 FORMAT(/1X,'A TOTAL OF',I4,' IMAGINARY FREQUENCY',
     *           ' POLARIZABILITIES ARE TO BE FOUND')
 9030 FORMAT(/1X,'COMPUTING DYNAMIC POLARIZABILTY FOR IMAGINARY',
     *           ' FREQUENCY',F12.7,'I (A.U.)')
 9040 FORMAT(1X,'USING ',A10,' SOLVER FOR LINEAR SYSTEM',
     *          ' OF SIZE',I10)
C
      END
C
C
C*MODULE LOCPOL    *DECK BCGSTB
      SUBROUTINE BCGSTB(H1,H2,X,RHS,ADIAG,DIAG,N,NITER,NWIDTH,
     *                  P,V,S,T,R0,R,RR,
     *                  NFT16,NFTH2,OMEGA,FREQI)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C
      LOGICAL DSKWRK,MASWRK,SVDSKW,GOPARR,OMEGA
      LOGICAL BRKDN
      EXTERNAL BRKDN
C
      DIMENSION H1(N,NWIDTH),H2(N,NWIDTH),X(N),RHS(N),
     *          ADIAG(N),DIAG(N),S(N),V(N),P(N),R(N),
     *          R0(N),T(N),RR(N)
C
C     ----- SOLVES A*X=Y BY USING ITERATIVE PRECONDITIONED ---------
C     STABILIZIED BICONLUGATE GRADIENT SOLVER (BCGSTAB).
C     MODIFIED FROM SPARSKIT2 SUBROUTINE BCGSTAB (AUTHOR YOUCEF SAAD)
C                 BY LYUDMILA SLIPCHENKO, 10/2005
C
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (FOUR=4.0D+00,EPS=1.0D-05)
C
C     ----- SOLVE (H2*H1 - FREQI^2)X = RHS -----
C     USING AN ITERATIVE BICONJUGATE GRADIENT STABILIZED ALGORITHM
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
C
      EPS2=N*EPS*EPS
C
C-------------- COMPUTE INITIAL RESIDUAL VECTOR --------------
C
C
C--H1*X_0
       CALL CALCAXS (H1,X,R,.FALSE.,N,NWIDTH,
     *             .TRUE.,NFT16,NFTH2)
C--H2*(H1*X_0)
       CALL CALCAXS (H2,R,RR,.FALSE.,N,NWIDTH,
     *             .FALSE.,NFT16,NFTH2)
C
       DO 130 J=1,N
          IF(OMEGA) THEN
             RR(J) = RR(J)+ FOUR*FREQI*FREQI*X(J)
             DIAG(J) = ADIAG(J) + FOUR*FREQI*FREQI
          END IF
          DIAG(J) = 1.0D+00 / DIAG(J)
C           DIAG(J) = 1.0D+00
C
          RR(J) = RHS(J)-RR(J)
C     PRECONDITIONING
          RR(J) = RR(J)*DIAG(J)
          R0(J) = RR(J)
          P(J) = RR(J)
C   FOR STOPPING CRITERIA
          R(J) = RHS(J)*DIAG(J)
  130  CONTINUE
C
C         WRITE(IW,*) 'DIAGONAL'
C         CALL PRSQ(DIAG,1,N,N)
C
C         WRITE(IW,*) 'RR'
C         CALL PRSQ(RR,1,N,N)
C
       RHO0 = DDOT(N, R0,1, RR,1)
       RHO=RHO0
       RESDU = SQRT(RHO)
C       RESDT = TOLREL*SQRT(DDOT(N, R,1,R,1)) + TOLABS
C
C       CONVERGENCE CRITERIA: ||CUR.RES.|| < THRESH*||RHS||
       RESDT = EPS2*SQRT(DDOT(N, R,1,R,1))
C         RESDT = TOLREL*RESD0 + TOLABS
       IF(RESDU.LE.RESDT) THEN
          WRITE(IW,*) 'INITIAL GUESS IS NOT GOOD'
          CALL ABRT
          STOP
       END IF
C
C  PRINT INITIAL VALUES
C
C       IF(MASWRK) WRITE(IW,*)
C     *       'RHO0 = ',RHO0,' STARTING RESIDUAL = ',RESD0,
C     *       ' TARGET RESIDUAL = ',RESDT
C       CALL FLSHBF(IW)
C
C   BEGINNING OF THE ITERATIONS
C
C
      DO 500 ITER=1,NITER
C
C     STEP (1), V = A*P
C
       CALL CALCAXS (H1,P,R,.FALSE.,N,NWIDTH,.TRUE.,NFT16,NFTH2)
       CALL CALCAXS (H2,R,V,.FALSE.,N,NWIDTH,.FALSE.,NFT16,NFTH2)
C
C  ADD IM. FREQUENCIES AND PRECONDITION
C
       DO 143 J=1,N
          IF(OMEGA)  V(J) = V(J)+FOUR*FREQI*FREQI*P(J)
           V(J) = V(J)*DIAG(J)
  143  CONTINUE
C
C
C    STEP (2)
C
       ALPHA = DDOT(N,R0,1,V,1)
       IF(BRKDN(ALPHA)) THEN
          WRITE(IW,*) 'ALPHA IS BAD'
          CALL ABRT
          STOP
       END IF
       ALPHA = RHO/ALPHA
C
C     STEP (3)
       DO 150 J=1,N
          S(J) = RR(J) - ALPHA*V(J)
  150  CONTINUE
C
C   STEP (4): THE SECOND MATVEC -- T = A*S
C
       CALL CALCAXS (H1,S,R,.FALSE.,N,NWIDTH,.TRUE.,NFT16,NFTH2)
       CALL CALCAXS (H2,R,T,.FALSE.,N,NWIDTH,.FALSE.,NFT16,NFTH2)
C
C   ADD IM. FREQUENCIES AND PRECONDITION
C
       DO 144 J=1,N
          IF(OMEGA)  T(J) = T(J)+FOUR*FREQI*FREQI*S(J)
          T(J) = T(J)*DIAG(J)
  144  CONTINUE
C
C     STEP (5)
       OMEG = DDOT(N,T,1,T,1)
       IF(BRKDN(OMEG)) THEN
          WRITE(IW,*) 'OMEG IS BAD'
          CALL ABRT
          STOP
       END IF
       OMEG = DDOT(N,T,1,S,1)/OMEG
       IF(BRKDN(OMEG)) THEN
          WRITE(IW,*) 'OMEG IS BAD'
          CALL ABRT
          STOP
       END IF
C
C     STEP (6) AND (7)
       DO 170 J=1,N
           R(J) = ALPHA*P(J) + OMEG*S(J)
           X(J) = X(J) + R(J)
           RR(J) = S(J) - OMEG*T(J)
  170  CONTINUE
C
       BETA = RHO
C
C   CONVERGENCE TEST
C
       RHO = DDOT(N,RR,1,RR,1)
       IF(BRKDN(RHO)) THEN
           WRITE(IW,*) 'RHO IS BAD'
           CALL ABRT
           STOP
       END IF
C
       RESDU = SQRT(RHO)
C
C  PRINT RESULT FOR THIS ITERATION
C
C       IF(MASWRK) WRITE(IW,9020) ITER, RESDU
C       CALL FLSHBF(IW)
C
       IF(RESDU.LE.RESDT) THEN
          IF(MASWRK)
     *       WRITE(IW,*) 'BICGTSTAB CONVERGED AFTER ',ITER,
     *       ' ITERATIONS'
          RETURN
       END IF
C
C     STEP (8): COMPUTING NEW P
C
       BETA = RHO*ALPHA/(OMEG*BETA)
C
       DO 180 J=1,N
          P(J) = RR(J)+BETA* (P(J) - OMEG*V(J))
  180  CONTINUE
C
C  END OF ITERATION
C
 500   CONTINUE
C
C
       IF(MASWRK) THEN
          WRITE(IW,*) ' '
          WRITE(IW,*) 'BCGSTAB ITERATIVE PROCEDURE DID NOT CONVERGED
     *                 IN', NITER,' ITERATIONS'
          WRITE(IW,*) 'GAMESS WILL PROCEED DESPITE THIS ...'
          WRITE(IW,*) ' '
       END IF
C
      DSKWRK = SVDSKW
       RETURN
C
C 9020  FORMAT(1X,'ITER',I4,' , MAX RESIDUE=',
C     *       1P,E15.7)
       END
C
C
C*MODULE LOCPOL    *DECK PGMRES
       SUBROUTINE PGMRES(H1, H2, RHS, SOL, N, IM, V, VV, WRK1,
     *                   MAXITS,
     *                   DIAG, DI, NWIDTH, NFT16, NFTH2, IERR,
     *                   OMEGA, FREQ)
C-----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DSKWRK,MASWRK,OMEGA,GOPARR
C
      PARAMETER (FOUR=4.0D+00,EPS=1.0D-06)
      PARAMETER (KMAX = 50)
C
      DIMENSION H1(N,NWIDTH),H2(N,NWIDTH),SOL(N),RHS(N),
     *          DIAG(N),DI(N),V(N),VV(N,*),WRK1(N),HH(KMAX+1,KMAX),
     *          C(KMAX),S(KMAX), RS(KMAX+1)
C
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C-------------------------------------------------------------
C ARNOLDI SIZE SHOULD NOT EXCEED KMAX=50 IN THIS VERSION.
C TO RESET MODIFY PARAMTER KMAX ACCORDINGLY.
C-------------------------------------------------------------
C
C     ----- SOLVE (H2*H1 + FREQI^2)X = RHS -----
C     USING AN ITERATIVE PRECONDITIONED GMRES ALGORITHM
C     MODIFIED FROM SPARSKIT PACKAGE BY LYUDMILA SLIPCHENKO. 09/20/2005
C
C----------------------------------------------------------------------*
C PARAMETERS
C-----------
C ON ENTRY:
C==========
C
C N     == INTEGER. THE DIMENSION OF THE MATRIX.
C IM    == SIZE OF KRYLOV SUBSPACE:  SHOULD NOT EXCEED 50 IN THIS
C          VERSION (CAN BE RESET BY CHANGING PARAMETER COMMAND FOR
C          KMAX BELOW)
C RHS   == REAL VECTOR OF LENGTH N CONTAINING THE RIGHT HAND SIDE.
C SOL   == REAL VECTOR OF LENGTH N CONTAINING AN INITIAL GUESS TO THE
C          SOLUTION ON INPUT. APPROXIMATE SOLUTION ON OUTPUT
C EPS   == TOLERANCE FOR STOPPING CRITERION. PROCESS IS STOPPED
C          AS SOON AS ( ||.|| IS THE EUCLIDEAN NORM):
C          || CURRENT RESIDUAL||/||RHS|| <= EPS
C MAXITS== MAXIMUM NUMBER OF ITERATIONS ALLOWED
C H1    == MATRIX H1
C H2    == MATRIX H2
C DIAG  == DIAGONAL PART OF H1H2 MATRIX
C NWIDTH== PART OF H1 OR H2 IN MEMORY
C NFT16 == STORAGE OF H1 ON DISK
C NFTH2 == STORAGE OF H2 ON DISK
C OMEGA
C FREQ
C
C                                                    *
C ON RETURN:
C==========
C SOL   == CONTAINS AN APPROXIMATE SOLUTION (UPON SUCCESSFUL RETURN).
C IERR  == INTEGER. ERROR MESSAGE WITH THE FOLLOWING MEANING.
C          IERR = 0 --> SUCCESSFUL RETURN.
C          IERR = 1 --> CONVERGENCE NOT ACHIEVED IN ITMAX ITERATIONS.
C          IERR =-1 --> THE INITIAL GUESS SEEMS TO BE THE EXACT
C                       SOLUTION (INITIAL RESIDUAL COMPUTED WAS ZERO)
C
C----------------------------------------------------------------------*
C
C WORK ARRAYS:
C=============
C VV    == WORK ARRAY OF LENGTH  N X (IM+1) (USED TO STORE THE ARNOLI
C          BASIS)
C WRK1  == REAL VECTOR OF LENGTH N
C DI    == RECIOPROCAL OF DIAGONAL PART OF (H1H2 + 4*FREQ^2) MATRIX
C----------------------------------------------------------------------*
       DATA EPSMAC/1.D-16/
C
      EPS2=N*EPS*EPS
      ITS = 0
      C(1)=0.0D+00
      S(1)=0.0D+00
C
      CALL SEQREW(NFT16)
      CALL SEQREW(NFTH2)
C
C      COMPUTE INITIAL RESIDUAL VECTOR
C  H1*X_0
       CALL CALCAXS (H1,SOL,V,.FALSE.,N,NWIDTH,.TRUE.,NFT16,NFTH2)
C  H2*(H1*X_0)
       CALL CALCAXS (H2,V,VV,.FALSE.,N,NWIDTH,.FALSE.,NFT16,NFTH2)
C
C  ADD FREQUENCIES IN HERE
         DO 130 J=1,N
            IF(OMEGA)
     *         VV(J,1) = VV(J,1)+ FOUR*FREQ*FREQ*SOL(J)
            VV(J,1) = RHS(J)-VV(J,1)
  130    CONTINUE
C
C
C  PREPARE CORRECT DIAGONAL
         DO 132 J=1,N
            DI(J) = 1.0D+00/(DIAG(J)+FOUR*FREQ*FREQ)
  132    CONTINUE
C
         BNORM = DNRM2(N, RHS, 1)
          EPS1=EPS2*BNORM
C
 20      RO = DNRM2(N, VV, 1)
C         IF (MASWRK .AND. ITS .EQ. 0)
C     *           WRITE(IW, 199) ITS, RO
C         CALL FLSHBF(IW)
         IF (RO .EQ. 0.0D+00) GOTO 999
         T = 1.0D+00/ RO
         DO 210 J=1, N
            VV(J,1) = VV(J,1)*T
 210     CONTINUE
C
C         IF (ITS .EQ. 0) EPS1=EPS*RO
C
C   INITIALIZE 1-ST TERM  OF RHS OF HESSENBERG SYSTEM..
        RS(1) = RO
        I = 0
 4      I=I+1
        ITS = ITS + 1
        I1 = I + 1
C
C    PRECONDITIONING
C
        DO 220 J=1, N
            WRK1(J) = VV(J,I)*DI(J)
 220    CONTINUE
C
C     H1*X_0
       CALL CALCAXS(H1,WRK1,V,.FALSE.,N,NWIDTH,.TRUE.,NFT16,NFTH2)
C     H2*(H1*X_0)
       CALL CALCAXS(H2,V,VV(1,I1),.FALSE.,N,NWIDTH,.FALSE.,NFT16,NFTH2)
C
C    ADD FREQUENCIES IN HERE
       DO 135 J=1,N
          IF(OMEGA)
     *        VV(J,I1) = VV(J,I1)+ FOUR*FREQ*FREQ*WRK1(J)
  135    CONTINUE
C
C
C-----------------------------------------
C     MODIFIED GRAM - SCHMIDT...
C-----------------------------------------
C       IF(MASWRK) WRITE(IW,*) 'GRAM-SCHMIDT START'
C       CALL FLSHBF(IW)
      DO 55 J=1, I
          T = DDOT(N, VV(1,J),1,VV(1,I1),1)
          HH(J,I) = T
          CALL DAXPY(N, -T, VV(1,J), 1, VV(1,I1), 1)
 55    CONTINUE
       T = DNRM2(N, VV(1,I1), 1)
       HH(I1,I) = T
       IF ( T .EQ. 0.0D+00) GOTO 58
       T = 1.0D+00/T
       DO 57  K=1,N
          VV(K,I1) = VV(K,I1)*T
 57    CONTINUE
C     IF(MASWRK) WRITE(IW,*) 'GRAM-SCHMIDT DONE'
C      CALL FLSHBF(IW)
C
C     DONE WITH MODIFIED GRAM SCHIMD AND ARNOLDI STEP..
C     NOW  UPDATE FACTORIZATION OF HH
C
 58    IF (I .EQ. 1) GOTO 121
C
C   PERFROM PREVIOUS TRANSFORMATIONS  ON I-TH COLUMN OF H
       DO 66 K=2,I
          K1 = K-1
          T = HH(K1,I)
          HH(K1,I) = C(K1)*T + S(K1)*HH(K,I)
          HH(K,I) = -S(K1)*T + C(K1)*HH(K,I)
 66    CONTINUE
 121   GAM = SQRT(HH(I,I)**2 + HH(I1,I)**2)
C
C     IF GAMMA IS ZERO THEN ANY SMALL VALUE WILL DO...
C     WILL AFFECT ONLY RESIDUAL ESTIMATE
C
       IF (GAM .EQ. 0.0D+00) GAM = EPSMAC
C
C     GET  NEXT PLANE ROTATION C
       C(I) = HH(I,I)/GAM
       S(I) = HH(I1,I)/GAM
       RS(I1) = -S(I)*RS(I)
       RS(I) =  C(I)*RS(I)
C
C     DETRERMINE RESIDUAL NORM AND TEST FOR CONVERGENCE- C
       HH(I,I) = C(I)*HH(I,I) + S(I)*HH(I1,I)
       RO = ABS(RS(I1))
C
C  UNCOMMENT FOR ITERATION RESULTS
C       IF(MASWRK) WRITE(IW,199) ITS, RO
C       CALL FLSHBF(IW)
C
       IF (I .LT. IM .AND. (RO .GT. EPS1))  GOTO 4
C
C     NOW COMPUTE SOLUTION. FIRST SOLVE UPPER TRIANGULAR SYSTEM.
C
       RS(I) = RS(I)/HH(I,I)
       DO 30 II=2,I
          K=I-II+1
          K1 = K+1
          T=RS(K)
          DO 40 J=K1,I
             T = T-HH(K,J)*RS(J)
 40       CONTINUE
          RS(K) = T/HH(K,K)
 30    CONTINUE
C
C     FORM LINEAR COMBINATION OF V(*,I)'S TO GET SOLUTION C
       T = RS(1)
       DO 15 K=1, N
          WRK1(K) = VV(K,1)*T
 15    CONTINUE
       DO 16 J=2, I
          T = RS(J)
          DO 161 K=1, N
             WRK1(K) = WRK1(K)+T*VV(K,J)
 161      CONTINUE
 16    CONTINUE
C
C     PRECONDITIONING
C
       DO 17 K=1, N
          WRK1(K) = WRK1(K)*DI(K)
          SOL(K) = SOL(K) + WRK1(K)
 17    CONTINUE
C
C     RESTART OUTER LOOP  WHEN NECESSARY
C
       IF (RO .LE. EPS1) GOTO 990
       IF (ITS .GE. MAXITS) GOTO 991
C
C     ELSE COMPUTE RESIDUAL VECTOR AND CONTINUE..
C
       DO 24 J=1,I
          JJ = I1-J+1
          RS(JJ-1) = -S(JJ-1)*RS(JJ)
          RS(JJ) = C(JJ-1)*RS(JJ)
 24    CONTINUE
       DO 25  J=1,I1
          T = RS(J)
          IF (J .EQ. 1)  T = T-1.0D+00
          CALL DAXPY (N, T, VV(1,J), 1,  VV, 1)
 25    CONTINUE
C
C 199   FORMAT('   ITS =', I4, ' RES. NORM =', D20.6)
C
C     RESTART OUTER LOOP.
       GOTO 20
C
 990   IERR = 0
       IF(MASWRK) WRITE(IW,*)'GMRES CONVERGED IN ',ITS,'ITERATIONS'
       CALL FLSHBF(IW)
       RETURN
 991   IERR = 1
       RETURN
 999   CONTINUE
       IERR = -1
       RETURN
C
C-----------------END OF PGMRES
       END
C
C*MODULE LOCPOL    *DECK CALCAXS
      SUBROUTINE CALCAXS(A,X,AX,INCORE,N,NWIDTH,IFH1,NFT16,NFTH2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL INCORE,IFH1
C
      DIMENSION A(N,NWIDTH),X(N),AX(N)
C
C     ----- THIS ROUTINE COMPUTES AX = A * X -----
C     SIMPLE A*X MULTIPLICATION, ONE VECTOR/ONE SOLUTION IN A TIME
C     A IS SYMMETRIC, SO ITS COLUMNS ARE ITS ROWS, NEEDED FOR
C     THE MULTIPLY, BUT COLUMNS HAVE THE ADVANTAGE OF UNIT STRIDE.
C     IF -INCORE-, ANY SYMMETRIC MATRIX -A- CAN BE GIVEN; OTHERWISE
C     ONLY MATRICES -H1- AND -H2- (WHICH ARE NEEDED FOR DYNAMICAL
C     POLARIZABILITY CALCULATIONS) CAN BE SUBMITTED. -IFH1- TELLS
C     WHETHER H1 (TRUE) OR H2 (FALSE) IS IN WORK.
C     LYUDMILA SLIPCHENKO 09/20/2005
C
      IF(INCORE) THEN
         DO 150 I=1,N
            AX(I) = DDOT(N,A(1,I),1,X(1),1)
  150    CONTINUE
      ELSE
         KOL2 = 0
  310    CONTINUE
            KOL0=KOL2
            KOL1=KOL2+1
            KOL2=KOL2+NWIDTH
            IF(KOL2.GT.N) KOL2=N
            NCOL = KOL2-KOL1+1
C
            LEN = N*NCOL
            IF(IFH1) THEN
                CALL SQREAD(NFT16,A,LEN)
            ELSE
                CALL SQREAD(NFTH2,A,LEN)
            END IF
C
            DO 350 I=1,NCOL
               AX(I+KOL0) = DDOT(N,A(1,I),1,X(1),1)
  350       CONTINUE
         IF(KOL2.LT.N) GO TO 310
C
             IF(IFH1) THEN
                CALL SEQREW(NFT16)
             ELSE
                CALL SEQREW(NFTH2)
             END IF
         END IF
C
      RETURN
      END
C
C
C-----------------------------------------------------------------------
      LOGICAL FUNCTION BRKDN(ALPHA)
      IMPLICIT  DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C-----------------------------------------------------------------------
C     TEST WHETHER ALPHA IS ZERO OR AN ABNORMAL NUMBER, IF YES,
C     THIS ROUTINE WILL RETURN .TRUE.
C-----------------------------------------------------------------------
      BRKDN = .FALSE.
      IF (ALPHA.GT.ZERO) THEN
         BETA = ONE / ALPHA
         IF (.NOT. BETA.GT.ZERO) THEN
            BRKDN = .TRUE.
         ENDIF
      ELSE IF (ALPHA.LT.ZERO) THEN
         BETA = ONE / ALPHA
         IF (.NOT. BETA.LT.ZERO) THEN
            BRKDN = .TRUE.
         ENDIF
      ELSE IF (ALPHA.EQ.ZERO) THEN
         BRKDN = .TRUE.
      ENDIF
      RETURN
      END
C-----END-OF-BRKDN
C
C*MODULE LOCPOL    *DECK SEMIIT
      SUBROUTINE SEMIIT(H1,H2,FREQ,NWIDTH,NFT16,NFTH2,
     *                  Y,NROT,QOUT,QIN,M0,M1Z,DIMAT,QREP,
     *                  TMP1,TMP2,IPVT,MXDII)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION  M0, M1Z
C
      DIMENSION H1(NROT,NWIDTH), H2(NROT,NWIDTH)
      DIMENSION Y(NROT),QIN(NROT),QOUT(NROT),M0(NROT),M1Z(NROT),
     *          DIMAT(MXDII+1,MXDII+1),QREP(*),
     *          TMP1(MXDII+1),TMP2(MXDII+1,MXDII+1),IPVT(MXDII+1)
C
      LOGICAL GOPARR, DSKWRK, MASWRK, SVDSKW
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
C     SOLVE A LINEAR EQUATION:
C        M*X = Y
C     VIA SEMI-ITERATIVE DIIS METHOD
C
C     HUI LI
C     IOWA STATE UNIVERSITY
C     SEPTEMBER 14, 2005
C
C     M : NROT*NROT SQUARE MATRIX
C     X : SOLUTION VECTOR --> QOUT
C     Y : A KNOWN VECTOR
C
C     1. PARTITION M INTO DIAGONAL M0 + OFF-DIAGONAL M1
C     2. X(1) = [Y - M1*X(0)]/M0
C     3. X(2) = [Y - M1*X(1)]/M0
C     4. X(3) = [Y - M1*Z(2)]/M0,  WITH Z(2) = DIIS[X(1), X(2)]
C     5. X(4) = [Y - M1*Z(3)]/M0,  WITH Z(3) = DIIS[X(1), X(2), X(3)]
C     6. X(N) = [Y - M1*Z(N-1)]/M0,  WITH Z(N-1) = DIIS[X(1),... X(N-1)]
C     7. CONVERGE IF SQRT{[X(N) - Z(N-1)]**2} < 1.0D-09
C
C
      DO I=1,NROT
         QIN(I) = Y(I)/(M0(I) + FREQ)
      ENDDO
C
C     ITERATIVE CYCLE
C
      MXITER=1000
      THRES = 1.0D-09
      DO ITER=1,MXITER
C
C        - COMPUTE M1*Z = H2*H1*Z
C
         DO I = 1, NROT
            M1Z(I) = 0.0D+00
         END DO
C
         SVDSKW = DSKWRK
         DSKWRK = .FALSE.
C
C        - FIRST H1*Z
         CALL SEQREW(NFT16)
         KOL2 = 0
  400    CONTINUE
         KOL0 = KOL2
         KOL2 = KOL2 + NWIDTH
         IF(KOL2.GT.NROT) KOL2=NROT
         NCOL = KOL2 - KOL0
         LEN = NROT*NCOL
         CALL SQREAD(NFT16,H1,LEN)
         DO 500 IROT=1,NROT
         DO 520 KK=1,NCOL
C  LVS
            M1Z(IROT) = M1Z(IROT) + H1(IROT,KK)*QIN(KOL0+KK)
  520    CONTINUE
  500    CONTINUE
         IF(KOL2.LT.NROT) GO TO 400
         CALL SEQREW(NFT16)
C
         DO KK = 1, NROT
            QOUT(KK) = M1Z(KK)
            M1Z(KK) = 0.0D+00
         END DO
C
C        - SECOND H2*(H1*Z)
         CALL SEQREW(NFTH2)
         LKOL2 = 0
  410    CONTINUE
         LKOL0 = LKOL2
         LKOL2 = LKOL2 + NWIDTH
         IF(LKOL2.GT.NROT) LKOL2=NROT
         LNCOL = LKOL2 - LKOL0
         LLEN = NROT*LNCOL
         CALL SQREAD(NFTH2,H2,LLEN)
         DO 510 LL=1,LNCOL
         DO 530 IROT=1,NROT
            M1Z(LKOL0+LL) = M1Z(LKOL0+LL) + H2(IROT,LL)*QOUT(IROT)
  530    CONTINUE
  510    CONTINUE
         IF(LKOL2.LT.NROT) GO TO 410
         CALL SEQREW(NFTH2)
C
         DSKWRK = SVDSKW
C
C
C        - COMPUTE THE NEXT ITERATE
         DO I=1,NROT
            AM0 = 1.0D+00/(M0(I)+FREQ)
            QOUT(I)=(Y(I)-M1Z(I))*AM0 + QIN(I) - AM0*FREQ*QIN(I)
         ENDDO
C
C        - TEST CONVERGENCE
         SDEV=0.0D+00
         DO I=1,NROT
            SDEV=SDEV+(QIN(I)-QOUT(I))**2
         ENDDO
         SDEV=SQRT(SDEV/NROT)
C
C        IF(MASWRK) WRITE(IW,*) ' THRES=',THRES,' SDEV=', SDEV
         IF (SDEV.LT.THRES) THEN
            IF(MASWRK) WRITE(IW,*)
     *                 'SEMIIT CONVERGED AFTER ',ITER,'ITERATIONS'
            CALL DCOPY(NROT,QIN,1,QOUT,1)
            RETURN
         END IF
C
C        DIIS INTERPOLATION
C
         IF(MXDII.GT.0) THEN
           NTSSAV = NTS
           NTS    = NROT
           CALL ASCDII(NROT,ITER,MXDII,QOUT,QIN,DIMAT,
     *                 QREP,TMP1,TMP2,IPVT)
           NTS    = NTSSAV
         END IF
C
C        COPY THE CHARGES
C
         IF(ITER.NE.MXITER) THEN
            DO I=1,NROT
               QIN(I)=QOUT(I)
            ENDDO
         END IF
      ENDDO
C
      IF(MASWRK) THEN
         WRITE(IW,*) ' '
         WRITE(IW,*) 'ITERATIVE PROCEDURE NOT CONVERGED IN',
     *                MXITER,' ITERATIONS'
         WRITE(IW,*) 'GAMESS WILL PROCEED DESPITE THIS ...'
         WRITE(IW,*) ' '
      END IF
      RETURN
      END
C
C*MODULE LOCPOL    *DECK ADDFR
      SUBROUTINE ADDFR(H21IM,NROT,FREQ,IFADD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION  H21IM(NROT,NROT)
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW,IFADD
      PARAMETER (FOUR=4.0D+00)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ---------ADD FREQENCIES TO H21--------------------------
C
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
C
        DO 20 MROT=1,NROT
          IF(IFADD) THEN
            H21IM(MROT,MROT) = H21IM(MROT,MROT) + FOUR*FREQ*FREQ
          ELSE
            H21IM(MROT,MROT) = H21IM(MROT,MROT) - FOUR*FREQ*FREQ
          END IF
  20    CONTINUE
        DSKWRK = SVDSKW
C
      RETURN
      END
