C 22 FEB 06 - LVS - VDWSEL ROUTINE MODIFIED TO BE USED FOR EFP SCREENING
C  2 FEB 06 - SK  - INPEF: MASWRK USED
C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - ADD 3D GRID ("CUBE FILE")
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  3 SEP 03 - SK  - INPEF: PRINT ONLY FROM MASTER NODE
C 24 JAN 02 - MWS - WFNDEN: RETURN TRUE SCF DENSITY FOR SCF STAGE OF CI
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 27 FEB 98 - MWS - ENGANL,WFNDEN: ALLOW FOR CITYP NE GUGA
C 20 NOV 97 - MWS - ENGANL: CORRECT MISTAKE OF 1 OCT STORING VTOT
C  1 OCT 97 - MWS - ENGANL: REFORMAT THE OUTPUT LINES
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C  8 AUG 97 - MWS - USE PARAMETER TO SET SIZE OF /POINTS/ COMMON
C 13 JUN 96 - VAG - CHANGES TO INTRODUCE CITYP VARIABLE
C 29 MAR 95 - SK  - ADD PIANL, PIANL1 ENERGY ANALYSIS
C 29 DEC 94 - TLW - ADD INPEF SUBROUTINE FOR EXTERNAL ELECTRIC FIELDS
C 17 NOV 94 - MS  - EXTENDED PDC IMPLEMENTATION
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 22 SEP 93 - MWS - FIX PRINTING OF PDC INPUT VALUES
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C 28 JUN 93 - MS  - CHKLOC,INPPDC,POINT,PROPIN,SPHERE: PDC-MEP CODE
C 11 MAR 93 - FJ  - ENGANL: SKIP PRINT IF EXTENDED VIRIAL ANALYSIS
C  2 APR 92 - MWS,TLW - COMMON ENRGYS MADE PURE FLOATING POINT
C 17 MAR 92 - TLW - INPPGS: PARALLEL I/O CHANGES
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  5 MAR 92 - MWS - CHANGE KINETIC ENERGY INTEGRAL DAF RECORD
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 28 NOV 91 - STE - INPPGS: USE MORE CHARACTER VARIABLES
C 17 OCT 90 - LJ  - GRID: COMPUTE CORNERS CORRECTLY
C  6 OCT 90 - MWS - INPPGS: FIX ERROR IN CONVERTING GRDSIZ TO BOHRS
C 25 JUL 90 - MWS - MOVE AIMPAC INTERFACE TO PARLEY MODULE
C 16 MAY 90 - MWS - AIMPAC: ADD SQRT3 FACTOR FOR DXY,DXZ,DYZ GAUSSIANS
C 20 MAR 90 - MK  - INPPGS: FIX CONVERSION TO BOHR.
C  8 MAR 90 - MWS - INPPGS,GRID: SIMPLIFY $GRID AND $POINTS INPUT
C 23 FEB 90 - MWS - MOVED PROPTY TO PRPPOP (DECSTN QUIRK)
C 23 OCT 89 - MWS - EXETYP=CHECK RUNS DON'T CALL AIMPAC
C 25 SEP 89 - MWS - ADD ROUTINES AIMMEM,AIMPAC
C 24 SEP 89 - MWS - NEW MODULE, CREATED FROM CODE IN PRPPOP AND PRPEL
C
C*MODULE PRPLIB  *DECK CHKLOC
      SUBROUTINE CHKLOC(DNAME,WHERE,IERR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION NUCLEI
C
      PARAMETER (MXPTPT=100)
C
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,MODGRD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PDCPAR/ CENTER(3),DPOLE(3),QPOLE(6),RMAX,DELR,CONSTR,
     *                PTSEL,VDWSCL,PTDENS,VDWINC,NFREQ,LAYER,NPDC
      COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
C
      DATA COMASS  /8HCOMASS  /,
     *     COCHRG  /8HCOCHARGE/,
     *     NUCLEI  /8HNUCLEI  /,
     *     POINTS  /8HPOINTS  /,
     *     GRID    /8HGRID    /,
     *     PDC     /8HPDC     /
C
C       CHECK TO SEE IF LOCATION SPECIFIED IN INPUT GROUP HAS BEEN
C       INCLUDED IN INPUT STREAM.  CALLING PARAMETERS INCLUDE -
C
C            DNAME  : NAME OF INPUT GROUP ($ELDENS, $ELMOM, ...)
C            WHERE  : LOCATION
C            IERR   : =1 IF LOCATION NOT INCLUDED
C
      IERR = 0
C
      IF(WHERE.NE.POINTS) GO TO 200
         IF(NPOINT.EQ.0 .AND. MASWRK) WRITE(IW,900) DNAME,WHERE
         IF(NPOINT.EQ.0) IERR=1
         RETURN
C
  200 CONTINUE
      IF(WHERE.NE.GRID) GO TO 300
         IF(NGRID.EQ.0 .AND. MASWRK) WRITE(IW,900) DNAME,WHERE
         IF(NGRID.EQ.0) IERR=1
         RETURN
C
  300 CONTINUE
      IF(WHERE.NE.PDC) GO TO 400
         IF(NPDC.EQ.0 .AND. MASWRK) WRITE(IW,900) DNAME,WHERE
         IF(NPDC.EQ.0) IERR=1
         RETURN
C
  400 CONTINUE
      IF(WHERE.EQ.COCHRG) RETURN
      IF(WHERE.EQ.COMASS) RETURN
      IF(WHERE.EQ.NUCLEI) RETURN
      IERR=1
      IF(MASWRK) WRITE(IW,910) DNAME,WHERE
      RETURN
  900 FORMAT(1X,'**** ERROR, YOUR $',A8,' REQUIRES A $',A8,' GROUP,'/
     *       1X,'WHICH IS EITHER MISSING OR BOGUS.')
  910 FORMAT(1X,'**** ERROR, $',A8,' - ILLEGAL WHERE=',A8)
      END
C*MODULE PRPLIB  *DECK ENGANL
      SUBROUTINE 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 /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 /VIRIAL/ SCALTE,SCALTT,GVIR,VTCONV,MAXVT,VTSCAL,VIROK,LVCLN
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,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     ----- EVERYTHING ELSE BY SUBTRACTION
C
      E2 = ETOT - E1 - 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
         WRITE(IW,9020) E1,E2,ENUCR,ETOT
         WRITE(IW,9030) VEE,VNE,VNN,VTOT,TKIN,VIRIAL
      END IF
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      RETURN
C
 9000 FORMAT(/10X,17(1H-)/10X,17HENERGY COMPONENTS/10X,17(1H-)/)
 9010 FORMAT( 1X,'        WAVEFUNCTION NORMALIZATION =',F19.10)
 9020 FORMAT(/1X,'               ONE ELECTRON ENERGY =',F19.10/
     *        1X,'               TWO ELECTRON ENERGY =',F19.10/
     *        1X,'          NUCLEAR REPULSION ENERGY =',F19.10/
     *       38X,18(1H-)/
     *        1X,'                      TOTAL ENERGY =',F19.10)
 9030 FORMAT(/1X,'ELECTRON-ELECTRON POTENTIAL ENERGY =',F19.10/
     *        1X,' NUCLEUS-ELECTRON POTENTIAL ENERGY =',F19.10/
     *        1X,'  NUCLEUS-NUCLEUS POTENTIAL ENERGY =',F19.10/
     *       38X,18(1H-)/
     *        1X,'            TOTAL POTENTIAL ENERGY =',F19.10/
     *        1X,'              TOTAL KINETIC ENERGY =',F19.10/
     *        1X,'                VIRIAL RATIO (V/T) =',F19.10)
      END
C*MODULE PRPLIB  *DECK GRID
      SUBROUTINE GRID
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION CORNER(3,4)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,MODGRD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (TM6=1.0D-06)
C
C       DETERMINE PARAMETERS TO GENERATE RECTANGULAR GRID OF POINTS.
C         UX(I),UY(I),UZ(I) : UNIT VECTORS DEFINING PLANE ORIENTATION
C         NXG,NYG           : NUMBER OF INCREMENTS ALONG GRID AXES
C
C       ON INPUT, WE KNOW
C         ORIGIN(I) : COORDINATES OF POINT DEFINING PLANE ORIGIN
C         XVEC(I)   : COORDINATES OF POINT DEFINING PLANE X-AXIS
C         YVEC(I)   : COORDINATES OF POINT DEFINING PLANE +Y DIRECTION
C         GRDSIZ    : GRID SIZE INCREMENT
C
C    DETERMINE VECTORS DEFINING GRID X-AXIS AND POSITIVE-Y DIRECTION
C
      DO 120 I=1,3
         UX(I) = XVEC(I) - ORIGIN(I)
         UY(I) = YVEC(I) - ORIGIN(I)
  120 CONTINUE
      IF(IAND(MODGRD,1).EQ.0) THEN
C
C                 DETERMINE VECTOR DEFINING GRID Z-AXIS
C
         UZ(1) = UX(2)*UY(3) - UX(3)*UY(2)
         UZ(2) = UX(3)*UY(1) - UX(1)*UY(3)
         UZ(3) = UX(1)*UY(2) - UX(2)*UY(1)
C
C                  CHECK FOR COLINEARITY OF DEFINING VECTORS
C
         UMOD  = SQRT(UZ(1)*UZ(1) + UZ(2)*UZ(2) + UZ(3)*UZ(3))
         IF(UMOD .LT. TM6) THEN
            IF (MASWRK) WRITE(IW,*)
     *         '$GRID X AND Y DIRECTIONS ARE COLLINEAR'
            CALL ABRT
         END IF
C
C                    NORMALIZE VECTOR DEFINING GRID Z-AXIS
C
         UZ(1) = UZ(1)/UMOD
         UZ(2) = UZ(2)/UMOD
         UZ(3) = UZ(3)/UMOD
C
C                    DETERMINE VECTOR DEFINING GRID Y-AXIS
C
         UY(1) = UZ(2)*UX(3) - UZ(3)*UX(2)
         UY(2) = UZ(3)*UX(1) - UZ(1)*UX(3)
         UY(3) = UZ(1)*UX(2) - UZ(2)*UX(1)
      ELSE
         DO I=1,3
            UZ(I) = ZVEC(I) - ORIGIN(I)
         ENDDO
         UMOD  = SQRT(UZ(1)*UZ(1) + UZ(2)*UZ(2) + UZ(3)*UZ(3))
         UZ(1) = UZ(1)/UMOD
         UZ(2) = UZ(2)/UMOD
         UZ(3) = UZ(3)/UMOD
      ENDIF
C
C                 NORMALIZE VECTOR DEFINING GRID Y-AXIS
C
      UMOD  = SQRT(UY(1)*UY(1) + UY(2)*UY(2) + UY(3)*UY(3))
      UY(1) = UY(1)/UMOD
      UY(2) = UY(2)/UMOD
      UY(3) = UY(3)/UMOD
C
C                 NORMALIZE VECTOR DEFINING GRID X-AXIS
C
      UMOD  = SQRT(UX(1)*UX(1) + UX(2)*UX(2) + UX(3)*UX(3))
      UX(1) = UX(1)/UMOD
      UX(2) = UX(2)/UMOD
      UX(3) = UX(3)/UMOD
C
C           DETERMINE NUMBER OF INCREMENTS IN X,Y DIRECTIONS
C
      XWID = (XVEC(1)-ORIGIN(1))*UX(1)
     *     + (XVEC(2)-ORIGIN(2))*UX(2)
     *     + (XVEC(3)-ORIGIN(3))*UX(3)
      YWID = (YVEC(1)-ORIGIN(1))*UY(1)
     *     + (YVEC(2)-ORIGIN(2))*UY(2)
     *     + (YVEC(3)-ORIGIN(3))*UY(3)
      ZWID = (ZVEC(1)-ORIGIN(1))*UZ(1)
     *     + (ZVEC(2)-ORIGIN(2))*UZ(2)
     *     + (ZVEC(3)-ORIGIN(3))*UZ(3)
      NXG  = INT(XWID/GRDSIZ) + 1
      NYG  = INT(YWID/GRDSIZ) + 1
      NZG  = INT(ZWID/GRDSIZ) + 1
      IF (MASWRK) THEN
         IF(NFG.EQ.0) THEN
            WRITE(IW,900) NXG,NYG
         ELSE
            WRITE(IW,1000) (UX(I),I=1,3),(UY(I),I=1,3),(UZ(I),I=1,3)
            WRITE(IW,905) NXG,NYG,NZG
         ENDIF
      ENDIF
C
C                   RENORMALIZE UNIT VECTORS TO GRDSIZ
C
      DO 180  I=1,3
         UX(I) = UX(I)*GRDSIZ
         UY(I) = UY(I)*GRDSIZ
         UZ(I) = UZ(I)*GRDSIZ
  180 CONTINUE
C
      DO 200 I=1,3
         CORNER(I,1) = ORIGIN(I)
         CORNER(I,2) = ORIGIN(I) + (NXG-1)*UX(I)
         CORNER(I,3) = ORIGIN(I)                 + (NYG-1)*UY(I)
         CORNER(I,4) = ORIGIN(I) + (NXG-1)*UX(I) + (NYG-1)*UY(I)
  200 CONTINUE
      IF (MASWRK) WRITE(IW,910) ((CORNER(I,J),I=1,3),J=1,4)
      RETURN
C
  900 FORMAT(1X,'THE GRID WILL CONTAIN',I6,' BY',I6,' MESH POINTS.')
  905 FORMAT(1X,'THE 3D GRID WILL CONTAIN',I6,' BY',I6,' BY',I6,
     *          ' MESH POINTS.')
  910 FORMAT(1X,'THE CORNERS OF THE GRID LIE AT (BOHR UNITS)'/
     *       1X,' LOWER LEFT=',3F20.10/
     *       1X,'LOWER RIGHT=',3F20.10/
     *       1X,' UPPER LEFT=',3F20.10/
     *       1X,'UPPER RIGHT=',3F20.10)
 1000 FORMAT(1X,'THE GRID VECTORS ARE:',
     *      /1X,'UX=',3F10.5,/1X,'UY=',3F10.5,/1X,'UZ=',3F10.5)
      END
C*MODULE PRPEL   *DECK INPEF
      SUBROUTINE INPEF(NOSYM,NZVAR,RUNTYP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (NNAM=2, ZERO=0.0D+00)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      LOGICAL SYM,GOPARR,DSKWRK,MASWRK,EFLDL
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA EFIELD  /8HEFIELD  /
      DATA QNAM    /8HEVEC    , 8HSYM     /
      DATA KQNAM   /33,0/
      DATA ENERGY  /8HENERGY  /, PROP /8HPROP    /
C
      JRET = 0
      EVEC(1) = ZERO
      EVEC(2) = ZERO
      EVEC(3) = ZERO
      SYM = .FALSE.
      EFLDL = .FALSE.
C
C                        READ NAMELIST $EFIELD
C
      CALL NAMEIO(IR,JRET,EFIELD,NNAM,QNAM,KQNAM,
     *            EVEC,SYM,
     *            0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF (JRET.EQ.1) RETURN
      IF (JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,910)
         CALL ABRT
      END IF
C
      IF ((EVEC(1).EQ.ZERO).AND.(EVEC(2).EQ.ZERO).AND.
     *    (EVEC(3).EQ.ZERO)) THEN
       IF(MASWRK) WRITE(IW,920)
       CALL ABRT
      END IF
C
      IF(MASWRK) WRITE(IW,930) EVEC(1),EVEC(2),EVEC(3),SYM
C
      IF (.NOT.SYM) THEN
       NOSYM = 1
       IF(MASWRK) WRITE(IW,940)
      END IF
C
C     CHECK THAT ONLY ENERGY AND PROPERTY RUNS USE AN
C     EXTERNAL ELECTRIC FIELD WITH INTERNALS
C
      IF ((NZVAR.NE.0).AND.(RUNTYP.NE.ENERGY.AND.RUNTYP.NE.PROP)) THEN
       IF(MASWRK) WRITE(IW,950) RUNTYP
       CALL ABRT
      END IF
C
      EFLDL = .TRUE.
      RETURN
  910 FORMAT(1H ,'*** $EFIELD - ILLEGAL INPUT, STOP')
  920 FORMAT(1H ,'A ZERO VECTOR HAS BEEN INPUT FOR',
     *           'THE ELECTRIC FIELD.',/
     *           'THE PROGRAM IS TERMINATING.')
  930 FORMAT(/5X,'$EFIELD OPTIONS '/5X,15(1H-)/
     *   5X,'EVEC=',3F8.4,'    SYM=',L8/)
  940 FORMAT(1H ,'SYMMETRY HAS BEEN TURNED OFF BY THE $EFIELD GROUP')
  950 FORMAT(1H ,'INTERNALS DO NOT WORK WITH THE $EFIELD FOR ',
     *          'RUNTYP = ',A8)
      END
C*MODULE PRPLIB  *DECK INPELD
      SUBROUTINE INPELD
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION NUCLEI
C
      PARAMETER (NNAM=5, MXATM=2000)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM),ELDENS(1)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      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
C
      DATA NUCLEI/8HNUCLEI   /
      DATA PUNCH/8HPUNCH   /, BOTH/8HBOTH    /, PAPER/8HPAPER   /
C
      DATA ELDENS  /8HELDENS  /
      DATA QNAM    /8HIEDEN   ,8HMORB    , 8HWHERE   , 8HOUTPUT  ,
     *              8HIEDINT  /
      DATA KQNAM   /1,1,5,5,1/
C
C       READ IN PARAMETERS FOR ELECTRON DENSITY CALCULATION, INCLUDING
C
C            MODENS - 0 > CALCULATE FOR ENTIRE WAVEFUNCTION
C                     K > CALCULATE FOR MOLECULAR ORBITAL K
C            WHERE  - LOCATIONS OF CALCULATIONS
C            OUTPUT - OUTPUT DESTINATION
C
C                        INITIALIZE PARAMETERS
C
      JRET=0
      OUTPUT=BOTH
      IEDEN  = 0
      MORB   = 0
      IEDOUT = 1
      IEDINT = 0
      WHERE  = NUCLEI
C
C                        READ NAMELIST $ELDENS
C
      CALL NAMEIO(IR,JRET,ELDENS,NNAM,QNAM,KQNAM,
     *            IEDEN,MORB,WHERE,OUTPUT,IEDINT,
     *            0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF (JRET .GT. 1) CALL ABRT
      IF(OUTPUT.EQ.PUNCH) IEDOUT=-1
      IF(OUTPUT.EQ.PAPER) IEDOUT=0
      IF(OUTPUT.EQ.BOTH ) IEDOUT=1
      IF((MORB.LT.0  .OR.  MORB.GT.NUM) .AND. MASWRK)
     *   WRITE(IW,900) MORB,NUM
      IF(MORB.LT.0  .OR.  MORB.GT.NUM) CALL ABRT
      MODENS=MORB
      IERR=0
      CALL CHKLOC(ELDENS,WHERE,IERR)
      IF(IERR .NE. 0) CALL ABRT
      ELDLOC=WHERE
      RETURN
C
  900 FORMAT(1H ,'*** $ELDENS - ILLEGAL VALUE, MORB = ',I3,
     *         ' NUM AO=',I4)
      END
C*MODULE PRPLIB  *DECK INPELF
      SUBROUTINE INPELF
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION NUCLEI
C
      PARAMETER (NNAM=4)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM),ELFLDG(1)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA PUNCH   /8HPUNCH   /, PAPER   /8HPAPER   /,
     *     BOTH    /8HBOTH    /, NUCLEI  /8HNUCLEI  /
C
      DATA ELFLDG  /8HELFLDG  /
      DATA QNAM    /8HIEFLD   , 8HWHERE   , 8HOUTPUT  ,8HIEFINT  /
      DATA KQNAM   /1,5,5,1/
C
C    READ IN PARAMETERS FOR ELECTROSTATIC FIELD/GRADIENT CALCULATION
C
C                        PROVIDE DEFAULTS
C
      JRET=0
      OUTPUT=BOTH
      IEFLD  = 0
      IEFOUT = 1
      IEFINT = 0
      WHERE  = NUCLEI
C
C                        READ NAMELIST $ELFLDG
C
      CALL NAMEIO(IR,JRET,ELFLDG,NNAM,QNAM,KQNAM,
     *            IEFLD,WHERE,OUTPUT,IEFINT,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF (JRET .GT. 1) CALL ABRT
      IF((IEFLD.LT.0  .OR.  IEFLD.GT.2) .AND. MASWRK)
     *   WRITE(IW,900) IEFLD
      IF(IEFLD.LT.0  .OR.  IEFLD.GT.2) CALL ABRT
      IERR=0
      CALL CHKLOC(ELFLDG,WHERE,IERR)
      IF(IERR .NE. 0) CALL ABRT
      ELFLOC = WHERE
      IF(OUTPUT.EQ.PUNCH) IEFOUT = -1
      IF(OUTPUT.EQ.PAPER) IEFOUT = 0
      IF(OUTPUT.EQ.BOTH ) IEFOUT = 1
      RETURN
  900 FORMAT(1H ,'*** $ELFLD - ILLEGAL VALUE, IEFLD = ',I4)
      END
C*MODULE PRPLIB  *DECK INPELM
      SUBROUTINE INPELM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (NNAM=4)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM),ELMOM(1)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA COMASS  /8HCOMASS  /,
     *     GRID    /8HGRID    /,
     *     PUNCH   /8HPUNCH   /,
     *     PAPER   /8HPAPER   /,
     *     BOTH    /8HBOTH    /
C
      DATA ELMOM   /8HELMOM   /
      DATA QNAM    /8HIEMOM   , 8HWHERE   , 8HOUTPUT  ,8HIEMINT  /
      DATA KQNAM   /1,5,5,1/
C
C       READ IN PARAMETERS FOR ELECTROSTATIC MOMENT CALCULATION
C
      JRET=0
      OUTPUT=BOTH
      IEMOM  = 1
      IEMOUT = 1
      IEMINT = 0
      WHERE  = COMASS
C
C                         READ NAMELIST $ELMOM
C
      CALL NAMEIO(IR,JRET,ELMOM,NNAM,QNAM,KQNAM,
     *            IEMOM,WHERE,OUTPUT,IEMINT,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,910)
         CALL ABRT
      END IF
C
      IF(IEMOM.GE.0 .AND. IEMOM.LE.3) GO TO 110
         IF (MASWRK) WRITE(IW,900) IEMOM
         CALL ABRT
  110 CONTINUE
      IF(WHERE.EQ.GRID) CALL ABRT
      IF(OUTPUT .EQ. PUNCH) IEMOUT=-1
      IF(OUTPUT .EQ. PAPER) IEMOUT=0
      IF(OUTPUT .EQ. BOTH) IEMOUT=1
      IERR=0
      CALL CHKLOC(ELMOM,WHERE,IERR)
      IF(IERR .NE. 0) CALL ABRT
      ELMLOC = WHERE
      RETURN
C
  900 FORMAT(1X,'ILLEGAL VALUE FOR IEMOM')
  910 FORMAT(1H ,'*** $ELMOM - ILLEGAL INPUT, STOP')
      END
C*MODULE PRPLIB  *DECK INPELP
      SUBROUTINE INPELP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION NUCLEI
C
      PARAMETER (NNAM=4)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM),ELPOT(1)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA NUCLEI  /8HNUCLEI  /, BOTH    /8HBOTH    /,
     *     PAPER   /8HPAPER   /, PUNCH   /8HPUNCH   /
C
      DATA ELPOT    /8HELPOT   /
      DATA QNAM    /8HIEPOT   , 8HWHERE   , 8HOUTPUT  ,8HIEPINT  /
      DATA KQNAM   /1,5,5,1/
C
C       READ IN PARAMETERS FOR ELECTROSTATIC POTENTIAL CALCULATION,
C           PROVIDE DEFAULTS
C
      JRET=0
      OUTPUT=BOTH
      IEPOT  = 0
      IEPOUT = 1
      IEPINT = 0
      WHERE  = NUCLEI
C
C                        READ NAMELIST $ELPOT
C
      CALL NAMEIO(IR,JRET,ELPOT,NNAM,QNAM,KQNAM,
     *            IEPOT,WHERE,OUTPUT,IEPINT,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF (JRET .GT. 1) CALL ABRT
      IF((IEPOT.LT.0   .OR.   IEPOT.GT.1) .AND. MASWRK)
     *   WRITE(IW,900)
      IF(IEPOT.LT.0   .OR.   IEPOT.GT.1) CALL ABRT
      IF(OUTPUT.EQ.PUNCH) IEPOUT= -1
      IF(OUTPUT.EQ.PAPER) IEPOUT=  0
      IF(OUTPUT.EQ.BOTH ) IEPOUT=  1
      IERR=0
      CALL CHKLOC(ELPOT,WHERE,IERR)
      IF(IERR .NE. 0) CALL ABRT
      ELPLOC = WHERE
      RETURN
  900 FORMAT(1X,'**** ERROR, $ELPOT HAS IEPOT=',I4)
      END
C*MODULE PRPLIB  *DECK INPPDC
      SUBROUTINE INPPDC
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MAKEFP
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (NNAM=14)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PDCPAR/ CENTER(3),DPOLE(3),QPOLE(6),RMAX,DELR,CONSTR,
     *                PTSEL,VDWSCL,PTDENS,VDWINC,NFREQ,LAYER,NPDC
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      PARAMETER (TOANGS=0.52917724924D+00, ONE=1.0D+00)
      PARAMETER (DEBYE=ONE/2.541766D+00, BUCK=ONE/1.345044D+00)
C
      DATA ANGS   /8HANGS    /, BOHR   /8HBOHR    /,
     *     CHARGE /8HCHARGE  /, DIPOLE /8HDIPOLE  /,
     *     QUPOLE /8HQUPOLE  /, RNONE  /8HNONE    /,
     *     CHELPG /8HCHELPG  /, GEODES /8HGEODESIC/, CONLLY /8HCONNOLLY/
      DATA MAKEFP/8HMAKEFP  /
C
C                 SET UP INPUT FOR $PDC
C
      DATA PDCWD/8HPDC     /
      DATA QNAM/8HRMAX    ,8HDELR    ,8HCONSTR  ,8HCENTER  ,8HDPOLE   ,
     *          8HQPOLE   ,8HPDUNIT  ,8HVDWSCL  ,8HPTSEL   ,8HNFREQ   ,
     *          8HPTDENS  ,8HVDWINC  ,8HLAYER   ,8HMAXPDC  /
      DATA KQNAM/3,3,5,33,33,63,5,3,5,1,3,3,1,1/
C
C     SPECIFY DEFAULT PARAMETERS
C     DEFAULT PTDENS IS DENSITY PER BOHR**2, WHICH IS 1.0 PTS PER ANG**2
C
      RMAX   = 3.0D+00
      DELR   = 0.8D+00
      CONSTR = CHARGE
      CALL VCLR(CENTER,1,3)
      CALL VCLR(DPOLE,1,3)
      CALL VCLR(QPOLE,1,6)
      PDUNIT = ANGS
      PTSEL  = GEODES
      NFREQ  = 30
      PTDENS = 0.28002830D+00
C       REAL DEFAULT TO BE SET VERY SOON, OR MUCH LATER...
      VDWSCL = 0.0D+00
C
      IF(RUNTYP.NE.MAKEFP) THEN
         VDWINC = 0.2D+00
         LAYER  = 4
         MAXPDC = 10000
      ELSE
         VDWINC = 0.1D+00
         LAYER  = 25
         MAXPDC = 100000
      END IF
C
C     READ IN PARAMETERS TO DETERMINE POINTS AT WHICH PROPERTIES
C     ARE TO BE CALCULATED.  THIS INCLUDES THE INPUT SPECIFIED IN
C          $PDC    - POINTS TO DETERMINE POTENTIAL-DERIVED CHARGES
C
      JRET=0
      CALL NAMEIO(IR,JRET,PDCWD,NNAM,QNAM,KQNAM,
     *            RMAX,DELR,CONSTR,CENTER,DPOLE,QPOLE,PDUNIT,VDWSCL,
     *            PTSEL,NFREQ,PTDENS,VDWINC,LAYER,MAXPDC,
     *    0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF(MASWRK) WRITE(IW,930) PDCWD
         CALL ABRT
      END IF
C
C     SET DEFAULT VDW SCALE FACTOR BASED ON PTSEL CHOICE.
C     MAKEFP JOBS HAVEN'T READ $STONE AT THIS POINT, AND THEREFORE
C     MUST POSTPONE SETTING THE DEFAULT UNTIL CHGPEN.SRC IS RUN.
C
      IF(VDWSCL.EQ.0.0D+00  .AND.  RUNTYP.NE.MAKEFP) THEN
        IF(PTSEL.EQ.CHELPG) VDWSCL=1.00D+00
        IF(PTSEL.EQ.GEODES) VDWSCL=1.40D+00
        IF(PTSEL.EQ.CONLLY) VDWSCL=1.40D+00
      END IF
C
C        IF THE GROUP WAS NOT FOUND, IT IS USUALLY IRRELEVANT TO PRINT
C
      IF(MASWRK  .AND.  (JRET.EQ.0  .OR.  RUNTYP.EQ.MAKEFP)) THEN
         WRITE(IW,950) VDWSCL,RMAX,DELR,VDWINC,PTDENS,LAYER,NFREQ,
     *                 MAXPDC,CONSTR,PDUNIT,PTSEL,CENTER,DPOLE,QPOLE
      END IF
C
      NERR=0
      IF(PDUNIT.NE.BOHR  .AND.  PDUNIT.NE.ANGS) THEN
         IF(MASWRK) WRITE(IW,*) 'ILLEGAL -PDUNIT- STRING GIVEN'
         NERR=NERR+1
      END IF
      IF(NFREQ.GT.99.OR.NFREQ.LT.-99) THEN
         IF(MASWRK) WRITE(IW,*) 'ILLEGAL VALUE FOR -NFREQ- GIVEN'
         NERR=NERR+1
      END IF
      IF(PTSEL.NE.CHELPG.AND.PTSEL.NE.GEODES.AND.PTSEL.NE.CONLLY) THEN
         IF(MASWRK) WRITE(IW,*) 'ILLEGAL -PTSEL- STRING GIVEN'
         NERR=NERR+1
      END IF
      IF(CONSTR.NE.RNONE  .AND.  CONSTR.NE.CHARGE  .AND.
     *   CONSTR.NE.DIPOLE .AND.  CONSTR.NE.QUPOLE) THEN
         IF(MASWRK) WRITE(IW,*) 'ILLEGAL -CONSTR- STRING GIVEN'
         NERR=NERR+1
      END IF
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,*) 'PLEASE FIX THE ABOVE ERRORS IN $PDC'
         CALL ABRT
      END IF
C
C     CONVERT CENTER AND MOMENTS TO ATOMIC UNITS FOR CALCULATIONS
C
      IF(PDUNIT.EQ.ANGS) THEN
         TOBOHR = ONE/TOANGS
         DO 20 I = 1,3
            CENTER(I) = CENTER(I)*TOBOHR
            DPOLE(I)  = DPOLE(I)*DEBYE
            QPOLE(2*I-1)=QPOLE(2*I-1)*BUCK
            QPOLE(2*I)=QPOLE(2*I)*BUCK
   20    CONTINUE
      END IF
C
      NPDC = MAXPDC
      RETURN
C
  930 FORMAT(1X,'**** ERROR IN $',A8,' INPUT')
  950 FORMAT(/5X,'INPUT FOR POTENTIAL-DERIVED CHARGES '/5X,35(1H-)/
     *   5X,'VDWSCL=',F8.2,'   RMAX=',F8.2,'   DELR=',F8.2/
     *   5X,'VDWINC=',F8.2,' PTDENS=',F8.2,/
     *   5X,'LAYER =',I8,  '  NFREQ=',I8,  ' MAXPDC=',I8/
     *   5X,'CONSTR=',A8,  '  UNITS=',A8,'  PTSEL=',A8/
     *   5X,'CENTER=',3F8.4/
     *   5X,'DPOLE =',3F8.4/
     *   5X,'QPOLE =',6F8.4/)
      END
C*MODULE PRPLIB  *DECK INPPGS
      SUBROUTINE INPPGS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (NNAM=7, MXPTPT=100)
C
      CHARACTER*8 WORD
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,MODGRD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
C
      PARAMETER (TOANGS=0.52917724924D+00, ZERO=0.0D+00, ONE=1.0D+00,
     *           QUARTR=2.5D-01)
C
      DATA ANGS    /8HANGS    /, BOHR    /8HBOHR    /
C
C                        SET UP INPUT FOR $GRID
C
      DATA GRIDWD/8HGRID    /
      DATA QNAM/8HORIGIN  ,8HXVEC    ,8HYVEC    ,8HSIZE    ,8HUNITS   ,
     *          8HZVEC    ,8HMODGRD  /
      DATA KQNAM/33,33,33,3,5,33,1/
C
C       READ IN PARAMETERS TO DETERMINE POINTS AT WHICH PROPERTIES
C       ARE TO BE CALCULATED.  THIS INCLUDES THE INPUT SPECIFIED IN
C            $POINTS - LIST OF POINT COORDINATES
C            $GRID   - PARAMETERS TO DETERMINE PLANAR GRID
C
      NPOINT = 0
      NGRID  = 0
      TOBOHR = ONE/TOANGS
C
C    POSITION INPUT TO $POINTS, AND READ THE POINTS
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $POINTS',IEOF)
      IF (IEOF.EQ.1) GO TO 200
      IERR=0
C
      CALL RDCARD('$POINTS ',IEOF)
      IF (IEOF.EQ.1) GO TO 200
      WORD = '        '
      LEN = -8
      CALL GSTRNG(WORD,LEN)
      NPOINT = IFIND('NPOINT  ',IERR)
C
      IF(NPOINT.GE.MXPTPT) THEN
         IF(MASWRK) WRITE(IW,920) MXPTPT
         CALL ABRT
         STOP
      END IF
C
      IPUNIT=0
      IF(WORD.EQ.'BOHR    ') IPUNIT= -1
      IF(WORD.EQ.'ANGS    ') IPUNIT=  1
      IF(IPUNIT.EQ.0) THEN
         IF (MASWRK) WRITE(IW,*) 'ILLEGAL CHOICE OF UNITS ',WORD
         CALL ABRT
      END IF
C
      IF (MASWRK) THEN
      DO 120 I=1,NPOINT
         READ(IR,*) X,Y,Z
         IF(IPUNIT .EQ. 1) THEN
            X = X*TOBOHR
            Y = Y*TOBOHR
            Z = Z*TOBOHR
         END IF
         XPOINT(I) = X
         YPOINT(I) = Y
         ZPOINT(I) = Z
  120 CONTINUE
      END IF
C
C SEND XPOINT, YPOINT, AND ZPOINT TO ALL SLAVES
C
C     IF (GOPARR) THEN
C        CALL DDI_BCAST(327,'F',XPOINT,NPOINT)
C        CALL DDI_BCAST(328,'F',YPOINT,NPOINT)
C        CALL DDI_BCAST(329,'F',ZPOINT,NPOINT)
C     END IF
C
C                      INITIALIZE GRID PARAMETERS
C
  200 CONTINUE
      DO 210 I = 1,3
         ORIGIN(I) = ZERO
         XVEC(I)   = ZERO
         YVEC(I)   = ZERO
         ZVEC(I)   = ZERO
  210 CONTINUE
      IF(NFG.NE.0) THEN
C        SET X,Y,Z UNIT VECTORS AS THE DEFAULT FOR FMO
         XVEC(1)=ONE
         YVEC(2)=ONE
         ZVEC(3)=ONE
      ENDIF
      GRDUNT = ANGS
      GRDSIZ = QUARTR
C
C     MODGRD=0 ORTHONORMALISE THE GRID VECTORS (XVEC,YVEC,ZVEC).
C     MODGRD=1 ONLY NORMALISE THE GRID VECTORS (XVEC,YVEC,ZVEC).
C
      MODGRD=0
C
C                         READ NAMELIST $GRID
C
      JRET=0
      CALL NAMEIO(IR,JRET,GRIDWD,NNAM,QNAM,KQNAM,
     *            ORIGIN,XVEC,YVEC,GRDSIZ,GRDUNT,
     *            ZVEC,MODGRD,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.1) RETURN
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,930) GRIDWD
         CALL ABRT
      END IF
      IF (MASWRK) THEN
        WRITE(IW,950) ORIGIN,XVEC,YVEC
        IF(NFG.NE.0) WRITE(IW,960) ZVEC
        WRITE(IW,970) GRDSIZ,GRDUNT
      ENDIF
C
      IF(GRDUNT.NE.BOHR  .AND.  GRDUNT.NE.ANGS) CALL ABRT
      IF(GRDUNT.EQ.BOHR) IGUNIT= -1
      IF(GRDUNT.EQ.ANGS) IGUNIT=  1
C
C                   CONVERT TO BOHR FOR CALCULATIONS
C
      IF(GRDUNT.EQ.ANGS) THEN
          GRDSIZ = GRDSIZ*TOBOHR
          DO 250 I = 1,3
             ORIGIN(I) = ORIGIN(I)*TOBOHR
             XVEC(I)   = XVEC(I)  *TOBOHR
             YVEC(I)   = YVEC(I)  *TOBOHR
             ZVEC(I)   = ZVEC(I)  *TOBOHR
  250     CONTINUE
      END IF
C
C          SET UP THE REST OF THE GRID VALUES
C
      CALL GRID
      NGRID = NXG*NYG
      RETURN
C
  920 FORMAT(1X,'**** ERROR, MAXIMUM NO. OF POINTS IN $POINTS IS',I5)
  930 FORMAT(1X,'**** ERROR IN $',A8,' INPUT')
  950 FORMAT(/5X,'INPUT DEFINING PLOTTING GRID'/5X,28(1H-)/
     *        1X,'ORIGIN=',3F20.10/
     *        1X,'  XVEC=',3F20.10/
     *        1X,'  YVEC=',3F20.10)
  960 FORMAT( 1X,'  ZVEC=',3F20.10)
  970 FORMAT( 1X,'  SIZE=',F20.10,'   UNITS=',A8)
      END
C*MODULE PRPLIB  *DECK PIANL
      SUBROUTINE PIANL(H,T,FA,FB,DA,DB,LABPI,L1,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION H(L2),T(L2),FA(L2),FB(L2),DA(L2),DB(L2),LABPI(L1)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, SMALL=1.0D-06)
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
      DATA RNONE/8HNONE    /
C
C        ---- SIGMA-PI ENERGY ANALYSIS ----
C        WRITTEN BY SHIRO KOSEKI, SPRING, 1995
C        THIS IS IMPLEMENTED ONLY FOR SPD BASIS FOR RHF,UHF,ROHF
C
      IF(NAT.LT.2) RETURN
      IF((SCFTYP.NE.RHF).AND.(SCFTYP.NE.UHF)
     *                  .AND.(SCFTYP.NE.ROHF)) RETURN
      IF(CITYP.NE.RNONE) RETURN
C
      CALL BASCHK(LMAX)
      IF(LMAX.GE.3) RETURN
C
C     IS THIS MOLECULE PLANAR OR LINEAR?
C
      XSUM = ZERO
      YSUM = ZERO
      ZSUM = ZERO
      DO 100 I=1,NAT
         XSUM = XSUM + ABS(C(1,I))
         YSUM = YSUM + ABS(C(2,I))
         ZSUM = ZSUM + ABS(C(3,I))
  100 CONTINUE
      KPIX = 0
      KPIY = 0
      KPIZ = 0
      IF(XSUM.LT.SMALL) KPIX = 1
      IF(YSUM.LT.SMALL) KPIY = 1
      IF(ZSUM.LT.SMALL) KPIZ = 1
      KPI = KPIX +KPIY +KPIZ
      IF(KPI.LT.1) RETURN
C
      IF(MASWRK) WRITE(IW,9000)
C
C     ASSIGN ATOMIC BASIS FUNCTIONS TO PI SPACE:
C
      DO 140 I=1,NUM
         LABPI(I)=0
  140 CONTINUE
      NAOTMP = 0
      DO 160 II = 1,NSHELL
         I = KATOM(II)
         MINI = KMIN(II)
         MAXI = KMAX(II)
C
         DO 170 I = MINI,MAXI
            NAOTMP = NAOTMP +1
C
C           S ORBITAL: 1
C           IF(I.EQ.1) LABPI(NAOTMP)=0
C
C           P ORBITAL: 2-4 (X,Y,Z)
            IF((I.EQ.2).AND.(KPIX.NE.0)) LABPI(NAOTMP)=1
            IF((I.EQ.3).AND.(KPIY.NE.0)) LABPI(NAOTMP)=1
            IF((I.EQ.4).AND.(KPIZ.NE.0)) LABPI(NAOTMP)=1
C
C           D ORBITAL: 5-10 (XX,YY,ZZ,XY,XZ,YZ)
C           IF(I.EQ.5) LABPI(NAOTMP)=0
C           IF(I.EQ.6) LABPI(NAOTMP)=0
C           IF(I.EQ.7) LABPI(NAOTMP)=0
            IF(I.EQ.8) THEN
               IF((KPIX.EQ.1).AND.(KPIY.EQ.0)) LABPI(NAOTMP)=1
               IF((KPIX.EQ.0).AND.(KPIY.EQ.1)) LABPI(NAOTMP)=1
               END IF
            IF(I.EQ.9) THEN
               IF((KPIX.EQ.1).AND.(KPIZ.EQ.0)) LABPI(NAOTMP)=1
               IF((KPIX.EQ.0).AND.(KPIZ.EQ.1)) LABPI(NAOTMP)=1
               END IF
            IF(I.EQ.10) THEN
               IF((KPIY.EQ.1).AND.(KPIZ.EQ.0)) LABPI(NAOTMP)=1
               IF((KPIY.EQ.0).AND.(KPIZ.EQ.1)) LABPI(NAOTMP)=1
               END IF
C
C           F ORBITAL: 11-20
C           G ORBITAL: 21-35
C
  170    CONTINUE
  160 CONTINUE
C
C     ----- READ BARE NUCLEUS HAMILTONIAN -----
C     ----- READ KINETIC ENERGY INTEGRALS -----
C     ----- READ FOCK MATRIX              -----
C     ----- READ DENSITY MATRIX           -----
C
      CALL DAREAD(IDAF,IODA,H ,L2,11,0)
      CALL DAREAD(IDAF,IODA,T ,L2,13,0)
      CALL DAREAD(IDAF,IODA,FA,L2,14,0)
      CALL DAREAD(IDAF,IODA,DA,L2,16,0)
      IF(SCFTYP.NE.RHF) THEN
         CALL DAREAD(IDAF,IODA,FB,L2,18,0)
         CALL DAREAD(IDAF,IODA,DB,L2,20,0)
      END IF
C
C     ----- CALCULATE THE ENERGY -----
C
      EHF1 = TRACEP(DA,FA,L1)
      EHF2 = TRACEP(DA,H ,L1)
      EHFK = TRACEP(DA,T ,L1)
C
      IF(SCFTYP.NE.RHF) THEN
         EHF1 = EHF1 + TRACEP(DB,FB,L1)
         EHF2 = EHF2 + TRACEP(DB,H ,L1)
         EHFK = EHFK + TRACEP(DB,T ,L1)
      END IF
C
      EHF = (EHF1+EHF2)/TWO
      ETOT  = EHF+EN
      IF(MASWRK) WRITE(IW,9010) EHF1,EHF2,EHF,EHFK,EN,ETOT
C
C     ----- DEVIDE THE ENERGY INTO SIGMA AND PI PARTS -----
C
C                         FOCK MATRIX...
C
      ASIGMA = ZERO
      APIPI  = ZERO
      AMIX   = ZERO
                        CALL PIANL1(DA,FA,LABPI,L1,ASIGMA,APIPI,AMIX,IW)
      IF(SCFTYP.NE.RHF) CALL PIANL1(DB,FB,LABPI,L1,ASIGMA,APIPI,AMIX,IW)
C
C                         BARE HAMITONIAN...
C
      BSIGMA = ZERO
      BPIPI  = ZERO
      BMIX   = ZERO
                        CALL PIANL1(DA,H,LABPI,L1,BSIGMA,BPIPI,BMIX,IW)
      IF(SCFTYP.NE.RHF) CALL PIANL1(DB,H,LABPI,L1,BSIGMA,BPIPI,BMIX,IW)
C
C                         KINETIC ENERGY PART...
C
      CSIGMA = ZERO
      CPIPI  = ZERO
      CMIX   = ZERO
                        CALL PIANL1(DA,T,LABPI,L1,CSIGMA,CPIPI,CMIX,IW)
      IF(SCFTYP.NE.RHF) CALL PIANL1(DB,T,LABPI,L1,CSIGMA,CPIPI,CMIX,IW)
C
C     TOTAL (1+2)
C
      ESIGMA = (ASIGMA+BSIGMA)/TWO
      EPIPI  = (APIPI +BPIPI )/TWO
      EMIX   = (AMIX  +BMIX  )/TWO
C
C     2 ELECTRON PART: (FOCK) - (BARE-H)
C
      FSIGMA = (ASIGMA-BSIGMA)/TWO
      FPIPI  = (APIPI -BPIPI )/TWO
      FMIX   = (AMIX  -BMIX  )/TWO
C
C     NUCLEAR-ELECTRON ATTRACTION: (BARE-H) - (KINETIC)
C
      GSIGMA = BSIGMA-CSIGMA
      GPIPI  = BPIPI -CPIPI
      GMIX   = BMIX  -CMIX
C
      ERROR = EHF -ESIGMA -EPIPI-EMIX
      SKELT = EN + ESIGMA
C
      IF(MASWRK) THEN
         WRITE(IW,9020) ESIGMA,CSIGMA,GSIGMA,FSIGMA,
     *                  EPIPI, CPIPI, GPIPI, FPIPI,
     *                  SKELT, ERROR,
     *                  EMIX,  CMIX,  GMIX,  FMIX
         WRITE(IW,FMT='('' ...... END OF PI ENERGY ANALYSIS ......'')')
      END IF
      RETURN
C
 9000 FORMAT(/1X,' ...... PI ENERGY ANALYSIS ......')
 9010 FORMAT(/1X,'ENERGY ANALYSIS:',
     *       /1X,'           FOCK ENERGY=',F18.10,
     *       /1X,'         BARE H ENERGY=',F18.10,
     *       /1X,'    ELECTRONIC ENERGY =',F18.10,
     *       /1X,'        KINETIC ENERGY=',F18.10,
     *       /1X,'         N-N REPULSION=',F18.10,
     *       /1X,'          TOTAL ENERGY=',F18.10)
 9020 FORMAT( 1X,'       SIGMA PART(1+2)=', F18.10,
     *       /1X,'              (K,V1,2)=',3F18.10,
     *       /1X,'          PI PART(1+2)=', F18.10,
     *       /1X,'              (K,V1,2)=',3F18.10,
     *       /1X,' SIGMA SKELETON, ERROR=',2F18.10,
     *       /1X,'            MIXED PART=',1P,4E12.5)
      END
C*MODULE PRPLIB  *DECK PIANL1
      SUBROUTINE PIANL1(A,B,LABPI,N,ESG,EPI,EMIX,IW)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),B(*), LABPI(*)
