C 20 AUG 07 - MWS - ADD STORAGE FOR DIPIT CALL
C 24 MAR 07 - MWS - PAD THE FRGINF,EFPPAR COMMON BLOCKS
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 17 JAN 06 - HL  - CHANGES FOR EFP INDUCED DIPOLES
C 14 NOV 05 - HL  - IMPLEMENT PCM/EFP GRADIENT
C 14 NOV 05 - DGF - DYNAM. ALLOCATE PCM ARRAYS; VARIOUS PCM IMPROVEMENTS
C 19 SEP 05 - HMN,IA - CHANGES FOR MOLECULAR DYNAMICS
C 15 JUL 05 - MWS - FIX DYNAMIC STORAGE FOR DIPIT
C  5 JUL 05 - HL  - ENUADD,CHGIND: FIX QM/EFP/PCM MATR. INVERSION
C  5 JUL 05 - MWS - ALLOCATE DYNAMIC MEMORY FOR DIPIT
C  1 JUN 05 - MWS - CHGIND: 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 FRGINF, EFPPAR
C  7 APR 04 - HL  - IEF-PCM/EFP FOR ALL TYPES OF WAVEFUNCTION
C  7 AUG 02 - CP,HL - ITERATIVE ISOTROPIC IEF-PCM/EFP
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C  8 OCT 01 - HL  - PARALLELIZE PCM
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - MWS - IEFFRG: DAF RECORD USED TO PASS IFR VALUES
C 11 OCT 00 - PB  - STATIC PART OF EFP CONTRIBUTION TO PCM
C
C*MODULE EFPCM  *DECK ENUADD
      SUBROUTINE ENUADD(DMATM1,VEC,VEC_I,VEC_TEST,
     *                  L2,XCTS,YCTS,ZCTS,AS,QSN,Q_FS,VEC_MUL)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      CHARACTER*8 FRGNME
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXATM=2000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION VEC(L2),VEC_I(L2),VEC_TEST(L2),
     *          DMATM1(NTS,NTS),XCTS(*),YCTS(*),ZCTS(*),
     *          AS(*),QSN(*),Q_FS(*),VEC_MUL(*)
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 /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 /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      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 /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,IKREP,IEF,IP_F,NFMOPCM
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO = 0.00D+00, ONE = 1.00D+00)
      PARAMETER (HALF = 0.5D+00)
      DATA RMC/8HMCSCF   /
      DATA GVB/8HGVB     /
C
C
      SCALE = (EPS-ONE)/EPS
      NUM2=(NUM*NUM+NUM)/2
      CALL VCLR(VEC_I,1,NUM2)
      CALL VCLR(Q_FS,1,NTS)
      CALL VCLR(VEC_MUL,1,NTS)
      P_NF=ZERO
      P_FF=ZERO
      P_FN=ZERO
C
      DO 300 ITS=1,NTS
         XI=XCTS(ITS)
         YI=YCTS(ITS)
         ZI=ZCTS(ITS)
C
         DO 500 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
 500    CONTINUE
 300  CONTINUE
C
      DO ITS=1,NTS
         DO JTS=1,NTS
            Q_FS(ITS)=Q_FS(ITS)+DMATM1(ITS,JTS)*VEC_MUL(JTS)
         ENDDO
         Q_FS(ITS)=Q_FS(ITS)*AS(ITS)
         IF(IEF.EQ.10) Q_FS(ITS)=Q_FS(ITS)*SCALE
      ENDDO
C
C     INTERACTIONS
C
      DO ITS=1,NTS
         CALL INTMEP(VEC,XCTS(ITS),YCTS(ITS),ZCTS(ITS))
         DO IBAS=1,NUM2
            VEC_I(IBAS)=VEC_I(IBAS)-VEC(IBAS)*Q_FS(ITS)
         ENDDO
C
         P_FF=P_FF+VEC_MUL(ITS)*Q_FS(ITS)
         P_FN=P_FN+VEC_MUL(ITS)*QSN(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)
            P_NF = P_NF + Q_FS(ITS)*ZAN(JATOM)/R
         ENDDO
      ENDDO
