C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
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 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C 23 JUL 04 - MWS - GRIDMAKE: REMOVE FTNCHEK WARNINGS
C 26 MAR 02 - KRG - USE ABRT CALL
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 10 APR 00 - MWS - ALLOCATE DYNAMIC MEMORY FOR NUMPRP TO USE
C 25 MAR 00 - KKB/LNB - AMOS DIST. MULTIPOLE ANALYSIS
C
C*MODULE DMULTI  *DECK NUMPROP
      SUBROUTINE NUMPROP
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      NATM = NAT
      CALL VALFM(LOADFM)
      LRADIUS   = LOADFM    + 1
      LATCHRGF  = LRADIUS   + NATM
      LATCHRG   = LATCHRGF  + NATM*3
      LATCHRGD  = LATCHRG   + NATM
      LATDIP    = LATCHRGD  + NATM*NATM
      LATDIPD   = LATDIP    + NATM*3
      LATDIPN   = LATDIPD   + NATM*NATM*3
      LATDIPM   = LATDIPN   + NATM*3*3
      LATCHRGN  = LATDIPM   + NATM*3*3
      LATQUADN  = LATCHRGN  + NATM*3
      LTATQUADN = LATQUADN  + NATM*3*3*3
      LATQUAD   = LTATQUADN + NATM*3*3*3
      LTATQUAD  = LATQUAD   + NATM*3*3
      LSATQUAD  = LTATQUAD  + NATM*3*3
      LATOCT    = LSATQUAD  + NATM*5
      LATHEX    = LATOCT    + NATM*7
      LATALPHA  = LATHEX    + NATM*9
      LAST      = LATALPHA  + NATM*3*3
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
      CALL NUMPRP(X(LRADIUS),X(LATCHRGF),X(LATCHRG),X(LATCHRGD),
     *            X(LATDIP),X(LATDIPD),X(LATDIPN),X(LATDIPM),
     *            X(LATCHRGN),X(LATQUADN),X(LTATQUADN),X(LATQUAD),
     *            X(LTATQUAD),X(LSATQUAD),X(LATOCT),X(LATHEX),
     *            X(LATALPHA),NATM)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE DMULTI  *DECK NUMPRP
      SUBROUTINE NUMPRP(RADIUS,ATCHRGF,ATCHRG,ATCHRGD,ATDIP,ATDIPD,
     *                  ATDIPN,ATDIPM,ATCHRGN,ATQUADN,TATQUADN,
     *                  ATQUAD,TATQUAD,SATQUAD,ATOCT,ATHEX,ATALPHA,NATM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*8 RTYPE
C
      PARAMETER (MXATM=2000, NPPA=1082)
C
      LOGICAL GDFT,KHNSHM,MOBAS,DFTGRD,LMOS
C
      DIMENSION RADIUS(NATM),ATCHRGF(NATM,3),ATCHRG(NATM),
     *          ATCHRGD(NATM,NATM),ATDIP(3*NATM),ATDIPD(3*NATM,NATM),
     *          ATDIPN(3*NATM,3),ATDIPM(3,NATM,3),ATCHRGN(NATM,3),
     *          ATQUADN(3,3,NATM,3),TATQUADN(3,3,NATM,3),
     *          ATQUAD(3,3,NATM),TATQUAD(3,3,NATM),SATQUAD(5,NATM),
     *          ATOCT(7,NATM),ATHEX(9,NATM),ATALPHA(3,NATM,3)
      DIMENSION BRAGG(36),QUADN(3,3,3),TQUADN(3,3,3),
     *          TOTCHGN(3),DIPXN(3),DIPYN(3),DIPZN(3)
C
      COMMON /COSMO1/ SE2,SECORR,ETOTS,CDUM,QVCOSMO,
     *                CSPOT(NPPA),ICORR,ITRIPO,ITRIP2,ITRIP3,ITRIP4,
     *                NATCOS,NQS,ITERC
      COMMON /DFTCWM/ GDFT,IQUAD,NRHX,NROX,NRSX,NTMAXX,
     *                NPHIX,IDFSYM,KHNSHM,MOBAS,DFTGRD,KSFUNC(6)
      COMMON /FMCOM / WORK(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PRPSEC/ ISECD1,ITYPD1,ISECD2,ITYPD2,LMOS,ISECV,ITYPV,
     *                ISECD3,ITYPD3
      COMMON /SQRS  / SQRP375,SQR3P75,SQR15,SQRP625,SQR516,
     *                SQR1P25,SQR4P37,SQR3564,SQR8P75
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
C****** BRAGG CONTAINS BRAGG-SLATER RADII FOR ATOMS. FOLLOWING
C****** BECKE, THE HYDROGEN RADIUS IS 0.35. FOR THE NOBLE GASES THE
C****** VALUES HAVE BEEN GUESSED
      DATA BRAGG/ 0.35D+00,0.35D+00,
     *  1.45D+00,1.05D+00,0.85D+00,0.70D+00,0.65D+00,0.60D+00,
     *  0.50D+00,0.45D+00,
     *  1.80D+00,1.50D+00,1.25D+00,1.10D+00,1.00D+00,1.00D+00,
     *  1.00D+00,1.00D+00,
     *  2.20D+00,1.80D+00,1.60D+00,1.40D+00,1.35D+00,1.40D+00,
     *  1.40D+00,1.40D+00,1.35D+00,1.35D+00,1.35D+00,1.35D+00,
     *  1.30D+00,1.25D+00,1.15D+00,1.15D+00,1.15D+00,1.15D+00/
C
C**** WRITTEN BY C. W. MURRAY AND D. J. TOZER.
C**   CHANGED FOR MY OWN PURPOSES; CAROLINE VAN HEUSDEN
C
C**** THIS ROUTINE WORKS OUT THE MEMORY ALLOCATIONS AND LENGTH OF
C**** VECTOR LOOP FOR THE QUADRATURE PROGRAM AND ALL TIME-INTENSIVE
C**** PARTS OF THE DFT CODE
C
      RTYPE='NUMPROP'
      IQUAD=1
C
      WRITE(6,*)
      WRITE(6,*)"************************"
      WRITE(6,*)"EVALUATING MULTIPOLES..."
      WRITE(6,*)"************************"
      WRITE(6,*)
C
C     ICORR IS TAKEN OUT AS IT ALWAYS SHOULD BE 1 AFTER A NEW CALL
C     TO CONSTS IS MADE IN RHFCL
C
C     ITRIP2=1 IS THE POINTER FOR SAVING THE RESULTS IN GRIDMAKE
C
      ITRIP2=1
C
C     COMPUTING SQUARE ROOTS
C
         SQRP375= SQRT(0.375D+00)
         SQR3P75= SQRT(3.75D+00)
         SQR15  = SQRT(15.0D+00)
         SQRP625= SQRT(0.625D+00)
         SQR516 = SQRT(5.0D+00/16.0D+00)
         SQR1P25= SQRT(1.25D+00)
         SQR4P37= SQRT(4.375D+00)
         SQR3564= SQRT(35.0D+00/64.0D+00)
         SQR8P75= SQRT(8.75D+00)
C
C       SOME MEMORY FOR GETTING DENSITY
C       NOTE THAT N2 IS NUMBER OF ELEMENTS -BELOW- DIAGONAL.
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      N2 = NATM*(NATM-1)/2
      N3 = NATM*NATM
C
      CALL VALFM(LOADFM)
      LA     = LOADFM + 1
      LB     = LA     + L2
      IPASC  = LB     + L2
      IRIJ   = IPASC  + 20
      IAIJ   = IRIJ   + N2
      ITHETA = IAIJ   + N2
      ITHEWT = ITHETA + 256
      IIPHI  = ITHEWT + 256
      IRBC   = IIPHI  + 512
      ITOP   = IRBC   + N3
C
C     MINDEN=.FALSE.
C     BRAG=.FALSE.
C
C****** IVECL IS THE NUMBER OF QUADRATURE POINTS TREATED AT ONCE
C     IVECL=(MEM-ITOP-2*NUM*NUM-3*L2)/(10+2*NUM)
C     IVECL=(MEM-ITOP+LOADFM-2*NUM*NUM-3*L2)/(10+2*NUM)
C     IF (IVECL.GT.256) IVECL=256
C
      IVECL=256
C
      IX     = ITOP
      IY     = IX     + IVECL
      IZ     = IY     + IVECL
      IWT    = IZ     + IVECL
      IPSMU  = IWT    + IVECL
      IRJ    = IPSMU  + NATM*IVECL
      IACCUM = IRJ    + NATM*IVECL
      IXMUN  = IACCUM + IVECL
      IXMU2  = IXMUN  + IVECL
      IWORK  = IXMU2  + IVECL
      IIPVT  = IWORK  + 10*IVECL+NUM*IVECL*2+NUM*NUM+L2
      IDELR  = IIPVT  + NATM
      LAST   = IDELR  + NATM*3
      NEED   = LAST - LOADFM -1
      CALL GETFM(NEED)
C
      ISDF=7
      ISDL=7
      IF(RTYPE.EQ.'NUMPOL')THEN
         ISDF=31
         ISDL=33
      ENDIF
      DO ISECD1=ISDF,ISDL
      IF (ISECD1.EQ.7.OR.ISECD1.EQ.33) THEN
      WRITE(6,19)
 19   FORMAT(/'NUMERICAL ANALYSIS WITH ORIGINAL BRAGG RADIUS')
      ENDIF
C
C         OBTAIN WAVEFUNCTION DENSITY
C
      CALL DENSGET(SCFTYP,WORK(LA),WORK(LB),L2)
C
C     BRAG=.TRUE.
      DO I=1,NATM
         RADIUS(I)=BRAGG(NINT(ZAN(I)))
         IF (ISECD1.EQ.7.OR.ISECD1.EQ.33) THEN
            WRITE(6,15) RADIUS(I)
 15         FORMAT('RADIUS ',' = ',F13.6)
         ENDIF
      ENDDO
C
      WRITE(6,*) "CALLING GRIDMAKE, IVECL =",IVECL
C
      CALL GRIDMAKE(RADIUS,ATCHRG,ATDIP,ATDIPN,TOTCHGN,
     *              DIPXN,DIPYN,DIPZN,ATCHRGN,QUADN,ATQUADN,TQUADN,
     *              TATQUADN,WORK(LA),WORK(IX),WORK(IY),WORK(IZ),
     *              WORK(IWT),WORK(IWORK),WORK(IPSMU),WORK(IRJ),
     *              WORK(IACCUM),WORK(IXMUN),WORK(IXMU2),WORK(IPASC),
     *              WORK(IRIJ),WORK(IAIJ),WORK(ITHETA),
     *              WORK(ITHEWT),WORK(IIPHI),IVECL,WORK(IRBC),
     *              .TRUE.,AVPOL,ATQUAD,TATQUAD,SATQUAD,ATOCT,ATHEX,
     *              ATALPHA,L1,L2,NATM)
C
       IF(RTYPE.EQ.'NUMPOL')THEN
          DO I=1,NATM
             ATCHRGF(I,ISECD1-30)=-ATCHRG(I)
          ENDDO
       ENDIF
      ENDDO
C
C-------------------------------------------------
C
      IF(RTYPE.EQ.'NUMPOL')THEN
      LMOS=.FALSE.
      ISECD1=7
C
      CALL DENSGET(SCFTYP,WORK(LA),WORK(LB),L2)
C
      STEP=0.01D+00
      DO NSTEP=1,2
      DO J=1,NATM
      DO I=1,NATM
         RADIUS(I)=BRAGG(NINT(ZAN(I)))
      ENDDO
      IF(NSTEP.EQ.1)THEN
         RADIUS(J)=RADIUS(J)+STEP
      ELSE
         RADIUS(J)=RADIUS(J)-STEP
      ENDIF
C
      CALL GRIDMAKE(RADIUS,ATCHRG,ATDIP,ATDIPM,TOTCHGN,
     *              DIPXN,DIPYN,DIPZN,ATCHRGN,QUADN,ATQUADN,TQUADN,
     *              TATQUADN,WORK(LA),WORK(IX),WORK(IY),WORK(IZ),
     *              WORK(IWT),WORK(IWORK),WORK(IPSMU),WORK(IRJ),
     *              WORK(IACCUM),WORK(IXMUN),WORK(IXMU2),
     *              WORK(IPASC),WORK(IRIJ),WORK(IAIJ),WORK(ITHETA),
     *              WORK(ITHEWT),WORK(IIPHI),IVECL,WORK(IRBC),
     *              .FALSE.,AVPOL,ATQUAD,TATQUAD,SATQUAD,ATOCT,ATHEX,
     *              ATALPHA,L1,L2,NATM)
C
      IF(NSTEP.EQ.1)THEN
         DO I=1,NATM
            ATCHRGD(I,J)=ATCHRG(I)
         ENDDO
         DO I=1,NATM*3
            ATDIPD(I,J)=ATDIP(I)
         ENDDO
      ELSE
         DO I=1,NATM
            ATCHRGD(I,J)=(ATCHRGD(I,J)-ATCHRG(I))/(STEP+STEP)
         ENDDO
         DO I=1,NATM*3
            ATDIPD(I,J)=(ATDIPD(I,J)-ATDIP(I))/(STEP+STEP)
         ENDDO
      ENDIF
      ENDDO
      ENDDO
C
C-----------------------------------------------
C     HAVE FOUND DERIVATIVES OF CHARGES AND DIPOLES
C     W.R.T CHANGES IN BRAGG RADII
C-------------------------------------------------
C
      WRITE(6,*)'DERIVATIVES OF ATOMIC CHARGES W.R.T. BRAGG RADII'
      DO I=1,NATM
         WRITE(6,99) (ATCHRGD(I,J),J=1,NATM)
      ENDDO
99    FORMAT(1X,6F10.5)
C
      CALL NEWPOL(ATCHRGF,ATCHRGD,ATDIPD,ATDIPN,WORK(IIPVT),
     *            WORK(IDELR),AVPOL,NATM)
      ENDIF
C
      CALL RETFM(NEED)
      RETURN
      END
C
C*MODULE DMULTI  *DECK NEWPOL
      SUBROUTINE NEWPOL(ATCHRGF,ATCHRGD,ATDIPD,ATDIPN,IPVT,DELR,
     *                  AVPOL,NATM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ATCHRGF(NATM,3),ATCHRGD(NATM,NATM),ATDIPD(3,NATM,NATM),
     *          ATDIPN(3,NATM,3),IPVT(NATM),DELR(NATM,3)
C
C     (C) R D AMOS , CAMBRIDGE , 1995
C
      WRITE(6,343)
  343 FORMAT(1X,'LOCAL ATOMIC DIPOLE POLARISABILITIES'
     1    ,/ 1X,'( TERM FROM CHANGE IN LOCAL DIPOLE)'
     1    ,/ 1X,'                    FX        FY        FZ')
      DO I=1,NATM
         DO ISECNR=1,3
            WRITE(6,3336) ATDIPN(1,I,ISECNR),
     *                    ATDIPN(2,I,ISECNR),ATDIPN(3,I,ISECNR)
         ENDDO
         WRITE(6,*)
      ENDDO
 3336 FORMAT(4X,2X,3F10.5)
C
      CALL F04AEF(ATCHRGD,NATM,ATCHRGF,NATM,NATM-1,3,DELR,
     *            NATM,IPVT,ATCHRGD,NATM,ATCHRGF,NATM,IFAIL)
      DO J=1,3
      DELR(NATM,J)=0.0D+00
      ENDDO
      DO K=1,3
      DO J=1,NATM
      DO L=1,3
      SUM=0.0D+00
      DO I=1,NATM
      SUM=SUM+ATDIPD(K,J,I)*DELR(I,L)
      ENDDO
      ATDIPN(K,J,L)=ATDIPN(K,J,L)+SUM
      ENDDO
      ENDDO
      ENDDO
C
  344 FORMAT(1X,'TOTAL ATOMIC DIPOLE POLARISABILITIES'
     1    ,/ 1X,'( INCLUDING RELAXATION OF BOUNDARIES)'
     1    ,/ 1X,'                    FX        FY        FZ')
      WRITE(6,344)
      DO I=1,NATM
      WRITE(6,*)'POLARISABILITIES'
      AV=0.0D+00
      DO ISECNR=1,3
      WRITE(6,3336) ATDIPN(1,I,ISECNR),
     & ATDIPN(2,I,ISECNR),ATDIPN(3,I,ISECNR)
C      WRITE(76,3336) ATDIPN(1,I,ISECNR),
C     & ATDIPN(2,I,ISECNR),ATDIPN(3,I,ISECNR)
      AV=AV+ATDIPN(ISECNR,I,ISECNR)
      ENDDO
      AV=AV/3.0D+00
2222  FORMAT(1X,'AVERAGE ATOMIC POLARISABILITY ',F10.5)
2223  FORMAT(1X,'FRACTION OF MOLECULAR POLARISABILITY ', F10.5)
      WRITE(6,2222)AV
      WRITE(6,2223) (AV/AVPOL)
      WRITE(6,*)
      ENDDO
      RETURN
      END
C*MODULE DMULTI  *DECK GRIDMAKE
      SUBROUTINE GRIDMAKE(RADIUS,ATCHRG,ATDIP,ATDIPN,TOTCHGN,
     *                    DIPXN,DIPYN,DIPZN,ATCHRGN,QUADN,ATQUADN,
     *                    TQUADN,TATQUADN,DENS,X,Y,Z,WT,WORK,PSMU,RJ,
     *                    ACCUM,XMUIJN,XMUIJ2,XPASC,RIJ,AIJ,THETA,
     *                    THETWT,IPHIDO,IVECL,RBC,WRTOUT,AVPOL,
     *                    ATQUAD,TATQUAD,SATQUAD,ATOCT,ATHEX,ATALPHA,
     *                    L1,L2,NATM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION RADIUS(NATM),ATCHRG(NATM),ATDIP(3,NATM),
     *          ATDIPN(3,NATM,3),TOTCHGN(3),DIPXN(3),DIPYN(3),DIPZN(3),
     *          ATCHRGN(NATM,3),QUADN(3,3,3),ATQUADN(3,3,NATM,3),
     *          TQUADN(3,3,3),TATQUADN(3,3,NATM,3),DENS(L2),
     *          X(IVECL),Y(IVECL),Z(IVECL),WT(IVECL),WORK(*),
     *          PSMU(IVECL,NATM),RJ(IVECL,NATM),ACCUM(IVECL),
     *          XMUIJN(IVECL),XMUIJ2(IVECL),XPASC(20),
     *          RIJ(NATM*(NATM-1)/2),AIJ(NATM*(NATM-1)/2),
     *          THETA(256),THETWT(256),IPHIDO(512),RBC(NATM,NATM),
     *          ATQUAD(3,3,NATM),TATQUAD(3,3,NATM),SATQUAD(5,NATM),
     *          ATOCT(7,NATM),ATHEX(9,NATM), ATALPHA(3,NATM,3)
      DIMENSION QUAD(3,3),TQUAD(3,3),SQUAD(5),SOCT(7),SHEX(9),
     *          A(3,3),WR(3),WI(3),ZV(3,3),FV1(3),IV1(3)
C
      LOGICAL GDFT,KHNSHM,MOBAS,DFTGRD,WRTOUT,LMOS
C
      PARAMETER (PI=3.1415926535898D+00)
      PARAMETER (MXATM=2000, MXSH=5000, MAXDEN=25*MXATM, LENABC=2000,
     *           NPPA=1082)
C
      COMMON /COSMO1/ SE2,SECORR,ETOTS,CDUM,QVCOSMO,
     *                CSPOT(NPPA),ICORR,ITRIPO,ITRIP2,ITRIP3,ITRIP4,
     *                NATCOS,NQS,ITERC
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /DFTCWM/ GDFT,IQUAD,NRHX,NROX,NRSX,NTMAXX,
     *                NPHIX,IDFSYM,KHNSHM,MOBAS,DFTGRD,KSFUNC(6)
      COMMON /PRPSEC/ ISECD1,ITYPD1,ISECD2,ITYPD2,LMOS,ISECV,ITYPV,
     *                ISECD3,ITYPD3
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     *                QDEN(MAXDEN),AR(LENABC),
     *                NSPA,NPS,NPS2,NDEN,NPSPHER
      COMMON /SQRS  / SQRP375,SQR3P75,SQR15,SQRP625,SQR516,
     *                SQR1P25,SQR4P37,SQR3564,SQR8P75
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C
C**** QUAD PERFORMS QUADRATURE FOR THE CALCULATION OF THE EXCHANGE
C**** CORRELATION ENERGY, FOCK MATRIX OR GRADIENT. IT DECIDES THE
C**** POINTS AND THE WEIGHTS AND SENDS A BATCH OF THEM TO DENXYZ
C**** FOR THE CONTRIBUTION.
C
      CALL VCLR(ATCHRG,1,NATM)
      CALL VCLR(ATDIP,1,3*NATM)
      CALL VCLR(ATQUAD,1,3*3*NATM)
      CALL VCLR(QUAD,1,3*3)
      CALL VCLR(TATQUAD,1,3*3*NATM)
      CALL VCLR(TQUAD,1,3*3)
      CALL VCLR(SQUAD,1,5)
      CALL VCLR(SATQUAD,1,5*NATM)
      CALL VCLR(SOCT,1,7)
      CALL VCLR(SHEX,1,9)
      CALL VCLR(ATOCT,1,7*NATM)
      CALL VCLR(ATHEX,1,9*NATM)
      TOTCHG=0.0D+00
      DIPX=0.0D+00
      DIPY=0.0D+00
      DIPZ=0.0D+00
C
      IF(ISECD1.EQ.7)THEN
      DO I=1,NATM
         DIPX=DIPX+ZAN(I)*C(1,I)
         DIPY=DIPY+ZAN(I)*C(2,I)
         DIPZ=DIPZ+ZAN(I)*C(3,I)
C
C MAKE SECOND MOMENTS TRACELESS
C
         TRACE=0.5D+00*ZAN(I)*(C(1,I)**2+C(2,I)**2+C(3,I)**2)
         TQUAD(1,1)=TQUAD(1,1)+1.5D+00*ZAN(I)*C(1,I)*C(1,I)-TRACE
         TQUAD(1,2)=TQUAD(1,2)+1.5D+00*ZAN(I)*C(1,I)*C(2,I)
         TQUAD(1,3)=TQUAD(1,3)+1.5D+00*ZAN(I)*C(1,I)*C(3,I)
         TQUAD(2,2)=TQUAD(2,2)+1.5D+00*ZAN(I)*C(2,I)*C(2,I)-TRACE
         TQUAD(2,3)=TQUAD(2,3)+1.5D+00*ZAN(I)*C(2,I)*C(3,I)
         TQUAD(3,3)=TQUAD(3,3)+1.5D+00*ZAN(I)*C(3,I)*C(3,I)-TRACE
         XX=C(1,I)*C(1,I)*ZAN(I)
         YY=C(2,I)*C(2,I)*ZAN(I)
         ZZ=C(3,I)*C(3,I)*ZAN(I)
         RR=XX+YY+ZZ
         SOCT(1)=SOCT(1)+2.5D+00*C(3,I)*ZZ-1.5D+00*C(3,I)*RR
         SOCT(2)=SOCT(2)+SQRP375*C(1,I)*(5.0D+00*ZZ-RR)
         SOCT(3)=SOCT(3)+SQRP375*C(2,I)*(5.0D+00*ZZ-RR)
         SOCT(4)=SOCT(4)+SQR3P75*C(3,I)*(XX-YY)
         SOCT(5)=SOCT(5)+SQR15*C(1,I)*C(2,I)*C(3,I)*ZAN(I)
         SOCT(6)=SOCT(6)+SQRP625*C(1,I)*(XX-3.0D+00*YY)
         SOCT(7)=SOCT(7)+SQRP625*C(2,I)*(3.0D+00*XX-YY)
         XX=C(1,I)*C(1,I)
         YY=C(2,I)*C(2,I)
         ZZ=C(3,I)*C(3,I)
         RR=XX+YY+ZZ
         SHEX(1)=SHEX(1)+0.125D+00*(35.0D+00*ZZ*ZZ-30.0D+00*ZZ*RR
     1      +3.0D+00*RR*RR)*ZAN(I)
         SHEX(2)=SHEX(2)+SQRP625*C(1,I)*C(3,I)*(7.0D+00*ZZ
     1     -3.0D+00*RR)*ZAN(I)
         SHEX(3)=SHEX(3)+SQRP625*C(2,I)*C(3,I)*(7.0D+00*ZZ
     1  -3.0D+00*RR)*ZAN(I)
         SHEX(4)=SHEX(4)+SQR516*(XX-YY)*(7.0D+00*ZZ-RR)
     1    *ZAN(I)
         SHEX(5)=SHEX(5)+SQR1P25*C(1,I)*C(2,I)*(7.0D+00*ZZ-RR)
     1    *ZAN(I)
         SHEX(6)=SHEX(6)+SQR4P37*C(1,I)*C(3,I)*(XX-3.0D+00*YY)
     1    *ZAN(I)
         SHEX(7)=SHEX(7)+SQR4P37*C(2,I)*C(3,I)*(3.0D+00*XX-YY)
     1    *ZAN(I)
         SHEX(8)=SHEX(8)+SQR3564*(XX*XX-6.0D+00*XX*YY
     1     +YY*YY)*ZAN(I)
         SHEX(9)=SHEX(9)+SQR8P75*C(1,I)*C(2,I)*(XX-YY)
     1    *ZAN(I)
      ENDDO
      ENDIF
C****** OF QUADRATURE
C
      NT=1
C     NTDFT=NT
C
C     COSMO CHANGE - KIM
C       WE CAN SET NT = 1 FOR OUR PURPOSES.
C       WE WANT TO ALWAYS HAVE IQUAD=1
C
C      IF (IDFSYM.EQ.1) NT=1
C
C      IF (NROX.NE.0) IQUAD=6
C
      IF (IQUAD.GT.0) THEN
          ILIM1=IQUAD
          ILIM2=IQUAD
      ELSE
          ILIM1=1
          ILIM2=5
      ENDIF
C
C NRTYPE DETERMINES THE POWER OF THE TRANSFORMATION IN EULER-MACAURIN
C RADIAL INTEGRATION. NRTYPE=1 IMPLIES THAT GAUSS-LAGUERRE INTEGRATION
C IS PERFORMED FOR HYDROGEN; NRTYPE.NE.1 MEANS EULER-MACLAURIN
C INTEGRATION IS PERFORMED. NTRANS DETERMINES THE POWER OF (1-U*U) USED
C IN THE FUZZY VORONOI POLYHEDRA SECTION OF CODE
C
C    EXPERIMENTAL VALUE OF NTRANS
      NTRANS=10
      BOHR=0.529177249D+00
C
C****** FORM PASCALS TRIANGLE IN XPASC FOR FUZZY VORONOI POLYHEDRA CODE
C
      ISIGN=-1
      DO 5 I=NTRANS,1,-1
      TEMP = NFAC(NTRANS)
      TEMP = TEMP/(NFAC(I)*NFAC(NTRANS-I))
      ISIGN=-ISIGN
      XPASC(I+1)=ISIGN*TEMP
      XPASC(1)=1.0D+00
5     CONTINUE
C
      APASC=0.0D+00
      DO 6 I=1,NTRANS+1
         XPASC(I)=XPASC(I)/(2*I-1)
         APASC=APASC+XPASC(I)
6     CONTINUE
      APASC=0.5D+00/APASC
      ITEMP=0
      CALL VCLR(RBC,1,NATM*NATM)
      DO 7 INA=1,NATM
      DO 7 JNA=1,INA-1
          ITEMP=ITEMP+1
          RIJ(ITEMP)=1.0D+00/(SQRT(
     1        (C(1,INA)-C(1,JNA))**2
     1       +(C(2,INA)-C(2,JNA))**2
     1       +(C(3,INA)-C(3,JNA))**2))
          RBC(INA,JNA)=1.0D+00/RIJ(ITEMP)
          RBC(JNA,INA)=1.0D+00/RIJ(ITEMP)
          CHI=RADIUS(INA)/RADIUS(JNA)
          TEMP=(CHI-1.0D+00)/(CHI+1.0D+00)
          AIJ(ITEMP)=TEMP/(TEMP*TEMP-1.0D+00)
          IF (AIJ(ITEMP).GT.0.5D+00) AIJ(ITEMP)=0.5D+00
          IF (AIJ(ITEMP).LT.-0.5D+00) AIJ(ITEMP)=-0.5D+00
7     CONTINUE
C
C***** NRTYP1 DETERMINES THE COARSENESS OF THE QUADRATURE GRID
C
      WRITE(6,*)
     *"NRTYP1 DETERMINES THE COARSENESS OF THE QUADRATURE GRID"
      DO 90 NRTYP1=ILIM1,ILIM2
      IF (NRTYP1.EQ.1) THEN
      NRO=64
      NRH=64
      NRS=74
      NTMAX=24
      NPMAX=48
      ENDIF
      IF (NRTYP1.EQ.2) THEN
      NRO=25
      NRH=20
      NRS=30
      NTMAX=8
      NPMAX=16
      ENDIF
      IF (NRTYP1.EQ.3) THEN
      NRO=25
      NRH=20
      NRS=30
      NTMAX=8
      NPMAX=16
      ENDIF
      IF (NRTYP1.EQ.4) THEN
      NRO=36
      NRH=31
      NRS=41
      NTMAX=12
      NPMAX=24
      ENDIF
      IF (NRTYP1.EQ.5) THEN
      NRO=54
      NRH=46
      NRS=60
      NTMAX=16
      NPMAX=32
      ENDIF
      IF (NRTYP1.EQ.6) THEN
      NRO=NROX
      NRH=NRHX
      NRS=NRSX
      NTMAX=NTMAXX
      NPMAX=NPHIX
      ENDIF
C
C****** LOOP OVER THE NUMBER OF ATOMS
C
      NPOINT=0
      DO 10 NATOM=1,NATM
         NATOMX = NATOM
C
C*********** PERFORM SYMETRY TEST
C
         SYMWT=1.0D+00
         IPHIS=1
         ITHES=1
         IPT=0
C         IF (NT.GT.1) CALL SYMDF(SYMWT,IPHIS,ITHES,NATOMX,NPMAX)
         IF (IPHIS.EQ.0) GOTO 10
C
         NRREAL=0
         BQEXP=RADIUS(NATOM)/BOHR
         NR=NRO
C HYDROGEN TREATED DIFFERENTLY, POTENTIALLY WITH GAUSS LAGUERRE POINTS
C FROM THE NAG ROUTINE BELOW
         IF (ZAN(NATOM).LT.1.5D+00) THEN
            NR=NRH
         ENDIF
         IF (ZAN(NATOM).GT.10.5D+00) NR=NRS
C
C********* RADIAL INTEGRATION
C
         DO 20 IR=1,NR
C************* EULER-MACLAURIN RADIAL INTEGRATION
                XNODE=IR
                XNODE=XNODE/(NR+1)
                RNODE=BQEXP*(XNODE/(1.0D+00-XNODE))**2
                IF (RNODE.GT.20.0D+00) GOTO 20
                RWT=2.0D+00*BQEXP**3*XNODE**5/(1.0D+00-XNODE)**7
                RWT=RWT/(NR+1)
            NRREAL=NRREAL+1
C
C****** ESTABLISH NUMBER OF THETA POINT ACCORDING TO RADIAL POSITION
C
            TEMP=RNODE/RADIUS(NATOM)*BOHR*NTMAX*5.0D+00
            NTHETA=NINT(TEMP)+1
            PHIZER=0.0D+00
            CALL SYMDF2(NTHETA,NTMAX,NPHI,NPMAX,
     *                  NATOMX,IPHIS,IPHIDO,PHIZER)
C
C****** FIND THETA WEIGHTS AND NODES BY GAUSS-LEGENDRE QUADRATURE
C
            CALL LEGWT(NTHETA,THETWT,THETA)
C
C****** THETA INTEGRATION
C
             DO 30 IQ=1,NTHETA
C***************** THE NEXT LINE IS THE THETA SYMMETRY TEST
                IF (ITHES.EQ.2.AND.THETA(IQ).LT.0) GOTO 30
                THETAN=ACOS(THETA(IQ))
C
C****** FIND TOTAL WEIGHT FOR THE PHI,THETA,R POINT.
                PHIWT=(2.0D+00*PI)/NPHI
                TOTWT=RWT*THETWT(IQ)*PHIWT*SYMWT
C
C****** PHI INTEGRATION
C
                DO 40 IPHI=1,NPHI
C***************** THE NEXT LINE IS THE PHI SYMMETRY TEST
                  IF (IPHIDO(IPHI).EQ.0) GOTO 40
                  IPT=IPT+1
                  PHINOD=PI*((2.0D+00*IPHI)/NPHI-1.0D+00)
                  PHINOD=PHINOD+PHIZER
                  XX=RNODE*SIN(THETAN)*COS(PHINOD)
                  YY=RNODE*SIN(THETAN)*SIN(PHINOD)
                  ZZ=RNODE*COS(THETAN)
                  X(IPT)=XX+C(1,NATOM)
                  Y(IPT)=YY+C(2,NATOM)
                  Z(IPT)=ZZ+C(3,NATOM)
                  WT(IPT)=TOTWT
                  IF (IPHIDO(IPHI).EQ.2) WT(IPT)=WT(IPT)*0.5D+00
C
                  IF (IPT.EQ.IVECL) THEN
C
C******************* HAVE A SUFFICIENTLY LARGE BATCH OF POINTS
C
      CALL WEIGHT(RJ,PSMU,RIJ,AIJ,APASC,XPASC,X,Y,Z,WT,ACCUM,
     *            XMUIJN,XMUIJ2,IVECL,IPT,NTRANS,NATOMX,NATM)
C
      CALL GRIDUSE(DENS,WORK,WORK(NUM*IPT+1),X,Y,Z,WT,IPT,NATOMX,
     *             ATCHRG,ATDIP,TATQUAD,TOTCHG,DIPX,DIPY,DIPZ,TQUAD,
     *             SOCT,SHEX,ATOCT,ATHEX,L1,L2,NATM)
      NPOINT=NPOINT+IPT
      IPT=0
      ENDIF
C
C*** PHI CLOSE
40             CONTINUE
C
C*** THETA CLOSE
30          CONTINUE
C
C*** FORM WHOLE CONTRIBUTION FROM THIS RADIAL POINT FOR THE RADIAL TEST
C*** RADIAL CLOSE
20       CONTINUE
C
C******************* HAVE SOME POINTS LEFT OVER!
C
      IF (IPT.NE.0) THEN
          CALL WEIGHT(RJ,PSMU,RIJ,AIJ,APASC,XPASC,X,Y,Z,WT,ACCUM,
     *                XMUIJN,XMUIJ2,IVECL,IPT,NTRANS,NATOM,NATM)
C
          CALL GRIDUSE(DENS,WORK,WORK(NUM*IPT+1),X,Y,Z,WT,IPT,NATOMX,
     *                 ATCHRG,ATDIP,TATQUAD,TOTCHG,DIPX,DIPY,DIPZ,TQUAD,
     *                 SOCT,SHEX,ATOCT,ATHEX,L1,L2,NATM)
          NPOINT=NPOINT+IPT
          IPT=0
      ENDIF
C
C*** ATOM CLOSE
   10 CONTINUE
C
   90 CONTINUE
C
C
C      COSMO CHANGE - KIM
C        WE CAN SET NT=1 FOR OUR PURPOSES.
C
       NT=1
C      NT=NTDFT
C
C WRITE OUT SECTION:
C
      IF(ISECD1.EQ.7) THEN
      DO I=1,NATM
      ATCHRG(I)=ZAN(I)+ATCHRG(I)
      ENDDO
      IF(WRTOUT)THEN
      WRITE(6,99) ISECD1
   99 FORMAT('ISECD1 =',I5)
C
      WRITE(6,321) TOTCHG
      WRITE(6,323)
      WRITE(6,3333) DIPX,DIPY,DIPZ
      WRITE(6,327)
      WRITE(6,3335) TQUAD(1,1),TQUAD(1,2),TQUAD(2,2),TQUAD(1,3),
     1  TQUAD(2,3),TQUAD(3,3)
      WRITE(6,328)
      ONETHIRD=1.0D+00/3.0D+00
      FOURTHIRD=4.0D+00/3.0D+00
      SQUAD(1)=TQUAD(3,3)
      SQUAD(2)=SQRT(FOURTHIRD)*TQUAD(1,3)
      SQUAD(3)=SQRT(FOURTHIRD)*TQUAD(2,3)
      SQUAD(4)=SQRT(ONETHIRD)*(TQUAD(1,1)-TQUAD(2,2))
      SQUAD(5)=SQRT(FOURTHIRD)*TQUAD(1,2)
C
      WRITE(6,3334)(SQUAD(I),I=1,5)
      WRITE(6,4444)(SOCT(I),I=1,7)
      WRITE(6,4445)(SHEX(I),I=1,9)
C
4444  FORMAT(1X,'MOLECULAR SPHERICAL HARMONIC OCTOPOLES'
     1 /6X,'Q30       Q31C      Q31S      Q32C      Q32S'
     1 ,'      Q33C      Q33S'/1X,7F10.5)
4445  FORMAT(1X,'MOLECULAR SPHERICAL HARMONIC HEXADECAPOLES'
     1 /6X,'Q40       Q41C      Q41S      Q42C      Q42S'
     1 ,'      Q43C      Q43S      Q44C      Q44S'/1X,9F10.5)
C
C      WRITE(76,*)"***********DEBUG************"
C      WRITE(76,*)'MOLECULAR MULTIPOLES'
C      WRITE(76,3333) TOTCHG
C      WRITE(76,3333) DIPZ,DIPX,DIPY
C      WRITE(76,3334) (SQUAD(I),I=1,5)
C      WRITE(76,4450) (SOCT(I),I=1,7)
C      WRITE(76,4451) (SHEX(I),I=1,9)
C      WRITE(76,*)"***********DEBUG************"
C
C
C    COSMO CHANGE - KIM
C
C
      IF(ITRIP2.EQ.1) THEN
         DO 3 IC=1,25*NATM
            QDEN(IC)=0.0D+00
 3       CONTINUE
         KSTRT=0
      ENDIF
C
      WRITE(6,331)
      DO I =1,NATM
      WRITE(6,333) I,ATCHRG(I)
      IF(ITRIP2.EQ.1) THEN
         QDEN(I+KSTRT)=ATCHRG(I)
C         WRITE(78,*)"I+KSTRT:",I+KSTRT
         KSTRT=KSTRT+24
      ENDIF
      ENDDO
C
      IF(ITRIP2.EQ.1) JSTRT=1
      WRITE(6,335)
      DO I=1,NATM
      WRITE(6,3336) I,ATDIP(1,I),ATDIP(2,I),ATDIP(3,I)
         IF(ITRIP2.EQ.1) THEN
            QDEN(I+JSTRT)=ATDIP(1,I)
            QDEN(I+JSTRT+1)=ATDIP(2,I)
            QDEN(I+JSTRT+2)=ATDIP(3,I)
C            WRITE(78,*)"I+JSTRT:",I+JSTRT
            JSTRT=JSTRT+24
         ENDIF
      ENDDO
C
C      IF(ITRIP2.EQ.1) ISTRT=2
C
      WRITE(6,3242)
      DO I=1,NATM
      WRITE(6,3338) I,TATQUAD(1,1,I),TATQUAD(1,2,I),
     1  TATQUAD(2,2,I),TATQUAD(1,3,I),TATQUAD(2,3,I),TATQUAD(3,3,I)
C      IF(ITRIP2.EQ.1) THEN
C
C  THE FOLLOWING CODE IS COMMENTED OUT, AS WE ARE NOW USING
C  ONLY SPHERICAL MULTIPOLES:
C
C THERE ARE DIFFERENT ONVENTIONS FOR QUADRUPOLE MOMENTS. ONE WAY IS TO
C WRITE THE QUADRUPOLE POTENTIAL AS
C
C        V(X,Y,Z)= ((3X**2-R**2)/2 * QXX+ ... + 3XY * QXY + ...)
C
C AND ANOTHER IS
C
C                = ( X**2 * Q'XX + ...+ 2XY * Q'XY)
C
C BOTH ARE IDENTICAL, IF THE TRACE OF THE TENSOR IS ZERO AND
C Q'=3/2 OF Q.
C
C IT APPEARS THAT THE CONVENTION IN THIS CODE IS THE LATTER, AND
C THUS WE NEED TO MULTIPLY THESE BY A FACTOR OF 2/3 TO WORK WITH
C COSMO.
C
C         QDEN(I+ISTRT)=TATQUAD(1,1,I)*TOTHIRDS
C         QDEN(I+ISTRT+3)=TATQUAD(2,2,I)*TOTHIRDS
C         QDEN(I+ISTRT+7)=TATQUAD(3,3,I)*TOTHIRDS
C         QDEN(I+ISTRT+2)=TATQUAD(1,2,I)*TOTHIRDS
C         QDEN(I+ISTRT+5)=TATQUAD(1,3,I)*TOTHIRDS
C         QDEN(I+ISTRT+6)=TATQUAD(2,3,I)*TOTHIRDS
C
C         ISTRT=ISTRT+16
C      ENDIF
C
      ENDDO
C
      IF(ITRIP2.EQ.1) ISTRT=4
      WRITE(6,3243)
      DO I=1,NATM
         SATQUAD(1,I)=TATQUAD(3,3,I)
         SATQUAD(2,I)=SQRT(FOURTHIRD)*TATQUAD(1,3,I)
         SATQUAD(3,I)=SQRT(FOURTHIRD)*TATQUAD(2,3,I)
         SATQUAD(4,I)=SQRT(ONETHIRD)*(TATQUAD(1,1,I)-TATQUAD(2,2,I))
         SATQUAD(5,I)=SQRT(FOURTHIRD)*TATQUAD(1,2,I)
         WRITE(6,3337) I,SATQUAD(1,I),SATQUAD(2,I),SATQUAD(3,I),
     &      SATQUAD(4,I),SATQUAD(5,I)
C
         FUDGE=1.00D+00
C
         IF(ITRIP2.EQ.1) THEN
            QDEN(I+ISTRT)=SATQUAD(1,I)*FUDGE
            QDEN(I+ISTRT+1)=SATQUAD(2,I)*FUDGE
            QDEN(I+ISTRT+2)=SATQUAD(3,I)*FUDGE
            QDEN(I+ISTRT+3)=SATQUAD(4,I)*FUDGE
            QDEN(I+ISTRT+4)=SATQUAD(5,I)*FUDGE
C            WRITE(78,*)"I+ISTRT:",I+ISTRT
C            WRITE(78,*)I,(SATQUAD(J,I),J=1,5)
            ISTRT=ISTRT+24
         ENDIF
C
      ENDDO
C
      FUDGE2=1.0D+00
      IF(ITRIP2.EQ.1) LSTRT=9
      WRITE(6,4446)
      DO I=1,NATM
         WRITE(6,4447) I,(ATOCT(J,I),J=1,7)
C         WRITE(78,*)"I+LSTRT:",I+LSTRT
C         WRITE(78,4447)I,(ATOCT(J,I),J=1,7)
         IF(ITRIP2.EQ.1) THEN
            QDEN(I+LSTRT)=ATOCT(1,I)*FUDGE2
            QDEN(I+LSTRT+1)=ATOCT(2,I)*FUDGE2
            QDEN(I+LSTRT+2)=ATOCT(3,I)*FUDGE2
            QDEN(I+LSTRT+3)=ATOCT(4,I)*FUDGE2
            QDEN(I+LSTRT+4)=ATOCT(5,I)*FUDGE2
            QDEN(I+LSTRT+5)=ATOCT(6,I)*FUDGE2
            QDEN(I+LSTRT+6)=ATOCT(7,I)*FUDGE2
            LSTRT=LSTRT+24
         ENDIF
      ENDDO
C
      FUDGE3=1.0D+00
      IF(ITRIP2.EQ.1) MSTRT=16
      WRITE(6,4448)
      DO I=1,NATM
         WRITE(6,4449) I,(ATHEX(J,I),J=1,9)
C        WRITE(78,*)"I+MSTRT:",I+MSTRT
         IF(ITRIP2.EQ.1) THEN
            QDEN(I+MSTRT)=ATHEX(1,I)*FUDGE3
            QDEN(I+MSTRT+1)=ATHEX(2,I)*FUDGE3
            QDEN(I+MSTRT+2)=ATHEX(3,I)*FUDGE3
            QDEN(I+MSTRT+3)=ATHEX(4,I)*FUDGE3
            QDEN(I+MSTRT+4)=ATHEX(5,I)*FUDGE3
            QDEN(I+MSTRT+5)=ATHEX(6,I)*FUDGE3
            QDEN(I+MSTRT+6)=ATHEX(7,I)*FUDGE3
            QDEN(I+MSTRT+7)=ATHEX(8,I)*FUDGE3
            QDEN(I+MSTRT+8)=ATHEX(9,I)*FUDGE3
            MSTRT=MSTRT+24
         ENDIF
      ENDDO
C
C      WRITE(76,*)"***********DEBUG************"
C      DO I=1,NATM
C         WRITE(76,*) 'POSITION'
C        WRITE(76,3333) C(1,I),C(2,I),C(3,I)
C        WRITE(76,*)'MULTIPOLES'
C         WRITE(76,3333) ATCHRG(I)
C         WRITE(76,3333) ATDIP(3,I),ATDIP(1,I),ATDIP(2,I)
C         WRITE(76,3334) (SATQUAD(J,I),J=1,5)
C         WRITE(76,4450) (ATOCT(J,I),J=1,7)
C         WRITE(76,4451) (ATHEX(J,I),J=1,9)
C      ENDDO
C      WRITE(76,*)"***********DEBUG************"
C 4450 FORMAT(1X,7F10.5)
C 4451 FORMAT(1X,9F10.5)
C
 4446 FORMAT(1X,'ATOMIC SPHERICAL HARMONIC OCTOPOLES'
     1 /11X,'Q30       Q31C      Q31S      Q32C      Q32S'
     1 ,'      Q33C      Q33S')
 4448 FORMAT(1X,'ATOMIC SPHERICAL HARMONIC HEXADECAPOLES'
     1 /11X,'Q40       Q41C      Q41S      Q42C      Q42S'
     1 ,'      Q43C      Q43S      Q44C      Q44S')
 4447 FORMAT(2X,I2,2X,7F10.5)
 4449 FORMAT(2X,I2,2X,9F10.5)
  321 FORMAT(1X,'TOTAL NUMBER OF ELECTRONS ',F10.5)
  323 FORMAT(1X,'DIPOLE         X          Y         Z ')
  328 FORMAT(1X,'SPHERICAL QUADRUPOLE ')
  327 FORMAT(1X,'QUADRUPOLE    XX        XY        YY        XZ',
     1'        YZ        ZZ')
 3333 FORMAT(10X,3F10.5)
 3334 FORMAT(10X,5F10.5)
 3335 FORMAT(10X,6F10.5)
  333 FORMAT(2X,I2,2X,F10.5)
  331 FORMAT(' ATOMIC CHARGES')
  335 FORMAT(' ATOMIC DIPOLES'/12X,'X        Y         Z')
 3338 FORMAT(2X,I2,2X,6F10.5)
 3336 FORMAT(2X,I2,2X,3F10.5)
 3242 FORMAT(' ATOMIC QUADRUPOLES'/11X,'XX        XY        YY',
     1 '        XZ        YZ        ZZ')
 3243 FORMAT(' ATOMIC SPHERICAL QUADRUPOLES')
 3337 FORMAT(2X,I2,2X,5F10.5)
C
      ENDIF
      ENDIF
C
C   ..... END OF SCF/DFT DENSITY ANALYSIS
C
      IF(ISECD1.EQ.7)ISECNR=1
      IF(ISECD1.EQ.31) ISECNR=1
      IF(ISECD1.EQ.32) ISECNR=2
      IF(ISECD1.EQ.33) ISECNR=3
      TOTCHGN(ISECNR)=TOTCHG
      DIPXN(ISECNR)=DIPX
      DIPYN(ISECNR)=DIPY
      DIPZN(ISECNR)=DIPZ
      DO J=1,3
      DO K=1,3
      QUADN(J,K,ISECNR)=QUAD(J,K)
      TQUADN(J,K,ISECNR)=TQUAD(J,K)
      ENDDO
      ENDDO
      DO I=1,NATM
      ATCHRGN(I,ISECNR)=ATCHRG(I)
      DO J=1,3
      ATDIPN(J,I,ISECNR)=ATDIP(J,I)
      DO K=1,3
      ATQUADN(J,K,I,ISECNR)=ATQUAD(J,K,I)
      TATQUADN(J,K,I,ISECNR)=TATQUAD(J,K,I)
      ENDDO
      ENDDO
      ENDDO
C
      IF(WRTOUT)THEN
      IF(ISECD1.EQ.33) THEN
      WRITE (88,421) (TOTCHGN(ISECNR),ISECNR=1,3)
      WRITE(6,329)
      DO ISECNR=1,3
      WRITE(6,3333) DIPXN(ISECNR),DIPYN(ISECNR),
     & DIPZN(ISECNR)
      A(ISECNR,1)=DIPXN(ISECNR)
      A(ISECNR,2)=DIPYN(ISECNR)
      A(ISECNR,3)=DIPZN(ISECNR)
      ENDDO
      WRITE(6,*)'MOLECULAR POLARISABILITY'
      DO ISECNR=1,3
      WRITE(6,3333) DIPXN(ISECNR),DIPYN(ISECNR),
     & DIPZN(ISECNR)
      ENDDO
      ND=3
      MATZ=1
      CALL RG(ND,ND,A,WR,WI,MATZ,ZV,IV1,FV1,IERR)
      AV=0.0D+00
      DO K=1,3
      IF(ABS(WI(K)).GT.1.0D-08)THEN
C      WRITE(6,*) 'COMPLEX EIGENVALUES'
      ELSE
      AV=AV+WR(K)
      ENDIF
      ENDDO
      AVPOL=AV/3.0D+00
      WRITE(6,*)'AVERAGE MOLECULAR POLARISABILITY ',AVPOL
      WRITE(6,341)
      DO ISECNR=1,3
      WRITE(6,3335) TQUADN(1,1,ISECNR),TQUADN(1,2,ISECNR),
     & TQUADN(2,2,ISECNR),TQUADN(1,3,ISECNR),
     & TQUADN(2,3,ISECNR),TQUADN(3,3,ISECNR)
      ENDDO
C
      WRITE(6,422)
      DO I =1,NATM
      WRITE(6,3336) I,ATCHRGN(I,1),
     & ATCHRGN(I,2),ATCHRGN(I,3)
      ENDDO
      WRITE(6,343)
      DO I=1,NATM
      DO ISECNR=1,3
      WRITE(6,3336) I,ATDIPN(1,I,ISECNR),
     & ATDIPN(2,I,ISECNR),ATDIPN(3,I,ISECNR)
      ENDDO
      WRITE(6,*)
      ENDDO
      WRITE(6,424)
      DO I=1,NATM
      DO ISECNR=1,3
      DO K=1,3
      ATALPHA(K,I,ISECNR)=ATCHRGN(I,ISECNR)*C(K,I)+
     & ATDIPN(K,I,ISECNR)
      ENDDO
      ENDDO
      DO ISECNR=1,3
      WRITE(6,3336) ATALPHA(1,I,ISECNR),
     & ATALPHA(2,I,ISECNR),ATALPHA(3,I,ISECNR)
      ENDDO
      WRITE(6,*)
      ENDDO
      ENDIF
      ENDIF
C
  421 FORMAT(1X,'CHARGE CONSERVATION'
     1    ,/ 1X,'               FX        FY        FZ',
     1     /10X,3F10.5)
  422 FORMAT(1X,'DERIVATIVES OF ATOMIC CHARGES'
     1    ,/ 1X,'                    FX        FY        FZ')
  424 FORMAT(1X,'TOTAL ATOMIC DIPOLE POLARISABILITIES',
     1    ' (ORIGINAL DEFINITION)'
     1    ,/ 1X, '(INCLUDING DERIVATIVES OF CHARGES)'
     1    ,/ 1X,'                    FX        FY        FZ')
  343 FORMAT(1X,'LOCAL ATOMIC DIPOLE POLARISABILITIES'
     1    ,/ 1X,'( TERM FROM CHANGE IN LOCAL DIPOLE)'
     1    ,/ 1X,'                    FX        FY        FZ')
  329 FORMAT(1X,'MOLECULAR DIPOLE POLARISABILITY'
     1    ,/ 1X,'               FX        FY        FZ')
  341 FORMAT(1X,'MOLECULAR DIPOLE-QUADRUPOLE POLARISABILITY'
     1     ,/1X,'              FXX       FXY       FYY       FXZ',
     1'        FYZ       FZZ')
C
      RETURN
      END
C*MODULE DMULTI  *DECK GRIDUSE
      SUBROUTINE GRIDUSE(DENS,Q,WORK,X,Y,Z,WT,IPT,NATOM,ATCHRG,ATDIP,
     *                   TATQUAD,TOTCHG,DIPX,DIPY,DIPZ,TQUAD,SOCT,
     *                   SHEX,ATOCT,ATHEX,L1,L2,NATM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION Q(IPT,L1),X(IPT),Y(IPT),Z(IPT),WT(IPT),DENS(L2),
     *          WORK(*),ATCHRG(NATM),ATDIP(3,NATM),TATQUAD(3,3,NATM),
     *          SOCT(7),SHEX(9),ATOCT(7,NATM),ATHEX(9,NATM)
C
C     (C) R D AMOS, CAMBRIDGE , 1995
C
      CALL RHOMAKE(Q,X,Y,Z,WORK,WORK(IPT*2+1),WORK(IPT*4+1),
     *             WORK(IPT*6+1),WORK(IPT*8+1),IPT,L1)
      I1=1
      I2=I1+IPT
      I3=I2+IPT*L1
      CALL RHOUSE(DENS,Q,X,Y,Z,WT,WORK(I1),WORK(I2),WORK(I3),
     *            IPT,L1,L2,NATM,NATOM,ATCHRG,ATDIP,TATQUAD,TOTCHG,
     *            DIPX,DIPY,DIPZ,TQUAD,C,SOCT,SHEX,ATOCT,ATHEX)
      RETURN
      END
C*MODULE DMULTI  *DECK RHOUSE
      SUBROUTINE RHOUSE(DENS,Q,X,Y,Z,WT,RHOAV,TEMP,DENSM2,IPT,L1,L2,
     *                  NATM,NATOM,ATCHRG,ATDIP,TATQUAD,TOTCHG,
     *                  DIPX,DIPY,DIPZ,TQUAD,C,SOCT,SHEX,ATOCT,ATHEX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DENS(L2),Q(IPT,L1),X(IPT),Y(IPT),Z(IPT),WT(IPT),
     *          RHOAV(IPT),TEMP(IPT,L1),DENSM2(L1,L1),
     *          ATCHRG(NATM),ATDIP(3,NATM),TATQUAD(3,3,NATM),
     *          TQUAD(3,3),C(3,NATM),SOCT(7),SHEX(9),
     *          ATOCT(7,NATM),ATHEX(9,NATM)
C
      COMMON /SQRS  / SQRP375,SQR3P75,SQR15,SQRP625,SQR516,
     *                SQR1P25,SQR4P37,SQR3564,SQR8P75
C          STATEMENT FUNCTION FOR TRIANGULAR INDEXING
      IKY(I)=(I*I-I)/2
C
      CALL VCLR(RHOAV,1,IPT)
      DO II=1,L1
         DO JJ=1,II
            IJ2=IKY(II)+JJ
            DENSM2(II,JJ)=DENS(IJ2)
            DENSM2(JJ,II)=DENS(IJ2)
         ENDDO
      ENDDO
C
C     AT THIS PLACE, SOME WORK HAS BEEN DONE TO GET THE FASTEST
C     MATRIX MULTIPLICATION ROUTINE. MXMA(PRPPOP.SRC), DGEMM(BLAS.SRC)
C     AND DGEMM(COSMO.SRC) ARE ALL SLOWER THAN MRARBR(MTHLIB,SRC)
C     WHEN COMPILED STANDALONE WITH -O3 ON SGI
C
C     CALL MXMA(Q,1,IPT,DENSM2,1,L1,TEMP,1,IPT,IPT,L1,L1)
C      CALL DGEMM
C    1('N','N',IPT,L1,L1,1.0D+00,Q,IPT,DENSM2,NUM,0.0D+00,TEMP,IPT)
C      CALL DGEMUL(Q,IPT,'N',DENSM2,L1,'N',TEMP,IPT,IPT,L1,L1)
C
      CALL MRARBR(Q,IPT,IPT,L1,DENSM2,L1,L1,TEMP,IPT)
C
      CALL VCLR(RHOAV,1,IPT)
C
      NUMB = (L1/6)*6
      DO II=1,NUMB,6
         DO M = 1,IPT
        RHOAV(M)=RHOAV(M)+TEMP(M,II)*Q(M,II)+TEMP(M,II+1)*Q(M,II+1)
        RHOAV(M)=RHOAV(M)+TEMP(M,II+2)*Q(M,II+2)+TEMP(M,II+3)*Q(M,II+3)
        RHOAV(M)=RHOAV(M)+TEMP(M,II+4)*Q(M,II+4)+TEMP(M,II+5)*Q(M,II+5)
         ENDDO
      ENDDO
C
      DO II=NUMB+1,L1
         DO M = 1,IPT
        RHOAV(M)=RHOAV(M)+TEMP(M,II)*Q(M,II)
         ENDDO
      ENDDO
C
      DO M=1,IPT
         RHOAV(M)=RHOAV(M)*WT(M)
      ENDDO
      DO M=1,IPT
         TOTCHG=TOTCHG-RHOAV(M)
         DIPX=DIPX-X(M)*RHOAV(M)
         DIPY=DIPY-Y(M)*RHOAV(M)
         DIPZ=DIPZ-Z(M)*RHOAV(M)
C
C MAKE SECOND MOMENT TENSOR TRACELESS
C
         TRACE=0.5D+00*RHOAV(M)*(X(M)**2+Y(M)**2+Z(M)**2)
         TQUAD(1,1)=TQUAD(1,1)-1.5D+00*RHOAV(M)*X(M)*X(M)+TRACE
         TQUAD(1,2)=TQUAD(1,2)-1.5D+00*RHOAV(M)*X(M)*Y(M)
         TQUAD(1,3)=TQUAD(1,3)-1.5D+00*RHOAV(M)*X(M)*Z(M)
         TQUAD(2,2)=TQUAD(2,2)-1.5D+00*RHOAV(M)*Y(M)*Y(M)+TRACE
         TQUAD(2,3)=TQUAD(2,3)-1.5D+00*RHOAV(M)*Y(M)*Z(M)
         TQUAD(3,3)=TQUAD(3,3)-1.5D+00*RHOAV(M)*Z(M)*Z(M)+TRACE
C
         XX=X(M)*X(M)*RHOAV(M)
         YY=Y(M)*Y(M)*RHOAV(M)
         ZZ=Z(M)*Z(M)*RHOAV(M)
         RR=XX+YY+ZZ
         SOCT(1)=SOCT(1)-2.5D+00*Z(M)*ZZ+1.5D+00*Z(M)*RR
         SOCT(2)=SOCT(2)-SQRP375*X(M)*(5.0D+00*ZZ-RR)
         SOCT(3)=SOCT(3)-SQRP375*Y(M)*(5.0D+00*ZZ-RR)
         SOCT(4)=SOCT(4)-SQR3P75*Z(M)*(XX-YY)
         SOCT(5)=SOCT(5)-SQR15*X(M)*Y(M)*Z(M)*RHOAV(M)
         SOCT(6)=SOCT(6)-SQRP625*X(M)*(XX-3.0D+00*YY)
         SOCT(7)=SOCT(7)-SQRP625*Y(M)*(3.0D+00*XX-YY)
         XX=X(M)*X(M)
         YY=Y(M)*Y(M)
         ZZ=Z(M)*Z(M)
         RR=XX+YY+ZZ
         SHEX(1)=SHEX(1)-0.125D+00*(35.0D+00*ZZ*ZZ-30.0D+00*ZZ*RR
     1   +3.0D+00*RR*RR)*RHOAV(M)
         SHEX(2)=SHEX(2)-SQRP625*X(M)*Z(M)*(7.0D+00*ZZ
     1  -3.0D+00*RR)*RHOAV(M)
         SHEX(3)=SHEX(3)-SQRP625*Y(M)*Z(M)*(7.0D+00*ZZ
     1  -3.0D+00*RR)*RHOAV(M)
         SHEX(4)=SHEX(4)-SQR516*(XX-YY)*(7.0D+00*ZZ-RR)
     1  *RHOAV(M)
         SHEX(5)=SHEX(5)-SQR1P25*X(M)*Y(M)*(7.0D+00*ZZ-RR)
     1  *RHOAV(M)
         SHEX(6)=SHEX(6)-SQR4P37*X(M)*Z(M)*(XX-3.0D+00*YY)
     1  *RHOAV(M)
         SHEX(7)=SHEX(7)-SQR4P37*Y(M)*Z(M)*(3.0D+00*XX-YY)
     1  *RHOAV(M)
         SHEX(8)=SHEX(8)-SQR3564*(XX*XX-6.0D+00*XX*YY
     1  +YY*YY)*RHOAV(M)
         SHEX(9)=SHEX(9)-SQR8P75*X(M)*Y(M)*(XX-YY)
     1  *RHOAV(M)
C
C ATOMIC PROPERTIES
C
         ATCHRG(NATOM)=ATCHRG(NATOM)-RHOAV(M)
         ATDIP(1,NATOM)=ATDIP(1,NATOM)-(X(M)-C(1,NATOM))*RHOAV(M)
         ATDIP(2,NATOM)=ATDIP(2,NATOM)-(Y(M)-C(2,NATOM))*RHOAV(M)
         ATDIP(3,NATOM)=ATDIP(3,NATOM)-(Z(M)-C(3,NATOM))*RHOAV(M)
         XAT=X(M)-C(1,NATOM)
         YAT=Y(M)-C(2,NATOM)
         ZAT=Z(M)-C(3,NATOM)
C
C MAKE SECOND MOMENT TENSOR TRACELESS
C
         TATRACE=0.5D+00*RHOAV(M)*(XAT**2+YAT**2+ZAT**2)
         TATQUAD(1,1,NATOM)=TATQUAD(1,1,NATOM)-
     *1.5D+00*XAT*XAT*RHOAV(M)
         TATQUAD(1,2,NATOM)=TATQUAD(1,2,NATOM)-
     *1.5D+00*XAT*YAT*RHOAV(M)
         TATQUAD(1,3,NATOM)=TATQUAD(1,3,NATOM)-
     *1.5D+00*XAT*ZAT*RHOAV(M)
         TATQUAD(2,2,NATOM)=TATQUAD(2,2,NATOM)-
     *1.5D+00*YAT*YAT*RHOAV(M)
         TATQUAD(2,3,NATOM)=TATQUAD(2,3,NATOM)-
     *1.5D+00*YAT*ZAT*RHOAV(M)
         TATQUAD(3,3,NATOM)=TATQUAD(3,3,NATOM)-
     *1.5D+00*ZAT*ZAT*RHOAV(M)
         TATQUAD(1,1,NATOM)=TATQUAD(1,1,NATOM)+TATRACE
         TATQUAD(2,2,NATOM)=TATQUAD(2,2,NATOM)+TATRACE
         TATQUAD(3,3,NATOM)=TATQUAD(3,3,NATOM)+TATRACE
C
C     OCTOPOLES (SPHERICAL HARMONIC FORM ONLY)
C
         XX=XAT*XAT*RHOAV(M)
         YY=YAT*YAT*RHOAV(M)
         ZZ=ZAT*ZAT*RHOAV(M)
         RR=XX+YY+ZZ
         ATOCT(1,NATOM)=ATOCT(1,NATOM)-
     *2.5D+00*ZAT*ZZ+1.5D+00*ZAT*RR
         ATOCT(2,NATOM)=ATOCT(2,NATOM)-
     *SQRP375*XAT*(5.0D+00*ZZ-RR)
         ATOCT(3,NATOM)=ATOCT(3,NATOM)-
     *SQRP375*YAT*(5.0D+00*ZZ-RR)
         ATOCT(4,NATOM)=ATOCT(4,NATOM)-
     *SQR3P75*ZAT*(XX-YY)
         ATOCT(5,NATOM)=ATOCT(5,NATOM)-
     *SQR15*XAT*YAT*ZAT*RHOAV(M)
         ATOCT(6,NATOM)=ATOCT(6,NATOM)-
     *SQRP625*XAT*(XX-3.0D+00*YY)
         ATOCT(7,NATOM)=ATOCT(7,NATOM)-
     *SQRP625*YAT*(3.0D+00*XX-YY)
C
C       DEBUG
C
C         WRITE(78,*)"XX,YY,ZZ,M,RHOAV,RR:"
C         WRITE(78,*)XX,YY,ZZ,M,RHOAV,RR
C         WRITE(78,*)"OCTOPOLES, AS COMPUTED IN RHOUSE:"
C         WRITE(78,*)(ATOCT(J,NATOM),J=1,7)
C
C     HEXADECAPOLES  (SPHERICAL HERMONIC FORM ONLY)
C
         XX=XAT*XAT
         YY=YAT*YAT
         ZZ=ZAT*ZAT
         RR=XX+YY+ZZ
         ATHEX(1,NATOM)=ATHEX(1,NATOM)-0.125D+00*(35.0D+00*ZZ*ZZ
     1   -30.0D+00*ZZ*RR + 3.0D+00*RR*RR)*RHOAV(M)
         ATHEX(2,NATOM)=ATHEX(2,NATOM)-SQRP625*XAT*ZAT*(7.0D+00*ZZ
     1  -3.0D+00*RR)*RHOAV(M)
         ATHEX(3,NATOM)=ATHEX(3,NATOM)-SQRP625*YAT*ZAT*(7.0D+00*ZZ
     1  -3.0D+00*RR)*RHOAV(M)
         ATHEX(4,NATOM)=ATHEX(4,NATOM)-SQR516*(XX-YY)*
     1   (7.0D+00*ZZ-RR)*RHOAV(M)
         ATHEX(5,NATOM)=ATHEX(5,NATOM)-SQR1P25*XAT*YAT*
     1   (7.0D+00*ZZ-RR)*RHOAV(M)
         ATHEX(6,NATOM)=ATHEX(6,NATOM)-SQR4P37*XAT*ZAT*
     1   (XX-3.0D+00*YY)*RHOAV(M)
         ATHEX(7,NATOM)=ATHEX(7,NATOM)-SQR4P37*YAT*ZAT*
     1   (3.0D+00*XX-YY)*RHOAV(M)
         ATHEX(8,NATOM)=ATHEX(8,NATOM)-SQR3564*
     1   (XX*XX-6.0D+00*XX*YY+YY*YY)*RHOAV(M)
         ATHEX(9,NATOM)=ATHEX(9,NATOM)-SQR8P75*XAT*YAT*(XX-YY)
     1  *RHOAV(M)
C
      ENDDO
      RETURN
      END
C*MODULE DMULTI  *DECK RHOMAKE
      SUBROUTINE RHOMAKE(P,X,Y,Z,DX,DY,DZ,AIGG,R2,NPT,L1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXGTOT=20000,MXSH=5000)
C
      DIMENSION P(NPT,L1),X(NPT),Y(NPT),Z(NPT),
     *          DX(NPT,2),DY(NPT,2),DZ(NPT,2),AIGG(NPT),R2(NPT)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
C
      DSQRT3=SQRT(3.0D+00)
      DSQRT5=SQRT(5.0D+00)
      DSQR15=SQRT(15.0D+00)
C
      L=0
      DO 501 II=1,NSHELL
      I=KATOM(II)
      XI=C(1,I)
      YI=C(2,I)
      ZI=C(3,I)
      I1=KSTART(II)
      I2=I1+KNG(II)-1
      MINK=KMIN(II)
      DO 400 M=1,NPT
      DX(M,1)=X(M)-XI
      DY(M,1)=Y(M)-YI
      DZ(M,1)=Z(M)-ZI
      DX(M,2)=DX(M,1)*DX(M,1)
      DY(M,2)=DY(M,1)*DY(M,1)
      DZ(M,2)=DZ(M,1)*DZ(M,1)
      R2(M)=DX(M,2)+DY(M,2)+DZ(M,2)
400   CONTINUE
      LL=MINK
      IF (MINK.EQ.5) LL=3
      IF (MINK.EQ.11) LL=4
      GO TO (521,522,523,524),LL
C
C......S BASIS FUNCTIONS
C
521   CONTINUE
      CALL VCLR(AIGG,1,NPT)
C$DIR SCALAR
        DO 432 IG=I1,I2
        DO 620 M=1,NPT
        AIGG(M)=AIGG(M)+CS(IG)*EXP(-EX(IG)*R2(M))
620     CONTINUE
432     CONTINUE
C
      L=L+1
C
      DO 600 M=1,NPT
      P(M,L)=AIGG(M)
600   CONTINUE
      IF (KMAX(II).EQ.1) GOTO 501
C
C....P BASIS FUNCTIONS
C
522   CONTINUE
      CALL VCLR(AIGG,1,NPT)
C$DIR SCALAR
        DO 442 IG=I1,I2
        DO 540 M=1,NPT
        AIGG(M)=AIGG(M)+CP(IG)*EXP(-EX(IG)*R2(M))
540     CONTINUE
442     CONTINUE
      DO 610 M=1,NPT
      P(M,L+1)=DX(M,1)*AIGG(M)
      P(M,L+2)=DY(M,1)*AIGG(M)
      P(M,L+3)=DZ(M,1)*AIGG(M)
610   CONTINUE
C
      L=L+3
      GO TO 501
C
C........D BASIS FUNCTIONS
C
523   CONTINUE
      CALL VCLR(AIGG,1,NPT)
C$DIR SCALAR
        DO 452 IG=I1,I2
        DO 660 M=1,NPT
        AIGG(M)=AIGG(M)+CD(IG)*EXP(-EX(IG)*R2(M))
660     CONTINUE
452     CONTINUE
      DO 621 M=1,NPT
      AIGGS=AIGG(M)*DSQRT3
      P(M,L+1)=DX(M,2)*AIGG(M)
      P(M,L+2)=DY(M,2)*AIGG(M)
      P(M,L+3)=DZ(M,2)*AIGG(M)
      P(M,L+4)=DX(M,1)*DY(M,1)*AIGGS
      P(M,L+5)=DX(M,1)*DZ(M,1)*AIGGS
      P(M,L+6)=DZ(M,1)*DY(M,1)*AIGGS
621   CONTINUE
      L=L+6
      GO TO 501
C
C......F BASIS FUNCTIONS
C
524   CONTINUE
      CALL VCLR(AIGG,1,NPT)
C$DIR SCALAR
        DO 462 IG=I1,I2
        DO 680 M=1,NPT
        AIGG(M)=AIGG(M)+CF(IG)*EXP(-EX(IG)*R2(M))
680     CONTINUE
462     CONTINUE
      DO 630 M=1,NPT
      DX3=DX(M,2)*DX(M,1)
      DY3=DY(M,2)*DY(M,1)
      DZ3=DZ(M,2)*DZ(M,1)
      DXYZ=DX(M,1)*DY(M,1)*DZ(M,1)
      AIGGS=AIGG(M)*DSQRT5
      AIGG15=AIGG(M)*DSQR15
      P(M,L+1)=DX3*AIGG(M)
      P(M,L+2)=DY3*AIGG(M)
      P(M,L+3)=DZ3*AIGG(M)
      P(M,L+4)=DX(M,2)*DY(M,1)*AIGGS
      P(M,L+5)=DX(M,2)*DZ(M,1)*AIGGS
      P(M,L+9)=DZ(M,2)*DY(M,1)*AIGGS
      P(M,L+6)=DX(M,1)*DY(M,2)*AIGGS
      P(M,L+8)=DX(M,1)*DZ(M,2)*AIGGS
      P(M,L+7)=DZ(M,1)*DY(M,2)*AIGGS
      P(M,L+10)=AIGG15*DXYZ
630   CONTINUE
      L=L+10
 501  CONTINUE
C
      RETURN
      END
C*MODULE DMULTI  *DECK SYMDF2
      SUBROUTINE SYMDF2(NTHETA,NTMAX,NPHI,NPMAX,NATOM,
     *                  IPHIS,IPHIDO,PHIZER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (PI=3.1415926535898D+00)
      PARAMETER (MXATM=2000, MXSH=5000)
C
      DIMENSION IPHIDO(128)
C
      LOGICAL GDFT,KHNSHM,MOBAS,DFTGRD
C
      COMMON /DFTCWM/ GDFT,IQUAD,NRHX,NROX,NRSX,NTMAXX,NPHIX,IDFSYM,
     *                KHNSHM,MOBAS,DFTGRD,KSFUNC(6)
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C****** WRITTEN BY C. W. MURRAY
C
C****** SYMDF DECIDES ON THE NUMBER OF ANGULAR POINTS IN THE QUADRATURE
C****** UP TO MAXIMUM DICTATED BY THE USER. FEWER ANGULAR POINTS ARE
C****** USED NEAR THE NUCLEI. IT ALSO ELIMINATES POINTS
C****** THAT ARE NOT SYMMETRY UNIQUE. SOME CODE RIPPED OFF FROM CADPAC
C
      SQ2=SQRT(2.0D+00)
C
C***** FIND THE VALUE OF NTHETA. IF ITHES=2 ADDITIONAL RESTRICTION IS
C***** THAT NTHETA IS ODD
C
      IF (NTHETA.GE.NTMAX) THEN
         NTHETA=NTMAX
         GOTO 36
      ENDIF
C
C     AMEND NTHETA TO ONE OF SEVERAL ALLOWED VALUES (FOR NAG ROUTINE)
C
      IF (NTHETA.LE.6) GOTO 36
      IF (NTHETA.EQ.7) NTHETA=8
      IF (NTHETA.EQ.9) NTHETA=10
      IF (NTHETA.EQ.11) NTHETA=12
      IF (NTHETA.EQ.13) NTHETA=14
      IF (NTHETA.EQ.15) NTHETA=16
      IF (NTHETA.LE.16) GOTO 36
      IF (NTHETA.LE.20) THEN
         NTHETA=20
         GOTO 36
      ENDIF
      IF (NTHETA.LE.24) THEN
         NTHETA=24
         GOTO 36
      ENDIF
      IF (NTHETA.LE.32) THEN
         NTHETA=32
         GOTO 36
      ENDIF
      IF (NTHETA.LE.48) THEN
         NTHETA=48
         GOTO 36
      ENDIF
      NTHETA=64
36    CONTINUE
C     IF (ITHES.EQ.2) NTHETA=((NTHETA+1)/2)*2
      NTHETA=((NTHETA+1)/2)*2
C
C********** DECIDE ON THE NUMBER OF PHI POINTS AND INITIALISE THE
C********** SYMMTERY TICKER IPHIDO
C
      NPHI=(NPMAX+IPHIS-1)/IPHIS*IPHIS
      IF (IPHIS.EQ.1) NPHI=NPMAX
      IF (NPMAX.EQ.(2*NTMAX)) NPHI=(2*NTHETA+IPHIS-1)/IPHIS*IPHIS
      DO 100 I=1,NPHI
100   IPHIDO(I)=1
      IF (NT.EQ.1.OR.IPHIS.LT.2.OR.DFTGRD) RETURN
C
C****** IF SYMMETRY IS TO BE USED ON PHI DETERMINE THE ARRAY IPHIDO
C
      IF (ABS(C(2,NATOM)).LT.1.0D-05.AND.ABS(C(1,NATOM))
     1.LT.1.0D-05) GOTO 101
      PHIZER=ATAN2(C(2,NATOM),C(1,NATOM))
      IF (PHIZER.LT.0.0D+00) PHIZER=PHIZER+2.0D+00*PI
101   DO 40 I=1,NPHI
          IF (IPHIDO(I).EQ.0) GOTO 40
          PHIOLD=PI*((2.0D+00*I)/NPHI-1.0D+00)
          PHIOLD=PHIOLD+PHIZER
          THEOLD=PI/SQ2
          ROLD=0.2D+00
          XOLD=ROLD*SIN(THEOLD)*COS(PHIOLD)
          YOLD=ROLD*SIN(THEOLD)*SIN(PHIOLD)
          ZOLD=ROLD*COS(THEOLD)
          XOLD=XOLD+C(1,NATOM)
          YOLD=YOLD+C(2,NATOM)
          ZOLD=ZOLD+C(3,NATOM)
C
         XS=U1*(XOLD-X0)+U2*(YOLD-Y0)+U3*(ZOLD-Z0)
         YS=V1*(XOLD-X0)+V2*(YOLD-Y0)+V3*(ZOLD-Z0)
         ZS=W1*(XOLD-X0)+W2*(YOLD-Y0)+W3*(ZOLD-Z0)
C
C  LOOP OVER SYMMETRY OPERATIONS AND GENERATE AND STORE RELEVANT POINTS
C
          DO 570 IT=2,NT
             NN=9*(IT-1)
             XNEW=XS*T(NN+1)+YS*T(NN+2)+ZS*T(NN+3)
             YNEW=XS*T(NN+4)+YS*T(NN+5)+ZS*T(NN+6)
             ZNEW=XS*T(NN+7)+YS*T(NN+8)+ZS*T(NN+9)
             XP=X0+U1*XNEW+V1*YNEW+W1*ZNEW
             YP=Y0+U2*XNEW+V2*YNEW+W2*ZNEW
             ZP=Z0+U3*XNEW+V3*YNEW+W3*ZNEW
             XP=XP-C(1,NATOM)
             YP=YP-C(2,NATOM)
             ZP=ZP-C(3,NATOM)
             RNEW=SQRT(XP*XP+YP*YP+ZP*ZP)
             IF (ABS(RNEW-ROLD).GT.1.0D-08) GOTO 570
             THENEW=ACOS(ZP/RNEW)
             PHINEW=ATAN2(YP,XP)
             IF (ABS(THENEW-THEOLD).GT.1.0D-08) GOTO 570
             PHITES=PHINEW-PHIZER
             IF (PHITES.LE.(-1.0D+00*PI+1.0D-08))
     *              PHITES=PHITES+2.0D+00*PI
             INEW=NINT((PHITES/PI+1.0D+00)*NPHI/2.0D+00)
             IF (INEW.EQ.I) THEN
C******** THIS POINT IS UNAFFECTED BY THIS SYMMETRY OPERATION, IT MUST
C******** THEREFORE BE ONE OF THE INTEGRATION LIMITS, AND MAY NEED TO
C******** TREATED DIFFERENTLY.
                IPHIDO(I)=2
             ELSE
                IPHIDO(INEW)=0
             ENDIF
570       CONTINUE
40    CONTINUE
C
C***** REPEAT PROCEDURE TO MAKE SURE THAT THERE IS EITHER 2 OR NO POINTS
C***** WITH IPHIDO=2
      IPDO2=0
      DO 41 I=1,NPHI
          IF (IPHIDO(I).EQ.0.OR.IPHIDO(I).EQ.1) GOTO 41
          IPDO2=IPDO2+1
          IPHI=I
41    CONTINUE
      IF (IPDO2.EQ.1) IPHIDO(IPHI)=1
      IF (IPDO2.GT.2) THEN
          WRITE(6,*) 'SOMETHING BAD HAPPENED IN NUMPROP:SYMDF2'
          CALL ABRT
      ENDIF
      RETURN
      END
C*MODULE DMULTI  *DECK LEGWT
      SUBROUTINE LEGWT(NTHETA,THETWT,THETA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION THETWT(64),THETA(64)
      DIMENSION THENAG(269),TWTNAG(269),NTHPT(64)
C
      DATA (THENAG(I),I=1,24)/ 0.0D+00,
     *  0.577350269189626D+00,-0.577350269189626D+00,
     *  0.774596669241483D+00,0.0D+00,-0.774596669241483D+00,
     *   0.861136311594052D+00,0.339981043584856D+00,
     *  -0.339981043584856D+00,-0.861136311594052D+00,
     *   0.906179845938664D+00,0.538469310105683D+00,0.0D+00,
     *  -0.538469310105683D+00,-0.906179845938664D+00,
     *   0.932469514203152D+00,0.661209386466264D+00,
     *   0.238619186083197D+00,-0.238619186083197D+00,
     *  -0.661209386466264D+00,-0.932469514203152D+00,
     *   0.960289856497536D+00,0.796666477413627D+00,
     *   0.525532409916329D+00/
      DATA (THENAG(I),I=25,51)/
     *   0.183434642495650D+00,-0.183434642495650D+00,
     *  -0.525532409916329D+00,-0.796666477413627D+00,
     *  -0.960289856497536D+00,0.973906528517172D+00,
     *  0.865063366688984D+00,0.679409568299024D+00,
     *  0.433395394129247D+00,0.148874338981631D+00,
     *  -0.148874338981631D+00,-0.433395394129247D+00,
     *  -0.679409568299024D+00,-0.865063366688984D+00,
     *  -0.973906528517172D+00,0.981560634246719D+00,
     *  0.904117256370475D+00,0.769902674194305D+00,
     *  0.587317954286617D+00,0.367831498998180D+00,
     *  0.125233408511469D+00,-0.125233408511469D+00,
     *  -0.367831498998180D+00,-0.587317954286617D+00,
     *  -0.769902674194305D+00,-0.904117256370475D+00,
     *  -0.981560634246719D+00/
      DATA (THENAG(I),I=52,80)/
     *   0.986283808696812D+00,0.928434883663573D+00,
     *   0.827201315069765D+00,0.687292904811685D+00,
     *   0.515248636358154D+00,0.319112368927890D+00,
     *   0.108054948707344D+00,-0.108054948707344D+00,
     *   -0.319112368927890D+00,-0.515248636358154D+00,
     *   -0.687292904811685D+00,-0.827201315069765D+00,
     *   -0.928434883663573D+00,-0.986283808696812D+00,
     *   0.989400934991650D+00,0.944575023073233D+00,
     *   0.865631202387832D+00,0.755404408355003D+00,
     *   0.617876244402644D+00,0.458016777657227D+00,
     *   0.281603550779259D+00,9.501250983763744D-02,
     *   -9.501250983763744D-02,-0.281603550779259D+00,
     *   -0.458016777657227D+00,-0.617876244402644D+00,
     *   -0.755404408355003D+00,-0.865631202387832D+00,
     *   -0.944575023073233D+00/
      DATA (THENAG(I),I=81,107)/
     *   -0.989400934991650D+00,0.993128599185095D+00,
     *    0.963971927277914D+00,0.912234428251326D+00,
     *    0.839116971822219D+00,0.746331906460151D+00,
     *    0.636053680726515D+00,0.510867001950827D+00,
     *    0.373706088715420D+00,0.227785851141645D+00,
     *    7.652652113349733D-02,-7.652652113349733D-02,
     *    -0.227785851141645D+00,-0.373706088715420D+00,
     *    -0.510867001950827D+00,-0.636053680726515D+00,
     *    -0.746331906460151D+00,-0.839116971822219D+00,
     *    -0.912234428251326D+00,-0.963971927277914D+00,
     *    -0.993128599185095D+00,0.995187219997021D+00,
     *    0.974728555971309D+00,0.938274552002733D+00,
     *    0.886415527004401D+00,0.820001985973903D+00,
     *    0.740124191578554D+00/
      DATA (THENAG(I),I=108,137)/
     *     0.648093651936976D+00,0.545421471388839D+00,
     *     0.433793507626045D+00,0.315042679696163D+00,
     *     0.191118867473616D+00,6.405689286260562D-02,
     *     -6.405689286260562D-02,-0.191118867473616D+00,
     *     -0.315042679696163D+00,-0.433793507626045D+00,
     *     -0.545421471388839D+00,-0.648093651936976D+00,
     *     -0.740124191578554D+00,-0.820001985973903D+00,
     *     -0.886415527004401D+00,-0.938274552002733D+00,
     *     -0.974728555971309D+00,-0.995187219997021D+00,
     *     0.997263861849481D+00,0.985611511545268D+00,
     *     0.964762255587506D+00,0.934906075937740D+00,
     *     0.896321155766052D+00,0.849367613732570D+00,
     *     0.794483795967942D+00,0.732182118740290D+00,
     *     0.663044266930215D+00,0.587715757240762D+00,
     *     0.506899908932229D+00,0.421351276130635D+00/
      DATA (THENAG(I),I=138,166)/
     *     0.331868602282128D+00,0.239287362252137D+00,
     *     0.144471961582796D+00,4.830766568773830D-02,
     *     -4.830766568773830D-02,-0.144471961582796D+00,
     *     -0.239287362252137D+00,-0.331868602282128D+00,
     *     -0.421351276130635D+00,-0.506899908932229D+00,
     *     -0.587715757240762D+00,-0.663044266930215D+00,
     *     -0.732182118740290D+00,-0.794483795967942D+00,
     *     -0.849367613732570D+00,-0.896321155766052D+00,
     *     -0.934906075937740D+00,-0.964762255587506D+00,
     *     -0.985611511545268D+00,-0.997263861849481D+00,
     *      0.998771007252426D+00,0.993530172266351D+00,
     *      0.984124583722827D+00,0.970591592546247D+00,
     *      0.952987703160431D+00,0.931386690706554D+00,
     *      0.905879136715570D+00,0.876572020274248D+00,
     *      0.843588261624393D+00/
      DATA (THENAG(I),I=167,196)/
     A0.807066204029443D+00,0.767159032515740D+00,0.724034130923815D+00,
     10.677872379632664D+00,0.628867396776513D+00,0.577224726083973D+00,
     20.523160974722233D+00,0.466902904750958D+00,0.408686481990717D+00,
     30.348755886292161D+00,0.287362487355456D+00,0.224763790394689D+00,
     40.161222356068892D+00,9.700469920946269D-02,3.238017096286936D-02,
     *   -3.238017096286936D-02,-9.700469920946269D-02,
     *   -0.161222356068892D+00,-0.224763790394689D+00,
     *   -0.287362487355456D+00,-0.348755886292161D+00,
     *   -0.408686481990717D+00,-0.466902904750958D+00,
     *   -0.523160974722233D+00,-0.577224726083973D+00,
     *   -0.628867396776513D+00,-0.677872379632664D+00,
     *   -0.724034130923815D+00,-0.767159032515740D+00,
     *   -0.807066204029443D+00/
      DATA (THENAG(I),I=197,226)/
     *   -0.843588261624393D+00,-0.876572020274248D+00,
     *   -0.905879136715570D+00,-0.931386690706554D+00,
     *   -0.952987703160431D+00,-0.970591592546247D+00,
     *   -0.984124583722827D+00,-0.993530172266351D+00,
     *   -0.998771007252426D+00,
     30.999305041735772D+00,0.996340116771955D+00,0.991013371476744D+00,
     40.983336253884626D+00,0.973326827789911D+00,0.961008799652054D+00,
     50.946411374858403D+00,0.929569172131939D+00,0.910522137078503D+00,
     60.889315445995114D+00,0.865999398154093D+00,0.840629296252580D+00,
     70.813265315122797D+00,0.783972358943341D+00,0.752819907260532D+00,
     80.719881850171611D+00,0.685236313054233D+00,0.648965471254657D+00,
     90.611155355172393D+00,0.571895646202634D+00,0.531279464019894D+00/
      DATA (THENAG(I),I=227,256)/
     *   0.489403145707053D+00,0.446366017253464D+00,
     *   0.402270157963992D+00,0.357220158337668D+00,
     *   0.311322871990211D+00,0.264687162208767D+00,
     *   0.217423643740007D+00,0.169644420423993D+00,
     *   0.121462819296121D+00, 7.299312178779903D-02,
     *   2.435029266342443D-02,-2.435029266342443D-02,
     *  -7.299312178779903D-02,-0.121462819296121D+00,
     *  -0.169644420423993D+00,-0.217423643740007D+00,
     *  -0.264687162208767D+00,-0.311322871990211D+00,
     *  -0.357220158337668D+00,-0.402270157963992D+00,
     *  -0.446366017253464D+00,-0.489403145707053D+00,
     *  -0.531279464019894D+00,-0.571895646202634D+00,
     *  -0.611155355172393D+00,-0.648965471254657D+00,
     *  -0.685236313054233D+00,-0.719881850171611D+00,
     *  -0.752819907260532D+00,-0.783972358943341D+00/
      DATA (THENAG(I),I=257,269)/
     *  -0.813265315122797D+00,-0.840629296252580D+00,
     *  -0.865999398154093D+00,-0.889315445995114D+00,
     *  -0.910522137078503D+00,-0.929569172131939D+00,
     *  -0.946411374858403D+00,-0.961008799652054D+00,
     *  -0.973326827789911D+00,-0.983336253884626D+00,
     *  -0.991013371476744D+00,-0.996340116771955D+00,
     *  -0.999305041735772D+00/
      DATA (TWTNAG(I),I=1,24)/
     A2.00000000000000D+00,
     11.00000000000000D+00,1.00000000000000D+00,
     20.555555555555556D+00,0.888888888888889D+00,0.555555555555556D+00,
     30.347854845137454D+00,0.652145154862546D+00,0.652145154862546D+00,
     40.347854845137454D+00,
     50.236926885056189D+00,0.478628670499366D+00,0.568888888888889D+00,
     60.478628670499366D+00,0.236926885056189D+00,
     70.171324492379170D+00,0.360761573048139D+00,0.467913934572691D+00,
     80.467913934572691D+00,0.360761573048139D+00,0.171324492379170D+00,
     90.101228536290376D+00,0.222381034453374D+00,0.313706645877887D+00/
      DATA (TWTNAG(I),I=25,51)/
     A0.362683783378362D+00,0.362683783378362D+00,0.313706645877887D+00,
     10.222381034453374D+00,0.101228536290376D+00,
     26.667134430868813D-02,0.149451349150581D+00,0.219086362515982D+00,
     30.269266719309996D+00,0.295524224714753D+00,0.295524224714753D+00,
     40.269266719309996D+00,0.219086362515982D+00,0.149451349150581D+00,
     56.667134430868813D-02,
     64.717533638651182D-02,0.106939325995318D+00,0.160078328543346D+00,
     70.203167426723066D+00,0.233492536538355D+00,0.249147045813403D+00,
     80.249147045813403D+00,0.233492536538355D+00,0.203167426723066D+00,
     90.160078328543346D+00,0.106939325995318D+00,4.717533638651182D-02/
      DATA (TWTNAG(I),I=52,80)/
     A3.511946033175186D-02,8.015808715976019D-02,0.121518570687903D+00,
     10.157203167158193D+00,0.185538397477938D+00,0.205198463721296D+00,
     20.215263853463158D+00,0.215263853463158D+00,0.205198463721296D+00,
     30.185538397477938D+00,0.157203167158193D+00,0.121518570687903D+00,
     48.015808715976019D-02,3.511946033175186D-02,
     52.715245941175409D-02,6.225352393864789D-02,9.515851168249277D-02,
     60.124628971255534D+00,0.149595988816577D+00,0.169156519395003D+00,
     70.182603415044924D+00,0.189450610455069D+00,0.189450610455069D+00,
     80.182603415044924D+00,0.169156519395003D+00,0.149595988816577D+00,
     90.124628971255534D+00,9.515851168249277D-02,6.225352393864789D-02/
      DATA (TWTNAG(I),I=81,107)/
     A2.715245941175409D-02,
     11.761400713915212D-02,4.060142980038694D-02,6.267204833410906D-02,
     28.327674157670474D-02,0.101930119817240D+00,0.118194531961518D+00,
     30.131688638449177D+00,0.142096109318382D+00,0.149172986472604D+00,
     40.152753387130726D+00,0.152753387130726D+00,0.149172986472604D+00,
     50.142096109318382D+00,0.131688638449177D+00,0.118194531961518D+00,
     60.101930119817240D+00,8.327674157670474D-02,6.267204833410906D-02,
     74.060142980038694D-02,1.761400713915212D-02,
     81.234122979998720D-02,2.853138862893366D-02,4.427743881741980D-02,
     95.929858491543678D-02,7.334648141108029D-02,8.619016153195327D-02/
      DATA (TWTNAG(I),I=108,137)/
     A9.761865210411388D-02,0.107444270115966D+00,0.115505668053726D+00,
     10.121670472927803D+00,0.125837456346828D+00,0.127938195346752D+00,
     20.127938195346752D+00,0.125837456346828D+00,0.121670472927803D+00,
     30.115505668053726D+00,0.107444270115966D+00,9.761865210411388D-02,
     48.619016153195327D-02,7.334648141108029D-02,5.929858491543678D-02,
     54.427743881741980D-02,2.853138862893366D-02,1.234122979998720D-02,
     67.018610009470096D-03,1.627439473090567D-02,2.539206530926205D-02,
     73.427386291302143D-02,4.283589802222667D-02,5.099805926237617D-02,
     85.868409347853554D-02,6.582222277636183D-02,7.234579410884849D-02,
     97.819389578707031D-02,8.331192422694674D-02,8.765209300440380D-02/
      DATA (TWTNAG(I),I=138,166)/
     A9.117387869576389D-02,9.384439908080456D-02,9.563872007927484D-02,
     19.654008851472779D-02,9.654008851472779D-02,9.563872007927484D-02,
     29.384439908080456D-02,9.117387869576389D-02,8.765209300440380D-02,
     38.331192422694674D-02,7.819389578707031D-02,7.234579410884849D-02,
     46.582222277636183D-02,5.868409347853554D-02,5.099805926237617D-02,
     54.283589802222667D-02,3.427386291302143D-02,2.539206530926205D-02,
     61.627439473090567D-02,7.018610009470096D-03,
     73.153346052305838D-03,7.327553901276261D-03,1.147723457923454D-02,
     81.557931572294385D-02,1.961616045735553D-02,2.357076083932438D-02,
     92.742650970835694D-02,3.116722783279808D-02,3.477722256477044D-02/
      DATA (TWTNAG(I),I=167,196)/
     A3.824135106583070D-02,4.154508294346475D-02,4.467456085669427D-02,
     14.761665849249046D-02,5.035903555385447D-02,5.289018948519367D-02,
     25.519950369998416D-02,5.727729210040321D-02,5.911483969839563D-02,
     36.070443916589388D-02,6.203942315989265D-02,6.311419228625402D-02,
     46.392423858464818D-02,6.446616443595008D-02,6.473769681268392D-02,
     56.473769681268392D-02,6.446616443595008D-02,6.392423858464818D-02,
     66.311419228625402D-02,6.203942315989265D-02,6.070443916589388D-02,
     75.911483969839563D-02,5.727729210040321D-02,5.519950369998416D-02,
     85.289018948519367D-02,5.035903555385447D-02,4.761665849249046D-02,
     94.467456085669427D-02,4.154508294346475D-02,3.824135106583070D-02/
      DATA (TWTNAG(I),I=197,226)/
     A3.477722256477044D-02,3.116722783279808D-02,2.742650970835694D-02,
     12.357076083932438D-02,1.961616045735553D-02,1.557931572294385D-02,
     21.147723457923454D-02,7.327553901276261D-03,3.153346052305838D-03,
     31.783280721696433D-03,4.147033260562467D-03,6.504457968978363D-03,
     48.846759826363947D-03,1.116813946013113D-02,1.346304789671864D-02,
     51.572603047602471D-02,1.795171577569734D-02,2.013482315353021D-02,
     62.227017380838325D-02,2.435270256871087D-02,2.637746971505466D-02,
     72.833967261425948D-02,3.023465707240248D-02,3.205792835485155D-02,
     83.380516183714160D-02,3.547221325688238D-02,3.705512854024004D-02,
     93.855015317861562D-02,3.995374113272033D-02,4.126256324262352D-02/
      DATA (TWTNAG(I),I=227,256)/
     A4.247351512365359D-02,4.358372452932345D-02,4.459055816375656D-02,
     14.549162792741814D-02,4.628479658131441D-02,4.696818281621001D-02,
     24.754016571483030D-02,4.799938859645830D-02,4.834476223480296D-02,
     34.857546744150342D-02,4.869095700913971D-02,4.869095700913971D-02,
     44.857546744150342D-02,4.834476223480296D-02,4.799938859645830D-02,
     54.754016571483030D-02,4.696818281621001D-02,4.628479658131441D-02,
     64.549162792741814D-02,4.459055816375656D-02,4.358372452932345D-02,
     74.247351512365359D-02,4.126256324262352D-02,3.995374113272033D-02,
     83.855015317861562D-02,3.705512854024004D-02,3.547221325688238D-02,
     93.380516183714160D-02,3.205792835485155D-02,3.023465707240248D-02/
      DATA (TWTNAG(I),I=257,269)/
     A2.833967261425948D-02,2.637746971505466D-02,2.435270256871087D-02,
     12.227017380838325D-02,2.013482315353021D-02,1.795171577569734D-02,
     21.572603047602471D-02,1.346304789671864D-02,1.116813946013113D-02,
     38.846759826363947D-03,6.504457968978363D-03,4.147033260562467D-03,
     41.783280721696433D-03/
      DATA NTHPT/
     A 1,2,4,7,11,16,0,22,0,30,0,40,
     1 0,52,0,66,0,0,0,82,0,0,0,102,
     2 0,0,0,0,0,0,0,126,0,0,0,0,
     3 0,0,0,0,0,0,0,0,0,0,0,158,
     4 0,0,0,0,0,0,0,0,0,0,0,0,
     5 0,0,0,206/
C
C.....JUST TO STOP THE PROGRAM CRASHING...
C
    5 CONTINUE
      NP=NTHPT(NTHETA)
      IF (NP.EQ.0.AND.NTHETA.NE.1) THEN
      NTHETA=NTHETA-1
      GOTO 5
      ENDIF
C
      IF (NP.EQ.0) THEN
         WRITE(6,*) 'DISALLOWED NUMBER OF THETA POINTS'
         CALL ABRT
      END IF
C
      DO 10 I=1,NTHETA
      THETA(I)=THENAG(NP)
      THETWT(I)=TWTNAG(NP)
      NP=NP+1
10    CONTINUE
      RETURN
      END
C*MODULE DMULTI  *DECK WEIGHT
      SUBROUTINE WEIGHT(RJ,PSMU,RIJ,AIJ,APASC,XPASC,X,Y,Z,WT,
     *                  ACCUM,XMUIJN,XMUIJ2,
     *                  IVECL,IPT,NTRANS,NATOM,NATM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      DIMENSION RJ(IVECL,NATM),PSMU(IVECL,NATM),
     *          RIJ(NATM*(NATM-1)/2),AIJ(NATM*(NATM-1)/2),
     *          XPASC(20),X(IVECL),Y(IVECL),Z(IVECL),WT(IVECL),
     *          ACCUM(IVECL),XMUIJN(IVECL),XMUIJ2(IVECL)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
C****** WRITTEN BY C. W. MURRAY
C
C***** WEIGHT CALCULATES THE WEIGHTS ASSOCIATED WITH
C***** THE BECKE PARTITIONING AMONGST THE ATOMS.
C
C******* WORK OUT THE WEIGHT FUNCTION BY BECKE PARTITIONING
C
                  DO 345 INA=1,NATM
                  DO 346 M=1,IPT
                  RJ(M,INA)=SQRT(
     1                    (C(1,INA)-X(M))**2
     1                   +(C(2,INA)-Y(M))**2
     1                   +(C(3,INA)-Z(M))**2)
                  PSMU(M,INA)=1.0D+00
346                CONTINUE
345                CONTINUE
C
                  ITEMP=0
                  DO 3350 INA=1,NATM
                  DO 3350 JNA=1,INA-1
                  ITEMP=ITEMP+1
                    DO 351 M=1,IPT
                    XMU=(RJ(M,INA)-RJ(M,JNA))*RIJ(ITEMP)
                    XMUIJ=XMU+AIJ(ITEMP)*(1-XMU*XMU)
                    XMUIJ2(M)=XMUIJ*XMUIJ
                    XMUIJN(M)=XMUIJ
                    ACCUM(M)=0.0D+00
351                 CONTINUE
                    DO 3200 I=1,NTRANS+1
                       DO 3201 M=1,IPT
                       ACCUM(M)=ACCUM(M)+XPASC(I)*XMUIJN(M)
                       XMUIJN(M)=XMUIJN(M)*XMUIJ2(M)
3201                   CONTINUE
3200                 CONTINUE
                    DO 3202 M=1,IPT
                    PSMU(M,INA)=PSMU(M,INA)*(0.5D+00-APASC*ACCUM(M))
                    PSMU(M,JNA)=PSMU(M,JNA)*(0.5D+00+APASC*ACCUM(M))
3202                CONTINUE
3350                CONTINUE
C
                  CALL VCLR(ACCUM,1,IPT)
                  DO 355 INA=1,NATM
                  DO 356 M=1,IPT
                  ACCUM(M)=ACCUM(M)+PSMU(M,INA)
356               CONTINUE
355               CONTINUE
C
C******** CONTRACT ACCUM INTO PSMU
C
                  DO 2670 INA=1,NATM
                  DO 2671 M=1,IPT
                  PSMU(M,INA)=PSMU(M,INA)/ACCUM(M)
2671              CONTINUE
2670              CONTINUE
C
                  DO 357 M=1,IPT
                  WT(M)=WT(M)*PSMU(M,NATOM)
357               CONTINUE
C
                  RETURN
                  END
C*MODULE DMULTI  *DECK NFAC
      INTEGER FUNCTION NFAC(N)
      NFAC=1
      DO 10 I=1,N
10    NFAC=NFAC*I
      RETURN
      END
C*MODULE DMULTI  *DECK RG
      SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR)
C
      INTEGER N,NM,IS1,IS2,IERR,MATZ
      DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,N),FV1(N)
      INTEGER IV1(N)
C
C     THIS ROUTINE CALLS THE RECOMMENDED SEQUENCE OF
C     ROUTINES FROM THE EIGENSYSTEM ROUTINE PACKAGE (EISPACK)
C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
C     OF A REAL GENERAL MATRIX.
C
C     ON INPUT
C
C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C        DIMENSION STATEMENT.
C
C        N  IS THE ORDER OF THE MATRIX  A.
C
C        A  CONTAINS THE REAL GENERAL MATRIX.
C
C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
C
C     ON OUTPUT
C
C        WR  AND  WI  CONTAIN THE REAL AND IMAGINARY PARTS,
C        RESPECTIVELY, OF THE EIGENVALUES.  COMPLEX CONJUGATE
C        PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE
C        EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST.
C
C        Z  CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS
C        IF MATZ IS NOT ZERO.  IF THE J-TH EIGENVALUE IS REAL, THE
C        J-TH COLUMN OF  Z  CONTAINS ITS EIGENVECTOR.  IF THE J-TH
C        EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE
C        J-TH AND (J+1)-TH COLUMNS OF  Z  CONTAIN THE REAL AND
C        IMAGINARY PARTS OF ITS EIGENVECTOR.  THE CONJUGATE OF THIS
C        VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE.
C
C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR
C           AND HQR2.  THE NORMAL COMPLETION CODE IS ZERO.
C
C        IV1  AND  FV1  ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IF (N .LE. NM) GO TO 10
      IERR = 10 * N
      GO TO 50
C
   10 CALL  BALANC(NM,N,A,IS1,IS2,FV1)
      CALL  ELMHES(NM,N,IS1,IS2,A,IV1)
      IF (MATZ .NE. 0) GO TO 20
C     .......... FIND EIGENVALUES ONLY ..........
      CALL  HQR(NM,N,IS1,IS2,A,WR,WI,IERR)
      GO TO 50
C     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
   20 CALL  ELTRAN(NM,N,IS1,IS2,A,IV1,Z)
      CALL  HQR2(NM,N,IS1,IS2,A,WR,WI,Z,IERR)
      IF (IERR .NE. 0) GO TO 50
      CALL  BALBAK(NM,N,IS1,IS2,FV1,N,Z)
   50 RETURN
      END
C*MODULE DMULTI  *DECK DENSGET
      SUBROUTINE DENSGET(SCFTYP,DA,DB,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION DA(L2),DB(L2)
C
      DATA RHF/8HRHF     /
      DATA RMC,RCI/8HMCSCF   ,8HCI      /
C
C        READ TOTAL DENSITY MATRIX
C
      CALL DAREAD(IDAF,IODA,DA,L2,16,0)
C      CALL PRTRIL(DA,NUM)
C
      IF(SCFTYP.EQ.RHF .OR. SCFTYP.EQ.RMC .OR. SCFTYP.EQ.RCI) RETURN
C
      CALL DAREAD(IDAF,IODA,DB,L2,20,0)
      CALL VADD(DA,1,DB,1,DA,1,L2)
C
      RETURN
      END
C*MODULE DMULTI  *DECK F04AEF
      SUBROUTINE F04AEF(A,IA,B,IB,N,M,C,IC,IPVT,AA,IAA,
     *                  BB,IBB,IFAIL)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(IA,N),B(IB,M),C(IC,M),IPVT(N),AA(IAA,N),BB(IBB,M)
      DO 10 I=1,N
      DO 10 J=1,N
10    AA(J,I)=A(J,I)
      DO 20 I=1,M
      DO 20 J=1,N
20    BB(J,I)=B(J,I)
      CALL LU(AA,IAA,N,IPVT)
      DO 30 I=1,M
30    CALL LUS(AA,IAA,N,IPVT,C(1,I),BB(1,I))
      IFAIL=0
      RETURN
      END
C*MODULE DMULTI  *DECK BALANC
      SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE)
C
      INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
      DOUBLE PRECISION A(NM,N),SCALE(N)
      DOUBLE PRECISION C,F,G,R,S,B2,RADIX
      LOGICAL NOCONV
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE,
C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C     THIS ROUTINE BALANCES A REAL MATRIX AND ISOLATES
C     EIGENVALUES WHENEVER POSSIBLE.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        A CONTAINS THE INPUT MATRIX TO BE BALANCED.
C
C     ON OUTPUT
C
C        A CONTAINS THE BALANCED MATRIX.
C
C        LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J)
C          IS EQUAL TO ZERO IF
C           (1) I IS GREATER THAN J AND
C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
C
C        SCALE CONTAINS INFORMATION DETERMINING THE
C           PERMUTATIONS AND SCALING FACTORS USED.
C
C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
C                 = D(J,J),      J = LOW,...,IGH
C                 = P(J)         J = IGH+1,...,N.
C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
C     THEN 1 TO LOW-1.
C
C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
C
C     THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN
C     BALANC  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
C     K,L HAVE BEEN REVERSED.)
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      RADIX = 16.0D+00
C
C   THESE 3 WILL BE SET IN LOOP 120 TO THEIR CORRECT INITIAL SETTINGS
C
      IEXC=0
      M = 0
      J = 0
C
      B2 = RADIX * RADIX
      K = 1
      L = N
      GO TO 100
C     .......... IN-LINE PROCEDURE FOR ROW AND
C                COLUMN EXCHANGE ..........
   20 SCALE(M) = J
      IF (J .EQ. M) GO TO 50
C
      DO 30 I = 1, L
         F = A(I,J)
         A(I,J) = A(I,M)
         A(I,M) = F
   30 CONTINUE
C
      DO 40 I = K, N
         F = A(J,I)
         A(J,I) = A(M,I)
         A(M,I) = F
   40 CONTINUE
C
   50 GO TO (80,130), IEXC
C     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
C                AND PUSH THEM DOWN ..........
   80 IF (L .EQ. 1) GO TO 280
      L = L - 1
C     .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
  100 DO 120 JJ = 1, L
         J = L + 1 - JJ
C
         DO 110 I = 1, L
            IF (I .EQ. J) GO TO 110
            IF (A(J,I) .NE. 0.0D+00) GO TO 120
  110    CONTINUE
C
         M = L
         IEXC = 1
         GO TO 20
  120 CONTINUE
C
      GO TO 140
C     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
C                AND PUSH THEM LEFT ..........
  130 K = K + 1
C
  140 DO 170 J = K, L
C
         DO 150 I = K, L
            IF (I .EQ. J) GO TO 150
            IF (A(I,J) .NE. 0.0D+00) GO TO 170
  150    CONTINUE
C
         M = K
         IEXC = 2
         GO TO 20
  170 CONTINUE
C     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
      DO 180 I = K, L
  180 SCALE(I) = 1.0D+00
C     .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
  190 NOCONV = .FALSE.
C
      DO 270 I = K, L
         C = 0.0D+00
         R = 0.0D+00
C
         DO 200 J = K, L
            IF (J .EQ. I) GO TO 200
            C = C + ABS(A(J,I))
            R = R + ABS(A(I,J))
  200    CONTINUE
C     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
         IF (C .EQ. 0.0D+00 .OR. R .EQ. 0.0D+00) GO TO 270
         G = R / RADIX
         F = 1.0D+00
         S = C + R
  210    IF (C .GE. G) GO TO 220
         F = F * RADIX
         C = C * B2
         GO TO 210
  220    G = R * RADIX
  230    IF (C .LT. G) GO TO 240
         F = F / RADIX
         C = C / B2
         GO TO 230
C     .......... NOW BALANCE ..........
  240    IF ((C + R) / F .GE. 0.95D+00 * S) GO TO 270
         G = 1.0D+00 / F
         SCALE(I) = SCALE(I) * F
         NOCONV = .TRUE.
C
         DO 250 J = K, N
  250    A(I,J) = A(I,J) * G
C
         DO 260 J = 1, L
  260    A(J,I) = A(J,I) * F
C
  270 CONTINUE
C
      IF (NOCONV) GO TO 190
C
  280 LOW = K
      IGH = L
      RETURN
      END
C*MODULE DMULTI  *DECK ELMHES
      SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT)
C
      INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
      DOUBLE PRECISION A(NM,N)
      DOUBLE PRECISION X,Y
      INTEGER INT(IGH)
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES,
C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C     GIVEN A REAL GENERAL MATRIX, THIS ROUTINE
C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
C     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C          ROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
C          SET LOW=1, IGH=N.
C
C        A CONTAINS THE INPUT MATRIX.
C
C     ON OUTPUT
C
C        A CONTAINS THE HESSENBERG MATRIX.  THE MULTIPLIERS
C          WHICH WERE USED IN THE REDUCTION ARE STORED IN THE
C          REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
C
C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
C          INTERCHANGED IN THE REDUCTION.
C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      LA = IGH - 1
      KP1 = LOW + 1
      IF (LA .LT. KP1) GO TO 200
C
      DO 180 M = KP1, LA
         MM1 = M - 1
         X = 0.0D+00
         I = M
C
         DO 100 J = M, IGH
            IF (ABS(A(J,MM1)) .LE. ABS(X)) GO TO 100
            X = A(J,MM1)
            I = J
  100    CONTINUE
C
         INT(M) = I
         IF (I .EQ. M) GO TO 130
C     .......... INTERCHANGE ROWS AND COLUMNS OF A ..........
         DO 110 J = MM1, N
            Y = A(I,J)
            A(I,J) = A(M,J)
            A(M,J) = Y
  110    CONTINUE
C
         DO 120 J = 1, IGH
            Y = A(J,I)
            A(J,I) = A(J,M)
            A(J,M) = Y
  120    CONTINUE
C     .......... END INTERCHANGE ..........
  130    IF (X .EQ. 0.0D+00) GO TO 180
         MP1 = M + 1
C
         DO 160 I = MP1, IGH
            Y = A(I,MM1)
            IF (Y .EQ. 0.0D+00) GO TO 160
            Y = Y / X
            A(I,MM1) = Y
C
            DO 140 J = M, N
  140       A(I,J) = A(I,J) - Y * A(M,J)
C
            DO 150 J = 1, IGH
  150       A(J,M) = A(J,M) + Y * A(J,I)
C
  160    CONTINUE
C
  180 CONTINUE
C
  200 RETURN
      END
C*MODULE DMULTI  *DECK HQR
      SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
C  RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG)
C
      INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR
      DOUBLE PRECISION H(NM,N),WR(N),WI(N)
      DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2
      LOGICAL NOTLAS
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR,
C     NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
C
C     THIS ROUTINE FINDS THE EIGENVALUES OF A REAL
C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C          ROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
C          SET LOW=1, IGH=N.
C
C        H CONTAINS THE UPPER HESSENBERG MATRIX.  INFORMATION ABOUT
C          THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG
C          FORM BY  ELMHES  OR  ORTHES, IF PERFORMED, IS STORED
C          IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
C
C     ON OUTPUT
C
C        H HAS BEEN DESTROYED.  THEREFORE, IT MUST BE SAVED
C          BEFORE CALLING  HQR  IF SUBSEQUENT CALCULATION AND
C          BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED.
C
C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C          FOR INDICES IERR+1,...,N.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED SEPTEMBER 1989.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      NORM = 0.0D+00
      K = 1
C     .......... STORE ROOTS ISOLATED BY BALANC
C                AND COMPUTE MATRIX NORM ..........
      DO 50 I = 1, N
C
         DO 40 J = K, N
   40    NORM = NORM + ABS(H(I,J))
C
         K = I
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
         WR(I) = H(I,I)
         WI(I) = 0.0D+00
   50 CONTINUE
C
      EN = IGH
      T = 0.0D+00
      ITN = 30*N
C     .......... SEARCH FOR NEXT EIGENVALUES ..........
   60 IF (EN .LT. LOW) GO TO 1001
      ITS = 0
      NA = EN - 1
      ENM2 = NA - 1
C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
   70 DO 80 LL = LOW, EN
         L = EN + LOW - LL
         IF (L .EQ. LOW) GO TO 100
         S = ABS(H(L-1,L-1)) + ABS(H(L,L))
         IF (S .EQ. 0.0D+00) S = NORM
         TST1 = S
         TST2 = TST1 + ABS(H(L,L-1))
         IF (TST2 .EQ. TST1) GO TO 100
   80 CONTINUE
C     .......... FORM SHIFT ..........
  100 X = H(EN,EN)
      IF (L .EQ. EN) GO TO 270
      Y = H(NA,NA)
      W = H(EN,NA) * H(NA,EN)
      IF (L .EQ. NA) GO TO 280
      IF (ITN .EQ. 0) GO TO 1000
      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C     .......... FORM EXCEPTIONAL SHIFT ..........
      T = T + X
C
      DO 120 I = LOW, EN
  120 H(I,I) = H(I,I) - X
C
      S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
      X = 0.75D+00 * S
      Y = X
      W = -0.4375D+00 * S * S
  130 ITS = ITS + 1
      ITN = ITN - 1
C     .......... LOOK FOR TWO CONSECUTIVE SMALL
C                SUB-DIAGONAL ELEMENTS.
C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
      DO 140 MM = L, ENM2
         M = ENM2 + L - MM
         ZZ = H(M,M)
         R = X - ZZ
         S = Y - ZZ
         P = (R * S - W) / H(M+1,M) + H(M,M+1)
         Q = H(M+1,M+1) - ZZ - R - S
         R = H(M+2,M+1)
         S = ABS(P) + ABS(Q) + ABS(R)
         P = P / S
         Q = Q / S
         R = R / S
         IF (M .EQ. L) GO TO 150
         TST1 = ABS(P)*(ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))
         TST2 = TST1 + ABS(H(M,M-1))*(ABS(Q) + ABS(R))
         IF (TST2 .EQ. TST1) GO TO 150
  140 CONTINUE
C
  150 MP2 = M + 2
C
      DO 160 I = MP2, EN
         H(I,I-2) = 0.0D+00
         IF (I .EQ. MP2) GO TO 160
         H(I,I-3) = 0.0D+00
  160 CONTINUE
C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C                COLUMNS M TO EN ..........
      DO 260 K = M, NA
         NOTLAS = K .NE. NA
         IF (K .EQ. M) GO TO 170
         P = H(K,K-1)
         Q = H(K+1,K-1)
         R = 0.0D+00
         IF (NOTLAS) R = H(K+2,K-1)
         X = ABS(P) + ABS(Q) + ABS(R)
         IF (X .EQ. 0.0D+00) GO TO 260
         P = P / X
         Q = Q / X
         R = R / X
  170    S = SIGN(SQRT(P*P+Q*Q+R*R),P)
         IF (K .EQ. M) GO TO 180
         H(K,K-1) = -S * X
         GO TO 190
  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
  190    P = P + S
         X = P / S
         Y = Q / S
         ZZ = R / S
         Q = Q / P
         R = R / P
         IF (NOTLAS) GO TO 225
C     .......... ROW MODIFICATION ..........
         DO 200 J = K, EN
            P = H(K,J) + Q * H(K+1,J)
            H(K,J) = H(K,J) - P * X
            H(K+1,J) = H(K+1,J) - P * Y
  200    CONTINUE
C
         J = MIN0(EN,K+3)
C     .......... COLUMN MODIFICATION ..........
         DO 210 I = L, J
            P = X * H(I,K) + Y * H(I,K+1)
            H(I,K) = H(I,K) - P
            H(I,K+1) = H(I,K+1) - P * Q
  210    CONTINUE
         GO TO 255
  225    CONTINUE
C     .......... ROW MODIFICATION ..........
         DO 230 J = K, EN
            P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
            H(K,J) = H(K,J) - P * X
            H(K+1,J) = H(K+1,J) - P * Y
            H(K+2,J) = H(K+2,J) - P * ZZ
  230    CONTINUE
C
         J = MIN0(EN,K+3)
C     .......... COLUMN MODIFICATION ..........
         DO 240 I = L, J
            P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
            H(I,K) = H(I,K) - P
            H(I,K+1) = H(I,K+1) - P * Q
            H(I,K+2) = H(I,K+2) - P * R
  240    CONTINUE
  255    CONTINUE
C
  260 CONTINUE
C
      GO TO 70
C     .......... ONE ROOT FOUND ..........
  270 WR(EN) = X + T
      WI(EN) = 0.0D+00
      EN = NA
      GO TO 60
C     .......... TWO ROOTS FOUND ..........
  280 P = (Y - X) / 2.0D+00
      Q = P * P + W
      ZZ = SQRT(ABS(Q))
      X = X + T
      IF (Q .LT. 0.0D+00) GO TO 320
C     .......... REAL PAIR ..........
      ZZ = P + SIGN(ZZ,P)
      WR(NA) = X + ZZ
      WR(EN) = WR(NA)
      IF (ZZ .NE. 0.0D+00) WR(EN) = X - W / ZZ
      WI(NA) = 0.0D+00
      WI(EN) = 0.0D+00
      GO TO 330
C     .......... COMPLEX PAIR ..........
  320 WR(NA) = X + P
      WR(EN) = X + P
      WI(NA) = ZZ
      WI(EN) = -ZZ
  330 EN = ENM2
      GO TO 60
C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C                CONVERGED AFTER 30*N ITERATIONS ..........
 1000 IERR = EN
 1001 RETURN
      END
C*MODULE DMULTI  *DECK ELTRAN
      SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z)
C
      INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
      DOUBLE PRECISION A(NM,IGH),Z(NM,N)
      INTEGER INT(IGH)
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C     THIS ROUTINE ACCUMULATES THE STABILIZED ELEMENTARY
C     SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A
C     REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY  ELMHES.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C          ROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
C          SET LOW=1, IGH=N.
C
C        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
C          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE
C          BELOW THE SUBDIAGONAL.
C
C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
C          INTERCHANGED IN THE REDUCTION BY  ELMHES.
C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
C
C     ON OUTPUT
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  ELMHES.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
C     .......... INITIALIZE Z TO IDENTITY MATRIX ..........
      DO 80 J = 1, N
C
         DO 60 I = 1, N
   60    Z(I,J) = 0.0D+00
C
         Z(J,J) = 1.0D+00
   80 CONTINUE
C
      KL = IGH - LOW - 1
      IF (KL .LT. 1) GO TO 200
C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
      DO 140 MM = 1, KL
         MP = IGH - MM
         MP1 = MP + 1
C
         DO 100 I = MP1, IGH
  100    Z(I,MP) = A(I,MP-1)
C
         I = INT(MP)
         IF (I .EQ. MP) GO TO 140
C
         DO 130 J = MP, IGH
            Z(MP,J) = Z(I,J)
            Z(I,J) = 0.0D+00
  130    CONTINUE
C
         Z(I,MP) = 1.0D+00
  140 CONTINUE
C
  200 RETURN
      END
C*MODULE DMULTI  *DECK HQR2
      SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR)
