C 11 APR 08 - JMM - STNMOM,STONE: MAKEFP OUTPUT TO UNIT NUMBER -NFTEFP-
C  4 MAR 08 - HN,JHS - PROPTY: INCLUDE GMCCI AS MCSCF CI STEP, NEO CALL
C  7 DEC 07 - JNM - STNMOM: ADD UNITS TO COORDINATE PRINTING LABEL
C 20 AUG 07 - DGF - SYNCHRONISE FMOOPT AND FMOPNT
C 21 MAY 07 - MWS - SPIND: ZAPT MUST EXIT W/O DOING PROPERTY
C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C 22 DEC 06 - DGF - ADD ISKPPRP OPTION, SYNCH. PSILVL, FIX BUG IN MOPOP
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 25 OCT 06 - SPC - PROPTY: NEO RUNS SHOULD JUST SKIP PROPERTIES
C  7 APR 06 - MWS - REMOVE NON-ARGUMENT ADJUSTABLE ARRAYS
C 29 MAR 06 - CHC - CHANGE COMPARISONS TO FORTRAN-STYLE .EQ. AND .GT.
C 22 FEB 06 - LVS - NUMERICAL PART FOR STONE ANALYSIS ADDED
C  2 FEB 06 - MWS - SPIND: PARALLELIZATION AND H+I
C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK HERMIT,WERMIT,NSHEL,FRGINF,EDCMP
C  5 FEB 05 - MWS - PAD COMMON MCINP
C 23 JUL 04 - MWS - PROPTY: FIX LOGIC FOR STATE-AVERAGED MCSCF PROPS
C 16 JAN 04 - MWS - STNMOM,STNXYZ: FIX F, ADD G, PARALLELIZE INTS,
C                   DETSPM: READ FILE 12 ONLY ON MASTER NODE
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  4 NOV 03 - MWS - STNMOM: ADD UNITS TO THE PRINT OUT
C 16 JUN 03 - HL  - STNMOM: BUG FIX FOR PRINT OUT
C 14 JAN 03 - JI  - ADD HOOKS FOR ORMAS CI STEP
C 12 DEC 02 - MWS - STNMOM: OUTPUT CHANGES
C  7 AUG 02 - HL,PM - PROPTY: INCLUDE IMGABI METHOD
C                   STNMOM,STNRD,STONE: END POINT MULTIPOLE CORRECTION
C                   METHOD THAT GENERATE INTEGER CHARGE EFP FRAGMENTS
C                   DENSITY MODIFICATION TO GENERATE INTEGER CHARGE EFP
C 20 JUN 02 - MWS - DETSPM: CHECK RUNS SHOULD ONLY CHECK
C 22 MAY 02 - GDF - DETSPM: AVOID PARALLEL I/O ON DET CI VECTOR FILE
C 17 APR 02 - PND - STNRD: ADD MASTER ONLY TO A PRINT LINE
C 16 FEB 02 - MWS - CORRECT THE CALLING LOGIC FOR DETSPM
C 24 JAN 02 - MVP - FIX THE PROBLEM WITH MCSCF SPIN DENSITY FOR C1 GROUP
C  6 SEP 01 - CHC - MODIFY STNXYZ FOR QFMM, ADD TX4,TY4,TZ4.
C  1 AUG 01 - JI  - TWEAKS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - MAF - STNMOM: DUMP A SECOND COPY OF THE PUNCH OUTPUT
C  3 MAY 00 - PM,HL - STNRD,STONE: OVERLAPPING REGIONS, CAP CHARGES
C 19 NOV 00 - RMM,HL - STNRD: PAD OUT EDCMP COMMON FOR POLAPP KEYWORD
C 26 OCT 00 - MWS - DETSPM,SPIND: EXIT QUICKLY IF MCSCF DETERMINANTS
C 11 OCT 00 - MWS - PAD COMMON BLOCK PCMPAR
C  1 MAY 00 - MWS - DETSPM: USE PURES OPTION TO SELECT CORRECT CI STATES
C 25 MAR 00 - KKB/LNB - CALL TO NUMPROP: AMOS DMA FOR COSMO
C 16 FEB 00 - MVP - MCSCF SPIN DENSITY IF DETERMINANTS
C 15 JAN 00 - VK  - STONE,STNRD,STNMOM: FIX NUCLEAR CHARGE BUG, TYPO IN
C                   STONE, CHANGES DUE TO PRESENCE OF QM/MM BUFFER
C 21 DEC 99 - MWS - MAKE GUGWFN COMMON CONSISTENT
C  9 APR 99 - MAF - MULKEN: ONLY MASTER NODE SHOULD CALL TO SHLPOP
C 13 FEB 99 - MAF - FIX VARIOUS POPULATION ANALYSIS BUGS
C 31 DEC 98 - MAF - SHLPOP: FROM HONDO, MULKEN;MOPOP: MOD FOR SHELL POPS
C 12 NOV 98 - MWS - KEEP COMMON MCINP CONSISTENT
C 12 SEP 98 - MWS - BNDORD: ADD SOME REFERENCES TO THE COMMENTS
C 13 APR 98 - MWS - PROPTY: SKIP PROPERTIES OF STATE AVERAGED RUNS
C 27 FEB 98 - MWS - BNDORD,PROPTY,LOWDIN,MULKEN,MOPOP,SPIND: CI NE GUGA
C 18 MAR 97 - MWS - PROPTY: SKIP PI ANALYSIS FOR PCM RUNS
C 18 DEC 96 - JHJ - STNRD: ADD 'BONDS' OPTION, CHANGES REG. MAKEFP
C 13 NOV 96 - MWS - PROPTY: SKIP PI ANALYSIS FOR PCM RUNS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 22 APR 96 - MWS - PROPTY: SKIP PI ANALYSIS FOR PROP RUNS
C  6 APR 95 - MWS - REMOVE ALL AMPERSAND CONTINUATIONS
C 29 MAR 95 - MWS - PROPTY: CALL SHIRO'S NEW PI ORBITAL ANALYSIS
C 28 NOV 94 - MWS - BNDORD: INCLUDE ISTVAN MAYER'S VALENCE COMPUTATION
C 11 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  9 MAR 94 - MWS - BNDORD: FIX FOR OPEN SHELL CASES
C 25 MAR 94 - MR  - BNDORD: INCLUDE GIAMBAGI/MAYER BOND ORDER ANALYSIS
C  9 DEC 93 - MWS - LOWDIN,MULKEN: CHANGE LOWDIN POP. DAF NUMBERS
C 11 MAR 93 - FJ  - PROPTY: ALWAYS DO PROPS IF EXTENDED VIRIAL ANALYSIS
C 20 APR 92 - DIS - STNXYZ: ADJUST LENGTHS FOR HERMIT,WERMIT
C  8 APR 92 - MWS - PROPTY: INSERT COMMENTED OUT CALL TO NBO
C 31 MAR 92 - MWS - STNRD: PARALLEL I/O CHANGES FOR $STONE INPUT
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  9 MAR 92 - JHJ - PROPTY:IF MPC WFN CALL MPCPRP.
C 10 JAN 92 - TLW - CHANGE REWIND TO CALL SEQREW
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 20 OCT 91 - TLW - SPIND: ADD F AND G CAPABILITIES
C  6 AUG 91 - JHJ - PROPTY:ADDED DRG'S CALL TO STONE FOR DMA.
C                   THE LAST 4 SUBROUTINES MUST APPEAR LAST OR DEC3100
C                   COMPILER YIELDS.
C 14 JAN 91 - DRG - ROUTINES STONE, STNRD, STNMOM AND STNXYZ ADDED TO DO
C                   DISTRIBUTED MULTIPOLAR ANALYSIS.
C 31 OCT 90 - MWS - FIX SPIND FOR DUMMY ATOMS, FIX MOPOP FOR ALL SCFTYP,
C                   BETTER POP PRINTOUT, USE /IJPAIR/ IN LOWDIN, MULKEN
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C 25 JUL 90 - MWS - PROPTY: MOVE AIMPAC INTERFACE CALL TO MAIN PROGRAM
C  2 JUN 90 - MWS - CHANGE FORMAT STATEMENT ABOVE MO POPULATION PRINTOUT
C 15 MAR 90 - TJP - ADD MO POPULATIONS
C 23 FEB 90 - MWS - BRING PROPTY BACK HERE (DECSTN QUIRK)
C 24 SEP 89 - MWS - MOVE CHKLOC,ENGANL,GRID,PROPTY,INPELD,INPELF
C                   INPELM,INPPGS,PROPIN,WFNDEN TO PRPLIB, USE MXATM
C  9 OCT 89 - MWS - SPIND NOW WORKS FOR HIGH SPIN GVB CASES.
C 16 DEC 88 - MWS - COMMENT OUT DAWRIT OF LOWDIN ORBS AND BOND ORDERS
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 30 MAY 88 - MWS - USE PARAMETERS TO DIMENSION COMMON
C 22 APR 88 - MWS - USE /GVBWFN/ IN MOPOP
C 14 APR 88 - MWS - FIX MULKEN FOR MORE ATOMS THAN AO-S (PETROS)
C 30 MAR 88 - MWS - FIX MULKEN TO DO ONLY TOTAL POPS FOR GVB,
C                   SUPPRESS ALL BUT TOTAL POP. PRINTOUT IN MULKEN
C 15 NOV 87 - STE - INSERT ABORTS FOR INPUT ERRORS AND USE EXETYP
C 30 APR 87 - STE - SYMTRN: REMOVE USE OF IA
C 12 APR 87 - LAM - PAD COMMON XYZPRP TO CORRECT LENGTH
C  4 NOV 86 - STE - USE PARAMETERS, TAKE NAV OUT OF USE
C  7 JUL 86 - JAB - CHANGE DEXP TO EXP, DFLOAT TO FLOAT
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS,
C                   USE GENERIC FUNCTIONS
C 26 FEB 86 - LAM - FIX MULKEN TO PRINT POPULATION WITH NORMAL PRINT
C  8 JAN 86 - STE - MODIFY COMMON RUNOPT, DELETE FROM MULKEN,SPIND
C 13 NOV 85 - LAM - PRINT LOTS IF NPRINT.GE.8 IN MOPOP AND MULKEN
C                   PRINT LOWDIN ANALYSIS ONLY IF NPRINT.GT.8
C  4 NOV 85 - STE - GRID: INSERT COMMON IOFILE TO DEFINE IW.
C                   INPPGS:MXPNT SHOULD BE MAXPNT.
C                   PROPTY,INPELD,INPELF,INPELM,INPELP,INPPGS,LOWDIN,
C                   MULKEN,SPIND,SYMTRN: CLEAN VARIABLE INITIALIZATION
C  9 SEP 85 - STE - REMOVE 600 LABEL IN PROPTY; 9260 FORMAT IN MULKEN
C 27 JUL 85 - MWS - RENAME COMMON XYZINT TO PRPINT TO AVOID CONFLICTS
C                   ON IBM MACHINES WITH ROUTINE XYZINT
C 14 JUL 85 - MWS - FIX FORMATS IN PROPIN,MULKEN
C  1 JUL 85 - MWS - IMPLEMENT ABBOTT LABS PROPERTIES PACKAGE
C  5 MAY 84 - STE - CHANGE ENERGY DATA IN DICT. RECORD 2 (ENGANL)
C  1 DEC 83 - STE - FIX INDEXING FOR NWDVAR=1 IN LOWDIN,MULKEN
C 17 NOV 83 - STE - SYNTAX CORRECTIONS IN LOWDIN,MULKEN
C  4 NOV 83 - STE - CHANGE GETFM(LOADFM) TO VALFM(LOADFM) IN LOWDIN
C                   CHANGE AOLIM,ATPOP,LATPOP TO ALLOW GHOST ATOMS
C                   DELETE ROUTINES TRANSD,TRNLIN,TRNMAT
C 15 JUN 83 - MWS - FIX BUG FOR SKIPPING PROPERTIES
C  6 MAY 83 - MWS - SKIP CI OR MC RUNTYP=PROP
C 18 MAR 83 - MWS - ADD ENGANL ROUTINE
C  3 OCT 82 - MWS,NDSU - CONVERT TO IBM
C 21 JUN 82 - MWS - NPRINT=8 FOR FULL OUTPUT IN MULKEN,MOPOP
C*MODULE PRPPOP  *DECK AOLIM
      SUBROUTINE AOLIM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /ATLIM / LIMLOW(MXATM),LIMSUP(MXATM)
      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
      LIMLOW(1) = 1
      LAT = 1
      J = 1
      DO 100 I = 1,NSHELL
         IAT = KATOM(I)
         IF (LAT .EQ. IAT) GO TO 100
         LAT = IAT
         LIMSUP(J) = KLOC(I)-1
         J = J+1
         LIMLOW(J) = KLOC(I)
  100 CONTINUE
      LIMSUP(J) = NUM
      IF(J.GE.NAT) GO TO 120
         JP1=J+1
         DO 110 J=JP1,NAT
            LIMLOW(J)=NUM
            LIMSUP(J)=1
  110    CONTINUE
  120 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK ATPOP
      SUBROUTINE ATPOP(A,IA,B,NAT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ATLIM/ LIMLOW(MXATM),LIMSUP(MXATM)
C
      DIMENSION A(*),B(*),IA(*)
C
      PARAMETER (ZERO=0.0D+00)
C
      DO 130 I = 1,NAT
         I1 = LIMLOW(I)
         I2 = LIMSUP(I)
         DO 120 J = 1,NAT
            J1 = LIMLOW(J)
            J2 = LIMSUP(J)
            DUM = ZERO
            IF(I1.GT.I2) GO TO 115
               DO 110 K = I1,I2
                  IF(J1.GT.J2) GO TO 105
                     DO 100 L = J1,J2
                        KK = MAX0(K,L)
                        LL = MIN0(K,L)
                        KL = IA(KK)+LL
                        DUM = DUM+A(KL)
  100                CONTINUE
  105             CONTINUE
  110          CONTINUE
  115       CONTINUE
            II = MAX0(I,J)
            JJ = MIN0(I,J)
            IJ = IA(II)+JJ
            B(IJ) = DUM
  120    CONTINUE
  130 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK BNDORM
      SUBROUTINE BNDORM(DETDEN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DETDEN
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C        ALLOCATE MEMORY FOR BOND ORDER ANALYSIS
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      CALL VALFM(LOADFM)
      LS     = LOADFM + 1
      LDA    = LS     + L2
      LDB    = LDA    + L2
      LQT    = LDB    + L2
      LQS    = LQT    + L3
      LBORD  = LQS    + L3
      LVT    = LBORD  + NAT*NAT
      LVB    = LVT    + NAT
      LAST   = LVB    + NAT
      NEED = LAST - LOADFM -1
C
      CALL GETFM(NEED)
      IF(EXETYP.NE.CHECK)
     *    CALL BNDORD(X(LS),X(LDA),X(LDB),X(LQT),X(LQS),
     *                X(LBORD),X(LVT),X(LVB),L1,L2,NAT,DETDEN)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PRPPOP  *DECK BNDORD
      SUBROUTINE BNDORD(S,DA,DB,QT,QS,BORD,VT,VB,L1,L2,NATM,DETDEN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,OPEN,DETDEN
C
      DIMENSION S(L2),DA(L2),DB(L2),QT(L1,L1),QS(L1,L1),
     *          BORD(NATM,NATM),VT(NATM),VB(NATM)
      DIMENSION IVALS(3),JVALS(3),RVALS(3),BVALS(3)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /ATLIM / LIMLOW(MXATM),LIMSUP(MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, BTHRSH=0.05D+00,
     *           UNIT=0.52917724924D+00)
C
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
      DATA GVB,RMC/8HGVB     ,8HMCSCF   /
      DATA RNONE/8HNONE    /
C
C     --- BOND ORDER ANALYSIS ---
C     WRITTEN BY MICHAEL RAMEK, GRAZ UNIVERSITY OF TECHNOLOGY, 1993.
C     EXTENDED TO COMPUTE THE VALENCES BY ISTVAN MAYER, BUDAPEST, 1994.
C
C     CLOSED SHELL FORMULA PRESENTED BY
C        M.GIAMBAGI, M.GIAMBAGI, D.R.GREMPEL, C.D.HEYMANN
C           J.CHIM.PHYS. 72, 15-22(1975)
C     AND REDISCOVERED BY
C        I.MAYER, CHEM.PHYS.LETT. 97, 270-274 (1983)
C     WITH AN EXTENSION TO OPEN SHELL SYSTEMS (UHF AND ROHF) IN ADDENDUM
C        I.MAYER, 117, 396(1985)
C     SEE AS WELL THE PAPERS
C        I.MAYER, INT.J.QUANTUM CHEM. 26, 151-154(1984)
C        M.S.GIAMBIAGI, M.GIAMBAGI, F.E.JORGE
C           Z.NATURFORSCH. 39A, 1259-73(1984)
C        I.MAYER THEORET.CHIM.ACTA 67, 315-322(1985)
C        I.MAYER INT.J.QUANTUM CHEM. 29, 73-84(1986)
C        I.MAYER INT.J.QUANTUM CHEM. 29, 477-483(1986)
C     THIS LATTER ONE BEING A BIT OF A SUMMARY.
C
C        THE METHOD REQUIRES SEPARATE ALPHA AND BETA DENSITIES,
C        WHICH ARE NOT AVAILABLE FOR GVB, OR GUGA-BASED CSFS.
C
      IF(SCFTYP.EQ.GVB) RETURN
      IF(SCFTYP.EQ.RMC  .AND. .NOT.DETDEN) RETURN
      IF(CITYP.NE.RNONE .AND. .NOT.DETDEN) RETURN
C
C        GET SEPARATE DENSITY MATRICES -DA- AND -DB-
C            NOTE, FOR RHF -DA- IS THE TOTAL DENSITY MATRIX.
C        GET OVERLAP MATRIX IN -S-
C
      OPEN = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF  .OR.  DETDEN
      CALL DAREAD(IDAF,IODA,DA,L2,16,0)
      IF(OPEN) CALL DAREAD(IDAF,IODA,DB,L2,20,0)
      CALL DAREAD(IDAF,IODA,S ,L2,12,0)
C
      IF(.NOT.MASWRK) RETURN
      WRITE(IW,9000) BTHRSH
C
      J = 0
      DO 100 I=1,L1
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
      DO 240 I=1,L1
         DO 230 J=1,L1
            SUMT = ZERO
            SUMS = ZERO
            DO 220 K=1,L1
               IK=IA(MAX(I,K)) + MIN(I,K)
               KJ=IA(MAX(K,J)) + MIN(K,J)
               IF(OPEN) THEN
                  SUMT = SUMT + (DA(IK)+DB(IK))*S(KJ)
                  SUMS = SUMS + (DA(IK)-DB(IK))*S(KJ)
               ELSE
                  SUMT = SUMT + DA(IK)*S(KJ)
               END IF
  220       CONTINUE
            IF(OPEN) THEN
               QT(I,J) = SUMT
               QS(I,J) = SUMS
            ELSE
               QT(I,J) = SUMT
            END IF
  230    CONTINUE
  240 CONTINUE
C
      CALL VCLR(VB,1,NAT)
C
      DO 390 I=1,NAT
         IP1 = I + 1
         BORD(I,I) = ZERO
         IB=LIMLOW(I)
         IE=LIMSUP(I)
         DO 330 J=IP1,NAT
            JB=LIMLOW(J)
            JE=LIMSUP(J)
            B=ZERO
            DO 320 K=IB,IE
               DO 310 L=JB,JE
                  IF(OPEN) THEN
                     B = B + QT(K,L)*QT(L,K) + QS(K,L)*QS(L,K)
                  ELSE
                     B = B + QT(K,L)*QT(L,K)
                  END IF
  310          CONTINUE
  320       CONTINUE
            BORD(I,J) = B
            BORD(J,I) = B
            VB(I) = VB(I) + B
            VB(J) = VB(J) + B
  330    CONTINUE
C
         QE = ZERO
         DO 350 MU=IB,IE
            QE = QE + QT(MU,MU)
  350    CONTINUE
         QE = QE + QE
         X = ZERO
         DO 370 MU=IB,IE
            DO 360 NU=IB,IE
               X = X+QT(MU,NU)*QT(NU,MU)
  360       CONTINUE
  370    CONTINUE
         VT(I) = QE - X
  390 CONTINUE
C
C        PRINT OUT RESULTS
C
      M=0
      DO 450 I=1,NAT
         DO 440 J=I,NAT
            B = BORD(I,J)
            IF(ABS(B).GT.BTHRSH) THEN
               M=M+1
               IVALS(M) = I
               JVALS(M) = J
               R = (C(1,I)-C(1,J))**2
     *           + (C(2,I)-C(2,J))**2
     *           + (C(3,I)-C(3,J))**2
               R = UNIT * SQRT(R)
               RVALS(M) = R
               BVALS(M) = B
               IF(M.EQ.3) THEN
                  WRITE(IW,9020) (IVALS(MM),JVALS(MM),
     *                            RVALS(MM),BVALS(MM),MM=1,3)
                  M=0
               END IF
            END IF
  440    CONTINUE
  450 CONTINUE
      IF(M.GT.0)
     *   WRITE(IW,9020) (IVALS(MM),JVALS(MM),RVALS(MM),BVALS(MM),MM=1,M)
C
      WRITE(IW,9030)
      DO 480 IAT=1,NAT
         VFREE = VT(IAT) - VB(IAT)
         WRITE(IW,9040) IAT,ANAM(IAT),BNAM(IAT),VT(IAT),VB(IAT),VFREE
  480 CONTINUE
      RETURN
C
 9000 FORMAT(/10X,31(1H-)/10X,'BOND ORDER AND VALENCE ANALYSIS',
     *                     5X,'BOND ORDER THRESHOLD=',F5.3/
     *        10X,31(1H-)//
     *   1X,3(18X,'BOND',5X)/
     *   1X,3(1X,'ATOM PAIR',1X,'DIST',2X,'ORDER',5X))
 9020 FORMAT(1X,2I4,2F7.3,5X,2I4,2F7.3,5X,2I4,2F7.3)
 9030 FORMAT(/23X,'TOTAL',7X,'BONDED',8X,'FREE'/
     *       6X,'ATOM',12X,'VALENCE',5X,'VALENCE',5X,'VALENCE')
 9040 FORMAT(1X,I4,1X,A8,A2,1X,3F12.3)
      END
C*MODULE PRPPOP  *DECK DETSPM
      SUBROUTINE DETSPM(DETDEN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DETDEN,DSKSAV
C
      PARAMETER (MXRT=100, MXATM=2000)
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORB,
     *                NA,NB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA AMCSCF/8HMCSCF   /, CHECK/8HCHECK   /
C
C      ALLOCATE MEMORY FOR SPIN DENSITY OF DETERMINANT WAVEFUNCTIONS
C
C         THIS CODE IS AGONIZINGLY SLOW FOR LARGE ACTIVE SPACES
C
      IF(NACT.GE.10) THEN
          DETDEN = .FALSE.
          RETURN
      END IF
C
      SPIN  = (MUL-1)/2.0D+00
      NOCC  = NACT+NCORSV
      NOCC2 = (NOCC*NOCC+NOCC)/2
      L1    = NUM
      L2    = (L1*L1+L1)/2
      L3    = L1*L1
C
      CALL VALFM(LOADFM)
      LCIVEC = LOADFM + 1
      LIFA   = LCIVEC + NCI
      LIAC   = LIFA   + (NACT+1)*(NACT+1)
      LIBC   = LIAC   + NA
      LIACNF = LIBC   + NA
      LIBCNF = LIACNF + NACT
      LACNF2 = LIBCNF + NACT
      LBCNF2 = LACNF2 + NACT
      LDENA  = LBCNF2 + NACT
      LDENB  = LDENA  + NOCC2
      LVEC   = LDENB  + NOCC2
      LWA    = LVEC   + L3
      LWB    = LWA    + L2
      LTA    = LWB    + L2
      LTB    = LTA    + L2
      LAST   = LTA    + L2
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL VCLR(X(LWA),1,L2)
      CALL VCLR(X(LWB),1,L2)
      CALL DAWRIT(IDAF,IODA,X(LWA),L2,16,0)
      CALL DAWRIT(IDAF,IODA,X(LWB),L2,20,0)
C
      IF(EXETYP.EQ.CHECK) GO TO 600
C
      CALL BINOM6(X(LIFA),NACT)
C
      DSKSAV = DSKWRK
      DSKWRK=.FALSE.
C
      CALL SEQREW(NFT12)
      IF (MASWRK) READ(NFT12) NSTATS,NDETS
      IF (GOPARR) CALL DDI_BCAST(2510,'I',NSTATS,1,MASTER)
      IF (GOPARR) CALL DDI_BCAST(2511,'I',NDETS ,1,MASTER)
C
      IF(NDETS.NE.NCI) THEN
         IF(MASWRK) WRITE(IW,*) 'NDETS.NE.NCI IN -DETSPM-'
         CALL ABRT
      END IF
C
      JJ = 0
      DO 77 J=1,NSTATS
         CALL SQREAD(NFT12,X(LCIVEC),NCI)
         IF (SCFTYP.EQ.AMCSCF) THEN
           IF(IPURES.EQ.1 .AND. ABS(SPINS(J)-SPIN).GT.0.1D+00) GO TO 77
           JJ=JJ+1
           IF (WSTATE(JJ).NE.ZERO) THEN
             CALL DETSPD(X(LCIVEC),X(LIAC),X(LIBC),X(LIFA),X(LIACNF),
     *                   X(LIBCNF),X(LACNF2),X(LBCNF2),X(LDENA),
     *                   X(LDENB),X(LVEC),X(LWA),X(LWB),X(LTA),X(LTB),
     *                   JJ,1,NA,NACT,NOCC,NOCC2,L1,L2,L3)
           END IF
         ELSE
           IF (J.EQ.IROOT) THEN
             CALL DETSPD(X(LCIVEC),X(LIAC),X(LIBC),X(LIFA),X(LIACNF),
     *                   X(LIBCNF),X(LACNF2),X(LBCNF2),X(LDENA),
     *                   X(LDENB),X(LVEC),X(LWA),X(LWB),X(LTA),X(LTB),
     *                   J,2,NA,NACT,NOCC,NOCC2,L1,L2,L3)
           END IF
         END IF
 77   CONTINUE
C
      CALL SEQREW(NFT12)
      DSKWRK = DSKSAV
C
  600 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PRPPOP  *DECK DETSPD
      SUBROUTINE DETSPD(CI,IACN,IBCN,IFA,IACNF,IBCNF,IACNF2,IBCNF2,
     *                  DNA,DNB,V,WRKA,WRKB,TMPA,TMPB,
     *                  NJ,NTYP,NA,NACT,NOCC,NOCC2,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXRT=100, MXAO=8192)
C
      DIMENSION CI(NCI),IACN(NA),IBCN(NA),IFA(0:NACT,0:NACT),
     *          IACNF(NACT),IBCNF(NACT),IACNF2(NACT),IBCNF2(NACT),
     *          DNA(NOCC2),DNB(NOCC2),V(L1,L1),
     *          WRKA(L2),WRKB(L2),TMPA(L2),TMPB(L2)
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NNACT,NORB,
     *                NNA,NNB,K,KST,IROOT,IPURES,MAXW1,NITER,MAXP,NCI,
     *                IGPDET,KSTSYM,NFTGCI
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
      COMMON /PSILVL/ IPSI,ISKPRP
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- ALPHA AND BETA DENSITY FOR DETERMINANT-BASED CI -----
C
      CALL VCLR(DNA,1,NOCC2)
      CALL VCLR(DNB,1,NOCC2)
C
      DO 11 II=1,NCORSV
         IL=IA(II)+II
         DNA(IL)=ONE
         DNB(IL)=ONE
  11  CONTINUE
      NALP = IFA(NACT,NNA)
      NBLP = IFA(NACT,NNB)
C
C     ***** DIAGONAL CONTRIBUTION ******
C
      DO 111 JJ=1,NCI
         CALL DETVAL(JJ,NNA,NNB,NACT,IACN,IBCN,IACNF,IBCNF,IFA)
         DO 22 II=1,NACT
            IL=IA(NCORSV+II)+(NCORSV+II)
            IF(IACNF(II).EQ.1) DNA(IL)=DNA(IL)+(CI(JJ)**2)
            IF(IBCNF(II).EQ.1) DNB(IL)=DNB(IL)+(CI(JJ)**2)
  22     CONTINUE
 111  CONTINUE
C
C     ***** NOW FOR THE NON-DIAGONAL PART *****
C        THE ALPHA PART FIRST
C
      DO 555 JJ=1,(NCI-NBLP+1),NBLP
         CALL DETVAL(JJ,NNA,NNB,NACT,IACN,IBCN,IACNF,IBCNF,IFA)
         DO 444 KK=JJ+NBLP,(NCI-NBLP+1),NBLP
            CALL DETVAL(KK,NNA,NNB,NACT,IACN,IBCN,IACNF2,IBCNF2,IFA)
            CALL DETCMP(NACT,M,N,IAB,NPR,IACNF,IBCNF,IACNF2,IBCNF2)
            IF (IAB.EQ.0) GO TO 444
            IL=IA(NCORSV+MAX(M,N))+(NCORSV+MIN(M,N))
            MNS=(-1)**NPR
            DNA(IL)=DNA(IL)+MNS*(CI(JJ)*CI(KK))
            DO 333 III=1,NBLP-1
               DNA(IL)=DNA(IL)+MNS*(CI(JJ+III)*CI(KK+III))
 333        CONTINUE
 444     CONTINUE
 555  CONTINUE
C
C        NOW GOES THE BETA PART
C
      DO 999 JJ=1,NBLP
         CALL DETVAL(JJ,NNA,NNB,NACT,IACN,IBCN,IACNF,IBCNF,IFA)
         DO 777 KK=JJ+1,NBLP
            CALL DETVAL(KK,NNA,NNB,NACT,IACN,IBCN,IACNF2,IBCNF2,IFA)
            CALL DETCMP(NACT,M,N,IAB,NPR,IACNF,IBCNF,IACNF2,IBCNF2)
            IF (IAB.EQ.0) GO TO 777
            IL=IA(NCORSV+MAX(M,N))+(NCORSV+MIN(M,N))
            MNS=(-1)**NPR
            DNB(IL)=DNB(IL)+MNS*(CI(JJ)*CI(KK))
            DO 666 III=1,NALP-1
               DNB(IL)=DNB(IL)+MNS*(CI((NBLP*III)+JJ)*CI((NBLP*III)+KK))
 666        CONTINUE
 777     CONTINUE
 999  CONTINUE
C
C     ----- TRANSFORM SPIN DENSITY FROM MO TO AO BASIS -----
C     WE MUST BE CAREFUL TO GET EXACTLY THE ORBITALS FOR WHICH
C     THE CI VECTOR WAS COMPUTED, I.E. AVOIDING ANY ACTIVE ORBITAL
C     PHASE CHANGES THAT MIGHT OCCUR IN THE CANONICALIZATION ROUTINE.
C
      IF(IPSI.EQ.0) THEN
         CALL DAREAD(IDAF,IODA,V,L3,325,0)
      ELSE
         CALL DAREAD(IDAF,IODA,V,L3, 15,0)
      END IF
      CALL TRPOSQ(V,L1)
      CALL TFTRI(WRKA,DNA,V,TMPA,L1,NOCC,L1)
      CALL TFTRI(WRKB,DNB,V,TMPB,L1,NOCC,L1)
C
C     FOR STATE AVERAGED MCSCF ONLY
C
      IF (NTYP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,TMPA,L2,16,0)
         CALL DAREAD(IDAF,IODA,TMPB,L2,20,0)
         DUM=WSTATE(NJ)
         DO 1111 J=1,L2
            WRKA(J)=TMPA(J)+DUM*WRKA(J)
            WRKB(J)=TMPB(J)+DUM*WRKB(J)
 1111    CONTINUE
      END IF
C
C     WRITE NEW SPIN DENSITY TO IDAF
C
      CALL DAWRIT(IDAF,IODA,WRKA,L2,16,0)
      CALL DAWRIT(IDAF,IODA,WRKB,L2,20,0)
      RETURN
      END
C*MODULE PRPPOP  *DECK DETCMP
      SUBROUTINE DETCMP(NACT,M,N,IAB,NPR,IACNF,IBCNF,IACNF2,IBCNF2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IACNF(NACT),IBCNF(NACT),IACNF2(NACT),IBCNF2(NACT)
C
      IAB =0
      NDIFA=0
      NDIFB=0
      ICA=0
      ICA2=0
      DO 111 J=1,NACT
         IF (IACNF(J).EQ.1) ICA=ICA+1
         IF (IACNF2(J).EQ.1) ICA2=ICA2+1
         IF (IACNF(J).EQ.IACNF2(J)) GO TO 111
         IF (NDIFA.GE.2) GO TO 666
         IF (NDIFA.EQ.0) THEN
            M = J
            IF (IACNF(J).EQ.1) MC=ICA
            IF (IACNF2(J).EQ.1) MC=ICA2
         END IF
         IF (NDIFA.EQ.1) THEN
            N = J
            IF (IACNF(J).EQ.1) NC=ICA
            IF (IACNF2(J).EQ.1) NC=ICA2
         END IF
         NDIFA = NDIFA + 1
 111  CONTINUE
      NPR=MC+NC
      ICB=0
      ICB2=0
      DO 222 J=1,NACT
         IF (IBCNF(J).EQ.1) ICB=ICB+1
         IF (IBCNF2(J).EQ.1) ICB2=ICB2+1
         IF (IBCNF(J).EQ.IBCNF2(J)) GO TO 222
         IF (NDIFA.NE.0) GO TO 666
         IF (NDIFB.GE.2) GO TO 666
         IF (NDIFB.EQ.0) THEN
            M = J
            IF (IBCNF(J).EQ.1) MC=ICB
            IF (IBCNF2(J).EQ.1) MC=ICB2
         END IF
         IF (NDIFB.EQ.1) THEN
            N = J
            IF (IBCNF(J).EQ.1) NC=ICB
            IF (IBCNF2(J).EQ.1) NC=ICB2
         END IF
         NDIFB = NDIFB + 1
 222  CONTINUE
      IF (NDIFB.EQ.2) NPR=MC+NC
      IF (NDIFA.EQ.2) IAB = 1
      IF (NDIFB.EQ.2) IAB = 2
 666  RETURN
      END
C*MODULE PRPPOP  *DECK DETRST
      SUBROUTINE DETRST
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,MA,MB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKO,IJKT,IDAF,NAV,IODA(950)
C
C         RESTORE THE TOTAL MCSCF DENSITY
C
      L2 = (NUM*NUM+NUM)/2
      CALL VALFM(LOADFM)
      LDA  = LOADFM + 1
      LDB  = LDA    + L2
      LAST = LDB    + L2
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL DAREAD(IDAF,IODA,X(LDA),L2,16,0)
      CALL DAREAD(IDAF,IODA,X(LDB),L2,20,0)
      CALL VADD(X(LDA),1,X(LDB),1,X(LDA),1,L2)
      CALL DAWRIT(IDAF,IODA,X(LDA),L2,16,0)
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PRPPOP  *DECK GROSSC
      SUBROUTINE GROSSC(A,B,IA,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(*),IA(N)
C
      PARAMETER (ZERO=0.0D+00)
C
      DO 120 I = 1,N
         DUM = ZERO
         DO 100 J = 1,N
            II = MAX0(I,J)
            JJ = MIN0(I,J)
            IJ = IA(II)+JJ
            DUM = DUM+A(IJ)
  100    CONTINUE
         B(I) = DUM
  120 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK PROPTY
      SUBROUTINE PROPTY(PRPTYP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION MAKEFP,METHMC
C
      CHARACTER*4 PRPTYP
C
      LOGICAL RHO,GOPARR,DSKWRK,MASWRK,VTSCAL,VIROK,LVCLN,
     *        CANONC,FCORE,FORS,NOCI,EKT,LINSER,ESAME,DETDEN,
     *        ISEPS,USEPS,MNMEDG,MNMEOP,UNVGSS,DGAPRX,RDCISV,
     *        TRIPLET,SG1T
C
      PARAMETER (MXATM=2000, MXAO=8192, MXRT=100, MXNORO=250,
     *           MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /CISPAR/ HAMTYP,DIAGZN,DAVCVG,PRTCIS,NSTCIS,ISTATE,MULCIS,
     &                MXV,NDAVIT,ICISPR,NACORE,NBCORE,NOA,NOB,NORBOC,
     &                NBFCIS,NGSVEC,MNMEDG,MNMEOP,ICLOBBR,UNVGSS,DGAPRX,
     &                RDCISV
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTDET,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /EOMPAR/ CVGCI,CVGEOM,GRPEOM,NSTEOM(8),NOACT,NUACT,
     *                MOACTCC(MXAO),MTHTRIP,MTHCI,MTHEOM,MTHINIT,
     *                MAXCI,MAXEOM,MICCI,MICEOM,IROOTCC(2),
     *                IPROPCC,IPROPCCE
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /ISEPS / ISEPS,USEPS
      COMMON /MCINP / METHMC,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /NEOJOB/ NEORUN,NELERM
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,NFMOPCM
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SVPOPT/ ISVP,NSVP
      COMMON /VIRIAL/ SCALTE,SCALTT,GVIR,VTCONV,MAXVT,VTSCAL,VIROK,LVCLN
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (ZERO=0.0D+00)
C
      DATA   NONE/4HNONE/
      DATA  RNONE/8HNONE    /
      DATA    GVB/8HGVB     /
      DATA AMCSCF/8HMCSCF   /
      DATA ALDET,GUGA,GENCI,ORMAS,GMCCI
     *      /8HALDET   ,8HGUGA    ,8HGENCI   ,8HORMAS   ,8HGMCCI   /
      DATA EOMCCSD/8HEOM-CCSD/
      DATA PROP  /8HPROP    /
      DATA  CHECK/8HCHECK   /
      DATA MAKEFP/8HMAKEFP  /
      DATA     C1/8HC1      /
C
C             ----- DRIVER FOR PROPERTY EVALUATION -----
C     IF STATE-SPECIFIC DENSITY MATRIX EXISTS, COMPUTE ALL PROPERTIES.
C
C     DO WE HAVE A DENSITY MATRIX FOR PROPERTIES?
C     1. FOR RUNTYP=PROP, WHERE ALL WE KNOW IS THE $VEC GIVEN AS INPUT,
C        SO ONLY VERY SIMPLE HF RUNS WILL HAVE THE DENSITY AVAILABLE.
C        WARNING: WE ASSUME GVB USER GAVE ANY NECESSARY -CICOEF- TOO.
C     2. THE CALLER HAS THE RESPONSIBILITY OF DECIDING WHAT KIND OF
C        DENSITY MATRIX EXISTS, SO WE CAN PRINT AN APPRPRIATE HEADER.
C
      RHO = .TRUE.
      IF(RUNTYP.EQ.PROP) THEN
         IF(SCFTYP.EQ.AMCSCF) RHO=.FALSE.
         IF(IPSI.GT.0)        RHO=.FALSE.
      END IF
C
C        SHOULD WE JUST EXIT?
C
      IF(.NOT.RHO) RETURN
      IF(NPRINT.EQ.-5  .AND.  .NOT.VTSCAL) RETURN
      IF(ISKPRP.NE.0) RETURN
C
C        IF A NUCLEAR/ELECTRON ORBITAL RUN, CALL NEO PROPERTIES
C
      IF(NEORUN.EQ.1) THEN
         CALL PRPNEO(PRPTYP)
         RETURN
      END IF
C
C        NON-VARIATIONAL METHODS MAY HAVE EXPECTATION DENSITY
C        OR RELAXED DENSITY MATRICES.
C
      IF(MASWRK) THEN
         IF(PRPTYP.EQ.'PROP') WRITE(IW,9210) SCFTYP
         IF(PRPTYP.EQ.'SCF ') WRITE(IW,9210) SCFTYP
         IF(PRPTYP.EQ.'DFT ') WRITE(IW,9220) DFTYPE,SCFTYP
         IF(PRPTYP.EQ.'VB  ') WRITE(IW,9230)
         IF(PRPTYP.EQ.'MP2 ') WRITE(IW,9240)
         IF(PRPTYP.EQ.'CI E') WRITE(IW,9250) CITYP,IROOT
         IF(PRPTYP.EQ.'CI R') WRITE(IW,9260) CITYP,IROOT
         IF(PRPTYP.EQ.'CISE') WRITE(IW,9270) ISTATE
         IF(PRPTYP.EQ.'CISR') WRITE(IW,9280) ISTATE
         IF(PRPTYP.EQ.'CCSD') WRITE(IW,9290)
         IF(PRPTYP.EQ.'EOMD') WRITE(IW,9300) IROOTCC(2),IROOTCC(1)
         IF(PRPTYP.EQ.'TDDF') WRITE(IW,9310) DFTYPE,NTHST
      END IF
C
C     ----- MOPAC USES STO BASIS FUNCTIONS, AND ITS OWN ROUTINE -----
C
      IF(MPCTYP.NE.NONE) THEN
         CALL MPCPRP(RHO,.TRUE.)
         RETURN
      END IF
C
C        SEE IF WE ARE RUNNING A SOLVATION MODEL
C
      NSOLV = 0
      IF(IZRF.EQ.1) NSOLV=NSOLV+1
      IF(IPCM.EQ.1) NSOLV=NSOLV+1
      IF(NFRG.GT.0) NSOLV=NSOLV+1
      IF(ISEPS)     NSOLV=NSOLV+1
      IF(ISVP.EQ.1) NSOLV=NSOLV+1
C
C        IF STATE-AVERAGED MCSCF, WE DON'T HAVE A STATE-SPECIFIC
C        DENSITY MATRIX AVAILABLE TO DO PROPERTIES, UNLESS WE
C        HAVE BEEN AVERAGING OVER DEGENERATE LEVELS
C
      IF(SCFTYP.EQ.AMCSCF) THEN
         IF(CISTEP.EQ.ALDET  .OR.
     *      CISTEP.EQ.GENCI  .OR.
     *      CISTEP.EQ.ORMAS  .OR.
     *      CISTEP.EQ.GMCCI) THEN
            NOCC = NCORSV + NACT
C
            NAVG=0
            DO 10 ISTAT=1,KDET
               IF(WSTATE(ISTAT).GT.0.0D+00) NAVG=NAVG+1
   10       CONTINUE
            IF(NAVG.EQ.1) GO TO 50
C
            NXTW=1
            NXTR=0
            ESAME=.TRUE.
            EMATCH= ZERO
            DO 20 ISTAT=1,KDET
               IF(IPURES.EQ.1  .AND.
     *            ABS(SPINS(ISTAT)-SDET).GT.0.03D+00) GO TO 20
               NXTR=NXTR+1
               IF(NXTR.EQ.IWTS(NXTW)) THEN
                  IF(EMATCH.EQ.ZERO) EMATCH = ESTATE(ISTAT)
                  ETEST = ESTATE(ISTAT)
                  ESAME = ESAME  .AND.  ABS(ETEST-EMATCH) .LT. 1.0D-07
                  NXTW=NXTW+1
               END IF
               IF(NXTW.GT.NAVG) GO TO 50
   20       CONTINUE
         END IF
C            WSTATE IS NOT SAVED TO MAKE ANALAGOUS TEST
         IF(CISTEP.EQ.GUGA) THEN
            NOCC = NMCC + (NDOC + NALP + NAOS + NBOS + NVAL)
            NAVG=1
         END IF
C
   50    CONTINUE
         IF(NAVG.GT.1) THEN
            IF(ESAME) THEN
               IF(MASWRK) WRITE(IW,9030)
            ELSE
               IF(MASWRK) WRITE(IW,9040) NAVG
               RETURN
C              ******
            END IF
         END IF
C            LIE ABOUT NUMBER OF E-, THIS MUST BE UNDONE BEFORE EXIT
         NASAVE = NA
         NA = NOCC
      END IF
C
C     ----- DETERMINE NUMBER OF BASIS FUNCTIONS PER ATOM -----
C
      CALL AOLIM
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
C
C        ----- ENERGY ANALYSIS -----
C     RUNTYP=PROP LACKS THE TOTAL ENERGY, SO WE CAN'T GET THE 2E- ENERGY
C     BY SUBTRACTING ALL 1-E ENERGIES (SO JUST SKIP).
C
C     CCSD AND EOM-CCSD ARE BOTH RUNTYP=EOM-CCSD AT THIS POINT, AND
C     HAVE THE APPROPRIATE TOTAL ENERGY STORED.  HOWEVER, CR-EOM OR
C     CR-CCL WILL HAVE A TRIPLES CORRECTION ADDED INTO -ETOT-, SO THE
C     CORRECT CCSD/EOM-CCSD ENERGY IS NOT AVAILABLE.
C
C     SOLVENT MODEL WOULD REQUIRE THINKING ABOUT THE VALIDITY OF IT ALL.
C
      IF(RUNTYP.EQ.PROP)  GO TO 101
      IF(CCTYP.NE.RNONE  .AND.  CCTYP.NE.EOMCCSD) GO TO 101
      IF(NSOLV.GT.0)      GO TO 101
C
      CALL VALFM(LOADFM)
      LDA    = LOADFM + 1
      LDB    = LDA    + L2
      LONEEI = LDB    + L2
      LLABPI = LONEEI + L2
      LFA    = LLABPI + L1
      LFB    = LFA    + L2
      LT     = LFB    + L2
      LAST   = LT     + L2
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 100
C
      CALL ENGANL(X(LDA),X(LDB),X(LONEEI),L1,L2)
C
      IF(IPSI.GT.0) GO TO 100
      IF(SCFTYP.EQ.GVB  .OR.  SCFTYP.EQ.AMCSCF) GO TO 100
C
      CALL PIANL(X(LONEEI),X(LT),X(LFA),X(LFB),
     *           X(LDA),X(LDB),X(LLABPI),L1,L2)
C
  100 CONTINUE
      CALL RETFM(NEED)
  101 CONTINUE
C
C        VIRIAL THEOREM SCALING MUST GET VIRIAL RATIO, EVEN IF
C        WE AREN'T PRINTING, BUT NEED NOT GO PAST -ENGANL-.
C
      IF (VTSCAL) RETURN
C
C            MWS: CONDITIONALLY CALL THE 2 POPULATION ANALYSIS CODES.
C            REMOVE THIS NEXT BIT, AND ITS MATCHING ENDIF LATER,
C            AFTER -SPHMAT- LEARNS ABOUT H+I.
    8 FORMAT(/1X,'SKIPPING POP.ANALYSIS DUE TO H+I AO''S IN THE BASIS,',
     *    ' AND ISPHER.GE.0')
      CALL BASCHK(LMAX)
      IF(LMAX.GE.5  .AND.  ISPHER.GE.0) THEN
         IF(MASWRK) WRITE(IW,8)
      ELSE
C
C        ----- SYMMETRICALLY TRANSFORM VECTORS AND DENSITY MATRIX -----
C
      CALL LOWDIN
C
C        ----- MULLIKEN AND LOWDIN POPULATION ANALYSIS -----
C
      CALL MULKEN
      END IF  ! END OF SKIPPING H+I POP.ANALYSIS
C
C        ----- IF DETERMINANTS USED, COMPUTE ALPHA/BETA DENSITY -----
C
      DETDEN = .FALSE.
      IF(IPSI.EQ.0) DETDEN= (SCFTYP.EQ.AMCSCF. AND. CISTEP.EQ.ALDET)
      IF(IPSI.EQ.1) DETDEN=   CITYP.EQ.ALDET
C------------PREVENT COMPUTATION OF SEPARATE ALPHA/BETA DENSITY---------
C    EXCEPT IN THE CASE WHERE THERE IS NO POINT GROUP SYMMETRY EXPLOITED
C    THIS IS DUE TO A C1 ASSUMPTION IN DETSPM'S FORMATION OF A+B DENSITY
      DETDEN = DETDEN .AND. GRPDET.EQ.C1
      IF (DETDEN) CALL DETSPM(DETDEN)
C
C        ----- ISTVAN MAYER'S (AND OTHER'S!) BOND ORDER ANALYSIS -----
C
      CALL BNDORM(DETDEN)
C
C        ----- ATOMIC SPIN DENSITY (AT ALL NUCLEI) -----
C
      CALL SPIND(DETDEN)
C
C        ----- RESTORE DETERMINANT MCSCF/CI'S TOTAL DENSITY -----
C
      IF (DETDEN) CALL DETRST
C
C        ----- SELECT DESIRED ELECTROSTATIC PROPERTIES -----
C        MULTIPOLE MOMENTS, MONOPOLE THROUGH OCTUPOLE
C        ELECTRON DENSITY, AT NUCLEI, USUALLY
C        ELECTRIC POTENTIAL, AT NUCLEI, OR TO FIT ATOMIC CHARGES
C        ELECTRIC FIELD/GRADIENT, 1ST/2ND DERIVATIVES OF POTENTIAL
C
      IF((IEMOM+IEDEN+IEFLD+IEPOT).GT.0   .AND.  NSOLV.GT.0) THEN
         IF(MASWRK) WRITE(IW,9050)
      END IF
C
      IF(IEMOM.NE.0) CALL ELMOMC
      IF(IEDEN.NE.0) CALL ELDENC
      IF(IEPOT.NE.0) CALL ELPOTC
      IF(IEFLD.NE.0) CALL ELFLDC
C
C        ----- ANTHONY STONE'S DISTRIBUTED MULTIPOLE ANALYSIS -----
C
      IF(RUNTYP.NE.MAKEFP) CALL STONE(IP,RNONE)
C
C        --- HUI LI'S DENSITY BASED DISTRIBUTED MULTIPOLE ANALYSIS ---
C
      IF(RUNTYP.EQ.MAKEFP) CALL DENMUL
C
C        ----- ROGER AMOS'S DISTRIBUTED MULTIPOLE ANALYSIS -----
C
      IF(ISEPS) CALL NUMPROP
C
      IF (MASWRK) WRITE(IW,9010)
      CALL TIMIT(1)
C
C     CALL THE NBO PROGRAM FROM FRANK WEINHOLD'S GROUP IN MADISON.
C     E-MAIL CONTACT IS WEINHOLD@CHEM.WISC.EDU
C     NOTE THAT NBO IS STRICTLY A PROGRAM FOR SEQUENTIAL EXECUTION.
C
C---  IF(NPROC.EQ.1) CALL RUNNBO
C---  IF(MASWRK) WRITE(IW,FMT='('' ......END OF NBO ANALYSIS......'')')
C---  CALL TIMIT(1)
C
      IF(SCFTYP.EQ.AMCSCF) NA = NASAVE
      RETURN
C
 9010 FORMAT(1X,'...... END OF PROPERTY EVALUATION ......')
 9030 FORMAT(/1X,'ALL STATES AVERAGED IN THE MCSCF ARE DEGENERATE,',
     *          ' THEREFORE'/
     *       1X,'PROPERTIES WILL BE COMPUTED USING THE AVERAGE',
     *          ' DENSITY MATRIX'/)
 9040 FORMAT(/1X,'THIS MCSCF AVERAGED',I4,' ROOTS OF DIFFERENT ENERGY.'/
     *       1X,'SINCE NO STATE-SPECIFIC DENSITY MATRIX IS AVAILABLE,'/
     *       1X,'PROPERTIES WILL NOT BE COMPUTED.'/)
 9050 FORMAT(/1X,'   **** A SOLVENT MODEL IS IN USE IN THIS RUN ****'/
     *   1X,'ELECTRICAL PROPERTIES',
     *      ' (MOMENTS, DENSITY, FIELD/GRADIENT, AND/OR POTENTIAL)'/
     *   1X,'ARE COMPUTED USING THE SOLVATED SYSTEM''S DENSITY MATRIX,',
     *      ' BUT'/
     *   1X,'ONLY WITH AO INTEGRAL CONTRIBUTIONS FROM THE QUANTUM',
     *      ' SOLUTE.')
C
 9210 FORMAT(/5X,64(1H-)/
     *        5X,'PROPERTY VALUES FOR THE ',A5,
     *           ' SELF-CONSISTENT FIELD WAVEFUNCTION'/
     *        5X,64(1H-))
 9220 FORMAT(/5X,69(1H-)/
     *        5X,'PROPERTIES FOR THE ',A8,' DFT FUNCTIONAL ',
     *           '(',A4,' TYPE) DENSITY MATRIX'/
     *        5X,69(1H-))
 9230 FORMAT(/5X,34(1H-)/
     *        5X,'PROPERTIES FOR THE VB WAVEFUNCTION'/
     *        5X,34(1H-))
 9240 FORMAT(/5X,49(1H-)/
     *        5X,'MP2 PROPERTIES...FOR THE FIRST ORDER WAVEFUNCTION'/
     *       10X,'USING THE RESPONSE (OR RELAXED) DENSITY'/
     *        5X,49(1H-))
 9250 FORMAT(/5X,56(1H-)/
     *        5X,A5,' CI PROPERTIES...FOR THE WAVEFUNCTION OF',
     *           ' STATE',I5/
     *       15X,'USING THE EXPECTATION VALUE DENSITY'/
     *        5X,56(1H-))
 9260 FORMAT(/5X,56(1H-)/
     *        5X,A5,' CI PROPERTIES...FOR THE WAVEFUNCTION OF',
     *           ' STATE',I5/
     *       13X,'USING THE RESPONSE (OR RELAXED) DENSITY'/
     *        5X,56(1H-))
 9270 FORMAT(/5X,59(1H-)/
     *        5X,'CIS PROPERTIES...FOR THE WAVEFUNCTION OF',
     *          ' EXCITED STATE',I5/
     *       17X,'USING THE EXPECTATION VALUE DENSITY'/
     *        5X,59(1H-))
 9280 FORMAT(/5X,59(1H-)/
     *        5X,'CIS PROPERTIES...FOR THE WAVEFUNCTION OF',
     *           ' EXCITED STATE',I5/
     *       15X,'USING THE RESPONSE (OR RELAXED) DENSITY'/
     *        5X,59(1H-))
 9290 FORMAT(/5X,38(1H-)/
     *        5X,'CCSD PROPERTIES...FOR THE GROUND STATE'/
     *        6X,'USING THE EXPECTATION VALUE DENSITY'/
     *        5X,38(1H-))
 9300 FORMAT(/5X,64(1H-)/
     *        5X,'EOM-CCSD PROPERTIES...FOR EXCITED STATE',I4,
     *           ' OF SPACE SYMMETRY',I3/
     *       20X,'USING THE EXPECTATION VALUE DENSITY'/
     *        5X,64(1H-))
 9310 FORMAT(/5X,42(1H-)/
     *        5X,'PROPERTIES FOR THE ',A8,' DFT FUNCTIONAL',
     *       /5X,' USING THE DENSITY OF EXCITED STATE',I5/
     *        5X,42(1H-))
      END
C*MODULE PRPPOP  *DECK LATPOP
      SUBROUTINE LATPOP(DD,DA,NAT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DD(*),DA(*)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ATLIM / LIMLOW(MXATM),LIMSUP(MXATM)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- COMPRESS FROM ORBITALS TO ATOMS FOR LOWDIN ANALYSIS -----
C
      DO 120 I = 1,NAT
         DA(I) = ZERO
         I1 = LIMLOW(I)
         I2 = LIMSUP(I)
         IF(I1.GT.I2) GO TO 110
            DO 100 J = I1,I2
               DA(I) = DA(I)+DD(J)
  100       CONTINUE
  110    CONTINUE
  120 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK LOWDIN
      SUBROUTINE LOWDIN
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT,BETA,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA ROHF/8HROHF    /, UHF/8HUHF     /, GVB/8HGVB     /
      DATA CHECK/8HCHECK   /, RNONE/8HNONE    /
C
C     ----- LOWDIN POPULATION ANALYSIS -----
C
      OUT = NPRINT.GT.8 .AND. MASWRK
      BETA = (SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF  .OR.  SCFTYP.EQ.GVB)
     *   .AND.   (CITYP.EQ.RNONE)
C
      L0 = NQMT
      L1 = NUM
      L2 = L1*(L1+1)/2
      L3 = L1*L1
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
      LOADFM=0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10+L2
      I21 = I10+L3
      I30 = I20+L2
      I31 = I20+L3
      I40 = I30+L2
      I41 = I30+L3
      I60 = I40+L2
      I70 = I60+L1
      I80 = I70 + 9*L1
C
C     ----- GET MEMORY
C
      LASTC = I80
      NEED = LASTC-I10
      CALL GETFM(NEED)
      IF ((NPRINT .EQ. 3) .AND. MASWRK) WRITE(IW,9008)
     *     I10,I20,I21,I30,I31,I40,I41,I60,LASTC
      IF (EXETYP .EQ. CHECK) GO TO 180
C
      DO 100 I = 1,L1
         IA(I) = I*(I-1)/2
  100 CONTINUE
C
C     ----- COMPUTE THE SYMMETRIC TRANSFORMATION MATRIX -----
C     FOR INFORMATION ON THE SYMMETRIC ORTHOGONALIZATION SEE
C              "ON THE NONORTHOGONALITY PROBLEM"
C          P.-O.-LOWDIN  ADV.QUANTUM CHEM. 5, 185(1970)
C
C     -S - AT X(I10)     OVERLAP MATRIX
C     -SV- AT X(I20)     EIGENVECTORS OF S
C     -SE- AT X(I31)     EIGENVALUES OF S
C     -SH- AT X(I10)     SYMMETRIC TRANSFORMATION MATRIX
C
      CALL DAREAD(IDAF,IODA,X(I10),L2,12,0)
      CALL SYMTRN(X(I10),X(I20),X(I31),X(I10),X(I70),
     *            L1,L1)
      CALL DAWRIT(IDAF,IODA,X(I10),L2,51,0)
C
C     ----- TRANSFORM THE VECTORS AND THE DENSITY MATRIX
C           TO THE SYMMETRIC ORTHOGONAL BASIS SET
C
C     -SH- AT X(I30)     SYMMETRIC TRANSFORMATION MATRIX
C     -V - AT X(I10)
C     -VO- AT X(I10)     VECTORS IN ORTHOGONAL BASIS SET
C     -VO2-AT X(I30)     ORBITAL POPULATION MATRIX
C                        VO2(I,J)= VO(I,J)**2*SIGN(VO(I,J))
C     -DO- AT X(I30)     DENSITY MATRIX IN THE ORTHOGONAL BASIS
C     -DOD-AT X(I60)     DIAGONAL OF DO
C
      NDAFV = 15
      NDAFD = 16
      NDAFE = 17
      NDAFL = 69
C
      IPASS = 0
  140 CONTINUE
      IF(OUT) THEN
         IF(BETA) THEN
            IF(IPASS.EQ.0) WRITE(IW,9028)
            IF(IPASS.EQ.1) WRITE(IW,9048)
         ELSE
            WRITE(IW,9068)
         END IF
      END IF
C
C     ----- CALCULATE ORBITAL POPULATIONS -----
C
      CALL DAREAD(IDAF,IODA,X(I30),L2,51,0)
      CALL DAREAD(IDAF,IODA,X(I10),L3,NDAFV,0)
      CALL SYMTRV(X(I30),X(I10),X(I10),X(I30),X(I60),L0,L1,L1)
      CALL DAREAD(IDAF,IODA,X(I60),L1,NDAFE,0)
      IF(OUT) CALL PREV(X(I30),X(I60),L1,L1,L1)
      CALL DAREAD(IDAF,IODA,X(I30),L2,51,0)
      CALL DAREAD(IDAF,IODA,X(I40),L2,NDAFD,0)
      CALL SYMTRD(X(I30),X(I40),X(I10),X(I40),X(I60),IA,L1,L1)
C
C     ----- PRINT THE LOWDIN BOND ORDER MATRIX -----
C
      IF(OUT) WRITE(IW,9088)
      IF(OUT) CALL PRTRIL(X(I40),L1)
      CALL DAWRIT(IDAF,IODA,X(I60),L1,NDAFL,0)
      IPASS = IPASS+1
      IF(BETA  .AND.  IPASS.LT.2) THEN
         NDAFV = NDAFV+4
         NDAFD = NDAFD+4
         NDAFE = NDAFE+4
         NDAFL = NDAFL+1
         GO TO 140
      END IF
C
  180 CONTINUE
      CALL RETFM(NEED)
      RETURN
C
 9008 FORMAT(' MEMORY ASSIGNEMENT'/31H I10, I20, I21, I30, I31, I40,
     *      ,'I41, I60 ='/8I8/8H LAST = ,I8)
 9028 FORMAT(/10X,58(1H-)/10X,23HLOWDIN ORBITAL (VO(I,J), 5H**2)*,
     *     30HSIGN(VO(I,J))   ALPHA ORBITALS/10X,58(1H-))
 9048 FORMAT(/10X,57(1H-)/10X,23HLOWDIN ORBITAL (VO(I,J), 5H**2)*,
     *     29HSIGN(VO(I,J))   BETA ORBITALS/10X,57(1H-))
 9068 FORMAT(/10X,41(1H-)/10X,23HLOWDIN ORBITAL (VO(I,J), 5H**2)*,
     *     13HSIGN(VO(I,J))/10X,41(1H-))
 9088 FORMAT(/10X,24(1H-),/10X,24HLOWDIN BOND ORDER MATRIX/10X,24(1H-))
      END
C*MODULE PRPPOP  *DECK MOPOP
      SUBROUTINE MOPOP(S,V,T,IA,L1,L2,L3,QMOAT,OCCNO,NATM,SPI)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION S(L2),V(L1,*),T(L1),IA(L1),QMOAT(NATM,L1),
     *          OCCNO(L1),SPI(L1,*)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /ATLIM/ LIMLOW(MXATM),LIMSUP(MXATM)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NMOGVB,NCONF(MXAO),NHAM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DATA RHF,UHF,ROHF,GVB/8HRHF     ,8HUHF     ,8HROHF    ,8HGVB     /
      DATA RMC,RNONE/8HMCSCF   ,8HNONE    /
      DATA ALPH,BET/8H   ALPHA,8H    BETA/
C
C     ----- MULLIKEN POPULATION FOR EACH OCCUPIED ORBITAL -----
C
      NORB = NA
      IF(SCFTYP.EQ.GVB) NORB=NMOGVB
C
C        THE ROUTINE IS SET UP FOR VARIOUS SCF FUNCTIONS ONLY
C
C     IF(MPLEVL.GT.0)    RETURN
      IF(MPLEVL.GT.0.AND.IPSI.GT.0) RETURN
      IF(CITYP.NE.RNONE) RETURN
      IF(CCTYP.NE.RNONE) RETURN
C
      IPASS = 1
      NBASE = 15
      IF(SCFTYP.EQ.RMC  .OR.  CITYP.NE.RNONE) NBASE=19
C
  140 CONTINUE
      RLAB = ALPH
      IF (IPASS.EQ.2) RLAB = BET
C
C         READ IN THE NATURAL ORBITAL
C
      CALL DAREAD(IDAF,IODA,V,L3,NBASE,0)
C
C         GET OR GENERATE THEIR OCCUPATION NUMBERS
C
      IF(CITYP.NE.RNONE  .OR.  SCFTYP.EQ.RMC) THEN
         CALL DAREAD(IDAF,IODA,OCCNO,L1,NBASE+2,0)
      ELSE
         DO 200 I=1,L1
            IF(SCFTYP.EQ.RHF)  OCCNO(I) = TWO
            IF(SCFTYP.EQ.UHF)  OCCNO(I) = ONE
            IF(SCFTYP.EQ.GVB)  OCCNO(I) = TWO*F(NCONF(I))
            IF(SCFTYP.EQ.ROHF) THEN
               IF(I.LE.NB) THEN
                  OCCNO(I) = TWO
               ELSE
                  OCCNO(I) = ONE
               END IF
            END IF
            IF(I.GT.NORB) OCCNO(I) = ZERO
  200    CONTINUE
      END IF
C
C     ----- TRANSFORM TO SPHERICAL HARMONICS IF NEEDED -----
C     -----    MULTIPLY BY   SPH**(-1)   ON THE LEFT   -----
C
      IF(ISPHER.GE.0) THEN
         DO J=1,NORB
            DO I=1,L1
               SUM=ZERO
               DO K=1,L1
                  SUM=SUM+SPI(I,K)*V(K,J)
               ENDDO
               T(I)=SUM
            ENDDO
            DO I=1,L1
               V(I,J)=T(I)
            ENDDO
         ENDDO
      END IF
C
C     COMPUTE MULLIKEN POP OF EACH AO IN EVERY OCCUPIED MO,
C     OVERWRITING THE MO'S COLUMN WITH THE POPULATION.
C
      IF (ISPHER.GE.0) THEN
      DO I=1,NORB
         DO J=1,L1
            SUM=ZERO
            DO K=1,L1
               JJ=MAX0(J,K)
               KK=MIN0(J,K)
               JK=IA(JJ)+KK
               SUM=SUM+V(K,I)*S(JK)*V(J,I)
            ENDDO
            T(J)=SUM*OCCNO(I)
         ENDDO
         DO J=1,L1
            V(J,I)=T(J)
         ENDDO
      ENDDO
      ELSE
      DO 340 I = 1,NORB
         DO 320 J = 1,L1
            SUM = ZERO
            DO 310 K = 1,L1
               JK = IA(J)+K
               IF (K .GT. J) JK = IA(K)+J
               SUM = SUM+V(K,I)*S(JK)*V(J,I)
  310       CONTINUE
            T(J) = SUM*OCCNO(I)
  320    CONTINUE
         DO 330 J = 1,L1
            V(J,I) = T(J)
  330    CONTINUE
  340 CONTINUE
      END IF
C
C        TREE KILLER OUTPUT
C
      IF((NPRINT.GE.8) .AND. MASWRK) THEN
         WRITE(IW,9000)
         IF(SCFTYP.EQ.UHF) WRITE(IW,9020) RLAB
         CALL PREVNL(V,OCCNO,NORB,L1,L1)
      END IF
C
C     ----- CONDENSE TO ATOMS -----
C
      DO 440 K=1,NORB
         DO 430 IAT=1,NAT
            SUM = ZERO
            IMIN = LIMLOW(IAT)
            IMAX = LIMSUP(IAT)
            DO 420 I=IMIN,IMAX
               SUM = SUM + V(I,K)
  420       CONTINUE
            QMOAT(IAT,K) = SUM
  430    CONTINUE
  440 CONTINUE
C
      IF (MASWRK) WRITE(IW,9010)
      IF((SCFTYP.EQ.UHF) .AND. MASWRK) WRITE(IW,9020) RLAB
      CALL PREVNL(QMOAT,OCCNO,NORB,NAT,NAT)
      IF((SCFTYP.EQ.RMC) .AND. MASWRK) WRITE(IW,9030)
C
      IF(SCFTYP.NE.UHF  .OR.  IPASS.EQ.2) RETURN
      IPASS = IPASS+1
      NBASE = NBASE+4
      NORB = NB
      GO TO 140
C
 9000 FORMAT(/5X,'AO MULLIKEN POPULATIONS IN EACH MOLECULAR ORBITAL')
 9010 FORMAT(/5X,'MULLIKEN ATOMIC POPULATION IN EACH MOLECULAR ORBITAL')
 9020 FORMAT(5X,A8,' ORBITALS')
 9030 FORMAT(1X,'WARNING! MCSCF POPULATIONS SHOWN ABOVE ARE FOR THE',
     *          ' NATURAL ORBITALS.'/
     *       1X,'IGNORE THE ABOVE DATA FOR MCSCF FUNCTIONS WHICH ARE',
     *          ' NOT OF -FORS- TYPE.'/
     *       1X,'THE FOLLOWING POPULATIONS ARE CORRECT FOR ANY MCSCF',
     *          ' WAVEFUNCTION.')
      END
C*MODULE PRPPOP  *DECK MULKEN
      SUBROUTINE MULKEN
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL LAST,DBG,LOTS,BETA,NRMPRT,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA ALPHA,BET,ALL /8H** ALPHA,8H*** BETA,8H**** ALL/
      DATA ROHF/8HROHF    /, UHF/8HUHF     /, GVB/8HGVB     /
      DATA CHECK/8HCHECK   /, RNONE/8HNONE    /
C
      BETA = (SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF)  .AND.
     *       (CITYP.EQ.RNONE)
      DBG = NPRINT .EQ. 3 .AND. MASWRK
      LOTS = NPRINT.GE.8 .AND. MASWRK
      NRMPRT = NPRINT .NE. -5 .AND. MASWRK
      IF (MASWRK) WRITE(IW,9020)
C
      L1 = NUM
      L2 = (L1*(L1+1))/2
      L3 = L1*L1
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
      LOADFM=0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I30 = I10+L2
      I40 = I30+(NAT*(NAT+1))/2
      I50 = I40+L2
      I60 = I50+L1
      I70 = I60+NAT
      I80 = I70+L1
      I90 = I80+NAT
      I41 = I90
      I51 = I41+L2
      I61 = I51+L1
      I71 = I61+NAT
      I72 = I71+L1
      I82 = I72
      I81 = I72+MAX(L1,NAT)
      I91 = I81+NAT
      I92 = I91+NAT*L1
      I93 = I92+L1
      I00 = I93
      I52 = I00+L3
      I53 = I52+L3
C
C     ----- GET MEMORY -----
C
      LASTC = I53
      NEED = LASTC-I10
      CALL GETFM(NEED)
      IF(DBG) WRITE(IW,9000) I10,I30,I40,I50,I60,I70,
     *                       I41,I51,I61,I71,LASTC
      IF (EXETYP .EQ. CHECK) GO TO 300
C
      MAXI = MAX(L1,NAT)
      DO 100 I = 1,MAXI
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     ----- READ IN OVERLAP MATRIX -----
C     ----- IF NEEDED, TRANSFORM TO SPHERICAL BASIS -----
C
C     -SPH- AT X(I00)
C     -S- AT X(I10)
C
      IF (ISPHER.GE.0) THEN
          IFINV=1
          IFTSP=0
          CALL SPHMAT(X(I00),L1,IFINV,IFTSP)
          CALL DAREAD(IDAF,IODA,X(I40),L2,12,0)
          CALL TFTRI(X(I10),X(I40),X(I00),X(I60),L1,L1,L1)
      ELSE
          CALL DAREAD(IDAF,IODA,X(I10),L2,12,0)
      END IF
C
C     ----- CALCULATE MO POPULATIONS -----
C     ----- TRANSFORM TO SPHERICAL HARMONICS IF NEEDED -----
C
C       -SPH**(-1) AT X(I00)
C              -S  AT X(I10)
C
      IF (ISPHER.GE.0) THEN
         IFINV=-1
         IFTSP= 0
         CALL SPHMAT(X(I00),L1,IFINV,IFTSP)
      END IF
      CALL MOPOP(X(I10),X(I40),X(I40+L3),IA,
     *           L1,L2,L3,X(I91),X(I92),NAT,X(I00))
C
C     ----- READ IN DENSITY MATRICES -----
C     ----- TRANSFORM TO SPHERICAL HARMONICS IF NEEDED -----
C
C     -SPH**(-1)(T)- AT X(I00)
C     -DA- AT X(I40)
C     -DB- AT X(I41)
C     SINCE ONLY THE SUM DA+DB IS CORRECT FOR GVB, WE DO ONLY
C     THE TOTAL POPULATION ANALYSIS FOR GVB FUNCTIONS.
C
      IF (ISPHER.GE.0) THEN
          IFINV=-1
          IFTSP= 1
          CALL SPHMAT(X(I00),L1,IFINV,IFTSP)
          CALL WFNDEN(X(I41),X(I52),L2)
          IF(SCFTYP.EQ.GVB) CALL VADD(X(I41),1,X(I52),1,X(I41),1,L2)
          CALL TFTRI(X(I40),X(I41),X(I00),X(I50),L1,L1,L1)
          IF(SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF)
     *       CALL TFTRI(X(I41),X(I52),X(I00),X(I50),L1,L1,L1)
      ELSE
          CALL WFNDEN(X(I40),X(I41),L2)
          IF(SCFTYP.EQ.GVB) CALL VADD(X(I40),1,X(I41),1,X(I40),1,L2)
      END IF
C
      IPASS = 1
      NDAFL = 69
      IDEN = I40
      IGOC = I50
      IGAC = I60
      ILOC = I70
      ILAC = I80
      ELEC=ALL
      NPASS=1
      LAST=.TRUE.
      IF(BETA) ELEC=ALPHA
      IF(BETA) NPASS=2
      IF(BETA) LAST=.FALSE.
C
C     ----- DO A MULLIKEN POPULATION ANALYSIS ----
C           CALCULATE OVERLAP POPULATION
C
  140 CONTINUE
      CALL OVLPOP(X(IDEN),X(I10),L2)
      IF(DBG) WRITE(IW,9035)
      IF(DBG) CALL PRTRIL(X(IDEN),L1)
C
C     ----- CALCULATE TOTAL GROSS POPULATION IN AO*S ----
C
  160 CONTINUE
      IF(SCFTYP.EQ.UHF  .AND.  (LAST.OR.LOTS) .AND. MASWRK)
     *        WRITE(IW,9030) ELEC
      CALL GROSSC(X(IDEN),X(IGOC),IA,L1)
C
C     ----- RECOVER PREVIOUSLY COMPUTED LOWDIN POPULATIONS -----
C     COMBINE "ALPHA" AND "BETA" LOWDIN POPULATIONS FOR GVB
C
      CALL DAREAD(IDAF,IODA,X(ILOC),L1,NDAFL,0)
      IF((SCFTYP.EQ.GVB .AND. CITYP.EQ.RNONE)
     *         .OR.  (BETA.AND.LAST)) THEN
         CALL DAREAD(IDAF,IODA,X(I71),L1,NDAFL+1,0)
         CALL VADD(X(ILOC),1,X(I71),1,X(ILOC),1,L1)
      END IF
C
C     ----- PRINT MULLIKEN AND LOWDIN POPS IN EACH AO -----
C
      IF((NRMPRT.AND.LAST) .OR. LOTS) THEN
         WRITE(IW,9040)
         WRITE(IW,9060) (I,BFLAB(I),X(I-1+IGOC),X(I-1+ILOC),I=1,L1)
      END IF
C
C     ----- COMPRESS POPS FROM ATOMIC ORBITALS TO ATOMS -----
C           FIRST MULLIKEN, AND THEN LOWDIN POPS.
C
      CALL ATPOP(X(IDEN),IA,X(I30),NAT)
      CALL LATPOP(X(ILOC),X(ILAC),NAT)
C
C     ----- PRINT MULLIKEN OVERLAP POPULATIONS -----
C
      IF(NRMPRT  .AND.  LAST) THEN
         IF (MASWRK) WRITE(IW,9035)
         CALL PRTRI(X(I30),NAT)
      END IF
C
C     ----- CALCULATE TOTAL GROSS POPULATION ON ATOMS -----
C
      CALL GROSSC(X(I30),X(IGAC),IA,NAT)
      IF(DBG) THEN
         WRITE(IW,9090)
         WRITE(IW,9100)
      END IF
      IF(LAST .AND. MASWRK) THEN
         WRITE(IW,9090)
         WRITE(IW,9120)
         WRITE(IP,9140)
      END IF
      DO 210 I=1,NAT
         PM = X(I-1+IGAC)
         QM = ZAN(I) - PM
         PL = X(I-1+ILAC)
         QL = ZAN(I) - PL
         IF(DBG) WRITE(IW,9160) I,ANAM(I),BNAM(I),PM,PL
         IF(LAST .AND. MASWRK) WRITE(IW,9180)
     *      I,ANAM(I),BNAM(I),PM,QM,PL,QL
         IF(LAST .AND. MASWRK) WRITE(IP,9200)
     *      ANAM(I),BNAM(I),PM,QM,PL,QL
  210 CONTINUE
C
C     ----- PRINT SHELL POPULATIONS -----
C
      IF(ISPHER.GE.0  .AND.  MASWRK) CALL SHLPOP(X(IGOC))
      IF (LAST) GO TO 300
C
C     ---- SET UP FOR BETA ANALYSIS FOR UHF OR ROHF ----
C
      IF (IPASS .EQ. NPASS) GO TO 220
      IPASS = IPASS+1
      NDAFL = 70
      IDEN = I41
      IGOC = I51
      IGAC = I61
      ILOC = I71
      ILAC = I81
      ELEC = BET
      GO TO 140
C
C     ----- CALCULATE ORBITAL AND ATOMIC SPIN DENSITIES -----
C     THESE ARE ALPHA-BETA POPULATIONS, FOR UHF AND ROHF ONLY
C
  220 CONTINUE
C
C           I50 = I50 - I51
C           I71 = I70 - I71   (COPY OF I71 IN I72)
C           I70 = I70 + I71
C
      CALL VSUB(X(I51),1,X(I50),1,X(I50),1,L1)
      CALL DCOPY(L1,X(I71),1,X(I72),1)
      CALL VSUB(X(I72),1,X(I70),1,X(I71),1,L1)
      CALL VADD(X(I72),1,X(I70),1,X(I70),1,L1)
C
C           I60 = I60 - I61
C           I81 = I80 - I81   (COPY OF I81 IN I82)
C           I80 = I80 + I81
C
      CALL VSUB(X(I61),1,X(I60),1,X(I60),1,NAT)
      CALL DCOPY(NAT,X(I81),1,X(I82),1)
      CALL VSUB(X(I82),1,X(I80),1,X(I81),1,NAT)
      CALL VADD(X(I82),1,X(I80),1,X(I80),1,NAT)
C
C     ----- PRINT THE AO AND ATOMIC SPIN DENSITY -----
C
      IF(LOTS) THEN
         WRITE(IW,9220)
         WRITE(IW,9060) (I,BFLAB(I),X(I-1+I50),X(I-1+I71),I = 1,NUM)
      END IF
      IF (MASWRK) THEN
         WRITE(IW,9240)
         WRITE(IW,9100)
         WRITE(IW,9160) (I,ANAM(I),BNAM(I),
     *                   X(I-1+I60),X(I-1+I81),I=1,NAT)
      END IF
C
C     ---- SET UP FOR ALL ELECTRON ANALYSIS OF ROHF AND UHF ----
C     THE VADD COMBINES THE ALPHA AND BETA AO POP. MATRIX
C
      LAST = .TRUE.
      ELEC = ALL
      CALL VADD(X(I40),1,X(I41),1,X(I40),1,L2)
      NDAFL = 69
      IDEN = I40
      IGOC = I50
      IGAC = I60
      ILOC = I70
      ILAC = I80
      GO TO 160
C
C     ----- RESET MEMORY -----
C
  300 CONTINUE
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(1X,'MEMORY ASSIGNMENT'/1X,'I10, I30, I40, I50,',
     *       ' I60, I70, I41, I51, I61, I71 ='/10I8/8H LAST = ,I8)
 9020 FORMAT(/10X,39(1H-),/10X,31HMULLIKEN AND LOWDIN POPULATION ,
     *     8HANALYSES/10X,39(1H-))
 9030 FORMAT(//10X,'*****',A8,' ELECTRONS ********')
 9035 FORMAT(/10X,'----- MULLIKEN ATOMIC OVERLAP POPULATIONS -----'/
     *        10X,'(OFF-DIAGONAL ELEMENTS NEED TO BE MULTIPLIED BY 2)')
 9040 FORMAT(/15X,'----- POPULATIONS IN EACH AO -----'/
     *        25X,'    MULLIKEN      LOWDIN')
 9060 FORMAT(10X,I5,2X,A8,2F12.5)
 9090 FORMAT(/10X,'TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS')
 9100 FORMAT(7X,'ATOM',9X,'MULL.POP.          ',10X,'LOW.POP.')
 9120 FORMAT(7X,'ATOM',9X,'MULL.POP.    CHARGE',10X,
     *           'LOW.POP.     CHARGE')
 9140 FORMAT(1X,'POPULATION ANALYSIS')
 9160 FORMAT(1X,I4,1X,A8,A2,F12.6,17X,F12.6)
 9180 FORMAT(1X,I4,1X,A8,A2,2F12.6,5X,2F12.6)
 9200 FORMAT(A8,A2,4F10.5)
 9220 FORMAT(/5X,'----- AO SPIN POPULATIONS -----'/
     *        25X,'    MULLIKEN      LOWDIN')
 9240 FORMAT(/5X,'ATOMIC SPIN POPULATION (ALPHA MINUS BETA)')
      END
C*MODULE PRPPOP  *DECK OVLPOP
      SUBROUTINE OVLPOP(A,B,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(N),B(N)
C
      DO 100 I = 1,N
         A(I) = A(I)*B(I)
  100 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK SHLPOP
      SUBROUTINE SHLPOP(AOPOP)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION AOPOP(*)
      DIMENSION SHL(7)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
      PARAMETER (ZERO=0.0D+00)
C
      WRITE(IW,9999)
C
C     ----- ON ENTRY, POPULATION MUST BE IN -SPH.HARM.- BASIS -----
C
      DO IAT=1,NAT
         DO ISHL=1,7
            SHL(ISHL)=ZERO
         ENDDO
         DO ISH=1,NSHELL
            IF(KATOM(ISH).EQ.IAT) THEN
               LIT=KTYPE(ISH)
               MIN=KMIN(ISH)
               MAX=KMAX(ISH)
               LOC=KLOC(ISH)-MIN
C
C                 THIS SUMS THE -S- PART OF AN -L- SHELL TO S POP.
               IF(LIT.EQ.2) THEN
                  IF(MIN.EQ.1) THEN
                     SHL(1)=SHL(1)+AOPOP( 1+LOC)
                     MIN=MIN+1
                  END IF
C
C        --- SUM CONTAMINANT POPULATIONS INTO LOWER -L- TYPES ---
C        THE CONTAMINANT POPULATIONS WILL NOT BE ZERO IF ISPHER=0,
C        BUT ARE EXPECTED TO BE VANISHING IF ISPHER=1 HAS BEEN USED.
C
C                    CARTESIAN D SHELL HAS S CONTAMINANT
               ELSE IF(LIT.EQ.3) THEN
                  SHL(1)=SHL(1)+AOPOP(10+LOC)
                  MAX=MAX-1
C                    CARTESIAN F SHELL HAS P CONTAMINANT
               ELSE IF(LIT.EQ.4) THEN
                  SHL(2)=SHL(2)+AOPOP(18+LOC)
     *                         +AOPOP(19+LOC)+AOPOP(20+LOC)
                  MAX=MAX-3
C                    CARTESIAN G SHELL HAS S,D CONTAMINANTS
               ELSE IF(LIT.EQ.5) THEN
                  SHL(1)=SHL(1)+AOPOP(30+LOC)
                  SHL(3)=SHL(3)+AOPOP(31+LOC)+AOPOP(32+LOC)
     *                         +AOPOP(33+LOC)+AOPOP(34+LOC)
     *                         +AOPOP(35+LOC)
                  MAX=MAX-6
C
C     MWS:     H+I ARE INCOMPLETE HERE.  REST OF THIS ROUTINE=OK,
C     SO AFTER SYMMETRY CODE SYMSLC,SYMORB IS UPDATED, SUCH THAT
C     ROUTINE -SPHMAT- KNOWS HOW TO DO H+I, FIX HERE WILL BE EASY.
C
C                    CARTESIAN H SHELL HAS P,F CONTAMINANTS
               ELSE IF(LIT.EQ.6) THEN
                  WRITE(IW,*) 'SHLPOP: H PARTITION MISSING'
                  CALL ABRT
C                    CARTESIAN I SHELL HAS S,D,G CONTAMINANTS
               ELSE IF(LIT.EQ.7) THEN
                  WRITE(IW,*) 'SHLPOP: I PARTITION MISSING'
                  CALL ABRT
               END IF
C                    SUM OVER THE SPHERICAL HARMONIC'S RANGE
               DO I=MIN,MAX
                  SHL(LIT)=SHL(LIT)+AOPOP(I+LOC)
               ENDDO
            END IF
         ENDDO
C
         SUM=SHL(1)+SHL(2)+SHL(3)+SHL(4)+SHL(5)+SHL(6)+SHL(7)
         WRITE(IW,9998) IAT,ANAM(IAT),BNAM(IAT),(SHL(I),I=1,7),SUM
      ENDDO
      RETURN
C
 9999 FORMAT(/10X,'MULLIKEN SPHERICAL HARMONIC POPULATIONS'/
     *       7X,'ATOM',11X,'S',7X,'P',7X,'D',6X,'F',6X,'G',
     *                         6X,'H',6X,'I',4X,'TOTAL')
 9998 FORMAT(1X,I4,1X,A8,A2,3F8.2,4F7.2,F8.2)
      END
C*MODULE PRPPOP  *DECK SPIND
      SUBROUTINE SPIND(DETDEN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DETDEN
      LOGICAL IANDJ,OUT,NORM,DOUBLE,GOPARR,DSKWRK,MASWRK
C
      DIMENSION XIN(49),YIN(49),ZIN(49),
     *          IX(84),IY(84),IZ(84),JX(84),JY(84),JZ(84),
     *          IJX(784),IJY(784),IJZ(784),SS(784),DIJ(784)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192)
C
      COMMON /FMCOM / X(1)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
      COMMON /XYZPRP/ XP,YP,ZP,PAD(35)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (SQRT3 =1.73205080756887729353D+00)
      PARAMETER (SQRT5 =2.23606797749978969641D+00)
      PARAMETER (SQRT7 =2.64575131106459059050D+00)
      PARAMETER (SQRT9 =3.00000000000000000000D+00)
      PARAMETER (SQRT11=3.31662479035539984911D+00)
C
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     *          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     *          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     *          2, 0, 2, 1, 1,
     *          5, 0, 0, 4, 4, 1, 0, 1, 0, 3,
     *          3, 2, 0, 2, 0, 3, 1, 1, 2, 2,
     *          1,
     *          6, 0, 0, 5, 5, 1, 0, 1, 0, 4,
     *          4, 2, 0, 2, 0, 4, 1, 1, 3, 3,
     *          0, 3, 3, 2, 1, 2, 1, 2/
      DATA IX / 1, 8, 1, 1,15, 1, 1, 8, 8, 1,
     *         22, 1, 1,15,15, 8, 1, 8, 1, 8,
     *         29, 1, 1,22,22, 8, 1, 8, 1,15,
     *         15, 1,15, 8, 8,
     *         36, 1, 1,29,29, 8, 1, 8, 1,22,
     *         22,15, 1,15, 1,22, 8, 8,15,15,
     *          8,
     *         43, 1, 1,36,36, 8, 1, 8, 1,29,
     *         29,15, 1,15, 1,29, 8, 8,22,22,
     *          1,22,22,15, 8,15, 8,15/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     *          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     *          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     *          0, 2, 1, 2, 1,
     *          0, 5, 0, 1, 0, 4, 4, 0, 1, 2,
     *          0, 3, 3, 0, 2, 1, 3, 1, 2, 1,
     *          2,
     *          0, 6, 0, 1, 0, 5, 5, 0, 1, 2,
     *          0, 4, 4, 0, 2, 1, 4, 1, 3, 0,
     *          3, 2, 1, 3, 3, 1, 2, 2/
      DATA IY / 1, 1, 8, 1, 1,15, 1, 8, 1, 8,
     *          1,22, 1, 8, 1,15,15, 1, 8, 8,
     *          1,29, 1, 8, 1,22,22, 1, 8,15,
     *          1,15, 8,15, 8,
     *          1,36, 1, 8, 1,29,29, 1, 8,15,
     *          1,22,22, 1,15, 8,22, 8,15, 8,
     *         15,
     *          1,43, 1, 8, 1,36,36, 1, 8,15,
     *          1,29,29, 1,15, 8,29, 8,22, 1,
     *         22,15, 8,22,22, 8,15,15/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     *          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     *          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     *          2, 2, 1, 1, 2,
     *          0, 0, 5, 0, 1, 0, 1, 4, 4, 0,
     *          2, 0, 2, 3, 3, 1, 1, 3, 1, 2,
     *          2,
     *          0, 0, 6, 0, 1, 0, 1, 5, 5, 0,
     *          2, 0, 2, 4, 4, 1, 1, 4, 0, 3,
     *          3, 1, 2, 1, 2, 3, 3, 2/
      DATA IZ / 1, 1, 1, 8, 1, 1,15, 1, 8, 8,
     *          1, 1,22, 1, 8, 1, 8,15,15, 8,
     *          1, 1,29, 1, 8, 1, 8,22,22, 1,
     *         15,15, 8, 8,15,
     *          1, 1,36, 1, 8, 1, 8,29,29, 1,
     *         15, 1,15,22,22, 8, 8,22, 8,15,
     *         15,
     *          1, 1,43, 1, 8, 1, 8,36,36, 1,
     *         15, 1,15,29,29, 8, 8,29, 1,22,
     *         22, 8,15, 8,15,22,22,15/
C
      DATA RHF/8HRHF     /, GVB/8HGVB     /, RMC/8HMCSCF   /,
     *     ROHF/8HROHF    /, RNONE/8HNONE    /
      DATA CHECK/8HCHECK   /
C
C    SPIN DENSITY FOR RHF OR PP-GVB IS IDENTICALLY ZERO.
C    GAMESS DOES NOT MAINTAIN SEPARATE A+B DENSITY MATRICES
C           FOR GUGA BASED CI AND MCSCF WAVEFUNCTIONS
C    GAMESS INCORRECTLY COMPUTES GVB DENSITY FOR LOW SPIN
C           SITUATIONS SO THAT PA+PB IS OK, BUT PA-PB IS WRONG.
C           HIGH SPIN COUPLING OF UNPAIRED ELECTRONS IS OK..
C    SOME OTHER CASES THAT MIGHT BE EXPECTED TO HAVE A
C    SEPARATE ALPHA AND BETA MATRIX, ACTUALLY DON'T.
C
C    NOTE, IPSI=0 MEANS SCF OR DFT,
C          IPSI=1 MEANS CI,
C          IPSI=2 MEANS MP,
C          IPSI=3 MEANS CC,
C          IPSI=4 MEANS VB
C
      IF(SCFTYP.EQ.RHF)                        RETURN
      IF(SCFTYP.EQ.RMC  .AND.  .NOT.DETDEN)    RETURN
      IF(CITYP.NE.RNONE .AND.  .NOT.DETDEN)    RETURN
      IF(SCFTYP.EQ.GVB  .AND.  MUL.NE.NOPEN+1) RETURN
      IF(SCFTYP.EQ.GVB  .AND.  NOPEN.EQ.0)     RETURN
      IF(SCFTYP.EQ.ROHF .AND.  IPSI.EQ.1)      RETURN
      IF(SCFTYP.EQ.ROHF .AND.  IPSI.EQ.2)      RETURN
      IF(SCFTYP.EQ.ROHF .AND.  IPSI.EQ.3)      RETURN
C
C        ASSUMING THAT WE HAVE NOT EXECUTED THE -IEFBR14- MANEUVER,
C        WE MUST HAVE THE REQUISITE DENSITIES TO EVALUATE THIS PROPERTY.
C
      TOL = RLN10*ITOL
      OUT = NPRINT .EQ. 3 .AND. MASWRK
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
      IF (MASWRK) WRITE(IW,9000)
C
      L1 = NUM
      L2 = (L1*(L1+1))/2
      DUM1=ZERO
      DUM2=ZERO
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
      CALL VALFM(LOADFM)
      LDA  = LOADFM + 1
      LDB  = LDA    + L2
      LINT = LDB    + L2
      LAT  = LINT   + L2
      LATA = LAT    + NAT
      LATB = LATA   + NAT
      LAST = LATB   + NAT
      NEED = LAST-LOADFM-1
      CALL GETFM(NEED)
C
      IF (EXETYP .EQ. CHECK) GO TO 700
C
C     ----- ZERO ATOMIC PROPERTY STORAGE
C
      IF(GOPARR) THEN
         CALL VCLR(X(LAT ),1,NAT)
         CALL VCLR(X(LATA),1,NAT)
         CALL VCLR(X(LATB),1,NAT)
      END IF
C
C     ----- GET DENSITY MATRICES -----
C
      CALL DAREAD(IDAF,IODA,X(LDA),L2,16,0)
      CALL DAREAD(IDAF,IODA,X(LDB),L2,20,0)
C
C     ----- INTIALIZE PARALLEL -----
C     PARALLELIZATION ON ATOMS LEADS TO LOW MESSAGING NEEDS,
C     BUT IS NOT GREAT IF THERE ARE ONLY A FEW ATOMS.
C
      IPCOUNT = ME - 1
C
C     ----- LOOP OVER ATOMS -----
C
      DO 600 IAT = 1,NAT
      XP = C(1,IAT)
      YP = C(2,IAT)
      ZP = C(3,IAT)
C
C     ----- GO PARALLEL! (STATIC LOAD BALANCING) -----
C
      IF (GOPARR) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 600
      END IF
C
C     ----- I SHELL
C
      DO 490 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
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
C
C     ----- J SHELL
C
      DO 480 JJ = 1,II
      J = KATOM(JJ)
      XJ = C(1,J)
      YJ = C(2,J)
      ZJ = C(3,J)
      J1 = KSTART(JJ)
      J2 = J1+KNG(JJ)-1
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      IJ = 0
      MAX = MAXJ
      DO 100 I = MINI,MAXI
      NX = IX(I)
      NY = IY(I)
      NZ = IZ(I)
      IF (IANDJ) MAX = I
      DO 100 J = MINJ,MAX
      IJ = IJ+1
      IJX(IJ) = NX+JX(J)
      IJY(IJ) = NY+JY(J)
      IJZ(IJ) = NZ+JZ(J)
  100 CONTINUE
C
      CALL VCLR(SS,1,IJ)
C
C     ----- I PRIMITIVE
C
      JGMAX = J2
      DO 440 IG = I1,I2
      AI = EX(IG)
      ARRI = AI*RR
      AXI = AI*XI
      AYI = AI*YI
      AZI = AI*ZI
      CSI = CS(IG)
      CPI = CP(IG)
      CDI = CD(IG)
      CFI = CF(IG)
      CGI = CG(IG)
      CHI = CH(IG)
      CII = CI(IG)
C
C     ----- J PRIMITIVE
C
      IF (IANDJ) JGMAX = IG
      DO 430 JG = J1,JGMAX
      AJ = EX(JG)
      AA = AI+AJ
      AA1 = ONE/AA
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ = CF(JG)
      CGJ = CG(JG)
      CHJ = CH(JG)
      CIJ = CI(JG)
C
      AX = (AXI+AJ*XJ)*AA1
      AY = (AYI+AJ*YJ)*AA1
      AZ = (AZI+AJ*ZJ)*AA1
      DUM = AJ*ARRI*AA1+AA*((XP-AX)**2+(YP-AY)**2+(ZP-AZ)**2)
      IF (DUM .GT. TOL) GO TO 430
      FAC = EXP(-DUM)
C
C     ----- DENSITY FACTOR
C
      DOUBLE=IANDJ.AND.IG.NE.JG
      MAX = MAXJ
      NN = 0
      DO 360 I = MINI,MAXI
         IF( I.EQ.1)            DUM1 = CSI*FAC
         IF( I.EQ.2)            DUM1 = CPI*FAC
         IF( I.EQ.5)            DUM1 = CDI*FAC
         IF( I.EQ.8  .AND.NORM) DUM1 = DUM1*SQRT3
C
         IF( I.EQ.11)           DUM1 = CFI*FAC
         IF((I.EQ.14).AND.NORM) DUM1 = DUM1*SQRT5
         IF((I.EQ.20).AND.NORM) DUM1 = DUM1*SQRT3
         IF( I.EQ.21)           DUM1 = CGI*FAC
C
         IF((I.EQ.24).AND.NORM) DUM1 = DUM1*SQRT7
         IF((I.EQ.30).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
         IF((I.EQ.33).AND.NORM) DUM1 = DUM1*SQRT3
C
         IF( I.EQ.36)           DUM1 = CHI*FAC
         IF((I.EQ.39).AND.NORM) DUM1 = DUM1*SQRT9
         IF((I.EQ.45).AND.NORM) DUM1 = DUM1*SQRT7/SQRT3
         IF((I.EQ.51).AND.NORM) DUM1 = DUM1*SQRT3
         IF((I.EQ.54).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
C
         IF( I.EQ.57)           DUM1 = CII*FAC
         IF((I.EQ.60).AND.NORM) DUM1 = DUM1*SQRT11
         IF((I.EQ.66).AND.NORM) DUM1 = DUM1*SQRT9/SQRT3
         IF((I.EQ.72).AND.NORM) DUM1 = DUM1*SQRT3
         IF((I.EQ.75).AND.NORM) DUM1 = DUM1*SQRT7/(SQRT3*SQRT5)
         IF((I.EQ.78).AND.NORM) DUM1 = DUM1*SQRT5
         IF((I.EQ.84).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
C
         IF (IANDJ) MAX = I
         DO 360 J = MINJ,MAX
C
            IF(J.EQ.1) THEN
              DUM2 = DUM1*CSJ
              IF(DOUBLE .AND. I.EQ.1) DUM2 = DUM2 + DUM2
              IF(DOUBLE .AND. I.GT.1) DUM2 = DUM2 + CSI*CPJ*FAC
C
            ELSE IF( J.EQ.2) THEN
              DUM2 = DUM1*CPJ
              IF(DOUBLE) DUM2 = DUM2 + DUM2
C
            ELSE IF( J.EQ.5) THEN
              DUM2 = DUM1*CDJ
              IF(DOUBLE) DUM2 = DUM2 + DUM2
            ELSE IF((J.EQ.8).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF (J.EQ.11) THEN
              DUM2 = DUM1*CFJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.14).AND.NORM) THEN
              DUM2 = DUM2*SQRT5
            ELSE IF((J.EQ.20).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF(J.EQ.21) THEN
              DUM2 = DUM1*CGJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.24).AND.NORM) THEN
              DUM2 = DUM2*SQRT7
            ELSE IF((J.EQ.30).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
            ELSE IF((J.EQ.33).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF( J.EQ.36) THEN
              DUM2 = DUM1*CHJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.39).AND.NORM) THEN
              DUM2 = DUM2*SQRT9
            ELSE IF((J.EQ.45).AND.NORM) THEN
              DUM2 = DUM2*SQRT7/SQRT3
            ELSE IF((J.EQ.51).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
            ELSE IF((J.EQ.54).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
C
            ELSE IF( J.EQ.57) THEN
              DUM2 = DUM1*CIJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.60).AND.NORM) THEN
              DUM2 = DUM2*SQRT11
            ELSE IF((J.EQ.66).AND.NORM) THEN
              DUM2 = DUM2*SQRT9/SQRT3
            ELSE IF((J.EQ.72).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
            ELSE IF((J.EQ.75).AND.NORM) THEN
              DUM2 = DUM2*SQRT7/(SQRT3*SQRT5)
            ELSE IF((J.EQ.78).AND.NORM) THEN
              DUM2 = DUM2*SQRT5
            ELSE IF((J.EQ.84).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
            END IF
C
         NN = NN+1
         DIJ(NN) = DUM2
  360 CONTINUE
C
C     ----- DENSITY INTEGRALS -----
C
      IN = -7
      DO 380 I = 1,LIT
         IN = IN+7
         NI = I
         DO 370 J = 1,LJT
            JN = IN+J
            NJ = J
            CALL DELINT
            XIN(JN) = XINT1
            YIN(JN) = YINT1
            ZIN(JN) = ZINT1
  370    CONTINUE
  380 CONTINUE
C
      DO 400 I = 1,IJ
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         SS(I) = SS(I)+DIJ(I)*XIN(NX)*YIN(NY)*ZIN(NZ)
  400 CONTINUE
C
C          THESE CLOSE PRIMITIVE LOOPS
  430 CONTINUE
  440 CONTINUE
C
C     ----- SET UP DENSITY INTEGRAL MATRIX -----
C
      MAX = MAXJ
      NN = 0
      DO I = MINI,MAXI
         LI = LOCI+I
         IN = (LI*(LI-1))/2
         IF (IANDJ) MAX = I
         DO J = MINJ,MAX
             LJ = LOCJ+J
             JN = LJ+IN
             NN = NN+1
             X(JN-1+LINT) = SS(NN)
         ENDDO
      ENDDO
C
C          THESE CLOSE SHELL LOOPS
  480    CONTINUE
  490 CONTINUE
C
C         DEBUG PRINT OF DELTA FUNCTION INTEGRALS.
C
      IF(OUT) THEN
         WRITE(IW,9010) IAT
         CALL PRTRIL(X(LINT),L1)
      END IF
C
C         COMPUTE THIS ATOM'S SPIN DENSITY.
C         ALL CONSTANTS MULTIPLYING THIS EXPECTATION VALUE ARE DELETED.
C
      X(IAT-1+LATA) = TRACEP(X(LDA),X(LINT),L1)
      X(IAT-1+LATB) = TRACEP(X(LDB),X(LINT),L1)
      X(IAT-1+LAT ) = X(IAT-1+LATA) - X(IAT-1+LATB)
  600 CONTINUE
C
C        END OF ATOM LOOP, OUTPUT RESULTS
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(917,X(LAT ),NAT)
         CALL DDI_GSUMF(918,X(LATA),NAT)
         CALL DDI_GSUMF(919,X(LATB),NAT)
      END IF
C
      IF(MASWRK) THEN
         WRITE(IW,9020)
         DO I = 1,NAT
            IF(ZAN(I).GT.ZERO) WRITE(IW,9030) I,ANAM(I),BNAM(I),
     *                ZAN(I),X(I-1+LAT),X(I-1+LATA),X(I-1+LATB)
         ENDDO
      END IF
C
  700 CONTINUE
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(/10X,41(1H-)/
     *        10X,'ATOMIC SPIN DENSITY AT THE NUCLEUS (A.U.)'/
     *        10X,41(1H-))
 9010 FORMAT(/10X,22(1H-)/10X,'SPIN DENSITY INTEGRALS',10X,' ATOM',I5/
     *        10X,22(1H-))
 9020 FORMAT(28X,'SPIN DENS',4X,'ALPHA DENS',5X,'BETA DENS')
 9030 FORMAT(1X,I4,2X,A8,A2,F6.1,F14.7,2F14.5)
      END
C*MODULE PRPPOP  *DECK SYMTRD
      SUBROUTINE SYMTRD(SH,D,DT,DO,DOD,IA,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION SH(*),D(*),DT(NDIM,*),DO(*),DOD(*),IA(*)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- TRANSFORM THE DENSITY MATRIX TO THE SYMMETRICALLY
C           ORTHOGONAL BASIS SET -----
C
      DO 140 I = 1,N
         DO 120 J = 1,N
            DTEMP = ZERO
            DO 100 K = 1,N
               IK = IA(I)+K
               IF (I .LT. K) IK = IA(K)+I
               JK = IA(J)+K
               IF (J .LT. K) JK = IA(K)+J
               DTEMP = DTEMP+SH(IK)*D(JK)
  100       CONTINUE
            DT(I,J) = DTEMP
  120    CONTINUE
  140 CONTINUE
C
      IJ = 0
      DO 200 I = 1,N
         DO 180 J = 1,I
            IJ = IJ+1
            DTEMP = ZERO
            DO 160 K = 1,N
               JK = IA(J)+K
               IF (J .LT. K) JK = IA(K)+J
               DTEMP = DTEMP+DT(I,K)*SH(JK)
  160       CONTINUE
            DO(IJ) = DTEMP
  180    CONTINUE
  200 CONTINUE
C
      II = 0
      DO 220 I = 1,N
         II = II+I
         DOD(I) = DO(II)
  220 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK SYMTRN
      SUBROUTINE SYMTRN(S,SV,SE,SH,B,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION S(*),SV(NDIM,*),SE(*),SH(*),B(NDIM,9)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- CALCULATE THE SYMMETRIC S **(+1/2) MATRIX
C
C     ----- DIAGONALIZE THE OVERLAP MATRIX -----
C
      IERR=0
      CALL GLDIAG(N,N,NDIM,S,B,SE,SV,IERR,B(1,9))
      IF (IERR .NE. 0) CALL ABRT
C
C     ----- TAKE THE SQUARE ROOT OF THE EIGENVALUES
C
      DO 100 I = 1,N
         SE(I) = SQRT(SE(I))
  100 CONTINUE
C
C     ----- FORM S **(+1/2) -----
C
      IJ = 0
      DO 160 I = 1,N
         DO 140 J = 1,I
            IJ = IJ+1
            SUM = ZERO
            DO 120 K = 1,N
               SUM = SUM+SV(I,K)*SE(K)*SV(J,K)
  120       CONTINUE
            SH(IJ) = SUM
  140    CONTINUE
  160 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK SYMTRV
      SUBROUTINE SYMTRV(SH,V,VO,VO2,T,M,N,NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION SH(*),V(NDIM,M),VO(NDIM,M),VO2(NDIM,M),T(*)
C
      PARAMETER (ZERO=0.0D+00)
C
C      TRANSFORM ORBITALS TO THE SYMMETRICALLY ORTHOGONALIZED BASIS SET
C      THERE ARE N AO-S AND M MO-S IN THE CALCULATION
C
      DO 180 J = 1,M
         IK = 0
         DO 140 I = 1,N
            IM1 = I-1
            DUM = ZERO
            VDUM = V(I,J)
            IF (IM1 .EQ. 0) GO TO 120
               DO 100 K = 1,IM1
                  T(K) = T(K)+SH(IK+K)*VDUM
                  DUM = DUM+SH(IK+K)*V(K,J)
  100          CONTINUE
  120       IK = IK+I
            T(I) = DUM+SH(IK)*VDUM
  140    CONTINUE
         DO 160 I = 1,N
            VO(I,J) = T(I)
  160    CONTINUE
  180 CONTINUE
C
      DO 220 I = 1,N
         DO 200 J = 1,M
            VT = VO(I,J)
            VO2(I,J) = VT*VT
            IF (VT .LT. ZERO) VO2(I,J) = -VO2(I,J)
  200    CONTINUE
  220 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK STNMOM
      SUBROUTINE STNMOM(DMAT,NUM2,CHG,X,Y,Z,XX,YY,ZZ,XY,XZ,YZ,
     *                  XXX,YYY,ZZZ,XXY,XXZ,XYY,YYZ,XZZ,YZZ,XYZ,
     *                  ZANSTN,ZMSSTN,KEQ,GRID,
     *                  TXYZ, LXYZ, W,WR,R, RHO, NSTART,NRAD,NANG,
     *                  MXDIM,NBECKE,NBCKRD,GRGRID,RRGRID,AAGRID,
     *                  SGRID,PPGRID,RADIUS,NFTEFP,FRGNME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MAKEFP
      CHARACTER*8 FRGNME
C
      DIMENSION DMAT(NUM2),CHG(NPTSTN),
     *          X(NPTSTN),Y(NPTSTN),Z(NPTSTN),
     *          XX(NPTSTN),YY(NPTSTN),ZZ(NPTSTN),
     *          XY(NPTSTN),XZ(NPTSTN),YZ(NPTSTN),
     *          XXX(NPTSTN),YYY(NPTSTN),ZZZ(NPTSTN),
     *          XXY(NPTSTN),XXZ(NPTSTN),XYY(NPTSTN),
     *          YYZ(NPTSTN),XZZ(NPTSTN),YZZ(NPTSTN),XYZ(NPTSTN),
     *          ZANSTN(NPTSTN),ZMSSTN(NPTSTN),
     *          KEQ(NPTSTN)
      DIMENSION GRID(4,NPTSTN*(NRAD-1)*NANG),
     *          TXYZ(3,NANG), LXYZ(3,MXDIM), W(NANG),WR(NRAD),R(NRAD),
     *          RHO(NPTSTN*(NRAD-1)*NANG),NSTART(NPTSTN+1),
     *          GRGRID(NPTSTN),RRGRID(NPTSTN,NPTSTN),
     *          AAGRID(NPTSTN,NPTSTN),SGRID(NPTSTN,NPTSTN),
     *          PPGRID(NPTSTN),RADIUS(NPTSTN)
C
      DIMENSION TIMSTR(3)
      DIMENSION X0(49),Y0(49),Z0(49),
     *          X1(49),Y1(49),Z1(49),
     *          X2(49),Y2(49),Z2(49),
     *          X3(49),Y3(49),Z3(49)
      DIMENSION JX(84),JY(84),JZ(84),IX(84),IY(84),IZ(84)
      DIMENSION DIJ(784),IJX(784),IJY(784),IJZ(784)
      DIMENSION XAI(0:5),YAI(0:5),ZAI(0:5),XBJ(0:5),YBJ(0:5),ZBJ(0:5)
C
      DATA XAI,YAI,ZAI,XBJ,YBJ,ZBJ/6*1.0,6*1.0,6*1.0,6*1.0,6*1.0,6*1.0/
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192,
     *           MXFGPT=12000)
C
      LOGICAL IANDJ,NORM,DOUBLE,GOPARR,DSKWRK,MASWRK
      LOGICAL DO_QUADRATURE
C
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /EFPPRT/ NOPRT(MXATM),MIDPRT(MXATM),NUMFFD(MXATM),
     *                INOPRT(2*MXATM),LSTGRP(20,MXATM+1),
     *                ISUM(20,MXFGPT),GRPSUM(20,4),NMDFFD
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
      COMMON /XYZMOM/ TX0,TY0,TZ0,TX1,TY1,TZ1,TX2,TY2,TZ2,TX3,TY3,TZ3,
     *                TX4,TY4,TZ4,T,PXI,PYI,PZI,PXJ,PYJ,PZJ,NI,NJ,CX,
     *                CY,CZ
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (THRESH=1.0D-06,RLN10=2.30258D+00,ETOL=36)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (SQRT3 =1.73205080756887729353D+00)
      PARAMETER (SQRT5 =2.23606797749978969641D+00)
      PARAMETER (SQRT7 =2.64575131106459059050D+00)
      PARAMETER (SQRT9 =3.00000000000000000000D+00)
      PARAMETER (SQRT11=3.31662479035539984911D+00)
C
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     *          3, 0, 0, 2, 2, 1, 0, 1, 0, 1,
     *          4, 0, 0, 3, 3, 1, 0, 1, 0, 2,
     *          2, 0, 2, 1, 1,
     *          5, 0, 0, 4, 4, 1, 0, 1, 0, 3,
     *          3, 2, 0, 2, 0, 3, 1, 1, 2, 2,
     *          1,
     *          6, 0, 0, 5, 5, 1, 0, 1, 0, 4,
     *          4, 2, 0, 2, 0, 4, 1, 1, 3, 3,
     *          0, 3, 3, 2, 1, 2, 1, 2/
      DATA IX / 1, 8, 1, 1,15, 1, 1, 8, 8, 1,
     *         22, 1, 1,15,15, 8, 1, 8, 1, 8,
     *         29, 1, 1,22,22, 8, 1, 8, 1,15,
     *         15, 1,15, 8, 8,
     *         36, 1, 1,29,29, 8, 1, 8, 1,22,
     *         22,15, 1,15, 1,22, 8, 8,15,15,
     *          8,
     *         43, 1, 1,36,36, 8, 1, 8, 1,29,
     *         29,15, 1,15, 1,29, 8, 8,22,22,
     *          1,22,22,15, 8,15, 8,15/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     *          0, 3, 0, 1, 0, 2, 2, 0, 1, 1,
     *          0, 4, 0, 1, 0, 3, 3, 0, 1, 2,
     *          0, 2, 1, 2, 1,
     *          0, 5, 0, 1, 0, 4, 4, 0, 1, 2,
     *          0, 3, 3, 0, 2, 1, 3, 1, 2, 1,
     *          2,
     *          0, 6, 0, 1, 0, 5, 5, 0, 1, 2,
     *          0, 4, 4, 0, 2, 1, 4, 1, 3, 0,
     *          3, 2, 1, 3, 3, 1, 2, 2/
      DATA IY / 1, 1, 8, 1, 1,15, 1, 8, 1, 8,
     *          1,22, 1, 8, 1,15,15, 1, 8, 8,
     *          1,29, 1, 8, 1,22,22, 1, 8,15,
     *          1,15, 8,15, 8,
     *          1,36, 1, 8, 1,29,29, 1, 8,15,
     *          1,22,22, 1,15, 8,22, 8,15, 8,
     *         15,
     *          1,43, 1, 8, 1,36,36, 1, 8,15,
     *          1,29,29, 1,15, 8,29, 8,22, 1,
     *         22,15, 8,22,22, 8,15,15/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     *          0, 0, 3, 0, 1, 0, 1, 2, 2, 1,
     *          0, 0, 4, 0, 1, 0, 1, 3, 3, 0,
     *          2, 2, 1, 1, 2,
     *          0, 0, 5, 0, 1, 0, 1, 4, 4, 0,
     *          2, 0, 2, 3, 3, 1, 1, 3, 1, 2,
     *          2,
     *          0, 0, 6, 0, 1, 0, 1, 5, 5, 0,
     *          2, 0, 2, 4, 4, 1, 1, 4, 0, 3,
     *          3, 1, 2, 1, 2, 3, 3, 2/
      DATA IZ / 1, 1, 1, 8, 1, 1,15, 1, 8, 8,
     *          1, 1,22, 1, 8, 1, 8,15,15, 8,
     *          1, 1,29, 1, 8, 1, 8,22,22, 1,
     *         15,15, 8, 8,15,
     *          1, 1,36, 1, 8, 1, 8,29,29, 1,
     *         15, 1,15,22,22, 8, 8,22, 8,15,
     *         15,
     *          1, 1,43, 1, 8, 1, 8,36,36, 1,
     *         15, 1,15,29,29, 8, 8,29, 1,22,
     *         22, 8,15, 8,15,22,22,15/
C
      DATA MAKEFP/8HMAKEFP  /, CHECK/8HCHECK   /
C
      DO_QUADRATURE = .FALSE.
      NG = NPTSTN*(NRAD-1)*NANG
C
      DO I=1,NPTSTN
         KEQ(I) = 0
         CHG(I) = 0.0D+00
         X(I)   = 0.0D+00
         Y(I)   = 0.0D+00
         Z(I)   = 0.0D+00
         XX(I)  = 0.0D+00
         YY(I)  = 0.0D+00
         ZZ(I)  = 0.0D+00
         XY(I)  = 0.0D+00
         XZ(I)  = 0.0D+00
         YZ(I)  = 0.0D+00
         XXX(I) = 0.0D+00
         YYY(I) = 0.0D+00
         ZZZ(I) = 0.0D+00
         XXY(I) = 0.0D+00
         XXZ(I) = 0.0D+00
         XYY(I) = 0.0D+00
         YYZ(I) = 0.0D+00
         XZZ(I) = 0.0D+00
         YZZ(I) = 0.0D+00
         XYZ(I) = 0.0D+00
      ENDDO
C
      IF(EXETYP.EQ.CHECK) GO TO 600
C
      TOL=RLN10*ITOL
      NORM=NORMF.NE.1.OR.NORMP.NE.1
C
C        PARALLELIZATION IS ON THE INTEGRAL SHELLS AS USUAL,
C        BUT GLOBAL SUM IS OVER THE COMPUTED PROPERTY VALUES.
C
      IPCOUNT = ME-1
C
C     FIRST SHELL LOOP.
C
      DO 9000 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
         LIT=KTYPE(II)
         MINI=KMIN(II)
         MAXI=KMAX(II)
         LOCI=KLOC(II)-MINI
C
C     SECOND SHELL LOOP.
C
         DO 8000 JJ=1,II
C
            IF (GOPARR) THEN
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 8000
            END IF
C
            J=KATOM(JJ)
            XJ=C(1,J)
            YJ=C(2,J)
            ZJ=C(3,J)
            J1=KSTART(JJ)
            J2=J1+KNG(JJ)-1
            LJT=KTYPE(JJ)
            MINJ=KMIN(JJ)
            MAXJ=KMAX(JJ)
            LOCJ=KLOC(JJ)-MINJ
            RR=(XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
            IANDJ=II.EQ.JJ
C
C     PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS.
C
            IJ=0
            MAX=MAXJ
            DO I=MINI,MAXI
               NX=IX(I)
               NY=IY(I)
               NZ=IZ(I)
               IF(IANDJ) MAX=I
               DO J=MINJ,MAX
                  IJ=IJ+1
                  IJX(IJ)=NX+JX(J)
                  IJY(IJ)=NY+JY(J)
                  IJZ(IJ)=NZ+JZ(J)
               ENDDO
            ENDDO
C
C     I PRIMITIVE.
C
            JGMAX=J2
            DO 7000 IG=I1,I2
               AI=EX(IG)
               ARRI=AI*RR
               AXI=AI*XI
               AYI=AI*YI
               AZI=AI*ZI
               CSI=CS(IG)
               CPI=CP(IG)
               CDI=CD(IG)
               CFI=CF(IG)
               CGI=CG(IG)
               CHI=CH(IG)
               CII=CI(IG)
C
C     J PRIMITIVE.
C
               IF(IANDJ) JGMAX=IG
               DO 6000 JG=J1,JGMAX
                  AJ=EX(JG)
                  AA=AI+AJ
                  DUM=AJ*ARRI/AA
                  IF(DUM.GT.TOL) GO TO 6000
                  FAC= EXP(-DUM)
                  CSJ=CS(JG)
                  CPJ=CP(JG)
                  CDJ=CD(JG)
                  CFJ=CF(JG)
                  CGJ=CG(JG)
                  CHJ=CH(JG)
                  CIJ=CI(JG)
C
                  AX=(AXI+AJ*XJ)/AA
                  AY=(AYI+AJ*YJ)/AA
                  AZ=(AZI+AJ*ZJ)/AA
C
C     DENSITY FACTOR.  MORE STANDARD STUFF.
C
      DOUBLE=IANDJ.AND.IG.NE.JG
      MAX=MAXJ
      NN=0
      DO 310 I=MINI,MAXI
         M=LOCI+I
         M2 = (M*(M-1))/2
         IF( I.EQ.1)            DUM1 = CSI*FAC
         IF( I.EQ.2)            DUM1 = CPI*FAC
         IF( I.EQ.5)            DUM1 = CDI*FAC
         IF( I.EQ.8  .AND.NORM) DUM1 = DUM1*SQRT3
C
         IF( I.EQ.11)           DUM1 = CFI*FAC
         IF((I.EQ.14).AND.NORM) DUM1 = DUM1*SQRT5
         IF((I.EQ.20).AND.NORM) DUM1 = DUM1*SQRT3
         IF( I.EQ.21)           DUM1 = CGI*FAC
C
         IF((I.EQ.24).AND.NORM) DUM1 = DUM1*SQRT7
         IF((I.EQ.30).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
         IF((I.EQ.33).AND.NORM) DUM1 = DUM1*SQRT3
C
         IF( I.EQ.36)           DUM1 = CHI*FAC
         IF((I.EQ.39).AND.NORM) DUM1 = DUM1*SQRT9
         IF((I.EQ.45).AND.NORM) DUM1 = DUM1*SQRT7/SQRT3
         IF((I.EQ.51).AND.NORM) DUM1 = DUM1*SQRT3
         IF((I.EQ.54).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
C
         IF( I.EQ.57)           DUM1 = CII*FAC
         IF((I.EQ.60).AND.NORM) DUM1 = DUM1*SQRT11
         IF((I.EQ.66).AND.NORM) DUM1 = DUM1*SQRT9/SQRT3
         IF((I.EQ.72).AND.NORM) DUM1 = DUM1*SQRT3
         IF((I.EQ.75).AND.NORM) DUM1 = DUM1*SQRT7/(SQRT3*SQRT5)
         IF((I.EQ.78).AND.NORM) DUM1 = DUM1*SQRT5
         IF((I.EQ.84).AND.NORM) DUM1 = DUM1*SQRT5/SQRT3
C
         IF(IANDJ) MAX=I
         DO 300 J=MINJ,MAX
C
            IF(J.EQ.1) THEN
              DUM2 = DUM1*CSJ
              IF(DOUBLE .AND. I.EQ.1) DUM2 = DUM2 + DUM2
              IF(DOUBLE .AND. I.GT.1) DUM2 = DUM2 + CSI*CPJ*FAC
C
            ELSE IF( J.EQ.2) THEN
              DUM2 = DUM1*CPJ
              IF(DOUBLE) DUM2 = DUM2 + DUM2
C
            ELSE IF( J.EQ.5) THEN
              DUM2 = DUM1*CDJ
              IF(DOUBLE) DUM2 = DUM2 + DUM2
            ELSE IF((J.EQ.8).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF (J.EQ.11) THEN
              DUM2 = DUM1*CFJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.14).AND.NORM) THEN
              DUM2 = DUM2*SQRT5
            ELSE IF((J.EQ.20).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF(J.EQ.21) THEN
              DUM2 = DUM1*CGJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.24).AND.NORM) THEN
              DUM2 = DUM2*SQRT7
            ELSE IF((J.EQ.30).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
            ELSE IF((J.EQ.33).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
C
            ELSE IF( J.EQ.36) THEN
              DUM2 = DUM1*CHJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.39).AND.NORM) THEN
              DUM2 = DUM2*SQRT9
            ELSE IF((J.EQ.45).AND.NORM) THEN
              DUM2 = DUM2*SQRT7/SQRT3
            ELSE IF((J.EQ.51).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
            ELSE IF((J.EQ.54).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
C
            ELSE IF( J.EQ.57) THEN
              DUM2 = DUM1*CIJ
              IF (DOUBLE) DUM2 = DUM2+DUM2
            ELSE IF((J.EQ.60).AND.NORM) THEN
              DUM2 = DUM2*SQRT11
            ELSE IF((J.EQ.66).AND.NORM) THEN
              DUM2 = DUM2*SQRT9/SQRT3
            ELSE IF((J.EQ.72).AND.NORM) THEN
              DUM2 = DUM2*SQRT3
            ELSE IF((J.EQ.75).AND.NORM) THEN
              DUM2 = DUM2*SQRT7/(SQRT3*SQRT5)
            ELSE IF((J.EQ.78).AND.NORM) THEN
              DUM2 = DUM2*SQRT5
            ELSE IF((J.EQ.84).AND.NORM) THEN
              DUM2 = DUM2*SQRT5/SQRT3
            END IF
C
            NN = NN+1
            DIJ(NN)=DMAT(M2+LOCJ+J)
            DIJ(NN) = DIJ(NN)*DUM2
  300    CONTINUE
  310 CONTINUE
C
C  LVS 2005
C  SWITCH BETWEEN ORIGINAL AND MOFIFIED STONE PROCEDURES
C  IF EXPONENTS ARE LARGE ENOUGH (NOT DIFFUSE) THEN CONTINUE WITH
C  USUAL STONE PROCEDURE, OTHERWISE DO NUMERICAL INTEGRATION
C  OVER GRID
C
                  IF (AA .GT. BIGEXP) THEN
C
C     ALL STANDARD STUFF TO THIS POINT.
C     FIND MINIMAL DISTANCE OF GAUSSIAN PRODUCT TO AN EXPANSION POINT.
C
                  DIST=1.0D+22
                  DO NPT=1,NPTSTN
                     TEST=(AX-STNPNT(2,NPT))**2
     *                   +(AY-STNPNT(3,NPT))**2
     *                   +(AZ-STNPNT(4,NPT))**2
                     IF(TEST.LT.DIST) DIST=TEST
                  ENDDO
C
C     FIND ALL EQUIVALENT MINIMALLY DISTANT EXPANSION POINTS.
C
                  IEQ=0
                  DO 67 NPT=1,NPTSTN
                     TEST=(AX-STNPNT(2,NPT))**2
     *                   +(AY-STNPNT(3,NPT))**2
     *                   +(AZ-STNPNT(4,NPT))**2
                     IF(TEST-DIST.GT.THRESH) GO TO 67
                     IEQ=IEQ+1
                     KEQ(IEQ)=NPT
   67             CONTINUE
                  IF(IEQ.EQ.0) GO TO 6000
                  CEQ=ONE/IEQ
C
C
C     CALCULATE THE MOMENT INTEGRALS.
C
                  T=SQRT(AA)
                  PXI=XI-AX
                  PYI=YI-AY
                  PZI=ZI-AZ
                  PXJ=XJ-AX
                  PYJ=YJ-AY
                  PZJ=ZJ-AZ
C
C     LOOP OVER EXPANSION POINT(S) FOR THESE PRIMITIVE PRODUCTS.
C
                  DO 500 NEQ=1,IEQ
                     IQ=KEQ(NEQ)
                     CX=STNPNT(2,IQ)-AX
                     CY=STNPNT(3,IQ)-AY
                     CZ=STNPNT(4,IQ)-AZ
                     IN=-7
                     DO 350 I=1,LIT
                        IN=IN+7
                        NI=I
                        DO 340 J=1,LJT
                           JN=IN+J
                           NJ=J
C
C     GAUSS-HERMITE QUADRATURE AROUND EXPANSION POINT.  WALT'S OLD
C     VERSION DID THE INTEGRALS AROUND THE GAUSSIAN PRODUCT CENTER
C     AND THEN SHIFTED MOMENTS TO THE EXPANSION POINT(S).
C
                           CALL STNXYZ
C
                           X0(JN)=TX0/T
                           Y0(JN)=TY0/T
                           Z0(JN)=TZ0/T
                           X1(JN)=TX1/T
                           Y1(JN)=TY1/T
                           Z1(JN)=TZ1/T
                           X2(JN)=TX2/T
                           Y2(JN)=TY2/T
                           Z2(JN)=TZ2/T
                           X3(JN)=TX3/T
                           Y3(JN)=TY3/T
                           Z3(JN)=TZ3/T
  340                   CONTINUE
  350                CONTINUE
C
                     MAX=MAXJ
                     NN=0
                     DO 490 I=MINI,MAXI
                        LI=LOCI+I
                        IN=(LI*(LI-1))/2
                        IF(IANDJ) MAX=I
                        DO 480 J=MINJ,MAX
                           LJ=LOCJ+J
                           JN=LJ+IN
                           NN=NN+1
                           NX=IJX(NN)
                           NY=IJY(NN)
                           NZ=IJZ(NN)
C
C     CONTRIBUTIONS TO MOMENTS FROM CURRENT PRIMITIVE PRODUCT.
C
                           DM     = -DIJ(NN)*CEQ
                           FCHG   = DM*X0(NX)*Y0(NY)*Z0(NZ)
                           CHG(IQ)= CHG(IQ)+FCHG
C
                           X(IQ)  =X(IQ)   + DM*X1(NX)*Y0(NY)*Z0(NZ)
                           Y(IQ)  =Y(IQ)   + DM*X0(NX)*Y1(NY)*Z0(NZ)
                           Z(IQ)  =Z(IQ)   + DM*X0(NX)*Y0(NY)*Z1(NZ)
C
                           XX(IQ) =XX(IQ)  + DM*X2(NX)*Y0(NY)*Z0(NZ)
                           YY(IQ) =YY(IQ)  + DM*X0(NX)*Y2(NY)*Z0(NZ)
                           ZZ(IQ) =ZZ(IQ)  + DM*X0(NX)*Y0(NY)*Z2(NZ)
                           XY(IQ) =XY(IQ)  + DM*X1(NX)*Y1(NY)*Z0(NZ)
                           XZ(IQ) =XZ(IQ)  + DM*X1(NX)*Y0(NY)*Z1(NZ)
                           YZ(IQ) =YZ(IQ)  + DM*X0(NX)*Y1(NY)*Z1(NZ)
C
                           XXX(IQ)=XXX(IQ) + DM*X3(NX)*Y0(NY)*Z0(NZ)
                           YYY(IQ)=YYY(IQ) + DM*X0(NX)*Y3(NY)*Z0(NZ)
                           ZZZ(IQ)=ZZZ(IQ) + DM*X0(NX)*Y0(NY)*Z3(NZ)
                           XXY(IQ)=XXY(IQ) + DM*X2(NX)*Y1(NY)*Z0(NZ)
                           XXZ(IQ)=XXZ(IQ) + DM*X2(NX)*Y0(NY)*Z1(NZ)
                           XYY(IQ)=XYY(IQ) + DM*X1(NX)*Y2(NY)*Z0(NZ)
                           YYZ(IQ)=YYZ(IQ) + DM*X0(NX)*Y2(NY)*Z1(NZ)
                           XZZ(IQ)=XZZ(IQ) + DM*X1(NX)*Y0(NY)*Z2(NZ)
                           YZZ(IQ)=YZZ(IQ) + DM*X0(NX)*Y1(NY)*Z2(NZ)
                           XYZ(IQ)=XYZ(IQ) + DM*X1(NX)*Y1(NY)*Z1(NZ)
  480                   CONTINUE
  490                CONTINUE
  500             CONTINUE
C
C           DO NUMERICAL INTEGRATION HERE
                  ELSE
                      IF (.NOT. DO_QUADRATURE) THEN
                        CALL DOGRID(ZANSTN,GRID,TXYZ,
     *                              LXYZ,W,WR,R,NSTART,
     *                              NRAD,NANG,MXDIM,NBECKE,NBCKRD,
     *                              GRGRID,RRGRID,AAGRID,SGRID,PPGRID,
     *                              RADIUS)
C
C              CLEAR GRID OF DENSITY VALUES
                        CALL VCLR(RHO,1,NG)
                        DO_QUADRATURE=.TRUE.
                      END IF
C
                      DO K=1,NG
                        XK=GRID(1,K)
                        YK=GRID(2,K)
                        ZK=GRID(3,K)
                        XA=XK-XI
                        YA=YK-YI
                        ZA=ZK-ZI
                        XB=XK-XJ
                        YB=YK-YJ
                        ZB=ZK-ZJ
                        E=AA*((XK-AX)**2+(YK-AY)**2+(ZK-AZ)**2)
                        IF (E .GT. ETOL) GOTO 220
                        E=EXP(-E)
                        CSI=CS(IG)
                        DO L=1,LIT
                          XAI(L)=XA*XAI(L-1)
                          YAI(L)=YA*YAI(L-1)
                          ZAI(L)=ZA*ZAI(L-1)
                        END DO
                        MAX = MAXJ
                        NN=0
                        DO IA=MINI,MAXI
C                         IF (IA .GT. 1 .AND. LIT .EQ. 1) CSI=CP(IG)
C                         IF (IA .GT. 4) CSI=CD(IG)
C                         CSJ=CS(JG)
                          DO L=1,LJT
                            XBJ(L)=XB*XBJ(L-1)
                            YBJ(L)=YB*YBJ(L-1)
                            ZBJ(L)=ZB*ZBJ(L-1)
                          END DO
                         IF(IANDJ) MAX=IA
                          DO JB=MINJ,MAX
                            NN = NN+1
                            RHO(K)=RHO(K)-E*DIJ(NN)*
     *                             XAI(JX(IA))*YAI(JY(IA))*ZAI(JZ(IA))*
     *                             XBJ(JX(JB))*YBJ(JY(JB))*ZBJ(JZ(JB))
                          END DO
                        END DO
 220                    CONTINUE
                      END DO
                    END IF
C
C  NOW WE HAVE ADDED IN THE OVERLAP CHARGE DENSITY FOR THIS PAIR OF
C  PRIMITIVES
C
C                       END LOOPS OVER EXPANSION POINTS
 6000          CONTINUE
 7000       CONTINUE
C                       END LOOPS OVER GAUSSIAN PRIMITIVES
 8000    CONTINUE
 9000 CONTINUE
C
C
      IF (DO_QUADRATURE) THEN
C  LOOP OVER SITES TO EVALUATE MULTIPOLE CONTRIBUTIONS FROM CHARGE
C  DENSITY
         X0(1)=ONE
         Y0(1)=ONE
         Z0(1)=ONE
         DO I=1,NPTSTN
            XI=STNPNT(2,I)
            YI=STNPNT(3,I)
            ZI=STNPNT(4,I)
C
C  WORK BACKWARDS (MOST DISTANT POINTS FIRST) TO AVOID LOSS OF
C  SIGNIFICANCE.
            KIN = NSTART(I+1)-1
            KF  = NSTART(I)
            DO K=KIN,KF,-1
C
                  XK=GRID(1,K)-XI
                  YK=GRID(2,K)-YI
                  ZK=GRID(3,K)-ZI
C
                  X1(1) = X0(1)*XK
                  Y1(1) = Y0(1)*YK
                  Z1(1) = Z0(1)*ZK
                  X2(1) = X1(1)*XK
                  Y2(1) = Y1(1)*YK
                  Z2(1) = Z1(1)*ZK
                  X3(1) = X2(1)*XK
                  Y3(1) = Y2(1)*YK
                  Z3(1) = Z2(1)*ZK
C
                  NX = 1
                  NY = 1
                  NZ = 1
                  IQ = I
C
                  DM     = GRID(4,K)*RHO(K)
                  CHG(IQ)= CHG(IQ)+ DM*X0(NX)*Y0(NY)*Z0(NZ)
C
                  X(IQ)  =X(IQ)   + DM*X1(NX)*Y0(NY)*Z0(NZ)
                  Y(IQ)  =Y(IQ)   + DM*X0(NX)*Y1(NY)*Z0(NZ)
                  Z(IQ)  =Z(IQ)   + DM*X0(NX)*Y0(NY)*Z1(NZ)
C
                  XX(IQ) =XX(IQ)  + DM*X2(NX)*Y0(NY)*Z0(NZ)
                  YY(IQ) =YY(IQ)  + DM*X0(NX)*Y2(NY)*Z0(NZ)
                  ZZ(IQ) =ZZ(IQ)  + DM*X0(NX)*Y0(NY)*Z2(NZ)
                  XY(IQ) =XY(IQ)  + DM*X1(NX)*Y1(NY)*Z0(NZ)
                  XZ(IQ) =XZ(IQ)  + DM*X1(NX)*Y0(NY)*Z1(NZ)
                  YZ(IQ) =YZ(IQ)  + DM*X0(NX)*Y1(NY)*Z1(NZ)
C
                  XXX(IQ)=XXX(IQ) + DM*X3(NX)*Y0(NY)*Z0(NZ)
                  YYY(IQ)=YYY(IQ) + DM*X0(NX)*Y3(NY)*Z0(NZ)
                  ZZZ(IQ)=ZZZ(IQ) + DM*X0(NX)*Y0(NY)*Z3(NZ)
                  XXY(IQ)=XXY(IQ) + DM*X2(NX)*Y1(NY)*Z0(NZ)
                  XXZ(IQ)=XXZ(IQ) + DM*X2(NX)*Y0(NY)*Z1(NZ)
                  XYY(IQ)=XYY(IQ) + DM*X1(NX)*Y2(NY)*Z0(NZ)
                  YYZ(IQ)=YYZ(IQ) + DM*X0(NX)*Y2(NY)*Z1(NZ)
                  XZZ(IQ)=XZZ(IQ) + DM*X1(NX)*Y0(NY)*Z2(NZ)
                  YZZ(IQ)=YZZ(IQ) + DM*X0(NX)*Y1(NY)*Z2(NZ)
                  XYZ(IQ)=XYZ(IQ) + DM*X1(NX)*Y1(NY)*Z1(NZ)
C
C 70            CONTINUE
            END DO
         END DO
      END IF
C
      CHARGE = ZERO
      DO I=1,NPTSTN
        CHARGE = CHARGE + CHG(I)
      END DO
      IF(MASWRK) WRITE(IW,9020) CHARGE
 9020 FORMAT(1X,'TOTAL ELECTRON CHARGE =',F18.10)
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(901,CHG,NPTSTN)
         CALL DDI_GSUMF(902,X  ,NPTSTN)
         CALL DDI_GSUMF(903,Y  ,NPTSTN)
         CALL DDI_GSUMF(904,Z  ,NPTSTN)
         CALL DDI_GSUMF(905,XX ,NPTSTN)
         CALL DDI_GSUMF(906,YY ,NPTSTN)
         CALL DDI_GSUMF(907,ZZ ,NPTSTN)
         CALL DDI_GSUMF(908,XY ,NPTSTN)
         CALL DDI_GSUMF(909,XZ ,NPTSTN)
         CALL DDI_GSUMF(910,YZ ,NPTSTN)
         CALL DDI_GSUMF(911,XXX,NPTSTN)
         CALL DDI_GSUMF(912,YYY,NPTSTN)
         CALL DDI_GSUMF(913,ZZZ,NPTSTN)
         CALL DDI_GSUMF(914,XXY,NPTSTN)
         CALL DDI_GSUMF(915,XXZ,NPTSTN)
         CALL DDI_GSUMF(916,XYY,NPTSTN)
         CALL DDI_GSUMF(917,YYZ,NPTSTN)
         CALL DDI_GSUMF(918,XZZ,NPTSTN)
         CALL DDI_GSUMF(919,YZZ,NPTSTN)
         CALL DDI_GSUMF(920,XYZ,NPTSTN)
      END IF
C                       END LOOPS OVER BASIS SET SHELLS
C
C          MODIFY THE MONOPOLE AND DIPOLE AT THE ATOM OR THE MIDPOINT
C          DEFINED IN LSTGRP(*,MXATM+1) SO THE
C          EFP HAS AN INTEGER CHARGE
C
      IF(NMDFFD.GT.0) THEN
        DO NMD = 1, NMDFFD
          GRPSUM(NMD,1)=0.0D+00
          DO I = 1, NPTSTN
            IF(ISUM(NMD,I).GT.0) THEN
              GRPSUM(NMD,1)=GRPSUM(NMD,1)+CHG(I)+ZANSTN(I)
            END IF
          ENDDO
          CHG(LSTGRP(NMD,MXATM+1)) =
     *           CHG(LSTGRP(NMD,MXATM+1)) + GRPSUM(NMD,1)
          GRPSUM(NMD,2)=GRPSUM(NMD,1)*0.5D+00*
     *                  (GRPSUM(NMD,2)-STNPNT(2,LSTGRP(NMD,MXATM+1)))
          GRPSUM(NMD,3)=GRPSUM(NMD,1)*0.5D+00*
     *                  (GRPSUM(NMD,3)-STNPNT(3,LSTGRP(NMD,MXATM+1)))
          GRPSUM(NMD,4)=GRPSUM(NMD,1)*0.5D+00*
     *                  (GRPSUM(NMD,4)-STNPNT(4,LSTGRP(NMD,MXATM+1)))
          X(LSTGRP(NMD,MXATM+1))=X(LSTGRP(NMD,MXATM+1))+GRPSUM(NMD,2)
          Y(LSTGRP(NMD,MXATM+1))=Y(LSTGRP(NMD,MXATM+1))+GRPSUM(NMD,3)
          Z(LSTGRP(NMD,MXATM+1))=Z(LSTGRP(NMD,MXATM+1))+GRPSUM(NMD,4)
          IF(MASWRK) WRITE(IW,*)' '
          IF(MASWRK) WRITE(IW,'(A/,1X,A/,A,2X,F14.10)')
     *' EFP/FORCE FIELD CONNECTION POINT  ',
     * STNPNT(1,LSTGRP(NMD,MXATM+1)),
     *' COLLECTED A MONOPOLE :',GRPSUM(NMD,1)
C
          IF(MASWRK) WRITE(IW,'(A,3(2X,F14.10),1X/,A)')
     *' AND A DIPOLE (X-Y-Z) :',(GRPSUM(NMD,1+LC),LC=1,3),
     *' FROM POINTS:'
C
          DO I = 1, NPTSTN
            IF(ISUM(NMD,I).GT.0  .AND. MASWRK)
     *        WRITE(IW,'(1X,A)') STNPNT(1,I)
          ENDDO
          IF(MASWRK)WRITE(IW,*)' '
        ENDDO
      END IF
C
C     ----- PRINT AND PUNCH THE RESULTS -----
C
  600 CONTINUE
C
      DO I=1,NPTSTN
        IF(INOPRT(I).GT.0 .AND. MASWRK) WRITE(IW,'(A,1X,A,1X,A)')
     *     ' POINT',STNPNT(1,I),' WILL NOT BE PRINTED/PUNCHED OUT'
      ENDDO
C
C              IF IMGABI, NO PRINTOUT
C
      IF(IPCMIT.EQ.1 .AND. IMGABI.EQ.1) GOTO 9950
C
C     A SECOND COPY OF THE "EFFECTIVE FRAGMENT INPUT FILE" BEING
C     WRITTEN TO XXX.EFP FOR KEEPS, IS ALSO WRITTEN TO XXX.IRC,
C     TO BE READ BY THE SCREENING CODE.  THE EXTRA COPY ON -IDMP-
C     WILL BE DELETED UPON ITS CLOSE BY THE SCREENING CODE, AS THIS
C     FILE IS JUST A KLUDGY WAY TO PASS STONE'S ANALYSIS TO IT.
C
C     THE SIGN OF -IDMP- TELLS IF WE ARE IN AN EFP GENERATING RUN!
C
      IF(RUNTYP.EQ.MAKEFP) THEN
         IDMP=4
         CALL SEQOPN(IDMP, 'IRCDATA',    'NEW',.FALSE.,'FORMATTED')
         CALL SEQOPN(NFTEFP,'MAKEFP','UNKNOWN',.FALSE.,'FORMATTED')
      ELSE
         IDMP=-1
      END IF
C
      IF (MASWRK  .AND.  INOPRT(I).LT.1) THEN
C
      WRITE(IW,9199)
 9199 FORMAT(/'          ---------------------------------------------',
     1       /'          DISTRIBUTED MULTIPOLAR ANALYSIS (A. J. STONE)',
     1       /'          ---------------------------------------------')
C                                     MONOPOLES
      WRITE(IW,9200)
      CALL TMDATE(TIMSTR)
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9050) FRGNME,TIMSTR,FRGNME,FRGNME,SCFTYP
         WRITE(  IDMP,9050) FRGNME,TIMSTR,FRGNME,FRGNME,SCFTYP
         WRITE(NFTEFP,9201)
         WRITE(  IDMP,9201)
      END IF
      DO I=1,NAT
         WRITE(IW,9211) I,ANAM(I),ZAN(I),(C(J,I),J=1,3)
      ENDDO
      DO I=1,NPTSTN
         WRITE(IW,9210) I+NAT,STNPNT(1,I),CHG(I),(STNPNT(J,I),J=2,4)
         IF(IDMP.GT.0) THEN
            WRITE(NFTEFP,9214) (STNPNT(J,I),J=1,4),ZMSSTN(I),ZANSTN(I)
            WRITE(  IDMP,9214) (STNPNT(J,I),J=1,4),ZMSSTN(I),ZANSTN(I)
         END IF
      ENDDO
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9206)
         WRITE(  IDMP,9206)
         WRITE(NFTEFP,9207)
         WRITE(  IDMP,9207)
         DO I=1,NPTSTN
            WRITE(NFTEFP,9209) STNPNT(1,I),CHG(I),ZANSTN(I)
            WRITE(  IDMP,9209) STNPNT(1,I),CHG(I),ZANSTN(I)
         ENDDO
         WRITE(NFTEFP,9206)
         WRITE(  IDMP,9206)
      END IF
 9200 FORMAT(/1X,'NET CHARGES AT POINTS                    COORDINATES',
     1       /1X,'---------------------                    -----------',
     1       /1X,'   #   NAME        CHARGE          X           Y',
     1        1X,'          Z')
 9050 FORMAT(10X,'RUNTYP=MAKEFP EFFECTIVE FRAGMENT POTENTIAL DATA',
     *           ' FOLLOWS...'/
     *       10X,A8,'EFP GENERATED AT ',3A8/
     *       ' $',A8/
     *       'EFP DATA FOR ',A8,' SCFTYP=',A8,
     *       '... GENERATED WITH BASIS SET=XXX')
 9201 FORMAT(' COORDINATES (BOHR)')
 9206 FORMAT(' STOP')
 9207 FORMAT(' MONOPOLES ')
 9209 FORMAT(A8,F15.10,F10.5)
 9210 FORMAT(1X,I4,3X,A8,2X,F9.5,1X,3(3X,F9.5))
 9211 FORMAT(1X,I4,3X,'Z',A8,1X,F9.5,1X,3(3X,F9.5))
 9214 FORMAT(A8,3F15.10,F12.7,F5.1)
C                                     DIPOLES
      WRITE(IW,9220)
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9221)
         WRITE(   IDMP,9221)
      END IF
      DO I=1,NPTSTN
         WRITE(IW,9230) STNPNT(1,I),X(I),Y(I),Z(I)
         IF(IDMP.GT.0) THEN
            WRITE(NFTEFP,9236) STNPNT(1,I),X(I),Y(I),Z(I)
            WRITE(  IDMP,9236) STNPNT(1,I),X(I),Y(I),Z(I)
         END IF
      ENDDO
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9206)
         WRITE(  IDMP,9206)
      END IF
 9220 FORMAT(/1X,'THE DISTRIBUTED MULTIPOLE ANALYSIS IS FOR THE',
     *           ' -ELECTRONIC- PART ONLY'/
     *        1X,'OF THE FOLLOWING MOMENTS.  ALL QUANTITIES ARE IN',
     *           ' ATOMIC UNITS'///
     *        1X,'FIRST MOMENTS AT POINTS'/
     *        1X,'-----------------------'/
     *        1X,'NAME',15X,'X',13X,'Y',13X,'Z')
 9221 FORMAT(' DIPOLES ')
 9230 FORMAT(1X,A8,3(2X,F12.5))
 9236 FORMAT(A8,3(1X,F15.10))
C                                     QUADRUPOLES
      WRITE(IW,9240)
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9241)
         WRITE(  IDMP,9241)
      END IF
      DO I=1,NPTSTN
            WRITE(IW,9250)     STNPNT(1,I),
     *                         XX(I),YY(I),ZZ(I),XY(I),XZ(I),YZ(I)
         IF(IDMP.GT.0) THEN
            WRITE(NFTEFP,9256) STNPNT(1,I),
     *                         XX(I),YY(I),ZZ(I),XY(I),XZ(I),YZ(I)
            WRITE(  IDMP,9256) STNPNT(1,I),
     *                         XX(I),YY(I),ZZ(I),XY(I),XZ(I),YZ(I)
         END IF
      ENDDO
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9206)
         WRITE(  IDMP,9206)
      END IF
 9240 FORMAT(/1X,'SECOND MOMENTS AT POINTS',
     1       /1X,'------------------------',
     1       /1X,'NAME',11X,'XX',9X,'YY',9X,'ZZ',9X,'XY',9X,
     1                      'XZ',9X,'YZ')
 9241 FORMAT(' QUADRUPOLES ')
 9250 FORMAT(1X,A8,6F11.5)
 9256 FORMAT(A8,4(1X,F15.10),' >',/8X,2(1X,F15.10))
C                                     OCTUPOLES
      WRITE(IW,9260)
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9261)
         WRITE(  IDMP,9261)
      END IF
      DO I=1,NPTSTN
         WRITE(IW,9270)        STNPNT(1,I),XXX(I),YYY(I),ZZZ(I),XXY(I),
     *                         XXZ(I),XYY(I),YYZ(I),XZZ(I),YZZ(I),XYZ(I)
         IF(IDMP.GT.0) THEN
            WRITE(NFTEFP,9276) STNPNT(1,I),XXX(I),YYY(I),ZZZ(I),XXY(I),
     *                         XXZ(I),XYY(I),YYZ(I),XZZ(I),YZZ(I),XYZ(I)
            WRITE(  IDMP,9276) STNPNT(1,I),XXX(I),YYY(I),ZZZ(I),XXY(I),
     *                         XXZ(I),XYY(I),YYZ(I),XZZ(I),YZZ(I),XYZ(I)
         END IF
      ENDDO
      IF(IDMP.GT.0) THEN
         WRITE(NFTEFP,9206)
         WRITE(  IDMP,9206)
      END IF
 9260 FORMAT(/1X,'THIRD MOMENTS AT POINTS',
     1       /1X,'-----------------------',
     1       /1X,'NAME',11X,'XXX',8X,'YYY',8X,'ZZZ',8X,'XXY',8X,
     1                     'XXZ',8X,'XYY'/
     *        16X,'YYZ',8X,'XZZ',8X,'YZZ',8X,'XYZ')
 9261 FORMAT(' OCTUPOLES  ')
 9270 FORMAT(1X,A8,6F11.5/9X,4F11.5)
 9276 FORMAT(A8,4(1X,F15.10),' >'/8X,4(1X,F15.10),' >'/8X,2(1X,F15.10))
C
      IF(IDMP.GT.0) THEN
         CALL FLSHBF(NFTEFP)
         CALL FLSHBF(IDMP)
         CALL SEQREW(IDMP)
      END IF
C
      END IF    ! END OF A GIGANTIC "IF(MASWRK)" BLOCK
C
C     FEED THE ABIIMG COMMON BLOCKS
C
 9950 CONTINUE
      IF(IPCMIT.EQ.1 .AND. IMGABI.EQ.1) THEN
      NUMIMG = NPTSTN
C NAME, COORDINATES, MASS AND NUCLEAR CHARGE
      DO I=1,NPTSTN
         NAMIMG(I)   = INT(STNPNT(1,I))
         CRDIMG(1,I) = STNPNT(2,I)
         CRDIMG(2,I) = STNPNT(3,I)
         CRDIMG(3,I) = STNPNT(4,I)
         MASIMG(I)   = INT(ZMSSTN(I))
         ZANIMG(I)   = ZANSTN(I)
      ENDDO
C MONOPOLES
      DO I=1,NPTSTN
         CHGIMG(1,I) = CHG(I)
         CHGIMG(2,I) = ZANSTN(I)
      ENDDO
C DIPOLES
      DO I=1,NPTSTN
         DIPIMG(1,I) = X(I)
         DIPIMG(2,I) = Y(I)
         DIPIMG(3,I) = Z(I)
      ENDDO
C QUADRUPOLES
      DO I=1,NPTSTN
         QADIMG(1,I) = XX(I)
         QADIMG(2,I) = YY(I)
         QADIMG(3,I) = ZZ(I)
         QADIMG(4,I) = XY(I)
         QADIMG(5,I) = XZ(I)
         QADIMG(6,I) = YZ(I)
      ENDDO
      END IF
C
      RETURN
      END
C*MODULE PRPPOP  *DECK STNRD
      SUBROUTINE STNRD(DBUG,ZANSTN,ZMSSTN,
     *                 IDMYAT,IPAIRS,IDIN,IDOUT,IDXTR,NATM,NNMO,
     *                 NRAD,NANG,NBECKE,NBCKRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION MAKEFP
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFGPT=12000)
      PARAMETER (NMO=500, MXMID=50)
      PARAMETER (PT5=0.5D+00, ZERO=0.0D+00, ZTOL=1.0D-05)
C
      CHARACTER*8 DUMMY, WORD2,WORD3
C
      LOGICAL DBUG,GOPARR,DSKWRK,MASWRK,MFRZ,DEEPIN,MOIDON,EDCOMP,
     *        DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,POLANG,POLAPP,KMIDPT,
     *        POLDYN
C
      DIMENSION IPAIRS(2,MXATM),ZANSTN(NATM),ZMSSTN(NATM),
     *          IDMYAT(NNMO),IDIN(NATM),IDOUT(NATM),IDXTR(NATM)
C
      COMMON /EDCMP / ZIJ(NMO),ZMO(5,NMO),OCCUP(NMO),DPFREQ(50),
     *                MOIDNO(5,NMO),IJMO(2,NMO),MOIJ(NMO),NMOIJ(NMO),
     *                NMOAT(NMO),NDPFREQ,IPROT(5),NPROT,
     *                MOIDON,EDCOMP,DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,
     *                POLANG,POLAPP,KMIDPT,POLDYN
      COMMON /EFPBF2/ QBORD(NMO)
      COMMON /EFPPRT/ NOPRT(MXATM),MIDPRT(MXATM),NUMFFD(MXATM),
     *                INOPRT(2*MXATM),LSTGRP(20,MXATM+1),
     *                ISUM(20,MXFGPT),GRPSUM(20,4),NMDFFD
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IS,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /MFRPAR/ MFRZ,NUMFRZ,IFRZ(MXAO)
      COMMON /NOPOL / NOPATM(MXATM),KPOINT(MXMID),NODENS(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
C
C
      DATA  ATOM/8HATOM    /, BOND/8HBOND    /, POINT/8HPOINT   /,
     *     CMASS/8HCMASS   /, ATOMS/8HATOMS   /
      DATA MAKEFP /8HMAKEFP  /, BONDS /8HBONDS   /,BUFMO/8HQMMMBUF /
C
      CALL PREFIN
C
C    POSITION INPUT TO $STONE, AND READ THE INPUT
C
      NPTSTN=0
      NBUFFM=0
      CALL VCLR(ZANSTN,1,MXATM)
      CALL VCLR(ZMSSTN,1,MXATM)
      CALL VICLR(IDMYAT,1,NMO)
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $STONE ',IEOF)
      IF(.NOT.(IPCMIT.EQ.1 .AND. IMGABI.EQ.1)) THEN
        IF (IEOF.EQ.1 .AND. RUNTYP.NE.MAKEFP) RETURN
      END IF
      ISTN = IEOF
C
      IF(MASWRK.AND.(.NOT.(IPCMIT.EQ.1.AND.IMGABI.EQ.1)))
     *                    WRITE(IW,900)
C
C         PROVIDE DEFAULTS
C
      DBUG =.FALSE.
      BIGEXP=ZERO
C      DEFAULT VALUE FOR # OF RADIAL AND ANGULAR  GRID POINTS
      NRAD=1
      NANG=0
      NBECKE = 0
      NBCKRD=-1
C
      IEOF = 0
      IERR = 0
C
      IF(ISTN.NE.1) THEN
      CALL RDCARD('$STONE1 ',IEOF)
      IF(IEOF.NE.0) CALL ABRT
C
C        THE FIRST LINE OF $STONE HAS A FEW GLOBAL OPTIONS
C
 12   CONTINUE
      WORD2='        '
      LEN=-8
      CALL GSTRNG(WORD2,LEN)
      IF(WORD2.EQ.'        ') GO TO 16
C
C         THE "DEBUG" KEYWORD IS NOT DOCUMENTED.
C
      IF(WORD2.EQ.'DEBUG   ') THEN
         WORD3='        '
         LGSTR=-8
         CALL GSTRNG(WORD3,LGSTR)
         DBUG=.FALSE.
         IF(WORD3.EQ.'.TRUE.  ' .OR. WORD3.EQ.'.T.     ') DBUG=.TRUE.
         GO TO 12
      END IF
C
C     LVS 2005
C     UPGRADED STONE PROCEDURE
C     (USING NUMERICAL INTEGRATION OVER DIFFUSE EXPONENTS):
C
      IF(WORD2.EQ.'BIGEXP') THEN
        BIGEXP=-2
        BIGEXP = IFIND('BIGEXP  ',IERR)
        IF(IERR.NE.0) CALL ABRT
        IF(BIGEXP.LT.ZERO) THEN
           IF (MASWRK) WRITE(IW,*)' BIGEXP HAS TO BE LARGER OR EQUAL ',
     *            'TO ZERO'
           CALL ABRT
           STOP
        END IF
        GO TO 12
      END IF
C
      IF(WORD2.EQ.'NRAD') THEN
        NRAD = -2
        NRAD = IFIND('NRAD    ',IERR)
        IF(IERR.NE.0) CALL ABRT
        IF(NRAD.LE.0) THEN
           IF (MASWRK) WRITE(IW,*)' NUMBER OF RADIAL GRID POINTS ',
     *            'SHOULD BE LARGER THAN ZERO'
           CALL ABRT
           STOP
        END IF
        GO TO 12
      END IF
C
      IF(WORD2.EQ.'NANG') THEN
        NANG=-2
        NANG = IFIND('NANG    ',IERR)
        IF(IERR.NE.0) CALL ABRT
        IF(NANG.LT.0) THEN
           IF (MASWRK) WRITE(IW,*)' NUMBER OF ANGULAR GRID POINTS ',
     *            'SHOULD BE LARGER THAN ZERO'
           CALL ABRT
           STOP
        END IF
        GO TO 12
      END IF
C
      IF(WORD2.EQ.'SMOOTH') THEN
         NBECKE=-2
         NBECKE = IFIND('SMOOTH  ',IERR)
        IF(IERR.NE.0) CALL ABRT
        IF(NBECKE.LT.0) THEN
           IF (MASWRK) WRITE(IW,*)' BECKE SMOOTHING PARAMETER  ',
     *            'SHOULD BE EQUAL OR LARGER THAN ZERO'
           CALL ABRT
           STOP
        END IF
        GO TO 12
      END IF
C
C   EFFECTIVE RADIUS FOR BECKE SMOOTHING
      IF(WORD2.EQ.'SMRAD') THEN
         NBCKRD=-1
         NBCKRD = IFIND('SMRAD   ',IERR)
         IF(IERR.NE.0) CALL ABRT
        IF(NBCKRD.NE.0 .AND. NBCKRD.NE.1) THEN
           IF (MASWRK) WRITE(IW,*)' SMRAD PARAMETER  ',
     *            'SHOULD BE SET TO 0 OR 1'
           CALL ABRT
           STOP
        END IF
        GO TO 12
      END IF
C
      IF(MASWRK) WRITE (IW,*) '$STONE: UNRECOGNIZED KEYWORD=',WORD2
      CALL ABRT
  16  CONTINUE
C
C         PROVIDE DEFAULTS IF KEYWORD WAS NOT READ
C
C      DEFAULT VALUES OF PARAMETERS FOR NUMERICAL STONE
      IF (BIGEXP.GT.ZERO) THEN
         IF(NRAD.EQ.1) NRAD = 100
         IF(NANG.EQ.0) NANG=590
         IF(NBECKE.EQ.0) NBECKE = 2
C    DEFAULT VALUE FOR SMOOTHING RADIUS:
C    USE CONSTANT RADII FOR SMOOTHING
C    OTHERWISE USE BRAGG-SLATER RADII OF ATOMS
         IF(NBCKRD.EQ.-1) NBCKRD=1
      END IF
C
C     PRINTING OUT
      IF(MASWRK) THEN
         WRITE(IW,9455) BIGEXP, NRAD, NANG, NBECKE
         IF(NBCKRD.EQ.0) THEN
            WRITE(IW,*) ' CONSTANT RADII WILL BE USED FOR SMOOTHING'
         ELSE IF (NBCKRD.EQ.1) THEN
            WRITE(IW,*)' BRAGG-SLATER RADII WILL BE USED FOR SMOOTHING'
         END IF
      END IF
C
      IERR = 0
      IEOF = 0
      CALL RDCARD('$STONE2 ',IEOF)
      IF(IEOF.NE.0) CALL ABRT
C
      DUMMY='        '
      LEN=-8
      CALL GSTRNG(DUMMY,LEN)
      IF (DUMMY.EQ.'$END   ') THEN
C           DO NOT RETURN STRAIGHT AWAY SINCE WANT TO OBTAIN
C           DEFAULT STONE ATOMS&MID-BOND POSITIONS FIRST
          ISTN=1
       ELSE
C            THERE IS SOMETHING ELSE TO READ IN $STONE GROUP...
          ISTN=0
          GOTO 501
       END IF
C
C         THIS IS END OF (ISTN.NE.1) IF-STATEMENT
      END IF
C
      ISTEP = 0
   50 CONTINUE
      ISTEP = ISTEP + 1
C
C        DEFAULT ATOMS AND BONDS FOR IMGABI METHOD
C
      IF ((ISTN.EQ.1 .AND. RUNTYP.EQ.MAKEFP)
     *.OR.(ISTN.EQ.1.AND.IMGABI.EQ.1.AND.IPCMIT.EQ.1)) THEN
         IF (ISTEP.EQ.1) THEN
            TOKEN = ATOMS
            GO TO 150
         END IF
         IF (ISTEP.EQ.2) THEN
            TOKEN = BONDS
            GO TO 150
         END IF
         IF (ISTEP.EQ.3) RETURN
      END IF
C
 500  CONTINUE
      CALL RDCARD('$STONE3 ',IEOF)
      IF(IEOF.NE.0) CALL ABRT
      DUMMY='        '
      LEN=-8
      CALL GSTRNG(DUMMY,LEN)
      IF(DUMMY.EQ.'$END    ') RETURN
C
 501  CONTINUE
      READ(UNIT=DUMMY, FMT='(A8)') TOKEN
      IF(MASWRK.AND.(.NOT.(IPCMIT.EQ.1.AND.IMGABI.EQ.1)))
     *                    WRITE(IW,910)TOKEN
C
C      IF(TOKEN.EQ.DEBUG .OR. TOKEN.EQ.SWITCH .OR. TOKEN.EQ.RADGRD
C     *   .OR. TOKEN.EQ.ANGGRD .OR. TOKEN.EQ.SMOOTH.OR.TOKEN.EQ.BCKRAD)
C     *   GO TO 50
C
C 'ATOM' TOKEN FOLLOWED BY ATOM NUMBER AND OPTIONAL NAME.
C
      IF(TOKEN.EQ.ATOM) THEN
        NPTSTN=NPTSTN+1
        IERR = 0
        LATM = IFIND('LATM    ',IERR)
        IF(IERR.NE.0) CALL ABRT
        LEN=-8
        CALL GSTRNG(DUMMY,LEN)
        IF(DUMMY.EQ.'        ') THEN
          STNPNT(1,NPTSTN) = ANAM(LATM)
        ELSE
          READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
        END IF
        DO 90 LC=1,3
         STNPNT(LC+1,NPTSTN)=C(LC,LATM)
 90     CONTINUE
        ZANSTN(NPTSTN)=ZAN(LATM)
        ZMSSTN(NPTSTN)=ZMASS(LATM)
        GO TO 500
      END IF
C
C       BOND MIDPOINT TOKEN WITH DEFINING ATOMS AND OPTIONAL NAME.
C
      IF(TOKEN.EQ.BOND) THEN
        NPTSTN=NPTSTN+1
        LATM1 = IFIND('LATM    ',IERR)
        IF(IERR.NE.0) CALL ABRT
        LATM2 = IFIND('LATM    ',IERR)
        IF(IERR.NE.0) CALL ABRT
        LEN=-8
        CALL GSTRNG(DUMMY,LEN)
        IF(DUMMY.EQ.'        ') THEN
          STNPNT(1,NPTSTN) = BOND
        ELSE
          READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
        END IF
        DO 120 LC=1,3
         STNPNT(LC+1,NPTSTN)=PT5*(C(LC,LATM1)+C(LC,LATM2))
 120    CONTINUE
        GO TO 500
      END IF
C
C 'POINT' TOKEN FOLLOWED BY COORDINATES AND OPTIONAL NAME.
C
      IF(TOKEN.EQ.POINT) THEN
        NPTSTN=NPTSTN+1
        DO 130 LC=2,4
          STNPNT(LC,NPTSTN) = RFIND('STNPNT  ',IERR)
          IF(IERR.NE.0) CALL ABRT
 130    CONTINUE
        LEN=-8
        CALL GSTRNG(DUMMY,LEN)
        IF(DUMMY.EQ.'        ') THEN
          STNPNT(1,NPTSTN) = POINT
        ELSE
          READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
        END IF
        GO TO 500
      END IF
C
C       CENTER-OF-MASS TOKEN FOLLOWED BY OPTIONAL NAME.
C
      IF(TOKEN.EQ.CMASS) THEN
        NPTSTN=NPTSTN+1
        LEN=-8
        CALL GSTRNG(DUMMY,LEN)
        IF(DUMMY.EQ.'        ') THEN
          STNPNT(1,NPTSTN) = CMASS
        ELSE
          READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
        END IF
        CALL CALCOM(CMX,CMY,CMZ)
        STNPNT(2,NPTSTN)=CMX
        STNPNT(3,NPTSTN)=CMY
        STNPNT(4,NPTSTN)=CMZ
        GO TO 500
      END IF
C
C       'ATOMS' TOKEN MEANING USE ALL ATOMS AS EXPANSION POINTS.
C
  150 CONTINUE
      IF(TOKEN.EQ.ATOMS) THEN
        DO 160 LATM=1,NAT
         NPTSTN=NPTSTN+1
         STNPNT(1,NPTSTN)=ANAM(LATM)
C
C         DETERMINE WHICH ATOM-POINTS SHOULD NOT BE PRINTED OUT
C         ACCORDING TO THE INFORMATION FROM NOPRT(I)
C
         DO IAT = 1, MXATM
           IF(NOPRT(IAT).EQ.LATM) INOPRT(NPTSTN)=1
         ENDDO
C
C
C         FIND THE ATOM-POINTS WHOSE CHARGES(MONOPOLE)
C         WILL BE COLLECTED AND ADDED TO THE ATOM OR THE MIDPOINT
C         DEFINED IN LSTGRP(*,MXATM+1)
C         THE DEFAULT LSTGRP(*,MXATM+1) IS THE ATOM LSTGRP(*,1)
C
          IF(NMDFFD.GT.0) THEN
            DO NMD = 1, NMDFFD
                DO K = 1, LSTGRP(NMD,2)
                  IF(LSTGRP(NMD,K+2).EQ.LATM)
     *                 ISUM(NMD,NPTSTN)=1
                ENDDO
                IF(LSTGRP(NMD,1).EQ.LATM)
     *                LSTGRP(NMD,MXATM+1)=NPTSTN
            ENDDO
          END IF
C
C
         DO 161 LC=1,3
          STNPNT(LC+1,NPTSTN)=C(LC,LATM)
 161    CONTINUE
        ZANSTN(NPTSTN)=ZAN(LATM)
        ZMSSTN(NPTSTN)=ZMASS(LATM)
 160    CONTINUE
        GO TO 50
      END IF
C
C       'BONDS' TOKEN MEANS USE ALL BOND MIDPOINTS
C
      IF(TOKEN.EQ.BONDS) THEN
         CALL GTBOND(IPAIRS,MXATM,NBONDS,0)
         DO I = 1,NBONDS
            NPTSTN=NPTSTN+1
            LATM1 = IPAIRS(1,I)
            LATM2 = IPAIRS(2,I)
C
C         DETERMINE WHICH BOND-POINTS SHOULD NOT BE PRINTED OUT
C         ACCORDING TO THE INFORMATION FROM NOPRT(I), MIDPRT(I)
C
            DO IAT = 1, MXATM
              IF(  (NOPRT(IAT).EQ.LATM1).OR.
     *             (NOPRT(IAT).EQ.LATM2)    )
     *                INOPRT(NPTSTN) =  1
            ENDDO
            DO IAT = 1, MXATM
              IF(  (MIDPRT(IAT).EQ.LATM1).OR.
     *             (MIDPRT(IAT).EQ.LATM2)     )
     *                INOPRT(NPTSTN) = -1
            ENDDO
C
C         FIND THE BOND-POINTS WHOSE CHARGES(MONOPOLE)
C         WILL BE COLLECTED AND ADDED TO THE ATOM OR THE MIDPOINT
C         DEFINED BY LSTGRP(*,MXATM+1)
C
          IF(NMDFFD.GT.0) THEN
            DO NMD = 1, NMDFFD
                DO K = 1, LSTGRP(NMD,2)
                  IF(LSTGRP(NMD,2+K).EQ.LATM1 .OR.
     *               LSTGRP(NMD,2+K).EQ.LATM2     )
     *                 ISUM(NMD,NPTSTN)= INOPRT(NPTSTN)
                ENDDO
C
C         HERE WE DEFINE LSTGRP(*,MXATM+1) AS THE THE MIDPOINT
C         IF THE MIDPOINT WILL BE PRINTED/PUNCHED OUT
C         ALSO FIND OUT THE COORDINATES OF THE CONNECTION POINT
C         AT THE DISCARDED SIDE AND STORE THEM AT GRPSUM(*,2-4).
C
                IF(LSTGRP(NMD,1).EQ.LATM1 .OR.
     *             LSTGRP(NMD,1).EQ.LATM2      ) THEN
                   DO K = 1, LSTGRP(NMD,2)
                     IF(LSTGRP(NMD,2+K).EQ.LATM1 .OR.
     *                  LSTGRP(NMD,2+K).EQ.LATM2     ) THEN
                          IF(INOPRT(NPTSTN).LT.0) THEN
                             IF(LSTGRP(NMD,1).EQ.LATM1) THEN
                               DO LC=1,3
                                  GRPSUM(NMD,1+LC)=C(LC,LATM2)
                               ENDDO
                             ELSE
                               DO LC=1,3
                                  GRPSUM(NMD,1+LC)=C(LC,LATM1)
                               ENDDO
                             END IF
                             LSTGRP(NMD,MXATM+1) = NPTSTN
                          ELSE
                             DO LC=1,3
                                GRPSUM(NMD,1+LC)=PT5*
     *                                (C(LC,LATM1)+C(LC,LATM2))
                             ENDDO
                          END IF
                     END IF
                   ENDDO
                END IF
            ENDDO
          END IF
C
            IF (LATM1.LT.10.AND.LATM2.LT.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I1,I1)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.GE.10.AND.LATM2.LT.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I2,I1)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.LT.10.AND.LATM2.GE.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I1,I2)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.GE.10.AND.LATM2.GE.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I2,I2)') 'BO',LATM1,LATM2
            END IF
            READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
C
            DO LC=1,3
               STNPNT(LC+1,NPTSTN)=PT5*(C(LC,LATM1)+C(LC,LATM2))
            END DO
         END DO
         GO TO 50
      END IF
C
C   QM/MM STUFF: KEYWORD RELATED TO THE PRESENCE OF THE BUFFER REGION
C   OR OVERLAPPING REGION. THE KEYWORD QMMMBUF MUST BE SPECIFIED
C   WHENEVER THERE IS AN OVERLAPPING REGION AND THE USER WANTS TO
C   DESELECT EXPANSION POINTS.
C
C   'BUFMO' TOKEN (QMMMBUF KEYWORD) :
C    1. SPECIFIES NUMBER OF BUFFER MO-S EQUAL TO THE NUMBER OF FROZEN
C       ORBITALS BY DEFAULT
C    2. PERFORMS THE STONE ANALYSIS FOR THE THE QM/MM CALCULATION.
C       THE POINTS WILL INCLUDE ATOMS AND BOND MIDPOINTS OUTSIDE THE
C       BUFFER REGION ONLY
C
      IF(TOKEN.EQ.BUFMO) THEN
        NBUFFM = IFIND('QMMMBUF ',IERR)
        IF(IERR.NE.0) CALL ABRT
        IF(NBUFFM.EQ.0) NBUFFM=NUMFRZ
        IF(MASWRK) WRITE(IW,9450)NBUFFM
        IF(NBUFFM.GT.NUMFRZ) THEN
           IF (MASWRK) WRITE(IW,*)' BUFMO HAS TO BE LESS OR EQUAL ',
     *            'TO THE NUMBER OF FROZEN MO-S IN $MOFRZ'
           CALL ABRT
           STOP
        END IF
C
C        SET MOIDON TO USE LOCALIZED MOS ARE USED INSTEAD OF CANONICAL
C
      MOIDON = .TRUE.
      IF(.NOT.MOIDON) CALL MOIDSM(15)
      NUMATB=0
C
C     COLLECT THE CENTERS WHICH ARE IN THE BUFFER ZONE INTO ARRAY IDMYAT
C
      CALL VICLR(IDMYAT,1,NMO)
      IJ=0
      DO 505 I=1,NAT
      IINB=0
        DO 510 J=1,NBUFFM
          DO 520 K=1,NMOAT(J)
            IF(I.EQ.MOIDNO(K,J)) IINB=1
  520     CONTINUE
  510   CONTINUE
      IF(IINB.EQ.1) THEN
         IJ=IJ+1
         IDMYAT(IJ)=I
      END IF
  505 CONTINUE
      NUMATB=IJ
      IF(MASWRK) WRITE(IW,9500) NUMATB,(IDMYAT(I),I=1,NUMATB)
C
C      COLLECT THE CENTERS WHICH ARE IN THE OVERLAPPING REGION (IDOUT)
C      COLLECT THE CENTERS THAT WILL STAY (IDIN)
C
      CALL VICLR(IDOUT,1,MXATM)
      CALL VICLR(IDIN,1,MXATM)
      IOUT=0
      JIN=0
      DO 600 JAT=1,NAT
         MIN=0
         DO 610 JEXP=1,MXATM
            IF (JAT.EQ.NOPATM(JEXP)) MIN=1
  610    CONTINUE
C
C COLLECT CENTERS THAT WILL BE EXCLUDED
C
         IF(MIN.EQ.1) THEN
            IOUT=IOUT+1
            IDOUT(IOUT)=JAT
         END IF
C
C COLLECT CENTERS THAT WILL STAY
C
         IF(MIN.EQ.0) THEN
            JIN=JIN+1
            IDIN(JIN)=JAT
         END IF
  600 CONTINUE
      NUMOUT=IOUT
      IF(MASWRK) WRITE(IW,9520) NUMOUT,(IDOUT(IOUT),IOUT=1,NUMOUT)
C
C      COLLECT THE CENTERS EXCLUDED FROM THE DENSITY
C
      CALL VICLR(IDXTR,1,MXATM)
      IXTR=0
      DO 640 JX=1,NAT
         MXTR=0
         DO 650 JXTR=1,MXATM
           IF (JX.EQ.NODENS(JXTR)) MXTR=1
  650    CONTINUE
C
         IF(MXTR.EQ.1) THEN
            IXTR=IXTR+1
            IDXTR(IXTR)=JX
         END IF
  640 CONTINUE
      NUMXTR=IXTR
      IF(MASWRK) WRITE(IW,9530) NUMXTR,(IDXTR(IXTR),IXTR=1,NUMXTR)
C
C   DETERMINE THE FORMAL CHARGE OF THE CAPPING ATOMS OUTSIDE THE BUFFER
C   REGION (IN OTHER WORDS, CHARGE IN THE BUFFER REGION IS ZEROED OUT)
C
      CALL VCLR(QBORD,1,NMO)
      DO 730 I=1,NUMATB
         IA=IDMYAT(I)
         PILEUP=ZERO
         DO 710 K=1,NBUFFM
            DO 720 J=1,NMOAT(K)
               IF(MOIDNO(J,K).EQ.IA) THEN
                  PILEUP=PILEUP+ZMO(J,K)
               END IF
  720       CONTINUE
  710    CONTINUE
C
C    DOUBLE CHECK CHARGES ON CAPPING ATOMS.
C
         PILE02=ZERO
         DO 711 K=1,NMO
            DO 721 J=1,NMOAT(K)
               IF(MOIDNO(J,K).EQ.IA) THEN
                  PILE02=PILE02+ZMO(J,K)
               END IF
  721       CONTINUE
  711    CONTINUE
C
         IF((ABS(PILE02-ZAN(IA)).GT.ZTOL) .AND.
     *      (PILE02-PILEUP).LT.1.5D+00)
     *          PILEUP=PILEUP-(PILE02-ZAN(IA))
         IF ( ABS(PILEUP-ZAN(IA) ) .GT. ZTOL ) QBORD(IA)=ZAN(IA)-PILEUP
         IF(MASWRK.AND.QBORD(IA).GT.ZTOL) WRITE(IW,9550)IA,QBORD(IA)
  730 CONTINUE
C
C    DETERMINE WHICH ATOMS WILL BECOME STONE EXPANSION POINTS -
C    EXCLUDING DUMMIES.
C
C COLLECT MIDPOINT MOS AND USE THEM TO DETERMINE
C BOUNDARY ATOM'S NUCLEAR CHARGE
C PABLO GOT RID OF THIS (UNNECESSARY)
C
      CALL GTBOND(IPAIRS,MXATM,NBONDS,0)
C
C    THIS REQUIRES TINKERING WITH THE NUCLEAR CHARGES
C
      DO 165 LATM=1,NAT
         DO 180 I=1,NUMATB
            IF (IDMYAT(I) .EQ. LATM .AND. QBORD(I).LT.ZTOL) THEN
               ZAN(LATM)=ZERO
               GO TO 165
            END IF
  180    CONTINUE
C
C PABLO MODIFY COMMENT
C      EXCLUDE "NOPATM" OVERLAPPING ATOMS FROM EXPANSION POINTS
C
         DO 185 J=1,NAT
            IF(IDOUT(J).EQ.LATM) THEN
               ZAN(LATM) = ZERO
               GO TO 165
            END IF
C PABLO
C ASSIGN CHARGES=0 FOR ATOMS THAT BELONG TO THE NODENS ARRAY
C BUT ARE NOT PART OF THE NOPATM ARRAY
C IN OTHER WORDS ZERO OUT THE CHARGES OF THE EXTRA EXPANSION POINTS
            IF (IDXTR(J).EQ.LATM) THEN
               ZAN(LATM) = ZERO
            END IF
C
  185    CONTINUE
C PABLO
C      ASSIGN CORRECT CHARGE TO BOUNDARING ATOMS OF OVERLAP
C         DO 190 JPT=1,NAT
C            IF (KPOINT(JPT).EQ.LATM) THEN
C               IF (KMIDPT)      ZAN(LATM)=ONE
C               IF (.NOT.KMIDPT) ZAN(LATM)=ZAN(LATM) - ONE
C            END IF
C  190    CONTINUE
C
         NPTSTN=NPTSTN+1
         STNPNT(1,NPTSTN)=ANAM(LATM)
         ZMSSTN(NPTSTN)=ZMASS(LATM)
C
C         DETERMINE WHICH ATOM-POINTS SHOULD NOT BE PRINTED OUT
C         ACCORDING TO THE INFORMATION FROM NOPRT(I)
C
         DO IAT = 1, MXATM
           IF(NOPRT(IAT).EQ.LATM) INOPRT(NPTSTN)=1
         ENDDO
C
C
C         FIND THE ATOM-POINTS WHOSE CHARGES(MONOPOLE)
C         WILL BE COLLECTED AND ADDED TO THE ATOM OR THE MIDPOINT
C         DEFINED IN LSTGRP(*,MXATM+1)
C         THE DEFAULT LSTGRP(*,MXATM+1) IS THE ATOM LSTGRP(*,1)
C
          IF(NMDFFD.GT.0) THEN
            DO NMD = 1, NMDFFD
                DO K = 1, LSTGRP(NMD,2)
                  IF(LSTGRP(NMD,K+2).EQ.LATM)
     *                 ISUM(NMD,NPTSTN)=1
                ENDDO
                IF(LSTGRP(NMD,1).EQ.LATM)
     *                LSTGRP(NMD,MXATM+1)=NPTSTN
            ENDDO
          END IF
C
         DO 170 LC=1,3
            STNPNT(LC+1,NPTSTN)=C(LC,LATM)
  170    CONTINUE
         IF( QBORD(LATM).GT.ZTOL) THEN
            ZAN(NPTSTN)=QBORD(LATM)
            ZANSTN(NPTSTN)=QBORD(LATM)
         ELSE
            ZANSTN(NPTSTN)=ZAN(LATM)
         END IF
  165 CONTINUE
C
C     DETERMINE WHICH BOND MIDPOINTS WILL BECOME STONE EXPANSION POINTS
C
        NUMBON=0
          DO 300 I = 1,NBONDS
             LATM1 = IPAIRS(1,I)
             LATM2 = IPAIRS(2,I)
             DO 310 J=1,NUMATB
             INBUF=IDMYAT(J)
             DEEPIN= QBORD(J) .LT. ZTOL
               IF((LATM1.EQ.INBUF .OR. LATM2.EQ.INBUF).AND. DEEPIN )
     *                     GO TO 300
 310         CONTINUE
C
C               EXCLUDE MIDPOINTS FROM FRAGMENT WHEN ATOMS ARE
C               ELEMENTS OF THE IDOUT ARRAY
C
             DO 320 K=1,NAT
                INOVL=IDOUT(K)
                IF(LATM1.EQ.INOVL .OR. LATM2.EQ.INOVL) GO TO 300
  320        CONTINUE
C
             NUMBON=NUMBON+1
             NPTSTN=NPTSTN+1
C
C         DETERMINE WHICH BOND-POINTS SHOULD NOT BE PRINTED OUT
C         ACCORDING TO THE INFORMATION FROM NOPRT(I), MIDPRT(I)
C
            DO IAT = 1, MXATM
              IF(  (NOPRT(IAT).EQ.LATM1).OR.
     *             (NOPRT(IAT).EQ.LATM2)    )
     *                INOPRT(NPTSTN) =  1
            ENDDO
            DO IAT = 1, MXATM
              IF(  (MIDPRT(IAT).EQ.LATM1).OR.
     *             (MIDPRT(IAT).EQ.LATM2)     )
     *                INOPRT(NPTSTN) = -1
            ENDDO
C
C         FIND THE BOND-POINTS WHOSE CHARGES(MONOPOLE)
C         WILL BE COLLECTED AND ADDED TO THE ATOM OR THE MIDPOINT
C         DEFINED BY LSTGRP(*,MXATM+1)
C
          IF(NMDFFD.GT.0) THEN
            DO NMD = 1, NMDFFD
                DO K = 1, LSTGRP(NMD,2)
                  IF(LSTGRP(NMD,2+K).EQ.LATM1 .OR.
     *               LSTGRP(NMD,2+K).EQ.LATM2     )
     *                 ISUM(NMD,NPTSTN)= INOPRT(NPTSTN)
                ENDDO
C
C         HERE WE DEFINE LSTGRP(*,MXATM+1) AS THE THE MIDPOINT
C         IF THE MIDPOINT WILL BE PRINTED/PUNCHED OUT
C         ALSO FIND OUT THE COORDINATES OF THE CONNECTION POINT
C         AT THE DISCARDED SIDE AND STORE THEM AT GRPSUM(*,2-4).
C
                IF(LSTGRP(NMD,1).EQ.LATM1 .OR.
     *             LSTGRP(NMD,1).EQ.LATM2      ) THEN
                   DO K = 1, LSTGRP(NMD,2)
                     IF(LSTGRP(NMD,2+K).EQ.LATM1 .OR.
     *                  LSTGRP(NMD,2+K).EQ.LATM2     ) THEN
                          IF(INOPRT(NPTSTN).LT.0) THEN
                             IF(LSTGRP(NMD,1).EQ.LATM1) THEN
                               DO LC=1,3
                                  GRPSUM(NMD,1+LC)=C(LC,LATM2)
                               ENDDO
                             ELSE
                               DO LC=1,3
                                  GRPSUM(NMD,1+LC)=C(LC,LATM1)
                               ENDDO
                             END IF
                             LSTGRP(NMD,MXATM+1) = NPTSTN
                          ELSE
                             DO LC=1,3
                                GRPSUM(NMD,1+LC)=PT5*
     *                                (C(LC,LATM1)+C(LC,LATM2))
                             ENDDO
                          END IF
                     END IF
                   ENDDO
                END IF
            ENDDO
          END IF
C
             IF (LATM1.LT.10.AND.LATM2.LT.10) THEN
                WRITE(UNIT=DUMMY,FMT='(A2,I1,I1)') 'BO',LATM1,LATM2
             ELSE IF (LATM1.GE.10.AND.LATM2.LT.10) THEN
                WRITE(UNIT=DUMMY,FMT='(A2,I2,I1)') 'BO',LATM1,LATM2
             ELSE IF (LATM1.LT.10.AND.LATM2.GE.10) THEN
                WRITE(UNIT=DUMMY,FMT='(A2,I1,I2)') 'BO',LATM1,LATM2
             ELSE IF (LATM1.GE.10.AND.LATM2.GE.10) THEN
                WRITE(UNIT=DUMMY,FMT='(A2,I2,I2)') 'BO',LATM1,LATM2
             END IF
             READ(UNIT=DUMMY, FMT='(A8)') STNPNT(1,NPTSTN)
             DO LC=1,3
                STNPNT(LC+1,NPTSTN)=PT5*(C(LC,LATM1)+C(LC,LATM2))
             END DO
 300      CONTINUE
          NBFBN = NBONDS-NUMBON
          IF(MASWRK) WRITE(IW,9600) NBONDS,NBFBN
        GO TO 500
      END IF
C
C        UNRECOGNIZED KEYWORD FOUND
C
      IF (MASWRK) WRITE(IW,920) TOKEN
      CALL ABRT
      STOP
C
C
  900 FORMAT(/1X,'PROCESSING $STONE INPUT...')
  910 FORMAT(1X,'STONE INPUT KEYWORD=',A8)
  920 FORMAT(1X,'ILLEGAL STONE ANALYSIS EXPANSION POINT TYPE=',A8)
 9450 FORMAT(/1X,'THE NUMBER OF BUFFER MO-S IS ',I3)
 9455 FORMAT(/1X,' BIGEXP ',F6.2/
     *      1X,' RADIAL GRID ',I4,'    ANGULAR GRID ',I4/
     *      1X,' BECKE SMOOTHING PARAMETER ',I2)
 9500 FORMAT(1X,'THE NUMBER OF ATOMS IN THE QM/MM BUFFER REGION IS ',
     *       I3,:/3X,'THEIR NUMBERS ARE:',(10I3))
 9520 FORMAT(1X,'THE NUMBER OF ATOMS IN THE OVERLAP REGION IS',
     *       I3,:/3X,'THEIR NUMBERS ARE:',(10I3))
 9530 FORMAT(1X,'NUMXTR=',I6/(1X,15I5/))
 9550 FORMAT(1X,'ATOM ',I3,' IS A CAPPING (FRONTIER) ATOM'
     *       /3X,' WITH A FORMAL NUCLEAR CHARGE IN MM REGION ',F4.2)
 9600 FORMAT(1X,'PROGRAM ORIGINALLY FOUND ',I4,' TWO-CENTER BONDS,'/
     *   3X,'OUT OF WHICH ',I4,' BELONG TO THE BUFFER REGION'/
     *   3X,'AND/OR THE OVERLAPPING REGION OF THE EFP FRAGMENT'/
     *   3X,'AND WILL BE EXCLUDED FROM STONE EXPANSION CENTERS LIST.'/)
      END
C*MODULE PRPPOP  *DECK STNXYZ
      SUBROUTINE STNXYZ
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /XYZMOM/ TX0,TY0,TZ0,TX1,TY1,TZ1,TX2,TY2,TZ2,TX3,TY3,TZ3,
     *                TX4,TY4,TZ4,T,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,CX,CY,CZ
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      DIMENSION H(36),W(36),MIN(8),MAX(8)
      DATA MIN /1,2,4, 7,11,16,22,29/
      DATA MAX /1,3,6,10,15,21,28,36/
      EQUIVALENCE (H(1),H1),(W(1),W1)
C
      DATA ZERO /0.0D+00/
C
C     ----- MULTIPOLE INTEGRALS (FOR STONE DMA) -----
C     THIS ROUTINE HANDLES S,P,D,F,G,H,I BASIS FUNCTIONS, AND
C     THIS ROUTINE HANDLES MULTIPOLES UP TO OCTUPOLE.
C
      TX0=ZERO
      TY0=ZERO
      TZ0=ZERO
      TX1=ZERO
      TY1=ZERO
      TZ1=ZERO
      TX2=ZERO
      TY2=ZERO
      TZ2=ZERO
      TX3=ZERO
      TY3=ZERO
      TZ3=ZERO
      TX4=ZERO
      TY4=ZERO
      TZ4=ZERO
C
C     GAUSS-HERMITE QUADRATURE USING MINIMUM POINT FORMULA.
C     THE 3 IN THE NEXT LINE IS MOMENTS THOUGH OCTUPOLE.
C
      NPTSTN=(NI+NJ-2+3)/2+1
      IMIN=MIN(NPTSTN)
      IMAX=MAX(NPTSTN)
C
      DO 200 I=IMIN,IMAX
         DUM=W(I)
         PX=DUM
         PY=DUM
         PZ=DUM
         DUM=H(I)/T
         AX=DUM-XI
         AY=DUM-YI
         AZ=DUM-ZI
         BX=DUM-XJ
         BY=DUM-YJ
         BZ=DUM-ZJ
         CTX=DUM-CX
         CTY=DUM-CY
         CTZ=DUM-CZ
C
         GO TO (1,2,3,4,5,6,7), NI
C
    7    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    6    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    5    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    4    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    3    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    2    PX=PX*AX
         PY=PY*AY
         PZ=PZ*AZ
    1    CONTINUE
C
         GO TO (11,12,13,14,15,16,17), NJ
C
   17    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   16    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   15    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   14    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   13    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   12    PX=PX*BX
         PY=PY*BY
         PZ=PZ*BZ
   11    CONTINUE
C
         TX0=TX0+PX
         TY0=TY0+PY
         TZ0=TZ0+PZ
         TX1=TX1+PX*CTX
         TY1=TY1+PY*CTY
         TZ1=TZ1+PZ*CTZ
         TX2=TX2+PX*CTX*CTX
         TY2=TY2+PY*CTY*CTY
         TZ2=TZ2+PZ*CTZ*CTZ
         TX3=TX3+PX*CTX*CTX*CTX
         TY3=TY3+PY*CTY*CTY*CTY
         TZ3=TZ3+PZ*CTZ*CTZ*CTZ
  200 CONTINUE
      RETURN
      END
C*MODULE PRPPOP  *DECK STONE
      SUBROUTINE STONE(NFTEFP,FRGNME)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, NMO=500, MXMID=50,
     *           MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MXDIM=5810)
      CHARACTER*8 FRGNME
C
      LOGICAL DBUG,GOPARR,DSKWRK,MASWRK,MOIDON,EDCOMP,DIPDCM,
     *        DEPRNT,QADDCM,ZDO,POLDCM,POLANG,POLAPP,KMIDPT,
     *        POLDYN
C
      COMMON /EDCMP / ZIJ(NMO),ZMO(5,NMO),OCCUP(NMO),DPFREQ(50),
     *                MOIDNO(5,NMO),IJMO(2,NMO),MOIJ(NMO),NMOIJ(NMO),
     *                NMOAT(NMO),NDPFREQ,IPROT(5),NPROT,
     *                MOIDON,EDCOMP,DIPDCM,DEPRNT,QADDCM,ZDO,POLDCM,
     *                POLANG,POLAPP,KMIDPT,POLDYN
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMOOPT/ ESPSCA(9),RESPAP(2),RESPPC(2),RESDIM,RESTRI(4),
     *                RCORSD,RESPCT,CONVFG,CNVDMP,COROFF,IXESP,MXITFG,
     *                NGUESS,NBSSE,MODORB,MODPAR,IRSTSTP,IRSTLAY,NPRFMO,
     *                NFMOPAL,MODPRP,MAXL1C,IPIEDA,MODGRD,MODESP,IVMUL
      COMMON /FMOPNT/ LICHFG,LMULFG,LIDMREC,LFRGNAM,LLAYFRG,LINDAT,
     *                LNCBS,LFMOZAN,LFMOC,LFMOMAS,LIZBAS,LIAGLOB,LIBDGH,
     *                LIABDFG,LJABDFG,LNCAO,LIDXCAO,LIAPRJO,LJAPRJO,
     *                LCOREAO,LOCCCOR,LSHIFTB,LIODFMO,LFMODA,LFMODB,
     *                LFMOESPA,LFMOESPB,LLOCFMO,LSCFFRG,LRIJ,LPOPMUL,
     *                LPOPMAT,LIALOC,LINDBD,LIATFRG,LINDFRG,LINDGFRG,
     *                LNATFRG,LNAT0FRG,LIANFRG,LZANFRG,LCFRG,LLIBISH,
     *                LLIBNSH,LLIBNG,LINDATG,LFMOBUF(3),LFMODE,LNUMFRG,
     *                LLOCTAT,LIAOGLOB,LLOADM,LFMOGE,LDGRID,LIODCFMO,
     *                LJOB2GRP,LFMOPG,LEMOCDR,LUNTXYZ,LUNTROT,LSTONEP,
     *                LMAPSU,LFRGMUL
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NOPOL / NOPATM(MXATM),KPOINT(MXMID),NODENS(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
C
C         THE FOLLOWING FOUR SHOULD BE ALLOCATED DYNAMICALLY SOMEDAY
C
      DIMENSION IDMYMO(NMO),OCC(NMO),IDZERO(NMO),IDMYMP(MXMID)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
C     SETUP CORE ALLOCATION HERE.
C     CHECK NPRINT AND DBUG OPERATION.
C
      DBUG=.FALSE. .AND. MASWRK
C
C                IF WE ARE USING IMGABI METHOD, SET NFRG=0
C
      NFRGBK = NFRG
      IF(IPCMIT.EQ.1 .AND. IMGABI.EQ.1) NFRG = 0
C
C     READ EXPANSION CENTERS AND GET CORE.
C         THE LENGTH OF -STNPNT- IS UNKNOWN BEFORE CALLING -STNRD-
C
      CALL VALFM(LOADFM)
      LIPAIRS = LOADFM  + 1
      LIDIN   = LIPAIRS + MXATM*2
      LIDOUT  = LIDIN   + MXATM
      LIDXTR  = LIDOUT  + MXATM
      LZANSTN = LIDXTR  + MXATM
      LZMSSTN = LZANSTN + MXATM
      LIDMYAT = LZMSSTN + MXATM
      LAST1   = LIDMYAT + NMO
C
      CALL STNRD(DBUG,X(LZANSTN),X(LZMSSTN),
     *           X(LIDMYAT),X(LIPAIRS),X(LIDIN),X(LIDOUT),X(LIDXTR),
     *           MXATM,NMO,NRAD,NANG,NBECKE,NBCKRD)
      IF(NPTSTN.EQ.0) THEN
         NFRG = NFRGBK
         RETURN
      END IF
C
C         CAREFUL HERE, THE 21*NPTSTN ARRAY IS ACTUALLY MANY ARRAYS,
C         SEE THE CALL TO STNMOM TO DOUBLE CHECK EVERYTHING.
C         THE EXTRA STORAGE FOR -NDM- IS A BIT OBSCURE TO ME.
C
      NUM2=(NUM*NUM+NUM)/2
      NUM3=NUM*NUM
      NPT41   = LAST1
      NDM     = NPT41   + NPTSTN*21
C  MEMORY ALLOCATION FOR NUMERICAL INTEGRATION
      LGRID   = NDM     + 3*NUM2+NUM3+NBUFFM
      LTXYZ   = LGRID   + 4*(NRAD-1)*NANG*NPTSTN
      LLXYZ   = LTXYZ   + NANG*3
      LW      = LLXYZ   + MXDIM*3
      LWR     = LW      + NANG
      LR      = LWR     + NRAD
      LRHO    = LR      + NRAD
      LSTART  = LRHO    + (NRAD-1)*NANG*NPTSTN
      LAST    = LSTART  + NPTSTN+1
C
      IF(BIGEXP.GT.ZERO) THEN
         LGRGRID = LAST
         LRRGRID = LGRGRID + NPTSTN
         LAAGRID = LRRGRID + NPTSTN*NPTSTN
         LSGRID  = LAAGRID + NPTSTN*NPTSTN
         LPPGRID = LSGRID  + NPTSTN*NPTSTN
         LRADIUS = LPPGRID + NPTSTN
         LAST    = LRADIUS + NPTSTN
      ELSE
         LGRGRID=LAST
         LRRGRID=LAST
         LAAGRID=LAST
         LSGRID =LAST
         LPPGRID=LAST
         LRADIUS=LAST
      END IF
C
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
C
C     COLLECT MIDPOINT (BOUNDARY ATOM CORE MO)
C
      CALL VICLR(IDMYMP,1,MXMID)
      INDMP=0
      DO 12 IMO=1,NA
         JMAX=NMOAT(IMO)
         NKMPT=0
         DO 14 IAT=1,MXMID
            IF((JMAX.EQ.1).AND.KMIDPT) THEN
               IF (KPOINT(IAT).EQ.MOIDNO(1,IMO)) NKMPT=1
            END IF
   14    CONTINUE
         IF(NKMPT.EQ.1) THEN
            INDMP=INDMP+1
            IDMYMP(INDMP)=IMO
         END IF
  12  CONTINUE
      NUMMPT=INDMP
      IF(MASWRK.AND.(.NOT.(IPCMIT.EQ.1.AND.IMGABI.EQ.1)))
     *      WRITE(IW,9530) NUMMPT,(IDMYMP(INDMP),INDMP=1,NUMMPT)
C
C              COLLECT MOS TO BE DELETED
C
         CALL VICLR(IDMYMO,1,NMO)
         INDMO=0
         DO 20 JMO=1,NA
            NDELMO=0
            DO 30 JAT=1,MXATM
               DO 40 KROW=1,NMOAT(JMO)
C PABLO CHANGE NOPATM TO DENSITY-RELATED ARRAY NODENS
                  IF (NODENS(JAT).EQ.MOIDNO(KROW,JMO)) NDELMO=1
   40          CONTINUE
   30       CONTINUE
C
C IDENTIFY MO BELONGING TO BOUNDARY-CORE,ADD TO ARRAY IF KMIDPT=.T.
C
            DO 50 IMP=1,NUMMPT
               IF (JMO.EQ.IDMYMP(IMP)) NDELMO=1
   50       CONTINUE
            IF(NDELMO.EQ.1) THEN
            INDMO=INDMO+1
            IDMYMO(INDMO)=JMO
            END IF
   20    CONTINUE
         NUMMO=INDMO
C
C NUMMO IS TOTAL NUMBER OF ORBITALS TO BE SUBTRACTED
C
C LOOP TO ASSIGN OCCUPANCY NUMBERS TO ORBITALS.  IF ORBITAL IS
C PART OF BUFFER OR OVERLAPPING REGION THE OCCUPANCY NUMBER IS ZERO
C
         CALL VICLR (OCC,1,NMO)
         CALL VICLR (IDZERO,1,NMO)
         INDZER=0
         DO 60 NORB=1,NA
            OCC(NORB)=TWO
            IF (NORB.LE.NBUFFM) OCC(NORB)=ZERO
            DO KOVLP=1,NUMMO
               IF (NORB.EQ.IDMYMO(KOVLP)) OCC(NORB)=ZERO
            ENDDO
            IF (OCC(NORB).LT.TWO) THEN
               INDZER=INDZER + 1
               IDZERO(INDZER) = NORB
            END IF
   60    CONTINUE
         NUMZER=INDZER
      IF(MASWRK.AND.(.NOT.(IPCMIT.EQ.1.AND.IMGABI.EQ.1)))
     *      WRITE(IW,9540) NUMZER,(IDZERO(INDZER),INDZER=1,NUMZER)
C
C     READ DENSITY MATRIX AND DOUBLE OFF-DIAGONAL ELEMENTS.
C
      IF(NBUFFM.EQ.0  .AND.  NUMMO.EQ.0) THEN
         CALL DENDD1(X(NDM),X(NDM+NUM2),NUM2)
C
C SUBTRACT DENSITY OF ORBITALS ASSOCIATED TO BUFFER OR OVERLAP
C
      ELSE
         CALL DAREAD(IDAF,IODA,X(NDM+3*NUM2),NUM3,71,0)
         CALL DMTX(X(NDM),X(NDM+3*NUM2),OCC,NA,NUM,NUM)
      END IF
C
      N=NDM-1
      DO I=1,NUM
         DO J=1,I
            N=N+1
            X(N)=TWO*X(N)
         ENDDO
         X(N)=X(N)/TWO
      ENDDO
C
C     ANALYSIS AND PRINTING.
C
      CALL STNMOM(X(NDM),NUM2,
     *        X(NPT41          ),X(NPT41+   NPTSTN),X(NPT41+ 2*NPTSTN),
     *        X(NPT41+ 3*NPTSTN),X(NPT41+ 4*NPTSTN),X(NPT41+ 5*NPTSTN),
     *        X(NPT41+ 6*NPTSTN),X(NPT41+ 7*NPTSTN),X(NPT41+ 8*NPTSTN),
     *        X(NPT41+ 9*NPTSTN),X(NPT41+10*NPTSTN),X(NPT41+11*NPTSTN),
     *        X(NPT41+12*NPTSTN),X(NPT41+13*NPTSTN),X(NPT41+14*NPTSTN),
     *        X(NPT41+15*NPTSTN),X(NPT41+16*NPTSTN),X(NPT41+17*NPTSTN),
     *        X(NPT41+18*NPTSTN),X(NPT41+19*NPTSTN),
     *        X(LZANSTN),X(LZMSSTN),
     *        X(NPT41+20*NPTSTN),X(LGRID),X(LTXYZ), X(LLXYZ), X(LW),
     *        X(LWR),X(LR), X(LRHO), X(LSTART),
     *        NRAD,NANG,MXDIM,NBECKE,NBCKRD,
     *        X(LGRGRID),X(LRRGRID),X(LAAGRID),X(LSGRID),X(LPPGRID),
     *        X(LRADIUS),NFTEFP,FRGNME)
C
C     FOR FMO WE ASSUME THAT ATOMIC CENTRES ARE ALWAYS USED.
C     FMO WANTS "POPULATIVE" CHARGES (CHARGES MEASURES IN ELECTRONS).
      IF(NFG.NE.0) THEN
        NFRGMUL=0
        IF(IVMUL.GT.1) NFRGMUL=NFRGMUL+3
        IF(IVMUL.GT.2) NFRGMUL=NFRGMUL+6
        IF(IVMUL.GT.3) NFRGMUL=NFRGMUL+10
C       HERE WE RELY ON NPT41 HAVING A CONTINUOUS FLOW OF MULTIPOLES.
        CALL DCOPY(NAT*(NFRGMUL+1),X(NPT41),1,X(LSTONEP),1)
        CALL DSCAL(NAT,-1.0D+00,X(LSTONEP),1)
C       CHANGE FROM CHARGES TO ATOMIC POPULATIONS.
      ENDIF
      CALL RETFM(NEED)
      IF(MASWRK.AND.(.NOT.(IPCMIT.EQ.1.AND.IMGABI.EQ.1))) CALL TIMIT(1)
      NFRG = NFRGBK
      RETURN
C
 9530 FORMAT(1X,'THE NUMBER OF "BOUNDARY" CORE MOS IS',
     *       I3,:/3X,'THEIR NUMBERS ARE:',(10I3))
 9540 FORMAT(1X,'THE NUMBER OF MOS WITH ZERO OCCUPATION NUMBER IS',
     *       I3,:/3X,'THEIR NUMBERS ARE:',(10I3))
      END
C*MODULE PRPPOP  *DECK DOGRID
      SUBROUTINE DOGRID(ZANSTN, GRID,XYZ, XYZLEB, W,WR,R, NSTART,
     *                  N_R,N_A,MXDIM,K_MU,IFBRRAD,
     *                  GR,RR,AA,S,PP,RADIUS)
C
C  RADIUS(K) IS ITS NOTIONAL RADIUS, ZANSTN(K) IS ITS ATOMIC
C  NUMBER (ZERO FOR AN OFF-ATOM SITE). THE RADIUS IS IGNORED UNLESS
C  THE SLATER FLAG IS TRUE, IN WHICH CASE IT CONTROLS THE BECKE
C  PARTITIONING OF SPACE BETWEEN SITES. THE ATOMIC NUMBER CONTROLS
C  THE SCALE OF THE RADIAL GRID.
C  ADOPTED FROM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C  NG IS THE TOTAL NUMBER OF GRID POINTS
C  N_R IS THE NUMBER OF RADIAL POINTS, N_A THE NUMBER OF ANGULAR POINTS.
C  K_MU IS THE VALUE OF THE BECKE SMOOTHING PARAMETER. M_R IS THE
C  PARAMETER FOR THE RADIAL QUADRATURE.
      PARAMETER (MXATM=2000,ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00,
     *      FOUR=4.0D+00,
     *      PI=3.14159265358979D+00,
     *      M_R=2)
C
      CHARACTER*8 DUMMY
      CHARACTER*6 CHTMP
C
      DIMENSION GRID(4,NPTSTN*(N_R-1)*N_A),GR(NPTSTN),RR(NPTSTN,NPTSTN),
     *          AA(NPTSTN,NPTSTN), S(NPTSTN,NPTSTN),
     *          PP(NPTSTN),NSTART(NPTSTN+1),
     *          ZANSTN(NPTSTN), RADIUS(NPTSTN),
     *          XYZ(N_A,3), XYZLEB(3,MXDIM), W(N_A),WR(N_R),R(N_R),
     *          SLTRAD(0:54)
C
C  BRAGG-SLATER RADII FROM SLATER, JCP (1964) 41, 3199. INERT GASES
C  ADDED WITH SAME RADIUS AS PRECEDING HALOGEN. HYDROGEN RADIUS IS
C  TWICE THE SLATER VALUE. THESE VALUES ARE IN ANGSTROM.
C  SLTRAD(0) : BOND MID-POINT
C      DATA (SLTRAD(I),I=0,54)/ 0.65D+00, 0.50D+00, 0.50D+00,
      DATA (SLTRAD(I),I=0,54)/ 0.5D+00, 0.35D+00, 0.50D+00,
     *  1.45D+00, 1.05D+00, 0.85D+00, 0.70D+00, 0.65D+00, 0.60D+00,
     *  0.50D+00, 0.50D+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, 2.35D+00, 2.00D+00, 1.80D+00, 1.55D+00,
     *  1.45D+00, 1.45D+00,
     *  1.35D+00, 1.30D+00, 1.35D+00, 1.40D+00, 1.60D+00, 1.55D+00,
     *  1.55D+00, 1.45D+00,
     *  1.45D+00, 1.40D+00, 1.40D+00, 1.40D+00/
C
C  THE RADIUS VALUES ARE SCALED UP BY THIS FACTOR FOR THE RADIAL
C INTEGRATION.
      PARAMETER (RSCALE=2.0D+00)
C
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
C
      NG = NPTSTN*(N_R-1)*N_A
C  DETERMINE RADII
      DO I=1,NPTSTN
      IF (IFBRRAD.EQ.1) THEN
        IF (ZANSTN(I).NE. ZERO) THEN
           IZNUC = INT(ZANSTN(I))
           RADIUS(I)=SLTRAD(IZNUC)
        ELSE
C            THIS IS FOR BOND MID-POINTS
           DUMMY='        '
           WRITE(UNIT=DUMMY,FMT='(A8)') STNPNT(1,I)
           CHTMP = DUMMY(3:8)
C           WRITE(IW,121) CHTMP
C 121       FORMAT('CHTMP=',A6)
           READ(UNIT=CHTMP, FMT='(I6)') IBOND
C           WRITE(IW,*)'IBOND=',IBOND
C            NOW IBOND CONTAINS THE ATOM UMBERS WHICH PRODUCE THE BOND
C             THIS IS A REDUNDANT CASE, JUST TO BE SURE
           IF(IBOND.EQ.0) THEN
              LATM1 = 0
              LATM2 =0
              IZNUC = INT(ZANSTN(I))
              RADIUS(I)=SLTRAD(IZNUC)
           ELSE
             IF (IBOND.LT.1000) THEN
                LATM1=IBOND/10
                LATM2=IBOND-(LATM1*10)
             ELSE IF(IBOND.GE.1000) THEN
                LATM1=IBOND/100
                LATM2=IBOND-(LATM1*100)
             END IF
            IZNUC1 = INT(ZANSTN(LATM1))
            IZNUC2 = INT(ZANSTN(LATM2))
            RADIUS(I)=(SLTRAD(IZNUC1)+SLTRAD(IZNUC2))/TWO
           END IF
           WRITE(IW,122) LATM1,LATM2,RADIUS(I)
 122       FORMAT(1X,'RADIUS FOR BOND MID-POINT BETWEEN ATOMS ',I2,
     *            ' AND ',I2,' IS ',F6.2)
        END IF
      ELSE
         RADIUS(I)=0.7D+00
      END IF
      END DO
C
      CALL FLSHBF(IW)
C  ANGULAR GRID
      CALL LEBPTW(XYZ,XYZLEB,W,N_A,MXDIM)
C
C  CONSTRUCT GRIDS AROUND EACH ATOM
      K=0
      DO N=1,NPTSTN
         IZNUC = INT(ZANSTN(N))
         ALPHA = RSCALE*SLTRAD(IZNUC)
         CALL RADGRD(ALPHA,WR,R,N_R,M_R)
         NSTART(N)=K+1
         DO I=1,N_R-1
            DO J=1,N_A
               K=K+1
               GRID(1,K)=R(I)*XYZLEB(1,J)+STNPNT(2,N)
               GRID(2,K)=R(I)*XYZLEB(2,J)+STNPNT(3,N)
               GRID(3,K)=R(I)*XYZLEB(3,J)+STNPNT(4,N)
               GRID(4,K)=FOUR*PI*WR(I)*W(J)
C               WRITE(IW,*) 'GRID(',K,') =',GRID(4,K)
            END DO
         END DO
      END DO
      NSTART(NPTSTN+1)=K+1
      IF (K.GT.NG) THEN
         WRITE(IW,*) 'NOT ENOUGH GRID POINTS ALLOCATED -- ', K,' NEEDED'
         CALL ABRT
      END IF
      NG=K
C
      IF (NPTSTN .GT. 1) THEN
C  ASSIGN WEIGHT TO EACH POINT ACCORDING TO BECKE FORMULA.
C  RR(M,N) IS THE DISTANCE BETWEEN NUCLEI M AND N.
         DO I=1,NPTSTN
            S(I,I)= ONE
            DO J=1,NPTSTN
               IF (J .NE. I) THEN
                   RR(I,J)=SQRT((STNPNT(2,J)-STNPNT(2,I))**2+
     *                     (STNPNT(3,J)-STNPNT(3,I))**2+
     *                     (STNPNT(4,J)-STNPNT(4,I))**2)
                   CHI=RADIUS(I)/RADIUS(J)
                   U=(CHI-ONE)/(CHI+ONE)
                   AA(I,J)=U/(U**2-ONE)
               END IF
            END DO
         END DO
C
C         WRITE(IW,*)'BECKE SMOOTHING PARAMETER = ', K_MU
         NATOM=0
         DO K=1,NG
            IF (K .GE. NSTART(NATOM+1)) NATOM=NATOM+1
C      THIS GRID POINT BELONGS TO ATOM NATOM
C      GR(M) IS DISTANCE FROM ATOM M
            DO M=1,NPTSTN
               GR(M)=SQRT((GRID(1,K)-STNPNT(2,M))**2+
     *              (GRID(2,K)-STNPNT(3,M))**2+
     *              (GRID(3,K)-STNPNT(4,M))**2)
C               WRITE(IW,*) 'GR=',GR(M)
            END DO
            DO MA=1,NPTSTN
               DO MB=1,NPTSTN
                  IF (MA .EQ. MB) GOTO 100
                  XMU_AB=(GR(MA)-GR(MB))/RR(MA,MB)
C     CHANGE VARIABLE TO ACCOUNT FOR ATOM RADII
                  XNU_AB=XMU_AB+AA(MA,MB)*(ONE-XMU_AB**2)
C     SMOOTHING FUNCTION
                  F=XNU_AB
C               WRITE(IW,*) 'F=',F
                  DO I=1,K_MU
                     F=F*(1.5D+00-0.5D+00*F*F)
                  END DO
                  S(MA,MB)=0.5D+00*(ONE-F)
 100           CONTINUE
               END DO
               PP(MA)=ONE
               DO II=1,NPTSTN
                  PP(MA)=PP(MA)*S(MA,II)
               END DO
C               WRITE(IW,*) 'S=',S
C               WRITE(IW,*) 'PP(',MA,')=',PP(MA)
            END DO
            SUMPP = ZERO
            DO II=1,NPTSTN
               SUMPP=SUMPP+PP(II)
            END DO
C            WRITE(IW,*) 'SUMPP=',SUMPP
            WEIGHT=PP(NATOM)/SUMPP
C            WRITE(IW,*) 'WEIGHT=',WEIGHT
            GRID(4,K)=GRID(4,K)*WEIGHT
C            IF (DEBUG) THEN
C               IF ((GRID(1,K)**2+GRID(2,K)**2).LE.1.0D-06) THEN
C                   WRITE(IW,*) NATOM, GRID(:,G), WEIGHT
C                   WRITE(IW,*) (S(I,:), I=1,NPTSTN), PP
C               END IF
C            END IF
         END DO
      ENDIF
C
C  SUPPRESS GRID POINTS WITH WEIGHTS BELOW CUTOFF
C Q=0
C DO G=1,NG
C   IF (G .GE. NSTART(P+1)) THEN
C     NSTART(P+1)=Q+1
C     P=P+1
C     IF (GRID(4,G).GT.CUTOFF) THEN
C       Q=Q+1
C       GRID(:,Q)=GRID(:,G)
C     END IF
C   END IF
C END DO
C NSTART(NPTSTN+1)=Q+1
C NG=Q
C
      RETURN
      END
C
C*MODULE PRPPOP  *DECK RADGRD
      SUBROUTINE RADGRD(ALPHA,W_R,R,N_R,M_R)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION W_R(N_R), R(N_R)
C
C     RADIAL GRID
C
      F=M_R*N_R*ALPHA**3
      DO I=1,N_R-1
         DUMMY = I
         DUMMY2= N_R-I
         R(I)=ALPHA*(DUMMY/DUMMY2)**M_R
         W_R(I)=F*(DUMMY**(3*M_R-1))/(DUMMY2**(3*M_R+1))
      END DO
C
      RETURN
      END