C
      UNZ = UNZ + (P_FF + P_FN + P_NF)*HALF
C
      IF(SCFTYP.EQ.RMC.OR.SCFTYP.EQ.GVB) THEN
         CALL DAREAD(IDAF,IODA,VEC_TEST,NUM2,87,0)
      ELSE
         CALL DAREAD(IDAF,IODA,VEC_TEST,NUM2,11,0)
      END IF
C
      DO IBAS = 1,NUM2
         VEC_TEST(IBAS) = VEC_TEST(IBAS) + VEC_I(IBAS)
      ENDDO
C
      IF(SCFTYP.EQ.RMC.OR.SCFTYP.EQ.GVB) THEN
         CALL DAWRIT(IDAF,IODA,VEC_TEST,NUM2,87,0)
      ELSE
         CALL DAWRIT(IDAF,IODA,VEC_TEST,NUM2,11,0)
      END IF
C
      RETURN
      END
C
C*MODULE EFPCM  *DECK CHGIND
      SUBROUTINE CHGIND(H1,FA,FB,DENTOT,DIP,ODIP,DMATM1,
     *                  XEFI,YEFI,ZEFI,EFLD,EFADD,ABFLD,VECMP,
     *                  VEC_2,VEC_S,L2,L1,
     *                  MCITER,XCTS,YCTS,ZCTS,AS,QSN,QSE,Q_FS,
     *                  Q_IND,AIND_PROJ,FIND_PROJ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXRT=100, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION H1(L2),FA(L2),FB(L2),DENTOT(L2),DIP(L2),ODIP(L2),
     *          DMATM1(NTS,NTS),XEFI(L2),YEFI(L2),ZEFI(L2),
     *          EFLD(3,NPTTPT),EFADD(3,NPTTPT),ABFLD(3,NPTTPT),
     *          VECMP(L2),VEC_2(L2),VEC_S(L2),
     *          XCTS(*),YCTS(*),ZCTS(*),AS(*),QSN(*),
     *          QSE(*),Q_FS(*),Q_IND(*),AIND_PROJ(*),FIND_PROJ(*)
      DIMENSION GSUMTMP(4)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      LOGICAL GOPARR,MASWRK,DSKWRK,SVDSKW,
     *        DIRSCF,FDIFF,DIRDIF,MINMEM
C
      COMMON /CONV  / ACURCY,ENUCR,AJUNK(3),DIFF,ITER,ICALCP,ICBET
      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 /ENRGYS/ ENUCRMC,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / XXX(1)
      COMMON /FRFD  / PPEFX(MXFGPT), PPEFY(MXFGPT), PPEFZ(MXFGPT)
      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 /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      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 /OPTSCF/ DIRSCF,FDIFF
      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 /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,LFIPRJ
     *                ,LPEL
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /PNUC  / PPNX(MXFGPT),PPNY(MXFGPT),PPNZ(MXFGPT),PNZRF(3)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ONE=1.0D+00, GLIM=1.0D-03)
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, DEBYE=2.541766D+00)
C
      DATA RHF,GVB/8HRHF     ,8HGVB     /
      DATA RMC/8HMCSCF   /
C
      SAVE EN0
C
C     THIS ROUTINE INDUCES THE COMBINED EFP AND PCM FIELD, AND THEN
C     MODIFIES THE 1E- PART OF THE FOCK MATRIX TO REFLECT THIS FIELD.
C
C     NOTE TO PROGRAMMERS:  THIS ROUTINE IS DESCENDED FROM -EFPCMP-
C     AND SHOULD BE KEPT AS SIMILAR AS POSSIBLE TO ITS PARENT!
C
C     CALCULATE ITERATIVELY THE INDUCED DIPOLES ON THE FRAGMENTS
C     UNTIL SELF-CONSISTENCY IS REACHED, ADD RESULTING FIELD DURING
C     SAME SCF ITERATION. ITERATION OF INDUCED DIPOLES ALSO COUPLED
C     TO SCF.
C
C     NOTE: GVB ITERATIONS START AT 1, RHF,UHF,ROHF AT 0.
C     EN0 IS USUAL NUCLEAR REPULSION, NO SOLVENT JUNK ADDED.
C
      SCALE = (EPS-ONE)/EPS
      CALL VCLR(VEC_2,1,L2)
      CALL VCLR(AIND_PROJ,1,NTS)
      CALL VCLR(FIND_PROJ,1,NTS)
      P_I = ZERO
      P_NUCC = ZERO