C
      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN,
     X        IGH,ITN,ITS,LOW,MP2,ENM2,IERR
      DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N)
      DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2
      LOGICAL NOTLAS
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2,
C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C     THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD.  THE
C     EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND
C     IF  ELMHES  AND  ELTRAN  OR  ORTHES  AND  ORTRAN  HAVE
C     BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM
C     AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C          ROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
C          SET LOW=1, IGH=N.
C
C        H CONTAINS THE UPPER HESSENBERG MATRIX.
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY  ELTRAN
C          AFTER THE REDUCTION BY  ELMHES, OR BY  ORTRAN  AFTER THE
C          REDUCTION BY  ORTHES, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE
C          IDENTITY MATRIX.
C
C     ON OUTPUT
C
C        H HAS BEEN DESTROYED.
C
C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C          FOR INDICES IERR+1,...,N.
C
C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
C          IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z
C          CONTAINS ITS EIGENVECTOR.  IF THE I-TH EIGENVALUE IS COMPLEX
C          WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH
C          COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS
C          EIGENVECTOR.  THE EIGENVECTORS ARE UNNORMALIZED.  IF AN
C          ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
C
C     CALLS CDIV FOR COMPLEX DIVISION.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      NORM = 0.0D+00
      K = 1
C     .......... STORE ROOTS ISOLATED BY BALANC
C                AND COMPUTE MATRIX NORM ..........
      DO 50 I = 1, N
