C 25 Mar 08 - JHS - ADD PBFLAB TO COMMON BLOCK NUCMON
C 02 AUG 07 - JHS - NEW ROUTINE TO CALCULATE NEO PROPERTIES
C
C*MODULE NEOPRP  *DECK PRPNEO
      SUBROUTINE PRPNEO(PRPTYP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*4 PRPTYP
C
      LOGICAL GOPARR,DSKWRK,MASWRK,VTSCAL,VIROK,LVCLN,DETDEN
      logical CANONC,FCORE,FORS,NOCI,EKT,LINSER
      DOUBLE PRECISION METHOD
C
      PARAMETER (MXATM=2000, MXRT=100, MXNORO=250)
C
      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 /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 /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /VIRIAL/ SCALTE,SCALTT,GVIR,VTCONV,MAXVT,VTSCAL,VIROK,LVCLN
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA AMCSCF/8HMCSCF   /
      DATA     C1/8HC1      /
      DATA  ALDET/8HALDET   /
      DATA   PROP/8HPROP    /
      DATA  CHECK/8HCHECK   /
C
C     THE FIRST HALF OF THIS ROUTINE TAKES CARE OF THE
C     ELECTRONIC PROPERTIES WITHIN THE NEO FRAMEWORK. THIS
C     PORTION IS A MODIFICATION OF THE SUBROUTINE PROPTYP.
C     THE SECOND HALF CALCULATES THE NUCLEAR PROPERTIES
C
C      PRESENTLY PROPERTIES NOT AVAILABLE FOR NEO-MCSCF
C
      IF(SCFTYP.EQ.AMCSCF) RETURN
      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.'MP2 ') WRITE(IW,9240)
         IF(PRPTYP.EQ.'CI E') WRITE(IW,9250) CITYP,IROOT
         IF(PRPTYP.EQ.'CI R') WRITE(IW,9260) CITYP,IROOT
      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
      IF(RUNTYP.EQ.PROP)  GO TO 101
C
      CALL VALFM(LOADFM)
      LDA    = LOADFM + 1
      LDB    = LDA    + L2
      LONEEI = LDB    + L2
      LAST   = LONEEI + L2
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 100
C
      CALL NEO_ENGANL(X(LDA),X(LDB),X(LONEEI),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 -NEO_ENGANL-.
C
      IF (VTSCAL) RETURN
C
C        ----- SYMMETRICALLY TRANSFORM VECTORS AND DENSITY MATRIX -----
C
      CALL LOWDIN
C
C        ----- ELECTRONIC MULLIKEN AND LOWDIN POPULATION ANALYSIS -----
C
      CALL MULKEN_NEO
C
C        ----- QM NUCLEAR MULIKEN AND LOWDIN POPULATION ANALYSIS -----
C
      CALL NUCPROP
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.NE.0) CALL ELMOMC_NEO
      IF(IEDEN.NE.0) CALL ELDENC
cjhs      IF(IEPOT.NE.0) CALL ELPOTC
cjhs      IF(IEFLD.NE.0) CALL ELFLDC
C
C     ----- GET EXPECTATION VALUE OF QM NUCLEI POSITIONS -----
C
      CALL NUCEXPCT
C
      IF (MASWRK) WRITE(IW,9010)
      CALL TIMIT(1)
      RETURN
C
 9010 FORMAT(1X,'...... END OF PROPERTY EVALUATION ......')
 9210 FORMAT(/2X,68(1H-)/
     *        2X,'PROPERTY VALUES FOR THE NEO-',A5,
     *           ' SELF-CONSISTENT FIELD WAVEFUNCTION'/
     *        2X,68(1H-))
 9220 FORMAT(/5X,69(1H-)/
     *        5X,'PROPERTIES FOR THE ',A8,' DFT FUNCTIONAL ',
     *           '(',A4,' TYPE) DENSITY MATRIX'/
     *        5X,69(1H-))
 9240 FORMAT(/5X,49(1H-)/
     *        5X,'NEO-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-))
      END
C
C*MODULE NEOPRP  *DECK NEO_ENGANL
      SUBROUTINE NEO_ENGANL(DA,DB,ONEEI,L1,L2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL BETA,GOPARR,DSKWRK,MASWRK,VTSCAL,VIROK,LVCLN
C
      DIMENSION DA(L2),DB(L2),ONEEI(L2)
C
      PARAMETER (MXATM=2000, MXRT=100)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VNE,VEE,VTOT,TKIN,ESTATE(MXRT),STATN,EDFT(2)
      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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NUCENG/ EELNU,ENUNU
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /VIRIAL/ SCALTE,SCALTT,GVIR,VTCONV,MAXVT,VTSCAL,VIROK,LVCLN
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA UHF,ROHF,GVB/8HUHF     ,8HROHF    ,8HGVB     /
      DATA RNONE/8HNONE    /
C
C     COMPUTE SELECTED ENERGY QUANTITIES, AND VIRIAL RATIO
C     M.W. SCHMIDT     NDSU   MAR 18,1983
C
      BETA = (SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF  .OR.  SCFTYP.EQ.GVB)
     *     .AND.  CITYP.EQ.RNONE
      CALL WFNDEN(DA,DB,L2)
      CALL DAREAD(IDAF,IODA,ENUCR,MXRT+15,2,0)
C
C     ----- E1 FROM BARE NUCLEUS HAMILTONIAN -----
C
      CALL DAREAD(IDAF,IODA,ONEEI,L2,11,0)
      E1 = TRACEP(DA,ONEEI,L1)
      IF(BETA) E1 = E1 + TRACEP(DB,ONEEI,L1)
C
C     ----- TKIN FROM KINETIC ENERGY INTEGRALS
C
      CALL DAREAD(IDAF,IODA,ONEEI,L2,13,0)
      TKIN = TRACEP(DA,ONEEI,L1)
      IF(BETA) TKIN = TKIN + TRACEP(DB,ONEEI,L1)
C
C     ----- WAVEFUNCTION NORMALIZATION FROM OVERLAP INTEGRALS
C
      CALL DAREAD(IDAF,IODA,ONEEI,L2,12,0)
      PSINRM = TRACEP(DA,ONEEI,L1)
      IF(BETA) PSINRM = PSINRM + TRACEP(DB,ONEEI,L1)
      PSINRM = PSINRM/NE
C
C     ----- GET MEMORY FOR NEO ENERGY ANALYSIS -----
C
      N1 =NUMNB
      N2 = (N1*N1+N1)/2

      CALL VALFM(LOADFM)
      LH1   = LOADFM + 1
      NDA   = LH1    + N2
      LTNUC = NDA    + N2
      NOVLP = LTNUC  + N2
      LAST  = NOVLP  + N2
      NEEDNUC = LAST - LOADFM - 1
C
      CALL GETFM(NEEDNUC)
C
C     ----- NUCLEAR WAVEFUNCTION NORMALIZATION FROM OVERLAP INTEGRALS
C
      CALL DAREAD(IDAF,IODA,X(NOVLP),N2,441,0)
      CALL DAREAD(IDAF,IODA,X(NDA),N2,445,0)
      PSINUCNRM = TRACEP(X(NDA),X(NOVLP),N1)
      PSINUCNRM = PSINUCNRM/NQMNUC
C
C     ----- 1NUC TERM -----
C
      CALL DAREAD(IDAF,IODA,X(LH1),N2,440,0)
      ENUC1 = TRACEP(X(NDA),X(LH1),N1)
C
C     ----- NUCLEAR KINETIC ENERGY INTEGRALS -----
C
      CALL DAREAD(IDAF,IODA,X(LTNUC),N2,442,0)
      TKNUC = TRACEP(X(NDA),X(LTNUC),N1)
C
C     ----- EVERYTHING ELSE BY SUBTRACTION -----
C
      EMIX  = EELNU
      ENUC2 = ENUNU
      E2 = ETOT - E1 - ENUC1 -ENUC2 -EMIX -ENUCR
      VNE = E1 - TKIN
      VNN = ENUCR
      VEE = E2
      VTOT = VNE + VNN + VEE
      VIRIAL = -VTOT/TKIN
C
C     ----- PRINT OUT THE RESULTS
C
      IF(MASWRK  .AND.  .NOT.VTSCAL) THEN
         WRITE(IW,9000)
         WRITE(IW,9010) PSINRM*PSINUCNRM
         WRITE(IW,9020) E1,E2,ENUC1,ENUC2,EMIX,ENUCR,ETOT
         WRITE(IW,9030) VEE,VNE,VNN,VTOT,TKIN,TKNUC,VIRIAL
      END IF
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      CALL RETFM(NEEDNUC)
C
      RETURN
C
 9000 FORMAT(/20X,17(1H-)/20X,17HENERGY COMPONENTS/20X,17(1H-)/)
 9010 FORMAT( 1X,'              NEO WAVEFUNCTION NORMALIZATION =',
     * F19.10)
 9020 FORMAT(/1X,'                         ONE ELECTRON ENERGY =',
     * F19.10/
     *        1X,'                         TWO ELECTRON ENERGY =',
     * F19.10/
     *        1X,'                          ONE NUCLEUS ENERGY =',
     * F19.10/
     *        1X,'                          TWO NUCLEUS ENERGY =',
     * F19.10/
     *        1X,'               MIXED ELECTRON-NUCLEUS ENERGY =',
     * F19.10/
     *        1X,'          CLASSICAL NUCLEAR REPULSION ENERGY =',
     * F19.10/
     *       38X,28(1H-)/
     *        1X,'                                TOTAL ENERGY =',
     * F19.10)
 9030 FORMAT(/1X,'          ELECTRON-ELECTRON POTENTIAL ENERGY =',
     * F19.10/
     *        1X,' CLASSICAL NUCLEUS-ELECTRON POTENTIAL ENERGY =',
     * F19.10/
     *        1X,'  CLASSICAL NUCLEUS-NUCLEUS POTENTIAL ENERGY =',
     * F19.10/
     *       38X,28(1H-)/
     *        1X,'           TOTAL ELECTRONIC POTENTIAL ENERGY =',
     * F19.10/
     *        1X,'             TOTAL ELECTRONIC KINETIC ENERGY =',
     * F19.10/
     *        1X,'                TOTAL NUCLEAR KINETIC ENERGY =',
     * F19.10/
     *        1X,'                          VIRIAL RATIO (V/T) =',
     * F19.10)
      END
C
C*MODULE NEOPRP  *DECK MULKEN_NEO
      SUBROUTINE MULKEN_NEO
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 /NUCPRP/ EPMUL(MXATM),EPLOW(MXATM)
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      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
      NTOTCN = NQMNUC + NTAUXB
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
      I54 = I53+NTOTCN
      I55 = I54+NTOTCN
C
C     ----- GET MEMORY -----
C
      LASTC = I55
      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)
C
C     ----- GET QM NUCLEAR CHARGES -----
C
      IQMN = I53
      IQLN = I54
      CALL NUC_CHARGE(X(IQMN),X(IQLN),NTOTCN,NUMNB)
C
      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)
         EPMUL(I) = PM
         QM = ZAN(I) - PM
         PL = X(I-1+ILAC)
         EPLOW(I) = PL
         QL = ZAN(I) - PL