C
      IFIRST = 0
      IF(SCFTYP.NE.GVB) IFIRST=1
      IF(SCFTYP.EQ.RMC) THEN
         IFIRST=2
         ITER=MCITER
      END IF
      DIRDIF = DIRSCF  .AND.  FDIFF  .AND.  SCFTYP.NE.GVB
C
C     -- FIELDS OF NUCLEI AND EFP MULTIPOLE POINTS --
C
      IF(ITER.EQ.IFIRST) THEN
         EN0=ENUCR
         IF(SCFTYP.EQ.RMC) EN0=ENUCRMC
         CALL EFNUCP(NAT,ZAN,C,ENUCP)
C
         CALL VCLR(Q_IND,1,NTS)
         CALL FRGFLD
C
         CALL VCLR(EFADD,1,3*NPTTPT)
      ELSE
         ENUCR=EN0
      END IF
C
C
C     -- ELECTRONIC FIELDS AT POL-POINTS
C        ELECTRIC FIELD INTEGRAL FILE -NFT25- IS A DISTRIBUTED FILE.
C        ALL NODES MUST DO STATIC LOAD BALANCE PROCESSING OF THIS.
C
      IF(GOPARR) THEN
         CALL VCLR(ABFLD,1,3*NPTTPT)
         CALL VCLR(EFLD ,1,3*NPTTPT)
      END IF
      NFT25 = 25
      SVDSKW = DSKWRK
      DSKWRK = .TRUE.
      CALL SEQREW(NFT25)
      IPCOUNT = ME-1
      DO 10 LEFP = 1, NPTTPT
         IF(GOPARR) THEN
            IPCOUNT = IPCOUNT + 1
            IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 10
         END IF
C
         CALL SQREAD(NFT25,XEFI,L2)
         CALL SQREAD(NFT25,YEFI,L2)
         CALL SQREAD(NFT25,ZEFI,L2)
         ABFLD(1,LEFP)=TRACEP(DENTOT,XEFI,L1)
         ABFLD(2,LEFP)=TRACEP(DENTOT,YEFI,L1)
         ABFLD(3,LEFP)=TRACEP(DENTOT,ZEFI,L1)
C        -- SOLVENT APPARENT CHARGE FIELD AT THE POINT ----
         PQ_X=ZERO
         PQ_Y=ZERO
         PQ_Z=ZERO
         DO ITS = 1,NTS
            QT = QSN(ITS) + QSE(ITS) + Q_FS(ITS) + Q_IND(ITS)
            RR = (EFP(1,LEFP)-XCTS(ITS))**2 +
     *            (EFP(2,LEFP)-YCTS(ITS))**2 +
     *            (EFP(3,LEFP)-ZCTS(ITS))**2
            IF(RR.GE.GLIM) THEN
               RR3 = RR * SQRT(RR)
               PQ_X = PQ_X + QT *(EFP(1,LEFP)-XCTS(ITS))/RR3
               PQ_Y = PQ_Y + QT *(EFP(2,LEFP)-YCTS(ITS))/RR3
               PQ_Z = PQ_Z + QT *(EFP(3,LEFP)-ZCTS(ITS))/RR3
            END IF
         ENDDO
C        -- TOTAL FIELD --
         EFLD(1,LEFP)=ABFLD(1,LEFP)+PPNX(LEFP)+PQ_X+PPEFX(LEFP)
         EFLD(2,LEFP)=ABFLD(2,LEFP)+PPNY(LEFP)+PQ_Y+PPEFY(LEFP)
         EFLD(3,LEFP)=ABFLD(3,LEFP)+PPNZ(LEFP)+PQ_Z+PPEFZ(LEFP)
 10   CONTINUE
      CALL SEQREW(NFT25)
      IF(GOPARR) THEN
         CALL DDI_GSUMF(667,EFLD, 3*NPTTPT)
         CALL DDI_GSUMF(668,ABFLD,3*NPTTPT)
      END IF