C
         DO 40 J = K, N
   40    NORM = NORM + ABS(H(I,J))
C
         K = I
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
         WR(I) = H(I,I)
         WI(I) = 0.0D+00
   50 CONTINUE
C
      EN = IGH
      T = 0.0D+00
      ITN = 30*N
C     .......... SEARCH FOR NEXT EIGENVALUES ..........
   60 IF (EN .LT. LOW) GO TO 340
      ITS = 0
      NA = EN - 1
      ENM2 = NA - 1
C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
   70 DO 80 LL = LOW, EN
         L = EN + LOW - LL
         IF (L .EQ. LOW) GO TO 100
         S = ABS(H(L-1,L-1)) + ABS(H(L,L))
         IF (S .EQ. 0.0D+00) S = NORM
         TST1 = S
         TST2 = TST1 + ABS(H(L,L-1))
         IF (TST2 .EQ. TST1) GO TO 100
   80 CONTINUE
C     .......... FORM SHIFT ..........
  100 X = H(EN,EN)
      IF (L .EQ. EN) GO TO 270
      Y = H(NA,NA)
      W = H(EN,NA) * H(NA,EN)
      IF (L .EQ. NA) GO TO 280
      IF (ITN .EQ. 0) GO TO 1000
      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C     .......... FORM EXCEPTIONAL SHIFT ..........
      T = T + X
