C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C 22 DEC 06 - DGF - FIX SPELLING IN COMMON BLOCK GRDPAR
C 22 FEB 06 - LVS - PDCPTS MODIFIED TO BE USED FOR EFP SCREENING
C  2 FEB 06 - MWS - INTEGRAL PARALLELIZATION, ADJUSTMENTS FOR H+I
C 14 NOV 05 - DGF - ADD EPOTEN TO COMPUTE PCM POTENTIAL
C 19 SEP 05 - HL  - ADD DENMUL TO DO DENSITY BASED MULTIPOLAR EXPANSION
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - DFG - EPCALC: FIX INDEXING FOR G FUNCTIONS
C 30 APR 05 - DGF - CHANGES FOR NEW FMO PROPERTIES
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK HERMIT, WERMIT, NSHEL, ROOT, FRGINF
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  9 DEC 03 - OQ  - SYNCH ARGUMENTS TO DIPDER
C 12 DEC 02 - AE  - USE MXPTPT TO DIMENSION COMMON
C 17 APR 02 - MWS - SYNCH UP FRGINF COMMON
C 25 JUN 01 - MWS - DELETE SOME SPURIOUS COMMENTS
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 10 APR 00 - MWS - REMOVE STATIC MEMORY FROM COSMO COMMONS
C 25 MAR 00 - KKB/LNB - EP FOR COSMO SOLVENT IN PRESENCE OF SOLUTE
C 12 SEP 98 - KA -  PDCCAL: SAVE/RESTORE CENTER AND MOMENTS
C  8 AUG 97 - MWS - USE PARAMETER TO SET SIZE OF /POINTS/ COMMON
C 18 DEC 96 - JHJ - EFLD2: 2 NEW ARGUMENTS IN DIPDER CALL
C 18 SEP 96 - MWS - ELPOTS: ALLOW FOR FRAGMENT BONDS
C 13 JUN 96 - VAG - CHANGES TO INTRODUCE CITYP INPUT VARIABLE
C 26 MAY 96 - MWS - PDCCAL: FIX BUG INTRODUCED BY 6 OCT CHANGE
C  3 JAN 96 - MWS - MOVE DIPINT TO INT1 MODULE, CHANGE EFLD1 ROUTINE
C  6 OCT 95 - MWS - ELPOTC: SAVE VALUE OF MAXIMUM NUMBER OF PDC POINTS
C 27 MAR 95 - SPW - CHANGE ARGUMENTS TO DIPDER CALL
C 29 DEC 94 - TLW - ADD EFLD1, EFLD2, AND DNUCEF FOR EXTERNAL EL. FIELD
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 15 JUL 94 - MR  - PDCCAL: FIX QUADRUPOLE CONSTRAINT ERROR
C 17 DEC 93 - MWS - ELPOTC: CHANGE GTBOND'S CALL
C  9 DEC 93 - MWS - ELMOMC: CHANGE DIPOLE DAF RECORD NUMBERS
C 29 NOV 93 - MS  - EPCALC: DON'T TRASH PROPERTY COORDINATES
C 28 JUN 93 - MS  - ELPOTC,PDCCAL: MEP CHANGES
C  4 JUN 92 - TLW - EPCALC,PRCALC: PARALLELIZE INTEGRAL CALCULATIONS
C 20 APR 92 - MK  - ELFLDC: DELETE XYZPRP TO AVOID TRASHING MOMENTS
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  7 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 20 OCT 91 - TLW - ADDED F AND G CAPABILITIES
C 15 AUG 90 - TLW - ADD H7(7) IN COMMON HERMIT AND W7(7) IN COMMON
C                   WERMIT.  ADD 7TH ELEMENT IN MINP, MAXP, MIN, AND
C                   MAX ARRAYS.
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C  8 MAR 90 - MWS - ELPOTC: CHANGE FORMAT OF MEPMAP'S INPUT FILE.
C 24 SEP 89 - MWS - MOVE POINT TO MODULE PRPLIB, /MASSES/ IN CALCOM
C 20 SEP 89 - MWS - USE MXATM PARAMETER, ADD EFG TENSOR COMMENTS
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 30 MAY 88 - MWS - USE PARAMETERS TO DIMENSION COMMONS
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C 14 FEB 88 - MWS - INCREASE RYS ROOTS TO 9 ROOTS
C 15 NOV 87 - STE - ELDENC,ELFLDC,ELMOMC,ELPOTC: CHANGE PRINT DEFAULTS
C  7 JUL 86 - MWS - SANITIZE FLOATING POINT CONSTANTS
C 12 JUN 86 - MWS - PUNCH ELECTROSTATIC POTENTIAL GRID DATA
C                   IN PLOTTING COMPATIBLE FORMAT IN ELPOTC
C 27 NOV 85 - STE - ELFLDC: GZY TO GXY; USE GENERIC ABS,EXP,SQRT
C                   ELMOMC: PRINT OXYY,OMXYY
C 27 JUL 85 - MWS - RENAME /XYZINT/ TO /PRPINT/ TO AVOID IBM
C                   CONFLICTS WITH THE XYZINT SUBROUTINE
C 16 JUL 85 - MWS - FIX COMPUTATION OF MONOPOLE FOR CI/MCSCF
C 15 JUL 85 - MWS - FIX INTEFL FOR D BASIS FUNCTIONS
C 11 JUL 85 - MWS - HOLLERITH DATA STATEMENTS
C  3 JUL 85 - MWS - FIX UNITS PRINTOUT FOR EL POTENTIAL
C  1 JUL 85 - MWS - IMPLEMENT ABBOTT LABS PROPERTIES IN GAMESS.
C
C*MODULE PRPEL   *DECK CALCMO
      SUBROUTINE CALCMO(MORB,VEXP,ORBS,IA,L1,VALUE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION ORBS(L1,2),VEXP(2),IA(2)
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                     RETRIEVE MOLECULAR ORBITALS
C
      L3=L1*L1
      CALL DAREAD(IDAF,IODA,ORBS,L3,15,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 PRPEL   *DECK CALCOC
      SUBROUTINE CALCOC(X,Y,Z)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      PARAMETER (ZERO=0.0D+00)
C
C       COMPUTE THE CENTER OF CHARGE.
C       NOTE, THIS ROUTINE DOES NOT ACCOUNT FOR EFP PARTICLES.
C
      X = ZERO
      Y = ZERO
      Z = ZERO
      TC = ZERO
      DO 120  I=1,NAT
         X = X + ZAN(I)*C(1,I)
         Y = Y + ZAN(I)*C(2,I)
         Z = Z + ZAN(I)*C(3,I)
         TC = TC + ZAN(I)
  120 CONTINUE
      X = X/TC
      Y = Y/TC
      Z = Z/TC
      RETURN
      END
C*MODULE PRPEL   *DECK CALCOM
      SUBROUTINE CALCOM(X,Y,Z)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MASSES/ ZMASS(MXATM)
C
      PARAMETER (ZERO=0.0D+00)
C
C       COMPUTE THE CENTER OF MASS.
C       NOTE, THIS ROUTINE DOES NOT ACCOUNT FOR EFP PARTICLES.
C
      X = ZERO
      Y = ZERO
      Z = ZERO
      TMASS = ZERO
      DO 120  I=1,NAT
         X = X + ZMASS(I)*C(1,I)
         Y = Y + ZMASS(I)*C(2,I)
         Z = Z + ZMASS(I)*C(3,I)
         TMASS = TMASS + ZMASS(I)
  120 CONTINUE
      X = X/TMASS
      Y = Y/TMASS
      Z = Z/TMASS
      RETURN
      END
C*MODULE PRPEL   *DECK DELINT
      SUBROUTINE DELINT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      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
C
      PARAMETER (ONE=1.0D+00)
C
C       EVALUATE DELTA FUNCTION INTEGRALS (USED FOR E- DENSITY).
C       FOR SPDFGHI BASIS, -NI- AND -NJ- WILL RUN UP TO 7
C
      XINT1 = ONE
      YINT1 = ONE
      ZINT1 = ONE
C
      PTXI = XP - XI
      PTYI = YP - YI
      PTZI = ZP - ZI
C
      PTXJ = XP - XJ
      PTYJ = YP - YJ
      PTZJ = ZP - ZJ
C
      GO TO (101,102,103,104,105,106,107), NI
C
  107 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  106 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  105 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  104 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  103 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  102 XINT1 = XINT1*PTXI
      YINT1 = YINT1*PTYI
      ZINT1 = ZINT1*PTZI
C
  101 CONTINUE
C
      GO TO (201,202,203,204,205,206,207), NJ
C
  207 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  206 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  205 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  204 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  203 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  202 XINT1 = XINT1*PTXJ
      YINT1 = YINT1*PTYJ
      ZINT1 = ZINT1*PTZJ
C
  201 CONTINUE
      RETURN
      END
C*MODULE PRPEL   *DECK DNUCEF
      SUBROUTINE DNUCEF(EG3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL EFLDL
C
      DIMENSION EG3(*),TFIND(3)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
C     GRADIENT NUCLEAR CONTRIBUTION FROM EXTERNAL ELECTRIC FIELD
C
      TFIND(1) = FIND(1)
      TFIND(2) = FIND(2)
      TFIND(3) = FIND(3)
C
C    MUST TEMPORARILY CHANGE THE SIGN OF THE FIELD.
C    IT IS CHANGED BACK IN EFLD2 WHEN THE GRADIENT
C    CONTRIBUTION FROM THE FIELD IS DONE.
C
      EVEC(1)  = -EVEC(1)
      EVEC(2)  = -EVEC(2)
      EVEC(3)  = -EVEC(3)
      FIND(1)  = EVEC(1)
      FIND(2)  = EVEC(2)
      FIND(3)  = EVEC(3)
C
      CALL DNUCZ(EG3)
C
      FIND(1)  = TFIND(1)
      FIND(2)  = TFIND(2)
      FIND(3)  = TFIND(3)
C
      RETURN
      END
C*MODULE PRPEL   *DECK EFLD1
      SUBROUTINE EFLD1(XMU,YMU,ZMU,H,L2,SOME)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL EFLDL,SOME
C
      DIMENSION H(L2),XMU(L2),YMU(L2),ZMU(L2)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- ADD ELECTRIC FIELD TO ONE ELECTRON HAMILTONIAN -----
C
      IF(SOME) WRITE(IW,9000) EVEC
C
C     DIPOLE INTEGRALS MUST BE EVALUATED WRT COORDINATE ORIGIN
C
      CALL DAREAD(IDAF,IODA,H  ,L2,11,0)
      CALL DAREAD(IDAF,IODA,XMU,L2,95,0)
      CALL DAREAD(IDAF,IODA,YMU,L2,96,0)
      CALL DAREAD(IDAF,IODA,ZMU,L2,97,0)
C
      IF(EVEC(1).NE.ZERO) CALL DAXPY(L2,EVEC(1),XMU,1,H,1)
      IF(EVEC(2).NE.ZERO) CALL DAXPY(L2,EVEC(2),YMU,1,H,1)
      IF(EVEC(3).NE.ZERO) CALL DAXPY(L2,EVEC(3),ZMU,1,H,1)
C
      CALL DAWRIT(IDAF,IODA,H,L2,11,0)
      RETURN
C
 9000 FORMAT(1X,'ELECTRIC FIELD =',3F10.5/
     *       1X,'HAS BEEN APPLIED TO ONE ELECTRON HAMILTONIAN')
      END
C*MODULE PRPEL   *DECK EFLD2
      SUBROUTINE EFLD2(DA,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL EFLDL
C
      DIMENSION TFIND(3)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      DATA EFL/8HEFIELD  /
C
C     ----- EXTERNAL ELECTRIC FIELD CONTRIBUTION TO THE GRADIENT -----
C
      TFIND(1) = FIND(1)
      TFIND(2) = FIND(2)
      TFIND(3) = FIND(3)
      FIND(1)  = EVEC(1)
      FIND(2)  = EVEC(2)
      FIND(3)  = EVEC(3)
      GNUCF = NE+ICH
      GNUCF = GNUCF/NE
C
      DUMMY = 0.0D+00
      CALL DIPDER(DA,DUMMY,DUMMY,DUMMY,L2,DUMMY,0,0,0,EFL)
C
      FIND(1)  = TFIND(1)
      FIND(2)  = TFIND(2)
      FIND(3)  = TFIND(3)
      EVEC(1)  = -EVEC(1)
      EVEC(2)  = -EVEC(2)
      EVEC(3)  = -EVEC(3)
      RETURN
      END
C*MODULE PRPEL   *DECK ELDENC
      SUBROUTINE ELDENC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PAPER,PUNCH,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000)
C
      DIMENSION ELDEN(1)
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 /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 /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 (ONE=1.0D+00)
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HELEDNC  /
      DATA ELDEN/8HELDEN   /
C
C       CALCULATE ELECTRON DENSITY FOR MOLECULAR ORBITAL OR
C       ENTIRE WAVE FUNCTION AT SPECIFIED POINTS.
C
      PAPER = IEDOUT.GE.0. AND. MASWRK
      IF(EXETYP.EQ.DEBUG)  PAPER=.TRUE.
      IF(EXETYP.EQ.DBUGME) PAPER=.TRUE.
      PAPER = PAPER .AND. MASWRK
C
      PUNCH = IEDOUT.NE.0
      PUNCH = PUNCH .AND. MASWRK
C
      IF(PAPER) WRITE(IW,920)
      IF(PAPER) WRITE(IW,930)
      IF(PUNCH) WRITE(IP,940)
C
C               SET POINTERS FOR PARTITIONING OF MEMORY
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      NVAL = 1
C
      CALL VALFM(LOADFM)
      IELD   = LOADFM + 1
      IEDW   = IELD   + L2
C
      IDENSA = IEDW   + 784
      IDENSB = IDENSA + L2
      LAST1  = IDENSB + L2
C
      IMOS   = IEDW   + 784
      INDEX  = IMOS   + L3
      LAST2  = INDEX  + L1
C
      LAST = MAX(LAST1,LAST2)
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      IF (EXETYP.EQ.CHECK) THEN
         CALL RETFM(NEED)
         RETURN
      END IF
C
C        OBTAIN TOTAL DENSITY MATRIX, IF NEEDED
C
      IF(MODENS.EQ.0) CALL DENDD1(X(IDENSA),X(IDENSB),L2)
C
C              LOOP OVER POINTS TO BE EVALUATED
C
      IPOINT = 0
  210 CONTINUE
      IPOINT = IPOINT + 1
      ISTAT = 0
      CALL POINT(ELDLOC,IPOINT,XP,YP,ZP,ISTAT)
      IF(ISTAT.LT.0) CALL ABRT
      IF(ISTAT.GT.0) THEN
         CALL RETFM(NEED)
         RETURN
C        ******
      END IF
C
C          CALCULATE DELTA FUNCTION INTEGRALS
C
      CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,GOPARR)
C
      IF(IEDINT.NE.0) THEN
         IF(GOPARR) CALL DDI_GSUMF(910,X(IELD),L2)
         IF(MASWRK) WRITE(IW,950)
         CALL PRTRIL(X(IELD),L1)
         IF(GOPARR) CALL DSCAL(L2,ONE/NPROC,X(IELD),1)
      END IF
C
C        CALCULATE ELECTRON DENSITY FOR ORBITAL OR WAVEFUNCTION
C
      IF(MODENS.EQ.0) THEN
         EDENS = TRACEP(X(IDENSA),X(IELD),L1)
      ELSE
         EDENS = 0.0D+00
         CALL CALCMO(MODENS,X(IELD),X(IMOS),X(INDEX),L1,EDENS)
      END IF
C
      IF(GOPARR) CALL DDI_GSUMF(911,EDENS,1)
C
C        OUTPUT THE RESULTS
C
      IF(PAPER) WRITE(IW,970) IPOINT,XP,YP,ZP,EDENS
      IF(PUNCH) WRITE(IP,980) IPOINT,XP,YP,ZP,EDENS
      GO TO 210
C
  920 FORMAT(/10X,16('-')/10X,'ELECTRON DENSITY'/10X,16('-'))
  930 FORMAT(2X,'POINT',6X,'X',9X,'Y',9X,'Z',13X,'DENSITY',
     *        /1H ,20X,'(BOHR)',20X,'(A.U.)')
  940 FORMAT(1X,'ELECTRON DENSITY, IPOINT,X,Y,Z,EDENS')
  950 FORMAT(1H ,10X,22('-'),/11X,'ELECTRON DENSITY INTEGRALS')
  970 FORMAT(1X,I5,2X,3F10.5,F17.6)
  980 FORMAT(1X,I3,2X,3F10.5,E15.6)
      END
C*MODULE PRPEL   *DECK ELFLDC
      SUBROUTINE ELFLDC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION NUCLEI
C
      LOGICAL PUNCH,PAPER,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ENAME(3),ENAM(3,2),TAGS(9),EFGEL(9)
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 /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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /XYZPRP/ XP,YP,ZP,DMY(35)
C
      PARAMETER (ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00, TM3=1.0D-03)
C
      DATA CHECK/8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HELFLDC  /
      DATA ENAM   /8HELECTRIC, 8H FIELD  , 8H        ,
     *             8HELECTRIC, 8H FIELD/G, 8HRADIENT /
      DATA ELFLD  /8HELFLD   /,
     *     NUCLEI /8HNUCLEI  /
      DATA TAGS/8H-X/R**3 ,8H-Y/R**3 ,8H-Z/R**3 ,8H-XX/R**5,
     *          8H-YY/R**5,8H-ZZ/R**5,8H-XY/R**5,8H-XZ/R**5,
     *          8H-YZ/R**5/
C
C     ----- CALCULATE ELECTRIC FIELD/GRADIENT -----
C
      IF(IEFLD.EQ.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ELFLDC: IMPOSSIBLE -IEFLD- VALUE'
         CALL ABRT
      END IF
C
      DO 10 I=1,3
         ENAME(I) = ENAM(I,IEFLD)
   10 CONTINUE
C
      PAPER = IEFOUT.GE.0
      IF(EXETYP.EQ.DEBUG)  PAPER=.TRUE.
      IF(EXETYP.EQ.DBUGME) PAPER=.TRUE.
      PAPER = PAPER .AND. MASWRK
C
      PUNCH = IEFOUT.NE.0
      PUNCH = PUNCH .AND. MASWRK
C
      IF(PUNCH) WRITE(IP,905) ENAME
      IF(PAPER) WRITE(IW,910) ENAME
      IF(PAPER  .AND.   ELFLOC.EQ.NUCLEI) WRITE(IW,920)
C
C          PARTITION MEMORY
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
                     NVAL=0
      IF(IEFLD.EQ.1) NVAL=3
      IF(IEFLD.EQ.2) NVAL=9
      IF (NVAL.EQ.0) THEN
         IF(MASWRK) WRITE(6,*) 'ELFLDC: IMPOSSIBLE -IEFLD- VALUE'
         CALL ABRT
      END IF
C
      CALL VALFM(LOADFM)
      IELF   = LOADFM + 1
      IEFW   = IELF   + NVAL*L2
      IDENSA = IEFW   + 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 EACH POINT TO BE COMPUTED -----
C
      IPOINT = 0
  210 CONTINUE
      IPOINT = IPOINT + 1
      ISTAT = 0
      CALL POINT(ELFLOC,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
      IF(PAPER) WRITE(IW,930) IPOINT,XP,YP,ZP
C
C             COMPUTE INTEGRALS OVER BASIS FUNCTIONS
C
      CALL PRCALC(ELFLD,X(IELF),X(IEFW),NVAL,L2,GOPARR)
C
      IF(IEFINT.EQ.0) GO TO 240
         MIN = 1
         MAX = 9
         IF(IEFINT.EQ.1) MAX=3
         IF(IEFINT.EQ.-2) MIN = 4
         IF(GOPARR) CALL DDI_GSUMF(912,X(IELF),MAX*L2)
         LOC=IELF+L2*(MIN-1)
         DO 230 I=MIN,MAX
            IF(MASWRK) WRITE(IW,980) TAGS(I)
            CALL PRTRI(X(LOC),L1)
            LOC=LOC+L2
  230    CONTINUE
         IF(GOPARR) CALL DSCAL(MAX*L2,ONE/NPROC,X(IELF),1)
  240 CONTINUE
C
C            EXPECTATION VALUE FOR ELECTRIC FIELD
C
      EFGEL(1) = TRACEP(X(IDENSA),X(IELF     ),L1)
      EFGEL(2) = TRACEP(X(IDENSA),X(IELF+  L2),L1)
      EFGEL(3) = TRACEP(X(IDENSA),X(IELF+2*L2),L1)
C
C            EXPECTATION VALUE FOR ELECTRIC FIELD GRADIENT
C
      IF(IEFLD.EQ.2) THEN
         EFGEL(4) = TRACEP(X(IDENSA),X(IELF+3*L2),L1)
         EFGEL(5) = TRACEP(X(IDENSA),X(IELF+4*L2),L1)
         EFGEL(6) = TRACEP(X(IDENSA),X(IELF+5*L2),L1)
         EFGEL(7) = TRACEP(X(IDENSA),X(IELF+6*L2),L1)
         EFGEL(8) = TRACEP(X(IDENSA),X(IELF+7*L2),L1)
         EFGEL(9) = TRACEP(X(IDENSA),X(IELF+8*L2),L1)
      END IF
C
      IF(GOPARR) CALL DDI_GSUMF(913,EFGEL,9)
C
      EFX  = EFGEL(1)
      EFY  = EFGEL(2)
      EFZ  = EFGEL(3)
      EFZ  = EFGEL(3)
      EGXX = EFGEL(4)
      EGYY = EFGEL(5)
      EGZZ = EFGEL(6)
      EGXY = EFGEL(7)
      EGXZ = EFGEL(8)
      EGYZ = EFGEL(9)
C
C          ADD NUCLEAR CONTRIBUTIONS
C
      DO 290  I=1,NAT
         ZNUC = ZAN(I)
         XN = C(1,I) - XP
         YN = C(2,I) - YP
         ZN = C(3,I) - ZP
         RR = XN*XN + YN*YN + ZN*ZN
         IF(RR.LT.TM3  .AND.  ELFLOC.NE.NUCLEI .AND. MASWRK)
     *               WRITE(IW,915) ENAME,XP,YP,ZP,I
         IF(RR.LT.TM3) GO TO 290
         RR3 = RR*SQRT(RR)
         EFX = EFX - ZNUC*XN/RR3
         EFY = EFY - ZNUC*YN/RR3
         EFZ = EFZ - ZNUC*ZN/RR3
         IF(IEFLD.EQ.2) THEN
            RR5  = RR*RR3
            EGXX = EGXX - ZNUC*XN*XN/RR5
            EGYY = EGYY - ZNUC*YN*YN/RR5
            EGZZ = EGZZ - ZNUC*ZN*ZN/RR5
            EGXY = EGXY - ZNUC*XN*YN/RR5
            EGXZ = EGXZ - ZNUC*XN*ZN/RR5
            EGYZ = EGYZ - ZNUC*YN*ZN/RR5
         END IF
  290 CONTINUE
C
      EFLD = SQRT(EFX*EFX + EFY*EFY + EFZ*EFZ)
C
C           FORM FIELD GRADIENT TENSOR
C     THIS TENSOR IS THE HESSIAN OF THE ELECTROSTATIC POTENTIAL, 1/R
C       EF(X) = D(1/R)/DX = -X/R**3
C     EFG(XX) = D(-X/R**3)/DX = 3X**2/R**5 - 1/R**3
C             = (2X**2 - Y**2 - Z**2)/R**5
C     EFG(XY) = 3XY/R**5, ETC.
C
      IF(IEFLD.EQ.2) THEN
         GXX = TWO*EGXX - EGYY - EGZZ
         GYY = TWO*EGYY - EGXX - EGZZ
         GZZ = TWO*EGZZ - EGXX - EGYY
         GXY = THREE*EGXY
         GXZ = THREE*EGXZ
         GYZ = THREE*EGYZ
      END IF
C
C          PRINT/PUNCH RESULTS
C
      IF(PAPER) WRITE(IW,935) EFX,EFY,EFZ,EFLD
      IF(PUNCH) WRITE(IP,950) IPOINT,XP,YP,ZP,EFX,EFY,EFZ,EFLD
      IF(PAPER .AND. IEFLD.EQ.2) WRITE(IW,940) GXX,GYY,GZZ,GXY,GXZ,GYZ
      IF(PUNCH .AND. IEFLD.EQ.2) WRITE(IP,960) GXX,GYY,GZZ,GXY,GXZ,GYZ
      GO TO 210
C
  905 FORMAT(2X,3A8)
  910 FORMAT(/10X,23('-'),/10X,3A8,/10X,23('-')/)
  915 FORMAT(/1X,'*** WARNING - ',3A8,'AT ',3F10.5,
     *       '.  CONTRIBUTION FROM NUCLEUS ',I3,' IGNORED')
  920 FORMAT(2X,'THE CONTRIBUTION FROM EACH NUCLEUS AT ITS ',
     *          'OWN POSITION IS NEGLECTED')
  930 FORMAT(/1X,'POINT',I4,11X,'X',11X,'Y',11X,'Z   (BOHR)'/
     *       13X,3F12.6)
  935 FORMAT(8X,'EFX',9X,'EFY',9X,'EFZ',9X,
     *       '/EF/      (A.U.)'/1X,4F12.6)
  940 FORMAT(7X,'EFGXX',7X,'EFGYY',7X,'EFGZZ',7X,'EFGXY',7X,
     *          'EFGXZ',7X,'EFGYZ'/1X,6F12.6)
  950 FORMAT(1X,'ELECTRIC FIELD/GRADIENT',I4,3F11.6/1X,4F12.7)
  960 FORMAT(1X,6F12.7)
  980 FORMAT(/10X,A8,' INTEGRALS'/10X,18(1H-))
      END
C*MODULE PRPEL   *DECK ELMOMC
      SUBROUTINE ELMOMC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PAPER,PUNCH,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ITAG(19),ELTEMP(16)
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 /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,QFAC=1.345044D+00 ,OFAC=7.117668D-01)
      PARAMETER (PT5=0.5D+00, ONE=1.0D+00, ONEPT5=1.5D+00,
     *           TWO=2.0D+00, TWOPT5=2.5D+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
      NSUM=0
      IF(IEMOM.GE.2) THEN
         ELTEMP( 1) = -TRACEP(X(IDENSA),X(IELM+ 3*L2),L1)
         ELTEMP( 2) = -TRACEP(X(IDENSA),X(IELM+ 4*L2),L1)
         ELTEMP( 3) = -TRACEP(X(IDENSA),X(IELM+ 5*L2),L1)
         ELTEMP( 4) = -TRACEP(X(IDENSA),X(IELM+ 6*L2),L1)
         ELTEMP( 5) = -TRACEP(X(IDENSA),X(IELM+ 7*L2),L1)
         ELTEMP( 6) = -TRACEP(X(IDENSA),X(IELM+ 8*L2),L1)
         NSUM = NSUM+6
      END IF
C                         OCTUPOLE
      IF(IEMOM.GE.3) THEN
         ELTEMP( 7) = -TRACEP(X(IDENSA),X(IELM+ 9*L2),L1)
         ELTEMP( 8) = -TRACEP(X(IDENSA),X(IELM+10*L2),L1)
         ELTEMP( 9) = -TRACEP(X(IDENSA),X(IELM+11*L2),L1)
         ELTEMP(10) = -TRACEP(X(IDENSA),X(IELM+12*L2),L1)
         ELTEMP(11) = -TRACEP(X(IDENSA),X(IELM+13*L2),L1)
         ELTEMP(12) = -TRACEP(X(IDENSA),X(IELM+14*L2),L1)
         ELTEMP(13) = -TRACEP(X(IDENSA),X(IELM+15*L2),L1)
         ELTEMP(14) = -TRACEP(X(IDENSA),X(IELM+16*L2),L1)
         ELTEMP(15) = -TRACEP(X(IDENSA),X(IELM+17*L2),L1)
         ELTEMP(16) = -TRACEP(X(IDENSA),X(IELM+18*L2),L1)
         NSUM = NSUM+10
      END IF
C
      IF(GOPARR  .AND.  NSUM.GT.0) CALL DDI_GSUMF(916,ELTEMP,NSUM)
C
      QMXX = ELTEMP(1)
      QMYY = ELTEMP(2)
      QMZZ = ELTEMP(3)
      QMXY = ELTEMP(4)
      QMXZ = ELTEMP(5)
      QMYZ = ELTEMP(6)
C
      OMXXX = ELTEMP( 7)
      OMYYY = ELTEMP( 8)
      OMZZZ = ELTEMP( 9)
      OMXXY = ELTEMP(10)
      OMXXZ = ELTEMP(11)
      OMXYY = ELTEMP(12)
      OMYYZ = ELTEMP(13)
      OMXZZ = ELTEMP(14)
      OMYZZ = ELTEMP(15)
      OMXYZ = ELTEMP(16)
C
C                 ADD 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
C
         IF(IEMOM .LT. 2) GO TO 280
         QMXX = QMXX + ZNUC*XN*XN
         QMYY = QMYY + ZNUC*YN*YN
         QMZZ = QMZZ + ZNUC*ZN*ZN
         QMXY = QMXY + ZNUC*XN*YN
         QMXZ = QMXZ + ZNUC*XN*ZN
         QMYZ = QMYZ + ZNUC*YN*ZN
C
         IF(IEMOM .LT. 3) GO TO 280
         OMXXX = OMXXX + ZNUC*XN*XN*XN
         OMYYY = OMYYY + ZNUC*YN*YN*YN
         OMZZZ = OMZZZ + ZNUC*ZN*ZN*ZN
         OMXXY = OMXXY + ZNUC*XN*XN*YN
         OMXXZ = OMXXZ + ZNUC*XN*XN*ZN
         OMXYY = OMXYY + ZNUC*XN*YN*YN
         OMYYZ = OMYYZ + ZNUC*YN*YN*ZN
         OMXZZ = OMXZZ + ZNUC*XN*ZN*ZN
         OMYZZ = OMYZZ + ZNUC*YN*ZN*ZN
         OMXYZ = OMXYZ + ZNUC*XN*YN*ZN
  280 CONTINUE
C
C     ----- CONVERT UNITS, FORM QUADRUPOLE,OCTUPOLE TENSORS -----
C
      DMX = DMX*DFAC
      DMY = DMY*DFAC
      DMZ = DMZ*DFAC
      DIP = SQRT(DMX*DMX+DMY*DMY+DMZ*DMZ)
C
      IF(IEMOM. LT. 2) GO TO 300
         QMXX = QMXX*QFAC
         QMYY = QMYY*QFAC
         QMZZ = QMZZ*QFAC
         QMXY = QMXY*QFAC
         QMXZ = QMXZ*QFAC
         QMYZ = QMYZ*QFAC
C
         QXX = PT5*(QMXX + QMXX - QMYY - QMZZ)
         QYY = PT5*(QMYY + QMYY - QMXX - QMZZ)
         QZZ = PT5*(QMZZ + QMZZ - QMXX - QMYY)
         QXY = ONEPT5*QMXY
         QXZ = ONEPT5*QMXZ
         QYZ = ONEPT5*QMYZ
C
      IF(IEMOM .LT. 3) GO TO 300
         OMXXX = OMXXX*OFAC
         OMYYY = OMYYY*OFAC
         OMZZZ = OMZZZ*OFAC
         OMXXY = OMXXY*OFAC
         OMXXZ = OMXXZ*OFAC
         OMXYY = OMXYY*OFAC
         OMYYZ = OMYYZ*OFAC
         OMXZZ = OMXZZ*OFAC
         OMYZZ = OMYZZ*OFAC
         OMXYZ = OMXYZ*OFAC
C
         OXXX =     OMXXX - ONEPT5*OMXYY - ONEPT5*OMXZZ
         OYYY =     OMYYY - ONEPT5*OMXXY - ONEPT5*OMYZZ
         OZZZ =     OMZZZ - ONEPT5*OMXXZ - ONEPT5*OMYYZ
         OXXY = TWO*OMXXY - PT5*OMYYY - PT5*OMYZZ
         OXXZ = TWO*OMXXZ - PT5*OMYYZ - PT5*OMZZZ
         OXYY = TWO*OMXYY - PT5*OMXXX - PT5*OMXZZ
         OYYZ = TWO*OMYYZ - PT5*OMXXZ - PT5*OMZZZ
         OXZZ = TWO*OMXZZ - PT5*OMXXX - PT5*OMXYY
         OYZZ = TWO*OMYZZ - PT5*OMXXY - PT5*OMYYY
         OXYZ = TWOPT5*OMXYZ
C
C     ---- OUTPUT RESULTS -----
C
  300 CONTINUE
      IF(.NOT.PAPER) GO TO 310
      WRITE(IW,930) IPOINT,XP,YP,ZP,CHARGE,DMX,DMY,DMZ,DIP
      IF(IEMOM.GT.1) WRITE(IW,940) QXX,QYY,QZZ,QXY,QXZ,QYZ
      IF(IEMOM.GT.2) WRITE(IW,950) OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                             OXZZ,OYZZ,OZZZ,OXYZ
  310 CONTINUE
      IF(.NOT.PUNCH) GO TO 210
      WRITE(IP,960) IPOINT,XP,YP,ZP,DMX,DMY,DMZ
      IF(IEMOM.GT.1) WRITE(IP,970) QXX,QYY,QZZ,QXY,QXZ,QYZ
      IF(IEMOM.GT.1) WRITE(IP,975) QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ
      IF(IEMOM.GT.2) WRITE(IP,980) OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                             OXZZ,OYZZ,OZZZ,OXYZ
      IF(IEMOM.GT.2) WRITE(IP,985) OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                             OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
      GO TO 210
C
  920 FORMAT(/10X,21('-'),/10X,'ELECTROSTATIC MOMENTS',
     *      /10X,21('-'),/)
  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)
  940 FORMAT(8X,'QXX',9X,'QYY',9X,'QZZ',9X,'QXY',9X,'QXZ',
     *       9X,'QYZ  (BUCKINGHAMS)'/1X,6F12.6)
  950 FORMAT(8X,'OXXX',8X,'OXXY',8X,'OXXZ',8X,'OXYY',8X,'OYYY',
     *       8X,'OYYZ'/1X,6F12.6/8X,'OXZZ',8X,'OYZZ',8X,'OZZZ',
     *       8X,'OXYZ  (X10**34 ESU-CM**3)'/1X,4F12.6)
  960 FORMAT(1X,'MOMENTS AT POINT',I5,' X,Y,Z=',3F10.6/
     *       1X,'DIPOLE     ',3F10.6)
  970 FORMAT(1X,'QUADRUPOLE ',6F10.6)
  975 FORMAT(1X,'2ND MOMENTS',6F10.6)
  980 FORMAT(1X,'OCTUPOLE   ',6F10.6/11X,4F10.6)
  985 FORMAT(1X,'3RD MOMENTS',6F10.6/11X,4F10.6)
  990 FORMAT(/10X,A4,' INTEGRALS'/11X,13(1H-))
      END
C*MODULE PRPEL   *DECK ELPOTC
      SUBROUTINE ELPOTC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION NUCLEI
C
      LOGICAL OUT,PAPER,PUNCH,NOGRID,NOPDC,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ROT(3,3)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXPTPT=100, MXAO=8192,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
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 /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GRDPAR/ ORIGIN(3),XVEC(3),YVEC(3),ZVEC(3),UX(3),UY(3),
     *                UZ(3),GRDSIZ,NGRID,IGUNIT,NXG,NYG,NZG,MODGRID
      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 /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TM3=1.0D-03)
C
      DATA CHECK/8HCHECK   /
      DATA COMASS/8HCOMASS  /,
     *     NUCLEI/8HNUCLEI  /, POINTS/8HPOINTS  /,
     *     GRID  /8HGRID    /,    PDC/8HPDC     /,
     *     CNONE /8HNONE    /
C
C       ----- CALCULATE ELECTROSTATIC POTENTIAL -----
C
C      IN ORDER TO INCREASE THE SPEED OF CALCULATION,
C      THE LOOP OVER POINTS IS PULLED THROUGH
C      THE LOOPS OVER SHELLS AND PRIMITIVES.
C      INSTEAD OF CALLING ROUTINE "PRCALC", "EPCALC" IS CALLED.
C
      NOGRID= ELPLOC.NE.GRID
      NOPDC = ELPLOC.NE.PDC
      OUT = NOPDC  .OR.  CONSTR.EQ.CNONE
C
      PAPER = IEPOUT.GE.0
      PAPER = PAPER .AND. MASWRK
C
      PUNCH = IEPOUT.NE.0
      PUNCH = PUNCH .AND. MASWRK
C
      IF(PAPER) WRITE(IW,910)
      IF(PAPER  .AND.  ELPLOC.EQ.NUCLEI) WRITE(IW,920)
      IF(PAPER  .AND.  OUT) WRITE(IW,930)
      IF(PUNCH) WRITE(IP,915)
C
C               SET POINTERS FOR PARTITIONING OF MEMORY
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
C
                           MAXPNT=0
      IF(ELPLOC.EQ.COMASS) MAXPNT=1
      IF(ELPLOC.EQ.NUCLEI) MAXPNT=NAT
      IF(ELPLOC.EQ.POINTS) MAXPNT=NPOINT
      IF(ELPLOC.EQ.GRID  ) MAXPNT=NGRID
      IF(ELPLOC.EQ.PDC   ) MAXPNT=NPDC
      IF(MAXPNT.EQ.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ELPOTC: IMPOSSIBLE NO. OF POINTS'
         CALL ABRT
      END IF
C
      CALL VALFM(LOADFM)
C
C        FIRST, WE MUST GENERATE THE POINTS FOR WHERE=PDC.
C        ALL OTHER LOCATIONS ARE KNOWN, OR GENERATED ON THE FLY BELOW.
C
      IF(ELPLOC.EQ.PDC) THEN
         IXYZ   = LOADFM + 1
         LRAD   = IXYZ   + 3*MAXPNT
         LCICOS = LRAD   + NAT
         LNLAYER = LCICOS + 3*2562
         LAST   = LNLAYER + LAYER + 1
         NEED   = LAST - LOADFM - 1
         CALL GETFM(NEED)
         NPDCSV = NPDC
         CALL PDCPTS(X(IXYZ),X(LRAD),X(LCICOS),MAXPNT,NAT,X(LNLAYER))
         MAXPNT = NPDC
         NPDC   = NPDCSV
         CALL RETFM(NEED)
      END IF
C
C     ALLOCATE MEMORY, NOTE -IXYZ- NEEDS TO BE FIRST IF -PDC- USED
C
      MAXBND = 5*NAT + 20*NFRG
C
      IXYZ   = LOADFM + 1
      IDENSA = IXYZ   + 3*MAXPNT
      IDENSB = IDENSA + L2
      IPOT   = IDENSB + L2
      IPAIRS = IPOT   + MAXPNT
      LAST   = IPAIRS + MAXBND
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      IF(EXETYP.EQ.CHECK) GO TO 700
C
C        GET TOTAL DENSITY MATRIX
C
      CALL DENDD1(X(IDENSA),X(IDENSB),L2)
C
C              SET UP ALL POINTS TO BE EVALUATED
C
      IF(ELPLOC.NE.PDC) THEN
         IPOINT = 0
         INDEX  = IXYZ - 3
  140    CONTINUE
         IPOINT = IPOINT + 1
         XPT = ZERO
         YPT = ZERO
         ZPT = ZERO
         ISTAT=0
         CALL POINT(ELPLOC,IPOINT,XPT,YPT,ZPT,ISTAT)
         IF(ISTAT.LT.0) CALL ABRT
         IF(ISTAT.GT.0) GO TO 170
         INDEX = INDEX + 3
         X(INDEX)   = XPT
         X(INDEX+1) = YPT
         X(INDEX+2) = ZPT
         GO TO 140
      END IF
C
C        CALCULATE ELECTRONIC CONTRIBUTION AT ALL POINTS
C
  170 CONTINUE
      CALL EPCALC(X(IXYZ),X(IPOT),X(IDENSA),MAXPNT,L2)
C
C        OUTPUT RESULTS OF ALL POINTS
C
      INDEX = IXYZ - 3
      DO 240  IPT=1,MAXPNT
         ELPOTE = -X(IPOT+IPT-1)
         INDEX  = INDEX + 3
         XPT = X(INDEX)
         YPT = X(INDEX+1)
         ZPT = X(INDEX+2)
C                 FORM NUCLEAR CONTRIBUTION
         ELPOTN = ZERO
         DO 210 I = 1,NAT
            XN = C(1,I) - XPT
            YN = C(2,I) - YPT
            ZN = C(3,I) - ZPT
            RR = SQRT(XN*XN + YN*YN + ZN*ZN)
            IF(RR.LT.TM3) THEN
              IF(ELPLOC.NE.NUCLEI .AND. MASWRK)
     *                  WRITE(IW,950) XPT,YPT,ZPT,I
              GO TO 210
            END IF
            ELPOTN = ELPOTN + ZAN(I)/RR
  210    CONTINUE
C
         ELPOTT = ELPOTE + ELPOTN
         X(IPOT+IPT-1) = ELPOTT
         IF(PAPER  .AND.  OUT)
     *      WRITE(IW,960) IPT,XPT,YPT,ZPT,ELPOTE,ELPOTN,ELPOTT
         IF(NOGRID  .AND.  PUNCH  .AND.  OUT)
     *      WRITE(IP,970) IPT,XPT,YPT,ZPT,ELPOTT
  240 CONTINUE
C
C        PUNCH THE ELECTROSTATIC POTENTIAL GRID WITH ALL THE
C        NECESSARY INFORMATION FOR A CONTOUR PLOT FROM -MEPMAP-
C
      IF(NOGRID  .OR.  .NOT.MASWRK) GO TO 700
C
      WRITE(IP,8000) 'START'
      CALL GTBOND(X(IPAIRS),MAXBND,NBONDS,0)
      WRITE(IP,8010) NAT,NBONDS,TITLE
      CALL GTBOND(X(IPAIRS),MAXBND,NBONDS,IP)
C
      WRITE(IP,8020)
      DO 300 I=1,NAT
         WRITE(IP,8030) (C(J,I),J=1,3)
  300 CONTINUE
C
      WRITE(IP,8040) (ORIGIN(I),I=1,3)
      FAC = ONE/GRDSIZ
      DO 310 I=1,3
         ROT(1,I) = FAC*UX(I)
         ROT(2,I) = FAC*UY(I)
         ROT(3,I) = FAC*UZ(I)
  310 CONTINUE
      WRITE(IP,8050) ((ROT(I,J),I=1,3),J=1,3)
      XMIN=ZERO
      XMAX=(NXG-1)*GRDSIZ
      YMIN=ZERO
      YMAX=(NYG-1)*GRDSIZ
      WRITE(IP,8060) XMIN,XMAX,YMIN,YMAX
      WRITE(IP,8070) NXG,NYG
      DO 330 I=1,NXG
         JMAX=0
         KARD=0
         LOC=(I-1)*NYG - 1
  320    CONTINUE
            JMIN=JMAX+1
            JMAX=JMAX+5
            KARD=KARD+1
            IF(JMAX.GT.NYG) JMAX=NYG
            WRITE(IP,8080) I,KARD,(X(IPOT+LOC+J),J=JMIN,JMAX)
         IF(JMAX.LT.NYG) GO TO 320
  330 CONTINUE
      WRITE(IP,8000) '  END'
C
  700 CONTINUE
      IF(ELPLOC.EQ.PDC  .AND.  CONSTR.NE.CNONE) THEN
         NEQ = NAT + 9
         CALL VALFM(LOADFM)
         LA    = LOADFM + 1
         LQ    = LA     + NEQ*NEQ
         LIPVT = LQ     + NEQ
         LWORK = LIPVT  + NEQ
         LAST  = LWORK  + NEQ
         NEED2 = LAST - LOADFM - 1
         CALL GETFM(NEED2)
         CALL PDCCAL(X(IXYZ),X(IPOT),X(LA),X(LQ),
     *               X(LIPVT),X(LWORK),MAXPNT,NEQ)
         CALL RETFM(NEED2)
      END IF
C
      CALL RETFM(NEED)
      RETURN
C
  910 FORMAT(/10X,23('-')/10X,'ELECTROSTATIC POTENTIAL',
     *          /10X,23('-'),/)
  915 FORMAT(1X,'ELECTRIC POTENTIAL, IPT,X,Y,Z,ELPOTT')
  920 FORMAT(2X,'THE CONTRIBUTION FROM EACH NUCLEUS AT ITS ',
     *          'OWN POSITION IS NEGLECTED'/)
  930 FORMAT(1X,'POINT',3X,'X',9X,'Y',9X,'Z',8X,
     *          'ELECTRONIC',5X,'NUCLEAR',7X,'TOTAL'/
     *          16X,'(BOHR)',31X,'(A.U.)')
  950 FORMAT(/1H ,'*** WARNING - ELECTROSTATIC POTENTIAL AT ',
     *      3F10.5,'. CONTRIBUTION FROM NUCLEUS ',I3,' IGNORED',/)
  960 FORMAT(1X,I4,F9.5,2F10.5,3F13.6)
  970 FORMAT(1X,I5,2X,3F10.5,E15.6)
C
 8000 FORMAT('----- ',A5,' OF INPUT FILE FOR -MEPMAP- PLOTTING',
     *       ' PROGRAM -----')
 8010 FORMAT('NATOMS=',I4,'  NBONDS=',I4,'  CUTOFF=60.5 DELTA=5.0'/
     *        10A8)
 8020 FORMAT('ATOMIC X,Y,Z COORDINATES (BOHRS)')
 8030 FORMAT(3F16.10)
 8040 FORMAT('TRANSLATION TO GRID ORIGIN'/3(1P,E17.8))
 8050 FORMAT('ROTATION TO PLOTTING PLANE'/
     *       3(1P,E17.8)/3(1P,E17.8)/3(1P,E17.8))
 8060 FORMAT('MINIMA AND MAXIMA ON PLOTTING PLANE'/4(1P,E17.8))
 8070 FORMAT(2I5,' X,Y MESH POINTS'/
     *   'THE ROWS OF THE ELECTROSTATIC POTENTIAL GRID (A.U.)')
 8080 FORMAT(I3,I2,5E15.7)
      END
C*MODULE PRPEL   *DECK EPCALC
      SUBROUTINE EPCALC(XYZ,VALUE,DENS,NP,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,NORM,DOUBLE
C
      DIMENSION XYZ(3,NP),VALUE(NP),DENS(L2)
      DIMENSION IX(84),IY(84),IZ(84),JX(84),JY(84),JZ(84),
     *          XIN(343),YIN(343),ZIN(343),
     *          IJX(784),IJY(784),IJZ(784),DIJ(784),WINT(784),WORK(784),
     *          HP(28),WP(28),MINP(7),MAXP(7)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      EQUIVALENCE (HP(1),H1),(WP(1),W1)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (PI212 =1.1283791670955D+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 MINP /1,2,4, 7,11,16,22/
      DATA MAXP /1,3,6,10,15,21,28/
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          EVALUATE THE MOLECULAR ELECTROSTATIC POTENTIAL
C       AT ALL -NP- POINTS GIVEN IN -XYZ- INTO ARRAY -VALUE-
C
      TOL  = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
C     ----- INITIALIZE PARALLEL -----
C     IT IS POSSIBLE THAT THIS ROUTINE IS CALLED FOR ONLY A FEW
C     POINTS, SUCH AS THE NUCLEI, BUT USUALLY IT IS CALLED FOR
C     AN ENTIRE 2D GRID OR SURFACE CAVITY OF UNITED SPHERES,
C     SO THE PARALLELIZATION IS OVER THE POINTS BEING EVALUATED.
C
      IPCOUNT = ME - 1
C
C           LOOP OVER PROPERTY COORDINATES (POINTS)
C
      CALL VCLR(VALUE,1,NP)
      DO 600 IPT=1,NP
C
C     ----- GO PARALLEL! (STATIC LOAD BALANCING) -----
C
      IF (GOPARR) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 600
      END IF
C
      XPP = XYZ(1,IPT)
      YPP = XYZ(2,IPT)
      ZPP = XYZ(3,IPT)
C
C                    LOOP OVER SHELLS II
C
      DO 510 II=1,NSHELL
      I    = KATOM(II)
      XI   = C(1,I)
      YI   = C(2,I)
      ZI   = C(3,I)
      I1   = KSTART(II)
      I2   = I1 + KNG(II) - 1
      LIT  = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II) - MINI
C
C                    LOOP OVER SHELLS JJ
C
      DO 500  JJ=1,II
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
      IANDJ  = II.EQ.JJ
      RR     = (XI-XJ)**2 + (YI-YJ)**2 + (ZI-ZJ)**2
      CALL VCLR(WORK,1,784)
C
C         FOR SPDFGHI BASIS, LIT AND LJT MAY RANGE UP TO 7
C
      NROOTS = (LIT + LJT - 2)/2 + 1
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
      FI  = PI212*AA1
C
      AAX = (AI*XI + AJ*XJ)
      AAY = (AI*YI + AJ*YJ)
      AAZ = (AI*ZI + AJ*ZJ)
C
      AX  = AAX*AA1
      AY  = AAY*AA1
      AZ  = AAZ*AA1
C
      DUM = AI*AJ*RR*AA1
      IF(DUM.GT.TOL) GO TO 400
      FAC = FI*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 200 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 180 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
  180    CONTINUE
  200 CONTINUE
C
C       CALCULATE POINTS AND WEIGHTS FOR RYS POLYNOMIAL
C
      XX  = AA * ((AX-XPP)**2 + (AY-YPP)**2 + (AZ-ZPP)**2)
      IF (NROOTS.LE.3) CALL RT123
      IF (NROOTS.EQ.4) CALL ROOT4
      IF (NROOTS.EQ.5) CALL ROOT5
      IF (NROOTS.GE.6) CALL ROOT6
C
C       LOOP OVER ROOTS OF RYS POLYNOMIAL TO CALCULATE INTEGRALS
C
      MM = 0
      DO 340  K=1,NROOTS
C
      UU = AA*U(K)
      WW = W(K)
      TT = ONE/(AA+UU)
      T  = SQRT(TT)
C
      X0 = (AAX + UU*XPP)*TT
      Y0 = (AAY + UU*YPP)*TT
      Z0 = (AAZ + UU*ZPP)*TT
C
C      CALCULATE 1-DIMENSIONAL INTEGRALS OVER ALL ANGULAR MOMENTA
C
      IN = -7+MM
      DO 320 I=1,LIT
      IN = IN+7
      NI = I
C
      DO 310 J=1,LJT
      JN = IN+J
      NJ = J
C
C       EVALUATE MOMENT INTEGRALS USING GAUSS-HERMITE QUADRATURE
C
      XINT0 = ZERO
      YINT0 = ZERO
      ZINT0 = ZERO
C
C       FOR SPDFGHI BASIS, -NI- AND -NJ- RANGE UP TO 7 EACH
C
      NPTS = (NI + NJ - 2)/2 + 1
      IMIN = MINP(NPTS)
      IMAX = MAXP(NPTS)
C
      DO 300  IROOT=IMIN,IMAX
C
      DUM = WP(IROOT)
      PX = DUM
      PY = DUM
      PZ = DUM
C
      DUM = HP(IROOT)*T
      PTX = DUM + X0
      PTY = DUM + Y0
      PTZ = DUM + Z0
C
      AXI = PTX - XI
      AYI = PTY - YI
      AZI = PTZ - ZI
C
      BXI = PTX - XJ
      BYI = PTY - YJ
      BZI = PTZ - ZJ
C
      GO TO (101,102,103,104,105,106,107), NI
C
  107 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  106 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  105 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  104 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  103 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  102 PX = PX*AXI
      PY = PY*AYI
      PZ = PZ*AZI
C
  101 CONTINUE
C
      GO TO (201,202,203,204,205,206,207), NJ
C
  207 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  206 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  205 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  204 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  203 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  202 PX = PX*BXI
      PY = PY*BYI
      PZ = PZ*BZI
C
  201 CONTINUE
C
      XINT0 = XINT0 + PX
      YINT0 = YINT0 + PY
      ZINT0 = ZINT0 + PZ
C
C        CLOSES LOOP OVER GAUSS-HERMITE NODES
  300 CONTINUE
C
      XIN(JN) = XINT0
      YIN(JN) = YINT0
      ZIN(JN) = ZINT0*WW
C
C        THESE TWO CLOSE LOOPS OVER ANGULAR MOMENTA
  310 CONTINUE
  320 CONTINUE
C
      MM = MM+49
C        CLOSES LOOP OVER RYS NODES
  340 CONTINUE
C
C                      LOOP OVER ORBITAL PRODUCTS
C
      DO 360 I=1,IJ
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         SUM = ZERO
         MM = 0
         DO 350 K=1,NROOTS
            SUM = SUM + XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
            MM = MM+49
  350    CONTINUE
         WINT(I) = SUM*DIJ(I)
  360 CONTINUE
C
C              SET UP EXPECTATION VALUE MATRICES
C
      INDEX = 0
      MAX = MAXJ
      NN2 = 0
      DO 380 I=MINI,MAXI
         IF (IANDJ) MAX = I
         DO 370 J=MINJ,MAX
            INDEX = INDEX + 1
            NN2 = NN2 + 1
            WORK(INDEX) = WORK(INDEX) + WINT(NN2)
  370    CONTINUE
  380 CONTINUE
C
C        THESE TWO CLOSE LOOPS OVER PRIMITIVES
  400 CONTINUE
  410 CONTINUE
C
      INDEX = 0
      VT = ZERO
      MAX = MAXJ
      DO 420  I=MINI,MAXI
         LI = LOCI + I
         IN = LI*(LI-1)/2
         IF (IANDJ) MAX = I
         DO 420  J=MINJ,MAX
            LJ = LOCJ + J
            JN = LJ + IN
            INDEX = INDEX + 1
            DW  = DENS(JN)*WORK(INDEX)
            IF(LI .NE. LJ) DW = DW + DW
            VT = VT + DW
  420 CONTINUE
      VALUE(IPT) = VALUE(IPT) + VT
C
C        THESE TWO END LOOPS OVER SHELLS
  500 CONTINUE
  510 CONTINUE
C
C        END OF LOOP OVER POINTS
  600 CONTINUE
C
      IF(GOPARR) CALL DDI_GSUMF(917,VALUE,NP)
C
      RETURN
      END
C*MODULE PRPEL   *DECK INTDEN
      SUBROUTINE INTDEN(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DIJ(*),WINT(*),IJX(*),IJY(*),IJZ(*)
      DIMENSION XIN(49),YIN(49),ZIN(49)
      COMMON /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
C
C     ----- CALCULATE INTEGRALS FOR ELECTRON DENSITY -----
C
C     CALCULATE 1-DIMENSIONAL INTEGRALS OVER ALL ANGULAR MOMENTA
C
      IN = -7
      DO 110 I=1,LIT
         IN = IN+7
         NI = I
         DO 100  J=1,LJT
            JN = IN+J
            NJ = J
            CALL DELINT
            XIN(JN) = XINT1
            YIN(JN) = YINT1
            ZIN(JN) = ZINT1
  100    CONTINUE
  110 CONTINUE
C
C                     FORM 3-DIMENSIONAL INTEGRALS
C
      DO 120 I = 1,IJ
         NX = IJX(I)
         NY = IJY(I)
         NZ = IJZ(I)
         WINT(I) = WINT(I) + DIJ(I)*XIN(NX)*YIN(NY)*ZIN(NZ)
  120 CONTINUE
      RETURN
      END
C*MODULE PRPEL   *DECK INTEFL
      SUBROUTINE INTEFL(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT,AA,AX,AY,AZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DIJ(*),WINT(*),IJX(*),IJY(*),IJZ(*)
      DIMENSION XIN(784,3),YIN(784,3),ZIN(784,3)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      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
      COMMON /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, THREE=3.0D+00)
      PARAMETER (PI212=1.1283791670955D+00)
C
C       CALCULATE INTEGRALS FOR ELECTROSTATIC FIELD AND GRADIENT.
C
      FAC = TWO*PI212/AA
      TWTHRD = TWO/THREE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
C           CALCULATE POINTS AND WEIGHTS OF RYS POLYNOMIAL
C
      XX  = AA*((AX-XP)**2 + (AY-YP)**2 + (AZ-ZP)**2)
      NROOTS = (LIT + LJT)/2 + 1
      IF(NROOTS.LE.3) CALL RT123
      IF(NROOTS.EQ.4) CALL ROOT4
      IF(NROOTS.EQ.5) CALL ROOT5
      IF(NROOTS.GE.6) CALL ROOT6
C
      NM = IEFLD
C
C                  LOOP OVER ROOTS OF RYS POLYNOMIAL
C
      MM = 0
      DO 120 K=1,NROOTS
C
      UU = AA*U(K)
      WW = W(K)
      TT = ONE/(AA+UU)
      T  = SQRT(TT)
C
      X0 = (AAX + UU*XP)*TT
      Y0 = (AAY + UU*YP)*TT
      Z0 = (AAZ + UU*ZP)*TT
C
C      CALCULATE 1-DIMENSIONAL INTEGRALS OVER ALL ANGULAR MOMENTA
C
      IN = -7+MM
      DO 110 I=1,LIT
      IN = IN+7
      NI = I
C
      DO 100 J=1,LJT
      JN = IN+J
      NJ = J
C
C        GAUSS-HERMITE QUADRATURE HERE FOR MOMENT INTEGRALS
C
      CALL MOMINT
C
      XIN(JN,1) = XINT0
      YIN(JN,1) = YINT0
      ZIN(JN,1) = ZINT0*WW
C
      XIN(JN,2) = XINT1
      YIN(JN,2) = YINT1
      ZIN(JN,2) = ZINT1*WW
C
      IF(NM.EQ.2) THEN
         XIN(JN,3) = XINT2
         YIN(JN,3) = YINT2
         ZIN(JN,3) = ZINT2*WW
      END IF
C
  100 CONTINUE
  110 CONTINUE
C
      MM = MM + 49
  120 CONTINUE
C
C                      LOOP OVER ORBITAL PRODUCTS
C
      DO 180  I=1,IJ
C
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
C
      SUMX  = ZERO
      SUMY  = ZERO
      SUMZ  = ZERO
      SUMXX = ZERO
      SUMYY = ZERO
      SUMZZ = ZERO
      SUMXY = ZERO
      SUMXZ = ZERO
      SUMYZ = ZERO
C
      MM = 0
      DO 160 K=1,NROOTS
C
      INDX = NX + MM
      INDY = NY + MM
      INDZ = NZ + MM
C
      XIN0 = XIN(INDX,1)
      YIN0 = YIN(INDY,1)
      ZIN0 = ZIN(INDZ,1)
C
      XIN1 = XIN(INDX,2)
      YIN1 = YIN(INDY,2)
      ZIN1 = ZIN(INDZ,2)
C
      UU  = AA*U(K)
C
      SUMX = SUMX + XIN1 * YIN0 * ZIN0 * UU
      SUMY = SUMY + XIN0 * YIN1 * ZIN0 * UU
      SUMZ = SUMZ + XIN0 * YIN0 * ZIN1 * UU
C
      IF(NM.EQ.2) THEN
C
         XIN2 = XIN(INDX,3)
         YIN2 = YIN(INDY,3)
         ZIN2 = ZIN(INDZ,3)
C
         UUU = UU*UU
C
         SUMXX = SUMXX + XIN2 * YIN0 * ZIN0 * UUU
         SUMYY = SUMYY + XIN0 * YIN2 * ZIN0 * UUU
         SUMZZ = SUMZZ + XIN0 * YIN0 * ZIN2 * UUU
         SUMXY = SUMXY + XIN1 * YIN1 * ZIN0 * UUU
         SUMXZ = SUMXZ + XIN1 * YIN0 * ZIN1 * UUU
         SUMYZ = SUMYZ + XIN0 * YIN1 * ZIN1 * UUU
      END IF
C
      MM = MM + 49
  160 CONTINUE
C
      DEN1 = FAC*DIJ(I)
      DEN2 = DEN1*TWTHRD
C
      INDEX = I
      WINT(INDEX) = WINT(INDEX) + SUMX*DEN1
      INDEX = INDEX + IJ
      WINT(INDEX) = WINT(INDEX) + SUMY*DEN1
      INDEX = INDEX + IJ
      WINT(INDEX) = WINT(INDEX) + SUMZ*DEN1
C
      IF(NM.EQ.2) THEN
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMXX*DEN2
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMYY*DEN2
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMZZ*DEN2
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMXY*DEN2
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMXZ*DEN2
         INDEX = INDEX + IJ
         WINT(INDEX) = WINT(INDEX) + SUMYZ*DEN2
      END IF
  180 CONTINUE
C
      RETURN
      END
C
C*MODULE PRPEL   *DECK INTMOM
      SUBROUTINE INTMOM(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT,AA,AX,AY,AZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DIJ(*),WINT(*),IJX(*),IJY(*),IJZ(*)
      DIMENSION XIN(49,4),YIN(49,4),ZIN(49,4)
C
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      COMMON /XYZORB/ T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ,NM
C
      PARAMETER (ONE=1.0D+00)
C
C       CALCULATE INTEGRALS FOR ELECTROSTATIC MOMENTS.
C       -NM- 1,2,3,4 MEANS MONOPOLE, DIPOLE, QUADRUPOLE, OCTUPOLE
C       THIS ROUTINE SHOULDN'T BE CALLED FOR LESS THAN DIPOLE.
C
      NM = IEMOM
C
      T  = ONE/SQRT(AA)
      X0 = AX
      Y0 = AY
      Z0 = AZ
C
C      CALCULATE 1-DIMENSIONAL INTEGRALS OVER ALL ANGULAR MOMENTA
C
      IN = -7
      DO 100  I=1,LIT
      IN = IN+7
      NI = I
C
      DO 100  J=1,LJT
      JN = IN+J
      NJ = J
C
      CALL MOMINT
C
      XIN(JN,1) = XINT0*T
      YIN(JN,1) = YINT0*T
      ZIN(JN,1) = ZINT0*T
C
      XIN(JN,2) = XINT1*T
      YIN(JN,2) = YINT1*T
      ZIN(JN,2) = ZINT1*T
C
      IF(NM.GE.2) THEN
         XIN(JN,3) = XINT2*T
         YIN(JN,3) = YINT2*T
         ZIN(JN,3) = ZINT2*T
      END IF
C
      IF(NM.GE.3) THEN
         XIN(JN,4) = XINT3*T
         YIN(JN,4) = YINT3*T
         ZIN(JN,4) = ZINT3*T
      END IF
C
  100 CONTINUE
C
C                      LOOP OVER ORBITAL PRODUCTS
C
      DO 180 I=1,IJ
C
      NX = IJX(I)
      NY = IJY(I)
      NZ = IJZ(I)
      DENS = DIJ(I)
C
      XIN0 = XIN(NX,1)
      YIN0 = YIN(NY,1)
      ZIN0 = ZIN(NZ,1)
C
      XIN1 = XIN(NX,2)
      YIN1 = YIN(NY,2)
      ZIN1 = ZIN(NZ,2)
C
      INDEX = I
      WINT(INDEX) = WINT(INDEX) + XIN1 * YIN0 * ZIN0 * DENS
      INDEX = INDEX + IJ
      WINT(INDEX) = WINT(INDEX) + XIN0 * YIN1 * ZIN0 * DENS
      INDEX = INDEX + IJ
      WINT(INDEX) = WINT(INDEX) + XIN0 * YIN0 * ZIN1 * DENS
C
      IF(NM.GE.2) THEN
C
        XIN2 = XIN(NX,3)
        YIN2 = YIN(NY,3)
        ZIN2 = ZIN(NZ,3)
C
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN2 * YIN0 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN2 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN0 * ZIN2 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN1 * YIN1 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN1 * YIN0 * ZIN1 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN1 * ZIN1 * DENS
      END IF
C
      IF(NM.GE.3) THEN
C
        XIN3 = XIN(NX,4)
        YIN3 = YIN(NY,4)
        ZIN3 = ZIN(NZ,4)
C
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN3 * YIN0 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN3 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN0 * ZIN3 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN2 * YIN1 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN2 * YIN0 * ZIN1 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN1 * YIN2 * ZIN0 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN2 * ZIN1 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN1 * YIN0 * ZIN2 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN0 * YIN1 * ZIN2 * DENS
        INDEX = INDEX + IJ
        WINT(INDEX) = WINT(INDEX) + XIN1 * YIN1 * ZIN1 * DENS
      END IF
C
  180 CONTINUE
C
      RETURN
      END
C*MODULE PRPEL   *DECK MOMINT
      SUBROUTINE MOMINT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      DIMENSION H(36),W(36),MIN(8),MAX(8)
C
      COMMON /PRPINT/ XINT0,XINT1,XINT2,XINT3,
     *                YINT0,YINT1,YINT2,YINT3,
     *                ZINT0,ZINT1,ZINT2,ZINT3
      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
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      EQUIVALENCE (H(1),H1),(W(1),W1)
C
      DATA MIN /1,2,4, 7,11,16,22,29/
      DATA MAX /1,3,6,10,15,21,28,36/
C
C       EVALUATE MOMENT INTEGRALS USING GAUSS-HERMITE QUADRATURE
C
      XINT0 = ZERO
      YINT0 = ZERO
      ZINT0 = ZERO
      XINT1 = ZERO
      YINT1 = ZERO
      ZINT1 = ZERO
      XINT2 = ZERO
      YINT2 = ZERO
      ZINT2 = ZERO
      XINT3 = ZERO
      YINT3 = ZERO
      ZINT3 = ZERO
C
C         -NI- AND -NJ- RUN UP TO 7 FOR SPDFGHI BASIS SETS,
C         -NM- IS BETWEEN 0 AND 2 FOR EP, ELEC. FIELD, E.F. GRADIENT,
C         -NM- IS BETWEEN 0 AND 3 FOR MOMENTS UP TO OCTUPOLE,
C         THEREFORE -NPTS- CANNOT EXCEED 8
C
      NPTS = (NI + NJ + NM - 2)/2 + 1
      IMIN = MIN(NPTS)
      IMAX = MAX(NPTS)
C
      DO 500 I = IMIN,IMAX
         DUM = W(I)
         PX = DUM
         PY = DUM
         PZ = DUM
C
         DUM = H(I)*T
         PTX = DUM + X0
         PTY = DUM + Y0
         PTZ = DUM + Z0
C
         AX = PTX - XI
         AY = PTY - YI
         AZ = PTZ - ZI
C
         BX = PTX - XJ
         BY = PTY - YJ
         BZ = PTZ - ZJ
C
         CX = PTX - XP
         CY = PTY - YP
         CZ = PTZ - ZP
C
         GO TO (101,102,103,104,105,106,107), NI
C
  107    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  106    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  105    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  104    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  103    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  102    PX = PX*AX
         PY = PY*AY
         PZ = PZ*AZ
C
  101    CONTINUE
C
         GO TO (201,202,203,204,205,206,207), NJ
C
  207    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  206    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  205    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  204    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  203    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  202    PX = PX*BX
         PY = PY*BY
         PZ = PZ*BZ
C
  201    CONTINUE
C
         NM1 = NM + 1
         GO TO (301,302,303,304),NM1
C
  304    XINT3 = XINT3 + PX*CX*CX*CX
         YINT3 = YINT3 + PY*CY*CY*CY
         ZINT3 = ZINT3 + PZ*CZ*CZ*CZ
C
  303    XINT2 = XINT2 + PX*CX*CX
         YINT2 = YINT2 + PY*CY*CY
         ZINT2 = ZINT2 + PZ*CZ*CZ
C
  302    XINT1 = XINT1 + PX*CX
         YINT1 = YINT1 + PY*CY
         ZINT1 = ZINT1 + PZ*CZ
C
  301    XINT0 = XINT0 + PX
         YINT0 = YINT0 + PY
         ZINT0 = ZINT0 + PZ
C
  500 CONTINUE
      RETURN
      END
C*MODULE PRPEL   *DECK PDCCAL
      SUBROUTINE PDCCAL(PTS,POT,A,Q,IPVT,WORK,NPTS,NEQ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION PTS(3,NPTS),POT(NPTS),
     *          A(NEQ,NEQ),Q(NEQ),IPVT(NEQ),WORK(NEQ)
      DIMENSION STATD(3),STATQ(6),KD(3),KQ(6),
     *          SVCENT(3),SVDIP(3),SVQUAD(6)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /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 /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /XYZPRP/ XP,YP,ZP,
     *                CDMX,CDMY,CDMZ,
     *                CQXX,CQYY,CQZZ,CQXY,CQXZ,CQYZ,
     *                CQMXX,CQMYY,CQMZZ,CQMXY,CQMXZ,CQMYZ,
     *                COXXX,COXXY,COXXZ,COXYY,COYYY,COYYZ,
     *                COXZZ,COYZZ,COZZZ,COXYZ,
     *                COMXXX,COMXXY,COMXXZ,COMXYY,COMYYY,
     *                COMYYZ,COMXZZ,COMYZZ,COMZZZ,COMXYZ
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, HALF=0.5D+00,
     *           PT75=0.75D+00, ONEPT5=1.5D+00, TWO=2.0D+00,
     *           TENP2=1.0D+02, TENM10=1.0D-10)
C
C        DEBYE = CONVERSION FROM DEBYES TO AU
C        BUCK  = CONVERSION FROM BUCKINGHAMS TO AU
C
      PARAMETER (DEBYE=ONE/2.541766D+00, BUCK=ONE/1.345044D+00,
     *           AU2CAL=627.51D+00)
C
      DATA DIPOLE/8HDIPOLE  /,QUPOLE/8HQUPOLE  /,
     *     FIXED /8HFIXED   /, FREE/8HFREE    /
C
C     THIS PROGRAM GENERATES ATOM CENTRED MONOPOLES
C     (NET ATOMIC CHARGES) BY A LEAST SQUARES FITTING OF
C     ELECTROSTATIC POTENTIALS AT SLECTED POINTS AROUND THE
C     MOLECULE, USING LAGRANGE MULTIPLIERS WITH THE
C     FOLLOWING OPTIONS FOR CONSTRAINTS:
C
C     CONSTR=NONE    ---  PRINT ELECTROSTATIC POTENTIAL AT POINTS BUT DO
C                         NOT PERFORM A FIT
C     CONSTR=CHARGE  ---  SUM OF CHARGES = TOTAL MOLECULAR CHARGE (DEFAU
C
C     CONSTR=DIPOLE  ---  CALCULATED DIPOLE COMPONENTS = INPUTTED VALUES
C                         (ALSO CONSTRAINS CHARGE)
C
C     CONSTR=QUPOLE  ---  CALCULATED QUADRUPOLE COMPONENTS = INPUTTED VA
C                         (ALSO CONSTRAINS DIPOLE AND CHARGE)
C
      CALL DCOPY(3,CENTER,1,SVCENT,1)
      CALL DCOPY(3,DPOLE ,1,SVDIP ,1)
      CALL DCOPY(6,QPOLE ,1,SVQUAD,1)
C
      TEST = DDOT(3,CENTER,1,CENTER,1)
      IF(TEST.LT.1.0D-02) THEN
         CENTER(1) = XP
         CENTER(2) = YP
         CENTER(3) = ZP
      END IF
C
      TEST = DDOT(3,DPOLE,1,DPOLE,1)
      IF(TEST.LT.1.0D-02) THEN
         DPOLE(1) = CDMX * DEBYE
         DPOLE(2) = CDMY * DEBYE
         DPOLE(3) = CDMZ * DEBYE
      END IF
C
      TEST = DDOT(6,QPOLE,1,QPOLE,1)
      IF(TEST.LT.1.0D-02) THEN
         QPOLE(1) = CQXX * BUCK
         QPOLE(2) = CQYY * BUCK
         QPOLE(3) = CQZZ * BUCK
         QPOLE(4) = CQXY * BUCK
         QPOLE(5) = CQXZ * BUCK
         QPOLE(6) = CQYZ * BUCK
      END IF
C
      NDIP  = 0
      NQUAD = 0
      STATC = FIXED
      DO 303,I=1,3
         KD(I)=0
         STATD(I)=FREE
 303  CONTINUE
      DO 302,I=1,6
         KQ(I)=0
         STATQ(I)=FREE
 302  CONTINUE
C
C     TEST FOR ALL X (OR Y OR Z) COORDINATES ZERO
C     IN WHICH CASE THE RELEVANT DIPOLE OR QUADRUPOLE MOMENT
C     IS NECESSARILY ZERO AND NEED NOT (CANNOT!) BE CONSTRAINED
C
      IF(CONSTR.EQ.DIPOLE.OR.CONSTR.EQ.QUPOLE) THEN
         DO 30 I=1,NAT
            DO 20 J=1,3
               IF(ABS(C(J,I)).GT.TENM10) KD(J)=1
   20       CONTINUE
            XREL=C(1,I)-CENTER(1)
            YREL=C(2,I)-CENTER(2)
            ZREL=C(3,1)-CENTER(3)
            XX=XREL*XREL-HALF*(YREL*YREL-ZREL*ZREL)
            YY=YREL*YREL-HALF*(ZREL*ZREL-XREL*XREL)
            ZZ=ZREL*ZREL-HALF*(XREL*XREL-YREL*YREL)
            XY=XREL*YREL
            XZ=XREL*ZREL
            YZ=YREL*ZREL
            IF(KQ(1).EQ.0.AND.ABS(XX).GT.TENM10) KQ(1)=1
            IF(KQ(2).EQ.0.AND.ABS(YY).GT.TENM10) KQ(2)=1
            IF(KQ(3).EQ.0.AND.ABS(ZZ).GT.TENM10) KQ(3)=1
            IF(KQ(4).EQ.0.AND.ABS(XY).GT.TENM10) KQ(4)=1
            IF(KQ(5).EQ.0.AND.ABS(XZ).GT.TENM10) KQ(5)=1
            IF(KQ(6).EQ.0.AND.ABS(YZ).GT.TENM10) KQ(6)=1
   30    CONTINUE
C
         DO 40 I=1,3
            IF(KD(I).NE.0) THEN
               NDIP=NDIP+1
               STATD(I)=FIXED
            END IF
   40    CONTINUE
      END IF
C
      IF(CONSTR.EQ.QUPOLE) THEN
         DO 50 I=1,6
            IF(KQ(1)+KQ(2)+KQ(3).EQ.3) KQ(1)=0
            IF(KQ(I).NE.0) THEN
               NQUAD=NQUAD+1
               STATQ(I)=FIXED
            END IF
   50    CONTINUE
      END IF
      NCNSTR=NDIP+NQUAD+1
C
C     SEE IF THERE ARE TOO MANY CONSTRAINTS
C
      IF(NCNSTR.GT.NAT) THEN
         IF(MASWRK) WRITE(IW,*)
     *      'NUMBER OF CONSTRAINTS CANNOT EXCEED NUMBER OF ATOMS'
         CALL ABRT
      END IF
C
C     SET UP MATRIX OF LINEAR COEFFICIENTS, A
C     BEGIN LOOP OVER ROWS
C
      DO 100 K=1,NAT
C
C     BEGIN LOOP OVER COLUMNS
C
      DO 200 MU=1,NAT
         SUM = ZERO
         DO 400 I=1,NPTS
            XPT=PTS(1,I)
            YPT=PTS(2,I)
            ZPT=PTS(3,I)
            RIK = (XPT-C(1,K))**2 + (YPT-C(2,K))**2 + (ZPT-C(3,K))**2
            RIK = SQRT(RIK)
            RIMU= (XPT-C(1,MU))**2 +(YPT-C(2,MU))**2 +(ZPT-C(3,MU))**2
            RIMU= SQRT(RIMU)
            SUM = SUM + ONE / (RIK * RIMU)
  400    CONTINUE
         A(K,MU) = SUM
  200 CONTINUE
  100 CONTINUE
C
      DO 500 I=1,NAT
         NFREE=1
C
C     FILL OUT COLUMNS CORRESPONDING TO LAGRANGE MULTIPLIERS
C     CHARGE FIRST:
C
         A(I,NAT+NFREE) = HALF
C
C     DX,DY,DZ
C
         IF(CONSTR.EQ.DIPOLE.OR.CONSTR.EQ.QUPOLE) THEN
            DO 60 M=1,3
               IF(KD(M).NE.ZERO) THEN
                  NFREE=NFREE+1
                  A(I,NAT+NFREE) = HALF*C(M,I)
               END IF
   60       CONTINUE
         END IF
C
C     QYY,QZZ,QXY,QXZ,QYZ   (N.B. SINCE QXX+QYY+QZZ=0,
C                            QXX IS NOT USED AS A CONSTRAINT)
C
         IF(CONSTR.EQ.QUPOLE) THEN
         XREL=C(1,I)-CENTER(1)
         YREL=C(2,I)-CENTER(2)
         ZREL=C(3,I)-CENTER(3)
         IF(KQ(1).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = HALF*(XREL*XREL-HALF*(YREL*YREL+ZREL*ZREL))
         END IF
         IF(KQ(2).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = HALF*(YREL*YREL-HALF*(ZREL*ZREL+XREL*XREL))
         END IF
         IF(KQ(3).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = HALF*(ZREL*ZREL-HALF*(XREL*XREL+YREL*YREL))
         END IF
         IF(KQ(4).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = PT75*(XREL*YREL)
         END IF
         IF(KQ(5).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = PT75*(XREL*ZREL)
         END IF
         IF(KQ(6).NE.ZERO) THEN
         NFREE=NFREE+1
         A(I,NAT+NFREE) = PT75*(YREL*ZREL)
         END IF
         END IF
C
C     FILL OUT THE ROWS CORRESPONDING TO CONSTRAINTS
C
         DO 995,M=1,NCNSTR
            A(NAT+M,I) = TWO*A(I,NAT+M)
  995    CONTINUE
  500 CONTINUE
C
C     FILL OUT THE BLOCK WHICH CONNECTS LAGRANGE MULTIPLIERS TO
C     CONSTRAINTS
C
      DO 600 K=NAT+1,NAT+NCNSTR
      DO 600 MU=NAT+1,NAT+NCNSTR
      A(K,MU) = ZERO
  600 CONTINUE
C
C     CONSTRUCT COLUMN VECTOR, Q
C
      DO 700 K=1,NAT
         SUM = ZERO
         DO 710 I=1,NPTS
            XPT = PTS(1,I)
            YPT = PTS(2,I)
            ZPT = PTS(3,I)
            RIK = (XPT-C(1,K))**2 + (YPT-C(2,K))**2 + (ZPT-C(3,K))**2
            RIK = SQRT(RIK)
            SUM = SUM + POT(I) / RIK
  710    CONTINUE
         Q(K) = SUM
  700 CONTINUE
C
C     CONSTRUCT THE PORTION CORRESPONDING TO LAGRANGE MULTIPLIER
C     IF THEY WERE NOT INPUT, TAKE CALCULATED DIPOLE AND QUADRUPOLE
C
      NFREE=1
      Q(NAT+NFREE) = ICH
C
      IF(CONSTR.EQ.DIPOLE  .OR.  CONSTR.EQ.QUPOLE) THEN
         DO 994 M=1,3
            IF(KD(M).NE.ZERO) THEN
               NFREE=NFREE+1
               Q(NAT+NFREE) = DPOLE(M)
            END IF
  994    CONTINUE
      END IF
C
      IF(CONSTR.EQ.QUPOLE) THEN
         DO 992 M=1,6
            IF(KQ(M).NE.ZERO) THEN
               NFREE=NFREE+1
               Q(NAT+NFREE) = QPOLE(M)
            END IF
  992    CONTINUE
      END IF
C
C     SOLVE MATRIX EQUATION AX = Q;
C     WHERE X = (Q1,Q2, ... QN,L1,L2,L3,L4)
C
      NSYST = NAT+NCNSTR
      CALL SLVLEQ(A,Q,IPVT,NSYST,NEQ,0,IERR)
C
C     COMPUTE INVERSE MATRIX FOR ERRORS IN CHARGES
C
      CALL DGEDI(A,NEQ,NSYST,IPVT,DET,WORK,1)
C
C     COMPUTE RMS DEVIATION AND RELATIVE RMS % DEVIATION
C
      RMS = ZERO
      RRMS = ZERO
      DO 800 I=1,NPTS
         XPT = PTS(1,I)
         YPT = PTS(2,I)
         ZPT = PTS(3,I)
         EQ = ZERO
         DO 810 J=1,NAT
            DIST = (XPT-C(1,J))**2 + (YPT-C(2,J))**2 + (ZPT-C(3,J))**2
            DIST = SQRT(DIST)
            EQ = EQ + Q(J) / DIST
  810    CONTINUE
         RMS = RMS + (POT(I) - EQ)**2
         RRMS = RRMS + (POT(I)**2)
  800 CONTINUE
C
C     CONVERT RMS TO KCAL/MOL
C
      CHISQ=SQRT(RMS/(NPTS-NAT-NCNSTR))
      RMS = SQRT(RMS/NPTS)*AU2CAL
      RRMS = SQRT(RRMS/NPTS)*AU2CAL
      RRMS = RMS/RRMS*TENP2
C
C     COMPUTE STANDARD ERRORS IN CHARGES - STORE IN WORK
C
      DO 900 I=1,NAT
  900 WORK(I)=SQRT(A(I,I))*CHISQ
C
C     COMMENCE OUPUT:
C        CALCULATE THE DIPOLE AND QUADRUPOLE MOMENTS
C        FROM THE FITTED CHARGES
C
      DIPX = ZERO
      DIPY = ZERO
      DIPZ = ZERO
      QMXX = ZERO
      QMYY = ZERO
      QMZZ = ZERO
      QMYZ = ZERO
      QMXY = ZERO
      QMXZ = ZERO
      SUM = ZERO
      DO 99 I=1,NAT
         DIPX=DIPX+(Q(I)*C(1,I))
         DIPY=DIPY+(Q(I)*C(2,I))
         DIPZ=DIPZ+(Q(I)*C(3,I))
         QMXX = QMXX + (C(1,I)-CENTER(1))**2*Q(I)
         QMYY = QMYY + (C(2,I)-CENTER(2))**2*Q(I)
         QMZZ = QMZZ + (C(3,I)-CENTER(3))**2*Q(I)
         QMXY = QMXY + (C(1,I)-CENTER(1))*(C(2,I)-CENTER(2))*Q(I)
         QMYZ = QMYZ + (C(2,I)-CENTER(2))*(C(3,I)-CENTER(3))*Q(I)
         QMXZ = QMXZ + (C(1,I)-CENTER(1))*(C(3,I)-CENTER(3))*Q(I)
         SUM = SUM + Q(I)
   99 CONTINUE
C
C     CONVERT ATOMIC UNITS BACK TO DEBYE AND BUCKINGHAM (UGH)
C
      QXX = (QMXX-HALF*(QMYY+QMZZ))/BUCK
      QYY = (QMYY-HALF*(QMZZ+QMXX))/BUCK
      QZZ = (QMZZ-HALF*(QMXX+QMYY))/BUCK
      QXY = ONEPT5*QMXY / BUCK
      QXZ = ONEPT5*QMXZ / BUCK
      QYZ = ONEPT5*QMYZ / BUCK
      QPOLE(1)=QPOLE(1)/BUCK
      QPOLE(2)=QPOLE(2)/BUCK
      QPOLE(3)=QPOLE(3)/BUCK
      QPOLE(4)=QPOLE(4)/BUCK
      QPOLE(5)=QPOLE(5)/BUCK
      QPOLE(6)=QPOLE(6)/BUCK
      DIPX=DIPX / DEBYE
      DIPY=DIPY / DEBYE
      DIPZ=DIPZ / DEBYE
      DPOLE(1)=DPOLE(1)/DEBYE
      DPOLE(2)=DPOLE(2)/DEBYE
      DPOLE(3)=DPOLE(3)/DEBYE
C
C     CALCULATE TOTAL DIPOLE MOMENT
C
      DIPTOT = SQRT(DIPX**2+DIPY**2+DIPZ**2)
      DPLSCF = SQRT(DPOLE(1)**2+DPOLE(2)**2+DPOLE(3)**2)
C
C     CREATE OUTPUT
C
      IF(MASWRK) THEN
         WRITE(IW,240)
         WRITE(IW,250)
         WRITE(IW,260) (ANAM(I),BNAM(I),Q(I),WORK(I),I=1,NAT)
         WRITE(IW,280) RMS
         WRITE(IW,285) RRMS
         WRITE(IW,286)
         CHG = ICH
         WRITE(IW,288) STATC,CHG,SUM
         WRITE(IW,289) DPLSCF,DIPTOT
         WRITE(IW,290) STATD(1),DPOLE(1),DIPX
         WRITE(IW,291) STATD(2),DPOLE(2),DIPY
         WRITE(IW,292) STATD(3),DPOLE(3),DIPZ
         WRITE(IW,293) STATQ(1),QPOLE(1),QXX
         WRITE(IW,294) STATQ(2),QPOLE(2),QYY
         WRITE(IW,295) STATQ(3),QPOLE(3),QZZ
         WRITE(IW,296) STATQ(4),QPOLE(4),QXY
         WRITE(IW,297) STATQ(5),QPOLE(5),QXZ
         WRITE(IW,298) STATQ(6),QPOLE(6),QYZ
      END IF
C
      CALL DCOPY(3,SVCENT,1,CENTER,1)
      CALL DCOPY(3,SVDIP ,1,DPOLE ,1)
      CALL DCOPY(6,SVQUAD,1,QPOLE ,1)
      RETURN
C
  240 FORMAT(/1X,'NET CHARGES:'/1X,37('-'))
  250 FORMAT(1X,'ATOM                CHARGE    E.S.D.'/1X,37('-'))
  260 FORMAT(1X,2A8,2F10.4)
  280 FORMAT(1X,37('-')//1X,'RMS DEVIATION IS         ',F8.3,
     *       ' KCAL/MOLE')
  285 FORMAT(1X,'RELATIVE RMS DEVIATION IS',F8.3,' %')
  286 FORMAT(/1X,'MULTIPOLE MOMENTS:',/1X,57('-')/
     *1X,'STATUS  ACTUAL VALUE       QUANTITY        FITTED VALUE '/
     $1X,57('-'))
  288 FORMAT(1X,A8,2X,F8.4,'       TOTAL CHARGE     ',F8.4)
  289 FORMAT(1X,10X,F8.4,' D','    DIPOLE MOMENT     ',F8.4,' D')
  290 FORMAT(1X,A8,2X,F8.4,' D','          DX          ',F8.4,' D')
  291 FORMAT(1X,A8,2X,F8.4,' D','          DY          ',F8.4,' D')
  292 FORMAT(1X,A8,2X,F8.4,' D','          DZ          ',F8.4,' D')
  293 FORMAT(1X,A8,2X,F8.4,' BU','         QXX         ',F8.4,' BU')
  294 FORMAT(1X,A8,2X,F8.4,' BU','         QYY         ',F8.4,' BU')
  295 FORMAT(1X,A8,2X,F8.4,' BU','         QZZ         ',F8.4,' BU')
  296 FORMAT(1X,A8,2X,F8.4,' BU','         QXY         ',F8.4,' BU')
  297 FORMAT(1X,A8,2X,F8.4,' BU','         QXZ         ',F8.4,' BU')
  298 FORMAT(1X,A8,2X,F8.4,' BU','         QYZ         ',F8.4,' BU'/
     *       1X,57('-'))
      END
C*MODULE PRPEL   *DECK PDCPTS
      SUBROUTINE PDCPTS(PTS,RADIUS,CICOS,MAXPNT,NATM,NLAYER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      DOUBLE PRECISION MAKEFP
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 /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 PTS(3,MAXPNT),RADIUS(NATM)
      DIMENSION CICOS(3,2562)
      DIMENSION NLAYER(LAYER+1)
C
      PARAMETER (ONE=1.0D+00, TOANGS=0.52917724924D+00)
C
      DATA CHELPG /8HCHELPG  /, GEODES /8HGEODESIC/, CONLLY /8HCONNOLLY/
      DATA PI/3.1415926535898D+00/
      DATA MAKEFP/8HMAKEFP  /
C
C        ---- GENERATE POINTS FOR PDC-MEP CALCULATION -----
C
      TOBOHR = ONE/TOANGS
      NPDC=0
C
C     GET VDW RADII
C
      CALL VDWSEL(RADIUS,NATM)
C
C          ----- CHELPG SELECTION SCHEME -----
C
      IF(PTSEL.EQ.CHELPG) THEN
C
C     CONVERT RMAX, DELR TO AU
C
      DELR = DELR * TOBOHR
      IF(RUNTYP.NE.MAKEFP) RMAX = RMAX * TOBOHR
C
C     GENERATE THE POINTS USING MODIFIED CHELPG ALGORITHM
C     POINTS ARE INITIALLY SELECTED IN A CUBE AROUND THE MOLECULE WHICH
C     IS SCALED TO THE SIZE OF THE MOLECULE+RMAX. THIS IS PRESENTLY AN I
C     PARAMETER.  POINTS ARE THEN EXCLUDED IF THEY FALL WITHIN THE INPUT
C     VDW RADIUS OF ANY OF THE ATOMS, OR, IF THEY FALL OUTSIDE
C     A DESIGNATED DISTANCE (RMAX) FROM ALL OF THE ATOMS.
C
C     FIND THE EXTREMA OF THE MOLECULAR GEOMETRY.
C
      XMAX=-1.0D+05
      XMIN= 1.0D+05
      YMAX=-1.0D+05
      YMIN= 1.0D+05
      ZMAX=-1.0D+05
      ZMIN= 1.0D+05
C
C        FOR MAKEFP SCREENING, USE ATOMS AND BOND MID-
C        POINTS AS CENTERS
      IF(RUNTYP.EQ.MAKEFP) THEN
         DO 39 I=1,NPTSTN
            IF(STNPNT(2,I).GT.XMAX) XMAX = STNPNT(2,I)
            IF(STNPNT(2,I).LT.XMIN) XMIN = STNPNT(2,I)
            IF(STNPNT(3,I).GT.YMAX) YMAX = STNPNT(3,I)
            IF(STNPNT(3,I).LT.YMIN) YMIN = STNPNT(3,I)
            IF(STNPNT(4,I).GT.ZMAX) ZMAX = STNPNT(4,I)
            IF(STNPNT(4,I).LT.ZMIN) ZMIN = STNPNT(4,I)
   39    CONTINUE
      ELSE
         DO 40 I=1,NAT
            IF(C(1,I).GT.XMAX) XMAX = C(1,I)
            IF(C(1,I).LT.XMIN) XMIN = C(1,I)
            IF(C(2,I).GT.YMAX) YMAX = C(2,I)
            IF(C(2,I).LT.YMIN) YMIN = C(2,I)
            IF(C(3,I).GT.ZMAX) ZMAX = C(3,I)
            IF(C(3,I).LT.ZMIN) ZMIN = C(3,I)
   40    CONTINUE
      END IF
C
C     DETERMINE THE MINIMUM CUBE DIMENSIONS REQUIRED TO CONTAIN
C     THE MOLECULE, INCLUDING THE MAXIMUM SELECTION RADIUS (RMAX)
C     ON BOTH SIDES.
C
      XRANGE = XMAX - XMIN + 2.0D+00 * RMAX
      YRANGE = YMAX - YMIN + 2.0D+00 * RMAX
      ZRANGE = ZMAX - ZMIN + 2.0D+00 * RMAX
C
      NXPTS = INT(XRANGE/DELR)
      NYPTS = INT(YRANGE/DELR)
      NZPTS = INT(ZRANGE/DELR)
      XMIN = XMIN+0.5D+00*(XRANGE-NXPTS*DELR)
      YMIN = YMIN+0.5D+00*(YRANGE-NYPTS*DELR)
      ZMIN = ZMIN+0.5D+00*(ZRANGE-NZPTS*DELR)
C
C     LOOP OVER POSSIBLE POINTS
C
      DO 50 II = 1,NXPTS + 1
      P1 = XMIN - RMAX + (II-1)*DELR
      DO 50 JJ = 1,NYPTS + 1
      P2 = YMIN - RMAX + (JJ-1)*DELR
      DO 50 KK = 1,NZPTS + 1
      P3 = ZMIN - RMAX + (KK-1)*DELR
C
C     IS THIS POINT WITHIN A VAN DER WAALS SPHERE OR OUTSIDE THE
C     RMAX DISTANCE FROM ALL ATOMS?
C
      RADMIN=50.0D+00
      IF(RUNTYP.EQ.MAKEFP) THEN
         INSIDE=0
         DO 59 I=1,NPTSTN
            VRAD = RADIUS(I)
            DIST = (P1 - STNPNT(2,I))**2 + (P2 - STNPNT(3,I))**2 +
     *           (P3 - STNPNT(4,I))**2
            DIST = SQRT(DIST)
            IF (DIST .LT. VDWSCL*VRAD) GO TO 70
            IF (DIST .LT. RMAX*VRAD) INSIDE =1
   59    CONTINUE
         IF (INSIDE .EQ. 0) GO TO 70
      ELSE
         DO 60 I=1,NAT
            VRAD = RADIUS(I)
            DIST = (P1-C(1,I))**2+(P2-C(2,I))**2+(P3-C(3,I))**2
            DIST = SQRT(DIST)
            IF (DIST .LT. VRAD) GO TO 70
            IF (DIST .LT. RADMIN) RADMIN = DIST
   60    CONTINUE
         IF (RADMIN .GT. RMAX) GO TO 70
      END IF
C
C     STORE POINTS (IN ATOMIC UNITS)
C
      NPDC = NPDC + 1
      IF(NPDC.GT.MAXPNT) THEN
         IF(MASWRK) WRITE(IW,970) MAXPNT
         CALL ABRT
      END IF
      PTS(1,NPDC) = P1
      PTS(2,NPDC) = P2
      PTS(3,NPDC) = P3
   70 CONTINUE
   50 CONTINUE
      END IF
C
C          ----- ICOSAHEDRAL GEODESIC SURFACE SELECTION SCHEME -----
C
C     GENERATE POINTS ON THE SCALED VDW SURFACE OF THE MOLECULE
C     FIRST CALL SPHERE TO GENERATE A SET OF (VERY CLOSELY) EVENLY
C     SPACED POINTS ON THE SURFACE OF A SPHERE, THEN TRANSLATE AND
C     SCALE ACCORDING TO THE ATOM TYPE.
C     NFREQ GOVERNS THE TESSELATION TYPE (SEE BELOW)
C
C          ----- CONNOLLY SURFACE POINT SELECTION SCHEME -----
C
C     THE SO-CALLED CONNOLLY SURFACE FAVOURED BY MERZ, KOLLMAN & SINGH
C     IS NOT A SMOOTHED VAN DER WAALS SURFACE AT ALL (CONTRARY TO ALL
C     IMPLICATIONS IN THEIR PUBLICATIONS).  RATHER, IT IS JUST THE
C     USUAL VDW SURFACE CONSISTING OF FUSED SPHERES - IN THIS CASE
C     THE DISTRIBUTION OF POINTS ON THE SURFACE OF THE SPHERES IS
C     DONE USING CONNOLLY'S GENUN ROUTINE.  THE RESULTS ARE IDENTICAL
C     TO THOSE FROM GAUSSIAN 9X ETC, WITH NO BEATING AROUND THE BUSH!
C
C     THE FOLLOWING CODE SUFFICES FOR BOTH GEODESIC AND "CONNOLLY"
C     SELECTION OF POINTS ON THE SURFACE, AND THE TWO CHOICES ARE
C     ALMOST INDISTINGUISHABLE (SEE SPACKMAN, J. COMP. CHEM. SUBMITTED)
C
      IF(PTSEL.EQ.GEODES.OR.PTSEL.EQ.CONLLY) THEN
C
C     FOR GEODESIC SURFACE GET TEMPLATE POINTS
C
      IF(PTSEL.EQ.GEODES) THEN
C
C     NFREQ HERE HAS THE FOLLOWING POSSIBLE MEANINGS:
C           = (10*NH + NK)        {3,5+}NH,NK TESSELATIONS
C           =-(10*NH + NK)        {5+,3}NH,NK TESSELATIONS
C     (OF COURSE BOTH NH AND NK MUST BE LESS THAN 10, SO
C     NFREQ MUST LIE WITHIN THE RANGE -99 TO 99)
C     THE DEFAULT VALUE IS NFREQ=30 (=03)
C
        ITYPE=1
        IF(NFREQ.LT.0) ITYPE=2
        NH=ABS(NFREQ)/10
        NK=MOD(ABS(NFREQ),10)
        CALL SPHGEO(ITYPE,NH,NK,CICOS,NPOINTS)
        IF(NFREQ.GT.0) WRITE(IW,910) NH,NK
        IF(NFREQ.LT.0) WRITE(IW,920) NH,NK
      END IF
      SCALE = VDWSCL
      DO 200 ILAYER = 1,LAYER
      IF(RUNTYP.EQ.MAKEFP) THEN
         NLAYER(ILAYER)=NPDC
         DO 160 K = 1,NPTSTN
         SCL = RADIUS(K)*SCALE
C
C        FOR "CONNOLLY" SURFACE GET TEMPLATE POINTS BASED ON PTDENS
C
         IF(PTSEL.EQ.CONLLY) THEN
            NPOINTS = INT(4.0D+00*PI*SCL*SCL*PTDENS)
            CALL SPHCNL(NPOINTS,CICOS)
         END IF
         DO 170 I = 1,NPOINTS
         XPT = CICOS(1,I)*SCL + STNPNT(2,K)
         YPT = CICOS(2,I)*SCL + STNPNT(3,K)
         ZPT = CICOS(3,I)*SCL + STNPNT(4,K)
C
C      DISCARD IF THIS POINT IS WITHIN SCALED VDW SPHERE OF ANOTHER ATOM
C
         DO 166 J = 1,NPTSTN
         IF(J.EQ.K) GO TO 166
         DISTJ = SQRT((XPT-STNPNT(2,J))**2+(YPT-STNPNT(3,J))**2+
     *          (ZPT-STNPNT(4,J))**2)
         IF(DISTJ.LT.RADIUS(J)*SCALE) GO TO 170
  166    CONTINUE
C
C        STORE POINTS (IN ATOMIC UNITS)
C
         NPDC = NPDC + 1
         IF(NPDC.GT.MAXPNT) THEN
            IF(MASWRK) WRITE(IW,970) MAXPNT
            CALL ABRT
         END IF
         PTS(1,NPDC) = XPT
         PTS(2,NPDC) = YPT
         PTS(3,NPDC) = ZPT
  170    CONTINUE
  160    CONTINUE
C
      ELSE
         DO 161 K = 1,NAT
         SCL = RADIUS(K)*SCALE
C
C        FOR "CONNOLLY" SURFACE GET TEMPLATE POINTS BASED ON PTDENS
C
         IF(PTSEL.EQ.CONLLY) THEN
            NPOINTS = INT(4.0D+00*PI*SCL*SCL*PTDENS)
            CALL SPHCNL(NPOINTS,CICOS)
         END IF
         DO 171 I = 1,NPOINTS
            XPT = CICOS(1,I)*SCL + C(1,K)
            YPT = CICOS(2,I)*SCL + C(2,K)
            ZPT = CICOS(3,I)*SCL + C(3,K)
C
C     DISCARD IF THIS POINT IS WITHIN SCALED VDW SPHERE OF ANOTHER ATOM
C
         DO 167 J = 1,NAT
         IF(J.EQ.K) GO TO 167
         DISTJ = SQRT((XPT-C(1,J))**2+(YPT-C(2,J))**2+(ZPT-C(3,J))**2)
         IF(DISTJ.LT.RADIUS(J)*SCALE) GO TO 171
  167    CONTINUE
C
C         STORE POINTS (IN ATOMIC UNITS)
C
         NPDC = NPDC + 1
         IF(NPDC.GT.MAXPNT) THEN
            IF(MASWRK) WRITE(IW,970) MAXPNT
            CALL ABRT
         END IF
         PTS(1,NPDC) = XPT
         PTS(2,NPDC) = YPT
         PTS(3,NPDC) = ZPT
  171    CONTINUE
  161    CONTINUE
      END IF
C
      SCALE = SCALE + VDWINC
  200 CONTINUE
      NLAYER(LAYER+1)=NPDC
      END IF
C
      IF(MASWRK) THEN
C        IF (RUNTYP.EQ.MAKEFP) WRITE(IW,*)'PDC NLAYER',NLAYER
         WRITE(IW,960) NPDC
      END IF
      RETURN
C
  910 FORMAT(' {3,5+}NH,NK TEMPLATE USED, NH=',I2,' NK=',I2)
  920 FORMAT(' {5+,3}NH,NK TEMPLATE USED, NH=',I2,' NK=',I2)
  960 FORMAT(' NUMBER OF POINTS SELECTED FOR FITTING =',I6)
  970 FORMAT(' MAXIMUM NUMBER OF POINTS MUST BE LESS THAN ',I6)
      END
C*MODULE PRPEL   *DECK PRCALC
      SUBROUTINE PRCALC(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)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      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
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (RLN10=2.30258D+00)
      PARAMETER (SQRT3 =1.73205080756887729353D+00)
      PARAMETER (SQRT5 =2.23606797749978969641D+00)
      PARAMETER (SQRT7 =2.64575131106459059050D+00)
      PARAMETER (SQRT9 =3.00000000000000000000D+00)
      PARAMETER (SQRT11=3.31662479035539984911D+00)
C
      DATA ELDEN  /8HELDEN   /,
     *     ELFLD  /8HELFLD   /,
     *     ELMOM  /8HELMOM   /
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=1,NSHELL
      I    = KATOM(II)
      XI   = C(1,I)
      YI   = C(2,I)
      ZI   = C(3,I)
      I1   = KSTART(II)
      I2   = I1 + KNG(II) - 1
      LIT  = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II) - MINI
C
C                         LOOP OVER SHELLS JJ
C
      DO 500  JJ=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
      IF(PRNAME.EQ.ELDEN)
     *  CALL INTDEN(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT)
C
      IF(PRNAME.EQ.ELFLD)
     *  CALL INTEFL(LIT,LJT,IJ,IJX,IJY,IJZ,DIJ,WINT,AA,AX,AY,AZ)
C
      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*MODULE PRPEL   *DECK COSPOT
      SUBROUTINE COSPOT(QVPOT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION NUCLEI
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MAXDEN=25*MXATM, LENABC=2000, NPPA=1082,
     *            MXPTPT=100)
C
      DIMENSION QVPOT(NPPA)
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 /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 /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
      COMMON /SOLV  / FEPSI,RDS,DISEX2,COSURF(3,LENABC),SRAD(MXATM),
     1                QDEN(MAXDEN),AR(LENABC),
     3                NSPA,NPS,NPS2,NDEN,NPSPHER
C
      PARAMETER (ZERO=0.0D+00, TM3=1.0D-03)
C
      DATA NUCLEI /8HNUCLEI  /
C
C       ----- CALCULATE ELECTROSTATIC POTENTIAL -----
C             FOR SOLUTE IN THE PRESENCE OF SOLVENT
C             AS CREATED BY COSMO, FOR CORRECTION
C
C      IN ORDER TO INCREASE THE SPEED OF CALCULATION,
C      THE LOOP OVER POINTS IS PULLED THROUGH
C      THE LOOPS OVER SHELLS AND PRIMITIVES.
C      INSTEAD OF CALLING ROUTINE "PRCALC", "EPCALC" IS CALLED.
C
C
C               SET POINTERS FOR PARTITIONING OF MEMORY
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      NPOINT=NPS
      MAXPNT=NPOINT
C
      CALL VALFM(LOADFM)
C
C     ALLOCATE MEMORY, NOTE -IXYZ- NEEDS TO BE FIRST IF -PDC- USED
C
      MAXBND = 5*NAT
C
      IXYZ   = LOADFM + 1
      IDENSA = IXYZ   + 3*MAXPNT
      IDENSB = IDENSA + L2
      IPOT   = IDENSB + L2
      IPAIRS = IPOT   + MAXPNT
      LAST   = IPAIRS + MAXBND
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
C        GET TOTAL DENSITY MATRIX
C
      CALL DENDD1(X(IDENSA),X(IDENSB),L2)
C
C              SET UP ALL POINTS TO BE EVALUATED
C
         IPOINT = 0
         INDEX  = IXYZ - 3
  140    CONTINUE
         IPOINT = IPOINT + 1
         XPT = ZERO
         YPT = ZERO
         ZPT = ZERO
         ISTAT=0
         CALL COSPTS(IPOINT,XPT,YPT,ZPT,ISTAT)
         IF(ISTAT.LT.0) CALL ABRT
         IF(ISTAT.GT.0) GO TO 170
         INDEX = INDEX + 3
         X(INDEX)   = XPT
         X(INDEX+1) = YPT
         X(INDEX+2) = ZPT
         GO TO 140
C
C        CALCULATE ELECTRONIC CONTRIBUTION AT ALL POINTS
C
  170 CONTINUE
      CALL EPCALC(X(IXYZ),X(IPOT),X(IDENSA),MAXPNT,L2)
C
C        OUTPUT RESULTS OF ALL POINTS
C
      INDEX  = IXYZ - 3
      DO 240  IPT=1,MAXPNT
         ELPOTE = -X(IPOT+IPT-1)
         INDEX  = INDEX + 3
         XPT = X(INDEX)
         YPT = X(INDEX+1)
         ZPT = X(INDEX+2)
C                 FORM NUCLEAR CONTRIBUTION
         ELPOTN = ZERO
         DO 210 I = 1,NAT
            XN = C(1,I) - XPT
            YN = C(2,I) - YPT
            ZN = C(3,I) - ZPT
            RR = SQRT(XN*XN + YN*YN + ZN*ZN)
            IF(RR.LT.TM3  .AND.  ELPLOC.NE.NUCLEI .AND. MASWRK)
     *                        WRITE(IW,950) XPT,YPT,ZPT,I
            IF(RR.LT.TM3) GO TO 210
               ELPOTN = ELPOTN + ZAN(I)/RR
  210    CONTINUE
C
         ELPOTT = ELPOTE + ELPOTN
         X(IPOT+IPT-1) = ELPOTT
         QVPOT(IPT) = ELPOTT
C         WRITE(95,960) IPT,XPT,YPT,ZPT,ELPOTE,ELPOTN,ELPOTT
  240 CONTINUE
C
      CALL RETFM(NEED)
      RETURN
C
  950 FORMAT(/1H ,'*** WARNING - ELECTROSTATIC POTENTIAL AT ',
     *      3F10.5,'. CONTRIBUTION FROM NUCLEUS ',I3,' IGNORED',/)
C
      END
C*MODULE PRPEL   *DECK COSPTS
      SUBROUTINE COSPTS(IPOINT,X,Y,Z,ISTAT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (NPPA=1082, MXPTPT=100)
C
      COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
      COMMON /COSMO3/ COSZAN(NPPA),CORZAN(3,NPPA)
C
C
C     ----- RETURN THE COORDINATES OF THE NEXT POINT
C           WHOSE PROPERTIES ARE TO BE COMPUTED -----
C
C     PRPLOC - SOLVENT ACCESSIBLE SURFACE
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(IPOINT.GT.NPOINT) ISTAT=1
         IF(IPOINT.GT.NPOINT) RETURN
C         X = COSURF(1,IPOINT)
C         Y = COSURF(2,IPOINT)
C         Z = COSURF(3,IPOINT)
C         X = C(1,IPOINT+NAT)*TOBOHR
C         Y = C(2,IPOINT+NAT)*TOBOHR
C         Z = C(3,IPOINT+NAT)*TOBOHR
         X = CORZAN(1,IPOINT)
         Y = CORZAN(2,IPOINT)
         Z = CORZAN(3,IPOINT)
         RETURN
C
      END
C*MODULE PRPEL   *DECK DENMUL
      SUBROUTINE DENMUL
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PAPER,GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 DUMMY
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION ELDEN(1), IPAIRS(2,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 /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 /MASSES/ ZMASS(MXATM)
      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 /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/8HELEDNC  /
      DATA ELDEN/8HELDEN   /
      DATA SIX/6.0D+00/, CMASS/8HCMASS   /
C
C     CALCULATE ELECTRON DENSITY FOR MOLECULAR ORBITAL OR ENTIRE
C     WAVE FUNCTION AT GRID POINTS THEN COLLECT THE CHARGES
C     ON ATOMIC CENTERS AND BOND MIDPOINTS AS MULTIPOLAR
C     EXPANSION POINTS IN EFP FORMAT.
C
      NATBAK = NAT
      XBAK   = C(1,1)
      YBAK   = C(2,1)
      ZBAK   = C(3,1)
      XPBAK  = XP
      YPBAK  = YP
      ZPBAK  = ZP
      ZANBAK = ZAN(1)
      ZMASSBK= ZMASS(1)
      ANAMBK = ANAM(1)
C
      PAPER = IEDOUT.GE.0
      IF(EXETYP.EQ.DEBUG)  PAPER=.TRUE.
      IF(EXETYP.EQ.DBUGME) PAPER=.TRUE.
      PAPER = PAPER .AND. MASWRK
C
      IF(PAPER) WRITE(IW,920) 'DENSITY BASED MULTIPOLAR EXPANSION'
C
C     -- SET POINTERS FOR PARTITIONING OF MEMORY
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
      NVAL = 1
      MXGRID = 100000 * NAT
C
C     -- ALLOCATE MEMORY
C
      CALL VALFM(LOADFM)
      LCRD = LOADFM + 1
      LCHG = LCRD   + MXGRID*3
      LMON = LCHG   + MXGRID
      LDIP = LMON   + 3*NAT
      LQUA = LDIP   + 3*3*NAT
      LOCT = LQUA   + 6*3*NAT
      IELD = LOCT   + 10*3*NAT
      IEDW = IELD   + L2
C
      IDENSA = IEDW   + 784
      IDENSB = IDENSA + L2
      LAST1  = IDENSB + L2 - 1
C
      IMOS   = IEDW   + 784
      INDEX  = IMOS   + L3
      LAST2  = INDEX  + L1 - 1
C
      LAST = MAX(LAST1,LAST2)
      NEED = LAST - LCRD + 1
      CALL GETFM(NEED)
      IF (EXETYP.EQ.CHECK) THEN
         CALL RETFM(NEED)
         RETURN
      END IF
C
C     -- OBTAIN TOTAL DENSITY MATRIX, IF NEEDED
C
      IF(MODENS.EQ.0) CALL DENDD1(X(IDENSA),X(IDENSB),L2)
C
C     -- DEFINE 3-D BOX FOR THE MOLECULE
C
      DO NCRD = 1, 3
        TMAX = C(NCRD,1)
        TMIN = C(NCRD,1)
        DO IAT = 2, NAT
          IF(C(NCRD,IAT).GT.TMAX) TMAX = C(NCRD,IAT)
          IF(C(NCRD,IAT).LT.TMIN) TMIN = C(NCRD,IAT)
        ENDDO
        IF(NCRD.EQ.1) XMAX = INT(TMAX + SIX)
        IF(NCRD.EQ.1) XMIN = INT(TMIN - SIX)
        IF(NCRD.EQ.2) YMAX = INT(TMAX + SIX)
        IF(NCRD.EQ.2) YMIN = INT(TMIN - SIX)
        IF(NCRD.EQ.3) ZMAX = INT(TMAX + SIX)
        IF(NCRD.EQ.3) ZMIN = INT(TMIN - SIX)
      ENDDO
      NGRIDX = INT(XMAX - XMIN)
      NGRIDY = INT(YMAX - YMIN)
      NGRIDZ = INT(ZMAX - ZMIN)
C
C     -- LOOP OVER GRID POINTS TO BE EVALUATED
C
      TOTAL = 0.0D+00
      IPOINT = 0
C
      DO 300 IX = 1, NGRIDX
      DO 300 IY = 1, NGRIDY
      DO 300 IZ = 1, NGRIDZ
      IPOINT = IPOINT + 1
      XP = XMIN - 0.50D+00 + IX
      YP = YMIN - 0.50D+00 + IY
      ZP = ZMIN - 0.50D+00 + IZ
C
C     -- CALCULATE DELTA FUNCTION INTEGRALS
C
      CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
C     -- CALCULATE ELECTRON DENSITY FOR ORBITAL OR WAVEFUNCTION
C
      IF(MODENS.EQ.0) THEN
         EDENS = TRACEP(X(IDENSA),X(IELD),L1)
      ELSE
         EDENS = 0.0D+00
         CALL CALCMO(MODENS,X(IELD),X(IMOS),X(INDEX),L1,EDENS)
      END IF
C
C     -- USE J-SUB-GRID
C
      IF (EDENS.GT.2.0D-03) THEN
         IPOINT = IPOINT - 1
         DO 310 JX = 0, 1
         DO 310 JY = 0, 1
         DO 310 JZ = 0, 1
         IPOINT = IPOINT + 1
         XP = XP - 0.25D+00 + 0.50D+00*JX
         YP = YP - 0.25D+00 + 0.50D+00*JY
         ZP = ZP - 0.25D+00 + 0.50D+00*JZ
C
         CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
         IF(MODENS.EQ.0) THEN
            EDENS = TRACEP(X(IDENSA),X(IELD),L1)
         ELSE
            EDENS = 0.0D+00
            CALL CALCMO(MODENS,X(IELD),
     *                  X(IMOS),X(INDEX),L1,EDENS)
         END IF
C
C        -- USE K-SUB-GRID
C
         IF (EDENS.GT.2.0D-02) THEN
            IPOINT = IPOINT - 1
            DO 320 KX = 0, 1
            DO 320 KY = 0, 1
            DO 320 KZ = 0, 1
            IPOINT = IPOINT + 1
            XP = XP - 0.125D+00 + 0.25D+00*KX
            YP = YP - 0.125D+00 + 0.25D+00*KY
            ZP = ZP - 0.125D+00 + 0.25D+00*KZ
C
            CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
            IF(MODENS.EQ.0) THEN
               EDENS = TRACEP(X(IDENSA),X(IELD),L1)
            ELSE
               EDENS = 0.0D+00
               CALL CALCMO(MODENS,X(IELD),
     *                     X(IMOS),X(INDEX),L1,EDENS)
            END IF
C
C           -- USE L-SUB-GRID
C
            IF (EDENS.GT.2.0D-01) THEN
               IPOINT = IPOINT - 1
               DO 330 LX = 0, 1
               DO 330 LY = 0, 1
               DO 330 LZ = 0, 1
               IPOINT = IPOINT + 1
               XP = XP - 0.0625D+00 + 0.125D+00*LX
               YP = YP - 0.0625D+00 + 0.125D+00*LY
               ZP = ZP - 0.0625D+00 + 0.125D+00*LZ
C
               CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
               IF(MODENS.EQ.0) THEN
                  EDENS = TRACEP(X(IDENSA),X(IELD),L1)
               ELSE
                  EDENS = 0.0D+00
                  CALL CALCMO(MODENS,X(IELD),
     *                        X(IMOS),X(INDEX),L1,EDENS)
               END IF
C
C              -- USE M-SUB-GRID
C
               IF (EDENS.GT.2.0D+00) THEN
                  IPOINT = IPOINT - 1
                  DO 340 MX = 0, 1
                  DO 340 MY = 0, 1
                  DO 340 MZ = 0, 1
                  IPOINT = IPOINT + 1
                  XP =XP - 0.03125D+00 + 0.0625D+00*MX
                  YP =YP - 0.03125D+00 + 0.0625D+00*MY
                  ZP =ZP - 0.03125D+00 + 0.0625D+00*MZ
C
                  CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
                  IF(MODENS.EQ.0) THEN
                     EDENS = TRACEP(X(IDENSA),X(IELD),L1)
                  ELSE
                     EDENS = 0.0D+00
                     CALL CALCMO(MODENS,X(IELD),
     *                           X(IMOS),X(INDEX),L1,EDENS)
                  END IF
C
C                 -- USE N-SUB-GRID
C
                  IF (EDENS.GT.2.0D+01) THEN
                     IPOINT = IPOINT - 1
                     DO 350 NX = 0, 1
                     DO 350 NY = 0, 1
                     DO 350 NZ = 0, 1
                     IPOINT = IPOINT + 1
                     XP =XP - 0.015625D+00 + 0.03125D+00*NX
                     YP =YP - 0.015625D+00 + 0.03125D+00*NY
                     ZP =ZP - 0.015625D+00 + 0.03125D+00*NZ
C
                     CALL PRCALC(ELDEN,X(IELD),X(IEDW),NVAL,L2,.FALSE.)
C
                     IF(MODENS.EQ.0) THEN
                        EDENS = TRACEP(X(IDENSA),X(IELD),L1)
                     ELSE
                        EDENS = 0.0D+00
                        CALL CALCMO(MODENS,X(IELD),
     *                              X(IMOS),X(INDEX),L1,EDENS)
                     END IF
C
C                    -- USE N2-SUB-GRID
C
                     IF (EDENS.GT.2.0D+02) THEN
                        IPOINT = IPOINT - 1
                        DO 360 N2X = 0, 1
                        DO 360 N2Y = 0, 1
                        DO 360 N2Z = 0, 1
                        IPOINT = IPOINT + 1
                        XP=XP- 0.0078125D+00 +0.015625D+00*N2X
                        YP=YP- 0.0078125D+00 +0.015625D+00*N2Y
                        ZP=ZP- 0.0078125D+00 +0.015625D+00*N2Z
C
                        CALL PRCALC(ELDEN,X(IELD),X(IEDW),
     *                              NVAL,L2,.FALSE.)
C
                        IF(MODENS.EQ.0) THEN
                           EDENS = TRACEP(X(IDENSA),X(IELD),L1)
                        ELSE
                           EDENS = 0.0D+00
                           CALL CALCMO(MODENS,X(IELD),
     *                                 X(IMOS),X(INDEX),L1,EDENS)
                        END IF
C
C                       -- USE N3-SUB-GRID
C
                        IF (EDENS.GT.2.0D+03) THEN
                           IPOINT = IPOINT - 1
                           DO 370 N3X = 0, 1
                           DO 370 N3Y = 0, 1
                           DO 370 N3Z = 0, 1
                           IPOINT = IPOINT + 1
                           XP=XP-3.90625D-03 +7.8125D-03*N3X
                           YP=YP-3.90625D-03 +7.8125D-03*N3Y
                           ZP=ZP-3.90625D-03 +7.8125D-03*N3Z
C
                           CALL PRCALC(ELDEN,X(IELD),X(IEDW),
     *                                 NVAL,L2,.FALSE.)
C
                           IF(MODENS.EQ.0) THEN
                              EDENS = TRACEP(X(IDENSA),X(IELD),L1)
                           ELSE
                              EDENS = 0.0D+00
                              CALL CALCMO(MODENS,X(IELD),
     *                                    X(IMOS),X(INDEX),L1,EDENS)
                           END IF
C
C                          -- USE N4-SUB-GRID
C
                           IF (EDENS.GT.2.0D+04) THEN
                              IPOINT = IPOINT - 1
                              DO 380 N4X = 0, 1
                              DO 380 N4Y = 0, 1
                              DO 380 N4Z = 0, 1
                              IPOINT = IPOINT + 1
                              XP=XP-1.953125D-03 +3.90625D-03*N4X
                              YP=YP-1.953125D-03 +3.90625D-03*N4Y
                              ZP=ZP-1.953125D-03 +3.90625D-03*N4Z
C
                              CALL PRCALC(ELDEN,X(IELD),X(IEDW),
     *                                    NVAL,L2,.FALSE.)
C
                              IF(MODENS.EQ.0) THEN
                                 EDENS = TRACEP(X(IDENSA),X(IELD),L1)
                              ELSE
                                 EDENS = 0.0D+00
                                 CALL CALCMO(MODENS,X(IELD),
     *                                       X(IMOS),X(INDEX),L1,EDENS)
                              END IF
                              X(LCRD+(IPOINT-1)*3  ) = XP
                              X(LCRD+(IPOINT-1)*3+1) = YP
                              X(LCRD+(IPOINT-1)*3+2) = ZP
                              DUM = EDENS*0.125D+00**8
                              X(LCHG-1+IPOINT) = DUM
                              TOTAL = TOTAL + DUM
                              XP=XP+1.953125D-03 -3.90625D-03*N4X
                              YP=YP+1.953125D-03 -3.90625D-03*N4Y
                              ZP=ZP+1.953125D-03 -3.90625D-03*N4Z
 380                          CONTINUE
                           ELSE
                              X(LCRD+(IPOINT-1)*3  ) = XP
                              X(LCRD+(IPOINT-1)*3+1) = YP
                              X(LCRD+(IPOINT-1)*3+2) = ZP
                              DUM = EDENS*0.125D+00**7
                              X(LCHG-1+IPOINT) = DUM
                              TOTAL = TOTAL + DUM
                           END IF
                           XP=XP+3.90625D-03 -7.8125D-03*N3X
                           YP=YP+3.90625D-03 -7.8125D-03*N3Y
                           ZP=ZP+3.90625D-03 -7.8125D-03*N3Z
 370                       CONTINUE
                        ELSE
                           X(LCRD+(IPOINT-1)*3  ) = XP
                           X(LCRD+(IPOINT-1)*3+1) = YP
                           X(LCRD+(IPOINT-1)*3+2) = ZP
                           DUM = EDENS*0.125D+00**6
                           X(LCHG-1+IPOINT) = DUM
                           TOTAL = TOTAL + DUM
                        END IF
                        XP=XP+ 0.0078125D+00 -0.015625D+00*N2X
                        YP=YP+ 0.0078125D+00 -0.015625D+00*N2Y
                        ZP=ZP+ 0.0078125D+00 -0.015625D+00*N2Z
 360                    CONTINUE
                     ELSE
                        X(LCRD+(IPOINT-1)*3  ) = XP
                        X(LCRD+(IPOINT-1)*3+1) = YP
                        X(LCRD+(IPOINT-1)*3+2) = ZP
                        DUM = EDENS*0.125D+00**5
                        X(LCHG-1+IPOINT) = DUM
                        TOTAL = TOTAL + DUM
                     END IF
                     XP =XP + 0.015625D+00 - 0.03125D+00*NX
                     YP =YP + 0.015625D+00 - 0.03125D+00*NY
                     ZP =ZP + 0.015625D+00 - 0.03125D+00*NZ
 350                 CONTINUE
                  ELSE
                     X(LCRD+(IPOINT-1)*3  ) = XP
                     X(LCRD+(IPOINT-1)*3+1) = YP
                     X(LCRD+(IPOINT-1)*3+2) = ZP
                     DUM = EDENS*0.125D+00**4
                     X(LCHG-1+IPOINT) = DUM
                     TOTAL = TOTAL + DUM
                  END IF
                  XP =XP + 0.03125D+00 - 0.0625D+00*MX
                  YP =YP + 0.03125D+00 - 0.0625D+00*MY
                  ZP =ZP + 0.03125D+00 - 0.0625D+00*MZ
 340              CONTINUE
               ELSE
                  X(LCRD+(IPOINT-1)*3  ) = XP
                  X(LCRD+(IPOINT-1)*3+1) = YP
                  X(LCRD+(IPOINT-1)*3+2) = ZP
                  DUM = EDENS*0.125D+00**3
                  X(LCHG-1+IPOINT) = DUM
                  TOTAL = TOTAL + DUM
               END IF
               XP = XP + 0.0625D+00 - 0.125D+00*LX
               YP = YP + 0.0625D+00 - 0.125D+00*LY
               ZP = ZP + 0.0625D+00 - 0.125D+00*LZ
 330           CONTINUE
            ELSE
               X(LCRD+(IPOINT-1)*3  ) = XP
               X(LCRD+(IPOINT-1)*3+1) = YP
               X(LCRD+(IPOINT-1)*3+2) = ZP
               DUM = EDENS*0.125D+00**2
               X(LCHG-1+IPOINT) = DUM
               TOTAL = TOTAL + DUM
            END IF
            XP = XP + 0.125D+00 - 0.25D+00*KX
            YP = YP + 0.125D+00 - 0.25D+00*KY
            ZP = ZP + 0.125D+00 - 0.25D+00*KZ
 320        CONTINUE
         ELSE
            X(LCRD+(IPOINT-1)*3  ) = XP
            X(LCRD+(IPOINT-1)*3+1) = YP
            X(LCRD+(IPOINT-1)*3+2) = ZP
            DUM = EDENS*0.125D+00
            X(LCHG-1+IPOINT) = DUM
            TOTAL = TOTAL + DUM
         END IF
         XP = XP + 0.25D+00 - 0.50D+00*JX
         YP = YP + 0.25D+00 - 0.50D+00*JY
         ZP = ZP + 0.25D+00 - 0.50D+00*JZ
 310     CONTINUE
      ELSE
         X(LCRD+(IPOINT-1)*3  ) = XP
         X(LCRD+(IPOINT-1)*3+1) = YP
         X(LCRD+(IPOINT-1)*3+2) = ZP
         DUM = EDENS
         X(LCHG-1+IPOINT) = DUM
         TOTAL = TOTAL + DUM
      END IF
 300  CONTINUE
C
      SCALE = NE/TOTAL
      DO JPT = 1, IPOINT
         X(LCHG-1+JPT) = X(LCHG-1+JPT) * SCALE
      END DO
      IF(PAPER) THEN
         WRITE(IW,*) IPOINT,'GRID POINTS ARE USED'
         WRITE(IW,*) 'TOTAL NUMBER OF ELECTRONS =', NE
         WRITE(IW,'(A,F9.5)')
     *     ' TOTAL CHARGE OBTAINED FROM GRID INTEGRAL =', TOTAL
         WRITE(IW,'(A,F7.5,A)')
     *     ' A SCALING FACTOR OF ', SCALE, ' IS APPLIED'
         WRITE(IW,*)
     *     'THE EXPANSION POINTS ARE PUNCHED OUT IN THE DAT FILE'
         WRITE(IW,*) ' '
      END IF
C
C     -- MULTIPOLAR EXPANSION
C        ATOMS, ATOM+MIDBOND OR MASS CENTER
C     METHOD = 1   :  ATOMS
C              2   :  ATOMS + BOND MIDPOINTS
C              3   :  MASS CENTER
      METHOD = 2
C
      CALL CALCOM(CMX,CMY,CMZ)
      NBONDS = 0
      IF(METHOD.EQ.2) THEN
         CALL GTBOND(IPAIRS,MXATM,NBONDS,0)
C        - NOTE: WE STORE THE BOND MIDPOINTS AS ATOMS
         DO I = 1,NBONDS
            NAT = NAT + 1
            LATM1 = IPAIRS(1,I)
            LATM2 = IPAIRS(2,I)
            DO LC=1,3
               C(LC,NAT)=0.50D+00*(C(LC,LATM1)+C(LC,LATM2))
            END DO
            IF (LATM1.LT.10.AND.LATM2.LT.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I1,I1)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.GE.10.AND.LATM2.LT.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I2,I1)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.LT.10.AND.LATM2.GE.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I1,I2)') 'BO',LATM1,LATM2
            ELSE IF (LATM1.GE.10.AND.LATM2.GE.10) THEN
               WRITE(UNIT=DUMMY,FMT='(A2,I2,I2)') 'BO',LATM1,LATM2
            END IF
            READ(UNIT=DUMMY, FMT='(A8)') ANAM(NAT)
         END DO
      END IF
C
      IF(METHOD.EQ.3) THEN
         JAT = 1
         NAT = 1
         C(1,1) = CMX
         C(2,1) = CMY
         C(3,1) = CMZ
         ZAN(1) = 0.0D+00
         ZMASS(1)= 0.0D+00
         ANAM(1)= CMASS
      END IF
C
      DO 400 JPT = 1, IPOINT
         XPT = X(LCRD+(JPT-1)*3  )
         YPT = X(LCRD+(JPT-1)*3+1)
         ZPT = X(LCRD+(JPT-1)*3+2)
         QPT = X(LCHG+JPT-1)
C
C        - FIND THE CLOSEST EXPANSION POINT
         IF(METHOD.EQ.1 .OR. METHOD.EQ.2) THEN
            TEMPR2 = 1.0D+10
            JAT = 1
            DO 410 IAT = 1, NAT
               R2 = (XPT-C(1,IAT))*(XPT-C(1,IAT))
     *             +(YPT-C(2,IAT))*(YPT-C(2,IAT))
     *             +(ZPT-C(3,IAT))*(ZPT-C(3,IAT))
               IF(R2.LT.TEMPR2) THEN
                  TEMPR2 = R2
                  JAT    = IAT
               END IF
 410        CONTINUE
         END IF
C
C        - DISTANCES
         DX = XPT - C(1,JAT)
         DY = YPT - C(2,JAT)
         DZ = ZPT - C(3,JAT)
C        - MONOPOLES
         X(LMON+JAT-1) = X(LMON+JAT-1) - QPT
C        - DIPOLES
         X(LDIP+(JAT-1)*3  ) = X(LDIP+(JAT-1)*3  ) - QPT*DX
         X(LDIP+(JAT-1)*3+1) = X(LDIP+(JAT-1)*3+1) - QPT*DY
         X(LDIP+(JAT-1)*3+2) = X(LDIP+(JAT-1)*3+2) - QPT*DZ
C        - QUADRUPOLES
         X(LQUA+(JAT-1)*6  ) = X(LQUA+(JAT-1)*6  ) - QPT*DX*DX
         X(LQUA+(JAT-1)*6+1) = X(LQUA+(JAT-1)*6+1) - QPT*DY*DY
         X(LQUA+(JAT-1)*6+2) = X(LQUA+(JAT-1)*6+2) - QPT*DZ*DZ
         X(LQUA+(JAT-1)*6+3) = X(LQUA+(JAT-1)*6+3) - QPT*DX*DY
         X(LQUA+(JAT-1)*6+4) = X(LQUA+(JAT-1)*6+4) - QPT*DX*DZ
         X(LQUA+(JAT-1)*6+5) = X(LQUA+(JAT-1)*6+5) - QPT*DY*DZ
C        - OCTUPOLES
         X(LOCT+(JAT-1)*10  ) = X(LOCT+(JAT-1)*10  ) - QPT*DX*DX*DX
         X(LOCT+(JAT-1)*10+1) = X(LOCT+(JAT-1)*10+1) - QPT*DY*DY*DY
         X(LOCT+(JAT-1)*10+2) = X(LOCT+(JAT-1)*10+2) - QPT*DZ*DZ*DZ
         X(LOCT+(JAT-1)*10+3) = X(LOCT+(JAT-1)*10+3) - QPT*DX*DX*DY
         X(LOCT+(JAT-1)*10+4) = X(LOCT+(JAT-1)*10+4) - QPT*DX*DX*DZ
         X(LOCT+(JAT-1)*10+5) = X(LOCT+(JAT-1)*10+5) - QPT*DX*DY*DY
         X(LOCT+(JAT-1)*10+6) = X(LOCT+(JAT-1)*10+6) - QPT*DY*DY*DZ
         X(LOCT+(JAT-1)*10+7) = X(LOCT+(JAT-1)*10+7) - QPT*DX*DZ*DZ
         X(LOCT+(JAT-1)*10+8) = X(LOCT+(JAT-1)*10+8) - QPT*DY*DZ*DZ
         X(LOCT+(JAT-1)*10+9) = X(LOCT+(JAT-1)*10+9) - QPT*DX*DY*DZ
 400  CONTINUE
C
      WRITE(IP,*) ' '
      WRITE(IP,'(A)') 'DENSITY BASED MULTIPOLAR EXPANSION'
      WRITE(IP,*) 'COORDINATES'
      DO I=1,NAT
         WRITE(IP  ,'(A8,3F15.10,F12.7,F5.1)')
     *     ANAM(I),(C(J,I),J=1,3),ZMASS(I),ZAN(I)
      ENDDO
      WRITE(IP,*) 'STOP'
      WRITE(IP,*) 'MONOPOLES'
      DO I=1,NAT
         WRITE(IP,'(A8,F15.10,F10.5)')
     *     ANAM(I),X(LMON+I-1),ZAN(I)
      ENDDO
      WRITE(IP,*) 'STOP'
      WRITE(IP,*) 'DIPOLES'
      DO I=1,NAT
         WRITE(IP,'(A8,3(1X,F15.10))')
     *     ANAM(I),X(LDIP+(I-1)*3),X(LDIP+(I-1)*3+1),
     *     X(LDIP+(I-1)*3+2)
      ENDDO
      WRITE(IP,*) 'STOP'
      WRITE(IP,*) 'QUADRUPOLES'
      DO I=1,NAT
         WRITE(IP,'(A8,4(1X,F15.10),A2,/2(1X,F15.10))')
     *     ANAM(I),X(LQUA+(I-1)*6),X(LQUA+(I-1)*6+1),
     *     X(LQUA+(I-1)*6+2),
     *     X(LQUA+(I-1)*6+3),' >',X(LQUA+(I-1)*6+4),
     *     X(LQUA+(I-1)*6+5)
      ENDDO
      WRITE(IP,*) 'STOP'
      WRITE(IP,*) 'OCTUPOLES'
      DO I=1,NAT
         WRITE(IP,'(A8,4(1X,F15.10),A2,/4(1X,F15.10),A2)')
     *     ANAM(I),X(LOCT+(I-1)*10),X(LOCT+(I-1)*10+1),
     *     X(LOCT+(I-1)*10+2),
     *     X(LOCT+(I-1)*10+3),' >',X(LOCT+(I-1)*10+4),
     *     X(LOCT+(I-1)*10+5),
     *     X(LOCT+(I-1)*10+6),X(LOCT+(I-1)*10+7),
     *     ' >',X(LOCT+(I-1)*10+8),
     *     X(LOCT+(I-1)*10+9)
      ENDDO
      WRITE(IP,*) 'STOP'
      WRITE(IP,*) '$END'
C
C
C     A SECOND COPY OF THE "EFFECTIVE FRAGMENT INPUT FILE"
C     CAN BE WRITTEN TO XXX.IRC TO
C     BE READ BY WALT STEVENS' CODE.
C
C      IF(RUNTYP.EQ.MAKEFP) THEN
C         IDMP =4
C         CALL SEQOPN(IDMP,'IRCDATA','NEW',.FALSE.,'FORMATTED')
C      ELSE
         IDMP=-1
C      END IF
C
      IF (IDMP.GT.0) THEN
        WRITE(IDMP,*) ' '
        WRITE(IDMP,'(A)') 'DENSITY BASED MULTIPOLAR EXPANSION'
        WRITE(IDMP,*) 'COORDINATES'
        DO I=1,NAT
         WRITE(IDMP  ,'(A8,3F15.10,F12.7,F5.1)')
     *     ANAM(I),(C(J,I),J=1,3),ZMASS(I),ZAN(I)
        ENDDO
        WRITE(IDMP,*) 'STOP'
        WRITE(IDMP,*) 'MONOPOLES'
        DO I=1,NAT
           WRITE(IDMP,'(A8,F15.10,F10.5)')
     *          ANAM(I),X(LMON+I-1),ZAN(I)
        ENDDO
        WRITE(IDMP,*) 'STOP'
        WRITE(IDMP,*) 'DIPOLES'
        DO I=1,NAT
           WRITE(IDMP,'(A8,3(1X,F15.10))')
     *          ANAM(I),X(LDIP+(I-1)*3),X(LDIP+(I-1)*3+1),
     *          X(LDIP+(I-1)*3+2)
        ENDDO
        WRITE(IDMP,*) 'STOP'
        WRITE(IDMP,*) 'QUADRUPOLES'
        DO I=1,NAT
           WRITE(IDMP,'(A8,4(1X,F15.10),A2,/2(1X,F15.10))')
     *          ANAM(I),X(LQUA+(I-1)*6),X(LQUA+(I-1)*6+1),
     *          X(LQUA+(I-1)*6+2),
     *          X(LQUA+(I-1)*6+3),' >',X(LQUA+(I-1)*6+4),
     *          X(LQUA+(I-1)*6+5)
        ENDDO
        WRITE(IDMP,*) 'STOP'
        WRITE(IDMP,*) 'OCTUPOLES'
        DO I=1,NAT
           WRITE(IDMP,'(A8,4(1X,F15.10),A2,/4(1X,F15.10),A2)')
     *          ANAM(I),X(LOCT+(I-1)*10),X(LOCT+(I-1)*10+1),
     *          X(LOCT+(I-1)*10+2),
     *          X(LOCT+(I-1)*10+3),' >',X(LOCT+(I-1)*10+4),
     *          X(LOCT+(I-1)*10+5),
     *          X(LOCT+(I-1)*10+6),X(LOCT+(I-1)*10+7),
     *          ' >',X(LOCT+(I-1)*10+8),
     *          X(LOCT+(I-1)*10+9)
        ENDDO
        WRITE(IDMP,*) 'STOP'
        WRITE(IDMP,*) '$END'
      END IF
C
      NAT    = NATBAK
      C(1,1) = XBAK
      C(2,1) = YBAK
      C(3,1) = ZBAK
      XP     = XPBAK
      YP     = YPBAK
      ZP     = ZPBAK
      ZAN(1) = ZANBAK
      ZMASS(1)= ZMASSBK
      ANAM(1)= ANAMBK
C
      CALL RETFM(NEED)
      RETURN
C
  920 FORMAT(/10X,34('-')/10X,A34,/10X,34('-'))
      END
C*MODULE PRPEL   *DECK EPOTEN
      SUBROUTINE EPOTEN(ACONST,XCOORD,YCOORD,ZCOORD,VALUE,DA,NP,L2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,NORM,DOUBLE,GOPARR,DSKWRK,MASWRK,NXT
C
      DIMENSION XCOORD(NP),YCOORD(NP),ZCOORD(NP),VALUE(NP),DA(L2)
      DIMENSION IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35),
     *          XIN(125),YIN(125),ZIN(125),
     *          DIJ(225),IJX(225),IJY(225),IJZ(225),
     *          HP(28),WP(28),MINP(7),MAXP(7)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /HERMIT/ H1,H2(2),H3(3),H4(4),H5(5),H6(6),H7(7),
     *                H8(8),H9(9),H10(10)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEFPCM,IP_F,NFMOPCM
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /WERMIT/ W1,W2(2),W3(3),W4(4),W5(5),W6(6),W7(7),
     *                W8(8),W9(9),W10(10)
C
      EQUIVALENCE (HP(1),H1),(WP(1),W1)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (PI212=1.1283791670955D+00)
      PARAMETER (SQRT3=1.73205080756888D+00)
      PARAMETER (SQRT5=2.23606797749979D+00)
      PARAMETER (SQRT7=2.64575131106459D+00)
      PARAMETER (RLN10=2.30258D+00)
      DATA MINP /1,2,4,7,11,16,22/
      DATA MAXP /1,3,6,10,15,21,28/
      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/
      DATA IX / 1, 6, 1, 1,11, 1, 1, 6, 6, 1,
     *         16, 1, 1,11,11, 6, 1, 6, 1, 6,
     *         21, 1, 1,16,16, 6, 1, 6, 1,11,
     *         11, 1,11, 6, 6/
      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/
      DATA IY / 1, 1, 6, 1, 1,11, 1, 6, 1, 6,
     *          1,16, 1, 6, 1,11,11, 1, 6, 6,
     *          1,21, 1, 6, 1,16,16, 1, 6,11,
     *          1,11, 6,11, 6/
      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/
      DATA IZ / 1, 1, 1, 6, 1, 1,11, 1, 6, 6,
     *          1, 1,16, 1, 6, 1, 6,11,11, 6,
     *          1, 1,21, 1, 6, 1, 6,16,16, 1,
     *         11,11, 6, 6,11/
C
C     EVALUATE THE VALUES OF THE ELECTROSTATIC POTENTIAL FOR ALL
C     POINTS GIVEN IN X,Y,ZCOORD.
C     AT PRESENT THERE ARE NO APPROXIMATIONS AND NO SCREENING OTHER
C     THAN THE USUAL INTEGRAL SCREENING WITH ITOL. THE INTEGRALS
C     ARE NEITHER ZEROED OUT AT THE BEGINNING NOR PARALLEL SUMMED
C     AT THE END.
C     THIS IS A CLONE OF EPCALC.
C
C     REDUCE THE ACCURACY 50 TIMES (E**4) FOR ITERATIVE SOLVER IN PCM
C     (AT PRESENT THIS ROUTINE IS ONLY CALLED IN CASE OF ITERATIVE PCM
C     SOLVER AND THE REASON FOR ACCURACY REDUCTION IS ITS LIMITED
C     NUMERIC ACCURACY).
C
C     THIS IS USED BY PCM AND FMO, AND IS LIMITED TO SPDFG ONLY.
C
      CALL BASCHK(LMAX)
      IF(LMAX.GT.4) THEN
         IF(MASWRK) WRITE(IW,*) 'EPOTEN: CODE IS SPDFG LIMITED'
         CALL ABRT
      END IF
C
      ITOLPCM=ITOL
      IF(IPCM.NE.0.AND.IPCMIT.NE.0) ITOLPCM=ITOLPCM-4
      TOL = RLN10*ITOLPCM
C     ITOLPCM=20
C     TOL  = RLN10*ITOLPCM
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NXT = IBTYP.EQ.1
      NEXT  = -1
      KOUNT = -1
C
C     PREPARE DENSITY FOR EASY TRACING: DOUBLE OFF-DIAGONAL ELEMENTS.
C
      L1=NUM
      LOOP=0
      DO I=1,L1
        DO J=1,I-1
          LOOP=LOOP+1
          DA(LOOP)=DA(LOOP)*TWO
        ENDDO
C       SKIP THE DIAGONAL
        LOOP=LOOP+1
      ENDDO
C
C                    LOOP OVER SHELLS II
C
      DO 510 II=1,NSHELL
        I    = KATOM(II)
        XI   = C(1,I)
        YI   = C(2,I)
        ZI   = C(3,I)
        I1   = KSTART(II)
        I2   = I1 + KNG(II) - 1
        LIT  = KTYPE(II)
        MINI = KMIN(II)
        MAXI = KMAX(II)
        LOCI = KLOC(II) - MINI
C
C                    LOOP OVER SHELLS JJ
C
        DO 500  JJ=1,II
C
          IF(GOPARR) THEN
            KOUNT=KOUNT+1
            IF(NXT) THEN
              IF(KOUNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
              IF(NEXT.NE.KOUNT) GOTO 500
            ELSE
              IF(MOD(KOUNT,NPROC).NE.ME) GOTO 500
            ENDIF
          ENDIF
          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
          NROOTS = (LIT + LJT - 2)/2 + 1
          IANDJ  = II.EQ.JJ
C
C             PREPARE INDICES FOR PAIRS OF (I,J) ORBITALS
C
          IJ = 0
          MAX = MAXJ
          DO 100  I=MINI,MAXI
            NX = IX(I)
            NY = IY(I)
            NZ = IZ(I)
            IF (IANDJ) MAX = I
            DO 100  J=MINJ,MAX
              IJ = IJ+1
              IJX(IJ) = NX+JX(J)
              IJY(IJ) = NY+JY(J)
              IJZ(IJ) = NZ+JZ(J)
  100       CONTINUE
C
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)
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)
C
              AA  = AI + AJ
              AA1 = ONE/AA
              FI  = PI212*AA1
C
              AAX = (AI*XI + AJ*XJ)
              AAY = (AI*YI + AJ*YJ)
              AAZ = (AI*ZI + AJ*ZJ)
C
              AX  = AAX*AA1
              AY  = AAY*AA1
              AZ  = AAZ*AA1
C
              DUM = AI*AJ*RR*AA1
              IF(DUM .GT. TOL) GO TO 400
              FAC = FI*EXP(-DUM)
C
C                       CALCULATE DENSITY FACTORS
C
              DOUBLE = IANDJ.AND.IG.NE.JG
              MAX = MAXJ
              NN  = 0
C
              DUM1 = ZERO
              DUM2 = ZERO
              DO 200 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
                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
                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
                IF(IANDJ) MAX = I
                DO 180 J = MINJ,MAX
                  NN = NN+1
                  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
C
                  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
C
                  ELSE IF ((J.EQ.14).AND.NORM) THEN
                    DUM2 = DUM2*SQRT5
C
                  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
C
                  ELSE IF ((J.EQ.24).AND.NORM) THEN
                    DUM2 = DUM2*SQRT7
C
                  ELSE IF ((J.EQ.30).AND.NORM) THEN
                    DUM2 = DUM2*SQRT5/SQRT3
C
                  ELSE IF ((J.EQ.33).AND.NORM) THEN
                    DUM2 = DUM2*SQRT3
C
                  END IF
C
                  DIJ(NN) = DUM2
  180           CONTINUE
  200         CONTINUE
C
C                    LOOP OVER PROPERTY COORDINATES
C
              DO 600 IPT=1,NP
                XPP = XCOORD(IPT)
                YPP = YCOORD(IPT)
                ZPP = ZCOORD(IPT)
C
C       CALCULATE POINTS AND WEIGHTS FOR RYS POLYNOMIAL
C
                XX  = AA * ((AX-XPP)**2 + (AY-YPP)**2 + (AZ-ZPP)**2)
                IF (NROOTS.LE.3) CALL RT123
                IF (NROOTS.EQ.4) CALL ROOT4
                IF (NROOTS.EQ.5) CALL ROOT5
                IF (NROOTS.GE.6) CALL ROOT6
C
C       LOOP OVER ROOTS OF RYS POLYNOMIAL TO CALCULATE INTEGRALS
C
                MM = 0
                DO 340  K=1,NROOTS
C
                  UU = AA*U(K)
                  WW = W(K)
                  TT = ONE/(AA+UU)
                  T  = SQRT(TT)
C
                  X0 = (AAX + UU*XPP)*TT
                  Y0 = (AAY + UU*YPP)*TT
                  Z0 = (AAZ + UU*ZPP)*TT
C
C      CALCULATE 1-DIMENSIONAL INTEGRALS OVER ALL ANGULAR MOMENTA
C
                  IN = -5+MM
                  DO 320  I=1,LIT
                    IN = IN+5
                    NI = I
C
                    DO 320  J=1,LJT
                      JN = IN+J
                      NJ = J
C
C       EVALUATE MOMENT INTEGRALS USING GAUSS-HERMITE QUADRATURE:
C
                      XINT0 = ZERO
                      YINT0 = ZERO
                      ZINT0 = ZERO
C
                      NPTS = (NI + NJ - 2)/2 + 1
                      IMIN = MINP(NPTS)
                      IMAX = MAXP(NPTS)
C
                      DO 310  IROOT=IMIN,IMAX
C
                        DUM = WP(IROOT)
                        PX = DUM
                        PY = DUM
                        PZ = DUM
C
                        DUM = HP(IROOT)*T
                        PTX = DUM + X0
                        PTY = DUM + Y0
                        PTZ = DUM + Z0
C
                        AXI = PTX - XI
                        AYI = PTY - YI
                        AZI = PTZ - ZI
C
                        BXI = PTX - XJ
                        BYI = PTY - YJ
                        BZI = PTZ - ZJ
C
                        GO TO (250,240,230,220,210),NI
C
  210                   PX = PX*AXI
                        PY = PY*AYI
                        PZ = PZ*AZI
C
  220                   PX = PX*AXI
                        PY = PY*AYI
                        PZ = PZ*AZI
C
  230                   PX = PX*AXI
                        PY = PY*AYI
                        PZ = PZ*AZI
C
  240                   PX = PX*AXI
                        PY = PY*AYI
                        PZ = PZ*AZI
C
  250                   CONTINUE
C
                        GO TO (300,290,280,270,260),NJ
C
  260                   PX = PX*BXI
                        PY = PY*BYI
                        PZ = PZ*BZI
C
  270                   PX = PX*BXI
                        PY = PY*BYI
                        PZ = PZ*BZI
C
  280                   PX = PX*BXI
                        PY = PY*BYI
                        PZ = PZ*BZI
C
  290                   PX = PX*BXI
                        PY = PY*BYI
                        PZ = PZ*BZI
C
  300                   CONTINUE
C
                        XINT0 = XINT0 + PX
                        YINT0 = YINT0 + PY
                        ZINT0 = ZINT0 + PZ
C
  310                 CONTINUE
C
                      XIN(JN) = XINT0
                      YIN(JN) = YINT0
                      ZIN(JN) = ZINT0*WW
C
  320               CONTINUE
C
                  MM = MM+25
  340           CONTINUE
C
                NN = 0
                MAX = MAXJ
                SUM = ZERO
C         COMPUTE V(I)=TR(D*V(I)), WHERE V(I) IS THE POTENTIAL
C         AT THE POINT IPT (==I) AND D IS THE ELECTRON DENSITY.
                DO 420  I=MINI,MAXI
                  LI = LOCI + I
                  IN = LI*(LI-1)/2
                  IF (IANDJ) MAX = I
                  DO 420  J=MINJ,MAX
                    LJ = LOCJ + J
                    JN = LJ + IN
                    NN = NN + 1
                    NX = IJX(NN)
                    NY = IJY(NN)
                    NZ = IJZ(NN)
                    MM = 0
                    SUM1 = ZERO
                    DO 350 K=1,NROOTS
                      SUM1=SUM1+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
                      MM = MM+25
  350               CONTINUE
                    SUM=SUM+DIJ(NN)*DA(JN)*SUM1
  420           CONTINUE
                VALUE(IPT) = VALUE(IPT) + ACONST*SUM
C
  600         CONTINUE
C
C        END OF LOOP OVER PROPERTY COORDINATES
C
  400       CONTINUE
  410    CONTINUE
C
C        END OF LOOPS OVER PRIMITIVES
C
C
  500   CONTINUE
  510 CONTINUE
C
C        END OF LOOPS OVER SHELLS
C
      IF(GOPARR.AND.NXT) CALL DDI_DLBRESET
C
C     RESTORE DENSITY IN COURTESY
C
      LOOP=0
      DO I=1,L1
        DO J=1,I-1
          LOOP=LOOP+1
          DA(LOOP)=DA(LOOP)/TWO
        ENDDO
C       SKIP THE DIAGONAL
        LOOP=LOOP+1
      ENDDO
C
      IF(NFMOPCM.NE.0) THEN
        IF(MASWRK) WRITE(IW,*) 'DONE PCM POT.'
        CALL TIMIT(1)
      ENDIF
C
      RETURN
      END