C
C     ALWAYS CALL TO DIPIT !
C
      CALL VALFM(LOADFM)
      LDIPNEW  = LOADFM   + 1
      LDIPNWD  = LDIPNEW  + NFRG
      LDINDC   = LDIPNWD  + NFRG
      LDINDDC  = LDINDC   + 3*NPTTPT
      LDIND2C  = LDINDDC  + 3*NPTTPT
      LDINDD2C = LDIND2C  + 3*NPTTPT
      LDIND_W  = LDINDD2C + 3*NPTTPT
      LDIND_C  = LDIND_W  + 6*NPTTPT
      LDIND2_W = LDIND_C  + 6*NPTTPT
      LDIND2_C = LDIND2_W + 6*NPTTPT
      LEWLD1   = LDIND2_C + 6*NPTTPT
      LEWLD2   = LEWLD1   + 3*NPTTPT
      LAST     = LEWLD2   + 3*NPTTPT
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL DIPIT(EFLD,EFADD,XXX(LDINDC),XXX(LDINDDC),XXX(LDIND2C),
     *           XXX(LDINDD2C),XXX(LDIND_W),XXX(LDIND_C),
     *           XXX(LDIND2_W),XXX(LDIND2_C),NPTTPT,XXX(LDIPNEW),
     *           XXX(LDIPNWD),XXX(LEWLD1),XXX(LEWLD2),NFRG)
      CALL RETFM(NEED)
C
C     *********************************
C     **    INSERT PCM CODES HERE    **
C     **  ONLY FOR MATRIX INVERSION  **
C     *********************************
C
      IF (IPCMIT.EQ.0) THEN
C
C        --  POTENTIAL DUE TO INDUCED DIPOLES --
C
         DO ITS=1,NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
            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
               AIND_PROJ(ITS) = AIND_PROJ(ITS)
     *           + (DIND(1,LEFP)*X+DIND(2,LEFP)*Y+DIND(3,LEFP)*Z)/R3
            ENDDO
         ENDDO
C
C        -- PCM CHARGES DUE TO INDUCED DIPOLES
C
         DO ITS=1,NTS
            Q = ZERO
            DO JTS=1,NTS
               Q=Q+DMATM1(ITS,JTS)*AIND_PROJ(JTS)
            ENDDO
            Q_IND(ITS)=Q*AS(ITS)
            IF(IEF.EQ.10) Q_IND(ITS)=Q_IND(ITS)*SCALE
         ENDDO
C
C        -- COMPUTE THE INTERACTIONS --
C
         PB_POL = ZERO
         DO ITS=1,NTS
            XI=XCTS(ITS)
            YI=YCTS(ITS)
            ZI=ZCTS(ITS)
C
            CALL INTMEP(VECMP,XI,YI,ZI)
C
C           - INTERACTION BETWEEN QM ELECTRONIC POTENTIAL AND Q_IND
C
            DO IBAS=1,L2
               VEC_2(IBAS)=VEC_2(IBAS)-VECMP(IBAS)*Q_IND(ITS)
            ENDDO
C
            PB_POL=PB_POL+XXX(LPEL+ITS-1)*Q_IND(ITS)
C
C           - INTERACTION BETWEEN QM NUCLEI POTENTIAL AND Q_IND
C
            DO JATOM = 1, NAT
               R2 = (C(1,JATOM)-XI)**2 +
     *              (C(2,JATOM)-YI)**2 +
     *              (C(3,JATOM)-ZI)**2
               R = SQRT(R2)
               P_I= P_I + Q_IND(ITS)*ZAN(JATOM) / R
            ENDDO
C
C           - INTERACTION BETWEEN EFP MULTIPOLES AND Q_IND
C
            P_NUCC = P_NUCC + Q_IND(ITS)*XXX(LVECMUL+ITS-1)
C
         ENDDO
C
         PB = PB + PB_POL
         ENUCR = ENUCR + PT5*(P_I + P_NUCC)