C
      DO 120 I = LOW, EN
  120 H(I,I) = H(I,I) - X
C
      S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
      X = 0.75D+00 * S
      Y = X
      W = -0.4375D+00 * S * S
  130 ITS = ITS + 1
      ITN = ITN - 1
C     .......... LOOK FOR TWO CONSECUTIVE SMALL
C                SUB-DIAGONAL ELEMENTS.
C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
      DO 140 MM = L, ENM2
         M = ENM2 + L - MM
         ZZ = H(M,M)
         R = X - ZZ
         S = Y - ZZ
         P = (R * S - W) / H(M+1,M) + H(M,M+1)
         Q = H(M+1,M+1) - ZZ - R - S
         R = H(M+2,M+1)
         S = ABS(P) + ABS(Q) + ABS(R)
         P = P / S
         Q = Q / S
         R = R / S
         IF (M .EQ. L) GO TO 150
         TST1 = ABS(P)*(ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))
         TST2 = TST1 + ABS(H(M,M-1))*(ABS(Q) + ABS(R))
         IF (TST2 .EQ. TST1) GO TO 150
  140 CONTINUE
C
  150 MP2 = M + 2
C
      DO 160 I = MP2, EN
         H(I,I-2) = 0.0D+00
         IF (I .EQ. MP2) GO TO 160
         H(I,I-3) = 0.0D+00
  160 CONTINUE