C
C       ----- ADD QM NUCLEAR CONTRIBUTION TO CHARGES -----
C
         DO J = 1,NTOTCN
           IF(I.EQ.IQMNUC(J)) THEN
            QM = QM + X(IQMN-1+J)
            QL = QL + X(IQLN-1+J)
           ENDIF
         ENDDO
C
         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
C*MODULE NEOPRP   *DECK NUCPROP
      SUBROUTINE NUCPROP
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MCSCF,NEOSCF,NEOCI
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NEOMTD/ NEOSCF,NEOCI
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DATA DIAGZN,RONHF/8HDIAGZN  ,8HRONHF   /
      DATA MCSCF/8HMCSCF   /
C
C   DETERMINE IF POPULATION DENSITIES CAN BE CALCULATED FOR RUN
C
         IF(NEOSCF.EQ.RONHF.OR.NEOSCF.EQ.DIAGZN.OR.NEOSCF.EQ.MCSCF) THEN
C
           L1 = NUMNB
           L2 = (L1*L1+L1)/2
           L3 = L1*L1
           NTOTCN = NQMNUC + NTAUXB
           NTOTCN2= (NTOTCN*NTOTCN+NTOTCN)/2
C
C     ----- GET MEMORY -----
C
           CALL VALFM(LOADFM)
           NOVER  = LOADFM + 1
           NVEC   = NOVER  + L2
           NT     = NVEC   + L3
           NMOATM = NT     + L1
           NOCC   = NMOATM + L1*NTOTCN
           NADEN  = NOCC   + L1
           NGPOPAO= NADEN  + L2
           LPOP   = NGPOPAO+ L1
           LSN2   = LPOP   + L1
           LAD2   = LSN2   + L3
           LV2    = LAD2   + L3
           LADTMP = LV2    + L3
           LVIND  = LADTMP + L2
           LVTRAN = LVIND  + L3
           LSEIGE = LVTRAN + L3
           LSON   = LSEIGE + L1
           LTELOW = LSON   + L3
           LFIN   = LTELOW + L3
           LWRK   = LFIN   + L3
           LIWRK  = LWRK   + L1*8
           LOVERD = LIWRK  + L1
           LOVD   = LOVERD + NTOTCN*NTOTCN
           LIMINB = LOVD   + NTOTCN2
           LIMAXB = LIMINB + NTOTCN
           LTOTM  = LIMAXB + NTOTCN
           LTOTL  = LTOTM  + NTOTCN
           LAST   = LTOTL  + NTOTCN
           NEEDNUP= LAST   - LOADFM - 1
C
        CALL GETFM(NEEDNUP)
C
C     ----- READ IN THE NUC MOLECULAR ORBITALS AND THE NUC DENSITY -----
C
        CALL DAREAD(IDAF,IODA,X(NVEC),L3,444,0)
        CALL DAREAD(IDAF,IODA,X(NOVER),L2,441,0)
C
C     ----- CALCULATE NUCLEAR POPULATION DENSITY -----
C
      CALL NUC_POPDEN(X(NOVER),X(NVEC),X(NT),X(NMOATM),X(NOCC),
     *                   X(NADEN),X(NGPOPAO),X(LPOP),L1,L2,NTOTCN,
     *                   X(LSN2),X(LAD2),X(LV2),X(LADTMP),
     *                   X(LVIND),X(LVTRAN),X(LSEIGE),X(LSON),X(LTELOW),
     *                   X(LFIN),X(LWRK),X(LIWRK),X(LOVERD),X(LOVD),
     *                   X(LIMINB),X(LIMAXB),X(LTOTM),X(LTOTL))