C
      N2 = (N*N+N)/2
      IJ = 0
      DO 10 I=1,N
         DO 20 J=1,I
            IJ = IJ+1
            IF(IJ.GT.N2) THEN
               WRITE(IW,200) IJ
               CALL ABRT
               END IF
C
            ABIJ = A(IJ)*B(IJ)
            IF(J.EQ.I) THEN
               IF(LABPI(I).EQ.1) THEN
                  EPI = EPI + ABIJ
                  ELSE
                  ESG = ESG + ABIJ
                  END IF
               ELSE
               IF((LABPI(I).EQ.1).AND.(LABPI(J).EQ.1)) THEN
                   EPI = EPI + ABIJ + ABIJ
                   GO TO 20
                   END IF
               IF((LABPI(I).EQ.0).AND.(LABPI(J).EQ.0)) THEN
                   ESG = ESG + ABIJ + ABIJ
                   GO TO 20
                   END IF
               EMIX = EMIX + ABIJ + ABIJ
               END IF
   20       CONTINUE
   10    CONTINUE
      RETURN
C
  200 FORMAT(/1X,'ERROR IN PIANL1;  IJ=',I10)
      END
C*MODULE PRPLIB  *DECK POINT
      SUBROUTINE POINT(PRPLOC,IPOINT,X,Y,Z,ISTAT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION NUCLEI
C
      PARAMETER (MXATM=2000, MXPTPT=100)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,MODGRD
C
      DATA COMASS  /8HCOMASS  /,
     *     COCHRG  /8HCOCHARGE/,
     *     NUCLEI  /8HNUCLEI  /,
     *     POINTS  /8HPOINTS  /,
     *     GRID    /8HGRID    /,
     *     PDC     /8HPDC     /
C
C     ----- RETURN THE COORDINATES OF THE NEXT POINT
C           WHOSE PROPERTIES ARE TO BE COMPUTED -----
C
C     PRPLOC - TYPE OF LOCATION, CHOOSE FROM
C              COCHRG, COMASS, NUCLEI, POINT, GRID, PDC
C     IPOINT - NUMBER OF POINT
C     X,Y,Z  - COORDINATE OF POINT
C     ISTAT  - STATUS : -1 = ERROR
C                        0 = CONTINUE
C                        1 = LAST POINT IN SEQUENCE ALREADY READ
C
      ISTAT = 0
C
      IF(PRPLOC.EQ.COCHRG) THEN
         IF(IPOINT.LE.1) CALL CALCOC(X,Y,Z)
         IF(IPOINT.GT.1) ISTAT=1
         RETURN
      ENDIF
C
      IF(PRPLOC.NE.COMASS) GO TO 100
         IF(IPOINT.LE.1) CALL CALCOM(X,Y,Z)
         IF(IPOINT.GT.1) ISTAT=1
         RETURN
C
  100 CONTINUE
      IF(PRPLOC.NE.NUCLEI) GO TO 200
         IF(IPOINT.GT.NAT) ISTAT=1
         IF(IPOINT.GT.NAT) RETURN
         X = C(1,IPOINT)
         Y = C(2,IPOINT)
         Z = C(3,IPOINT)
         RETURN
C
  200 CONTINUE
      IF(PRPLOC.NE.POINTS) GO TO 300
         IF(IPOINT.GT.NPOINT) ISTAT=1
         IF(IPOINT.GT.NPOINT) RETURN
         X = XPOINT(IPOINT)
         Y = YPOINT(IPOINT)
         Z = ZPOINT(IPOINT)
         RETURN
C
  300 CONTINUE
      IF(PRPLOC.NE.GRID) GO TO 400
         IF(IPOINT.GT.NGRID) ISTAT=1
         IF(IPOINT.GT.NGRID) RETURN
         IX = (IPOINT - 1)/NYG
         IY = IPOINT - IX*NYG - 1
         X  = ORIGIN(1) + IX*UX(1) + IY*UY(1)
         Y  = ORIGIN(2) + IX*UX(2) + IY*UY(2)
         Z  = ORIGIN(3) + IX*UX(3) + IY*UY(3)
         RETURN
C
  400 CONTINUE
      IF(PRPLOC.NE.PDC) GO TO 500
         CALL ABRT
         STOP
C
  500 CONTINUE
      CALL ABRT
      END
C*MODULE PRPLIB  *DECK PROPIN
      SUBROUTINE PROPIN
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION OUT(3)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA OUT/8HPUNCH   ,8HPAPER   ,8HBOTH    /
C
C     READ INPUT DEFINING POINTS, GRIDS, OR SURFACES.
C
      CALL INPPGS
      CALL INPPDC
C
C     READ INPUT FOR ELECTRON DENSITY, ELECTROSTATIC FIELD/GRADIENT,
C     ELECTOSTATIC MOMENTS, ELECTROSTATIC POTENTIAL
C
      CALL INPELD
      CALL INPELF
      CALL INPELM
      CALL INPELP
      IF (MASWRK)
     *WRITE(IW,900) IEMOM,IEFLD,IEPOT,IEDEN,
     *              ELMLOC,ELFLOC,ELPLOC,ELDLOC,
     *              OUT(IEMOUT+2),OUT(IEFOUT+2),OUT(IEPOUT+2),
     *              OUT(IEDOUT+2),IEMINT,IEFINT,IEDINT,MODENS
      RETURN
C
  900 FORMAT(/10X,16(1H-)/10X,'PROPERTIES INPUT'/10X,16(1H-)//
     *   1X,4X,'MOMENTS',12X,'FIELD',11X,'POTENTIAL',10X,'DENSITY'/
     *   1X,'IEMOM =',I8,3X,'IEFLD =',I8,3X,'IEPOT =',I8,3X,'IEDEN =',
     *   I8/1X,4('WHERE =',A8,3X)/1X,4('OUTPUT=',A8,3X)/
     *   1X,'IEMINT=',I8,3X,'IEFINT=',I8,21X,'IEDINT=',I8/
     *   55X,'MORB  =',I8)
      END
C*MODULE PRPLIB  *DECK SPHCNL
      SUBROUTINE SPHCNL(N,U)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION U(3,N)
C
      DATA TWO/2.0D+00/, EPS/1.0D-10/, PI/3.1415926535898D+00/
C
C     GENERATE UNIT VECTORS OVER SPHERE FOR "CONNOLLY" SURFACE
C     THIS IS THE OLD GENUN ROUTINE FROM MICHAEL CONNOLLY'S MS
C
      NEQUAT = INT(SQRT(PI*N))
      NVERT = NEQUAT/2
      NU = 0
      DO 200 I=1,NVERT+1
        FI = (PI*(I-1))/NVERT
        Z = COS(FI)
        XY = SIN(FI)
        NHOR = INT(NEQUAT*XY+EPS)
        IF(NHOR.LT.1) NHOR = 1
        DO 100 J=1,NHOR
          FJ = (TWO*PI*(J-1))/NHOR
          X = COS(FJ)*XY
          Y = SIN(FJ)*XY
          IF(NU.GE.N) GO TO 300
          NU = NU+1
          U(1,NU) = X
          U(2,NU) = Y
  100     U(3,NU) = Z
  200   CONTINUE
  300 N = NU
      RETURN
      END
C*MODULE PRPLIB  *DECK SPHGEO
      SUBROUTINE SPHGEO(ITYPE,NH,NK,C,NPOINTS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION CIC(3,12),C(3,10000)
      DIMENSION IPS(20),JPS(20),KPS(20)
      DIMENSION BARA(3,406),BARB(3,406),BAR(3,406)
C
      DATA CIC/
     *  0.0D+00, 1.0D+00, 1.61803398875D+00,
     *  0.0D+00,-1.0D+00, 1.61803398875D+00,
     *  0.0D+00, 1.0D+00,-1.61803398875D+00,
     *  0.0D+00,-1.0D+00,-1.61803398875D+00,
     *  1.61803398875D+00,0.0D+00, 1.0D+00,
     *  1.61803398875D+00,0.0D+00,-1.0D+00,
     * -1.61803398875D+00,0.0D+00, 1.0D+00,
     * -1.61803398875D+00,0.0D+00,-1.0D+00,
     *  1.0D+00, 1.61803398875D+00,0.0D+00,
     *  1.0D+00,-1.61803398875D+00,0.0D+00,
     * -1.0D+00, 1.61803398875D+00,0.0D+00,
     * -1.0D+00,-1.61803398875D+00,0.0D+00/
      DATA IPS/5,5,9,5,6,6,6,6,11,9,3,3,11,11,7,12,4,12,12,2/
      DATA JPS/1,2,1,6,5,10,4,3,1,3,8,4,7,8,8,8,10,2,10,1/
      DATA KPS/2,10,5,9,10,4,3,9,9,11,11,8,1,7,12,4,12,7,2,7/
C
C     GENERATES SETS OF COORDINATES ON SURFACE OF SPHERE
C     STARTING FROM 12 POINTS OF THE ICOSAHEDRON AND TESSELATING
C     EACH FACE ACCORDING TO THE NH, NK INDICES
C     THE RESULT IS A TRIANGULAR {3,5+}H,K TESSELATION FOR ITYPE = 1
C     AND A HEXAGONAL {5+,3}H,K TESSELATION FOR ITYPE = 2
C
      IF(ITYPE.EQ.1) THEN
        CALL BARGEN(NH,NK,BAR,NBAR)
      END IF
      IF(ITYPE.EQ.2) THEN
        IF(NH.GT.NK) THEN
          CALL BARGEN(NH-NK,NH+2*NK,BARA,NBARA)
          CALL BARGEN(NH,NK,BARB,NBARB)
        ELSE
          CALL BARGEN(NK+2*NH,NK-NH,BARA,NBARA)
          CALL BARGEN(NH,NK,BARB,NBARB)
        END IF
        CALL BARDIF(BARA,NBARA,BARB,NBARB,BAR,NBAR)
      END IF
C
C     NOW DECORATE EACH OF THE 20 FACES OF THE ICOSAHEDRON
C
      NPT=0
      DO 20 NFACE=1,20
        XA=CIC(1,IPS(NFACE))
        YA=CIC(2,IPS(NFACE))
        ZA=CIC(3,IPS(NFACE))
        XB=CIC(1,JPS(NFACE))
        YB=CIC(2,JPS(NFACE))
        ZB=CIC(3,JPS(NFACE))
        XC=CIC(1,KPS(NFACE))
        YC=CIC(2,KPS(NFACE))
        ZC=CIC(3,KPS(NFACE))
C
C     AND FOR EACH FACE LOOP OVER THE BARYCENTRIC COORDS
C
        DO 30 IB=1,NBAR
          NPT=NPT+1
          C(1,NPT)=BAR(3,IB)*XA+BAR(1,IB)*XB+BAR(2,IB)*XC
          C(2,NPT)=BAR(3,IB)*YA+BAR(1,IB)*YB+BAR(2,IB)*YC
          C(3,NPT)=BAR(3,IB)*ZA+BAR(1,IB)*ZB+BAR(2,IB)*ZC
C
C     IF THIS POINT IS ON AN EDGE OR CORNER THEN
C     CHECK TO SEE IF THIS IS THE SAME AS A PREVIOUS POINT
C
          IF(BAR(1,IB)*BAR(2,IB)*BAR(3,IB).LT.1.0D-08) THEN
            LOOP=NPT-1
            DO 35 K=1,LOOP
              IF(ABS(C(1,K)-C(1,LOOP+1)).GT.1.0D-05) GO TO 35
              IF(ABS(C(2,K)-C(2,LOOP+1)).GT.1.0D-05) GO TO 35
              IF(ABS(C(3,K)-C(3,LOOP+1)).GT.1.0D-05) GO TO 35
              NPT=NPT-1
   35       CONTINUE
          END IF
   30   CONTINUE
   20 CONTINUE
C
C     NORMALISE ALL POINTS TO UNIT SPHERE
C
      NPOINTS=NPT
      DO 50 I=1,NPT
        R=SQRT(C(1,I)*C(1,I)+C(2,I)*C(2,I)+C(3,I)*C(3,I))
        C(1,I)=C(1,I)/R
        C(2,I)=C(2,I)/R
   50   C(3,I)=C(3,I)/R
      RETURN
      END
C*MODULE PRPLIB  *DECK BARDIF
      SUBROUTINE BARDIF(A,NPTA,B,NPTB,C,NPTC)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(3,406),B(3,406),C(3,406)
      DATA DELTA/1.0D-08/
C
      NPTC=0
      DO 100 IA=1,NPTA
        DO 90 IB=1,NPTB
C
C     IF THE IA'TH POINT ISN'T A DUPLICATE, KEEP IT
C
          IF(ABS(A(1,IA)-B(1,IB)).LT.DELTA.AND.
     $       ABS(A(2,IA)-B(2,IB)).LT.DELTA.AND.
     $       ABS(A(3,IA)-B(3,IB)).LT.DELTA) GO TO 100
   90   CONTINUE
        NPTC=NPTC+1
        C(1,NPTC)=A(1,IA)
        C(2,NPTC)=A(2,IA)
        C(3,NPTC)=A(3,IA)
  100 CONTINUE
      RETURN
      END
C*MODULE PRPLIB  *DECK BARGEN
      SUBROUTINE BARGEN(NH,NK,BAR,NBAR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IACC
C
      DIMENSION BAR(3,406)
      DATA HALF/0.5D+00/, TWO/2.0D+00/, RT3/1.73205080756888D+00/
C
C     CREATES BARYCENTRIC COORDINATES FOR A SINGLE FACE
C     LOCATED IN XY-PLANE FOR CONVENIENCE
C     USED BY SPHGEO IN CONJUNCTION WITH BARYCO
C
      RH=NH
      RK=NK
      TT=TWO*(RH*RH+RH*RK+RK*RK)
C
      NBAR=0
      V1X=(RH-RK)/TT
      V1Y=-RT3*(RH+RK)/TT
      V2X=(TWO*RH+RK)/TT
      V2Y=-RT3*RK/TT
      V3X=-(RH+TWO*RK)/TT
      V3Y=-RT3*RH/TT
      XO=HALF
      YO=RT3*HALF
      DO 100 I1=1,NH+NK+1
      CALL BARYCO(XO,YO,NBAR,BAR,IACC)
      IF(.NOT.IACC) GO TO 101
      XP=XO
      YP=YO
      DO 80 I2=1,NK+1
      XP=XP+V2X
      YP=YP+V2Y
      CALL BARYCO(XP,YP,NBAR,BAR,IACC)
      IF(.NOT.IACC) GO TO 81
   80 CONTINUE
   81 XQ=XO
      YQ=YO
      DO 90 I3=1,NH+1
      XQ=XQ+V3X
      YQ=YQ+V3Y
      CALL BARYCO(XQ,YQ,NBAR,BAR,IACC)
      IF(.NOT.IACC) GO TO 91
   90 CONTINUE
   91 XO=XO+V1X
      YO=YO+V1Y
  100 CONTINUE
  101 CONTINUE
      RETURN
      END
C*MODULE PRPLIB  *DECK BARYCO
      SUBROUTINE BARYCO(XR,YR,NBAR,BAR,IACC)
C
C     TESTS WHETHER (XR,YR) IS INSIDE TRIANGLE; IF SO
C     INCREMENTS NBAR AND RETURNS BARYCENTRIC COORDS IN BAR
C     IACC FLAGS WHETHER THE POINT WAS ACCEPTED (T) OR NOT (F)
C     USED BY SPHGEO IN CONJUNCTION WITH BARGEN
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL IACC
C
      DIMENSION BAR(3,406)
      DATA HALF/0.5D+00/,ONE/1.0D+00/,TWO/2.0D+00/,THREE/3.0D+00/,
     *     FOUR/4.0D+00/, ZERO/0.0D+00/,DELTA/1.0D-06/
      DATA RT3/1.73205080756888D+00/,ABC/0.4330127018922192D+00/
C
C     TEST TO SEE IF (XR,YR) LIES OUTSIDE THE EQUILATERAL TRIANGLE
C
      IF(YR.GT.(RT3*XR+DELTA).OR.YR.LT.-DELTA.OR.YR.GT.RT3*(ONE-XR)
     $+DELTA) THEN
      IACC=.FALSE.
      ELSE
C
C     IF INSIDE THEN DETERMINE BARYCENTRIC COORDINATES
C     AR, BR AND CR ARE LENGTHS FROM (XR,YR) TO VERTICES
C     ARC, BRA AND CRB ARE AREAS OF TRIANGLES INSIDE THE
C     EQUILATERAL TRIANGULAR FACE
C
      IACC=.TRUE.
      AR=SQRT(XR*XR+YR*YR)
      BR=SQRT(MAX(ZERO,XR*XR-TWO*XR+ONE+YR*YR))
      CR=SQRT(MAX(ZERO,XR*XR-XR+HALF*HALF+YR*YR-RT3*YR+THREE/FOUR))
      S=HALF*(ONE+AR+CR)
      ARC=SQRT(MAX(ZERO,S*(S-ONE)*(S-AR)*(S-CR)))
      S=HALF*(ONE+BR+AR)
      BRA=SQRT(MAX(ZERO,S*(S-ONE)*(S-BR)*(S-AR)))
      S=HALF*(ONE+CR+BR)
      CRB=SQRT(MAX(ZERO,S*(S-ONE)*(S-CR)*(S-BR)))
      NBAR=NBAR+1
      BAR(1,NBAR)=ARC/ABC
      BAR(2,NBAR)=BRA/ABC
      BAR(3,NBAR)=CRB/ABC
      END IF
      RETURN
      END
C*MODULE PRPLIB  *DECK WFNDEN
      SUBROUTINE WFNDEN(DA,DB,L2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DA(L2),DB(L2)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PSILVL/ IPSI
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA RHF/8HRHF     /, RNONE/8HNONE    /, RMC/8HMCSCF   /
C
C     ----- OBTAIN 1-PARTICLE DENSITY MATRIX FROM DAF -----
C     -WFNDED- DIFFERS FROM -DENDD1- IN THAT IT DOES NOT ADD DA+DB.
C
C            FOR RHF, MCSCF, AND CI, DA IS THE TOTAL DENSITY.
C            FOR UHF AND ROHF, DA AND DB ARE ALPHA AND BETA DENSITY.
C            FOR GVB, DA AND DB ARE CORRECT ONLY FOR HIGH SPIN CASE,
C                     AND OPEN SHELL SINGLET.
C                     DTOT=DA+DB IS CORRECT, BUT DA-DB IS NOT.
C
      CALL DAREAD(IDAF,IODA,DA,L2,16,0)
      IF(SCFTYP.EQ.RHF) RETURN
      IF(SCFTYP.EQ.RMC) RETURN
      IF(CITYP.NE.RNONE  .AND.  IPSI.EQ.1) RETURN
      CALL DAREAD(IDAF,IODA,DB,L2,20,0)
      RETURN
      END
C
C*MODULE PRPLIB  *DECK VDWSEL
      SUBROUTINE VDWSEL(RADIUS,NATM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     LOAD SELECTED VAN DER WAALS RADII DPENDING UPON PTSEL:
C     CONLLY ... MERZ AND KOLLMAN VALUES FOR CONNOLLY SURFACE
C     GEODES ... RECOMMENDED BY GAVEZZOTTI + SPACKMAN
C     CHELPG ... RECOMMENDED BY BRENEMAN AND WIBERG FOR CHELPG.
C
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      CHARACTER*8 DUMMY
      CHARACTER*6 CHTMP
      DOUBLE PRECISION MAKEFP
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      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 /PDCPAR/ CENTER(3),DPOLE(3),QPOLE(6),RMAX,DELR,CONSTR,
     *                PTSEL,VDWSCL,PTDENS,VDWINC,NFREQ,LAYER,NPDC
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
C
      DIMENSION RADIUS(NATM),VANDER(36,3)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, 
     *          TOANGS=0.52917724924D+00)
C
      DATA CHELPG /8HCHELPG  /, GEODES /8HGEODESIC/, CONLLY /8HCONNOLLY/
      DATA MAKEFP/8HMAKEFP  /
C
C     VDW RADII FROM MERZ/KOLLMAN/SINGH
C
      DATA (VANDER(I,1),I=1,36)/1.20D+00,1.20D+00,1.37D+00,1.45D+00,
     *     1.45D+00,1.50D+00,1.50D+00,1.40D+00,1.35D+00,1.30D+00,
     *     1.57D+00,1.36D+00,1.24D+00,1.17D+00,1.80D+00,1.75D+00,
     *     1.70D+00,19*0.0D+00/
C
C     VDW RADII FROM GAVEZZOTTI (J.AM.CHEM.SOC. 105, 5220 (1983))
C     SUPPLEMENTED BY VALUES FOR BORON AND THIRD AND FOURTH ROWS
C     FROM M.A.SPACKMAN (J.CHEM.PHYS. 85, 6579 (1986))
C
      DATA (VANDER(I,2),I=1,36)/1.20D+00,0.00D+00,0.00D+00,0.00D+00,
     *     1.85D+00,1.50D+00,1.50D+00,1.40D+00,1.35D+00,0.00D+00,
     *     0.00D+00,0.00D+00,2.07D+00,2.05D+00,1.96D+00,1.89D+00,
     *     1.80D+00,19*0.0D+00/
C
C     VDW RADII FROM BRENEMAN & WIBERG
C
      DATA (VANDER(I,3),I=1,36)/1.45D+00,1.45D+00,1.50D+00,1.50D+00,
     *     1.50D+00,1.50D+00,1.70D+00,1.70D+00,1.70D+00,1.70D+00,
     *     2.00D+00,2.00D+00,2.00D+00,2.00D+00,2.00D+00,2.00D+00,
     *     2.00D+00,19*0.0D+00/
C
C     SELECT VDW RADII AND CONVERT TO BOHR
C
      TOBOHR = ONE/TOANGS
      IF(PTSEL.EQ.GEODES) THEN
        IVDW=2
        WRITE(IW,1010)
      END IF
      IF(PTSEL.EQ.CHELPG) THEN
         IF(RUNTYP.EQ.MAKEFP) THEN
           IVDW=2
           WRITE(IW,1010)
         ELSE
           IVDW=3
           WRITE(IW,1020)
         END IF
      END IF
      IF(PTSEL.EQ.CONLLY) THEN
        IVDW=1
        WRITE(IW,1000)
      END IF
C
C        FOR MAKEFP SCREENING
      IF(RUNTYP.EQ.MAKEFP) THEN
         IF(MASWRK) WRITE(IW,915) 
      DO 29 I=1,NPTSTN
         ZNUC = ZAN(I) + IZCORE(I)
         IZ = INT(ZNUC)
         RADIUS(I)=ZERO
         IF(IZ.LT.0) CALL ABRT
C           BOND MID-POINT         
         IF(IZ.EQ.0) THEN
            DUMMY='        '
            WRITE(UNIT=DUMMY,FMT='(A8)') STNPNT(1,I)
            CHTMP = DUMMY(3:8)
            READ(UNIT=CHTMP, FMT='(I6)') IBOND
C            WRITE(IW,*)'IBOND=',IBOND
C            NOW IBOND CONTAINS THE ATOM UMBERS WHICH PRODUCE THE BOND           
C
C             THIS IS A REDUNDANT CASE, JUST TO BE SURE
            IF(IBOND.EQ.0) THEN 
               LATM1 = 0
               LATM2 =0
               RADIUS(I)=1.40D+00
            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
              ZNUC1=ZAN(LATM1) + IZCORE(LATM1)
              IZ1=INT(ZNUC1)
              ZNUC2=ZAN(LATM2) + IZCORE(LATM2)
              IZ2=INT(ZNUC2)
              RADIUS(I)=(VANDER(IZ1,IVDW)+VANDER(IZ2,IVDW))/TWO
            END IF
            WRITE(IW,122) LATM1,LATM2
 122        FORMAT(1X,'RADIUS FOR BOND MID-POINT BETWEEN ATOMS ',I2,
     *            ' AND ',I2)          
         END IF
         IF(IZ.GT.0 .AND. IZ.LE.36) RADIUS(I)=VANDER(IZ,IVDW)
         IF(RADIUS(I).EQ.ZERO) THEN
            RADIUS(I)=1.8D+00
            IF(MASWRK) WRITE(IW,920) IZ,RADIUS(I)
         END IF
         RADIUS(I)=RADIUS(I)*TOBOHR
         IF(MASWRK) WRITE(IW,930) STNPNT(1,I),RADIUS(I)
   29 CONTINUE
C
C        OTHER CASES
      ELSE
      DO 30 I=1,NAT
         ZNUC = ZAN(I) + IZCORE(I)
         IZ = INT(ZNUC)
         RADIUS(I)=ZERO
         IF(IZ.LE.0) CALL ABRT
         IF(IZ.LE.36) RADIUS(I)=VANDER(IZ,IVDW)
         IF(RADIUS(I).EQ.ZERO) THEN
            RADIUS(I)=1.8D+00
            IF(MASWRK) WRITE(IW,920) IZ,RADIUS(I)
         END IF
         RADIUS(I)=RADIUS(I)*TOBOHR
   30 CONTINUE
      END IF
C
  915 FORMAT(10X,'EFP CENTER      VDW RADIUS')
  920 FORMAT(10X,'THE VDW RADIUS OF ATOMIC NUMBER',I3,
     *           'IS UNKNOWN AND HAS BEEN SET TO',F4.2,' ANGSTROMS')
  930 FORMAT(10X,A8,8X,F4.2,' BOHR')
 1000 FORMAT(' MERZ-KOLLMAN RADII USED FOR CHARGE FITTING')
 1010 FORMAT(' GAVEZZOTTI RADII USED FOR CHARGE FITTING')
 1020 FORMAT(' BRENEMAN (CHELPG) RADII USED FOR CHARGE FITTING')
      RETURN
      END