C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C                COLUMNS M TO EN ..........
      DO 260 K = M, NA
         NOTLAS = K .NE. NA
         IF (K .EQ. M) GO TO 170
         P = H(K,K-1)
         Q = H(K+1,K-1)
         R = 0.0D+00
         IF (NOTLAS) R = H(K+2,K-1)
         X = ABS(P) + ABS(Q) + ABS(R)
         IF (X .EQ. 0.0D+00) GO TO 260
         P = P / X
         Q = Q / X
         R = R / X
  170    S = SIGN(SQRT(P*P+Q*Q+R*R),P)
         IF (K .EQ. M) GO TO 180
         H(K,K-1) = -S * X
         GO TO 190
  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
  190    P = P + S
         X = P / S
         Y = Q / S
         ZZ = R / S
         Q = Q / P
         R = R / P
         IF (NOTLAS) GO TO 225
C     .......... ROW MODIFICATION ..........
         DO 200 J = K, N
            P = H(K,J) + Q * H(K+1,J)
            H(K,J) = H(K,J) - P * X
            H(K+1,J) = H(K+1,J) - P * Y
  200    CONTINUE
C
         J = MIN0(EN,K+3)
C     .......... COLUMN MODIFICATION ..........
         DO 210 I = 1, J
            P = X * H(I,K) + Y * H(I,K+1)
            H(I,K) = H(I,K) - P
            H(I,K+1) = H(I,K+1) - P * Q
  210    CONTINUE