C
         CALL RETFM(NEEDNUP)
C
      WRITE(IW,9000)
         ELSE
C
      WRITE(IW,*)"NUCLEAR PROP DOES NOT SUPPORT THE NEOSCF DECLARED"
             CALL ABRT
C
         ENDIF
C
       RETURN
 9000 FORMAT(/1X,'..... END OF NUCLEAR POPULATION ANALYSES .....'/)
       END
C
C*MODULE NEOPRP   *DECK NUC_POPDEN
      SUBROUTINE NUC_POPDEN(SNUC,VNUC,TNUC,OATNUC,OCCNUC,ADEN,GPOPAO,
     *                      POPLOW,N1,N2,NTOTCN,SNUC2D,ADEN2D,
     *                      VNUC2,ADENTEMP,VINDIAG,VTRAN2,
     *                      SEIGEN,SONEH,TEMPLOW,FINLOW,WRK,IWRK,
     *                      OVERDEN,OVERD,IMINB,IMAXB,TOTMUL,TOTLOW)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 QNUN,QNN
      CHARACTER*10 PBFLAB
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
      PARAMETER (MXAO=8192, MXATM=2000)
C
      DIMENSION IA(MXAO),IWRK(N1),IMINB(NTOTCN),IMAXB(NTOTCN)
      DIMENSION SNUC(N2),VNUC(N1,N1),TNUC(N1),OATNUC(NTOTCN,N1),
     *          OCCNUC(N1),ADEN(N2),GPOPAO(N1),POPLOW(N1),
     *          SNUC2D(N1,N1),ADEN2D(N1,N1),
     *          VNUC2(N1,N1),ADENTEMP(N2),VINDIAG(N1,N1),
     *          VTRAN2(N1,N1),SEIGEN(N1),SONEH(N1,N1),TEMPLOW(N1,N1),
     *          FINLOW(N1,N1),WRK(N1,8),OVERDEN(NTOTCN,NTOTCN),
     *          OVERD((NTOTCN*NTOTCN+NTOTCN)/2),TOTMUL(NTOTCN),
     *          TOTLOW(NTOTCN)
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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NUCMON/ QNUN(20),QNN(20),PBFLAB(MXAO)
      COMMON /NUCPRP/ EPMUL(MXATM),EPLOW(MXATM)
      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)
C
      NOCCNORB=NNA  ! NUMBER OF OCCUPIED NUCLEAR ORBS
      NBAS_ATOM = NUMNB/NTOTCN  ! NUMBER OF BASIS FUNC PER NUCLEUS CENTER
C
C    POPULATION ANALYSIS ROUTINE FOR QUANTUM PROTONS
C
      DO 210 I = 1,NTOTCN
         OCCNUC(I) = ONE
  210 CONTINUE
C
      DO 220 I = NTOTCN+1,NUMNB
         OCCNUC(I) = ZERO
  220 CONTINUE
C
      DO 100 I = 1,N1
         IA(I) = (I*I-I)/2
  100 CONTINUE
C
C     I LOOPS OVER ALL THE OCCUPIED NUCLEAR ORBITALS
C     J LOOPS OVER ALL THE NUCLEAR MO's
C
         DO 340 I = 1,NOCCNORB
         DO 320 J = 1,N1
            SUM = ZERO
            DO 310 K = 1,N1
               JK = IA(J)+K
               IF (K .GT. J) JK = IA(K)+J
               SUM = SUM+VNUC(K,I)*SNUC(JK)*VNUC(J,I)
  310       CONTINUE
            TNUC(J) = SUM*OCCNUC(I)
  320    CONTINUE
         DO 330 J = 1,N1
            VNUC2(J,I) = TNUC(J)
  330    CONTINUE
  340 CONTINUE
C
      WRITE(IW,9000)
      IF((NPRINT.GE.8) .AND. MASWRK) THEN
      WRITE(IW,9005)
      CALL PREVNL(VNUC2,OCCNUC,NOCCNORB,N1,N1)
      ENDIF
C
C     ----- CONDENSE TO ATOMS -----
C
      DO 440 K=1,NOCCNORB
            IMIN = 0
            IMAX = 0
         DO 430 IAT=1,NTOTCN
            SUM = ZERO
            IMIN = IMIN + NBAS_ATOM
            IF(IAT.EQ.1) IMIN =1
            IMAX = IMAX + NBAS_ATOM
C
            DO 420 I=IMIN,IMAX
               SUM = SUM + VNUC2(I,K)
  420       CONTINUE
            OATNUC(IAT,K) = SUM
  430    CONTINUE
  440 CONTINUE
C
C     ----- PRINT ATOM POP IN EACH MO -----
C
       WRITE(IW,9010)
       CALL PREVNL_NU(OATNUC,OCCNUC,IQMNUC,NOCCNORB,NTOTCN,NTOTCN)
C
C     ----- GET ALPHA DENSITY MATRIX -----
C
        CALL DAREAD(IDAF,IODA,ADEN,N2,445,0)
        CALL CPYTSQ(SNUC,SNUC2D,N1,1)
        CALL CPYTSQ(ADEN,ADEN2D,N1,1)
C
        DO I = 1,N2
          ADENTEMP(I) = ADEN(I)
        ENDDO
C
C     ----- GET OVERLAP POPULATION -----
C
        CALL OVLPOP(ADENTEMP,SNUC,N2)
C
C     ----- CALCULATE TOTAL GROSS POPULATION IN AO*S -----
C
      CALL GROSSC(ADENTEMP,GPOPAO,IA,N1)
C
C     ----- CALCULATE LOWDIN POP -----
C
      CALL NUC_LOWDIN(SNUC,POPLOW,ADEN2D,VINDIAG,VTRAN2,SEIGEN,
     *                SONEH,TEMPLOW,FINLOW,WRK,IWRK,N1,N2)
C
C     ----- CALCULATE OVERLAP -----
C
         CALL CPYTSQ(ADENTEMP,ADEN2D,N1,1)
         DO I = 1,NTOTCN
           IF (I.EQ.1) THEN
           IMINB(I)=1
           IMAXB(I)=NBAS_ATOM
C
           ELSE
           IMINB(I)=IMINB(I-1) + NBAS_ATOM
           IMAXB(I)=IMAXB(I-1) + NBAS_ATOM
           ENDIF
C
         ENDDO
C
       DO 540 K=1,NTOTCN
         DO 530 IAT=1,NTOTCN
            SUM = ZERO
C
               DO I = IMINB(K),IMAXB(K)
                 DO J = IMINB(IAT),IMAXB(IAT)
                   SUM = SUM + ADEN2D(I,J)
                 ENDDO
               ENDDO
C
            OVERDEN(IAT,K) = SUM
  530    CONTINUE
  540 CONTINUE
      CALL CPYSQT(OVERDEN,OVERD,NTOTCN,1)
C
         WRITE(IW,9040)
         WRITE(IW,9060) (I,PBFLAB(I+NUM),GPOPAO(I),POPLOW(I),I=1,N1)
         WRITE(IW,9035)
         CALL PRTRI_NU(OVERD,IQMNUC,NTOTCN)
