C 11 APR 08 - MWS - IMABIX: SYNCHRONIZE CALL TO STONE
C 28 AUG 07 - HL  - IEFBEM: USE NEW CAV+DIS+REP GRADIENTS
C 24 MAR 07 - MWS - PAD THE FRGINF,EFPPAR COMMON BLOCKS
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C  7 APR 06 - MWS - REMOVE ICENT/NESFF/INA/INF TRACES
C  2 FEB 06 - MWS - UPDATE STORAGE FOR PRCALC CALL
C 19 NOV 05 - HL,DGF - IMPLEMENT PCM+EFP GRADIENTS, FMO+PCM CHANGES
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C 15 JUL 05 - MWS - FIX DYNAMIC STORAGE FOR DIPIT
C 13 JUL 05 - HL  - ADDPOT: PASS FORMAL DIMENSIONING ARGUMENT
C  5 JUL 05 - DGF - IXMATV: AVOID COPYING DATA BEFORE INITIALIZED
C  1 JUN 05 - MWS - ADDPOT: ADJUST TO HMN'S PARALLEL ELEC FIELD INTS
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, ROOT, FRGINF, EFPPAR
C  3 JUN 04 - HL  - ADDPOT: COMMENT OUT SOME UNWANTED CODE
C  7 APR 04 - HL  - IEF-PCM/EFP FOR ALL TYPES OF WAVEFUNCTION
C 12 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 16 JUN 03 - HL  - CHANGES FOR PCM GRADIENTS. ADD C-PCM
C  7 AUG 02 - HL,CP - ITERATIVE ISOTROPIC IEF-PCM/EFP
C                   IXMATV: INCLUDE PB_MUL,PB_POL,PC_MUL,PC_POL
C  8 OCT 01 - HL  - PARALLELIZE PCM AND NEW ROUTINE IEFCMV1
C 20 FEB 01 - BM  - IEFBEM,OVER: SPHERE LINKAGE GRADIENT FIXES
C 29 DEC 00 - BM  - CORRECTED THE GRADIENT BUG
C 26 OCT 00 - BM  - ALL IEF ROUTINES ARE HERE
C
C*MODULE PCMIEF  *DECK IEFCMM
      SUBROUTINE IEFCMM(DMATM1,SE,DE,SI,DI,VERT,CENTR,WORK,IPVT,
     *                  XCTS,YCTS,ZCTS,AS,NVERT,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXVRT=20)
C
      DIMENSION DMATM1(NTS,NTS),SE(NTS,NTS),DE(NTS,NTS),
     *          SI(NTS,NTS),DI(NTS,NTS),VERT(NTS*MXVRT*3),
     *          CENTR(NTS*MXVRT*3),WORK(NTS),IPVT(NTS),XCTS(*),YCTS(*),
     *          ZCTS(*),AS(*),NVERT(*),XE(*),YE(*),ZE(*),RE(*),ISPHE(*)
      DIMENSION DET(2)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /ANIDAT/ EPSM,EPSXX,EPSXY,EPSXZ,EPSYY,EPSYZ,EPSZZ,
     *                EPSM1XX,EPSM1XY,EPSM1XZ,EPSM1YY,EPSM1YZ,EPSM1ZZ,
     *                EPS1,EPS2,EPS3,ROT(3,3)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /IONDAT/ EPSI,DK2,DALP
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      DATA TWO/2.0D+00/
      DATA FPI/12.56637061D+00/
C
C        COMPUTE THE C AND THEN C-INVERSE MATRIX
C
      IF(IPRINT.EQ.1. AND .MASWRK) THEN
         WRITE(IW,*) ' -------------------'
         WRITE(IW,*) ' -- MATRIX C^(-1) --'
         WRITE(IW,*) ' -------------------'
      END IF
C
      IF(IEF.EQ.1) DETEPS=SQRT(EPS1*EPS2*EPS3)
       IF(MASWRK) WRITE(6,*) 'DETEPS',DETEPS,EPSM1XX,EPSM1XY
       IF(MASWRK) WRITE(6,*) EPSM1XZ,EPSM1YY,EPSM1YZ,EPSM1ZZ
C
C         COMPUTE THE C MATRIX, STORED AT -DMATM1-
C
      DO 140 I=1,NTS
C
C -- CALCULATION OF THE DIAGONAL TERMS
C
        LI=ISPHE(I)
        SSI=AS(I)*1.07D+00*SQRT(FPI*AS(I))/FPI
        DDI=-SSI/(2.0D+00*RE(LI))
        SI(I,I)=SSI
        DI(I,I)=DDI
C
        IF(IEF.EQ.1) THEN
C
C   ANISOTROPIC DIELECTRICS
C
          CALL DIAGAN(I,SSE,DDE,VERT,CENTR,XCTS,YCTS,ZCTS,AS,NVERT,
     *                XE,YE,ZE,RE,ISPHE)
          SE(I,I)=SSE/DETEPS
          DE(I,I)=DDE/DETEPS
        ELSE
C
C   IONIC SOLUTIONS
C
          CALL DIAGION(I,SSE,DDE,VERT,CENTR,XCTS,YCTS,ZCTS,AS,
     *                 NVERT,XE,YE,ZE,RE,ISPHE)
          SE(I,I)=SSE
          DE(I,I)=DDE
        END IF
C
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
C
C -- CALCULATION OF THE OFF-DIAGONAL TERMS
C
        DO 150 J=1,NTS
C
          IF (J.EQ.I) GO TO 150
          LJ=ISPHE(J)
C
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          CXJ=(XJ-XE(LJ))/RE(LJ)
          CYJ=(YJ-YE(LJ))/RE(LJ)
          CZJ=(ZJ-ZE(LJ))/RE(LJ)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          DRIJ3=DRIJ**3
C
          SSI=AS(I)*AS(J)/(FPI*DRIJ)
          DDI=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
C
          SI(I,J)=SSI
          DI(I,J)=DDI
        IF(IEF.EQ.1) THEN
          VECEPSX=EPSM1XX*(XI-XJ)+EPSM1XY*(YI-YJ)+EPSM1XZ*(ZI-ZJ)
          VECEPSY=EPSM1XY*(XI-XJ)+EPSM1YY*(YI-YJ)+EPSM1YZ*(ZI-ZJ)
          VECEPSZ=EPSM1XZ*(XI-XJ)+EPSM1YZ*(YI-YJ)+EPSM1ZZ*(ZI-ZJ)
          DRIJEPS=SQRT((XI-XJ)*VECEPSX+(YI-YJ)*VECEPSY+(ZI-ZJ)*VECEPSZ)
          DRIJEPS3=DRIJEPS**3
C
          SSE=AS(I)*AS(J)/(FPI*DRIJEPS)
          DDE=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJEPS3)
C
          SE(I,J)=SSE/DETEPS
          DE(I,J)=DDE/DETEPS
        ELSE
          SSE=AS(I)*AS(J)*EXP(-DALP*DRIJ)/(FPI*DRIJ*EPSI)
          DDE=AS(I)*AS(J)*EXP(-DALP*DRIJ)*(1.0D+00+DALP*DRIJ)
     *            *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
C
         SE(I,J)=SSE
         DE(I,J)=DDE
        END IF
C
  150   CONTINUE
  140 CONTINUE
C
      DO I=1,NTS
      DO K=1,NTS
        DMATP=0.0D+00
        DO J=1,NTS
          DELTAIJ=0.0D+00
          DELTAJK=0.0D+00
          IF (I.EQ.J) DELTAIJ=AS(I)
          IF (J.EQ.K) DELTAJK=AS(K)
          DMATP=DMATP
     *          +(DELTAIJ/TWO-DE(I,J))*(1.0D+00/AS(J))*SI(J,K)
     *          +SE(I,J)*(1.0D+00/AS(J))*(DELTAJK/TWO+DI(K,J))
        ENDDO
        DMATM1(I,K)=FPI*DMATP
      ENDDO
      ENDDO
C
      IF(IPRINT.EQ.1. AND .MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      INFO=0
      CALL DGEFA(DMATM1,NTS,NTS,IPVT,INFO)
      IF(INFO.NE.0.) THEN
         IF(MASWRK) WRITE(IW,*) ' THE A MATRIX IS SINGULAR'
         CALL ABRT
      END IF
C
      CALL DGEDI(DMATM1,NTS,NTS,IPVT,DET,WORK,01)
C
      DO I=1,NTS
      DO K=1,NTS
        DMATM1(I,K)=-DMATM1(I,K)
      ENDDO
      ENDDO
C      IF(MASWRK) WRITE(6,333)(DMATM1(1,NN),NN=1,5)
C      IF(MASWRK) WRITE(6,333)(DMATM1(2,NN),NN=1,5)
C      IF(MASWRK) WRITE(6,333)(DMATM1(3,NN),NN=1,5)
C 333  FORMAT(5(F8.3,1X))
C
      IF(IPRINT.EQ.1. AND .MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A INVERSE MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      RETURN
      END
C*MODULE PCMIEF  *DECK ICVEBF
      SUBROUTINE ICVEBF(DMATM1,SE,DE,BEMCHG,VEC,
     *                  VPOT,VCAM,FLD,FLW,L2,NFT27,SOME,
     *                  XCTS,YCTS,ZCTS,AS,QSN,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMATM1(NTS,NTS),BEMCHG(L2),
     *          SE(NTS,NTS),DE(NTS,NTS),VEC(NTS),
     *          VPOT(NTS),VCAM(NTS),FLD(L2,3),FLW(784,3),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),XE(*),YE(*),ZE(*),
     *          RE(*),ISPHE(*)
C
      LOGICAL SOME
      LOGICAL GOPARR, DSKWRK, MASWRK
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 /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
      COMMON /XYZPRP/ XP,YP,ZP,DMY(35)
C
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      DATA ZERO, ONE, TWO  /0.0D+00, 1.0D+00, 2.0D+00/
      DATA FPI/12.56637061D+00/
      DATA ELFLD  /8HELFLD   /
C
C     AFTER PROJECTION ONTO THE NORMAL TO THE TESSERA,
C     THE ELECTRIC FIELD INTEGRALS ARE WRITTEN TO DISK FILE -NFT27-.
C
C  1) CALCULATION OF < CHI(MU) | ALPHA/R | CHI(NU) > INTEGRALS
C     (ELECT. FIELD INTEGRALS) ON REPRESENTATIVE POINTS OF TESSERAE.
C
      IEFLDOLD=IEFLD
      IEFLD=1
      CALL SEQREW(NFT27)
      DO ITS = 1, NTS
         XP=XCTS(ITS)
         YP=YCTS(ITS)
         ZP=ZCTS(ITS)
         L=ISPHE(ITS)
         CNX=(XP-XE(L))/RE(L)
         CNY=(YP-YE(L))/RE(L)
         CNZ=(ZP-ZE(L))/RE(L)
         CALL PRCALC(ELFLD,FLD,FLW,3,L2,.FALSE.)
         DO IBAS = 1, L2
            BEMCHG(IBAS) = FLD(IBAS,1)*CNX
     *                   + FLD(IBAS,2)*CNY
     *                   + FLD(IBAS,3)*CNZ
         ENDDO
         CALL SQWRIT(NFT27,BEMCHG,L2)
      ENDDO
      CALL SEQREW(NFT27)
      IEFLD=IEFLDOLD
C
C  2) CALCULATION OF APPARENT CHARGES GENERATED BY THE SOLUTE'S NUCLEI.
C
      DO  ITS = 1, NTS
        XI=XCTS(ITS)
        YI=YCTS(ITS)
        ZI=ZCTS(ITS)
        L=ISPHE(ITS)
        CX=(XI-XE(L))/RE(L)
        CY=(YI-YE(L))/RE(L)
        CZ=(ZI-ZE(L))/RE(L)
        VPOT(ITS) = ZERO
        VCAM(ITS) = ZERO
         DO  JATOM = 1, NAT
          XL=C(1,JATOM)
          YL=C(2,JATOM)
          ZL=C(3,JATOM)
          RIL=SQRT((XI-XL)**2+(YI-YL)**2+(ZI-ZL)**2)
          XM=((XI-XL)*CX+(YI-YL)*CY+(ZI-ZL)*CZ)/RIL**3
          VPOT(ITS) = VPOT(ITS) - ZAN(JATOM)/RIL
          VCAM(ITS) = VCAM(ITS) + XM*ZAN(JATOM)
         ENDDO
C
      ENDDO
C
      DO I=1,NTS
        VECGP=0.0D+00
        DO J=1,NTS
          DELTAIJ=0.0D+00
          IF (I.EQ.J) DELTAIJ=AS(I)
          VECGP=VECGP-(DELTAIJ/TWO-DE(I,J))*VPOT(J)
     *               -SE(I,J)*VCAM(J)
        ENDDO
        VEC(I)=VECGP
      ENDDO
C
      QNTOT = ZERO
      QNUC2 = ZERO
      DO ITS = 1, NTS
         QSN(ITS) = ZERO
         DO JTS = 1, NTS
            QSN(ITS) = QSN(ITS)+DMATM1(ITS,JTS)*VEC(JTS)
         ENDDO
         QSN(ITS) = QSN(ITS)*AS(ITS)
         QNUC2 = QNUC2 + VCAM(ITS)*AS(ITS)
         QNTOT = QNTOT + QSN(ITS)
      ENDDO
      QNUC2 = QNUC2/FPI
      IF(SOME. AND .MASWRK) WRITE(IW,1200) QNUC2
C
C  4) RENORMALIZATION OF `NUCLEAR' APPARENT CHARGES.
C     IF ICOMP=1 THE CORRECTION IS PROPORTIONAL TO THE TESSERA'S AREA.
C     IF ICOMP=2 OR ICOMP=3 THE CORRECTION IS EQUAL FOR EACH TESSERA.
C
      CHG = ZERO
      DO JATOM = 1, NAT
         CHG = CHG + ZAN(JATOM)
      ENDDO
C
      TCH = - CHG * (EPS - ONE) / EPS
      QNTOTN = ZERO
C
C     OPTION 1:
C
      FN = ONE
      IF(ICOMP.EQ.1) THEN
         SUPTOT = STOT * ANTOAU * ANTOAU
         DIFF = TCH - QNTOT
         DO ITS = 1, NTS
            QSN(ITS) = QSN(ITS) + DIFF * AS(ITS) / SUPTOT
            QNTOTN = QNTOTN + QSN(ITS)
         ENDDO
C
C     OPTION 2:
C
      ELSE IF(ICOMP.EQ.2) THEN
         FN = TCH / QNTOT
         DO ITS = 1, NTS
            QSN(ITS) = QSN(ITS) * FN
            QNTOTN = QNTOTN + QSN(ITS)
         ENDDO
      END IF
C
C     TOTAL APPARENT CHARGE BEFORE AND AFTER RENORMALIZATION.
C
      IF(SOME. AND .MASWRK) THEN
         IF(ICOMP.NE.0) WRITE(IW,1000) QNTOT, TCH, QNTOTN
         IF(ICOMP.EQ.0) WRITE(IW,1100) QNTOT, TCH
      END IF
      RETURN
C
 1000 FORMAT(1X,'NUCLEAR APPARENT CHARGE',F10.5,' THEORETICAL',
     *          F10.5,' RENORMALIZED',F10.5)
 1100 FORMAT(1X,'NUCLEAR APPARENT CHARGE',F10.5,' THEORETICAL',
     *          F10.5,' NOT RENORMALIZED',F10.5)
 1200 FORMAT(/1X,'ESTIMATE OF NUCLEAR CHARGE',F15.5)
      END
C*MODULE PCMIEF  *DECK IXMATBF
      SUBROUTINE IXMATBF(D,XX,Q,SCR,DMATM1,SE,DE,QPOT,QCHG,SCRPOT,QET,
     *                   QETN,TCH,NUM2,NFT27,XCTS,YCTS,ZCTS,AS,QSN,QSE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      DIMENSION D(NUM2),XX(NUM2),Q(NTS),SCR(NUM2),DMATM1(NTS,NTS),
     *          SE(NTS,NTS),DE(NTS,NTS),QPOT(NTS),QCHG(NTS),
     *          SCRPOT(NUM2),XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),QSE(*)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
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 /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
C
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      DATA ZERO, ONE, TWO/0.0D+00, 1.0D+00, 2.0D+00/
C
C     ROUTINE FOR THE CALCULATION OF MATRIX XX (INTERACTION BETWEEN
C     `ELECTRONIC' APPARENT CHARGES AND ELECTRONIC POTENTIAL OF THE
C     SOLUTE):
C         XX(M,N) = SUM_I[ V(M,N;I)QSE(I) ]
C         QSE(I) = A(I)SUM_J[ DMATM1(I,J){SUM_M,N[ D(M,N)DV(M,N;I)DN ]}
C     WHICH IS USED IN THE SCF ROUTINE IN ORDER TO CORRECT TWO-ELECTRON
C     INTEGRALS.
C
C     DISK FILE NFT27 STORES THE NORMAL COMPONENTS OF THE ELECTRIC
C     FIELD -DV(M,N;I)/DN.
C
      CALL SEQREW(NFT27)
C
      DO ITS = 1, NTS
        XI=XCTS(ITS)
        YI=YCTS(ITS)
        ZI=ZCTS(ITS)
        CALL INTMEP(SCRPOT,XI,YI,ZI)
        CALL SQREAD(NFT27,SCR,NUM2)
        QPOT(ITS) = TRACEP(D,SCRPOT,NUM)
        QCHG(ITS) = TRACEP(D,SCR,NUM)
      ENDDO
C
      DO I=1,NTS
        VECGP =ZERO
C
        DO J=1,NTS
          DELTAIJ=ZERO
          IF(I.EQ.J)DELTAIJ=AS(I)
C         IF(I.EQ.J)DELTAIJ=1.0D+00
          VECGP=VECGP-(DELTAIJ/TWO-DE(I,J))*QPOT(J)
     *               - SE(I,J)*QCHG(J)
        ENDDO
C
        Q(I) =VECGP
      ENDDO
C
      CALL SEQREW(NFT27)
C
      QET = ZERO
      DO ITS = 1, NTS
         QSE(ITS) = ZERO
         DO ITSJ = 1, NTS
            QSE(ITS) = QSE(ITS)+DMATM1(ITS,ITSJ)*Q(ITSJ)
         ENDDO
         QSE(ITS) = QSE(ITS)*AS(ITS)
         QET = QET + QSE(ITS)
      ENDDO
C
C  4) NORMALIZZAZIONE DELLE CARICHE VIRTUALI ELEMENTARI ELETTRONICHE.
C     PER ICOMP=1 LA CORREZIONE E' DISTRIBUITA IN PROPORZIONE
C     ALL'AREA DELLA TESSERA.
C     PER ICOMP=2 LA CORREZIONE E' EFFETTUATA CON FATTORI COSTANTI.
C
      DNE = NE
      TCH = DNE * (EPS - ONE) / EPS
C
C     OPTION 1
C
      QETN = ZERO
      FE = ONE
      IF(ICOMP.EQ.1) THEN
         SUPTOT = STOT * ANTOAU * ANTOAU
         DIFF = TCH - QET
         DO ITS = 1, NTS
            QSE(ITS) = QSE(ITS) + DIFF * AS(ITS) / SUPTOT
            QETN = QETN + QSE(ITS)
        ENDDO
C
C     OPTION 2
C
      ELSE IF(ICOMP.EQ.2) THEN
         FE = TCH / QET
         DO ITS = 1, NTS
            QSE(ITS) = QSE(ITS) * FE
            QETN = QETN + QSE(ITS)
         ENDDO
      END IF
C
C     PRINT TOTAL INDUCED CHARGE BEFORE AND AFTER RENORMALIZATION
C
      IF(NPRINT.EQ.5 .AND. MASWRK) THEN
         IF(ICOMP.NE.0) WRITE(IW,1000) QET, TCH, QETN
         IF(ICOMP.EQ.0) WRITE(IW,1100) QET, TCH
      END IF
C
C   CALCULATION OF INTERACTION MATRICES TO BE ADDED TO THE FOCK MATRIX
C
      CALL VCLR(XX,1,NUM2)
      PB = ZERO
      PC = ZERO
      PX = ZERO
C
      DO ITS = 1, NTS
        CALL INTMEP(SCR,XCTS(ITS),YCTS(ITS),ZCTS(ITS))
        DO IBAS = 1, NUM2
          XX(IBAS) = XX(IBAS) - SCR(IBAS) * QSE(ITS)
        ENDDO
C
C     PB = INTERACTION ELECTRONS-NUCLEAR INDUCED CHARGES
C     PX = INTERACTION ELECTRONS-ELECTRONIC INDUCED CHARGES
C
        VEL =  - TRACEP(D,SCR,NUM)
        PB = PB + VEL * QSN(ITS)
        PX = PX + VEL * QSE(ITS)
C
C     PC = INTERACTION NUCLEI-ELECTRONIC INDUCED CHARGES
C
        DO JATOM = 1, NAT
          R2 = (C(1,JATOM)-XCTS(ITS))**2 + (C(2,JATOM)-
     *         YCTS(ITS))**2 + (C(3,JATOM)-ZCTS(ITS))**2
          R = SQRT(R2)
          PC = PC + QSE(ITS) * ZAN(JATOM) / R
        ENDDO
      ENDDO
      RETURN
C
 1000 FORMAT(1X,'ELECTRONIC APPARENT CHARGE',F10.5,' THEORETICAL ',
     * F10.5,' RENORMALIZED',F10.5)
 1100 FORMAT(1X,'ELECTRONIC APPARENT CHARGE',F10.5,' THEORETICAL ',
     * F10.5,' NOT RENORMALIZED',F10.5)
      END
C*MODULE PCMIEF  *DECK DIAGAN
      SUBROUTINE DIAGAN(I,SSE,DDE,VERT,CENTR,XCTS,YCTS,ZCTS,AS,
     *                  NVERT,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXVRT=20)