C     .......... ACCUMULATE TRANSFORMATIONS ..........
         DO 220 I = LOW, IGH
            P = X * Z(I,K) + Y * Z(I,K+1)
            Z(I,K) = Z(I,K) - P
            Z(I,K+1) = Z(I,K+1) - P * Q
  220    CONTINUE
         GO TO 255
  225    CONTINUE
C     .......... ROW MODIFICATION ..........
         DO 230 J = K, N
            P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
            H(K,J) = H(K,J) - P * X
            H(K+1,J) = H(K+1,J) - P * Y
            H(K+2,J) = H(K+2,J) - P * ZZ
  230    CONTINUE
C
         J = MIN0(EN,K+3)
C     .......... COLUMN MODIFICATION ..........
         DO 240 I = 1, J
            P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
            H(I,K) = H(I,K) - P
            H(I,K+1) = H(I,K+1) - P * Q
            H(I,K+2) = H(I,K+2) - P * R
  240    CONTINUE
C     .......... ACCUMULATE TRANSFORMATIONS ..........
         DO 250 I = LOW, IGH
            P = X * Z(I,K) + Y * Z(I,K+1) + ZZ * Z(I,K+2)
            Z(I,K) = Z(I,K) - P
            Z(I,K+1) = Z(I,K+1) - P * Q
            Z(I,K+2) = Z(I,K+2) - P * R
  250    CONTINUE
  255    CONTINUE