C
C        -- ADD VEC_2 TO THE HAMILTONIAN
C
         CALL VADD(H1,1,VEC_2,1,H1,1,L2)
C
C        -- ADD DIFFERENCE IN PERTURBATION WHEN DIRSCF AND FDIFF -----
C           SINCE VEC_S IS NOT USED, WE USE IT TO SAVE VEC_2
C
         IF(SCFTYP .EQ. RHF) THEN
            IF(DIRSCF  .AND.  FDIFF) THEN
               IF(ITER.EQ.1) CALL VCLR(VEC_S,1,L2)
               CALL VSUB(VEC_S,1,VEC_2,1,VEC_S,1,L2)
               CALL VADD(FA,1,VEC_S,1,FA,1,L2)
               CALL DCOPY(L2,VEC_2,1,VEC_S,1)
            END IF
         ELSE
            IF(DIRSCF  .AND.  FDIFF .AND. SCFTYP .NE. GVB) THEN
               IF(ITER.EQ.1) CALL VCLR(VEC_S,1,L2)
               CALL VSUB(VEC_S,1,VEC_2,1,VEC_S,1,L2)
               CALL VADD(FA,1,VEC_S,1,FA,1,L2)
               CALL VADD(FB,1,VEC_S,1,FB,1,L2)
               CALL DCOPY(L2,VEC_2,1,VEC_S,1)
            END IF
         END IF
      END IF
C
C     *********************************
C     ***      END OF PCM CODES     ***
C     *********************************
C
C
C
      IF(GOPARR) THEN
         ENUCR = ENUCR/NPROC
         IF(.NOT.DIRDIF) CALL DSCAL(L2,ONE/NPROC,H1,1)
      END IF
      IPCOUNT = ME-1
C
      IF(DIRDIF) CALL VCLR(DIP,1,L2)
C
      TMUX=ZERO
      TMUY=ZERO
      TMUZ=ZERO
      DO 25 LEFP=1,NPTTPT
C
C     ----- GO PARALLEL! -----
C
         IF (GOPARR) THEN
            IPCOUNT = IPCOUNT + 1
            IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 25
         END IF
C
         PNX=PPNX(LEFP)
         PNY=PPNY(LEFP)
         PNZ=PPNZ(LEFP)
         PMUX = DIND(1,LEFP)
         PMUY = DIND(2,LEFP)
         PMUZ = DIND(3,LEFP)
         PMUXD = DINDD(1,LEFP)
         PMUYD = DINDD(2,LEFP)
         PMUZD = DINDD(3,LEFP)
         PEX=ABFLD(1,LEFP)
         PEY=ABFLD(2,LEFP)
         PEZ=ABFLD(3,LEFP)
C
C        --- CALCULATE ALPHA*FMU VALUES THAT ARE USED IN LOOPS BELOW ---
C
         EFX = EFPOL(1,LEFP)*EFADD(1,LEFP)
     *        +EFPOL(7,LEFP)*EFADD(2,LEFP)
     *        +EFPOL(8,LEFP)*EFADD(3,LEFP)
         EFY = EFPOL(4,LEFP)*EFADD(1,LEFP)
     *        +EFPOL(2,LEFP)*EFADD(2,LEFP)
     *        +EFPOL(9,LEFP)*EFADD(3,LEFP)
         EFZ = EFPOL(5,LEFP)*EFADD(1,LEFP)
     *        +EFPOL(6,LEFP)*EFADD(2,LEFP)
     *        +EFPOL(3,LEFP)*EFADD(3,LEFP)
C
         TOTMUX = PMUX + PMUXD
         TOTMUY = PMUY + PMUYD
         TOTMUZ = PMUZ + PMUZD
C
C        ----------- NUCLEAR ENERGY MODIFICATION -----------
C        ---- PLUS TOTAL ENERGY CORRECTION ACCORDING TO ----
C        +1/2*MU*[<SI|F|SI>-FEFP] - 0.5*ALPHA*FMU*[<SI|F|SI>]
C
         ENUCR = ENUCR + ((-PMUX*PNX - PMUY*PNY - PMUZ*PNZ)
     *                    + PMUXD*PEX + PMUYD*PEY + PMUZD*PEZ
     *                    - PMUX*PPEFX(LEFP)
     *                    - PMUY*PPEFY(LEFP)
     *                    - PMUZ*PPEFZ(LEFP)
     *                    - (EFX*PEX) - (EFY*PEY) - (EFZ*PEZ))*PT5