C
      DIMENSION VERT(NTS*MXVRT*3),CENTR(NTS*MXVRT*3)
      DIMENSION NUMB(MXVRT),THETA(MXVRT),PHI(MXVRT),PHINUMB(MXVRT),
     *          XGP16PTS(8),WGP16PTS(8), XGP64PTS(32),WGP64PTS(32),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),NVERT(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
C
      COMMON /ANIDAT/ EPSM,EPSXX,EPSXY,EPSXZ,EPSYY,EPSYZ,EPSZZ,
     *                EPSM1XX,EPSM1XY,EPSM1XZ,EPSM1YY,EPSM1YZ,EPSM1ZZ,
     *                EPS1,EPS2,EPS3,ROT(3,3)
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      DATA ZERO,TWO/0.0D+00,2.0D+00/
      DATA TPI,FPI/6.28318531D+00,12.56637061D+00/
C
      DATA NG16PTS/8/
      DATA XGP16PTS/9.894009350D-01,9.445750231D-01,
     *              8.656312024D-01,7.554044084D-01,
     *              6.178762444D-01,4.580167777D-01,
     *              2.816035508D-01,9.501250984D-02/
      DATA WGP16PTS/2.715245941D-02,6.225352394D-02,
     *              9.515851168D-02,1.246289713D-01,
     *              1.495959888D-01,1.691565194D-01,
     *              1.826034150D-01,1.894506105D-01/
C
      DATA NG64PTS/32/
      DATA XGP64PTS/9.993050417D-01,9.963401168D-01,
     *              9.910133715D-01,9.833362539D-01,
     *              9.733268278D-01,9.610087996D-01,
     *              9.464113748D-01,9.295691721D-01,
     *              9.105221371D-01,8.893154460D-01,
     *              8.659993982D-01,8.406292963D-01,
     *              8.132653151D-01,7.839723589D-01,
     *              7.528199073D-01,7.198818502D-01,
     *              6.852363131D-01,6.489654713D-01,
     *              6.111553552D-01,5.718956462D-01,
     *              5.312794640D-01,4.894031457D-01,
     *              4.463660173D-01,4.022701580D-01,
     *              3.572201583D-01,3.113228720D-01,
     *              2.646871622D-01,2.174236437D-01,
     *              1.696444204D-01,1.214628193D-01,
     *              7.299312179D-02,2.435029266D-02/
      DATA WGP64PTS/1.783280722D-03,4.147033261D-03,
     *              6.504457969D-03,8.846759826D-03,
     *              1.116813946D-02,1.346304790D-02,
     *              1.572603048D-02,1.795171578D-02,
     *              2.013482315D-02,2.227017381D-02,
     *              2.435270257D-02,2.637746972D-02,
     *              2.833967261D-02,3.023465707D-02,
     *              3.205792835D-02,3.380516184D-02,
     *              3.547221326D-02,3.705512854D-02,
     *              3.855015318D-02,3.995374113D-02,
     *              4.126256324D-02,4.247351512D-02,
     *              4.358372453D-02,4.459055816D-02,
     *              4.549162793D-02,4.628479658D-02,
     *              4.696818282D-02,4.754016571D-02,
     *              4.799938860D-02,4.834476223D-02,
     *              4.857546744D-02,4.869095701D-02/
C
C
      SSE=ZERO
      DDE=ZERO
      LI=ISPHE(I)
      NV=30*(I-1)
C
C -- LOOP ON INTEGRATION POINTS
C
        SSEP=0.0D+00
        DDEP=0.0D+00
        XZ=(XCTS(I)-XE(LI))/RE(LI)
        YZ=(YCTS(I)-YE(LI))/RE(LI)
        ZZ=(ZCTS(I)-ZE(LI))/RE(LI)
        RMIN=0.99D+00
        IF (ABS(XZ).LE.RMIN) THEN
          RMIN=ABS(XZ)
          XX=0.0D+00
          YX=-ZZ/SQRT(1.0D+00-XZ*XZ)
          ZX=YZ/SQRT(1.0D+00-XZ*XZ)
        END IF
        IF (ABS(YZ).LE.RMIN) THEN
          RMIN=ABS(YZ)
          XX=ZZ/SQRT(1.0D+00-YZ*YZ)
          YX=0.0D+00
          ZX=-XZ/SQRT(1.0D+00-YZ*YZ)
        END IF
        IF (ABS(ZZ).LE.RMIN) THEN
          XX=YZ/SQRT(1.0D+00-ZZ*ZZ)
          YX=-XZ/SQRT(1.0D+00-ZZ*ZZ)
          ZX=0.0D+00
        END IF
        XY=YZ*ZX-YX*ZZ
        YY=ZZ*XX-ZX*XZ
        ZY=XZ*YX-XX*YZ
C
        DO K=1,NVERT(I)
          VECX=VERT(NV+3*(K-1)+1)-XE(LI)
          VECY=VERT(NV+3*(K-1)+2)-YE(LI)
          VECZ=VERT(NV+3*(K-1)+3)-ZE(LI)
          DC=(VECX*XZ+VECY*YZ+VECZ*ZZ)/RE(LI)
          IF(DC.GE.1.0D+00)DC=1.0D+00
          IF(DC.LE.-1.0D+00)DC=-1.0D+00
          THETA(K)=ACOS(DC)
          DCS=(VECX*XX+VECY*YX+VECZ*ZX)
     *        /(RE(LI)*SIN(THETA(K)))
          IF(DCS.GE.1.0D+00)DCS=1.0D+00
          IF(DCS.LE.-1.0D+00)DCS=-1.0D+00
          PHI(K)=ACOS(DCS)
          SNPHI=(VECX*XY+VECY*YY+VECZ*ZY)/(RE(LI)*SIN(THETA(K)))
          IF (SNPHI.LE.0.0D+00) PHI(K)=TPI-PHI(K)
        ENDDO
C
        DO K=2,NVERT(I)
          PHI(K)=PHI(K)-PHI(1)
          IF (PHI(K).LT.0.0D+00) PHI(K)=TPI+PHI(K)
        ENDDO
C
        XX=XX*COS(PHI(1))+XY*SIN(PHI(1))
        YX=YX*COS(PHI(1))+YY*SIN(PHI(1))
        ZX=ZX*COS(PHI(1))+ZY*SIN(PHI(1))
        XY=YZ*ZX-YX*ZZ
        YY=ZZ*XX-ZX*XZ
        ZY=XZ*YX-XX*YZ
C
        PHI(1)=0.0D+00
        NUMB(1)=1
        NUMB(2)=2
        PHINUMB(1)=PHI(1)
        PHINUMB(2)=PHI(2)
C
        DO 210 K=3,NVERT(I)
          DO L=2,K-1
            IF (PHI(K).LT.PHINUMB(L)) THEN
              DO M=1,K-L
                NUMB(K-M+1)=NUMB(K-M)
                PHINUMB(K-M+1)=PHINUMB(K-M)
              ENDDO
              NUMB(L)=K
              PHINUMB(L)=PHI(K)
              GO TO 210
            END IF
          ENDDO
          NUMB(K)=K
          PHINUMB(K)=PHI(K)
  210   CONTINUE
        NUMB(NVERT(I)+1)=NUMB(1)
        PHINUMB(NVERT(I)+1)=TPI
C
C -- LOOP ON GAUSS POINTS
C
        DO NEDGE=1,NVERT(I)
C
          APH=(PHINUMB(NEDGE+1)-PHINUMB(NEDGE))/2.0D+00
          BPH=(PHINUMB(NEDGE+1)+PHINUMB(NEDGE))/2.0D+00
          THA=THETA(NUMB(NEDGE))
          THB=THETA(NUMB(NEDGE+1))
          PHA=PHINUMB(NEDGE)
          PHB=PHINUMB(NEDGE+1)
          OCX=(CENTR(NV+3*(NEDGE-1)+1)-XE(LI))/RE(LI)
          OCY=(CENTR(NV+3*(NEDGE-1)+2)-YE(LI))/RE(LI)
          OCZ=(CENTR(NV+3*(NEDGE-1)+3)-ZE(LI))/RE(LI)
          RCC=OCX**2+OCY**2+OCZ**2
C
C
          DO NPH=1,NG64PTS
          DO NPHSGN=0,1
C
           PH=(2*NPHSGN-1)*APH*XGP64PTS(NPH)+BPH
           COSPH=COS(PH)
           SINPH=SIN(PH)
C
          IF (RCC.LT.1.0D-07) THEN
           COTGTHMAX=(SIN(PH-PHA)/TAN(THB)+SIN(PHB-PH)/TAN(THA))
     *               /SIN(PHB-PHA)
           THMAX=ATAN(1.0D+00/COTGTHMAX)
          ELSE
           XCC=XX*OCX+YX*OCY+ZX*OCZ
           YCC=XY*OCX+YY*OCY+ZY*OCZ
           ZCC=XZ*OCX+YZ*OCY+ZZ*OCZ
           AA=(XCC*COSPH+YCC*SINPH)**2+ZCC**2
           BB=-ZCC*RCC
           CC=RCC**2-(XCC*COSPH+YCC*SINPH)**2
           DS=BB**2-AA*CC
           IF(DS.LT.ZERO)DS=ZERO
           CS=(-BB+SQRT(DS))/AA
           IF(CS.GT.1.0D+00)CS=1.0D+00
           IF(CS.LT.-1.0D+00)CS=1.0D+00
           THMAX=ACOS(CS)
          END IF
           IF(THMAX.LT.1.0D-08) GO TO 33
C
           ATH=THMAX/TWO
C
           SSEPP=0.0D+00
           DDEPP=0.0D+00
C
           DO NTH=1,NG16PTS
           DO NTHSGN=0,1
C
            TH=(2*NTHSGN-1)*ATH*XGP16PTS(NTH)+ATH
            COSTH=COS(TH)
            SINTH=SIN(TH)
            VX=XX*SINTH*COSPH+XY*SINTH*SINPH+XZ*(COSTH-1.0D+00)
            VY=YX*SINTH*COSPH+YY*SINTH*SINPH+YZ*(COSTH-1.0D+00)
            VZ=ZX*SINTH*COSPH+ZY*SINTH*SINPH+ZZ*(COSTH-1.0D+00)
            RTH=SQRT(2*(1-COSTH))
            RTHEPS=SQRT(EPSM1XX*VX*VX+2*EPSM1XY*VX*VY+2*EPSM1XZ*VX*VZ
     *                               +  EPSM1YY*VY*VY+2*EPSM1YZ*VY*VZ
     *                                               +  EPSM1ZZ*VZ*VZ)
            SSEPP=SSEPP+(RE(LI)/(FPI*RTHEPS))*SINTH*ATH*WGP16PTS(NTH)
            DDEPP=DDEPP-(RTH**2/(2*FPI*RTHEPS**3))*SINTH
     *                   *ATH*WGP16PTS(NTH)
C
           ENDDO
           ENDDO
C
           SSEP=SSEP+SSEPP*APH*WGP64PTS(NPH)
           DDEP=DDEP+DDEPP*APH*WGP64PTS(NPH)
  33       CONTINUE
C
          ENDDO
          ENDDO
C
        ENDDO
C
        SSE=SSEP*AS(I)
        DDE=DDEP*AS(I)
C
      RETURN
      END
C*MODULE PCMIEF  *DECK DIAGION
      SUBROUTINE DIAGION(I,SSE,DDE,VERT,CENTR,XCTS,YCTS,ZCTS,AS,
     *                   NVERT,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXVRT=20)
C
      DIMENSION VERT(NTS*MXVRT*3),CENTR(NTS*MXVRT*3)
      DIMENSION NUMB(MXVRT),THETA(MXVRT),PHI(MXVRT),PHINUMB(MXVRT),
     *          XGP16PTS(8),WGP16PTS(8),XGP64PTS(32),WGP64PTS(32),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),NVERT(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
C
      COMMON /IONDAT/ EPSI,DK2,DALP
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      DATA ZERO,TWO/0.0D+00,2.0D+00/
      DATA TPI,FPI/6.28318531D+00,12.56637061D+00/
C
      DATA NG16PTS/8/
      DATA XGP16PTS/9.894009350D-01,9.445750231D-01,
     *              8.656312024D-01,7.554044084D-01,
     *              6.178762444D-01,4.580167777D-01,
     *              2.816035508D-01,9.501250984D-02/
      DATA WGP16PTS/2.715245941D-02,6.225352394D-02,
     *              9.515851168D-02,1.246289713D-01,
     *              1.495959888D-01,1.691565194D-01,
     *              1.826034150D-01,1.894506105D-01/
C
      DATA NG64PTS/32/
      DATA XGP64PTS/9.993050417D-01,9.963401168D-01,
     *              9.910133715D-01,9.833362539D-01,
     *              9.733268278D-01,9.610087996D-01,
     *              9.464113748D-01,9.295691721D-01,
     *              9.105221371D-01,8.893154460D-01,
     *              8.659993982D-01,8.406292963D-01,
     *              8.132653151D-01,7.839723589D-01,
     *              7.528199073D-01,7.198818502D-01,
     *              6.852363131D-01,6.489654713D-01,
     *              6.111553552D-01,5.718956462D-01,
     *              5.312794640D-01,4.894031457D-01,
     *              4.463660173D-01,4.022701580D-01,
     *              3.572201583D-01,3.113228720D-01,
     *              2.646871622D-01,2.174236437D-01,
     *              1.696444204D-01,1.214628193D-01,
     *              7.299312179D-02,2.435029266D-02/
      DATA WGP64PTS/1.783280722D-03,4.147033261D-03,
     *              6.504457969D-03,8.846759826D-03,
     *              1.116813946D-02,1.346304790D-02,
     *              1.572603048D-02,1.795171578D-02,
     *              2.013482315D-02,2.227017381D-02,
     *              2.435270257D-02,2.637746972D-02,
     *              2.833967261D-02,3.023465707D-02,
     *              3.205792835D-02,3.380516184D-02,
     *              3.547221326D-02,3.705512854D-02,
     *              3.855015318D-02,3.995374113D-02,
     *              4.126256324D-02,4.247351512D-02,
     *              4.358372453D-02,4.459055816D-02,
     *              4.549162793D-02,4.628479658D-02,
     *              4.696818282D-02,4.754016571D-02,
     *              4.799938860D-02,4.834476223D-02,
     *              4.857546744D-02,4.869095701D-02/
C
      SSE=ZERO
      DDE=ZERO
      LI=ISPHE(I)
      NV=30*(I-1)
C
C -- LOOP ON INTEGRATION POINTS
C
        SSEP=0.0D+00
        DDEP=0.0D+00
        XZ=(XCTS(I)-XE(LI))/RE(LI)
        YZ=(YCTS(I)-YE(LI))/RE(LI)
        ZZ=(ZCTS(I)-ZE(LI))/RE(LI)
        RMIN=0.99D+00
        IF (ABS(XZ).LE.RMIN) THEN
          RMIN=ABS(XZ)
          XX=0.0D+00
          YX=-ZZ/SQRT(1.0D+00-XZ*XZ)
          ZX=YZ/SQRT(1.0D+00-XZ*XZ)
        END IF
        IF (ABS(YZ).LE.RMIN) THEN
          RMIN=ABS(YZ)
          XX=ZZ/SQRT(1.0D+00-YZ*YZ)
          YX=0.0D+00
          ZX=-XZ/SQRT(1.0D+00-YZ*YZ)
        END IF
        IF (ABS(ZZ).LE.RMIN) THEN
          XX=YZ/SQRT(1.0D+00-ZZ*ZZ)
          YX=-XZ/SQRT(1.0D+00-ZZ*ZZ)
          ZX=0.0D+00
        END IF
        XY=YZ*ZX-YX*ZZ
        YY=ZZ*XX-ZX*XZ
        ZY=XZ*YX-XX*YZ
C
        DO K=1,NVERT(I)
          VECX=VERT(NV+3*(K-1)+1)-XE(LI)
          VECY=VERT(NV+3*(K-1)+2)-YE(LI)
          VECZ=VERT(NV+3*(K-1)+3)-ZE(LI)
          DC=(VECX*XZ+VECY*YZ+VECZ*ZZ)/RE(LI)
          IF(DC.GE.1.0D+00)DC=1.0D+00
          IF(DC.LE.-1.0D+00)DC=-1.0D+00
          THETA(K)=ACOS(DC)
          DCS=(VECX*XX+VECY*YX+VECZ*ZX)
     *        /(RE(LI)*SIN(THETA(K)))
          IF(DCS.GE.1.0D+00)DCS=1.0D+00
          IF(DCS.LE.-1.0D+00)DCS=-1.0D+00
          PHI(K)=ACOS(DCS)
          SNPHI=(VECX*XY+VECY*YY+VECZ*ZY)/(RE(LI)*SIN(THETA(K)))
          IF (SNPHI.LE.0.0D+00) PHI(K)=TPI-PHI(K)
        ENDDO
C
        DO K=2,NVERT(I)
          PHI(K)=PHI(K)-PHI(1)
          IF (PHI(K).LT.0.0D+00) PHI(K)=TPI+PHI(K)
        ENDDO
C
        XX=XX*COS(PHI(1))+XY*SIN(PHI(1))
        YX=YX*COS(PHI(1))+YY*SIN(PHI(1))
        ZX=ZX*COS(PHI(1))+ZY*SIN(PHI(1))
        XY=YZ*ZX-YX*ZZ
        YY=ZZ*XX-ZX*XZ
        ZY=XZ*YX-XX*YZ
C
        PHI(1)=0.0D+00
        NUMB(1)=1
        NUMB(2)=2
        PHINUMB(1)=PHI(1)
        PHINUMB(2)=PHI(2)
C
        DO 210 K=3,NVERT(I)
          DO L=2,K-1
            IF (PHI(K).LT.PHINUMB(L)) THEN
              DO M=1,K-L
                NUMB(K-M+1)=NUMB(K-M)
                PHINUMB(K-M+1)=PHINUMB(K-M)
              ENDDO
              NUMB(L)=K
              PHINUMB(L)=PHI(K)
              GO TO 210
            END IF
          ENDDO
          NUMB(K)=K
          PHINUMB(K)=PHI(K)
  210   CONTINUE
        NUMB(NVERT(I)+1)=NUMB(1)
        PHINUMB(NVERT(I)+1)=TPI
C
C -- LOOP ON GAUSS POINTS
C
        DO NEDGE=1,NVERT(I)
C
          APH=(PHINUMB(NEDGE+1)-PHINUMB(NEDGE))/2.0D+00
          BPH=(PHINUMB(NEDGE+1)+PHINUMB(NEDGE))/2.0D+00
          THA=THETA(NUMB(NEDGE))
          THB=THETA(NUMB(NEDGE+1))
          PHA=PHINUMB(NEDGE)
          PHB=PHINUMB(NEDGE+1)
          OCX=(CENTR(NV+3*(NEDGE-1)+1)-XE(LI))/RE(LI)
          OCY=(CENTR(NV+3*(NEDGE-1)+2)-YE(LI))/RE(LI)
          OCZ=(CENTR(NV+3*(NEDGE-1)+3)-ZE(LI))/RE(LI)
          RCC=OCX**2+OCY**2+OCZ**2
C
C
          DO NPH=1,NG64PTS
          DO NPHSGN=0,1
C
           PH=(2*NPHSGN-1)*APH*XGP64PTS(NPH)+BPH
           COSPH=COS(PH)
           SINPH=SIN(PH)
C
          IF (RCC.LT.1.0D-07) THEN
           COTGTHMAX=(SIN(PH-PHA)/TAN(THB)+SIN(PHB-PH)/TAN(THA))
     *               /SIN(PHB-PHA)
           THMAX=ATAN(1.0D+00/COTGTHMAX)
          ELSE
           XCC=XX*OCX+YX*OCY+ZX*OCZ
           YCC=XY*OCX+YY*OCY+ZY*OCZ
           ZCC=XZ*OCX+YZ*OCY+ZZ*OCZ
           AA=(XCC*COSPH+YCC*SINPH)**2+ZCC**2
           BB=-ZCC*RCC
           CC=RCC**2-(XCC*COSPH+YCC*SINPH)**2
           DS=BB**2-AA*CC
           IF(DS.LT.ZERO)DS=ZERO
           CS=(-BB+SQRT(DS))/AA
           IF(CS.GT.1.0D+00)CS=1.0D+00
           IF(CS.LT.-1.0D+00)CS=1.0D+00
           THMAX=ACOS(CS)
          END IF
           IF(THMAX.LT.1.0D-08) GO TO 33
C
           ATH=THMAX/TWO
C
           SSEPP=0.0D+00
           DDEPP=0.0D+00
C
           DO NTH=1,NG16PTS
           DO NTHSGN=0,1
C
            TH=(2*NTHSGN-1)*ATH*XGP16PTS(NTH)+ATH
            COSTH=COS(TH)
            SINTH=SIN(TH)
            RTH=SQRT(2*(1-COSTH))
            SSEPP=SSEPP+(RE(LI)*EXP(-DALP*RE(LI)*RTH)
     *                   /(FPI*RTH*EPSI))*SINTH*ATH*WGP16PTS(NTH)
            DDEPP=DDEPP-(RTH**2*EXP(-DALP*RE(LI)*RTH)
     *                   *(1.0D+00+DALP*RE(LI)*RTH)
     *                   /(2*FPI*RTH**3))*SINTH*ATH*WGP16PTS(NTH)
C
           ENDDO
           ENDDO
C
           SSEP=SSEP+SSEPP*APH*WGP64PTS(NPH)
           DDEP=DDEP+DDEPP*APH*WGP64PTS(NPH)
  33       CONTINUE
C
          ENDDO
          ENDDO
C
        ENDDO
C
        SSE=SSEP*AS(I)
        DDE=DDEP*AS(I)
C
      RETURN
      END
C*MODULE PCMIEF  *DECK NEWC
      SUBROUTINE NEWC(DMFR,IPVT,WORK,SE,DE,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,
     *                RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMFR(NTS,NTS),IPVT(NTS),WORK(NTS),SE(NTS,NTS),
     *          DE(NTS,NTS),XCTS(*),YCTS(*),ZCTS(*),AS(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
      DATA TWO/2.0D+00/
      DATA FPI/12.56637061D+00/
C
      IF(MASWRK) WRITE(IW,*) '---------------------------'
      IF(MASWRK) WRITE(IW,*) '-- MATRIX C^(-1)[EPS(W)] --'
      IF(MASWRK) WRITE(IW,*) '---------------------------'
C
      IF(IEF.LT.3) THEN
C
      DO I=1,NTS
      DO K=1,NTS
        DMATP=0.0D+00
        DO J=1,NTS
          DELTAIJ=0.0D+00
          DELTAJK=0.0D+00
          SEIJ=SE(I,J)*EPS/EPSINF
          SIJK=SE(J,K)*EPS
          IF (I.EQ.J) DELTAIJ=AS(I)
          IF (J.EQ.K) DELTAJK=AS(K)
          DMATP=DMATP
     *          +(DELTAIJ/TWO-DE(I,J))*(1.0D+00/AS(J))*SIJK
     *          +SEIJ*(1.0D+00/AS(J))*(DELTAJK/TWO+DE(K,J))
        ENDDO
        DMFR(I,K)=FPI*DMATP
      ENDDO
      ENDDO
C
      ELSE
C
      DO 140 I=1,NTS
        LI=ISPHE(I)
        SSI=AS(I)*1.07D+00*SQRT(FPI*AS(I))/FPI
        DDI=-SSI/(2.0D+00*RE(LI))
        SE(I,I)=SSI
        DE(I,I)=DDI
C
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
C
        DO 150 J=1,NTS
C
          IF (J.EQ.I) GO TO 150
          LJ=ISPHE(J)
C
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          CXJ=(XJ-XE(LJ))/RE(LJ)
          CYJ=(YJ-YE(LJ))/RE(LJ)
          CZJ=(ZJ-ZE(LJ))/RE(LJ)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          DRIJ3=DRIJ**3
C
          SSI=AS(I)*AS(J)/(FPI*DRIJ)
          DDI=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
C
          SE(I,J)=SSI
          DE(I,J)=DDI
C
  150   CONTINUE
  140 CONTINUE
C
      FACT=(EPSINF+1.0D+00)/(EPSINF-1.0D+00)
      DO I=1,NTS
      DO K=1,NTS
        DMATP=0.0D+00
        DO J=1,NTS
          DELTAIJ=0.0D+00
          IF (I.EQ.J) DELTAIJ=AS(I)*FACT
          DMATP=DMATP
     *          +(DELTAIJ/TWO-DE(I,J))*(1.0D+00/AS(J))*SE(J,K)
        ENDDO
        DMFR(I,K)=FPI*DMATP
      ENDDO
      ENDDO
C
      END IF
C
      INFO=0
      CALL DGEFA(DMFR,NTS,NTS,IPVT,INFO)
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'C MATRIX IS SINGULAR'
         CALL ABRT
      END IF
C
      CALL DGEDI(DMFR,NTS,NTS,IPVT,DET,WORK,01)
C
      IF(IEF.LT.3) THEN
C
      DO I=1,NTS
      DO K=1,NTS
        DMFR(I,K)=-DMFR(I,K)
      ENDDO
      ENDDO
C
      ELSE
C
      DO I=1,NTS
      DO K=1,NTS
        DMATP=0.0D+00
        DO J=1,NTS
          DELTAJK=0.0D+00
          IF (K.EQ.J) DELTAJK=AS(J)
          DMATP=DMATP
     *          -DMFR(I,J)*(DELTAJK/TWO-DE(J,K))
        ENDDO
        SE(I,K)=DMATP
      ENDDO
      ENDDO
C
      DO I=1,NTS
      DO K=1,NTS
        DMFR(I,K)= SE(I,K)
      ENDDO
      ENDDO
C
      END IF
C
      IF(MASWRK) WRITE(IW,*)
     *       '..... DONE GENERATING NEW C INVERSE MATRIX .....'
      CALL TIMIT(1)
      RETURN
      END
C*MODULE PCMIEF  *DECK DMQEXTF
      SUBROUTINE DMQEXTF(DMATM1,DE,IWFLA,IPVT,WORK,XCTS,YCTS,ZCTS,
     *                   AS,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMATM1(NTS,NTS),DE(NTS,NTS),IPVT(NTS),
     *          WORK(NTS),XCTS(*),YCTS(*),ZCTS(*),AS(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
      DIMENSION DET(2)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
      DATA ONE,TWO/1.0D+00,2.0D+00/
      DATA FPI/12.56637061D+00/
C
      IF(IWFLA.EQ.0) THEN
        IF(MASWRK) WRITE(IW,*) 'CASO STATICO: EPS',EPS
        FACT=(EPS+ONE)/(EPS-ONE)
      ELSE
        IF(MASWRK) WRITE(IW,*) 'CASO DINAMICO: EPS',EPSINF
        FACT=(EPSINF+ONE)/(EPSINF-ONE)
      END IF
C
C         COMPUTE THE DE MATRIX
C
      IF(IEF.EQ.3) THEN
C
      DO 140 I=1,NTS
C
C -- CALCULATION OF THE DIAGONAL TERMS
C
        LI=ISPHE(I)
        SSI=AS(I)*1.07D+00*SQRT(FPI*AS(I))/FPI
        DDI=-SSI/(2.0D+00*RE(LI))
        DE(I,I)=DDI
C
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
C
C -- CALCULATION OF THE OFF-DIAGONAL TERMS
C
        DO 150 J=1,NTS
C
          IF (J.EQ.I) GO TO 150
          LJ=ISPHE(J)
C
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          CXJ=(XJ-XE(LJ))/RE(LJ)
          CYJ=(YJ-YE(LJ))/RE(LJ)
          CZJ=(ZJ-ZE(LJ))/RE(LJ)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          DRIJ3=DRIJ**3
C
          DDI=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
C
         DE(I,J)=DDI
C
  150   CONTINUE
  140 CONTINUE
C
      END IF
C
      DO I=1,NTS
        DO J=1,NTS
         DELTAIJ=0.0D+00
         IF (I.EQ.J) DELTAIJ=AS(I)*FACT
         DMATM1(I,J)=FPI*(DELTAIJ/TWO-DE(J,I))*(1.0D+00/AS(I))
        ENDDO
      ENDDO
C
      INFO=0
      CALL DGEFA(DMATM1,NTS,NTS,IPVT,INFO)
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) ' THE A MATRIX IS SINGULAR'
         CALL ABRT
      END IF
C
      CALL DGEDI(DMATM1,NTS,NTS,IPVT,DET,WORK,01)
C
      DO I=1,NTS
      DO K=1,NTS
        DMATM1(I,K)=-DMATM1(I,K)
      ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE PCMIEF  *DECK IEFCMMV
      SUBROUTINE IEFCMMV(DMATM1,SI,DI,WORK,IPVT,DDIMT1,DDIMT2,
     *                   XCTS,YCTS,ZCTS,AS,VAD,QOR,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMATM1(NTS,NTS),SI(NTS,NTS),DI(NTS,NTS),
     *          WORK(NTS),IPVT(NTS),DDIMT1(NTS),DDIMT2(NTS),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),VAD(*),QOR(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
C
      DIMENSION DET(2)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
      DATA ZERO,TWO/0.0D+00,2.0D+00/
      DATA FPI/12.56637061D+00/
C
C        COMPUTE THE C AND THEN C-INVERSE MATRIX
C
         IF(MASWRK) WRITE(IW,*) ' -------------------'
         IF(MASWRK) WRITE(IW,*) ' -- MATRIX C^(-1) --'
         IF(MASWRK) WRITE(IW,*) ' -------------------'
C
C   COMPUTE THE C MATRIX, STORED AT -DMATM1-
C   ISOTROPIC DIELECTRICS
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 1: COMPUTING MATRIX SI AND DI'
C
      DO 140 I=1,NTS
C
C -- CALCULATION OF THE DIAGONAL TERMS
C
        LI=ISPHE(I)
        SSI=AS(I)*1.07D+00*SQRT(FPI*AS(I))/FPI
        DDI=-SSI/(2.0D+00*RE(LI))
        SI(I,I)=SSI
        DI(I,I)=DDI
C
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
C
C -- CALCULATION OF THE OFF-DIAGONAL TERMS
C
        DO 150 J=1,NTS
C
          IF (J.EQ.I) GO TO 150
          LJ=ISPHE(J)
C
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          CXJ=(XJ-XE(LJ))/RE(LJ)
          CYJ=(YJ-YE(LJ))/RE(LJ)
          CZJ=(ZJ-ZE(LJ))/RE(LJ)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          DRIJ3=DRIJ**3
C
          SSI=AS(I)*AS(J)/(FPI*DRIJ)
          DDI=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
C
          SI(I,J)=SSI
          DI(I,J)=DDI
C
  150   CONTINUE
  140 CONTINUE
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 2: COMPUTING MATRIX A '
C
      FACT=(EPS+1.0D+00)/(EPS-1.0D+00)
      DO 200 I=1,NTS
C
C           INITIALIZE PARALLEL
      IPCOUNT = ME - 1
      CALL VCLR(DDIMT1,1,NTS)
C
      DO 210 K=1,NTS
        DMATP=0.0D+00
C            GO PARALLEL!
        IF(GOPARR) THEN
          IPCOUNT = IPCOUNT + 1
          IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 210
        END IF
C
        DO 220 J=1,NTS
          DELTAIJ=0.0D+00
          IF (I.EQ.J) DELTAIJ=AS(I)*FACT
          DMATP=DMATP
     *          +(DELTAIJ/TWO-DI(I,J))*(1.0D+00/AS(J))*SI(J,K)
 220  CONTINUE
C        DMATM1(I,K)=FPI*DMATP
         DDIMT1(K)=FPI*DMATP
 210  CONTINUE
C          SUM UP
        IF(GOPARR) THEN
          CALL DDI_GSUMF(2405,DDIMT1,NTS)
        END IF
           DO 250 J=1,NTS
           DMATM1(I,J)=DDIMT1(J)
 250  CONTINUE
C
 200  CONTINUE
C
      IF(IPRINT.EQ.1 .AND. MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 3: COMPUTING MATRIX A^(-1) '
C
      INFO=0
      CALL DGEFA(DMATM1,NTS,NTS,IPVT,INFO)
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) ' THE A MATRIX IS SINGULAR'
         CALL ABRT
      END IF
C
      CALL DGEDI(DMATM1,NTS,NTS,IPVT,DET,WORK,01)
C
      IF(IPRINT.EQ.1 .AND. MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A INVERSE MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 4: COMPUTING MATRIX C^(-1)'
      IF(MASWRK) WRITE(IW,*) ' '
C
      DO 400 I=1,NTS
C
C           INITIALIZE PARALLEL
      IPCOUNT = ME - 1
      CALL VCLR(DDIMT2,1,NTS)
C
      DO 410 K=1,NTS
        DMATP=0.0D+00
C            GO PARALLEL!
        IF(GOPARR) THEN
          IPCOUNT = IPCOUNT + 1
          IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 410
        END IF
C
        DO 420 J=1,NTS
          DELTAJK=0.0D+00
          IF (K.EQ.J) DELTAJK=AS(J)
          DMATP=DMATP
     *          -DMATM1(I,J)*(DELTAJK/TWO-DI(J,K))
 420  CONTINUE
C       DMATM1(I,K)=-DMATM1(I,K)
C        SI(I,K)=DMATP
         DDIMT2(K)=DMATP
 410  CONTINUE
C          SUM UP
        IF(GOPARR) THEN
          CALL DDI_GSUMF(2406,DDIMT2,NTS)
        END IF
           DO 450 J=1,NTS
           SI(I,J)=DDIMT2(J)
 450  CONTINUE
 400  CONTINUE
C
      DO I=1,NTS
      DO K=1,NTS
        DMATM1(I,K)= SI(I,K)
      ENDDO
      ENDDO
      IF(IEF.EQ.5.OR.IEF.EQ.8) THEN
C
C     THE ORIENTATIONAL APPARENT CHARGES ARE READ FROM DISK FILE -NFT29-
C     ASCOR - NFT29 - ORIENTATIONAL APPARENT CHARGES.
C
      NFT29 = 29
      CALL SEQOPN(NFT29,'ASCOR','OLD',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT29)
C
      CALL SQREAD(NFT29,QOR,NTS)
      LEN=1
      CALL SQREAD(NFT29,PORT,LEN)
      CALL SEQREW(NFT29)
C
C         POTENTIAL AND FIELD DUE TO THE ADDITIONAL ORIENTATIONAL
C         SURFACE CHARGE
C
       DO I=1,NTS
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
        SSI=ZERO
        DO J=1,NTS
         IF (J.EQ.I) THEN
          SS0=1.07D+00*SQRT(FPI/AS(I))
          SSI=SSI+QOR(I)*SS0
         ELSE
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          SSI=SSI+QOR(J)/DRIJ
         END IF
        ENDDO
        VAD(I)=SSI
       ENDDO
      END IF
C
      RETURN
      END
C*MODULE PCMIEF  *DECK IEFCMV1
      SUBROUTINE IEFCMV1(DMATM1,WORK,IPVT,DDIMT1,DDIMT2,XCTS,YCTS,
     *                   ZCTS,AS,VAD,QOR,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMATM1(NTS,NTS),
     *          WORK(NTS),IPVT(NTS),DDIMT1(NTS),DDIMT2(NTS),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),VAD(*),QOR(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
C
      DIMENSION DET(2)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
      DATA ZERO,TWO/0.0D+00,2.0D+00/
      DATA FPI/12.56637061D+00/
C
C        COMPUTE THE C AND THEN C-INVERSE MATRIX
C
         IF(MASWRK) WRITE(IW,*) ' -------------------'
         IF(MASWRK) WRITE(IW,*) ' -- MATRIX C^(-1) --'
         IF(MASWRK) WRITE(IW,*) ' -------------------'
C
C   COMPUTE THE C MATRIX, STORED AT -DMATM1-
C   ISOTROPIC DIELECTRICS
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 1: COMPUTING MATRIX A '
C
      FACT=(EPS+1.0D+00)/(EPS-1.0D+00)
      DO 200 I=1,NTS
C
C           INITIALIZE PARALLEL
      IPCOUNT = ME - 1
      CALL VCLR(DDIMT1,1,NTS)
C
      DO 210 K=1,NTS
        DMATP=0.0D+00
C            GO PARALLEL!
        IF(GOPARR) THEN
          IPCOUNT = IPCOUNT + 1
          IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 210
        END IF
C
        DO 220 J=1,NTS
C
C       -- COMPUTE DDIJ --
C
        IF(J.EQ.I) THEN
          LI=ISPHE(I)
          SSI=AS(I)*1.07D+00*SQRT(FPI*AS(I))/FPI
          DDIJ=-SSI/(2.0D+00*RE(LI))
        ELSE
          XI=XCTS(I)
          YI=YCTS(I)
          ZI=ZCTS(I)
          LJ=ISPHE(J)
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          CXJ=(XJ-XE(LJ))/RE(LJ)
          CYJ=(YJ-YE(LJ))/RE(LJ)
          CZJ=(ZJ-ZE(LJ))/RE(LJ)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          DRIJ3=DRIJ**3
C
          DDIJ=AS(I)*AS(J)
     *        *((XI-XJ)*CXJ+(YI-YJ)*CYJ+(ZI-ZJ)*CZJ)/(FPI*DRIJ3)
        END IF
C
C       -- COMPUTE SSJK --
C
        IF (K.EQ.J) THEN
          SSJK=AS(J)*1.07D+00*SQRT(FPI*AS(J))/FPI
        ELSE
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          XK=XCTS(K)
          YK=YCTS(K)
          ZK=ZCTS(K)
          DRJK=SQRT((XJ-XK)**2+(YJ-YK)**2+(ZJ-ZK)**2)
C
          SSJK=AS(J)*AS(K)/(FPI*DRJK)
C
        END IF
C
          DELTAIJ=0.0D+00
          IF (I.EQ.J) DELTAIJ=AS(I)*FACT
          DMATP=DMATP
     *          +(DELTAIJ/TWO-DDIJ)*(1.0D+00/AS(J))*SSJK
 220  CONTINUE
C
         DDIMT1(K)=FPI*DMATP
 210  CONTINUE
C          SUM UP
        IF(GOPARR) THEN
          CALL DDI_GSUMF(2405,DDIMT1,NTS)
        END IF
           DO 250 J=1,NTS
           DMATM1(I,J)=DDIMT1(J)
 250  CONTINUE
C
 200  CONTINUE
C
      IF(IPRINT.EQ.1 .AND. MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 2: COMPUTING MATRIX A^(-1) '
C
      INFO=0
      CALL DGEFA(DMATM1,NTS,NTS,IPVT,INFO)
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'THE A MATRIX IS SINGULAR'
         CALL ABRT
      END IF
C
      CALL DGEDI(DMATM1,NTS,NTS,IPVT,DET,WORK,01)
C
      IF(IPRINT.EQ.1 .AND. MASWRK) THEN
         WRITE(IW,*) ' IDMM: THE A INVERSE MATRIX IS'
         CALL PRSQ(DMATM1,NTS,NTS,NTS)
      END IF
C
      IF(MASWRK) WRITE(IW,*) ' '
      IF(MASWRK) CALL TIMIT(1)
      IF(MASWRK) WRITE(IW,*) 'STEP 3: COMPUTING MATRIX C^(-1)'
      IF(MASWRK) WRITE(IW,*) ' '
C
      DO 400 I=1,NTS
C
C           INITIALIZE PARALLEL
      IPCOUNT = ME - 1
      CALL VCLR(DDIMT2,1,NTS)
C
      DO 410 K=1,NTS
        DMATP=0.0D+00
C            GO PARALLEL!
        IF(GOPARR) THEN
          IPCOUNT = IPCOUNT + 1
          IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 410
        END IF
C
        DO 420 J=1,NTS
C
C      -- COMPUTE DDJK --
C
        IF(J.EQ.K) THEN
          LJ=ISPHE(J)
          SSJ=AS(J)*1.07D+00*SQRT(FPI*AS(J))/FPI
          DDJK=-SSJ/(2.0D+00*RE(LJ))
        ELSE
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          LK=ISPHE(K)
          XK=XCTS(K)
          YK=YCTS(K)
          ZK=ZCTS(K)
          CXK=(XK-XE(LK))/RE(LK)
          CYK=(YK-YE(LK))/RE(LK)
          CZK=(ZK-ZE(LK))/RE(LK)
          DRJK=SQRT((XJ-XK)**2+(YJ-YK)**2+(ZJ-ZK)**2)
          DRJK3=DRJK**3
C
          DDJK=AS(J)*AS(K)
     *        *((XJ-XK)*CXK+(YJ-YK)*CYK+(ZJ-ZK)*CZK)/(FPI*DRJK3)
        END IF
C
          DELTAJK=0.0D+00
          IF (K.EQ.J) DELTAJK=AS(J)
          DMATP=DMATP
     *          -DMATM1(I,J)*(DELTAJK/TWO-DDJK)
 420  CONTINUE
C
         DDIMT2(K)=DMATP
 410  CONTINUE
C          SUM UP
        IF(GOPARR) THEN
          CALL DDI_GSUMF(2406,DDIMT2,NTS)
        END IF
           DO 450 J=1,NTS
           DMATM1(I,J)=DDIMT2(J)
 450  CONTINUE
 400  CONTINUE
C
      IF(IEF.EQ.5.OR.IEF.EQ.8) THEN
C
C     THE ORIENTATIONAL APPARENT CHARGES ARE READ FROM DISK FILE -NFT29-
C     ASCOR - NFT29 - ORIENTATIONAL APPARENT CHARGES.
C
      NFT29 = 29
      CALL SEQOPN(NFT29,'ASCOR','OLD',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT29)
C
      CALL SQREAD(NFT29,QOR,NTS)
      LEN=1
      CALL SQREAD(NFT29,PORT,LEN)
      CALL SEQREW(NFT29)
C
C         POTENTIAL AND FIELD DUE TO THE ADDITIONAL ORIENTATIONAL
C         SURFACE CHARGE
C
       DO I=1,NTS
        XI=XCTS(I)
        YI=YCTS(I)
        ZI=ZCTS(I)
        SSI=ZERO
        DO J=1,NTS
         IF (J.EQ.I) THEN
          SS0=1.07D+00*SQRT(FPI/AS(I))
          SSI=SSI+QOR(I)*SS0
         ELSE
          XJ=XCTS(J)
          YJ=YCTS(J)
          ZJ=ZCTS(J)
          DRIJ=SQRT((XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2)
          SSI=SSI+QOR(J)/DRIJ
         END IF
        ENDDO
        VAD(I)=SSI
       ENDDO
      END IF
C
      RETURN
      END
C*MODULE PCMIEF  *DECK ICVEV
      SUBROUTINE ICVEV(DMATM1,BEMPOT,VPOT,VCAM,L2,NFT27,SOME,XCTS,
     *                 YCTS,ZCTS,AS,QSN,VAD,QOR,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DMATM1(NTS,NTS),BEMPOT(L2),VPOT(NTS),VCAM(NTS),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),VAD(*),QOR(*),
     *          XE(*),YE(*),ZE(*),RE(*),ISPHE(*)
C
      LOGICAL SOME
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      PARAMETER (MXATM=2000)
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 /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,N26,N27,IREP,IEF,IP_F,NFMOPCM
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
C
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      DATA ZERO, ONE/0.0D+00, 1.0D+00/
      DATA FPI/12.56637061D+00/
C
C-- FOR RPAC  7/11/97 ----------------------------------------------
C       PCMINTS - NFT28 - ELEC POT INTS AT THE  TESSERAE'S CENTRE.
C------------------------------------------------------------------
C
C-    IF(RPAC) THEN
C-     NFT28=28
C-     CALL SEQREW(NFT28)
C-    END IF
C
C  2) CALCULATION OF APPARENT CHARGES GENERATED BY THE SOLUTE'S NUCLEI.
C
      CALL SEQREW(NFT27)
      DO  ITS = 1, NTS
        XI=XCTS(ITS)
        YI=YCTS(ITS)
        ZI=ZCTS(ITS)
         CALL INTMEP(BEMPOT,XI,YI,ZI)
C-       IF(RPAC) CALL SQWRIT(NFT28,BEMPOT,L2)
         CALL SQWRIT(NFT27,BEMPOT,L2)
        L=ISPHE(ITS)
        CX=(XI-XE(L))/RE(L)
        CY=(YI-YE(L))/RE(L)
        CZ=(ZI-ZE(L))/RE(L)
        VCAM(ITS) = ZERO
        VPOT(ITS) = ZERO
         DO  JATOM = 1, NAT
          XL=C(1,JATOM)
          YL=C(2,JATOM)
          ZL=C(3,JATOM)
          RIL=SQRT((XI-XL)**2+(YI-YL)**2+(ZI-ZL)**2)
          XM=((XI-XL)*CX+(YI-YL)*CY+(ZI-ZL)*CZ)/RIL**3
          VPOT(ITS) = VPOT(ITS) - ZAN(JATOM)/RIL
          VCAM(ITS) = VCAM(ITS) + XM*ZAN(JATOM)
         ENDDO
      ENDDO
C-    IF(RPAC) CALL SEQREW(NFT28)
      CALL SEQREW(NFT27)
C
      QNTOT = ZERO
      QNUC2 = ZERO
      DO ITS = 1, NTS
         QSN(ITS) = ZERO
         DO JTS = 1, NTS
          VP=VPOT(JTS)
          IF(IEF.EQ.5.OR.IEF.EQ.8)VP = VP - VAD(JTS)
          QSN(ITS) = QSN(ITS)-DMATM1(ITS,JTS)*VP
         ENDDO
         QSN(ITS) = QSN(ITS)*AS(ITS)
         QNUC2 = QNUC2 + VCAM(ITS)*AS(ITS)
         QNTOT = QNTOT + QSN(ITS)
      ENDDO
      QNUC2 = QNUC2/FPI
C
C     SCALE C-PCM CHARGE BY (EPS-1)/EPS
C
      IF(IEF.EQ.10) THEN
        SCALE  = (EPS-ONE)/EPS
        DO ITS = 1, NTS
           QSN(ITS) = QSN(ITS)*SCALE
        END DO
        QNUC2 = QNUC2*SCALE
        QNTOT = QNTOT*SCALE
      END IF
C
      IF(SOME .AND. MASWRK) WRITE(IW,1200) QNUC2
C
C  4) RENORMALIZATION OF `NUCLEAR' APPARENT CHARGES.
C     IF ICOMP=1 THE CORRECTION IS PROPORTIONAL TO THE TESSERA'S AREA.
C     IF ICOMP=2 OR ICOMP=3 THE CORRECTION IS EQUAL FOR EACH TESSERA.
C
      CHG = ZERO
      DO JATOM = 1, NAT
         CHG = CHG + ZAN(JATOM)
      ENDDO
C
      IF(IEF.EQ.5.OR.IEF.EQ.8) THEN
       QTOR=ZERO
       DO ITS=1,NTS
        QTOR=QTOR+QOR(ITS)
       ENDDO
       CHG=CHG+QTOR
      END IF
C
      TCH = - CHG * (EPS - ONE) / EPS
      QNTOTN = ZERO
C
C     OPTION 1:
C
      FN = ONE
      IF(ICOMP.EQ.1) THEN
         SUPTOT = STOT * ANTOAU * ANTOAU
         DIFF = TCH - QNTOT
         DO ITS = 1, NTS
            QSN(ITS) = QSN(ITS) + DIFF * AS(ITS) / SUPTOT
            QNTOTN = QNTOTN + QSN(ITS)
         ENDDO
C
C     OPTION 2:
C
      ELSE IF(ICOMP.EQ.2) THEN
         FN = TCH / QNTOT
         DO ITS = 1, NTS
            QSN(ITS) = QSN(ITS) * FN
            QNTOTN = QNTOTN + QSN(ITS)
         ENDDO
      END IF
C
C     TOTAL APPARENT CHARGE BEFORE AND AFTER RENORMALIZATION.
C
      IF(SOME .AND. MASWRK) THEN
         IF(ICOMP.NE.0) WRITE(IW,1000) QNTOT, TCH, QNTOTN
         IF(ICOMP.EQ.0) WRITE(IW,1100) QNTOT, TCH
      END IF
      RETURN
C
 1000 FORMAT(1X,'NUCLEAR APPARENT CHARGE',F10.5,' THEORETICAL',
     *          F10.5,' RENORMALIZED',F10.5)
 1100 FORMAT(1X,'NUCLEAR APPARENT CHARGE',F10.5,' THEORETICAL',
     *          F10.5,' NOT RENORMALIZED',F10.5)
 1200 FORMAT(/1X,'ESTIMATE OF NUCLEAR CHARGE',F15.5)
      END
C*MODULE PCMIEF  *DECK IXMATV
      SUBROUTINE IXMATV(D,XX,SCR,DMATM1,QPOT,SCRPOT,
     *                  Q0,Q1,Q2,Q3,D0,QA,DIMAT,QREP,TMP,
     *                  TMP1,TMP2,IPVT,POTTMP,RMUL,QET,TCH,QESC,QETN,
     *                  NFT27,IPCFP,L2,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,
     *                  AS,QSN,QSE,Q_FS,VEC_MUL,
     *                  CHG2,CHG2N,PEL,XE,YE,ZE,RE,ISPHE,LIST,
     *                  ASCCRD,ASCCHG,ASCDIP,ASCQAD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL NOTRF
C
      DIMENSION D(L2),XX(L2),SCR(L2),DMATM1(NTS,NTS),QPOT(NTS),
     *          SCRPOT(L2),Q0(NTS),Q1(NTS),Q2(NTS),Q3(NTS),
     *          D0(NTS),QA(NTS),DIMAT(MXDII1,MXDII1),
     *          QREP(NTSPAR,MXDIIS,2),TMP(NTS,3),TMP1(MXDII1),
     *          TMP2(MXDII1,MXDII1),IPVT(MXDII1),POTTMP(NTS),
     *          RMUL(MXSP,10),XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),
     *          QSE(*),Q_FS(*),VEC_MUL(*),
     *          CHG2(*),CHG2N(*),PEL(*),XE(*),YE(*),ZE(*),
     *          RE(*),ISPHE(*),LIST(*),ASCCRD(3,*),ASCCHG(*),
     *          ASCDIP(3,*),ASCQAD(6,*)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR, DSKWRK, MASWRK, USEIMG, IPCFP
C
      COMMON /CONV  / ACURCY,ENUCR,AJUNK(3),DIFF,ITERRO,ICALCP,ICBET
      COMMON /FMCOM / WORK(1)
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /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 /NTNOPT/ ENERGY,ENERG0,DEMAX,SQCDF,ITERMC,MICRO,NOTRF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMCHG/ PB,PX,PC,UNZ,QNUC,FN,FE,UNZOR,PORT
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMEFP/ PB_MUL,PC_MUL,PB_POL,PC_POL,P_FF,P_NF,P_FN,P_I,
     *                P_IBIS,P_J,P_NUCC,P_NUCCBIS,ENPCM
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,MFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      DATA ZERO, ONE/0.0D+00, 1.0D+00/
      DATA RMC/8HMCSCF   /
      DATA GVB/8HGVB     /
C
C     ROUTINE FOR THE CALCULATION OF MATRIX XX (INTERACTION BETWEEN
C     `ELECTRONIC' APPARENT CHARGES AND ELECTRONIC POTENTIAL OF THE
C     SOLUTE):
C         XX(M,N) = SUM_I[ V(M,N;I)QSE(I) ]
C     WHICH IS USED IN THE SCF ROUTINE IN ORDER TO CORRECT TWO-ELECTRON
C     INTEGRALS.
C
C
      IF(SCFTYP.EQ.RMC) THEN
        ITER=ITERMC-1
      ELSE IF(SCFTYP.EQ.GVB) THEN
        ITER=ITERRO+1
      ELSE
        ITER=ITERRO
      END IF
C
C     --- GENERATE AN IMAGE OF THE AB INITIO REGION
C
      IF(IPCMIT.EQ.1 .AND. IMGABI.EQ.1) CALL IMABIX
C
C     --- COMPUTE THE ELECTRONIC POTENTIAL AT EACH TESSERA
C
C     IF(MASWRK.AND.ITER.EQ.1) WRITE(IW,*) ' '
C     IF(MASWRK.AND.ITER.EQ.1) CALL TIMIT(1)
C     IF(MASWRK.AND.ITER.EQ.1) WRITE(IW,*) ' '
C     IF(MASWRK.AND.ITER.EQ.1) WRITE(IW,*)
C    *   'CALLING QPOT(ITS) SQREAD(NFT27) ...'
C
      IF(IPCMIT.EQ.1 .AND. IDIRCT.EQ.1) THEN
        IF(IMGABI.EQ.1) THEN
C         INITIALIZE PARALLEL
          IPCOUNT = ME - 1
          CALL VCLR(QPOT,1,NTS)
          KTIMES = 0
C
          DO 90 ITS = 1, NTS
C           GO PARALLEL!
            IF(GOPARR) THEN
              IPCOUNT = IPCOUNT + 1
              IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 90
            END IF
C
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
C           INTRODUCE IMGABI
            USEIMG=XI.GE.XMAX+RABI*ANTOAU
     *      .OR.   XI.LE.XMIN-RABI*ANTOAU
     *      .OR.   YI.GE.YMAX+RABI*ANTOAU
     *      .OR.   YI.LE.YMIN-RABI*ANTOAU
     *      .OR.   ZI.GE.ZMAX+RABI*ANTOAU
     *      .OR.   ZI.LE.ZMIN-RABI*ANTOAU
            IF(USEIMG) THEN
              CALL POTIMG(QADD,XI,YI,ZI)
              QPOT(ITS) = - QADD
              KTIMES = KTIMES + 1
            ELSE
              CALL INTMEP(SCRPOT,XI,YI,ZI)
              QPOT(ITS) = TRACEP(D,SCRPOT,NUM)
            END IF
  90      CONTINUE
C               SUM UP
          IF(GOPARR) THEN
            CALL DDI_GSUMF(2451,QPOT,NTS)
            CALL DDI_GSUMI(2460,KTIMES,1)
          END IF
C         IF(MASWRK.AND.ITER.EQ.1 .AND. IMGABI.EQ.1)
C    *             WRITE(IW,*) 'IMGABI USED',KTIMES,'TIMES'
C
        ELSE
          CALL VCLR(QPOT,1,NTS)
          CALL EPOTEN(ONE,XCTS,YCTS,ZCTS,QPOT,D,NTS,L2)
C         EPOTEN IS PARALLELISED, BUT NO GLOBAL SUM INSIDE.
          IF(GOPARR) CALL DDI_GSUMF(2451,QPOT,NTS)
        ENDIF
      ELSE
C
        CALL SEQREW(NFT27)
        DO ITS = 1, NTS
          CALL SQREAD(NFT27,SCRPOT,L2)
          QPOT(ITS) = TRACEP(D,SCRPOT,NUM)
        ENDDO
        CALL SEQREW(NFT27)
C
      END IF
C
C     IF(MASWRK.AND.ITER.EQ.1) CALL TIMIT(1)
C     IF(MASWRK.AND.ITER.EQ.1) WRITE(IW,*)
C    *    'QPOT(ITS) SQREAD(NFT27) IS DONE'
C     IF(MASWRK.AND.ITER.EQ.1) WRITE(IW,*) ' '
C
C     --- STORE THE POTENTIAL
C
      CALL DCOPY(NTS,QPOT,1,PEL,1)
      CALL DSCAL(NTS,-ONE,PEL,1)
C
C
C    --- COMPUTE THE PCM SURFACE CHARGE
C
C     BRANCH BETWEEN MATRIX INVERSION AND ITERATIVE CALCULATIONS
      ENPCM = ZERO
C
C        -- MATRIX INVERSION --
      IF (IPCMIT.EQ.0) THEN
C
        QET = ZERO
        DO ITS = 1, NTS
          QSE(ITS) = ZERO
          DO ITSJ = 1, NTS
            QSE(ITS) = QSE(ITS)-DMATM1(ITS,ITSJ)*QPOT(ITSJ)
          ENDDO
          QSE(ITS) = QSE(ITS)*AS(ITS)
          QET = QET + QSE(ITS)
        ENDDO
C
C       SCALE C-PCM CHARGE BY (EPS-1)/EPS
C
        IF(IEF.EQ.10) THEN
          SCALE  = (EPS-ONE)/EPS
          DO ITS = 1, NTS
             QSE(ITS) = QSE(ITS)*SCALE
          END DO
          QET = QET*SCALE
        END IF
C
      ELSE
C
C       -- ITERATIVE --
C
C           -- COMPUTE THE TOTAL POTENTIAL AT EACH TESSERA
C
C  USE SEPARATE CHARGE FOR ICOMP=2
C
        IF(ICOMP.EQ.2 .AND.(.NOT.IPCFP)) THEN
C          DO ITS = 1, NTS
C            QSN(ITS)=POTTMP(ITS)
C          END DO
           IF(ITER.EQ.1) THEN
               CALL VCLR(QSN,1,NTS)
               DO ITS=1,NTS
                  DO IAT=1,NAT
                     DINVR=1.0D+00/SQRT((XCTS(ITS)-C(1,IAT))**2+
     *                                  (YCTS(ITS)-C(2,IAT))**2+
     *                                  (ZCTS(ITS)-C(3,IAT))**2)
                     QSN(ITS)=QSN(ITS)-ZAN(IAT)*DINVR
                  ENDDO
               ENDDO
C       NUCLEAR CHARGES ARE GOT ONLY ONCE FOR THE CHARGE COMPENSATION.
C       TRY TO GIVE MORE ITERATIONS TO ENSURE CONVERGENCE.
               MXISAV1=MXITR1
               MXISAV2=MXITR2
               MXITR1=MAX(MXITR1,2000)
               MXITR2=MAX(MXITR2,2000)
               CALL ASCIT(ITER,QSN,Q0,Q1,Q2,Q3,D0,QA,DIMAT,QREP,TMP,TMP1
     *                   ,TMP2,IPVT,RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,
     *                    AS,QSE,XE,YE,ZE,RE,ISPHE,LIST)
               MXITR1=MXISAV1
               MXITR2=MXISAV2
               DO ITS = 1, NTS
                 CHG2N(ITS) = Q1(ITS)
                 SCALE = ONE
                 IF(IEF.EQ.10) SCALE=(EPS-ONE)/EPS
                 QSN(ITS)=QSE(ITS)*SCALE
               END DO
           ELSE
               DO ITS = 1, NTS
                 QSN(ITS)=POTTMP(ITS)
               END DO
           END IF
        END IF
C
        CALL ADDPOT(QPOT,VEC_MUL,IPCFP,XCTS,YCTS,ZCTS)
C
C           -- COMPUTE THE SURFACE CHARGE ON EACH TESSERA
C
        THRSBK=THRES
        IF(SCFTYP.NE.RMC.AND.SCFTYP.NE.GVB) THEN
        IF(THRSLS.GT.THRES .AND. (DIFF.GT.DENSLS.OR.ITER.EQ.1))
     *         THRES=THRSLS
        IF(THRES.EQ.THRSBK .AND. DENSLS.LT.99.0D+00) THEN
            IF(MASWRK) WRITE(IW,*)
     *  '          * * *   PCM SWITCHES TO FULL ACCURACY * * *'
            DENSLS=100.0D+00
        END IF
        END IF
        CALL ASCIT(ITER,QPOT,Q0,Q1,Q2,Q3,D0,QA,DIMAT,QREP,TMP,TMP1,TMP2,
     *             IPVT,RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,QSE,XE,YE,
     *             ZE,RE,ISPHE,LIST)
        THRES=THRSBK
C
C  SCALE C-PCM CHARGE BY (EPS-1)/EPS
C
        IF(IEF.EQ.10) THEN
          DO ITS = 1, NTS
             QSE(ITS) = QSE(ITS)*(EPS-ONE)/EPS
          END DO
        END IF
C
C  USE SEPARATE CHARGE FOR ICOMP=2
C
        QTOTUN=ZERO
        IF(ICOMP.EQ.2 .AND.(.NOT.IPCFP)) THEN
          QNT =0.0D+00
          QET =0.0D+00
          DO ITS = 1, NTS
            QNT =QNT +QSN(ITS)
            QET =QET + (QSE(ITS)-QSN(ITS))
          END DO
          FN = (NE + ICH)*(EPS-ONE)/EPS/ABS(QNT)
          FE = NE*(EPS-ONE)/EPS/ABS(QET)
C         IF(MASWRK) WRITE(IW,*) 'FN=',FN,'FE=',FE
          DO ITS = 1, NTS
            CHG2(ITS)=CHG2N(ITS)*FN+
     *                (Q1(ITS)-CHG2N(ITS))*FE
C           V_ELE(ITS)=QSE(ITS)-QSN(ITS)
            QSE(ITS)=QSN(ITS)*FN +(QSE(ITS)-QSN(ITS))*FE
            POTTMP(ITS)=QSN(ITS)
            QSN(ITS)=0.0D+00
          ENDDO
C
          QTOTUN=QNT+QET
        END IF
C
C
C           -- COMPUTE TOTAL PCM INTERACTION
C              ALSO COMPUTE THE AVERAGE ABSOLUTE CHARGE
C
        QET = ZERO
        AVEASC=0.0D+00
        DO ITS=1,NTS
          QET=QET+QSE(ITS)
          AVEASC=AVEASC+ABS(QSE(ITS))
          VNUC = ZERO
          DO IAT=1,NAT
             DINVR=1.0D+00/SQRT(
     *            (XCTS(ITS)-C(1,IAT))**2+
     *            (YCTS(ITS)-C(2,IAT))**2+
     *            (ZCTS(ITS)-C(3,IAT))**2)
             VNUC=VNUC+ZAN(IAT)*DINVR
          ENDDO
          ENPCM = ENPCM + (PEL(ITS)+VNUC+VEC_MUL(ITS))*QSE(ITS)
        ENDDO
        IF(ICOMP.NE.2.OR.IPCFP) QTOTUN=QET
C
C       AVEASC=AVEASC/REAL(NTS)
C       IF(MASWRK) WRITE(IW,*) 'AVEASC=',AVEASC
C       IF(MASWRK) WRITE(IW,*) 'QET=',QET
C
      END IF
C
C  SKIP NORMALIZATION IF IPCMIT=1
      IF(IPCMIT.EQ.1) GO TO 100
C
C
C
C     --- NORMALIZATION ---
C
C  4) NORMALIZZAZIONE DELLE CARICHE VIRTUALI ELEMENTARI ELETTRONICHE.
C     PER ICOMP=1 LA CORREZIONE E' DISTRIBUITA IN PROPORZIONE
C     ALL'AREA DELLA TESSERA.
C     PER ICOMP=2 LA CORREZIONE E' EFFETTUATA CON FATTORI COSTANTI.
C
      TCH = NE * (EPS - ONE) / EPS
C
C     OPTION 1
C
      QETN = ZERO
      FE = ONE
      IF(ICOMP.EQ.1) THEN
         SUPTOT = STOT * ANTOAU * ANTOAU
         DIFF = TCH - QET
         DO ITS = 1, NTS
            QSE(ITS) = QSE(ITS) + DIFF * AS(ITS) / SUPTOT
            QETN = QETN + QSE(ITS)
        ENDDO
C
C     OPTION 2
C
      ELSE IF(ICOMP.EQ.2) THEN
         FE = TCH / QET
         DO ITS = 1, NTS
            QSE(ITS) = QSE(ITS) * FE
            QETN = QETN + QSE(ITS)
         ENDDO
      END IF
C
C     PRINT TOTAL INDUCED CHARGE BEFORE AND AFTER RENORMALIZATION
C
      IF(NPRINT.EQ.5 .AND. MASWRK) THEN
         IF(ICOMP.NE.0) WRITE(IW,1000) QET, TCH, QETN
         IF(ICOMP.EQ.0) WRITE(IW,1100) QET, TCH
      END IF
C
 100  CONTINUE
C
C     CALCULATION OF INTERACTION MATRICES TO BE ADDED TO THE FOCK MATRIX
C
C     PB = INTERACTION ELECTRONS-NUCLEAR INDUCED CHARGES
C     PC = INTERACTION NUCLEI-ELECTRONIC INDUCED CHARGES
C     PX = INTERACTION ELECTRONS-ELECTRONIC INDUCED CHARGES
C     PB_MUL = INTERACTION ELECTRONS-EFP MULTIPOLE INDUCED CHRG
C     PC_MUL = INTERACTION EFP MUL POTENTIAL-ELECTRON INDUCED CHRG
C     PB_POL = INTERACTION ELECTRONS-EFP POLARIZABLE INDUCED CHRG
C     PC_POL = INTERACTION EFP POL POTENTIAL-ELECTRON INDUCED CHRG
C
      CALL VCLR(XX,1,L2)
      PB = ZERO
      PB_MUL=ZERO
      PB_POL=ZERO
      PC = ZERO
      PC_MUL=ZERO
      PC_POL=ZERO
      PX = ZERO
C
C         PRINT CPU TIME
C     IF(MASWRK.AND.ITER.EQ.1) THEN
C     WRITE(IW,*) ' '
C     CALL TIMIT(1)
C     WRITE(IW,*) ' '
C     WRITE(IW,*) 'CALLING XX(IBAS) SQREAD(NFT27) ...'
C     END IF
C
C
C         -- MATRIX INVERSION
C
      IF(IPCMIT.EQ.0) THEN
        CALL SEQREW(NFT27)
        DO 110 ITS = 1, NTS
          CALL SQREAD(NFT27,SCR,L2)
          DO IBAS = 1, L2
            XX(IBAS) = XX(IBAS) - SCR(IBAS) * QSE(ITS)
          ENDDO
C
          PB = PB + PEL(ITS) * QSN(ITS)
          PX = PX + PEL(ITS) * QSE(ITS)
          PB_MUL=PB_MUL+PEL(ITS) * Q_FS(ITS)
          PC_MUL=PC_MUL+VEC_MUL(ITS)*QSE(ITS)
C
C         NOT HERE, MUST BE IN DECK CHGIND
C         PB_POL=PB_POL+PEL(ITS) * Q_IND(ITS)
C
C             PRUNED AWAY, NOVEMBER 2005 BY HUI
C         PC_POL=PC_POL+(AIND_PROJ(ITS)+
C    *           FIND_PROJ(ITS))*QSE(ITS)
C
          DO JATOM = 1, NAT
            R2 = (C(1,JATOM)-XCTS(ITS))**2 + (C(2,JATOM)-
     *           YCTS(ITS))**2 + (C(3,JATOM)-ZCTS(ITS))**2
            R  = SQRT(R2)
            PC = PC + QSE(ITS) * ZAN(JATOM) / R
          ENDDO
C
 110    CONTINUE
        CALL SEQREW(NFT27)
      END IF
C
C
C         -- ITERATIVE
C
      IF(IPCMIT.EQ.1) THEN
C
C             -- NO USE IMGASC --
C
        IF(IMGASC.EQ.0) THEN
C
          IF(IDIRCT.EQ.1) THEN
C                   INITIALIZE PARALLEL
            IPCOUNT = ME - 1
C
            DO 120 ITS = 1, NTS
C                     GO PARALLEL!
              IF(GOPARR) THEN
                IPCOUNT = IPCOUNT + 1
                IF(MOD(IPCOUNT,NPROC).NE.0) GO TO 120
              END IF
C
              CALL INTMEP(SCR,XCTS(ITS),YCTS(ITS),ZCTS(ITS))
              DO IBAS = 1, L2
                XX(IBAS) = XX(IBAS) - SCR(IBAS) * QSE(ITS)
              ENDDO
 120        CONTINUE
C                   SUM UP
            IF(GOPARR) THEN
              CALL DDI_GSUMF(2452,XX,L2)
            END IF
          ELSE
C
            CALL SEQREW(NFT27)
            DO 130 ITS = 1, NTS
              CALL SQREAD(NFT27,SCR,L2)
              DO IBAS = 1, L2
                XX(IBAS) = XX(IBAS) - SCR(IBAS) * QSE(ITS)
              ENDDO
 130        CONTINUE
            CALL SEQREW(NFT27)
          END IF
C
C
C              -- USE IMGASC --
        ELSE
C
C         GENERATE PCM MULTIPOLE IMAGE
C
          CALL VALFM(LOADFM)
          LTMP   = LOADFM + 1
          NEED   = 3*NTS
          CALL GETFM(NEED)
          CALL IMASCX(WORK(LTMP),XCTS,YCTS,ZCTS,QSE,XE,YE,ZE,ISPHE,LIST,
     *                ASCCRD,ASCCHG,ASCDIP,ASCQAD)
          CALL RETFM(NEED)
C
          CALL CHGINT(XX,XCTS,YCTS,ZCTS,QSE,LIST,ASCCRD,ASCCHG)
C         IF(MASWRK) WRITE(IW,*) ' XX(1)=',XX(1),' XX(2)=',XX(2)
C         IF(MASWRK) WRITE(IW,*) 'CALLED CHGINT'
          CALL DCOPY(L2,XX,1,SCR,1)
C
          CALL DPLINT(XX,ASCCRD,ASCDIP)
C         IF(MASWRK) WRITE(IW,*) ' XX(1)=',XX(1),' XX(2)=',XX(2)
C         IF(MASWRK) WRITE(IW,*) 'CALLED DPLINT'
          CALL VADD(SCR,1,XX,1,SCR,1,L2)
C
          CALL QADINT(XX,ASCCRD,ASCQAD)
C         IF(MASWRK) WRITE(IW,*) ' XX(1)=',XX(1),' XX(2)=',XX(2)
C         IF(MASWRK) WRITE(IW,*) 'CALLED QADINT'
          CALL VADD(SCR,1,XX,1,XX,1,L2)
C
        END IF
C
          DO 140 ITS = 1, NTS
            PX = PX + PEL(ITS) * QSE(ITS)
 140      CONTINUE
C
      END IF
C
C     IF(MASWRK) WRITE(IW,*) ' XX(1)=',XX(1),' XX(2)=',XX(2)
C     IF(MASWRK) WRITE(IW,*) ' XX(3)=',XX(3),' XX(4)=',XX(4)
C
C     IF(MASWRK.AND.ITER.EQ.1) THEN
C     CALL TIMIT(1)
C     WRITE(IW,*) 'XX(IBAS) SQREAD(NFT27) IS DONE'
C     WRITE(IW,*) ' '
C     END IF
C
      PB=PB+PB_MUL+PB_POL
      PC=PC+PC_MUL+PC_POL
C     NOTE PX HAS DIFFERENT MEANING IF IPCMIT
      IF(IPCMIT.NE.0) THEN
        PX=ENPCM - PX
        P_FF=ZERO
        P_NF=ZERO
        P_FN=ZERO
        P_I=ZERO
        P_IBIS=ZERO
        P_J=ZERO
        P_NUCC=ZERO
        P_NUCCBIS=ZERO
        PB_MUL=ZERO
        PB_POL=ZERO
        PC_MUL=ZERO
        PC_POL=ZERO
        PB=ZERO
        PC=ZERO
      END IF
C     QET=QTOTUN
C     QESC=QTOTUN-QET
      QESC=QTOTUN
C     THE MEANING OF QESC IS DIFFERENT COMPARED TO ICOMP=3.  HERE
C     IT HOLDS THE UNNORMALISED TOTAL CHARGE FOR PRINTING PURPOSES.
C
      RETURN
C
 1000 FORMAT(1X,'ELECTRONIC APPARENT CHARGE',F10.5,' THEORETICAL ',
     * F10.5,' RENORMALIZED',F10.5)
 1100 FORMAT(1X,'ELECTRONIC APPARENT CHARGE',F10.5,' THEORETICAL ',
     * F10.5,' NOT RENORMALIZED',F10.5)
      END
C*MODULE PCMIEF  *DECK DERIEFX
      SUBROUTINE DERIEFX(XCTS,YCTS,ZCTS,AS,QSN,QSE,Q_FS,Q_IND,XE,YE,ZE,
     *                   RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
      DIMENSION XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),QSE(*),
     *          Q_FS(*),Q_IND(*),XE(*),YE(*),ZE(*),RE(*),ISPHE(*)
C
      COMMON /FMCOM / XX(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 /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
      COMMON /PCMPLY/ STOT,VOL,CCX,CCY,CCZ,RDIF,NESF,NESFP,NESFF,I_NESF,
     *                L_AST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
C     ----- POLARIZABLE CONTINUUM MODEL GRADIENT CORRECTIONS -----
C
C        REGENERATE ORIGINAL CAVITY IF NECESSARY
C
C     IF(ICAV.EQ.1  .OR.  IDISP.EQ.1) CALL PEDRAM
C
      L2 = (NUM*NUM+NUM)/2
      NATM = NAT
      NESFT=NESF
C
      CALL VALFM(LOADFM)
      LD     = LOADFM + 1
      LDRSLV = LD     + L2
      LQNDER = LDRSLV + 3*NATM
      LQEDER = LQNDER + 3*NATM
      LDRCNT = LQEDER + 3*NATM
      LDRRAD = LDRCNT + NESFT*NATM*3*3
      LQT    = LDRRAD + NESFT*NATM*3
      LAST   = LQT    + NTS
      NEED   = LAST   - LOADFM -1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 600
C     IF(NESF.GT.NESFP) THEN
C
C        READ PREVIOUSLY COMPUTED QUANTITIES SAVED ON DISK:
C        DERCENTR, DERRAD
C
      CALL SEQREW(NFT26)
      CALL SEQADV(NFT26)
      CALL SEQADV(NFT26)
      LEN = NESFT*NATM*3*3
      CALL SQREAD(NFT26,XX(LDRCNT),LEN)
      LEN = NESFT*NATM*3
      CALL SQREAD(NFT26,XX(LDRRAD),LEN)
      CALL SEQREW(NFT26)
C     END IF
C
      CALL IEFBEM(XX(LD),XX(LDRSLV),XX(LQNDER),XX(LQEDER),XX(LDRCNT),
     *            XX(LDRRAD),XX(LQT),L2,NATM,NESFT,XCTS,YCTS,ZCTS,AS,
     *            QSN,QSE,Q_FS,Q_IND,XE,YE,ZE,RE,ISPHE)
C
  600 CONTINUE
      CALL RETFM(NEED)
C
      IF(MASWRK) WRITE(IW,*)
      IF(MASWRK) WRITE(IW,*)
     *      '.... DONE WITH IEF CONTRIBUTION TO GRADIENT ....'
C
      CALL TIMIT(1)
      RETURN
      END
C*MODULE PCMIEF  *DECK IEFBEM
      SUBROUTINE IEFBEM(DEN,DERSOLV,QNDER,QEDER,DERCENTR,DERRAD,QT,L2,
     *                  NATM,NESFT,XCTS,YCTS,ZCTS,AS,QSN,QSE,Q_FS,Q_IND,
     *                  XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050)
C
      DIMENSION DEN(L2),DERSOLV(3,NATM),QNDER(3,NATM),QEDER(3,NATM),
     *          DERCENTR(NESFT,NATM,3,3),DERRAD(NESFT,NATM,3),QT(NTS),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),QSE(*),Q_FS(*),
     *          Q_IND(*),XE(*),YE(*),ZE(*),RE(*),ISPHE(*)
C
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
      COMMON /PCMTMP/ IDSP(MXATM+MXFRG*5),DPCMFRG(6,MXFRG),
     *                EGPCM(3,MXATM),PCMCDR(3,MXATM)
C
      DATA ZERO, PT5, ONE /0.0D+00, 0.5D+00, 1.0D+00/
      DATA FPI/12.56637062D+00/
C---------------------------------------------------------------------
C
C     CALCOLA LA DERIVATA DEGLI ELEMENTI DI MATRICE BEM RISPETTO ALLE
C     COORDINATE (ICOORD) DEGLI ATOMI DI SOLUTO (NSJ). INDICANDO CON
C     G(X) LA DERIVATA DI G RISPETTO A X, SI HA {JCP, 101, 3888 (1994)}:
C     G(X) = 1/2 TR[PH'(X)] + 1/2 TR[PG'(X)(P)] - TR[S(X)W(P)] + V'(X)NN
C
C---------------------------------------------------------------------
C
      FACT=FPI*EPS/(EPS-ONE)
C
C     CONTROLLA LE DIMENSIONI DELLA BASE
C
C     LEGGE LA MATRICE DENSITA'
C
      CALL DAREAD(IDAF,IODA,DEN,L2,16,0)
C
C     1) PRELIMINARE
C     GRADIENTE DELL'ENERGIA EFC CON LE CARICHE VIRTUALI NUCLEARI
C     FISSE: CALCOLATO CON UNA ROUTINE (MODIFICATA) DI HONDO8,
C     USATO PER LA DERIVATA DI J
C
      IF(IP_F.EQ.1) THEN
        DO ITS=1,NTS
         QT(ITS)=QSN(ITS)+Q_FS(ITS)+Q_IND(ITS)
        ENDDO
        INDQ=1
        CALL CHGDER(DEN,QT,QNDER,QEDER,L2,NATM,INDQ,XCTS,YCTS,ZCTS)
      ELSE
        INDQ=1
        CALL CHGDER(DEN,QSN,QNDER,QEDER,L2,NATM,INDQ,XCTS,YCTS,ZCTS)
      END IF
C
C     2) PRELIMINARE
C     GRADIENTE DELL'ENERGIA EFC CON LE CARICHE VIRTUALI ELETTRONICHE
C     FISSE: CALCOLATO CON UNA ROUTINE (MODIFICATA) DI HONDO8,
C     USATO PER LA DERIVATA DI X
C
      INDQ=2
      CALL CHGDER(DEN,QSE,QNDER,QEDER,L2,NATM,INDQ,XCTS,YCTS,ZCTS)
C
C     LOOP SUGLI ATOMI E SULLE COORDINATE
C
      DO 100 NSJ = 1, NAT
        NSJR=NSJ
        DO 100 ICOORD = 1, 3
          SESE=ZERO
          SNSN=ZERO
          SESN=ZERO
          IF(IP_F.EQ.1) THEN
             CALL OVER(NSJ,NSJR,ICOORD,NATM,NESFT,QSE,QT,SESE,SNSN,SESN,
     *              DERCENTR,DERRAD,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,RE,ISPHE)
          ELSE
             CALL OVER(NSJ,NSJR,ICOORD,NATM,NESFT,QSE,QSN,SESE,SNSN,SESN
     *             ,DERCENTR,DERRAD,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,RE,ISPHE)
          END IF
C
          DERJ = QNDER(ICOORD,NSJ)
          DERX = QEDER(ICOORD,NSJ)
C
          DERY = ZERO
          DERU = ZERO
C
          DO 200 ITS = 1, NTS
C
C           PUNTO RAPPRESENTATIVO
C
            XI = XCTS(ITS)
            YI = YCTS(ITS)
            ZI = ZCTS(ITS)
C
              XN = C(1,NSJ)
              YN = C(2,NSJ)
              ZN = C(3,NSJ)
              DIST = SQRT( (XN-XI)**2+(YN-YI)**2+(ZN-ZI)**2 )
              IF(ICOORD.EQ.1)PROD=XI-XN
              IF(ICOORD.EQ.2)PROD=YI-YN
              IF(ICOORD.EQ.3)PROD=ZI-ZN
              DVNUC = ZAN(NSJ) * PROD / DIST**3
C
            DERY = DERY+DVNUC * QSE(ITS)
C
            IF(IP_F.EQ.1) THEN
             DERU = DERU+DVNUC * QT(ITS)
            ELSE
             DERU = DERU+DVNUC * QSN(ITS)
            END IF
C
 200      CONTINUE
C
         DOVERS = FACT*(PT5*(SESE+SNSN)+SESN)
C
        DEDN   = DERJ+DERY
        DERSOLV(ICOORD,NSJ) = DEDN+DERX+DERU + DOVERS
C       IF(MASWRK) WRITE(IW,33) NSJ,ICOORD,DEDN+DERX+DERU,DOVERS,
C    *               DERSOLV(ICOORD,NSJ)
C33     FORMAT(2I3,3F10.6)
C
 100  CONTINUE
C
C     IF(MASWRK) CALL TIMIT(1)
C
C
C     -- CAVITATION, DISPERSION & REPULSION CONTRIBUTIONS
C
      DO IAT=1,NAT
      DO IXYZ=1,3
         EG(IXYZ,IAT)=EG(IXYZ,IAT)+DERSOLV(IXYZ,IAT)
         IF(ICAV.EQ.1 .OR. IDISP.EQ.1)
     *        EG(IXYZ,IAT)=EG(IXYZ,IAT)+PCMCDR(IXYZ,IAT)
      ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE PCMIEF  *DECK OVER
      SUBROUTINE OVER(NSJ,NSJR,IC,NATM,NESFT,QSE,QSN,SESE,SNSN,SESN,
     *              DERCENTR,DERRAD,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,RE,ISPHE)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION DERCENTR(NESFT,NATM,3,3),DERRAD(NESFT,NATM,3),
     *          QSE(*),QSN(*),XCTS(*),YCTS(*),ZCTS(*),AS(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*)
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      DATA ZERO/0.0D+00/
C
      SESE=ZERO
      SNSN=ZERO
      SESN=ZERO
C
      DO ITS=1,NTS
       L=ISPHE(ITS)
       XNI = - (XE(L) - XCTS(ITS)) / RE(L)
       YNI = - (YE(L) - YCTS(ITS)) / RE(L)
       ZNI = - (ZE(L) - ZCTS(ITS)) / RE(L)
       IF(L.EQ.NSJR) THEN
         IF(IC.EQ.1)DN=XNI
         IF(IC.EQ.2)DN=YNI
         IF(IC.EQ.3)DN=ZNI
         SESE=SESE+DN*(QSE(ITS)**2)/AS(ITS)
         SNSN=SNSN+DN*(QSN(ITS)**2)/AS(ITS)
         SESN=SESN+DN*QSE(ITS)*QSN(ITS)/AS(ITS)
       ELSE
        DCENTN=XNI*DERCENTR(L,NSJ,IC,1)+
     *         YNI*DERCENTR(L,NSJ,IC,2)+
     *         ZNI*DERCENTR(L,NSJ,IC,3)
        DN=DERRAD(L,NSJ,IC)+DCENTN
        SESE=SESE+DN*(QSE(ITS)**2)/AS(ITS)
        SNSN=SNSN+DN*(QSN(ITS)**2)/AS(ITS)
        SESN=SESN+DN*QSE(ITS)*QSN(ITS)/AS(ITS)
       END IF
      ENDDO
C
      RETURN
      END
C*MODULE PCMIEF    *DECK ASCIT
      SUBROUTINE ASCIT(ITER,V,Q0,Q1,Q2,Q3,D0,QA,DIMAT,QREP,TMP,TMP1,
     *                 TMP2,IPVT,RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,
     *                 QSE,XE,YE,ZE,RE,ISPHE,LIST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION V(NTS),Q0(NTS),Q1(NTS),Q2(NTS),Q3(NTS),D0(NTS),QA(NTS),
     *          DIMAT(MXDII1,MXDII1),QREP(NTSPAR,MXDIIS,2),TMP(NTS,3),
     *          TMP1(MXDII1),TMP2(MXDII1,MXDII1),IPVT(MXDII1),
     *          RMUL(MXSP,10),XCTS(*),YCTS(*),ZCTS(*),AS(*),QSE(*),
     *          XE(*),YE(*),ZE(*),RE(*),ISPHE(*),LIST(*)
C
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
C     THIS IS THE MAIN DRIVER FOR THE ASC ITERATIVE CALCULATION
C     V: MOLECULAR FIELD (OR POTENTIAL)
C
C     MOLECULAR REGION INITIALIZATION
C
      IF(ITER.EQ.1) CALL REGINI(ISPHE,LIST)
C
C     ITERATIVE CYCLE (STEP 1)
C
      CALL ASCCYC(1,ITER,V,Q0,Q1,D0,QA,DIMAT,QREP,TMP,TMP1,TMP2,IPVT,
     *            RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,RE,
     *            ISPHE,LIST)
C
C  C-PCM ONLY NEEDS STEP 1.
      IF(IEF.EQ.10) THEN
        CALL DCOPY(NTS,Q1,1,QSE,1)
        CALL DCOPY(NTS,Q1,1,Q0,1)
        RETURN
      END IF
C
C     ITERATIVE CYCLE (STEP 2)
C
C     WRITE THE INTERMEDIATE POTENTIAL, THE STEP 2 INITIAL GUESS
C     AND COPY THE INTERMEDIATE CHARGES
C
      FACIP=-1.0D+00/(EPS-1.0D+00)
      FACIG=-1.0D+00/EPS
      DO ITS=1,NTS
         Q0(ITS)=      Q2(ITS)
         Q2(ITS)=FACIP*Q1(ITS)
         Q3(ITS)=      Q1(ITS)
         IF(ITER.EQ.1) Q0(ITS)=FACIG*Q1(ITS)
      ENDDO
C
      CALL ASCCYC(2,ITER,Q2,Q0,Q1,D0,QA,DIMAT,QREP,TMP,TMP1,TMP2,IPVT,
     *            RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,RE,
     *            ISPHE,LIST)
C
C     COPY THE CHARGES
C
      DO ITS=1,NTS
         QSE(ITS)=Q3(ITS)+Q1(ITS)
      ENDDO
C
C     USE Q1 AND Q3 FOR NEXT GUESSES
C
      CALL DCOPY(NTS,Q3,1,Q0,1)
      CALL DCOPY(NTS,Q1,1,Q2,1)
CONTROLLO
      RETURN
      END
C*MODULE PCMIEF  *DECK REGINI
      SUBROUTINE REGINI(ISPHE,LIST)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION ISPHE(*),LIST(*)
C
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
C
C     REGION INITIALIZATION
C
C     REGION INDECES
C
      K=0
      L=0
      DO I=1,NTS
         J=ISPHE(I)
         IF (J.NE.K) THEN
  10         CONTINUE
             IF (J.EQ.(K+1)) THEN
            K=K+1
            L=L+1
            LIST(L)=I
             ELSE
            K=K+1
                GO TO 10
             END IF
         END IF
      ENDDO
      LIST(L+1)=NTS+1
      NREG=L
C     IF(MASWRK) WRITE(IW,*) ' NREG=',NREG
C     WRITE(*,*) ' NREG=',NREG
      RETURN
      END
C*MODULE PCMIEF  *DECK ASCCYC
      SUBROUTINE ASCCYC(ISTEP,ICALL,V,QIN,QOUT,D0,QA,DIMAT,QREP,TMP,TMP1
     *                 ,TMP2,IPVT,RMUL,MXDII1,NTSPAR,XCTS,YCTS,ZCTS,AS,
     *                  XE,YE,ZE,RE,ISPHE,LIST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION V(NTS),QIN(NTS),QOUT(NTS),D0(NTS),QA(NTS),
     *          DIMAT(MXDII1,MXDII1),QREP(NTSPAR,MXDIIS,2),TMP(NTS,3),
     *          TMP1(MXDII1),TMP2(MXDII1,MXDII1),IPVT(MXDII1),
     *          RMUL(MXSP,10),XCTS(*),YCTS(*),ZCTS(*),AS(*),XE(*),YE(*),
     *          ZE(*),RE(*),ISPHE(*),LIST(*)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      INTEGER ISTEP
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IREP,IEF,IP_F,NFMOPCM
C
      DATA FPI,EPI/12.56637061D+00,25.13274122D+00/
C
C     COMMON /CONV  / ACURCY,ENUCR,AJUNK(3),DIFF,ITERRO,ICALCP,ICBET
C     COMMON /NTNOPT/ ENERGY,ENERG0,DEMAX,SQCDF,ITERMC,MICRO,NOTRF
C     COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
C    *                MPLEVL,MPCTYP
C     DATA RMC/8HMCSCF   /, GVB/8HGVB     /
C     IF(SCFTYP.EQ.RMC) THEN
C       ITER=ITERMC-1
C     ELSE IF(SCFTYP.EQ.GVB) THEN
C       ITER=ITERRO+1
C     ELSE
C       ITER=ITERRO
C     END IF
C
C     INITIAL MATRIX
C
      COST=1.07D+00
      FACEPS=(EPS+1.0D+00)/(2.0D+00*(EPS-1.0D+00))
      IF (ISTEP.EQ.1) THEN
         DO ITS=1,NTS
            D0(ITS) = 1.0D+00/(COST*SQRT(FPI/AS(ITS)))
         ENDDO
      ELSE
         DO ITS=1,NTS
            L=ISPHE(ITS)
            D0(ITS) = 1.0D+00/
     *         (FACEPS + COST*SQRT(FPI*AS(ITS))/(EPI*RE(L)))
         ENDDO
      END IF
C
C     INITIAL GUESS:AT THE FIRST CALL OF THE FIRST STEP
C     IN THE OTHER CASES IS EXTERNAL PROVIDED IN QIN
C
C
      IF(ISTEP.EQ.1) THEN
          IF (ICALL.EQ.1) THEN
            DO ITS=1,NTS
              QIN(ITS)=0.050D+00*V(ITS)*D0(ITS)
            ENDDO
          END IF
      END IF
C
C     ITERATIVE CYCLE
C
      IF(ISTEP.EQ.1)MXITER=MXITR1
      IF(ISTEP.EQ.2)MXITER=MXITR2
      DO I=1,MXITER
C
C     CHARGE-CHARGE INTERACTION
C
C         PRINT CPU TIME
C     IF(MASWRK.AND.I.EQ.1.AND.ITER.EQ.1) THEN
C     WRITE(IW,*) ' '
C     CALL TIMIT(1)
C     WRITE(IW,*) ' '
C     WRITE(IW,*) 'CALLING ASCPOT ...'
C     END IF
C
         CALL ASCPOT(ISTEP,QIN,QA,RMUL,TMP,XCTS,YCTS,ZCTS,AS,XE,YE,ZE,
     *               RE,ISPHE,LIST)
C
C     IF(MASWRK.AND.I.EQ.1.AND.ITER.EQ.1) THEN
C     CALL TIMIT(1)
C     WRITE(IW,*) 'ASCPOT IS DONE'
C     WRITE(IW,*) ' '
C     END IF
C
C     NEW CHARGES
C
         DO ITS=1,NTS
            QOUT(ITS)=(V(ITS)-QA(ITS))*D0(ITS)
         ENDDO
C
C     TEST CONVERGENCE
C
         SDEV=0.0D+00
         DO ITS=1,NTS
            SDEV=SDEV+(QIN(ITS)-QOUT(ITS))**2
         ENDDO
         SDEV=SQRT(SDEV/NTS)
         IF(NFMOPCM.NE.0.AND.MASWRK) WRITE(IW,9100) I,SDEV
         IF (SDEV.LT.THRES) THEN
            IF(NFMOPCM.NE.0.AND.MASWRK) WRITE(IW,9200) ISTEP,I
            CALL DCOPY(NTS,QIN,1,QOUT,1)
            RETURN
         END IF
C
C     INTERPOLATION
C
         IF(MXDIIS.GT.0)
     *    CALL ASCDII(NTSPAR,I,MXDIIS,QOUT,QIN,
     *                DIMAT,QREP,TMP1,TMP2,IPVT)
C
C     COPY THE CHARGES
C
         IF(I.NE.MXITER) CALL DCOPY(NTS,QOUT,1,QIN,1)
      ENDDO
C
      IF(MASWRK) THEN
         WRITE(IW,*) ' '
         WRITE(IW,*) '! ! ! ! WARNING ! ! ! !'
         WRITE(IW,*) 'PCM ITERATIVE PROCEDURE NOT CONVERGED IN',
     *                MXITER,' ITERATIONS'
C        WRITE(IW,*) 'CONVERGED TO ',SDEV,' EXPECTED TO ',THRES
         WRITE(IW,*) 'GAMESS WILL PROCEED DESPITE THIS ...'
         WRITE(IW,*) ' '
      END IF
C     CALL ABRT
      RETURN
 9100 FORMAT(1X,'ASC ITERATION=',I4,' CONVERGED TO',E14.6)
 9200 FORMAT(/1X,'THE ASC EQUATIONS FOR PHASE=',I1,' CONVERGED IN',
     *           I4,' ITERATIONS.'/)
      END
C*MODULE PCMIEF  *DECK ASCPOT
      SUBROUTINE ASCPOT(ISTEP,QIN,POT,RMUL,TMP,XCTS,YCTS,ZCTS,AS,XE,YE,
     *                  ZE,RE,ISPHE,LIST)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR, DSKWRK, MASWRK, NXT
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      DIMENSION QIN(NTS),POT(NTS),RMUL(MXSP,10),TMP(NTS,3),XCTS(*),
     *          YCTS(*),ZCTS(*),AS(*),XE(*),YE(*),ZE(*),RE(*),ISPHE(*),
     *          LIST(*)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
C
      DATA FPI/12.56637061D+00/,ZERO/0.0D+00/
C
C     CLEAN POT
C
      CALL VCLR(POT,1,NTS)
C
C     CALCULATE REGION MULTIPOLES
C
      CALL REGMUL(RMUL,QIN,TMP,XCTS,YCTS,ZCTS,XE,YE,ZE,ISPHE,LIST)
C
      IF(ISTEP.NE.1) THEN
         DO ITS=1,NTS
            L=ISPHE(ITS)
            TMP(ITS,1)=(XCTS(ITS)-XE(L))/RE(L)
            TMP(ITS,2)=(YCTS(ITS)-YE(L))/RE(L)
            TMP(ITS,3)=(ZCTS(ITS)-ZE(L))/RE(L)
         ENDDO
      ENDIF
C
C     CHARGE-CHARGE INTERACTION VIA REGION SCHEME
C
      NXT = IBTYP.EQ.1
      NEXT  = -1
      IPCOUNT = -1
C
      R1=(RCUT(1)*ANTOAU)**2
      R2=(RCUT(2)*ANTOAU)**2
C
      DO 100 I=1,NREG
         IF(GOPARR) THEN
            IPCOUNT = IPCOUNT + 1
            IF (NXT) THEN
               IF(IPCOUNT.GT.NEXT) CALL DDI_DLBNEXT(NEXT)
               IF (NEXT.NE.IPCOUNT) GO TO 100
            ELSE
               IF(MOD(IPCOUNT,NPROC).NE.ME) GO TO 100
            ENDIF
         END IF
         LISTI=LIST(I)
         LISTIP=LIST(I+1)-1
C
C     --- SELF-REGION INTERACTION (WITHIN REGION)
C
         IF (ISTEP.EQ.1) THEN
            DO 110 J=LISTI,LISTIP
               DO 120 K=J+1,LISTIP
                  RJKM1=SQRT((XCTS(J)-XCTS(K))**2+
     &                       (YCTS(J)-YCTS(K))**2+(ZCTS(J)-ZCTS(K))**2)
                  POT(J)=POT(J)+QIN(K)/RJKM1
                  POT(K)=POT(K)+QIN(J)/RJKM1
 120           CONTINUE
 110        CONTINUE
         ELSE
            DO J=LISTI,LISTIP
               DO K=J+1,LISTIP
                  XD=XCTS(J)-XCTS(K)
                  YD=YCTS(J)-YCTS(K)
                  ZD=ZCTS(J)-ZCTS(K)
                  RJKM1=SQRT(XD**2+YD**2+ZD**2)
                  RJKM3=RJKM1**3
                  SCALJ=XD*TMP(J,1)+YD*TMP(J,2)+ZD*TMP(J,3)
                  SCALK=XD*TMP(K,1)+YD*TMP(K,2)+ZD*TMP(K,3)
                  POT(J)=POT(J)-QIN(K)/RJKM3*SCALJ
                  POT(K)=POT(K)+QIN(J)/RJKM3*SCALK
               ENDDO
            ENDDO
         END IF
         XEI=XE(ISPHE(LIST(I)))
         YEI=YE(ISPHE(LIST(I)))
         ZEI=ZE(ISPHE(LIST(I)))
C
C        NO I,J PERMUTATION SYMMETRY IS USED IN THE LOOPS. SHOULD IT BE
C        POSSIBLE TO HALVE WORK HERE?
C
         DO 200 J=1,NREG
            IF (I.EQ.J) GOTO 200
            XEJ=XE(ISPHE(LIST(J)))
            YEJ=YE(ISPHE(LIST(J)))
            ZEJ=ZE(ISPHE(LIST(J)))
            RIJ=(XEI-XEJ)**2+(YEI-YEJ)**2+(ZEI-ZEJ)**2
            LISTJ=LIST(J)
            LISTJP=LIST(J+1)-1
C
C           --- REGION-REGION FULL INTERACTION (NEAR RANGE)
C
            IF (RIJ.LT.R1) THEN
               DO II=LISTI,LISTIP
                  POTI=ZERO
                  IF (ISTEP.EQ.1) THEN
                     DO KK=LISTJ,LISTJP
                        RJKM1=SQRT((XCTS(II)-XCTS(KK))**2+
     &                    (YCTS(II)-YCTS(KK))**2+(ZCTS(II)-ZCTS(KK))**2)
                        POTI=POTI+QIN(KK)/RJKM1
                     ENDDO
                  ELSE
                     DO KK=LISTJ,LISTJP
                        XD=XCTS(II)-XCTS(KK)
                        YD=YCTS(II)-YCTS(KK)
                        ZD=ZCTS(II)-ZCTS(KK)
                        SCAL=XD*TMP(II,1)+YD*TMP(II,2)+ZD*TMP(II,3)
                        RJKM1=SQRT(XD**2+YD**2+ZD**2)
                        POTI=POTI-QIN(KK)*SCAL/RJKM1**3
                     ENDDO
                  END IF
                  POT(II)=POT(II)+POTI
               ENDDO
C
C     --- REGION-REGION MULTIPOLE INTERACTION (MIDDLE AND LONG RANGE)
C
            ELSE
C
               IF (IMUL.LE.0) GO TO 50
C
               DO 40 II=LISTI,LISTIP
                  XD=XCTS(II)-XEJ
                  YD=YCTS(II)-YEJ
                  ZD=ZCTS(II)-ZEJ
                  R=SQRT(XD**2+YD**2+ZD**2)
                  DINVR=1.0D+00/R
                  POTI=ZERO
C
C                 CHARGE - REGION MONOPOLE
C
                  IF (ISTEP.EQ.1) THEN
                     POTI=POTI+RMUL(J,1)*DINVR
                  ELSE
                     SCAL=XD*TMP(II,1)+YD*TMP(II,2)+ZD*TMP(II,3)
                     POTI=POTI-RMUL(J,1)*SCAL*DINVR**3
                  ENDIF
C
                  IF (IMUL.LE.1.OR.RIJ.GE.R2) GO TO 30
C
C                 CHARGE - REGION DIPOLE
C
                  IF (ISTEP.EQ.1) THEN
                     POTI=POTI+
     &                 (RMUL(J,2)*XD+RMUL(J,3)*YD+RMUL(J,4)*ZD)*DINVR**3
                  ELSE
                     SCAL1=RMUL(J,2)*XD+RMUL(J,3)*YD+RMUL(J,4)*ZD
                     SCAL2=RMUL(J,2)*TMP(II,1)+RMUL(J,3)*TMP(II,2)+
     &                     RMUL(J,4)*TMP(II,3)
                     POTI=POTI-(3.0D+00*SCAL1*SCAL-SCAL2*R*2)*DINVR**5
                  ENDIF
C
                  IF (IMUL.LE.2) GO TO 30
C
C                CHARGE - REGION QUADRUPOLE
C
                  IF (ISTEP.EQ.1) THEN
                     POTI=POTI + DINVR**5 * (
     &                ( RMUL(J, 5)*XD**2+
     &                  RMUL(J, 8)*YD**2+
     &                  RMUL(J,10)*ZD**2 )/6.0D+00  +
     &                ( RMUL(J, 6)*XD*YD+
     &                  RMUL(J, 7)*XD*ZD+
     &                  RMUL(J, 9)*YD*ZD )/3.0D+00   )
                  ELSE
                     EQX=RMUL(J,5)*XD+RMUL(J,6)*YD+RMUL(J,7)*ZD
                     EQY=RMUL(J,6)*XD+RMUL(J,8)*YD+RMUL(J,9)*ZD
                     EQZ=RMUL(J,7)*XD+RMUL(J,9)*YD+RMUL(J,10)*ZD
                     POTI=POTI-2.5D+00*SCAL*DINVR**7+
     &                 (EQX*TMP(II,1)+EQY*TMP(II,2)+EQZ*TMP(II,3))
     &                 *DINVR**5
                  ENDIF
 30               CONTINUE
                  POT(II)=POT(II)+POTI
 40            CONTINUE
 50            CONTINUE
            ENDIF
 200     CONTINUE
 100  CONTINUE
C
      IF(GOPARR) THEN
         IF(NXT) CALL DDI_DLBRESET
         CALL DDI_GSUMF(2470,POT,NTS)
      ENDIF
C
C     SURFACE AREA WEIGHTING (ONLY FOR STEP 2)
C
      IF (ISTEP.EQ.2) THEN
         DINVFPI=1.0D+00/FPI
         DO ITS=1,NTS
            POT(ITS)=-AS(ITS)*DINVFPI*POT(ITS) !!!!!!!!
         ENDDO
      END IF
      RETURN
      END
C*MODULE PCMIEF  *DECK REGMUL
      SUBROUTINE REGMUL(RMUL,QIN,TMP,XCTS,YCTS,ZCTS,XE,YE,ZE,ISPHE,LIST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION RMUL(MXSP,10),TMP(NTS,3),QIN(NTS),XCTS(*),YCTS(*),
     *          ZCTS(*),XE(*),YE(*),ZE(*),ISPHE(*),LIST(*)
C
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      DATA ZERO/0.0D+00/
C
C     REGION MULTIPOLES UP TO QUADRUPOLES
C
      CALL VCLR(RMUL,1,MXSP*10)
C
      IF (IMUL.LE.0) GO TO 8000
C
C
C
      DO I=1,NREG
         SUM=ZERO
         DO J=LIST(I),LIST(I+1)-1
            SUM=SUM+QIN(J)
         ENDDO
         RMUL(I,1)=SUM
      ENDDO
C
      IF (IMUL.LE.1) GO TO 8000
C
C     DISTANCES RELATIVES TO THE REGION CENTER
C
      DO ITS=1,NTS
         L=ISPHE(ITS)
         TMP(ITS,1)=XCTS(ITS)-XE(L)
         TMP(ITS,2)=YCTS(ITS)-YE(L)
         TMP(ITS,3)=ZCTS(ITS)-ZE(L)
      ENDDO
C
C     DIPOLES
C
      DO I=1,NREG
         SUMX=ZERO
         SUMY=ZERO
         SUMZ=ZERO
         DO J=LIST(I),LIST(I+1)-1
            SUMX=SUMX+QIN(J)*TMP(J,1) !DX
            SUMY=SUMY+QIN(J)*TMP(J,2) !DY
            SUMZ=SUMZ+QIN(J)*TMP(J,3) !DZ
         ENDDO
         RMUL(I,2)=SUMX
         RMUL(I,3)=SUMY
         RMUL(I,4)=SUMZ
      ENDDO
C
      IF (IMUL.LE.2) GO TO 8000
C
C     QUADRUPOLES
C
      DO I=1,NREG
         SUMXX=ZERO
         SUMXY=ZERO
         SUMXZ=ZERO
         SUMYY=ZERO
         SUMYZ=ZERO
         SUMZZ=ZERO
         DO J=LIST(I),LIST(I+1)-1
            R2=TMP(J,1)**2+TMP(J,2)**2+TMP(J,3)**2
            QJ=QIN(J)
            SUMXX=SUMXX+QJ*(3.0D+00*TMP(J,1)**2-R2)  !QXX
            SUMXY=SUMXY+QJ*TMP(J,1)*TMP(J,2)         !QXY
            SUMXZ=SUMXZ+QJ*TMP(J,1)*TMP(J,3)         !QXZ
            SUMYY=SUMYY+QJ*(3.0D+00*TMP(J,2)**2-R2)  !QYY
            SUMYZ=SUMYZ+QJ*TMP(J,2)*TMP(J,3)         !QYZ
            SUMZZ=SUMZZ+QJ*(3.0D+00*TMP(J,3)**2-R2)  !QZZ
         ENDDO
         RMUL(I, 5)=SUMXX
         RMUL(I, 6)=SUMXY*3.0D+00
         RMUL(I, 7)=SUMXZ*3.0D+00
         RMUL(I, 8)=SUMYY
         RMUL(I, 9)=SUMYZ*3.0D+00
         RMUL(I,10)=SUMZZ
      ENDDO
C
 8000 CONTINUE
C        IF(MASWRK) WRITE(*,*) 'PCM MULTIPOLAR EXPANSION'
C     DO I=1,NREG
C        IF(MASWRK) WRITE(*,9000) I,(RMUL(I,N),N=1,4)
C        IF(MASWRK) WRITE(*,9001) (RMUL(I,N),N=5,10)
C     ENDDO
C
      RETURN
      END
C*MODULE PCMIEF  *DECK ASCDII
      SUBROUTINE ASCDII(NTSPAR,NIT,MXDIIS,QOUT,QIN,
     *                  DIMAT,QREP,TMP,TMPMAT,IPVT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DIMAT(MXDIIS+1,MXDIIS+1),QIN(NTS)
      DIMENSION QOUT(NTS),QREP(NTSPAR,MXDIIS,2)
      DIMENSION TMP(MXDIIS+1),TMPMAT(MXDIIS+1,MXDIIS+1)
      DIMENSION IPVT(MXDIIS+1)
C
      PARAMETER (ONE=1.0D+00)
C
      LOGICAL GOPARR, DSKWRK, MASWRK
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
C     QREP ELEMENTS ARE NO LONGER EXTRAVAGANTLY COPIED.
C     INSTEAD THEY ARE FILLED IN AND USED ON THE ROTARY BASIS
C     (LIKE A REVOLVER).
C     THE SAME COULD BE BUT IS NOT DONE WITH DIMAT, BECAUSE IT IS SMALL
C     AND ONE CAN CREATE SPACE BY BRUTE FORCE SHIFTING.
C     NIT0 POINTS TO THE CURRENT SLOT IN QREP.
C     NITMAX IS THE NUMBER OF FILLED SLOTS IN QREP.
C     NIT0 AND NITMAX ARE EQUAL TO NIT, UNLESS NIT>MXDIIS, IN THE
C     LATTER CASE NITMAX IS MXDIIS AND NIT0 IS THE SPOT WITHIN
C     MXDIIS (SEE BELOW).  I0 TAKES CARE OF CONVERTING I (LATEST
C     DIIS ITERATION INDEX) INTO THE PROPER SLOT IN QREP.
C
      NITMAX=MIN(NIT,MXDIIS)
      NIT0=MOD(NIT-1,MXDIIS)+1
C
      IMIN=ME*NTSPAR+1
      IMAX=MIN((ME+1)*NTSPAR,NTS)
      NTSME=IMAX-IMIN+1
C
C     STORE THE CHARGES
C     QREP(,,1)=QOUT
C     QREP(,,2)=QOUT-QIN
C
      CALL DCOPY(NTSME,QOUT(IMIN),1,QREP(1,NIT0,1),1)
      CALL DCOPY(NTSME,QOUT(IMIN),1,QREP(1,NIT0,2),1)
      CALL DAXPY(NTSME,-ONE,QIN(IMIN),1,QREP(1,NIT0,2),1)
C
C     UPGRADE THE INTERPOLATION MATRIX
C
      IF(NIT.GT.MXDIIS) THEN
         DO I=1,MXDIIS-1
            DO J=1,MXDIIS-1
               DIMAT(I+1,J+1)=DIMAT(I+2,J+2)
            ENDDO
         ENDDO
      ENDIF
C
      I0=NIT0
      DO I=NITMAX,1,-1
         TMP(I)=DDOT(NTSME,QREP(1,NIT0,2),1,QREP(1,I0,2),1)
         I0=I0-1
         IF(I0.EQ.0) I0=I0+MXDIIS
      ENDDO
      IF(GOPARR) CALL DDI_GSUMF(2452,TMP,NITMAX)
C
      DIMAT(NITMAX+1,1)=-1.0D+00
      DIMAT(1,NITMAX+1)=-1.0D+00
      DO I=NITMAX,1,-1
         DIMAT(NITMAX+1,I+1)=TMP(I)
         DIMAT(I+1,NITMAX+1)=TMP(I)
      ENDDO
C
C     AT THE FIRST ITERATION NO INTERPOLATION ONLY MATRIX INITIALIZATION
C
      IF (NIT.EQ.1) THEN
         DIMAT(1,1)=0.0D+00
         RETURN
      END IF
C
C     VECTOR INITIALIZATION
C
      TMP(1)=-1.0D+00
      CALL VCLR(TMP(2),1,NITMAX)
C
C     SOLVE THE LINEAR SYSTEM
C
C     COPY THE MATRIX (SHOULD BE DESTROYED)
C
      DO I=1,NITMAX+1
         DO J=1,NITMAX+1
            TMPMAT(I,J)=DIMAT(I,J)
         ENDDO
      ENDDO
C
C     SOLVE THE LINEAR SYSTEM
C
      CALL DGEFA(TMPMAT,MXDIIS+1,NITMAX+1,IPVT,INFO)
      IF (INFO.NE.0) THEN
         WRITE(6,*) 'SINGULAR MATRIX IN ASCDII'
         CALL ABRT
      END IF
      CALL DGESL(TMPMAT,MXDIIS+1,NITMAX+1,IPVT,TMP,0)
C
C     INTERPOLATE
C
      CALL VCLR(QOUT,1,NTS)
      I0=NIT0
      DO I=NITMAX,1,-1
         CALL DAXPY(NTSME,TMP(I+1),QREP(1,I0,1),1,QOUT(IMIN),1)
         I0=I0-1
         IF(I0.EQ.0) I0=I0+MXDIIS
      ENDDO
      IF(GOPARR) CALL DDI_GSUMF(2452,QOUT,NTS)
      RETURN
      END
C
C*MODULE PCMIEF  *DECK ADDPOT
      SUBROUTINE ADDPOT(QPOT,VEC_MUL,IPCFP,XCTS,YCTS,ZCTS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO = 0.00D+00, ONE = 1.00D+00)
C
      DIMENSION QPOT(NTS),VEC_MUL(*),XCTS(*),YCTS(*),ZCTS(*)
C
      CHARACTER*8 FRGNME,POLNAM,DPOLNAM
C
      LOGICAL IPCFP
C
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /EFPPAR/ EFP(3,MXFGPT),EFPOL(9,MXFGPT),
     *                EFDP(3,MXIFRQ*MXFGPT),EFDPOL(9,MXIFRQ*MXFGPT),
     *                ENO,DIND(3,MXFGPT),DINDD(3,MXFGPT),
     *                POLNAM(MXFGPT),DPOLNAM(MXIFRQ*MXFGPT),
     *                POLSCR(MXFGPT),IPLSCR,INDSCR
      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 /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
C     THIS ROUTINE MERGES TOGETHER THE ELECTRONIC, NUCLEAR, EFFECTIVE
C     FRAGMENT POTENTIALS (BOTH FIX AND POLARIZABLE)
C
C     -- NUCLEAR  POTENTIAL
C
      DO ITS=1,NTS
         DO IAT=1,NAT
            DINVR=1.0D+00/SQRT(
     *           (XCTS(ITS)-C(1,IAT))**2+
     *           (YCTS(ITS)-C(2,IAT))**2+
     *           (ZCTS(ITS)-C(3,IAT))**2)
            QPOT(ITS)=QPOT(ITS)-ZAN(IAT)*DINVR
         ENDDO
      ENDDO
C
C     -- EFP STATIC MULTIPOLE POTENTIAL
C
      IF(IEFC+IEFD+IEFQ.GT.0) THEN
         CALL VCLR(VEC_MUL,1,NTS)
         DO ITS = 1, NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
            DO I=1,NMTTPT
              ALFA=EFATRM(I)
              BETA=EFBTRM(I)
              XX=EFQAD(1,I)
              YY=EFQAD(2,I)
              ZZ=EFQAD(3,I)
              XY=EFQAD(4,I)
              XZ=EFQAD(5,I)
              YZ=EFQAD(6,I)
              DUM=XX+YY+ZZ
              QXX=(3.0D+00*XX-DUM)/2.0D+00
              QYY=(3.0D+00*YY-DUM)/2.0D+00
              QZZ=(3.0D+00*ZZ-DUM)/2.0D+00
              QXY=1.5D+00*XY
              QXZ=1.5D+00*XZ
              QYZ=1.5D+00*YZ
              DX = XCTS(ITS)-EFC(1,I)
              DY = YCTS(ITS)-EFC(2,I)
              DZ = ZCTS(ITS)-EFC(3,I)
              R2=(EFC(1,I)-XI)**2+(EFC(2,I)-YI)**2+(EFC(3,I)-ZI)**2
              R = SQRT(R2)
              R3= R*R2
              R5= R2*R3
              FACT=ONE-(BETA*EXP(-ALFA*R2))
              VEC_MUL(ITS)=VEC_MUL(ITS)
     *                   + (FACT*EFCHG(1,I) + EFCHG(2,I))/R
     *                   + FACT*(EFDIP(1,I)*DX
     *                         + EFDIP(2,I)*DY
     *                         + EFDIP(3,I)*DZ)/R3
     *                   + FACT*(QXX*DX*DX+QYY*DY*DY+QZZ*DZ*DZ
     *                          +QXY*2.0D+00*DX*DY
     *                          +QXZ*2.0D+00*DX*DZ
     *                          +QYZ*2.0D+00*DY*DZ)/R5
            ENDDO
            QPOT(ITS)=QPOT(ITS)-VEC_MUL(ITS)
         ENDDO
      END IF
C
C     -- POTENTIAL DUE TO EFP INDUCED DIPOLES --
C
      IF (IPCFP) THEN
         DO ITS=1,NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
            DUM=ZERO
            DO LEFP=1,NPTTPT
               XL=EFP(1,LEFP)
               YL=EFP(2,LEFP)
               ZL=EFP(3,LEFP)
               X = XI - XL
               Y = YI - YL
               Z = ZI - ZL
               R2=X*X + Y*Y + Z*Z
               R = SQRT(R2)
               R3 = R*R*R
               DUM= DUM
     *           + (DIND(1,LEFP)*X+DIND(2,LEFP)*Y+DIND(3,LEFP)*Z)/R3
            ENDDO
            QPOT(ITS)=QPOT(ITS)-DUM
         ENDDO
      END IF
      RETURN
C
      END
C*MODULE PCMIEF  *DECK ITIEFIN
      SUBROUTINE ITIEFIN
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      PARAMETER (NNAM=16)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA SOLWD/8HPCMITR  /
      DATA QNAM /8HMXDIIS  ,8HIMUL    ,8HRCUT1   ,8HRCUT2   ,
     *           8HTHRES   ,8HIDIRCT  ,8HMXITR1  ,8HMXITR2  ,
     *           8HIMGABI  ,8HRABI    ,8HIMGASC  ,8HRASC    ,
     *           8HIEFPOL  ,8HREFPOL  ,8HTHRSLS  ,8HDENSLS  /
      DATA KQNAM/1,1,3,3,3,1,1,1,1,3,1,3,1,3,3,3/
C
C     IEF ITERATIVE PARAMETERS: DEFAULTS
C
C     MXDIIS: MAXIMUM SIZE OF THE DIIS LINEAR EQUATIONS,
C     THE VALUE ESTABILISH THE AMOUNT OF MEMORY USED
C
C     MEM= 5*NTS + 2*MXDIIS*NTS + 2*(MXDIIS+1)*(MXDIIS+1) + MXDIIS +1
C                  ************
C
      MXDIIS=50
C
C     MXITR1: MAXIMUM NUMBER OF ITERATIONS IN PHASE 1
C
      MXITR1=50
C
C     MXITR2: MAXIMUM NUMBER OF ITERATIONS IN PHASE 2
C
      MXITR2=50
C
      IF(NFG.NE.0) THEN
        MXITR1=200
        MXITR2=200
      ENDIF
C
C     IMUL: MULTIPOLE REGION EXPANSION ORDER:
C
C      O: NEGLECTED (ONLY FOR TEST POURPOSES)
C      1: MONOPOLE
C      2: MONOPOLE+DIPOLE
C     >2: MONOPOLE+DIPOLE+QUADRUPOLE
C
      IMUL=3
C
C     RCUT: CUTOFF RADII, VALUES IN ANGSTROMS
C     RCUT1: CUTOFF FOR MID RANGE ITERATIONS
C     RCUT2: CUTOFF FOR LONG RANGE ITERATIONS
C
      RCUT1=15.0D+00
      RCUT2=30.0D+00
      IF(NFG.NE.0) THEN
        RCUT1=20.0D+00
        RCUT2=50.0D+00
      ENDIF
C
C     THRES: CONVERGENCE THRESHOLD FOR THE CHARGE DIIS PROCEDURE
C
      THRES=1.0D-08
      IF(NFG.NE.0) THRES=1.0D-09
C
C     IDIRCT: =1, COMPUTE ELECTRIC POTENTIAL EACH TIME, NO DISC WORK
C
      IDIRCT=1
C
C     IMGABI: =1, USE IMAGE OF THE AB INITIO REGION
C
      IMGABI=0
C
C     RABI: THE DISTANCE IMGABI WORKS(ANGSTROMS)
C
      RABI=4.0D+00
C
C     IMGASC: =1, USE IMAGE OF THE PCM ASC
C
      IMGASC=0
C
C     RASC: THE DISTANCE IMGASC WORKS(ANGSTROMS)
C
      RASC=20.0D+00
C
C     IEFPOL: =1, USE PCM ASC TO INDUCE EFP DIPOLES
C
      IEFPOL=1
C
C     REFPOL: THE DISTANCE ASC WORKS(ANGSTROMS) TO INDUCE DIPOLES
C
      REFPOL=0.0D+00
C
C     THRSLS: THE LOOSE CRITERION FOR EARLY SCF CYCLES.
C             IF THRSLS < THRES , NO LOOSE CRITERIA WILL BE USED
C
      THRSLS=5.0D-04
C
C     DENSLS: THE DENSITY CHANGE THRESH FOR LOOSE CRITERIA.
C             IN EARLY SCF CYCLES, IF DENSITY CHANGE IS ABOVE
C             DENSLS, THRSLS WILL BE USED INSTEAD OF THRES.
C             IF DENSLS IS VERY LARGE, ALL THE SCF CYCLES WILL
C             USE THRES, NO LOOSE CRITERIA WILL BE USED
C
      DENSLS=0.010D+00
C
C     READ THE PCMITR NAMELIST
C
      CALL NAMEIO(IR,JRET,SOLWD,NNAM,QNAM,KQNAM,
     *            MXDIIS, IMUL, RCUT1, RCUT2, THRES,
     *            IDIRCT,MXITR1,MXITR2,
     *            IMGABI,RABI,IMGASC,RASC,IEFPOL,REFPOL,
     *            THRSLS,DENSLS,
     *       0,0,0,    0,0,0,0,0,
     *   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0,
     *   0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
C
      RCUT(1)=RCUT1
      RCUT(2)=RCUT2
      IF(MXDIIS.GT.MAX(MXITR1,MXITR2))MXDIIS=MAX(MXITR1,MXITR2)
      RETURN
      END
C*MODULE PCMIEF  *DECK ELEPOT
      SUBROUTINE ELEPOT(BEMPOT,L2,NFT27,XCTS,YCTS,ZCTS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION BEMPOT(L2),XCTS(*),YCTS(*),ZCTS(*)
C
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      CALL SEQREW(NFT27)
      DO  ITS = 1, NTS
        XI=XCTS(ITS)
        YI=YCTS(ITS)
        ZI=ZCTS(ITS)
         CALL INTMEP(BEMPOT,XI,YI,ZI)
         CALL SQWRIT(NFT27,BEMPOT,L2)
      ENDDO
      CALL SEQREW(NFT27)
      RETURN
      END
C*MODULE PCMIEF  *DECK POTIMG
      SUBROUTINE POTIMG(QADD,XI,YI,ZI)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
C
C     COMPUTE THE POTENTIAL AT (XI,YI,ZI) INDUCED
C     BY THE IMAGE OF THE AB INITIO REGION(INCLUDE BUFFER)
C
            QADD=0.0D+00
            DO  I=1,NUMIMG
               XX=QADIMG(1,I)
               YY=QADIMG(2,I)
               ZZ=QADIMG(3,I)
               XY=QADIMG(4,I)
               XZ=QADIMG(5,I)
               YZ=QADIMG(6,I)
               DUM=XX+YY+ZZ
               QXX=(3.0D+00*XX-DUM)/2.0D+00
               QYY=(3.0D+00*YY-DUM)/2.0D+00
               QZZ=(3.0D+00*ZZ-DUM)/2.0D+00
               QXY=1.5D+00*XY
               QXZ=1.5D+00*XZ
               QYZ=1.5D+00*YZ
               DX = XI-CRDIMG(1,I)
               DY = YI-CRDIMG(2,I)
               DZ = ZI-CRDIMG(3,I)
               R2=DX**2 + DY**2 + DZ**2
               R = SQRT(R2)
C    -- MONOPOLE (ELECTRON ONLY)   --
               QADD = QADD + CHGIMG(1,I)/R
C    -- DIPOLE     --
               QADD = QADD +
     *                    (DIPIMG(1,I)*DX+
     *                    DIPIMG(2,I)*DY+
     *                    DIPIMG(3,I)*DZ)/(R*R2)
C    -- QUADROPOLE --
               QADD = QADD + (QXX*3.0D+00*DX*DX
     *                     + QYY*3.0D+00*DY*DY
     *                     + QZZ*3.0D+00*DZ*DZ
     *                     + QXY*6.0D+00*DX*DY
     *                     + QXZ*6.0D+00*DX*DZ
     *                     + QYZ*6.0D+00*DY*DZ)
     *                       /(3.0D+00*R2*R2*R)
            ENDDO
      RETURN
      END
C*MODULE PCMIEF   *DECK IMABIX
      SUBROUTINE IMABIX
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DATA RNONE/8HNONE    /
C
C             MAKE AN IMAGE OF THE AB INITIO REGION(INCLUDE BUFFER)
C
C     FIND THE BOUNDARY OF THE AB INITIO REGION
C
      DO N = 1, 3
        TMAX = C(N,1)
        TMIN = C(N,1)
        DO IAT = 2, NAT
          IF(C(N,IAT).GT.TMAX) TMAX = C(N,IAT)
          IF(C(N,IAT).LT.TMIN) TMIN = C(N,IAT)
        ENDDO
        IF(N.EQ.1) XMAX = TMAX
        IF(N.EQ.1) XMIN = TMIN
        IF(N.EQ.2) YMAX = TMAX
        IF(N.EQ.2) YMIN = TMIN
        IF(N.EQ.3) ZMAX = TMAX
        IF(N.EQ.3) ZMIN = TMIN
      ENDDO
C
C     --- STONE MULTIPOLE EXPANSION
C
      CALL STONE(IP,RNONE)
C
      RETURN
      END
C*MODULE PCMIEF  *DECK CHGINT
      SUBROUTINE CHGINT(QQ,XCTS,YCTS,ZCTS,QSE,LIST,ASCCRD,ASCCHG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,OUT,NORM,DOUBLE
      LOGICAL GOPARR,MASWRK,DSKWRK,USEIMG
C
      DIMENSION QQ(*),XCTS(*),YCTS(*),ZCTS(*),QSE(*),LIST(*),
     *          ASCCRD(3,*),ASCCHG(*)
C
      DIMENSION DIJ(225),XIN(125),YIN(125),ZIN(125),
     *          IX(35),IY(35),IZ(35),JX(35),JY(35),JZ(35),
     *          IJX(225),IJY(225),IJZ(225)
      DIMENSION CHCINT(225)
      DIMENSION FIJ(225)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /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 /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STV   / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
      COMMON /SYMIND/ TOL,II,JJ,LIT,LJT,MINI,MINJ,MAXI,MAXJ,IANDJ
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,
     *           PI212=1.1283791670955D+00, SQRT3=1.73205080756888D+00,
     *           SQRT5=2.23606797749979D+00, SQRT7=2.64575131106459D+00,
     *           RLN10=2.30258D+00)
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
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/
      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
      DATA DBUGME/8HINT1    /
C
C    CALCULATE THE CHARGE-CHARGE CONTRIBUTION TO THE ONE ELECTRON
C    INTEGRALS FOR THE PCM (REGION WISE)
C
C     CHECK ON BASIS SET, FOR THIS ROUTINE IS LIMITED TO SPDFG
C
      CALL BASCHK(LMAX)
      IF(LMAX.GT.4) THEN
         IF(MASWRK) WRITE(IW,*) 'CHGINT: CODE IS SPDFG LIMITED'
         CALL ABRT
      END IF
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
C
      NEFC = NREG
      TOL = RLN10*ITOL
      OUT = NPRINT.EQ.3  .OR.  EXETYP.EQ.DBUGME
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
C
      L1 = NUM
      L2 = (NUM*NUM+NUM)/2
      CALL VCLR(QQ,1,L2)
C
C     ----- ISHELL
C
      DO 720 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     ----- JSHELL
C
         DO 700 JJ = 1,II
C
C     ----- GO PARALLEL! -----
C
            IF (GOPARR) THEN
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 690
            END IF
            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
            NROOTS = (LIT+LJT-2)/2+1
            RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
            IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
            IJ = 0
            MAX = MAXJ
            DO 160 I = MINI,MAXI
               NX = IX(I)
               NY = IY(I)
               NZ = IZ(I)
               IF (IANDJ) MAX = I
               DO 140 J = MINJ,MAX
                  IJ = IJ+1
                  IJX(IJ) = NX+JX(J)
                  IJY(IJ) = NY+JY(J)
                  IJZ(IJ) = NZ+JZ(J)
  140          CONTINUE
  160       CONTINUE
            DO 180 I = 1,IJ
               CHCINT(I) = ZERO
  180       CONTINUE
C
C     ----- I PRIMITIVE
C
            JGMAX = J2
            DO 520 IG = I1,I2
               AI = EX(IG)
               ARRI = AI*RR
               AXI = AI*XI
               AYI = AI*YI
               AZI = AI*ZI
               CSI = CS(IG)
               CPI = CP(IG)
               CDI = CD(IG)
               CFI = CF(IG)
               CGI = CG(IG)
C
C
C     ----- J PRIMITIVE
C
               IF (IANDJ) JGMAX = IG
               DO 500 JG = J1,JGMAX
                  AJ = EX(JG)
                  AA = AI+AJ
                  AA1 = ONE/AA
                  DUM = AJ*ARRI*AA1
                  IF (DUM .GT. TOL) GO TO 500
                  FAC = EXP(-DUM)
                  CSJ = CS(JG)
                  CPJ = CP(JG)
                  CDJ = CD(JG)
                  CFJ = CF(JG)
                  CGJ = CG(JG)
                  AX = (AXI+AJ*XJ)*AA1
                  AY = (AYI+AJ*YJ)*AA1
                  AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
                  DOUBLE=IANDJ.AND.IG.NE.JG
                  MAX = MAXJ
                  NN = 0
                  DUM1 = ZERO
                  DUM2 = ZERO
                  DO 220 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 200 J = MINJ,MAX
                        IF (J.EQ.1) THEN
                           DUM2=DUM1*CSJ
                           IF (DOUBLE) THEN
                              IF (I.LE.1) THEN
                                 DUM2=DUM2+DUM2
                              ELSE
                                 DUM2=DUM2+CSI*CPJ*FAC
                              END IF
                           END IF
                        ELSE IF (J.EQ.2) THEN
                           DUM2=DUM1*CPJ
                           IF (DOUBLE) DUM2=DUM2+DUM2
                        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
                        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
                        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
                        END IF
                        NN = NN+1
                        DIJ(NN) = DUM2
  200                CONTINUE
  220             CONTINUE
C
C     PCM REGIONAL CHARGES INTEGRALS.
C     IF THE REGION IS CLOSE TO AB INITIO ATOMS, DO NOT USE IMGASC
C     IF THE REGION IS FAR AWAY FROM AB INITIO ATOMS, USE IMGASC
C
                  DUM = PI212*AA1
                  DO 400 I=1,IJ
                     FIJ(I) = DIJ(I)*DUM
  400             CONTINUE
                  AAX = AA*AX
                  AAY = AA*AY
                  AAZ = AA*AZ
                  KTIMES = 0
                  DO 492 IC = 1,NEFC
                     CX = ASCCRD(1,IC)
                     CY = ASCCRD(2,IC)
                     CZ = ASCCRD(3,IC)
                     USEIMG = .FALSE.
                     USEIMG = CX.GE.XMAX+RASC*ANTOAU
     *               .OR.     CX.LE.XMIN-RASC*ANTOAU
     *               .OR.     CY.GE.YMAX+RASC*ANTOAU
     *               .OR.     CY.LE.YMIN-RASC*ANTOAU
     *               .OR.     CZ.GE.ZMAX+RASC*ANTOAU
     *               .OR.     CZ.LE.ZMIN-RASC*ANTOAU
C
                  IF(USEIMG) THEN
C                 --- USE IMGASC ---
C
                     KTIMES = KTIMES + 1
                     ZNUC =  -ASCCHG(IC)
                     XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
                     IF (NROOTS.LE.3) CALL RT123
                     IF (NROOTS.EQ.4) CALL ROOT4
                     IF (NROOTS.EQ.5) CALL ROOT5
                     MM = 0
                     DO 477 K = 1,NROOTS
                        UU = AA*U(K)
                        WW = W(K)*ZNUC
                        TT = ONE/(AA+UU)
                        T = SQRT(TT)
                        X0 = (AAX+UU*CX)*TT
                        Y0 = (AAY+UU*CY)*TT
                        Z0 = (AAZ+UU*CZ)*TT
                        IN = -5+MM
                        DO 476 I = 1,LIT
                           IN = IN+5
                           NI = I
                           DO 475 J = 1,LJT
                              JN = IN+J
                              NJ = J
                              CALL STVINT
                              XIN(JN) = XINT
                              YIN(JN) = YINT
                              ZIN(JN) = ZINT*WW
  475                      CONTINUE
  476                   CONTINUE
                        MM = MM+25
  477                CONTINUE
                     DO 481 I = 1,IJ
                        NX = IJX(I)
                        NY = IJY(I)
                        NZ = IJZ(I)
                        DUM = ZERO
                        MM = 0
                        DO 479 K = 1,NROOTS
                           DUM = DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
                           MM = MM+25
  479                   CONTINUE
                        CHCINT(I) = CHCINT(I)+DUM*FIJ(I)
  481                CONTINUE
                  ELSE
C
C                 --- USE EVEREY TESSERA IN EACH REGION ---
C
                  DO 592 JC=LIST(IC),LIST(IC+1)-1
                     CX = XCTS(JC)
                     CY = YCTS(JC)
                     CZ = ZCTS(JC)
                     ZNUC =  -QSE(JC)
C
                     XX = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
                     IF (NROOTS.LE.3) CALL RT123
                     IF (NROOTS.EQ.4) CALL ROOT4
                     IF (NROOTS.EQ.5) CALL ROOT5
                     MM = 0
                     DO 577 K = 1,NROOTS
                        UU = AA*U(K)
                        WW = W(K)*ZNUC
                        TT = ONE/(AA+UU)
                        T = SQRT(TT)
                        X0 = (AAX+UU*CX)*TT
                        Y0 = (AAY+UU*CY)*TT
                        Z0 = (AAZ+UU*CZ)*TT
                        IN = -5+MM
                        DO 576 I = 1,LIT
                           IN = IN+5
                           NI = I
                           DO 575 J = 1,LJT
                              JN = IN+J
                              NJ = J
                              CALL STVINT
                              XIN(JN) = XINT
                              YIN(JN) = YINT
                              ZIN(JN) = ZINT*WW
  575                      CONTINUE
  576                   CONTINUE
                        MM = MM+25
  577                CONTINUE
                     DO 581 I = 1,IJ
                        NX = IJX(I)
                        NY = IJY(I)
                        NZ = IJZ(I)
                        DUM = ZERO
                        MM = 0
                        DO 579 K = 1,NROOTS
                           DUM = DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
                           MM = MM+25
  579                   CONTINUE
                        CHCINT(I) = CHCINT(I)+DUM*FIJ(I)
  581                CONTINUE
  592             CONTINUE
                  END IF
  492             CONTINUE
C                 IF(MASWRK.AND.II.EQ.1.AND.JJ.EQ.1.AND.IG.EQ.I1
C    *               .AND.ITER.EQ.1)
C    *            WRITE(IW,*) 'IMGASC USED',KTIMES,'TIMES'
C
C     ----- END PRIMITIVES -----
C
  500          CONTINUE
  520       CONTINUE
C
            MAX = MAXJ
            NN = 0
            DO 620 I = MINI,MAXI
               LI = LOCI+I
               IN = (LI*(LI-1))/2
               IF (IANDJ) MAX = I
               DO 600 J = MINJ,MAX
                  LJ = LOCJ+J
                  JN = LJ+IN
                  NN = NN+1
                  QQ(JN) = QQ(JN) + CHCINT(NN)
  600          CONTINUE
  620       CONTINUE
C
C     ----- END PARALLEL
C
  690    CONTINUE
C
C     ----- END SHELLS -----
C
  700    CONTINUE
  720 CONTINUE
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) CALL DDI_GSUMF(920,QQ,L2)
C
C     ----- DONE WITH INTEGRALS -----
C
      IF(MASWRK.AND.OUT) THEN
         WRITE(IW,*) 'EFC INTEGRALS'
         CALL PRTRIL(QQ,L1)
         WRITE(IW,9100)
         CALL TIMIT(1)
      END IF
      RETURN
C
 9100 FORMAT(1X,'...... END OF ASCCHG INTEGRALS ......')
      END
C*MODULE PCMIEF  *DECK DPLINT
      SUBROUTINE DPLINT(QQ,ASCCRD,ASCDIP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,MASWRK,DSKWRK,USEIMG
      LOGICAL IANDJ,NORM,DOUBLE
      LOGICAL OUT
C
      DIMENSION QQ(*),ASCCRD(3,*),ASCDIP(3,*)
C
      DIMENSION CHDINT(100)
      DIMENSION DIJ(100),XIN(432),YIN(432),ZIN(432)
      DIMENSION FIJ(100)
      DIMENSION IX(20),IY(20),IZ(20),JX(20),JY(20),JZ(20)
      DIMENSION IJX(100),IJY(100),IJZ(100)
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      COMMON /CSSTV / CX,CY,CZ
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /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 /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STV   / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
C
      DATA PI212 /1.1283791670955D+00/
      DATA SQRT3 /1.73205080756888D+00/
      DATA SQRT5 /2.23606797749979D+00/
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 1/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1/
      DATA IX / 1, 5, 1, 1, 9, 1, 1, 5, 5, 1,
     1         13, 1, 1, 9, 9, 5, 1, 5, 1, 5/
      DATA IY / 1, 1, 5, 1, 1, 9, 1, 5, 1, 5,
     1          1,13, 1, 5, 1, 9, 9, 1, 5, 5/
      DATA IZ / 1, 1, 1, 5, 1, 1, 9, 1, 5, 5,
     1          1, 1,13, 1, 5, 1, 5, 9, 9, 5/
C
      DATA DBUGME/8HINT1    /
      DATA ZERO,ONE/0.0D+00,1.0D+00/
      DATA RLN10 /2.30258D+00/
C
C    CALCULATE THE CHARGE-DIPOLE CONTRIBUTION TO THE ONE ELECTRON
C    INTEGRALS FOR THE PCM (REGION WISE)
C
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
      NEFD = NREG
      OUT = NPRINT.EQ.3  .OR.  EXETYP.EQ.DBUGME
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NUM2=(NUM*(NUM+1))/2
      DO 10 I=1,NUM2
      QQ(I) = ZERO
 10   CONTINUE
C
C
C     ----- ISHELL
C
      DO 600 II = 1,NSHELL
      IAT= KATOM(II)
      XI = C(1,IAT)
      YI = C(2,IAT)
      ZI = C(3,IAT)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
C
C     ----- JSHELL
C
      DO 580 JJ = 1,II
C
C     ----- GO PARALLEL! -----
C
      IF (GOPARR) THEN
         IPCOUNT = IPCOUNT + 1
         IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 570
      END IF
      JAT= KATOM(JJ)
      XJ = C(1,JAT)
      YJ = C(2,JAT)
      ZJ = C(3,JAT)
      J1 = KSTART(JJ)
      J2 = J1+KNG(JJ)-1
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      NROOTS = (LIT+LJT+1-2)/2 + 1
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      IJ = 0
      MAX = MAXJ
      DO 50 I = MINI,MAXI
      NX = IX(I)
      NY = IY(I)
      NZ = IZ(I)
      IF (IANDJ) MAX = I
      DO 50 J = MINJ,MAX
      IJ = IJ+1
      IJX(IJ) = NX+JX(J)
      IJY(IJ) = NY+JY(J)
      IJZ(IJ) = NZ+JZ(J)
  50  CONTINUE
      DO 60 I=1,IJ
  60  CHDINT(I) = ZERO
C
C     ----- I PRIMITIVE
C
      JGMAX = J2
      DO 520 IG = I1,I2
      AI = EX(IG)
      ARRI = AI*RR
      AXI = AI*XI
      AYI = AI*YI
      AZI = AI*ZI
      CSI = CS(IG)
      CPI = CP(IG)
      CDI = CD(IG)
      CFI=CF(IG)
C
C
C     ----- J PRIMTIVE
C
      IF (IANDJ) JGMAX = IG
      DO 500 JG = J1,JGMAX
      AJ = EX(JG)
      AA = AI+AJ
      AA1 = ONE/AA
      DUM = AJ*ARRI*AA1
      IF (DUM .GT. TOL) GO TO 500
      FAC =  EXP(-DUM)
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ=CF(JG)
      AX = (AXI+AJ*XJ)*AA1
      AY = (AYI+AJ*YJ)*AA1
      AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
      DOUBLE=IANDJ.AND.IG.NE.JG
      MAX = MAXJ
      NN = 0
      DO 170 I=MINI,MAXI
      GO TO (70,80,110,110,90,110,110,100,110,110,
     1       102,110,110,104,110,110,110,110,110,106),I
   70 DUM1=CSI*FAC
      GO TO 110
   80 DUM1=CPI*FAC
      GO TO 110
   90 DUM1=CDI*FAC
      GO TO 110
  100 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 110
  102 DUM1=CFI*FAC
      GO TO 110
  104 IF(NORM) DUM1=DUM1*SQRT5
      GO TO 110
  106 IF(NORM) DUM1=DUM1*SQRT3
  110 IF(IANDJ) MAX=I
      DO 170 J=MINJ,MAX
      GO TO (125,130,160,160,140,160,160,150,160,160,
     1       152,160,160,154,160,160,160,160,160,156),J
  125 DUM2=DUM1*CSJ
      IF(.NOT.DOUBLE) GO TO 160
      IF(I.GT.1) GO TO 126
      DUM2=DUM2+DUM2
      GO TO 160
  126 DUM2=DUM2+CSI*CPJ*FAC
      GO TO 160
  130 DUM2=DUM1*CPJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  140 DUM2=DUM1*CDJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  150 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 160
  152 DUM2=DUM1*CFJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  154 IF(NORM) DUM2=DUM2*SQRT5
      GO TO 160
  156 IF(NORM) DUM2=DUM2*SQRT3
  160 NN=NN+1
  170 DIJ(NN)=DUM2
C
C     ..... HELLMANN-FEYNMAN TERM .....
C
      DUM = PI212*AA1
      DUM=DUM+DUM
      DO 380 I = 1,IJ
      FIJ(I) = DIJ(I)*DUM
  380 CONTINUE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
      DO 480 IC = 1,NEFD
      CX   = ASCCRD(1,IC)
      CY   = ASCCRD(2,IC)
      CZ   = ASCCRD(3,IC)
      USEIMG = .FALSE.
      USEIMG = (CX.GE.XMAX+RASC*ANTOAU
     *   .AND. CX.LE.XMAX+RASC*ANTOAU*3.0D+00)
     *   .OR.  (CX.LE.XMIN-RASC*ANTOAU
     *   .AND. CX.GE.XMIN-RASC*ANTOAU*3.0D+00)
     *   .OR.  (CY.GE.YMAX+RASC*ANTOAU
     *   .AND. CY.LE.YMAX+RASC*ANTOAU*3.0D+00)
     *   .OR.  (CY.LE.YMIN-RASC*ANTOAU
     *   .AND. CY.GE.YMIN-RASC*ANTOAU*3.0D+00)
     *   .OR.  (CZ.GE.ZMAX+RASC*ANTOAU
     *   .AND. CZ.LE.ZMAX+RASC*ANTOAU*3.0D+00)
     *   .OR.  (CZ.LE.ZMIN-RASC*ANTOAU
     *   .AND. CZ.GE.ZMIN-RASC*ANTOAU*3.0D+00)
      IF(.NOT.USEIMG) GO TO 480
C
      XMU  = ASCDIP(1,IC)
      YMU  = ASCDIP(2,IC)
      ZMU  = ASCDIP(3,IC)
      XX   = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
      IF(NROOTS.LE.3) CALL RT123
      IF(NROOTS.EQ.4) CALL ROOT4
      IF(NROOTS.EQ.5) CALL ROOT5
      IF(NROOTS.EQ.6) CALL ROOT6
      IF(NROOTS.GE.7) THEN
         IF (MASWRK) WRITE(IW,9008)
         CALL ABRT
      END IF
      MM = 0
      DO 401 K = 1,NROOTS
      UU = AA*U(K)
      WW = W(K)
      WW=WW*UU
      TT = ONE/(AA+UU)
      T  =  SQRT(TT)
      X0 = (AAX+UU*CX)*TT
      Y0 = (AAY+UU*CY)*TT
      Z0 = (AAZ+UU*CZ)*TT
      IN = -4+MM
      DO 400 I = 1,LIT
      IN = IN+4
      NI = I
      DO 400 J = 1,LJT
      JN = IN+J
      NJ = J
      CALL STVINT
      XIN(JN   ) = XINT
      YIN(JN   ) = YINT
      ZIN(JN   ) = ZINT*WW
      CALL POLXYZ
      XIN(JN+125) = XINT
      YIN(JN+125) = YINT
      ZIN(JN+125) = ZINT*WW
  400 CONTINUE
  401 MM = MM+16
      DO 403 I = 1,IJ
      NX    = IJX(I)
      NY    = IJY(I)
      NZ    = IJZ(I)
      DUMX = ZERO
      DUMY = ZERO
      DUMZ = ZERO
      MM    = 0
      DO 402 K = 1,NROOTS
      DUMX= DUMX+XIN(NX+MM+125)*YIN(NY+MM    )*ZIN(NZ+MM    )
      DUMY= DUMY+XIN(NX+MM    )*YIN(NY+MM+125)*ZIN(NZ+MM    )
      DUMZ= DUMZ+XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM+125)
  402 MM  = MM+16
      DUM = FIJ(I)
      CHDINT(I)=CHDINT(I)-DUM*(DUMX*XMU+DUMY*YMU+DUMZ*ZMU)
 403  CONTINUE
 480  CONTINUE
C
C-------- END OF ASCDIP CENTERS LOOP
C
 500  CONTINUE
 520  CONTINUE
C
C     ----- END OF *PRIMITIVE* LOOPS -----
C     ----- SET QQ MATRIX
C
      MAX=MAXJ
      NN=0
      DO 550 I=MINI,MAXI
      LI=LOCI+I
      IN = (LI*(LI-1))/2
      IF(IANDJ) MAX=I
      DO 550 J=MINJ,MAX
      LJ=LOCJ+J
      JN=LJ+IN
      NN=NN+1
      QQ(JN)=QQ(JN) + CHDINT(NN)
  550 CONTINUE
C
C     ----- END PARALLEL
C
  570    CONTINUE
C
C     ----- END SHELLS -----
C
  580 CONTINUE
  600 CONTINUE
C
C     ----- END OF *SHELL* LOOPS -----
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) CALL DDI_GSUMF(921,QQ,NUM2)
C
C     ----- PRINTING SECTION ----
C
      IF(MASWRK.AND.OUT) THEN
         WRITE(6,9380)
         CALL PRTRIL(QQ,NUM)
         WRITE(6,9128)
         CALL TIMIT(1)
      END IF
      RETURN
C
 9008 FORMAT(/' NUMBER OF POLYNOMIAL ROOTS NEEDED (NROOTS) IS GREATER',
     *        ' THAN 6 IN EFDINT.  CALL A PROGRAMMER/QUANTUM CHEMIST.')
 9128 FORMAT(/,' ...... END OF ASCDIP INTEGRALS....')
 9380 FORMAT(10X,14(1H-),/,10X,14H   EFD  MATRIX,' (CHARGE DIPOLE)',
     * /,10X,14(1H-))
      END
C*MODULE PCMIEF  *DECK QADINT
      SUBROUTINE QADINT(QQ,ASCCRD,ASCQAD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION QQ(*),ASCCRD(3,*),ASCQAD(6,*)
      COMMON /FMCOM / XX(1)
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
C
      CALL VALFM(LOADFM)
      NEED=6*MXSP
      CALL GETFM(NEED)
      CALL QADIT(QQ,XX(LOADFM+1),ASCCRD,ASCQAD)
      CALL RETFM(NEED)
C
      RETURN
      END
C*MODULE PCMIEF  *DECK QADIT
      SUBROUTINE QADIT(QQ,QUAD,ASCCRD,ASCQAD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL IANDJ,NORM,DOUBLE
      LOGICAL OUT
      LOGICAL GOPARR,MASWRK,DSKWRK,USEIMG
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
      PARAMETER (TOANGS=0.52917724924D+00, ANTOAU=1.0D+00/TOANGS)
C
      COMMON /CSSTV / CX,CY,CZ
      COMMON /ABIIMG/ XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,
     *                NAMIMG(3*MXATM),
     *                MASIMG(3*MXATM),ZANIMG(3*MXATM),
     *                CRDIMG(3,3*MXATM),CHGIMG(2,3*MXATM),
     *                DIPIMG(3,3*MXATM),QADIMG(6,3*MXATM),NUMIMG
      COMMON /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 /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /STV   / XINT,YINT,ZINT,T,X0,Y0,Z0,XI,YI,ZI,XJ,YJ,ZJ,NI,NJ
C
      DIMENSION QUAD(6,*),CHQINT(100),QQ(1),ASCCRD(3,*),ASCQAD(6,*)
      DIMENSION DIJ(100),XIN(432),YIN(432),ZIN(432)
      DIMENSION FIJ(100)
      DIMENSION IX(20),IY(20),IZ(20),JX(20),JY(20),JZ(20)
      DIMENSION IJX(100),IJY(100),IJZ(100),W2(5),W4(5)
C
      DATA PI212 /1.1283791670955D+00/
      DATA SQRT3 /1.73205080756888D+00/
      DATA SQRT5 /2.23606797749979D+00/
      DATA JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 1/
      DATA JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1/
      DATA IX / 1, 5, 1, 1, 9, 1, 1, 5, 5, 1,
     1         13, 1, 1, 9, 9, 5, 1, 5, 1, 5/
      DATA IY / 1, 1, 5, 1, 1, 9, 1, 5, 1, 5,
     1          1,13, 1, 5, 1, 9, 9, 1, 5, 5/
      DATA IZ / 1, 1, 1, 5, 1, 1, 9, 1, 5, 5,
     1          1, 1,13, 1, 5, 1, 5, 9, 9, 5/
      DATA ZERO,ONE,TWO,FOUR/0.0D+00,1.0D+00,2.0D+00,4.0D+00/
      DATA PT5,ONEPT5,THREE /0.5D+00,1.5D+00,3.0D+00/
      DATA RLN10 /2.30258D+00/
      DATA DBUGME/8HINT1    /
C
C    CALCULATE THE CHARGE-QUADRUPOLE CONTRIBUTION TO THE ONE ELECTRON
C    INTEGRALS FOR THE PCM (REGION WISE)
C
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
      NEFQ = NREG
      OUT = NPRINT.EQ.3  .OR.  EXETYP.EQ.DBUGME
      TOL = RLN10*ITOL
      NORM = NORMF .NE. 1 .OR. NORMP .NE. 1
      NUM2=(NUM*(NUM+1))/2
      DO 10 I=1,NUM2
      QQ(I) = ZERO
 10   CONTINUE
C
C  REPLACE SECOND MOMENTS IN ARRAY SM WITH ELECTRIC QUADRUPOLE TENSORS
C
      DO 20 I=1,NEFQ
      XX =  ASCQAD(1,I)
      YY =  ASCQAD(2,I)
      ZZ =  ASCQAD(3,I)
      XY =  ASCQAD(4,I)
      XZ =  ASCQAD(5,I)
      YZ =  ASCQAD(6,I)
      DUM = XX + YY + ZZ
      QUAD(1,I) = (THREE * XX - DUM) * PT5
      QUAD(2,I) = (THREE * YY - DUM) * PT5
      QUAD(3,I) = (THREE * ZZ - DUM) * PT5
      QUAD(4,I) = ONEPT5 * XY
      QUAD(5,I) = ONEPT5 * XZ
      QUAD(6,I) = ONEPT5 * YZ
 20   CONTINUE
C
C     ----- ISHELL
C
      DO 600 II = 1,NSHELL
      IAT= KATOM(II)
      XI = C(1,IAT)
      YI = C(2,IAT)
      ZI = C(3,IAT)
      I1 = KSTART(II)
      I2 = I1+KNG(II)-1
      LIT = KTYPE(II)
      MINI = KMIN(II)
      MAXI = KMAX(II)
      LOCI = KLOC(II)-MINI
C
C     ----- JSHELL
C
      DO 580 JJ = 1,II
C
C     ----- GO PARALLEL! -----
C
            IF (GOPARR) THEN
               IPCOUNT = IPCOUNT + 1
               IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 570
            END IF
      JAT= KATOM(JJ)
      XJ = C(1,JAT)
      YJ = C(2,JAT)
      ZJ = C(3,JAT)
      J1 = KSTART(JJ)
      J2 = J1+KNG(JJ)-1
      LJT = KTYPE(JJ)
      MINJ = KMIN(JJ)
      MAXJ = KMAX(JJ)
      LOCJ = KLOC(JJ)-MINJ
      NROOTS = (LIT+LJT+2-2)/2 + 1
      RR = (XI-XJ)**2+(YI-YJ)**2+(ZI-ZJ)**2
      IANDJ = II .EQ. JJ
C
C     ----- PREPARE INDICES FOR PAIRS OF (I,J) FUNCTIONS
C
      IJ = 0
      MAX = MAXJ
      DO 30 I = MINI,MAXI
      NX = IX(I)
      NY = IY(I)
      NZ = IZ(I)
      IF (IANDJ) MAX = I
      DO 30 J = MINJ,MAX
      IJ = IJ+1
      IJX(IJ) = NX+JX(J)
      IJY(IJ) = NY+JY(J)
      IJZ(IJ) = NZ+JZ(J)
  30  CONTINUE
      DO 60 I=1,IJ
  60  CHQINT(I) = ZERO
C
C     ----- I PRIMITIVE
C
      JGMAX = J2
      DO 520 IG = I1,I2
      AI = EX(IG)
      ARRI = AI*RR
      AXI = AI*XI
      AYI = AI*YI
      AZI = AI*ZI
      CSI = CS(IG)
      CPI = CP(IG)
      CDI = CD(IG)
      CFI=CF(IG)
C
C
C     ----- J PRIMTIVE
C
      IF (IANDJ) JGMAX = IG
      DO 500 JG = J1,JGMAX
      AJ = EX(JG)
      AA = AI+AJ
      AA1 = ONE/AA
      DUM = AJ*ARRI*AA1
      IF (DUM .GT. TOL) GO TO 500
      FAC =  EXP(-DUM)
      CSJ = CS(JG)
      CPJ = CP(JG)
      CDJ = CD(JG)
      CFJ=CF(JG)
      AX = (AXI+AJ*XJ)*AA1
      AY = (AYI+AJ*YJ)*AA1
      AZ = (AZI+AJ*ZJ)*AA1
C
C     ----- DENSITY FACTOR
C
      DOUBLE=IANDJ.AND.IG.NE.JG
      MAX = MAXJ
      NN = 0
      DO 170 I=MINI,MAXI
      GO TO (70,80,110,110,90,110,110,100,110,110,
     1       102,110,110,104,110,110,110,110,110,106),I
   70 DUM1=CSI*FAC
      GO TO 110
   80 DUM1=CPI*FAC
      GO TO 110
   90 DUM1=CDI*FAC
      GO TO 110
  100 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 110
  102 DUM1=CFI*FAC
      GO TO 110
  104 IF(NORM) DUM1=DUM1*SQRT5
      GO TO 110
  106 IF(NORM) DUM1=DUM1*SQRT3
  110 IF(IANDJ) MAX=I
      DO 170 J=MINJ,MAX
      GO TO (125,130,160,160,140,160,160,150,160,160,
     1       152,160,160,154,160,160,160,160,160,156),J
  125 DUM2=DUM1*CSJ
      IF(.NOT.DOUBLE) GO TO 160
      IF(I.GT.1) GO TO 126
      DUM2=DUM2+DUM2
      GO TO 160
  126 DUM2=DUM2+CSI*CPJ*FAC
      GO TO 160
  130 DUM2=DUM1*CPJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  140 DUM2=DUM1*CDJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  150 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 160
  152 DUM2=DUM1*CFJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 160
  154 IF(NORM) DUM2=DUM2*SQRT5
      GO TO 160
  156 IF(NORM) DUM2=DUM2*SQRT3
  160 NN=NN+1
  170 DIJ(NN)=DUM2
C
C     ..... HELLMANN-FEYNMAN TERM .....
C
      DUM = PI212*AA1/THREE
      DO 380 I = 1,IJ
      FIJ(I) = DIJ(I)*DUM
  380 CONTINUE
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
      DO 480 IC = 1,NEFQ
      CX   = ASCCRD(1,IC)
      CY   = ASCCRD(2,IC)
      CZ   = ASCCRD(3,IC)
      USEIMG = .FALSE.
      USEIMG = (CX.GE.XMAX+RASC*ANTOAU
     *   .AND. CX.LE.XMAX+RASC*ANTOAU*2.0D+00)
     *   .OR.  (CX.LE.XMIN-RASC*ANTOAU
     *   .AND. CX.GE.XMIN-RASC*ANTOAU*2.0D+00)
     *   .OR.  (CY.GE.YMAX+RASC*ANTOAU
     *   .AND. CY.LE.YMAX+RASC*ANTOAU*2.0D+00)
     *   .OR.  (CY.LE.YMIN-RASC*ANTOAU
     *   .AND. CY.GE.YMIN-RASC*ANTOAU*2.0D+00)
     *   .OR.  (CZ.GE.ZMAX+RASC*ANTOAU
     *   .AND. CZ.LE.ZMAX+RASC*ANTOAU*2.0D+00)
     *   .OR.  (CZ.LE.ZMIN-RASC*ANTOAU
     *   .AND. CZ.GE.ZMIN-RASC*ANTOAU*2.0D+00)
      IF(.NOT.USEIMG) GO TO 480
C
      XX   = AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
      IF(NROOTS .LE. 3) CALL RT123
      IF(NROOTS .EQ. 4) CALL ROOT4
      IF(NROOTS.EQ.5) CALL ROOT5
      IF(NROOTS.EQ.6) CALL ROOT6
      IF(NROOTS.GE.7) THEN
       IF (MASWRK) WRITE(IW,9008)
       CALL ABRT
      END IF
      MM = 0
      DO 401 K = 1,NROOTS
      U2 = AA*U(K)
      U4 = U2*U2
      WW = W(K)
      W2(K) = TWO *WW*U2
      W4(K) = FOUR*WW*U4
      TT = ONE/(AA+U2)
      T  =  SQRT(TT)
      X0 = (AAX+U2*CX)*TT
      Y0 = (AAY+U2*CY)*TT
      Z0 = (AAZ+U2*CZ)*TT
      IN = -4+MM
      DO 400 I = 1,LIT
      IN = IN+4
      NI = I
      DO 400 J = 1,LJT
      JN = IN+J
      NJ = J
      CALL STVINT
      XIN(JN   ) = XINT
      YIN(JN   ) = YINT
      ZIN(JN   ) = ZINT
      CALL POLXYZ
      XIN(JN+144) = XINT
      YIN(JN+144) = YINT
      ZIN(JN+144) = ZINT
      CALL EFQXYZ
      XIN(JN+288) = XINT
      YIN(JN+288) = YINT
      ZIN(JN+288) = ZINT
  400 CONTINUE
  401 MM = MM+16
      DO 403 I = 1,IJ
      NX    = IJX(I)
      NY    = IJY(I)
      NZ    = IJZ(I)
      DUMXX = ZERO
      DUMYY = ZERO
      DUMZZ = ZERO
      DUMXY = ZERO
      DUMXZ = ZERO
      DUMYZ = ZERO
      MM    = 0
      DO 402 K = 1,NROOTS
      DUMXX= DUMXX-XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM    )*W2(K)
     1            +XIN(NX+MM+288)*YIN(NY+MM    )*ZIN(NZ+MM    )*W4(K)
      DUMYY= DUMYY-XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM    )*W2(K)
     1            +XIN(NX+MM    )*YIN(NY+MM+288)*ZIN(NZ+MM    )*W4(K)
      DUMZZ= DUMZZ-XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM    )*W2(K)
     1            +XIN(NX+MM    )*YIN(NY+MM    )*ZIN(NZ+MM+288)*W4(K)
      DUMXY= DUMXY+XIN(NX+MM+144)*YIN(NY+MM+144)*ZIN(NZ+MM    )*W4(K)
      DUMXZ= DUMXZ+XIN(NX+MM+144)*YIN(NY+MM    )*ZIN(NZ+MM+144)*W4(K)
      DUMYZ= DUMYZ+XIN(NX+MM    )*YIN(NY+MM+144)*ZIN(NZ+MM+144)*W4(K)
  402 MM  = MM+16
      DUM = FIJ(I)
      CHQINT(I) = CHQINT(I) -
     1 DUM * ( DUMXX * QUAD(1,IC)
     2       + DUMYY * QUAD(2,IC)
     3       + DUMZZ * QUAD(3,IC)
     4 +TWO* ( DUMXY * QUAD(4,IC)
     5       + DUMXZ * QUAD(5,IC)
     6       + DUMYZ * QUAD(6,IC) ) )
 403  CONTINUE
 480  CONTINUE
C
C-------- END OF ASCQAD CENTERS LOOP
C
 500  CONTINUE
 520  CONTINUE
C
C     ----- END OF *PRIMITIVE* LOOPS -----
C     ----- SET QQ MATRIX
C
      MAX=MAXJ
      NN=0
      DO 550 I=MINI,MAXI
      LI=LOCI+I
      IN= (LI*(LI-1))/2
      IF(IANDJ) MAX=I
      DO 550 J=MINJ,MAX
      LJ=LOCJ+J
      JN=LJ+IN
      NN=NN+1
      QQ(JN)=QQ(JN) + CHQINT(NN)
  550 CONTINUE
C
C     ----- END PARALLEL
C
  570    CONTINUE
C
C     ----- END SHELLS -----
C
  580 CONTINUE
  600 CONTINUE
C
C     ----- END OF *SHELL* LOOPS -----
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) CALL DDI_GSUMF(923,QQ,NUM2)
C
C     ----- PRINTING SECTION ----
C
      IF(MASWRK.AND.OUT) THEN
         WRITE(IW,9380)
         CALL PRTRIL(QQ,NUM)
         WRITE(IW,9128)
         CALL TIMIT(1)
      END IF
      RETURN
C
 9008 FORMAT(/' NUMBER OF POLYNOMIAL ROOTS NEEDED (NROOTS) IS GREATER',
     *        ' THAN 6 IN EFQINT.  CALL A PROGRAMMER/QUANTUM CHEMIST.')
 9128 FORMAT(/,' ...... END OF ASCQAD INTEGRALS....')
 9380 FORMAT(10X,14(1H-),/,10X,14H   EFQ  MATRIX,' (CHARGE QUARUPOLE)',
     * /,10X,14(1H-))
      END
C*MODULE PCMIEF  *DECK IMASCX
      SUBROUTINE IMASCX(TMP,XCTS,YCTS,ZCTS,QSE,XE,YE,ZE,ISPHE,LIST,
     *                  ASCCRD,ASCCHG,ASCDIP,ASCQAD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
C
      DIMENSION TMP(NTS,3),XCTS(*),YCTS(*),ZCTS(*),QSE(*),XE(*),YE(*),
     *          ZE(*),ISPHE(*),LIST(*),ASCCRD(3,*),ASCCHG(*),
     *          ASCDIP(3,*),ASCQAD(6,*)
C
C
C     REGION MULTIPOLES UP TO QUADRUPOLES
C
C     -- ZERO OUT ASCIMG
C
      DO I=1,MXSP
            ASCCHG(I)   = 0.0D+00
         DO J=1,3
            ASCCRD(J,I) = 0.0D+00
            ASCDIP(J,I) = 0.0D+00
         ENDDO
         DO J=1,6
            ASCQAD(J,I) = 0.0D+00
         ENDDO
      ENDDO
C
C     -- COORDINATES
C
      DO I=1,NREG
        ASCCRD(1,I) = XE(ISPHE(LIST(I)))
        ASCCRD(2,I) = YE(ISPHE(LIST(I)))
        ASCCRD(3,I) = ZE(ISPHE(LIST(I)))
      ENDDO
C
C
C     -- MONOPOLES
C
      DO I=1,NREG
         DO J=LIST(I),LIST(I+1)-1
            ASCCHG(I)=ASCCHG(I)+QSE(J)
         ENDDO
      ENDDO
C
C
C     DISTANCES RELATIVES TO THE REGION CENTER
C
      DO ITS=1,NTS
         L=ISPHE(ITS)
         TMP(ITS,1)=XCTS(ITS)-XE(L)
         TMP(ITS,2)=YCTS(ITS)-YE(L)
         TMP(ITS,3)=ZCTS(ITS)-ZE(L)
      ENDDO
C
C     -- DIPOLES
C
      DO I=1,NREG
         DO J=LIST(I),LIST(I+1)-1
            ASCDIP(1,I)=ASCDIP(1,I)+QSE(J)*TMP(J,1) !DX
            ASCDIP(2,I)=ASCDIP(2,I)+QSE(J)*TMP(J,2) !DY
            ASCDIP(3,I)=ASCDIP(3,I)+QSE(J)*TMP(J,3) !DZ
         ENDDO
      ENDDO
C
C
C     -- QUADRUPOLES
C          SINCE 'ASCQAD' IS CALLED BY 'QADINT', WHICH FOLLOWS
C          STONE'S MULTIPOLAR EXPANSION FORMALISM, HERE
C          WE FOLLOW STONE'S FORMALISM TOO.
C          SEE: REGMUL, WHICH FOLLOWS THE MOST POPULAR FORMALISM
C
      DO I=1,NREG
       DO J=LIST(I),LIST(I+1)-1
        ASCQAD(1,I)=ASCQAD(1,I)+QSE(J)*TMP(J,1)*TMP(J,1)  !QXX
        ASCQAD(2,I)=ASCQAD(2,I)+QSE(J)*TMP(J,2)*TMP(J,2)  !QYY
        ASCQAD(3,I)=ASCQAD(3,I)+QSE(J)*TMP(J,3)*TMP(J,3)  !QZZ
        ASCQAD(4,I)=ASCQAD(4,I)+QSE(J)*TMP(J,1)*TMP(J,2)  !QXY
        ASCQAD(5,I)=ASCQAD(5,I)+QSE(J)*TMP(J,1)*TMP(J,3)  !QXZ
        ASCQAD(6,I)=ASCQAD(6,I)+QSE(J)*TMP(J,2)*TMP(J,3)  !QYZ
       ENDDO
      ENDDO
C
      RETURN
      END