C
  260 CONTINUE
C
      GO TO 70
C     .......... ONE ROOT FOUND ..........
  270 H(EN,EN) = X + T
      WR(EN) = H(EN,EN)
      WI(EN) = 0.0D+00
      EN = NA
      GO TO 60
C     .......... TWO ROOTS FOUND ..........
  280 P = (Y - X) / 2.0D+00
      Q = P * P + W
      ZZ = SQRT(ABS(Q))
      H(EN,EN) = X + T
      X = H(EN,EN)
      H(NA,NA) = Y + T
      IF (Q .LT. 0.0D+00) GO TO 320
C     .......... REAL PAIR ..........
      ZZ = P + SIGN(ZZ,P)
      WR(NA) = X + ZZ
      WR(EN) = WR(NA)
      IF (ZZ .NE. 0.0D+00) WR(EN) = X - W / ZZ
      WI(NA) = 0.0D+00
      WI(EN) = 0.0D+00
      X = H(EN,NA)
      S = ABS(X) + ABS(ZZ)
      P = X / S
      Q = ZZ / S
      R = SQRT(P*P+Q*Q)
      P = P / R
      Q = Q / R
C     .......... ROW MODIFICATION ..........
      DO 290 J = NA, N
         ZZ = H(NA,J)
         H(NA,J) = Q * ZZ + P * H(EN,J)
         H(EN,J) = Q * H(EN,J) - P * ZZ
  290 CONTINUE