C
C
C        --- PREPARE POLARIZED ENVIRONMENT FOR NEXT SCF ITERATION -----
C
         CALL SQREAD(NFT25,XEFI,L2)
         CALL SQREAD(NFT25,YEFI,L2)
         CALL SQREAD(NFT25,ZEFI,L2)
C
C        --- CALCULATE PERTURBATION DUE TO INDUCED DIPOLES -----
         IF(DIRDIF) THEN
            DO K=1,L2
               DIP(K) = DIP(K)
     *         +(-TOTMUX*XEFI(K) - TOTMUY*YEFI(K) - TOTMUZ*ZEFI(K)
     *         + EFX*XEFI(K) + EFY*YEFI(K) + EFZ*ZEFI(K))*PT5
            ENDDO
         ELSE
            DO K=1,L2
               H1(K)=H1(K)
     *         +(-TOTMUX*XEFI(K) - TOTMUY*YEFI(K) - TOTMUZ*ZEFI(K)
     *         + EFX*XEFI(K) + EFY*YEFI(K) + EFZ*ZEFI(K))*PT5
            ENDDO
         END IF
C
         TMUX=TMUX+PMUX
         TMUY=TMUY+PMUY
         TMUZ=TMUZ+PMUZ
 25   CONTINUE
C
      CALL SEQREW(NFT25)
      DSKWRK = SVDSKW
C
      IF(GOPARR) THEN
         GSUMTMP(1) = TMUX
         GSUMTMP(2) = TMUY
         GSUMTMP(3) = TMUZ
         GSUMTMP(4) = ENUCR
         CALL DDI_GSUMF(685,GSUMTMP,4)
         TMUX = GSUMTMP(1)
         TMUY = GSUMTMP(2)
         TMUZ = GSUMTMP(3)
         ENUCR= GSUMTMP(4)
         IF(DIRDIF) THEN
            CALL DDI_GSUMF(686,DIP,L2)
         ELSE
            CALL DDI_GSUMF(687, H1,L2)
         END IF
      END IF
C
      ENUCRMC=ENUCR
C
      TMUXD=TMUX*DEBYE
      TMUYD=TMUY*DEBYE
      TMUZD=TMUZ*DEBYE
C
C     ----- ADD DIFFERENCE IN PERTURBATION WHEN DIRSCF AND FDIFF -----
C
      IF(SCFTYP .EQ. RHF) THEN
         IF(DIRSCF  .AND.  FDIFF) THEN
            IF(ITER.EQ.1) CALL VCLR(ODIP,1,L2)
            CALL VADD(H1,1,DIP,1,H1,1,L2)
            CALL VSUB(ODIP,1,DIP,1,ODIP,1,L2)
            CALL VADD(FA,1,ODIP,1,FA,1,L2)
            CALL DCOPY(L2,DIP,1,ODIP,1)
         END IF
      ELSE
         IF(DIRDIF) THEN
            IF(ITER .EQ. 1) CALL VCLR(ODIP,1,L2)
            CALL VADD(H1,1,DIP,1,H1,1,L2)
            CALL VSUB(ODIP,1,DIP,1,ODIP,1,L2)
            CALL VADD(FA,1,ODIP,1,FA,1,L2)
            CALL VADD(FB,1,ODIP,1,FB,1,L2)
            CALL DCOPY(L2,DIP,1,ODIP,1)
         END IF
      END IF
C
C     --- FOR MINMEM AND GVB H1+PERTURBATION MUST BE WRITTEN TO DAF ---
C
      CALL DAWRIT(IDAF,IODA,H1,L2,80,0)
      IF (MINMEM .OR. SCFTYP.EQ.GVB .OR. SCFTYP.EQ.RMC) THEN
         CALL DAWRIT(IDAF,IODA,H1,L2,11,0)
      END IF
      RETURN
      END