C
C     ----- TOTAL ATOMIC POPULATIONS -----
C
             L = 0
          DO I = 1,NTOTCN
                SUM1 = ZERO
                SUM2 = ZERO
            DO J = L + 1,L + NBAS_ATOM
                SUM1 = SUM1 +GPOPAO(J)
                SUM2 = SUM2 + POPLOW(J)
            ENDDO
          TOTMUL(I) = SUM1
          TOTLOW(I) = SUM2
          L = L + NBAS_ATOM
          ENDDO
       WRITE(IW,9090)
       WRITE(IW,9120)
       DO I=1,NTOTCN
         PM = TOTMUL(I)
         PL = TOTLOW(I)
         ILOC = IQMNUC(I)
         QM = PM - EPMUL(ILOC)
         QL = PL - EPLOW(ILOC)
         WRITE(IW,9180) ILOC,QNN(I),BNAM(ILOC),
     *                 PM,QM,PL,QL
       ENDDO
C
 9000 FORMAT(//1X,'..... NUCLEAR POPULATION ANALYSES .....')
 9005 FORMAT(/5X,'AO MULLIKEN POPULATIONS IN EACH MOLECULAR ORBITAL')
 9010 FORMAT(/5X,'MULLIKEN ATOMIC POPULATION IN EACH NUC MOLECULAR',
     *       ' ORBITAL')
 9040 FORMAT(/15X,'----- POPULATIONS IN EACH NUC AO -----'/
     *        25X,'      MULLIKEN      LOWDIN')
 9035 FORMAT(/10X,'----- MULLIKEN NUCLEI OVERLAP POPULATIONS -----'/
     *        10X,'(OFF-DIAGONAL ELEMENTS NEED TO BE MULTIPLIED BY 2)')
 9060 FORMAT(10X,I5,2X,A10,2F12.5)
 9090 FORMAT(/10X,'TOTAL MULLIKEN AND LOWDIN NUC ATOMIC POPULATIONS')
 9120 FORMAT(7X,'NUC ',9X,'MULL.POP.    CHARGE',10X,
     *           'LOW.POP.     CHARGE')
 9180 FORMAT(1X,I4,1X,A8,A2,2F12.6,5X,2F12.6)
C
      RETURN
C
      END
C
C*MODULE NEOPRP   *DECK NUC_LOWDIN
      SUBROUTINE NUC_LOWDIN(SNLOW,POPLOW,DEN,VINDIAG,VTRAN2,
     *               SEIGEN,SONEH,TEMPLOW,FINLOW,WRK,IWRK,
     *               N1,N2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION IWRK(N1)
      DIMENSION SNLOW(N2),POPLOW(N1),DEN(N1,N1),VINDIAG(N1,N1),
     *          VTRAN2(N1,N1),SEIGEN(N1),SONEH(N1,N1),TEMPLOW(N1,N1),
     *          FINLOW(N1,N1),WRK(N1,8)
C
      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
C
C  ---- DIAGONALIZE OVERLAP MATRIX USING JACOBI (kdiag =3) -----
C
      KDIAGSAV = KDIAG
      KDIAG = 3
      CALL GLDIAG(N1,N1,N1,SNLOW,WRK,SEIGEN,VINDIAG,
     *               IERR,IWRK)
      IF(IERR .NE. 0) CALL ABRT
      KDIAG = KDIAGSAV
C
C  ---- TAKE SQUARE ROOT OF S EIGENVALUES ----
C
         DO  I = 1,n1
             SEIGEN(I) = DSQRT(SEIGEN(I))
         ENDDO
C
C  ------ TRANSPOSE THE S EIGENVECTOR --------
C
      DO I =1,N1
         DO J = 1,N1
           VTRAN2(I,J)=VINDIAG(J,I)
         ENDDO
      ENDDO
C
C  ----- FORM SYMMETRIC  ORTHONORMAL ORBITALS -----
C              X = S^1/2 = U*(s^1/2)
C
      DO 410 J=1,N1
         EJ = SEIGEN(J)
         DO 400 I=1,N1
            QIJ=VINDIAG(I,J)
            QIJ=QIJ*EJ
           VINDIAG(I,J)=QIJ
  400    CONTINUE
  410  CONTINUE
C
        CALL MRARBR(VINDIAG,N1,N1,N1,VTRAN2,N1,
     *               N1,SONEH,N1)
C
C  ---- MULTIPLY S^1/2*P*S^1/2 ----
C
         CALL MRARBR(SONEH,N1,N1,N1,DEN,N1,
     *               N1,TEMPLOW,N1)
         CALL MRARBR(TEMPLOW,N1,N1,N1,SONEH,N1,
     *               N1,FINLOW,N1)
C
         DO I =1,N1
           POPLOW(I) = FINLOW(I,I)
         ENDDO
C
       RETURN
       END
C
C*MODULE NEOPRP   *DECK ELMOMC_NEO
      SUBROUTINE ELMOMC_NEO
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PAPER,PUNCH,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ITAG(19)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /FMCOM / X(1)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      PARAMETER (DFAC=2.541766D+00)
      PARAMETER (ONE=1.0D+00)
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HELMOMC  /
      DATA ELMOM/8HELMOM   /
      DATA ITAG/4H   X,4H   Y,4H   Z,4H  XX,4H  YY,4H  ZZ,4H  XY,
     *          4H  XZ,4H  YZ,4H XXX,4H YYY,4H ZZZ,4H XXY,4H XXZ,
     *          4H XYY,4H YYZ,4H XZZ,4H YZZ,4H XYZ/
C
C     ----- CALCULATE ELECTROSTATIC MOMENTS -----
C
      PAPER = IEMOUT.GE.0
      IF(EXETYP.EQ.DEBUG)  PAPER=.TRUE.
      IF(EXETYP.EQ.DBUGME) PAPER=.TRUE.
      PAPER = PAPER .AND. MASWRK
C
      PUNCH = IEMOUT.NE.0 .AND. NFG.EQ.0
      PUNCH = PUNCH .AND. MASWRK
C
      IF(PAPER) WRITE(IW,920)
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
C
                     NVAL=0
      IF(IEMOM.EQ.1) NVAL=3
      IF(IEMOM.EQ.2) NVAL=9
      IF(IEMOM.EQ.3) NVAL=19
      IF(NVAL.EQ.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ELMOMC: IMPOSSIBLE -IEMOM- VALUE'
         CALL ABRT
      END IF
C
C     ----- GROW MEMORY -----
C
      CALL VALFM(LOADFM)
      IELM   = LOADFM + 1
      IEMW   = IELM   + NVAL*L2
      IDENSA = IEMW   + NVAL*784
      IDENSB = IDENSA + L2
      LAST   = IDENSB + L2
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      IF (EXETYP.EQ.CHECK) THEN
         CALL RETFM(NEED)
         RETURN
      END IF
C
C     ----- GET TOTAL DENSITY MATRIX -----
C
      CALL DENDD1(X(IDENSA),X(IDENSB),L2)
C
C     ----- LOOP OVER POINTS FOR MULTIPOLE EXPANSION -----
C
      IPOINT = 0
  210 IPOINT = IPOINT + 1
      ISTAT = 0
      CALL POINT(ELMLOC,IPOINT,XP,YP,ZP,ISTAT)
      IF(ISTAT.LT.0) CALL ABRT
      IF(ISTAT.GT.0) THEN
         CALL RETFM(NEED)
         RETURN
C        ******
C
      END IF
C
C     ----- COMPUTE MULTIPOLE INTEGRALS OVER BASIS FUNCTIONS -----
C
      CALL PRCALC(ELMOM,X(IELM),X(IEMW),NVAL,L2,GOPARR)
C
      IF(GOPARR) CALL DDI_GSUMF(914,X(IELM),3*L2)
C
      IF(IEMINT.EQ.0 .OR. .NOT.PAPER) GO TO 240
                              IMAX=0
         IF(ABS(IEMINT).EQ.1) IMAX=3
         IF(ABS(IEMINT).EQ.2) IMAX=9
         IF(ABS(IEMINT).EQ.3) IMAX=19
                          IMIN=1
         IF(IEMINT.EQ.-2) IMIN=4
         IF(IEMINT.EQ.-3) IMIN=10
         IF(GOPARR) CALL DDI_GSUMF(915,X(IELM+3*L2),(IMAX-3)*L2)
         IF(IMAX.EQ.0) THEN
            IF(MASWRK) WRITE(IW,*) 'ELMOMC: IMPOSSIBLE -IEMINT- VALUE'
            CALL ABRT
         END IF
         INDEX = IELM+(IMIN-1)*L2
         DO 230 I=IMIN,IMAX
            IF(MASWRK) WRITE(IW,990) ITAG(I)
            CALL PRTRIL(X(INDEX),L1)
            INDEX = INDEX + L2
  230    CONTINUE
         IF(GOPARR) CALL DSCAL((IMAX-3)*L2,ONE/NPROC,X(IELM+3*L2),1)
  240 CONTINUE
C
C     ---- STORE DIPOLE INTEGRALS ON DAF -----
C
      CALL DAWRIT(IDAF,IODA,X(IELM      ),L2,95,0)
      CALL DAWRIT(IDAF,IODA,X(IELM+   L2),L2,96,0)
      CALL DAWRIT(IDAF,IODA,X(IELM+ 2*L2),L2,97,0)
C
C                 CALCULATE ELECTROSTATIC MOMENTS
C
C                         MONOPOLE  (FROM OVERLAP INTEGRALS)
      CALL DAREAD(IDAF,IODA,X(IDENSB),L2,12,0)
      CHARGE = -TRACEP(X(IDENSA),X(IDENSB),L1)
C                         DIPOLE
      DMX = -TRACEP(X(IDENSA),X(IELM     ),L1)
      DMY = -TRACEP(X(IDENSA),X(IELM+  L2),L1)
      DMZ = -TRACEP(X(IDENSA),X(IELM+2*L2),L1)
C                         QUADRUPOLE
      IF(IEMOM.GE.2) THEN
          WRITE(IW,925)
      ENDIF
C         ELTEMP( 1) = -TRACEP(X(IDENSA),X(IELM+ 3*L2),L1)
C         ELTEMP( 2) = -TRACEP(X(IDENSA),X(IELM+ 4*L2),L1)
C         ELTEMP( 3) = -TRACEP(X(IDENSA),X(IELM+ 5*L2),L1)
C         ELTEMP( 4) = -TRACEP(X(IDENSA),X(IELM+ 6*L2),L1)
C         ELTEMP( 5) = -TRACEP(X(IDENSA),X(IELM+ 7*L2),L1)
C         ELTEMP( 6) = -TRACEP(X(IDENSA),X(IELM+ 8*L2),L1)
C         NSUM = NSUM+6
C     END IF
C                         OCTUPOLE
C      IF(IEMOM.GE.3) THEN
C        ELTEMP( 7) = -TRACEP(X(IDENSA),X(IELM+ 9*L2),L1)
C        ELTEMP( 8) = -TRACEP(X(IDENSA),X(IELM+10*L2),L1)
C        ELTEMP( 9) = -TRACEP(X(IDENSA),X(IELM+11*L2),L1)
C        ELTEMP(10) = -TRACEP(X(IDENSA),X(IELM+12*L2),L1)
C        ELTEMP(11) = -TRACEP(X(IDENSA),X(IELM+13*L2),L1)
C        ELTEMP(12) = -TRACEP(X(IDENSA),X(IELM+14*L2),L1)
C        ELTEMP(13) = -TRACEP(X(IDENSA),X(IELM+15*L2),L1)
C        ELTEMP(14) = -TRACEP(X(IDENSA),X(IELM+16*L2),L1)
C        ELTEMP(15) = -TRACEP(X(IDENSA),X(IELM+17*L2),L1)
C        ELTEMP(16) = -TRACEP(X(IDENSA),X(IELM+18*L2),L1)
C        NSUM = NSUM+10
C     END IF
C
C     IF(GOPARR  .AND.  NSUM.GT.0) CALL DDI_GSUMF(916,ELTEMP,NSUM)
C
C     QMXX = ELTEMP(1)
C     QMYY = ELTEMP(2)
C     QMZZ = ELTEMP(3)
C     QMXY = ELTEMP(4)
C     QMXZ = ELTEMP(5)
C     QMYZ = ELTEMP(6)
C
C     OMXXX = ELTEMP( 7)
C     OMYYY = ELTEMP( 8)
C     OMZZZ = ELTEMP( 9)
C     OMXXY = ELTEMP(10)
C     OMXXZ = ELTEMP(11)
C     OMXYY = ELTEMP(12)
C     OMYYZ = ELTEMP(13)
C     OMXZZ = ELTEMP(14)
C     OMYZZ = ELTEMP(15)
C     OMXYZ = ELTEMP(16)
C
C
C                ADD QM NUCLEAR CONTRIBUTION TO MOMENTS
C
      N1=NUMNB
      N2=(NUMNB*NUMNB+NUMNB)/2
      NVAL=3
C
C     ----- GROW MEMORY -----
C
      CALL VALFM(LOADFM)
      IDENNUC= LOADFM + 1
      INUCOV = IDENNUC+ N2
      INLM   = INUCOV + N2
      INMW   = INLM   + NVAL*N2
      LAST   = INMW   + NVAL*784
      NEEDNUC = LAST - LOADFM - 1
      CALL GETFM(NEEDNUC)
C
C     ---- GET NUCLEAR DENSITY AND OVERLAP MATRICES ----
C
          CALL DAREAD(IDAF,IODA,X(IDENNUC),N2,445,0)
          CALL DAREAD(IDAF,IODA,X(INUCOV),N2,441,0)
C
C                         MONOPOLE
C
       CHARGE = CHARGE + TRACEP(X(IDENNUC),X(INUCOV),N1)
C
C                         DIPOLE
C
      CALL PRCALC_NEO(ELMOM,X(INLM),X(INMW),NVAL,N2,GOPARR)
C
C     ---- STORE QM NUCLEAR DIPOLE INTEGRALS ON DAF -----
C
      CALL DAWRIT(IDAF,IODA,X(INLM      ),N2,464,0)
      CALL DAWRIT(IDAF,IODA,X(INLM+   N2),N2,465,0)
      CALL DAWRIT(IDAF,IODA,X(INLM+ 2*N2),N2,466,0)
C
      DNX = TRACEP(X(IDENNUC),X(INLM),N1)
      DMX = DMX +DNX
      DNY = TRACEP(X(IDENNUC),X(INLM + N2),N1)
      DMY = DMY+DNY
      DNZ = TRACEP(X(IDENNUC),X(INLM + 2*N2),N1)
      DMZ = DMZ +DNZ
C
       CALL RETFM(NEEDNUC)
C
C                 ADD CLASSICAL NUCLEAR CONTRIBUTION TO MOMENTS
C
      DO 280  I=1,NAT
         ZNUC = ZAN(I)
         CHARGE = CHARGE + ZNUC
         XN = C(1,I) - XP
         YN = C(2,I) - YP
         ZN = C(3,I) - ZP
         DMX = DMX + ZNUC*XN
         DMY = DMY + ZNUC*YN
         DMZ = DMZ + ZNUC*ZN
  280 CONTINUE
C
C     ----- CONVERT UNITS -----
C
      DMX = DMX*DFAC
      DMY = DMY*DFAC
      DMZ = DMZ*DFAC
      DIP = SQRT(DMX*DMX+DMY*DMY+DMZ*DMZ)
C
C     ---- OUTPUT RESULTS -----
C
      IF(.NOT.PAPER) GO TO 310
      WRITE(IW,930) IPOINT,XP,YP,ZP,CHARGE,DMX,DMY,DMZ,DIP
  310 CONTINUE
      IF(.NOT.PUNCH) GO TO 210
      WRITE(IP,960) IPOINT,XP,YP,ZP,DMX,DMY,DMZ
      GO TO 210
C
  920 FORMAT(/10X,21('-'),/10X,'ELECTROSTATIC MOMENTS',
     *      /10X,21('-'),/)
  925 FORMAT(/1X,5('*'),'QUADRUPOLE AND OCTUPOLE MOMENTS ARE NOT ',
     *          'AVAILABLE FOR NEO RUNS',5('*')/)
  930 FORMAT(1X,'POINT',I4,11X,'X',11X,'Y',11X,'Z (BOHR)',4X,
     *       'CHARGE'/13X,3F12.6,6X,F6.2,1X,'(A.U.)'/
     *       9X,'DX',10X,'DY',10X,'DZ',9X,'/D/  (DEBYE)'/
     *       1X,4F12.6)
  960 FORMAT(1X,'MOMENTS AT POINT',I5,' X,Y,Z=',3F10.6/
     *       1X,'DIPOLE     ',3F10.6)
  990 FORMAT(/10X,A4,' INTEGRALS'/11X,13(1H-))
      END
C
C*MODULE NEOPRP   *DECK PRCALC_NEO
      SUBROUTINE PRCALC_NEO(PRNAME,XVAL,WINT,NVAL,L2,GPINTS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,NORM,DOUBLE,GPINTS
C
C            EXTERNAL STORAGE REQUIREMENT FOR WINT IS NVAL*784
      DIMENSION XVAL(NVAL*L2),WINT(*)
C
      DIMENSION IX(84),IY(84),IZ(84),JX(84),JY(84),JZ(84),
     *          DIJ(784),IJX(784),IJY(784),IJZ(784)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)

      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NEOSHL/ NGAUSS,NGAUSE,NGAUSN,NTSHEL,NNUCSH
      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 /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ

      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 ELDEN  /8HELDEN   /
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
C       CALCULATE INTEGRALS NECESSARY FOR EVALUATION OF PROPERTIES.
C
C       THE PROPERTY REQUESTED BY -PRNAME- CAN BE (HOLLERITH!)
C             'ELDEN   ', 'ELMOM   ', OR 'ELFLD   '
C       THE NUMBER OF INTEGRAL TYPES IN -NVAL- DEPENDS ON THE VALUE OF
C             IEDEN         IEMOM         IEFLD
C       RESPECTIVELY (SEE $ELDENS, $ELMOM, $ELFLD) BUT TOPS OUT AT
C          IEDEN  NVAL    IEMOM NVAL    IEFLD NVAL
C            1     1        1     3       1    3
C                           2     9       2    9
C                           3    19
C
      NIJ = 784*NVAL
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
C     ----- INITIALIZE PARALLEL -----
C     PARALLELIZATION CONTROLLED BY CALLING ARG -GPINTS- NOT -GOPARR-.
C
C     PARALLELIZATION IS ON THE INTEGRAL EVALUATION, BUT THERE IS NO
C     GLOBAL SUM PERFORMED AT THE BOTTOM OF THE ROUTINE.  THE AMOUNT
C     OF DATA SUMMED OVERWHELMS THE PARALLELIZATION,
C
C        FOR 86 ATOMS AND 500 AOS, 6-31G(D), INTEGRAL TIME ONLY,
C        ON THE "INFERNO" AXP UNIPROCESSOR/MYRINET CLUSTER, FOR
C        DOING MOMENTS TO OCTUPOLE AT CENTER OF MASS (1 TIME) AND
C        DOING DENSITY/FIELD GRADIENT/MEP AT ALL NUCLEI (86 TIMES):
C                         P=1   P=4, NO GSUMF   P=4, W/GSUMF
C          MOMENTS       1.89        0.53          1.63  SEC.
C          DENSITY       2.82        0.83          4.29
C          FIELD GRAD   38.86       10.42         45.20
C          M.ELEC.POT.   6.44        1.61          1.59
C        FOOTNOTE: THE LAST IS NOT DONE IN THIS ROUTINE, SEE EPCALC.
C        THE SAME RUN DOES MULLIKEN ANALYSIS FOR P=1 IN 20.76, AND
C        FOR P=4 IN 21.09, ITS MATRIX WORK HAS NO PARALLEL CONTENT.
C
C     IF THE CALLING ROUTINE CAN SUM A SMALLER AMOUNT OF DATA, SUCH
C     AS TRACING THE DENSITY AGAINST THE INTEGRAL SUBSET ON EACH
C     NODE, THEN SUMMING ONLY THE PROPERTY VALUES, SET -GPINTS- ON.
C     OTHERWISE RUN THIS ROUTINE IN SERIAL WITH -GPINTS- OFF.
C
      IPCOUNT = ME - 1
      IF(GPINTS) CALL VCLR(XVAL,1,NVAL*L2)
C
C                         LOOP OVER SHELLS II
C
      DO 510  II=NSHELL+1,NTSHEL
      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                         LOOP OVER SHELLS JJ
C
          DO 500  JJ=NSHELL+1,II
C
C     ----- GO PARALLEL! (STATIC LOAD BALANCING) -----
C
      IF (GPINTS) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 500
      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
C
      RR     = (XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2
      IANDJ  = II.EQ.JJ
C
      CALL VCLR(WINT,1,NIJ)
C
C             PREPARE INDICES FOR PAIRS OF (I,J) ORBITALS
C
      IJ = 0
      MAX = MAXJ
      DO 110  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
  110 CONTINUE
C
C                     LOOP OVER PRIMITIVES IG
C
      JGMAX = J2
      DO 410 IG=I1,I2
      AI  = EX(IG)
      CSI = CS(IG)
      CPI = CP(IG)
      CDI = CD(IG)
      CFI = CF(IG)
      CGI = CG(IG)
      CHI = CH(IG)
      CII = CI(IG)
C
C                     LOOP OVER PRIMITIVES JG
C
      IF (IANDJ) JGMAX = IG
      DO 400 JG=J1,JGMAX
      AJ  = EX(JG)
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ = CF(JG)
      CGJ = CG(JG)
      CHJ = CH(JG)
      CIJ = CI(JG)
C
      AA  = AI + AJ
      AA1 = ONE/AA
      AX  = (AI*XI + AJ*XJ)*AA1
      AY  = (AI*YI + AJ*YJ)*AA1
      AZ  = (AI*ZI + AJ*ZJ)*AA1
C
      DUM = AI*AJ*RR*AA1
C
      IF(PRNAME.EQ.ELDEN)
     *     DUM = DUM + AA*((XP-AX)**2 + (YP-AY)**2 + (ZP-AZ)**2)
      IF(DUM .GT. TOL) GO TO 400
C
      FAC = EXP(-DUM)
C
C                   CALCULATE DENSITY FACTORS
C
      DOUBLE = IANDJ.AND.IG.NE.JG
      MAX = MAXJ
      NN  = 0
C
      DUM1 = ZERO
      DUM2 = ZERO
C
      DO 380 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 380 J = MINJ,MAX
            NN = NN+1
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
            DIJ(NN) = DUM2
  380 CONTINUE
C
C             CALCULATE INTEGRALS FOR SPECIFIC PROPERTIES
C
cjhs      IF(PRNAME.EQ.ELDEN)
cjhs     *  CALL INTDEN(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT)
C
cjhs      IF(PRNAME.EQ.ELFLD)
cjhs     *  CALL INTEFL(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT,AA,AX,AY,AZ)
C
cjhs      IF(PRNAME.EQ.ELMOM)
        CALL INTMOM(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT,AA,AX,AY,AZ)
C
C                  NEXT TWO CLOSE LOOPS OVER GAUSSIAN PRIMITIVES
  400 CONTINUE
  410 CONTINUE
C
C                  STORE THIS INTEGRAL SHELL SET IN OUTPUT MATRIX
C
      MAX = MAXJ
      DO 460  K=1,NVAL
         NL2 = (K-1)*L2
         NN  = (K-1)*IJ
         DO 450  I=MINI,MAXI
            LI = LOCI + I
            IN = LI*(LI-1)/2 + NL2
            IF (IANDJ) MAX = I
            DO 440  J=MINJ,MAX
               LJ = LOCJ + J
               JN = LJ + IN
               NN = NN+1
               XVAL(JN) = WINT(NN)
  440       CONTINUE
  450    CONTINUE
  460 CONTINUE
C
C         NEXT TWO CLOSE SHELL LOOPS
  500 CONTINUE
  510 CONTINUE
C
C        TO COMPLETE PARALLEL HERE, DDI_GSUMF(918,XVAL,NVAL*L2)
C        SEE NOTES AT TOP OF THIS ROUTINE ABOUT WHY THIS IS NOT DONE.
C
      RETURN
      END
C
C*MODULE NEOPRP   *DECK NUC_CHARGE
      SUBROUTINE NUC_CHARGE(QMNUC,QLNUC,NTOTCN,N1)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXAO=8192)
c
      DIMENSION QMNUC(NTOTCN),QLNUC(NTOTCN)
C
      COMMON /FMCOM / X(1)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
      N2 = (N1*N1+N1)/2
      N3 = N1*N1
      NBAS_ATOM =  N1/NTOTCN
C
C     GET MEMORY FOR CALCULATING QM NUCLEAR CHARGE
C
      CALL VALFM(LOADFM)
      ISNUC    = LOADFM + 1
      IADEN    = ISNUC    + N2
      IADEN2D  = IADEN    + N2
      IGPOPAO  = IADEN2D  + N3
      IPOPLOW  = IGPOPAO  + N1
      IVINDIAG = IPOPLOW  + N1
      IVTRAN2  = IVINDIAG + N3
      ISEIGEN  = IVTRAN2  + N3
      ISONEH   = ISEIGEN  + N1
      ITEMPLOW = ISONEH   + N3
      IFINLOW  = ITEMPLOW + N3
      IWRK     = IFINLOW  + N3
      IIWRK    = IWRK     + N1*8
      LAST     = IIWRK    + N1
C
      NEEDQNUC = LAST - LOADFM -1
      CALL GETFM(NEEDQNUC)
C
      CALL DAREAD(IDAF,IODA,X(ISNUC),N2,441,0)
      CALL DAREAD(IDAF,IODA,X(IADEN),N2,445,0)
      CALL CPYTSQ(X(IADEN),X(IADEN2D),N1,1)
      CALL OVLPOP(X(IADEN),X(ISNUC),N2)
      CALL GROSSC(X(IADEN),X(IGPOPAO),IA,N1)
C
      CALL NUC_LOWDIN(X(ISNUC),X(IPOPLOW),X(IADEN2D),
     *                X(IVINDIAG),X(IVTRAN2),X(ISEIGEN),X(ISONEH),
     *                X(ITEMPLOW),X(IFINLOW),X(IWRK),X(IIWRK),
     *                N1,N2)
C
C    ---- TOTAL ATOMIC POPULATIONS ----
C
             L = 0
          DO I = 1,NTOTCN
                SUM1 = ZERO
                SUM2 = ZERO
            DO J = L + 1,L + NBAS_ATOM
                SUM1 = SUM1 + X(IGPOPAO-1+J)
                SUM2 = SUM2 + X(IPOPLOW-1+J)
            ENDDO
          QMNUC(I) = SUM1
          QLNUC(I) = SUM2
          L = L + NBAS_ATOM
          ENDDO
C
      CALL RETFM(NEEDQNUC)
C
      RETURN
      END
C
C*MODULE NEOPRP   *DECK NUCEXPCT
      SUBROUTINE NUCEXPCT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT,PAPER,GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 QNUN,QNN
      CHARACTER*10 PBFLAB
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION ITAG(19)
      DIMENSION CTMP(3,MXATM)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /NUCMON/ QNUN(20),QNN(20),PBFLAB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HELMOMC  /
      DATA ELMOM/8HELMOM   /
      DATA ITAG/4H   X,4H   Y,4H   Z,4H  XX,4H  YY,4H  ZZ,4H  XY,
     *          4H  XZ,4H  YZ,4H XXX,4H YYY,4H ZZZ,4H XXY,4H XXZ,
     *          4H XYY,4H YYZ,4H XZZ,4H YZZ,4H XYZ/
C
      OUT   = NPRINT.EQ.3 .AND. MASWRK
      PAPER = (IEMOUT.GE.0  .AND. NPRINT.NE.-5) .AND. MASWRK
      IF (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME) THEN
        OUT   = .TRUE. .AND. MASWRK
        PAPER = .TRUE. .AND. MASWRK
      END IF
C
C     ----- CALCULATE EXPECTATION VALUE OF THE QM NUCLEAR COORDINATES -----
C
      L1 = NUMNB
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      NVAL = 3
C
C     ----- GROW MEMORY -----
C
      LOADFM = 0
      CALL VALFM(LOADFM)
      IELM   = LOADFM + 1
      IMOS   = IELM   + NVAL*L2
      IEMW   = IMOS   + L3
      IAINDX = IEMW   + NVAL*784
      LAST   = IAINDX + L1
      NEED = LAST - LOADFM -1
C
      IF(OUT) WRITE(IW,900) IELM,IEMW,LAST
      CALL GETFM(NEED)
      IF (EXETYP.EQ.CHECK) THEN
         CALL RETFM(NEED)
         RETURN
      END IF
      IF(PAPER) WRITE(IW,910)
C
C     ----- CALCULATE THE EXPECTATION VALUE OF THE QM NUCLEI COORDINATES
C           RELATIVE TO THE ORIGIN OF THE COORDINATE SYSTEM -----
C
      XPTMP = XP
      YPTMP = YP
      ZPTMP = ZP
      XP = 0.0D+00
      YP = 0.0D+00
      ZP = 0.0D+00
C
C     END IF
C
C     ----- COMPUTE DIPOLE INTEGRALS OVER BASIS FUNCTIONS -----
C
      CALL PRCALC_NEO(ELMOM,X(IELM),X(IEMW),NVAL,L2,GOPARR)
      IF(IEMINT.EQ.0 .OR. .NOT.PAPER) GO TO 240
         IMIN=1
         IMAX = 0
         IF(ABS(IEMINT).EQ.1) IMAX=3
         IF(ABS(IEMINT).EQ.2) IMAX=9
         IF(ABS(IEMINT).EQ.3) IMAX=19
         IF(IEMINT.EQ.-2) IMIN=4
         IF(IEMINT.EQ.-3) IMIN=10
         IF (IMAX .EQ. 0) CALL ABRT
         INDEX = IELM+(IMIN-1)*L2
         DO 230 I=IMIN,IMAX
            IF(MASWRK) WRITE(IW,930) ITAG(I)
            CALL PRTRI(X(INDEX),L1)
            INDEX = INDEX + L2
  230    CONTINUE
  240 CONTINUE
C
C     ----- LOOP OVER THE OCCUPIED NUCLEAR ORBITALS AND COORDINATES -----
C
        INDEX=0
        DO 250 IMONB=1, NNA
          INDEX=IELM
          DO I=1,3
            IAT=IQMNUC(IMONB)
            CALL CALCMO_NEO(IMONB,X(INDEX),X(IMOS),X(IAINDX),L1,VALL)
            INDEX=INDEX+L2
            CTMP(I,IAT)=VALL
          ENDDO
  250 CONTINUE
C
C     ----- PRINT THE EXPECTATION VALUES OF THE NUCLEAR POSITIONS -----
C
      WRITE (IW,940)
      DO 270 IQN = 1,NNA
        IATM = IQMNUC(IQN)
        WRITE (IW,950) QNN(IQN),
     *  CTMP(1,IATM),CTMP(2,IATM),CTMP(3,IATM)
  270 CONTINUE
C
      IF(NQMNUC.GT.1) THEN
        WRITE(IW,960)
      ENDIF
      XP = XPTMP
      YP = YPTMP
      ZP = ZPTMP
      CALL RETFM(NEED)
      RETURN
C
  900 FORMAT(1H ,'MEMORY ASSIGNMENT (EXPECTATION VALUE XYZ):',/1H ,
     *       10X,'IELM   =',I8,4X,'IEMW   =',I8,4X,'LAST   =',I8)
  910 FORMAT(/10X,47('-'),/10X,
     * 'EXPECTATION VALUE OF THE QM NUCLEAR COORDINATES', /10X,47('-'))
  930 FORMAT(/10X,A4,' INTEGRALS'/11X,13(1H-))
  940 FORMAT(/6X,'NUC ',27X,'COORDINATES (BOHR)'/
     *         24X,'X',19X,'Y',19X,'Z')
  950 FORMAT(1X,4X,1X,A8,F17.10,2F20.10/)
  960 FORMAT(10X,55('*'),/
     * 10X,'EXPECTATION VALUE OF EACH PROTON POSITION IS CALCULATED',/
     * 10X,'FROM ITS CORRESPONDING MO AND NOT THE TOTAL NUCLEAR WFN',/
     * 10X,'THE PROTONS MAY NOT BE IN THE SAME ORDER AS LISTED IN  ',/
     * 10X,'THE INPUT PROTON COORDINATE LIST AND MAY INSTEAD BE IN ',/
     * 10X,'ORDER OF LOCALIZED NUCLEAR MOS                         ',/
     * 10X,'IF THE NUCLEAR MOS ARE DELOCALIZED THESE VALUES WILL ',/
     * 10X,'NOT CORRESPOND TO THE AVERAGE POSITION OF EACH PROTON  ',/
     * 10X,55('*')/)
      END
C
C*MODULE NEOPRP   *DECK CALCMO_NEO
      SUBROUTINE CALCMO_NEO(MORB,VEXP,ORBS,IA,L1,VALUE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION ORBS(L1,2),VEXP(*),IA(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00)
C
C       CALCULATE EXPECTATION VALUE OF PROPERTY WITH RESPECT TO A
C       SINGLE MOLECULAR ORBITAL.  CALLING PARAMETERS INCLUDE -
C
C         MORB  : MOLECULAR ORBITAL NUMBER
C         VEXP  : MATRIX OF EXPECTATION VALUE WITH RESPECT TO BASIS
C         ORBS  : MATRIX OF MOLECULAR ORBITALS
C         IA    : VECTOR CONTAINING ARRAY INDICES
C         L1    : DIMENSION OF MATRICES
C         VALUE : EXPECTATION VALUE OF PROPERTY
C
C                  SET UP THE TRIANGULAR INDEX ARRAY
C
      VALUE = ZERO
C
      II = 0
      DO 100  I=1,L1
         IA(I) = II
         II = II + I
  100 CONTINUE
C
C     ----- GET NUCLEAR MOLECULAR ORBITALS -----
C
      L3=L1*L1
      CALL DAREAD(IDAF,IODA,ORBS,L3,444,0)
C
C     ----- CALCULATE EXPECTATION VALUE -----
C
      DO 120  K=1,L1
         DO 120  L=1,K
            KL = IA(K) + L
            VALUE = VALUE + ORBS(K,MORB)*VEXP(KL)*ORBS(L,MORB)
  120 CONTINUE
C
      VALUE = TWO*VALUE
C
      DO 140  K=1,L1
         KL = IA(K) + K
         VALUE = VALUE - ORBS(K,MORB)*VEXP(KL)*ORBS(K,MORB)
  140 CONTINUE
      RETURN
      END
C*MODULE NEOPRP  *DECK PRTRI_NU
      SUBROUTINE PRTRI_NU(D,IA,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION D(*), IA(*)
C
      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
C
C     ----- PRINT SYMMETRIC MATRIX -D- OF DIMENSION -N- -----
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      MM1 = MAX-1
      DO 120 I0=1,N,MAX
         IL = MIN(N,I0+MM1)
         WRITE(IW,9008)
         WRITE(IW,9028) (IA(I),I=I0,IL)
         WRITE(IW,9008)
         IL = -1
         DO 100 I=I0,N
            IL=IL+1
            J0=I0+(I*I-I)/2
            JL=J0+MIN(IL,MM1)
            WRITE(IW,9048) IA(I),(D(J),J=J0,JL)
  100    CONTINUE
  120 CONTINUE
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(6X,10(4X,I4,4X))
 9048 FORMAT(I5,1X,10F12.7)
      END
C*MODULE NEOPRP  *DECK PREVNL_NU
      SUBROUTINE PREVNL_NU(V,E,IA,M,N,LDV)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION V(LDV,M),E(M),IA(M)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PRINT OUT EIGENDATA (VECTORS AND VALUES) -----
C     THE ROWS ARE NOT LABELED WITH BASIS TAGS.
C     -V- IS N X M, WITH TRUE LEADING DIMENSION -LDV-
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      IMAX = 0
  100 IMIN = IMAX+1
      IMAX = IMAX+MAX
      IF (IMAX .GT. M) IMAX = M
      WRITE (IW,9008)
      WRITE (IW,9028) (I,I = IMIN,IMAX)
      WRITE (IW,9008)
      WRITE (IW,9068) (E(I),I = IMIN,IMAX)
      WRITE (IW,9008)
      DO 120 J = 1,N
  120 WRITE (IW,9048) IA(J),(V(J,I),I = IMIN,IMAX)
      IF (IMAX .LT. M) GO TO 100
      END IF
      RETURN
C
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,10X,10F11.6)
 9068 FORMAT(15X,10F11.6)
      END