C     .......... COLUMN MODIFICATION ..........
      DO 300 I = 1, EN
         ZZ = H(I,NA)
         H(I,NA) = Q * ZZ + P * H(I,EN)
         H(I,EN) = Q * H(I,EN) - P * ZZ
  300 CONTINUE
C     .......... ACCUMULATE TRANSFORMATIONS ..........
      DO 310 I = LOW, IGH
         ZZ = Z(I,NA)
         Z(I,NA) = Q * ZZ + P * Z(I,EN)
         Z(I,EN) = Q * Z(I,EN) - P * ZZ
  310 CONTINUE
C
      GO TO 330
C     .......... COMPLEX PAIR ..........
  320 WR(NA) = X + P
      WR(EN) = X + P
      WI(NA) = ZZ
      WI(EN) = -ZZ
  330 EN = ENM2
      GO TO 60
C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
C                VECTORS OF UPPER TRIANGULAR FORM ..........
  340 IF (NORM .EQ. 0.0D+00) GO TO 1001
C     .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
      DO 800 NN = 1, N
         EN = N + 1 - NN
         P = WR(EN)
         Q = WI(EN)
         NA = EN - 1
         IF (Q) 710, 600, 800
C     .......... REAL VECTOR ..........
  600    M = EN
         H(EN,EN) = 1.0D+00
         IF (NA .EQ. 0) GO TO 800
C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
         DO 700 II = 1, NA
            I = EN - II
            W = H(I,I) - P
            R = 0.0D+00
C
            DO 610 J = M, EN
  610       R = R + H(I,J) * H(J,EN)
C
            IF (WI(I) .GE. 0.0D+00) GO TO 630
            ZZ = W
            S = R
            GO TO 700
  630       M = I
            IF (WI(I) .NE. 0.0D+00) GO TO 640
            T = W
            IF (T .NE. 0.0D+00) GO TO 635
               TST1 = NORM
               T = TST1
  632          T = 0.01D+00 * T
               TST2 = NORM + T
               IF (TST2 .GT. TST1) GO TO 632
  635       H(I,EN) = -R / T
            GO TO 680
C     .......... SOLVE REAL EQUATIONS ..........
  640       X = H(I,I+1)
            Y = H(I+1,I)
            Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I)
            T = (X * S - ZZ * R) / Q
            H(I,EN) = T
            IF (ABS(X) .LE. ABS(ZZ)) GO TO 650
            H(I+1,EN) = (-R - W * T) / X
            GO TO 680
  650       H(I+1,EN) = (-S - Y * T) / ZZ
C
C     .......... OVERFLOW CONTROL ..........
  680       T = ABS(H(I,EN))
            IF (T .EQ. 0.0D+00) GO TO 700
            TST1 = T
            TST2 = TST1 + 1.0D+00/TST1
            IF (TST2 .GT. TST1) GO TO 700
            DO 690 J = I, EN
               H(J,EN) = H(J,EN)/T
  690       CONTINUE
C
  700    CONTINUE
C     .......... END REAL VECTOR ..........
         GO TO 800
C     .......... COMPLEX VECTOR ..........
  710    M = NA
C     .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
C                EIGENVECTOR MATRIX IS TRIANGULAR ..........
         IF (ABS(H(EN,NA)) .LE. ABS(H(NA,EN))) GO TO 720
         H(NA,NA) = Q / H(EN,NA)
         H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA)
         GO TO 730
  720    CALL CDIV(0.0D+00,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN))
  730    H(EN,NA) = 0.0D+00
         H(EN,EN) = 1.0D+00
         ENM2 = NA - 1
         IF (ENM2 .EQ. 0) GO TO 800
C     .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- ..........
         DO 795 II = 1, ENM2
            I = NA - II
            W = H(I,I) - P
            RA = 0.0D+00
            SA = 0.0D+00
C
            DO 760 J = M, EN
               RA = RA + H(I,J) * H(J,NA)
               SA = SA + H(I,J) * H(J,EN)
  760       CONTINUE
C
            IF (WI(I) .GE. 0.0D+00) GO TO 770
            ZZ = W
            R = RA
            S = SA
            GO TO 795
  770       M = I
            IF (WI(I) .NE. 0.0D+00) GO TO 780
            CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN))
            GO TO 790
C     .......... SOLVE COMPLEX EQUATIONS ..........
  780       X = H(I,I+1)
            Y = H(I+1,I)
            VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q
            VI = (WR(I) - P) * 2.0D+00 * Q
            IF (VR .NE. 0.0D+00 .OR. VI .NE. 0.0D+00) GO TO 784
               TST1 = NORM * (ABS(W) + ABS(Q) + ABS(X)
     X                      + ABS(Y) + ABS(ZZ))
               VR = TST1
  783          VR = 0.01D+00 * VR
               TST2 = TST1 + VR
               IF (TST2 .GT. TST1) GO TO 783
  784       CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,
     X                H(I,NA),H(I,EN))
            IF (ABS(X) .LE. ABS(ZZ) + ABS(Q)) GO TO 785
            H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X
            H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X
            GO TO 790
  785       CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q,
     X                H(I+1,NA),H(I+1,EN))
C
C     .......... OVERFLOW CONTROL ..........
  790       T = MAX(ABS(H(I,NA)), ABS(H(I,EN)))
            IF (T .EQ. 0.0D+00) GO TO 795
            TST1 = T
            TST2 = TST1 + 1.0D+00/TST1
            IF (TST2 .GT. TST1) GO TO 795
            DO 792 J = I, EN
               H(J,NA) = H(J,NA)/T
               H(J,EN) = H(J,EN)/T
  792       CONTINUE
C
  795    CONTINUE
C     .......... END COMPLEX VECTOR ..........
  800 CONTINUE
C     .......... END BACK SUBSTITUTION.
C                VECTORS OF ISOLATED ROOTS ..........
      DO 840 I = 1, N
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
C
         DO 820 J = I, N
  820    Z(I,J) = H(I,J)
C
  840 CONTINUE
C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
C                VECTORS OF ORIGINAL FULL MATRIX.
C                FOR J=N STEP -1 UNTIL LOW DO -- ..........
      DO 880 JJ = LOW, N
         J = N + LOW - JJ
         M = MIN0(J,IGH)
C
         DO 880 I = LOW, IGH
            ZZ = 0.0D+00
C
            DO 860 K = LOW, M
  860       ZZ = ZZ + Z(I,K) * H(K,J)
C
            Z(I,J) = ZZ
  880 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C                CONVERGED AFTER 30*N ITERATIONS ..........
 1000 IERR = EN
 1001 RETURN
      END
C*MODULE DMULTI  *DECK BALBAK
      SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z)
      INTEGER I,J,K,M,N,II,NM,IGH,LOW
      DOUBLE PRECISION SCALE(N),Z(NM,M)
      DOUBLE PRECISION S
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK,
C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C     THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C     BALANCED MATRIX DETERMINED BY  BALANC.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        LOW AND IGH ARE INTEGERS DETERMINED BY  BALANC.
C
C        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
C          AND SCALING FACTORS USED BY  BALANC.
C
C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
C
C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
C
C     ON OUTPUT
C
C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IF (M .EQ. 0) GO TO 200
      IF (IGH .EQ. LOW) GO TO 120
C
      DO 110 I = LOW, IGH
         S = SCALE(I)
C     .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED
C                IF THE FOREGOING STATEMENT IS REPLACED BY
C                S=1.0D+00/SCALE(I). ..........
         DO 100 J = 1, M
  100    Z(I,J) = Z(I,J) * S
C
  110 CONTINUE
C     ......... FOR I=LOW-1 STEP -1 UNTIL 1,
C               IGH+1 STEP 1 UNTIL N DO -- ..........
  120 DO 140 II = 1, N
         I = II
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140
         IF (I .LT. LOW) I = LOW - II
         K = INT(SCALE(I))
         IF (K .EQ. I) GO TO 140
C
         DO 130 J = 1, M
            S = Z(I,J)
            Z(I,J) = Z(K,J)
            Z(K,J) = S
  130    CONTINUE
C
  140 CONTINUE
C
  200 RETURN
      END
C*MODULE DMULTI  *DECK CDIV
      SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI)
      DOUBLE PRECISION AR,AI,BR,BI,CR,CI
C
C     COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI)
C
      DOUBLE PRECISION S,ARS,AIS,BRS,BIS
      S = ABS(BR) + ABS(BI)
      ARS = AR/S
      AIS = AI/S
      BRS = BR/S
      BIS = BI/S
      S = BRS**2 + BIS**2
      CR = (ARS*BRS + AIS*BIS)/S
      CI = (AIS*BRS - ARS*BIS)/S
      RETURN
      END
