C  4 MAR 08 - MWS - RAISE EFP PROCESSOR LIMIT FROM 64 TO 1024
C  7 DEC 07 - LVS - DIPIT: INCREASE NUMBER OF ITERATIONS ALLOWED
C 20 AUG 07 - LVS - PBC-EFP CLEANED
C 24 MAR 07 - LVS - POLARIZATION SCREENING ADDED
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 17 JAN 06 - HL  - SWITCHING FUNCTIONS FOR MD, EFPCMP: SIMPLIFICATION
C                   DIPIT: REMOVE METHOD=3, VARIOUS SIMPLIFICATIONS
C                   ADD FLDIND, VARIOUS CHANGES IN PARALLELIZATION
C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C 15 JUL 05 - MWS - DIPIT,EFPCMP: LIMIT VCLR RANGES, FIX DIPIT STORAGE
C 13 JUL 05 - SPW - EFPCMP: REMOVE WORK FROM INSIDE L2 LOOPS
C  5 JUL 05 - MWS - ALLOCATE DYNAMIC MEMORY FOR DIPIT
C 27 JUN 05 - HL,IA - PERIODIC BOUNDARY COND. FOR EFP
C  1 JUN 05 - MWS - DIPIT: FIX PARALLEL MP2/EFP1; DIPIT,EFPCMP: CLEAN UP
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 22 FEB 05 - MWS - DIPIT: USE INTEGER ZERO CALL
C 13 FEB 05 - HMN - EFPCMP,FRGFLD,POLINT,DIPIT: EFP PARALLELIZATION,
C                   ADD DUMMY EFP MD HOOKS INTO DIPIT AND FRGFLD
C 13 FEB 05 - JS  - ADD EFP MP2 CALLS TO DIPIT AND FRGFLD
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 JUL 03 - JMM - SUPPRESS PRINTING FOR MONTE CARLO JOBS
C  7 AUG 02 - IA  - FRGFLD: CHANGE CHARGE SCREENING, DIPIT: TIGHTEN CONV
C  7 AUG 02 - HL  - EFPCMP: ZERO OFF LARGE INDUCED DIPOLES
C 17 APR 02 - PND - EFPCMP: ONLY MASTER PRINTS A LINE
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 01 - HL  - TRIVIAL FORMAT CHANGE
C 20 FEB 01 - PND - DIPIT: PROVIDE ESCAPE HATCH IF NOT CONVERGENT
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C 16 FEB 00 - VK  - EFNUCP: SKIP OVER BUFFER NUCLEI WHEN DOING QM/MM
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 18 DEC 96 - JHJ - EFPCMP: WRITE 1-E HAILTONIAN TO DAF 80
C 30 OCT 96 - SPW - EFPCMP: CHANGES FOR MCSCF RUNS
C 29 SEP 96 - MWS - EFPCMP: REMOVE RESETTING OF NPRINT VALUE
C 18 SEP 96 - MWS - POLINT: REMOVE TIMING CALL
C 10 SEP 96 - MWS - INCLUDE THIS CODE IN DISTRIBUTION VERSION
C 13 JUN 96 - MWS - INTRODUCE CI TYPE VARIABLE, DELETE ZDPMLT
C 24 MAY 96 - WC  - POLINT: REMOVE DAMPING, AXE EFPEX IN EFPPAR
C  9 JAN 96 - WC  - REMOVE EXTRANEOUS OUTPUT
C 12 SEP 95 - WC  - FIX BUGS IN DIPIT AND FRGFLD
C  5 JAN 95 - PND - CHANGES FOR UNSYMMETRIC POLARIZABILITY TENSORS
C 23 NOV 94 - MWS - REMOVE ALL FTNCHEK ERRORS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C  1 AUG 94 - SPW - CHANGE ENERGY CORRECTION AND PERTURBATION TO H1
C 27 MAY 94 - SPW - EFPCMP: ADD IND. DIPOLE EACH SCF ITERATION
C 20 NOV 91 - JHJ - POLINT,REPINT,ZRFINT REPLACED IA ARRAY
C 13 NOV 91 - JHJ - ECPCMP: ENUCR=EN.
C 31 OCT 91 - JHJ - EFPCMP: CHANGED /CONV/, E PRINTOUT.
C 28 OCT 91 - JHJ - EFPCMP: DAF FILES REPLACE FILE 25 FOR ZRF.
C 23 OCT 91 - JHJ - EFPCMP: REPLACED ENERGY WITH ENRGY.
C                   INTRODUCED DIELEC AND ICALC TO COMMON /ZRFPAR/.
C 22 OCT 91 - JHJ - REPLACED OPTMIZ WITH GRAD THROUGHOUT.
C 23 JUL 91 - JHJ - ENUCPO:ADDED DRG'S NUCLEAR-INDUCED DIPOLE INTERAC.
C 19 JUL 91 - JHJ - ZRFINP: ADDED SCRF OPTION.
C 17 JUL 91 - JHJ - ZRFINP (FROM DRG'S EFPOTI).
C 16 JUL 91 - JHJ - ZDPMLT:CALC DIPOLES DUE TO FRG. MULTIPOLES.
C                   DIPIT:ITERATE INDUCED DIPOLE TO SELF CONSISTENCY.
C 15 JUL 91 - JHJ - FRGFLD:CALLED IN EFPCMP TO CALCULATE ELECTRIC
C                   FIELD DUE TO FRAGMENTS.
C 11 JUL 91 - JHJ - REPINT:INTRODUCED DRG'S REP. INTEGRAL CODE.
C  9 JUL 91 - JHJ - GATHERED EFPCMP AND POLINT TO FORM EFINTB.
C 20 FEB 91 - DRG - ZERNER REACTION FIELD ADDITION TO EFPCMP AND ENUCP.
C                   READ $ZRF IN EFPOTI FOR ZERNER REACTION FIELD SPEC.
C 18 FEB 91 - DRG - ZRFINT, INTEGRAL ROUTINE FOR ZERNER REACTION FIELD.
C  5 JUN 90 - DRG - DIPOLE FIELD INTEGRALS ADDED TO STANDV, CALLING
C  4 JUN 90 - DRG - MODIFIED POLINT FOR FIXED INDUCED MOMENTS.
C 26 MAY 90 - DRG - POLINT FOR EFP INTEGRALS.
C
C*MODULE EFINTB  *DECK DIPIT
      SUBROUTINE DIPIT(EFLD,EFADD,DINDC,DINDDC,DIND2C,DINDD2C,
     *                 DIND_W,DIND_C,DIND2_W,DIND2_C,N,
     *                 DIPNEW,DIPNWD,EWIFLD1,EWIFLD2,NXFRG)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      CHARACTER*8 POLNAM,DPOLNAM,DNAME,DPSNAM
C
      LOGICAL IFEWLD,TNFOIL,EWLDPL,DEBUGOPT
      LOGICAL GOPARR,MASWRK,DSKWRK
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MXCPUEFP=1024)
C
C        NOTE THAT -N- IS -NPTTPT-, AND -NXFRG- IS -NFRG-
C
      DIMENSION EFLD(3,N),EFADD(3,N),
     *          DINDC(3,N),DINDDC(3,N),DIND2C(3,N),DINDD2C(3,N),
     *          DIND_W(6,N),DIND_C(6,N),DIND2_W(6,N),DIND2_C(6,N),
     *          DIPNEW(NXFRG),DIPNWD(NXFRG),
     *          EWIFLD1(3,N),EWIFLD2(3,N)
C        NEXT ARE A STATIC ALLOCATION,
C        NOTE INITIALIZATION AND -SAVE- STATEMENT BELOW
      DIMENSION DIPLD(MXFRG),DIPDD(MXFRG)
C
      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 /EFPPBC/ XBOX,YBOX,ZBOX,XCUT,YCUT,ZCUT,RCUT
     *                ,SWR1,SWR2,ISWF
      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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
C
      COMMON /EFPAR / LOCNFG,NLOW,NUP
      COMMON /EFPAR2/ NFRGAL(MXCPUEFP),NLOWAL(MXCPUEFP),NUPAL(MXCPUEFP),
     *                NPPTAL(MXCPUEFP)
      COMMON /EWALD / IFEWLD,TNFOIL,BETA,ELCNST,EPSLN,CUTOFF,CUTLST,
     *                KMAX,NKVEC,LEVEL,EWLDPL
      COMMON /EWFLD / RCPFLD(3,MXFGPT),RCPIIFLD1(3,MXFGPT),
     *                RCPIIFLD2(3,MXFGPT), TNFLD(3),
     *                TNIIFLD1(3),TNIIFLD2(3),
     *                SELFFLD1(3,MXFGPT),SELFFLD2(3,MXFGPT),
     *                EXCLFLD(3,MXFGPT),
     *                EXCLIIFLD1(3,MXFGPT),EXCLIIFLD2(3,MXFGPT)
      COMMON /SCRPAR/ SCROFF
      INTEGER DDI_NP, DDI_ME
C
      PARAMETER (ZERO=0.00D+00, THREE=3.00D+00)
      PARAMETER (ONE=1.00D+00, TWO=2.00D+00)

C
      INTEGER RIGHT_N,LEFT_N,BUFF_PAS,BUFF_PAR
C
      SAVE DIPLD,DIPDD
      DATA DIPLD/MXFRG*0.0D+00/,DIPDD/MXFRG*0.0D+00/
      PI=4*ATAN(1.0D+00)
C
C     -- THIS ROUTINE IS LIFTED FROM "ELECT" WRITTEN BY DR. W.STEVENS.
C
C     --- EFLD  - TOTAL EXTERNAL FIELD EXCEPT FOR INDUCED DIPOLE ---
C     IDPUNC=0/1 MEANS DIPOLE ITERATION CONVERGED/UNCONVERGED
C
C
C     DIND2/DINDD2 ARE COMPUTED FOR MP2 LEVEL EFP1, AND
C     ARE IDENTICAL BECAUSE THE POLARIZABILITY TENSOR IS DIAGONAL
C     AND SYMMETRIC.
C
C     DIND AND DINDD ARE DIFFERENT FOR EFP1RHF AND EFP1DFT
C     BECAUSE THE POLARIZABILITY TENSORS ARE NOT SYMMETRIC
C
      CALL DDI_NPROC( DDI_NP, DDI_ME )
C
      RIGHT_N = MOD(DDI_ME+1,DDI_NP)
      LEFT_N = MOD(DDI_ME-1+DDI_NP,DDI_NP)
C
      NIN1 = 1
      DO I = 1,NLOW-1
        NIN1 = NIN1 + NPPTS(I)
      ENDDO
      NIN2=0
      DO I = NLOW,NUP
        NIN2=NPPTS(I)+NIN2
      ENDDO
C
      IF(GOPARR) THEN
        CALL VICLR(NPPTAL,1,MXCPUEFP)
        NPPTAL(DDI_ME+1) = NIN2
        CALL DDI_GSUMI(664,NPPTAL,DDI_NP)
      ELSE
        NPPTAL(1) = NIN2
      END IF
C
      IDPUNC = 0
C
      IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
         CALL VCLR(DIND2_C,1,6*NPTTPT)
         CALL VCLR(DIND2_W,1,6*NPTTPT)
      ELSE
         CALL VCLR(DIND_C,1,6*NPTTPT)
         CALL VCLR(DIND_W,1,6*NPTTPT)
      END IF
C
C          TOP OF BIG LOOP OVER SELF CONSISTENCY
C
C     PARALLEL AND SERIAL ALGORITHMS CONVERGE TO THE SAME RESULTS
C     IF AND ONLY IF WE USE A TIGHT CRITERION = 1.0D-10
C
C     PARALLEL ALGORITHM REQUIRES MORE ITERS (INCREASED 30->60 12/2007)
C
      ITLIM= 60
      CONV = 1.0D-10
      ITER = 0
  390 CONTINUE
      ITER = ITER + 1
C
      IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
         CALL VCLR(DIND2C, 1,3*NPTTPT)
         CALL VCLR(DINDD2C,1,3*NPTTPT)
      ELSE
         CALL VCLR(DINDC, 1,3*NPTTPT)
         CALL VCLR(DINDDC,1,3*NPTTPT)
      END IF
      CALL VCLR(DIPNEW,1,NFRG)
      CALL VCLR(DIPNWD,1,NFRG)
      KP = NIN1
C
C   ----- CONTRINUTION FROM EWALD TERMS
      IF (IFEWLD .AND. EWLDPL) THEN 
         CALL VCLR (EWIFLD1,1,3*NPTTPT)
         CALL VCLR (EWIFLD2,1,3*NPTTPT)
         CALL VCLR (RCPIIFLD1,1,3*NPTTPT)
         CALL VCLR (RCPIIFLD2,1,3*NPTTPT)
         CALL VCLR (TNIIFLD1,1,3*NPTTPT)
         CALL VCLR (TNIIFLD2,1,3*NPTTPT)
         CALL VCLR (SELFFLD1,1,3*NPTTPT)
         CALL VCLR (SELFFLD2,1,3*NPTTPT)
         CALL VCLR (EXCLIIFLD1,1,3*NPTTPT)
         CALL VCLR (EXCLIIFLD2,1,3*NPTTPT)
C    ---  DEBUGGING FUNCTION: CALCULATES INDUCED FIELD EXPLICITELY
         DEBUGOPT=.FALSE.
         IF (DEBUGOPT) CALL DIRECTSUM_IND
         CALL RECIPR_IIFLD
C   ---- SELF- AND EXCLUSION CONTRIBUTIONS TO THE INDUCED FIELD
         CALL SELF_FIELD
         CALL TINFOIL_FIELD
         DO K=1,NPTTPT
            DO I=1,3
            EWIFLD1(I,K)= RCPIIFLD1(I,K)+TNIIFLD1(I)
     *                   +SELFFLD1(I,K)
     *                   +EXCLIIFLD1(I,K)       

            EWIFLD2(I,K)= RCPIIFLD2(I,K)+TNIIFLD2(I)
     *                   +SELFFLD2(I,K)
     *                   +EXCLIIFLD2(I,K)     
            END DO
         END DO
C         WRITE(6,*)'TOTAL FIELD 1',(EWIFLD1(I,1),I=1,3),
C     *              (EWIFLD1(I,2),I=1,3),(EWIFLD1(I,3),I=1,3)     
C         WRITE(6,*)'TOTAL FIELD 2',(EWIFLD2(I,1),I=1,3),
C     *              (EWIFLD2(I,2),I=1,3),(EWIFLD2(I,3),I=1,3)              
      END IF
      
      IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
         DO I = 1,NIN2
            DIND2_C(1,I) = DIND2(1,KP)
            DIND2_C(2,I) = DIND2(2,KP)
            DIND2_C(3,I) = DIND2(3,KP)
            DIND2_C(4,I) = DINDD2(1,KP)
            DIND2_C(5,I) = DINDD2(2,KP)
            DIND2_C(6,I) = DINDD2(3,KP)
            KP = KP + 1
         ENDDO
      ELSE
         DO I = 1,NIN2
            DIND_C(1,I) = DIND(1,KP)
            DIND_C(2,I) = DIND(2,KP)
            DIND_C(3,I) = DIND(3,KP)
            DIND_C(4,I) = DINDD(1,KP)
            DIND_C(5,I) = DINDD(2,KP)
            DIND_C(6,I) = DINDD(3,KP)
            KP = KP + 1
         ENDDO
      END IF
C
      BUFF_PAS = DDI_ME
      DO 340 I = 1,DDI_NP
C
         IF (GOPARR .AND. I .NE. DDI_NP) THEN
            BUFF_PAR = MOD(DDI_ME-I+DDI_NP,DDI_NP)
            IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
               CALL DDI_IRECV(DIND2_W,6*NPPTAL(BUFF_PAR+1)*8,
     *                        LEFT_N,IREQRCV)
               CALL DDI_ISEND(DIND2_C,6*NPPTAL(BUFF_PAS+1)*8,
     *                        RIGHT_N,IREQSND)
            ELSE
               CALL DDI_IRECV(DIND_W,6*NPPTAL(BUFF_PAR+1)*8,
     *                        LEFT_N,IREQRCV)
               CALL DDI_ISEND(DIND_C,6*NPPTAL(BUFF_PAS+1)*8,
     *                        RIGHT_N,IREQSND)
            END IF
         END IF
C
         INDEX = NIN1
C
C      PARALLELIZATION OF ORIGINAL LOOP, IM = 1,NFRG
C
         DO 330 IM = NLOW,NUP
             COMX1 = EFCENT(1,IM)
             COMY1 = EFCENT(2,IM)
             COMZ1 = EFCENT(3,IM)
C
            I1 = INDEX
            INDEX = INDEX + NPPTS(IM)
            I2 = I1 + NPPTS(IM) - 1
C
            IF(I.EQ.1) THEN
               DIPX  = ZERO
               DIPY  = ZERO
               DIPZ  = ZERO
               DDIPX = ZERO
               DDIPY = ZERO
               DDIPZ = ZERO
            END IF
C
            DO 320 IPX = I1,I2
               IF(BUFF_PAS .EQ. DDI_ME) THEN
                  IF (IFEWLD .AND. EWLDPL) THEN 
                     EFX  = EFLD(1,IPX)+EWIFLD1(1,IPX)
                     EFY  = EFLD(2,IPX)+EWIFLD1(2,IPX)
                     EFZ  = EFLD(3,IPX)+EWIFLD1(3,IPX)
                     EFXD = EFLD(1,IPX)+EWIFLD2(1,IPX)
                     EFYD = EFLD(2,IPX)+EWIFLD2(2,IPX)
                     EFZD = EFLD(3,IPX)+EWIFLD2(3,IPX)
                  ELSE
                     EFX  = EFLD(1,IPX)
                     EFY  = EFLD(2,IPX)
                     EFZ  = EFLD(3,IPX)
                     EFXD = EFLD(1,IPX)
                     EFYD = EFLD(2,IPX)
                     EFZD = EFLD(3,IPX)
                  END IF
               ELSE
                  EFX  =ZERO
                  EFY  =ZERO
                  EFZ  =ZERO
                  EFXD =ZERO
                  EFYD =ZERO
                  EFZD =ZERO
               END IF
C               WRITE(6,*)'EFX1=',IPX,EFX,EFY,EFZ
C               WRITE(6,*)'EFX2=',IPX,EFXD,EFYD,EFZD
               X1 = EFP(1,IPX)
               Y1 = EFP(2,IPX)
               Z1 = EFP(3,IPX)
C
               INDEX2 = 1
               K1 = 1
               DO J = 1,NLOWAL(BUFF_PAS+1)-1
                 INDEX2 = INDEX2 + NPPTS(J)
               ENDDO
C                   PARALLELIZE ORIGINAL LOOP JM = 1,NFRG
               DO 310 JM = NLOWAL(BUFF_PAS+1),NUPAL(BUFF_PAS+1)
                  IF (IM.EQ.JM) THEN
                     K1 = K1 + NPPTS(JM)
                     GO TO 309
                  END IF
C
                  COMX2 = EFCENT(1,JM)
                  COMY2 = EFCENT(2,JM)
                  COMZ2 = EFCENT(3,JM)
C
                  COMX = COMX1 - COMX2
                  COMY = COMY1 - COMY2
                  COMZ = COMZ1 - COMZ2
C
                  PBCX = ZERO
                  PBCY = ZERO
                  PBCZ = ZERO
               IF (XBOX .NE. ZERO) PBCX = XBOX * ANINT(COMX/XBOX)
               IF (YBOX .NE. ZERO) PBCY = YBOX * ANINT(COMY/YBOX)
               IF (ZBOX .NE. ZERO) PBCZ = ZBOX * ANINT(COMZ/ZBOX)
C
C   ----  DISTANCE BETWEEN FRAGMENTS
               X12  = ABS(COMX - PBCX)
               Y12  = ABS(COMY - PBCY)
               Z12  = ABS(COMZ - PBCZ)
               R12  = X12*X12+Y12*Y12+Z12*Z12
C
               IF(X12.GT.XCUT.OR.Y12.GT.YCUT.OR.Z12.GT.ZCUT) THEN
                 K1 = K1 + NPPTS(JM)
                 GOTO 309
               END IF
               R12  = X12*X12+Y12*Y12+Z12*Z12
               IF(R12.GT.RCUT) THEN
                 K1 = K1 + NPPTS(JM)
                 GOTO 309
               END IF
C
               SWF   = 1.0D+00
               IF(R12.GT.SWR2*SWR2) THEN
                 K1 = K1 + NPPTS(JM)
                 GOTO 309
               ELSE IF(R12.GT.SWR1*SWR1) THEN
                 IF(ISWF.EQ.1) THEN
                    DUMY1 = 1.0D+00/(SWR2*SWR2 - SWR1*SWR1)
                    DUMY2 = DUMY1*DUMY1
                    DUMY3 = DUMY2*DUMY1
                    DUMY4 = DUMY2*DUMY2
                    DUMY5 = DUMY2*DUMY3
                    VUMY1 = R12 - SWR1*SWR1
                    VUMY2 = VUMY1*VUMY1
                    VUMY3 = VUMY2*VUMY1
                    VUMY4 = VUMY2*VUMY2
                    VUMY5 = VUMY2*VUMY3
                    DSWF  = -30.0D+00*DUMY3*VUMY2
     *                      +60.0D+00*DUMY4*VUMY3
     *                      -30.0D+00*DUMY5*VUMY4
                    DSWF  = 2.0D+00*DSWF
                    SWF   = 1.0D+00 - 10.0D+00*DUMY3*VUMY3
     *                              + 15.0D+00*DUMY4*VUMY4
     *                              -  6.0D+00*DUMY5*VUMY5
                 ELSE IF(ISWF.EQ.2) THEN
                    R12   = SQRT(R12)
                    PIL   = 3.1415926535897932D+00/(SWR2-SWR1)
                    DUMY  = (R12 - SWR1) * PIL
                    SWF   = 0.50D+00 * (COS(DUMY) + 1.0D+00)
                 END IF
               END IF
C
                  J1 = INDEX2
                  J2 = J1 + NPPTS(JM) - 1
                  DO 300 JP = J1,J2
                     X2 = EFP(1,JP)
                     Y2 = EFP(2,JP)
                     Z2 = EFP(3,JP)
                     X = X1 - X2
                     Y = Y1 - Y2
                     Z = Z1 - Z2
C
                     X = X - PBCX
                     Y = Y - PBCY
                     Z = Z - PBCZ
C
                     R2 = X * X + Y * Y + Z * Z
                     R = SQRT(R2)
                     R3 = R * R * R
                     R5 = R * R * R3
C
C  ----  LVS: ELECTROSTATIC FRAG-FRAG SCREENING
C  ----  ZERO OUT ALL PENETRATION TERMS FIRST
                     P1 = ZERO
                     P1D = ZERO
                     P2 = ZERO
C  ----  IF DO FRAGMENT-FRAGMENT SCREENING
                     IF(IPLSCR.EQ.1) THEN
C                     IF(1.EQ.0) THEN
                       ALPHAB = POLSCR(JP)
                       ALPHAA = POLSCR(IPX)
C
                       IF(R.LE.SCROFF) THEN
                         AP=ALPHAA*R
                         BP=ALPHAB*R
                         AA2=ALPHAA*ALPHAA
                         AB2=ALPHAB*ALPHAB
                         DIFF=ABS(ALPHAA-ALPHAB)
C
C  ----  NOW CALCULATE ALL FRAG-FRAG CHARGE-PENETRATION TERMS
                         CALL EPENCHDIP(AP,BP,AA2,AB2,DIFF,
     *                                  P1,P1_N1,P1_N2)
                         CALL EPENCHQU(AP,BP,AA2,AB2,DIFF,R2,
     *                            P1D,P1D_N1,P1D_N2)
C      LVS DIFFERENT VARIANTS OF SCREENING
C                         CALL SCR5CHDIP(AP,BP,AA2,AB2,DIFF,
C     *                                  P1,P1_N1,P1_N2)
C                         CALL SCR5CHQU(AP,BP,AA2,AB2,DIFF,R2,
C     *                            P1D,P1D_N1,P1D_N2)
                         P2 =  -P1D/R3
                       END IF
                     END IF
C   ---- END OF ELECTROSTATIC SCREENING PART
C
                     IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
C
                        IF(GOPARR) THEN
                           DOT = DIND2_C(1,K1)*X
     *                         + DIND2_C(2,K1)*Y
     *                         + DIND2_C(3,K1)*Z
                           DOT = THREE * DOT / R5
                           EFX = EFX- SWF*DIND2_C(1,K1)/R3 + SWF*DOT*X
                           EFY = EFY- SWF*DIND2_C(2,K1)/R3 + SWF*DOT*Y
                           EFZ = EFZ- SWF*DIND2_C(3,K1)/R3 + SWF*DOT*Z
C
                           DOTD = DIND2_C(4,K1)*X
     *                          + DIND2_C(5,K1)*Y
     *                          + DIND2_C(6,K1)*Z
                           DOTD = THREE * DOTD / R5
                           EFXD =EFXD-SWF*DIND2_C(4,K1)/R3+SWF*DOTD*X
                           EFYD =EFYD-SWF*DIND2_C(5,K1)/R3+SWF*DOTD*Y
                           EFZD =EFZD-SWF*DIND2_C(6,K1)/R3+SWF*DOTD*Z
                           K1 = K1 + 1
                        ELSE
                           DOT = DIND2(1,JP)*X
     *                         + DIND2(2,JP)*Y
     *                         + DIND2(3,JP)*Z
                           DOT = THREE * DOT / R5
                           EFX = EFX - SWF*DIND2(1,JP)/R3 + SWF*DOT*X
                           EFY = EFY - SWF*DIND2(2,JP)/R3 + SWF*DOT*Y
                           EFZ = EFZ - SWF*DIND2(3,JP)/R3 + SWF*DOT*Z
C
                           DOTD = DINDD2(1,JP)*X
     *                          + DINDD2(2,JP)*Y
     *                          + DINDD2(3,JP)*Z
                           DOTD = THREE * DOTD / R5
                           EFXD = EFXD-SWF*DINDD2(1,JP)/R3+SWF*DOTD*X
                           EFYD = EFYD-SWF*DINDD2(2,JP)/R3+SWF*DOTD*Y
                           EFZD = EFZD-SWF*DINDD2(3,JP)/R3+SWF*DOTD*Z
                        END IF
                     ELSE
C
C  ----  LVS: SCREENING FOR FRAGMENT-FRAGMENT PART ADDED
C-------- EWALD PART
                        IF (IFEWLD .AND. EWLDPL) THEN 
C
                        F0 = ERFC ( BETA * R )
                        EXPBR = EXP(-BETA * BETA * R2)
                        F1 = TWO*BETA*R/SQRT(PI)*EXPBR + F0
                        F2 = 4.0D+00*BETA*BETA*BETA/SQRT(PI)*EXPBR/R2
C
                        IF(GOPARR) THEN
                           DOT1 = DIND_C(1,K1)*X
     *                          + DIND_C(2,K1)*Y
     *                          + DIND_C(3,K1)*Z
                           DOT = THREE * DOT1 / R5
                           EFX = EFX+ (F1+P1)*(-DIND_C(1,K1)/R3+DOT*X)
     *                               +(P2+F2)*DOT1*X
                           EFY = EFY+ (F1+P1)*(-DIND_C(2,K1)/R3+DOT*Y)
     *                               +(F2+P2)*DOT1*Y
                           EFZ = EFZ+ (F1+P1)*(-DIND_C(3,K1)/R3+DOT*Z)
     *                               +(F2+P2)*DOT1*Z
C
                           DOTD1 = DIND_C(4,K1)*X
     *                           + DIND_C(5,K1)*Y
     *                           + DIND_C(6,K1)*Z
                           DOTD = THREE * DOTD1 / R5
                           EFXD = EFXD+(F1+P1)*(-DIND_C(4,K1)/R3+DOTD*X)       
     *                                 +(F2+P2)*DOTD1*X
                           EFYD = EFYD+(F1+P1)*(-DIND_C(5,K1)/R3+DOTD*Y)
     *                                 +(F2+P2)*DOTD1*Y
                           EFZD = EFZD+(F1+P1)*(-DIND_C(6,K1)/R3+DOTD*Z)
     *                                 +(F2+P2)*DOTD1*Z
                           K1 = K1 + 1
                        ELSE
                           DOT1 = DIND(1,JP)*X
     *                          + DIND(2,JP)*Y
     *                          + DIND(3,JP)*Z
                           DOT = THREE * DOT1 / R5
                           EFX = EFX+ (F1+P1)*
     *                          (-DIND(1,JP)/R3 + DOT*X)+(F2+P2)*DOT1*X
                           EFY = EFY+ (F1+P1)*
     *                          (-DIND(2,JP)/R3 + DOT*Y)+(F2+P2)*DOT1*Y
                           EFZ = EFZ+ (F1+P1)*
     *                          (-DIND(3,JP)/R3 + DOT*Z)+(F2+P2)*DOT1*Z
C
                           DOTD1 = DINDD(1,JP)*X
     *                           + DINDD(2,JP)*Y
     *                           + DINDD(3,JP)*Z
                           DOTD = THREE * DOTD1 / R5
                           EFXD = EFXD+ (F1+P1)*
     *                       (-DINDD(1,JP)/R3 + DOTD*X)+(F2+P2)*DOTD1*X
                           EFYD = EFYD+ (F1+P1)*
     *                       (-DINDD(2,JP)/R3 + DOTD*Y)+(F2+P2)*DOTD1*Y
                           EFZD = EFZD+ (F1+P1)*
     *                       (-DINDD(3,JP)/R3 + DOTD*Z)+(F2+P2)*DOTD1*Z
                        END IF
C------- END OF EWALD PART
                        ELSE
                        IF(GOPARR) THEN
                           DOT1 = DIND_C(1,K1)*X
     *                          + DIND_C(2,K1)*Y
     *                          + DIND_C(3,K1)*Z
                           DOT = THREE * DOT1 / R5
                           EFX = EFX+ SWF*((ONE+P1)*
     *                          (-DIND_C(1,K1)/R3 + DOT*X)+P2*DOT1*X)
                           EFY = EFY+ SWF*((ONE+P1)*
     *                          (-DIND_C(2,K1)/R3 + DOT*Y)+P2*DOT1*Y)
                           EFZ = EFZ+ SWF*((ONE+P1)*
     *                          (-DIND_C(3,K1)/R3 + DOT*Z)+P2*DOT1*Z)
C
                           DOTD1 = DIND_C(4,K1)*X
     *                           + DIND_C(5,K1)*Y
     *                           + DIND_C(6,K1)*Z
                           DOTD = THREE * DOTD1 / R5
                           EFXD = EFXD+ SWF*((ONE+P1)*
     *                          (-DIND_C(4,K1)/R3 + DOTD*X)+P2*DOTD1*X)
                           EFYD = EFYD+ SWF*((ONE+P1)*
     *                          (-DIND_C(5,K1)/R3 + DOTD*Y)+P2*DOTD1*Y)
                           EFZD = EFZD+ SWF*((ONE+P1)*
     *                          (-DIND_C(6,K1)/R3 + DOTD*Z)+P2*DOTD1*Z)
                           K1 = K1 + 1
                        ELSE
                           DOT1 = DIND(1,JP)*X
     *                          + DIND(2,JP)*Y
     *                          + DIND(3,JP)*Z
                           DOT = THREE * DOT1 / R5
                           EFX = EFX+ SWF*((ONE+P1)*
     *                          (-DIND(1,JP)/R3 + DOT*X)+P2*DOT1*X)
                           EFY = EFY+ SWF*((ONE+P1)*
     *                          (-DIND(2,JP)/R3 + DOT*Y)+P2*DOT1*Y)
                           EFZ = EFZ+ SWF*((ONE+P1)*
     *                          (-DIND(3,JP)/R3 + DOT*Z)+P2*DOT1*Z)
C
                           DOTD1 = DINDD(1,JP)*X
     *                           + DINDD(2,JP)*Y
     *                           + DINDD(3,JP)*Z
                           DOTD = THREE * DOTD1 / R5
                           EFXD = EFXD+ SWF*((ONE+P1)*
     *                          (-DINDD(1,JP)/R3 + DOTD*X)+P2*DOTD1*X)
                           EFYD = EFYD+ SWF*((ONE+P1)*
     *                          (-DINDD(2,JP)/R3 + DOTD*Y)+P2*DOTD1*Y)
                           EFZD = EFZD+ SWF*((ONE+P1)*
     *                          (-DINDD(3,JP)/R3 + DOTD*Z)+P2*DOTD1*Z)
                        END IF
                        END IF
                     END IF
  300             CONTINUE
  309             CONTINUE
                  INDEX2 = INDEX2 + NPPTS(JM)
  310          CONTINUE
C
C               WRITE(6,*)'FINAL EFX1=',IPX,EFX,EFY,EFZ
C               WRITE(6,*)'FINAL EFX2=',IPX,EFXD,EFYD,EFZD
C
               IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
                  IF(IPX.EQ.I2) THEN
                     IF(GOPARR) THEN
                        DIND2C(1,IPX) =DIND2C(1,IPX) +EFX*EFPOL2(1,IM)
     *                                               +EFY*EFPOL2(4,IM)
     *                                               +EFZ*EFPOL2(5,IM)
                        DIND2C(2,IPX) =DIND2C(2,IPX) +EFX*EFPOL2(7,IM)
     *                                               +EFY*EFPOL2(2,IM)
     *                                               +EFZ*EFPOL2(6,IM)
                        DIND2C(3,IPX) =DIND2C(3,IPX) +EFX*EFPOL2(8,IM)
     *                                               +EFY*EFPOL2(9,IM)
     *                                               +EFZ*EFPOL2(3,IM)
                        DINDD2C(1,IPX)=DINDD2C(1,IPX)+EFXD*EFPOL2(1,IM)
     *                                               +EFYD*EFPOL2(7,IM)
     *                                               +EFZD*EFPOL2(8,IM)
                        DINDD2C(2,IPX)=DINDD2C(2,IPX)+EFXD*EFPOL2(4,IM)
     *                                               +EFYD*EFPOL2(2,IM)
     *                                               +EFZD*EFPOL2(9,IM)
                        DINDD2C(3,IPX)=DINDD2C(3,IPX)+EFXD*EFPOL2(5,IM)
     *                                               +EFYD*EFPOL2(6,IM)
     *                                               +EFZD*EFPOL2(3,IM)
                     ELSE
                        DIND2(1,IPX) = EFX*EFPOL2(1,IM)
     *                               + EFY*EFPOL2(4,IM)
     *                               + EFZ*EFPOL2(5,IM)
                        DIND2(2,IPX) = EFX*EFPOL2(7,IM)
     *                               + EFY*EFPOL2(2,IM)
     *                               + EFZ*EFPOL2(6,IM)
                        DIND2(3,IPX) = EFX*EFPOL2(8,IM)
     *                               + EFY*EFPOL2(9,IM)
     *                               + EFZ*EFPOL2(3,IM)
                        DINDD2(1,IPX)= EFXD*EFPOL2(1,IM)
     *                               + EFYD*EFPOL2(7,IM)
     *                               + EFZD*EFPOL2(8,IM)
                        DINDD2(2,IPX)= EFXD*EFPOL2(4,IM)
     *                               + EFYD*EFPOL2(2,IM)
     *                               + EFZD*EFPOL2(9,IM)
                        DINDD2(3,IPX)= EFXD*EFPOL2(5,IM)
     *                               + EFYD*EFPOL2(6,IM)
     *                               + EFZD*EFPOL2(3,IM)
                     END IF
                     IF(I.EQ.1) THEN
                        DIPX  = DIPX + DIND2(1,IPX)
                        DIPY  = DIPY + DIND2(2,IPX)
                        DIPZ  = DIPZ + DIND2(3,IPX)
                        DDIPX = DDIPX + DINDD2(1,IPX)
                        DDIPY = DDIPY + DINDD2(2,IPX)
                        DDIPZ = DDIPZ + DINDD2(3,IPX)
                     END IF
                  END IF
               ELSE
                  IF(GOPARR) THEN
                     DINDC(1,IPX) = DINDC(1,IPX)+ EFX*EFPOL(1,IPX)
     *                                          + EFY*EFPOL(4,IPX)
     *                                          + EFZ*EFPOL(5,IPX)
                     DINDC(2,IPX) = DINDC(2,IPX)+ EFX*EFPOL(7,IPX)
     *                                          + EFY*EFPOL(2,IPX)
     *                                          + EFZ*EFPOL(6,IPX)
                     DINDC(3,IPX) = DINDC(3,IPX)+ EFX*EFPOL(8,IPX)
     *                                          + EFY*EFPOL(9,IPX)
     *                                          + EFZ*EFPOL(3,IPX)
                     DINDDC(1,IPX)=DINDDC(1,IPX)+ EFXD*EFPOL(1,IPX)
     *                                          + EFYD*EFPOL(7,IPX)
     *                                          + EFZD*EFPOL(8,IPX)
                     DINDDC(2,IPX)=DINDDC(2,IPX)+ EFXD*EFPOL(4,IPX)
     *                                          + EFYD*EFPOL(2,IPX)
     *                                          + EFZD*EFPOL(9,IPX)
                     DINDDC(3,IPX)=DINDDC(3,IPX)+ EFXD*EFPOL(5,IPX)
     *                                          + EFYD*EFPOL(6,IPX)
     *                                          + EFZD*EFPOL(3,IPX)
                  ELSE
                     DIND(1,IPX) = EFX*EFPOL(1,IPX)
     *                           + EFY*EFPOL(4,IPX)
     *                           + EFZ*EFPOL(5,IPX)
                     DIND(2,IPX) = EFX*EFPOL(7,IPX)
     *                           + EFY*EFPOL(2,IPX)
     *                           + EFZ*EFPOL(6,IPX)
                     DIND(3,IPX) = EFX*EFPOL(8,IPX)
     *                           + EFY*EFPOL(9,IPX)
     *                           + EFZ*EFPOL(3,IPX)
                     DINDD(1,IPX)= EFXD*EFPOL(1,IPX)
     *                           + EFYD*EFPOL(7,IPX)
     *                           + EFZD*EFPOL(8,IPX)
                     DINDD(2,IPX)= EFXD*EFPOL(4,IPX)
     *                           + EFYD*EFPOL(2,IPX)
     *                           + EFZD*EFPOL(9,IPX)
                     DINDD(3,IPX)= EFXD*EFPOL(5,IPX)
     *                           + EFYD*EFPOL(6,IPX)
     *                           + EFZD*EFPOL(3,IPX)
                  END IF
                  IF(I.EQ.1) THEN
                     DIPX  = DIPX + DIND(1,IPX)
                     DIPY  = DIPY + DIND(2,IPX)
                     DIPZ  = DIPZ + DIND(3,IPX)
                     DDIPX = DDIPX + DINDD(1,IPX)
                     DDIPY = DDIPY + DINDD(2,IPX)
                     DDIPZ = DDIPZ + DINDD(3,IPX)
                  END IF
               END IF
  320       CONTINUE
C
            IF(I.EQ.1) THEN
               DIPNEW(IM) = DIPX * DIPX + DIPY * DIPY + DIPZ * DIPZ
               DIPNEW(IM) = SQRT(DIPNEW(IM))
               DIPNWD(IM) = DDIPX*DDIPX + DDIPY*DDIPY + DDIPZ*DDIPZ
               DIPNWD(IM) = SQRT(DIPNWD(IM))
            END IF
  330    CONTINUE
C
         IF(GOPARR .AND. I .NE. DDI_NP) THEN
            CALL DDI_WAIT(IREQRCV)
            CALL DDI_WAIT(IREQSND)
C
            IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
               DO IJ = 1,NPPTAL(BUFF_PAR+1)
                  DIND2_C(1,IJ) = DIND2_W(1,IJ)
                  DIND2_C(2,IJ) = DIND2_W(2,IJ)
                  DIND2_C(3,IJ) = DIND2_W(3,IJ)
                  DIND2_C(4,IJ) = DIND2_W(4,IJ)
                  DIND2_C(5,IJ) = DIND2_W(5,IJ)
                  DIND2_C(6,IJ) = DIND2_W(6,IJ)
               ENDDO
            ELSE
               DO IJ = 1,NPPTAL(BUFF_PAR+1)
                  DIND_C(1,IJ) = DIND_W(1,IJ)
                  DIND_C(2,IJ) = DIND_W(2,IJ)
                  DIND_C(3,IJ) = DIND_W(3,IJ)
                  DIND_C(4,IJ) = DIND_W(4,IJ)
                  DIND_C(5,IJ) = DIND_W(5,IJ)
                  DIND_C(6,IJ) = DIND_W(6,IJ)
               ENDDO
            END IF
C
            BUFF_PAS = BUFF_PAR
         END IF
C
  340 CONTINUE   ! END LOOP OVER ALL PROCESSORS DDI_NP
C
C     ----- CHECK INDUCED DIPOLE AGAINST CONVERGENCE CRITERIA -----
C
      ICONV = 0
C
C               PARALLELIZATION OF ORIGINAL LOOP, I = 1,NFRG
      DO 350 I = NLOW,NUP
         DIFF1 = DIPLD(I) - DIPNEW(I)
         DIFF2 = DIPDD(I) - DIPNWD(I)
         DIPLD(I) = DIPNEW(I)
         DIPDD(I) = DIPNWD(I)
         DIFF1 = ABS(DIFF1)
         DIFF2 = ABS(DIFF2)
         IF (DIFF1.LT.CONV) ICONV = ICONV + 1
         IF (DIFF2.LT.CONV) ICONV = ICONV + 1
  350 CONTINUE
C
      CALL DDI_GSUMI(668,ICONV,1)
C
C     HUI: ITER MUST BE > 2
      IF(ITER.LT.2) ICONV = 0
C
      IF (ICONV.NE.2*NFRG .AND. ITER.LE.ITLIM) THEN
         IF(GOPARR) THEN
            IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
               CALL VCLR(DIND2, 1,3*NPTTPT)
               CALL VCLR(DINDD2,1,3*NPTTPT)
               DO IPX=NIN1,NIN1+NIN2-1
                  DIND2(1,IPX) = DIND2C(1,IPX)
                  DIND2(2,IPX) = DIND2C(2,IPX)
                  DIND2(3,IPX) = DIND2C(3,IPX)
                  DINDD2(1,IPX) = DINDD2C(1,IPX)
                  DINDD2(2,IPX) = DINDD2C(2,IPX)
                  DINDD2(3,IPX) = DINDD2C(3,IPX)
               ENDDO
            ELSE
               CALL VCLR(DIND, 1,3*NPTTPT)
               CALL VCLR(DINDD,1,3*NPTTPT)
               DO IPX=NIN1,NIN1+NIN2-1
                  DIND(1,IPX) = DINDC(1,IPX)
                  DIND(2,IPX) = DINDC(2,IPX)
                  DIND(3,IPX) = DINDC(3,IPX)
                  DINDD(1,IPX) = DINDDC(1,IPX)
                  DINDD(2,IPX) = DINDDC(2,IPX)
                  DINDD(3,IPX) = DINDDC(3,IPX)
               ENDDO
            END IF  ! IEFPOL2=1/0 BLOCKS
         END IF  ! GOPARR
C
C             BRANCH BACK TO BEGINNING OF INDUCED DIPOLE ITERATIONS
C
         GO TO 390
      END IF
C
      IF (ICONV.NE.2*NFRG) THEN
         IF(MASWRK) WRITE(6,9010) ITER,NACC
         IF(NACC.EQ.0) CALL ABRT
         IDPUNC=1
         RETURN
      END IF
C
C
      IF(GOPARR) THEN
         IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
            CALL DDI_GSUMF(770,DINDD2,3*NPTTPT)
            CALL DDI_GSUMF(771,DIND2, 3*NPTTPT)
         ELSE
            CALL DDI_GSUMF(770,DINDD ,3*NPTTPT)
            CALL DDI_GSUMF(771,DIND , 3*NPTTPT)
         END IF
      END IF
C      WRITE(6,*)'DIND: ',(DIND(J,I),I=1,NPTTPT)
C      WRITE(6,*)'DINDD:',(DINDD(J,I),I=1,NPTTPT)
C
C     -- USE CONVERGED INDUCED DIPOLES TO CALCULATE THEIR FIELDS
C        ONLY WHEN NAT>0, EFADD WILL BE COMPUTED IN FLDIND.
C        EFADD IS USED ONLY FOR QM/EFP
C
      CALL FLDIND(EFADD)
C
      RETURN
C
 9010 FORMAT(1X,'DIPIT FAILED TO CONVERGE IN',I5,'ITERATIONS, NACC=',I4)
      END
C*MODULE EFINTB  *DECK EFPCMP
      SUBROUTINE EFPCMP(H1,FA,FB,DENTOT,DIP,ODIP,
     *                  XEFI,YEFI,ZEFI,EFLD,EFADD,ABFLD,L2,L1,MCITER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000,
     *           MXRT=100, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00, ONE=1.0D+00,
     *           DEBYE=2.541766D+00)
C
      DIMENSION H1(L2),FA(L2),FB(L2),DENTOT(L2),
     *          DIP(L2),ODIP(L2),XEFI(L2),YEFI(L2),ZEFI(L2),
     *          EFLD(3,NPTTPT),EFADD(3,NPTTPT),ABFLD(3,NPTTPT)
      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 / XX(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 /PNUC  / PPNX(MXFGPT),PPNY(MXFGPT),PPNZ(MXFGPT),PNZRF(3)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA RHF,GVB/8HRHF     ,8HGVB     /
      DATA RMC/8HMCSCF   /
C
      SAVE EN0
C
C     THIS ROUTINE INDUCES THE EFP FIELD, AND MODIFIES THE 1E- PART
C     OF THE FOCK MATRIX TO REFLECT THE PRESENCE OF EFP'S.
C
C     NOTE TO PROGRAMMERS:  THIS ROUTINE HAS BEEN COPIED TO
C     MAKE SMALL ADDITIONS.  IF YOU CHANGE SOMETHING IN THE PARENT
C     ROUTINE -EFPCMP-, DO THE SAME THING TO -CHGIND- AND CONVERSELY.
C
C     --- EFLD  - TOTAL EXTERNAL FIELD EXCEPT FOR INDUCED DIPOLE ---
C
C     CAUTION, GVB ITERATIONS START AT 1, BUT RHF,UHF,ROHF AT 0!
C     EN0 IS USUAL TERM
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)
         CALL FRGFLD
      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
      NFT25 = 25
      SVDSKW = DSKWRK
      DSKWRK = .TRUE.
      CALL SEQREW(NFT25)
      IPCOUNT = ME-1
      DO 10 LEFP = 1, NPTTPT
         IF(GOPARR) THEN
            ABFLD(1,LEFP)=ZERO
            ABFLD(2,LEFP)=ZERO
            ABFLD(3,LEFP)=ZERO
            IPCOUNT = IPCOUNT + 1
            IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 10
         END IF
         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)
 10   CONTINUE
      CALL SEQREW(NFT25)
      IF(GOPARR) CALL DDI_GSUMF(668,ABFLD,3*NPTTPT)
C
C     -- TOTAL FIELDS --
C
      DO LEFP = 1, NPTTPT
         EFLD(1,LEFP)=ABFLD(1,LEFP) + PPNX(LEFP) + PPEFX(LEFP)
         EFLD(2,LEFP)=ABFLD(2,LEFP) + PPNY(LEFP) + PPEFY(LEFP)
         EFLD(3,LEFP)=ABFLD(3,LEFP) + PPNZ(LEFP) + PPEFZ(LEFP)
      ENDDO
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,XX(LDINDC),XX(LDINDDC),XX(LDIND2C),
     *           XX(LDINDD2C),XX(LDIND_W),XX(LDIND_C),XX(LDIND2_W),
     *           XX(LDIND2_C),NPTTPT,XX(LDIPNEW),XX(LDIPNWD),
     *           XX(LEWLD1),XX(LEWLD2),NFRG)
      CALL RETFM(NEED)
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        - WARNING IF INDUCED DIPOLE IS TOO LARGE
C          WHY NOT SCALE IT BACK TO THE MAXIMUM, ASKS MIKE?
C          SOOHAENG RAISED THIS FROM 0.2 TO 0.3 IN AUGUST 2007
C
         TOTALP=SQRT(PMUX*PMUX+PMUY*PMUY+PMUZ*PMUZ)
         IF(TOTALP.GT.0.3D+00) THEN
            IF(MASWRK) WRITE(IW,9999) LEFP,TOTALP
            PMUX =ZERO
            PMUY =ZERO
            PMUZ =ZERO
            PMUXD=ZERO
            PMUYD=ZERO
            PMUZD=ZERO
            TOTALP=ZERO
         END IF
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 ITER -----
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
C
 9999 FORMAT(1X,'POINT ',I5,' HAS A LARGE INDUCED DIPOLE',
     *          ' TOTAL: ',F9.4,' THUS IS ZEROED OFF')
      END
C*MODULE EFINTB  *DECK EFNUCP
      SUBROUTINE EFNUCP(N,Z,C,ENUCP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION Z(N),C(3,N)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      PARAMETER (MXAO=8192, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /EFPBUF/ POLCHG(10),NBUFMO,LBUFF(MXAO),LBF,NAPOL,IAPOL(10)
      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 /PNUC  / PNX(MXFGPT),PNY(MXFGPT),PNZ(MXFGPT),PNZRF(3)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, GLIM=1.0D-03)
C
C     CALCULATES NUCLEI INTERACTING WITH EXTERNAL POINT
C     POLARIZABILITIES
C
      NEFP = NPTTPT
C
C.... CHANGED TO ADD NUCLEAR DAMPING .... WJS 2/2
C.... TO AGREE WITH MEYER, NUCLEAR TERM SHOULD NOT BE DAMPED!!! 12/22
C.... MAKE UP YOUR MIND, ALREADY .... DRG 8/90
C
      ENFP=ZERO
      DO 340 I=1,NEFP
      PNX(I) = ZERO
      PNY(I) = ZERO
      PNZ(I) = ZERO
C
C     ----- CALCULATE FIELD AT THE POINT ---
C
      IF (NBUFMO .EQ. 0) THEN
      DO 330 J = 1,N
      RR = ZERO
      DO 320 K = 1,3
 320  RR = RR + (EFP(K,I) - C(K,J))**2
      IF(RR.LT.GLIM) GO TO 330
      FACT=ONE
      RR3 = RR * SQRT(RR)
      PNX(I) = PNX(I) + FACT*Z(J) *(EFP(1,I)-C(1,J))/RR3
      PNY(I) = PNY(I) + FACT*Z(J) *(EFP(2,I)-C(2,J))/RR3
      PNZ(I) = PNZ(I) + FACT*Z(J) *(EFP(3,I)-C(3,J))/RR3
 330  CONTINUE
      ELSE
C     ---TAKE INTO ACCOUNT POSSIBLE PRESENCE OF BUFFER IN QM/MM ----
C     INTERACTIONS ONLY WITH THE AB INITIO ZONE NUCLEI AND
C     THE 'BORDER' ATOMS WITH A SPECIAL CHARGE ARE CONSIDERED ONLY
      NATAB=N-LBF
      DO 430 JJ = 1,NATAB+NAPOL
      IF(JJ.GT.NATAB) THEN
         NUMP=JJ-NATAB
         J=IAPOL(NUMP)
         CHARGE=POLCHG(NUMP)
      ELSE
         J=JJ
         CHARGE=Z(J)
      END IF
      RR = ZERO
      DO 420 K = 1,3
 420  RR = RR + (EFP(K,I) - C(K,J))**2
      IF(RR.LT.GLIM) GO TO 430
      FACT=ONE
      RR3 = RR * SQRT(RR)
      PNX(I) = PNX(I) + FACT*CHARGE *(EFP(1,I)-C(1,J))/RR3
      PNY(I) = PNY(I) + FACT*CHARGE *(EFP(2,I)-C(2,J))/RR3
      PNZ(I) = PNZ(I) + FACT*CHARGE *(EFP(3,I)-C(3,J))/RR3
 430  CONTINUE
      END IF
C
C     ----- CALCULATE INTERACTION OF POLARIZABILITY TENSOR WITH FIELD --
C
      POLXX=-EFPOL(1,I)
      POLYY=-EFPOL(2,I)
      POLZZ=-EFPOL(3,I)
      POLXY=-EFPOL(4,I)
      POLXZ=-EFPOL(5,I)
      POLYZ=-EFPOL(6,I)
      POLYX=-EFPOL(7,I)
      POLZX=-EFPOL(8,I)
      POLZY=-EFPOL(9,I)
      ENFP = ENFP+(POLXX*PNX(I)**2+POLYY*PNY(I)**2+POLZZ*PNZ(I)**2
     1           +(POLXY+POLYX)*PNX(I)*PNY(I)
     2           +(POLXZ+POLZX)*PNX(I)*PNZ(I)
     3           +(POLYZ+POLZY)*PNY(I)*PNZ(I))
 340  CONTINUE
      ENUCP = ENFP
      RETURN
      END
C*MODULE EFINTB  *DECK FRGFLD
      SUBROUTINE FRGFLD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL DOMONO,DODIPO,DOQUAD,DOOCTU
      LOGICAL GOPARR,MASWRK,DSKWRK
      LOGICAL IFEWLD,TNFOIL,EWLDPL
C
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION QUA(6,MXFGPT)
C
      CHARACTER*8 FRGNME,POLNAM,DPOLNAM,DNAME,DPSNAM
C
      COMMON /DOMULT/ DOMONO(MXFGPT),DODIPO(MXFGPT),DOQUAD(MXFGPT),
     *                DOOCTU(MXFGPT)
      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 /EFPPBC/ XBOX,YBOX,ZBOX,XCUT,YCUT,ZCUT,RCUT
     *                ,SWR1,SWR2,ISWF
      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)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPAR / LOCNFG,NLOW,NUP
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
C
C
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
      COMMON /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /SCRPAR/ SCROFF
      COMMON /EWALD / IFEWLD,TNFOIL,BETA,ELCNST,EPSLN,CUTOFF,CUTLST,
     *                KMAX,NKVEC,LEVEL,EWLDPL
C
      PARAMETER (ZERO = 0.0D+00, FIVE=5.0D+00, TWO=2.0D+00,ONE=1.0D+00)
      PARAMETER (THREE=3.0D+00, PT5=0.5D+00, ONEPT5=1.5D+00)
      PI=4*ATAN(1.0D+00)
C
C     THIS CODE WAS ORIGINALLY WRITTEN BY DR. W.STEVENS, AND REWRITTEN
C     FOR GAMESS BY JHJ 7/15/91.  IT CALCULATES THE ELECTRIC FIELD DUE
C     TO POINT CHARGES, DIPOLES, AND QUADRUPOLES.
C
      ISCR = INT(ICHGP/2.0D+00)
      CALL VCLR(PPEFX,1,NPTTPT)
      CALL VCLR(PPEFY,1,NPTTPT)
      CALL VCLR(PPEFZ,1,NPTTPT)
      IF(NFRG.LE.1) RETURN
C
      DO 210 I = 1,NMTTPT
        IF(.NOT.DOQUAD(I)) GO TO 210
C
          IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
            XX = QADHF(1,I)
            YY = QADHF(2,I)
            ZZ = QADHF(3,I)
            XY = QADHF(4,I)
            XZ = QADHF(5,I)
            YZ = QADHF(6,I)
          ELSE
            XX = EFQAD(1,I)
            YY = EFQAD(2,I)
            ZZ = EFQAD(3,I)
            XY = EFQAD(4,I)
            XZ = EFQAD(5,I)
            YZ = EFQAD(6,I)
          END IF
C
        DUM = XX + YY + ZZ
        QUA(1,I) = PT5 * (THREE * XX - DUM)
        QUA(2,I) = PT5 * (THREE * YY - DUM)
        QUA(3,I) = PT5 * (THREE * ZZ - DUM)
        QUA(4,I) = ONEPT5 * XY
        QUA(5,I) = ONEPT5 * XZ
        QUA(6,I) = ONEPT5 * YZ
210   CONTINUE
      INDEX = 1
C
      DO I = 1,NLOW-1
       INDEX = INDEX + NPPTS(I)
      ENDDO
C      DO 300 IM = 1,NFRG
      DO 300 IM =NLOW,NUP
         COMX1 = EFCENT(1,IM)
         COMY1 = EFCENT(2,IM)
         COMZ1 = EFCENT(3,IM)
C
        I1 = INDEX
        INDEX = INDEX + NPPTS(IM)
        I2 = I1 + NPPTS(IM) - 1
        DO 250 IPX=I1,I2
          EFX = ZERO
          EFY = ZERO
          EFZ = ZERO
          X1 = EFP(1,IPX)
          Y1 = EFP(2,IPX)
          Z1 = EFP(3,IPX)
          INDEX2 = 1
          DO 200 JM = 1,NFRG
            J1 = INDEX2
            INDEX2 = INDEX2 + NMPTS(JM)
            J2 = J1 + NMPTS(JM) - 1
            IF(IM.EQ.JM) GO TO 200
C
               COMX2 = EFCENT(1,JM)
               COMY2 = EFCENT(2,JM)
               COMZ2 = EFCENT(3,JM)
C
               COMX = COMX1 - COMX2
               COMY = COMY1 - COMY2
               COMZ = COMZ1 - COMZ2
C
               PBCX = ZERO
               PBCY = ZERO
               PBCZ = ZERO
               IF (XBOX .NE. ZERO) PBCX = XBOX * ANINT(COMX/XBOX)
               IF (YBOX .NE. ZERO) PBCY = YBOX * ANINT(COMY/YBOX)
               IF (ZBOX .NE. ZERO) PBCZ = ZBOX * ANINT(COMZ/ZBOX)
C
C   ----  DISTANCE BETWEEN FRAGMENTS
               X12  = ABS(COMX - PBCX)
               Y12  = ABS(COMY - PBCY)
               Z12  = ABS(COMZ - PBCZ)
               R12  = X12*X12+Y12*Y12+Z12*Z12
C
               IF(X12.GT.XCUT.OR.Y12.GT.YCUT.OR.Z12.GT.ZCUT) THEN
                 GOTO 200
               END IF
               R12  = X12*X12+Y12*Y12+Z12*Z12
               IF(R12.GT.RCUT) THEN
                 GOTO 200
               END IF
C
               SWF   = 1.0D+00
               IF(R12.GT.SWR2*SWR2) THEN
                 GOTO 200
               ELSE IF(R12.GT.SWR1*SWR1) THEN
                 IF(ISWF.EQ.1) THEN
                    DUMY1 = 1.0D+00/(SWR2*SWR2 - SWR1*SWR1)
                    DUMY2 = DUMY1*DUMY1
                    DUMY3 = DUMY2*DUMY1
                    DUMY4 = DUMY2*DUMY2
                    DUMY5 = DUMY2*DUMY3
                    VUMY1 = R12 - SWR1*SWR1
                    VUMY2 = VUMY1*VUMY1
                    VUMY3 = VUMY2*VUMY1
                    VUMY4 = VUMY2*VUMY2
                    VUMY5 = VUMY2*VUMY3
                    DSWF  = -30.0D+00*DUMY3*VUMY2
     *                      +60.0D+00*DUMY4*VUMY3
     *                      -30.0D+00*DUMY5*VUMY4
                    DSWF  = 2.0D+00*DSWF
                    SWF   = 1.0D+00 - 10.0D+00*DUMY3*VUMY3
     *                              + 15.0D+00*DUMY4*VUMY4
     *                              -  6.0D+00*DUMY5*VUMY5
                 ELSE IF(ISWF.EQ.2) THEN
                    R12   = SQRT(R12)
                    PIL   = 3.1415926535897932D+00/(SWR2-SWR1)
                    DUMY  = (R12 - SWR1) * PIL
                    SWF   = 0.50D+00 * (COS(DUMY) + 1.0D+00)
                 END IF
               END IF
C
            DO 100 JP = J1,J2
C       WRITE(6,*)'IN FRGFLD,IM,JM,IPX,JP',IM,JM,IPX,JP
              X2 = EFC(1,JP)
              Y2 = EFC(2,JP)
              Z2 = EFC(3,JP)
C
              X = X1 - X2
              Y = Y1 - Y2
              Z = Z1 - Z2
              X = X - PBCX
              Y = Y - PBCY
              Z = Z - PBCZ
C
              R2 = X*X + Y*Y + Z*Z
              R = SQRT(R2)
              R3 = R*R*R
              R5 = R3 * R2
              R7 = R5 * R2
C
C  ---- INFORMATION ABOUT ELECTROSTATIC SCREENING
C  ----  ZERO OUT ALL FRAG-FRAG PENETRATION TERMS
              P1_E = ZERO
              P1_N2 = ZERO
              P1D_E = ZERO
              P2_E = ZERO
              P3_E = ZERO
C
C  ----  AB INITIO-FRAGMENT SCREENING
              IF(ISCR.EQ.1) THEN
C
                IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
                  CF = BMULTHF(JP)
                  EX = AMULTHF(JP)
                ELSE
                  CF = EFBTRM2(JP)
                  EX = EFATRM2(JP)
                END IF
C
              ELSE
                 CF = ZERO
                 EX = ZERO
C
                 IF (EFBTRM2(JP) .NE. ZERO) THEN
                    ALPHAB = EFATRM2(JP)
                 ELSE
                    ALPHAB = 1.0D+20
                 END IF
                 ALPHAA = POLSCR(IPX)
C
C  ----  FRAGMENT-FRAGMENT SCREENING
                 IF(IPLSCR.EQ.1 .AND.ISCR.EQ.3.AND.INDSCR.EQ.0) THEN
C                 IF(1.EQ.0) THEN
C
                    IF(R.LE.SCROFF) THEN
                       AP=ALPHAA*R
                       BP=ALPHAB*R
                       AA2=ALPHAA*ALPHAA
                       AB2=ALPHAB*ALPHAB
                       DIFF=ABS(ALPHAA-ALPHAB)
C
C  ----  NOW CALCULATE ALL FRAG-FRAG CHARGE-PENETRATION TERMS
                    CALL EPENCHDIP(AP,BP,AA2,AB2,DIFF,P1_E,P1_N1,P1_N2)
                    CALL EPENCHQU(AP,BP,AA2,AB2,DIFF,R2,
     *                            P1D_E,P1D_N1,P1D_N2)
                    P2_E =  -P1D_E/R3
                    CALL EPENCHOCT(AP,BP,AA2,AB2,DIFF,R7,
     *                             P3_E,P3_N1,P3_N2)
                    END IF
                 END IF
              END IF
C   ---- END OF ELECTROSTATIC-POLARIZATION SCREENING PART
C
C   ---- EWALD SUMS IF ANY
              F0 = ONE
              F1 = ONE
              F2 = ZERO
              IF (IFEWLD .AND. EWLDPL) THEN
                 F0 = ERFC(BETA*R)
                 EXPBR = EXP(-BETA*BETA*R2)
                 F1 = TWO*BETA*R/SQRT(PI)*EXPBR + F0
                 F2 = 4.0D+00*BETA*BETA*BETA/SQRT(PI)*EXPBR/R2
              END IF    
C
C  FIELD DUE TO TOTAL CHARGE Q (=ELECTRONIC+NUCLEAR)
C
          IF(DOMONO(JP)) THEN
C
            IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
              Q = CHGHF(1,JP) + CHGHF(2,JP)
            ELSE
              Q = EFCHG(1,JP) + EFCHG(2,JP)
              Q_E = EFCHG(1,JP)
              Q_N = EFCHG(2,JP)
            END IF
C
              IF (IFEWLD .AND. EWLDPL) THEN
                 IF(ISCR.EQ.1) THEN
                    EX2=EXP(-EX*R2)
                    SCREEN=1.0D+00-CF*EX2
                    EFX = EFX + Q*X*((F1-ONE)+SCREEN)/R3
     *                     - 2.0D+00*Q*CF*X*EX*EX2/R
                    EFY = EFY + Q*Y*((F1-ONE)+SCREEN)/R3
     *                     - 2.0D+00*Q*CF*Y*EX*EX2/R
                    EFZ = EFZ + Q*Z*((F1-ONE)+SCREEN)/R3
     *                     - 2.0D+00*Q*CF*Z*EX*EX2/R
C
                 ELSE
C     PENCD=0 IF THERE IS NO FRAG-FRAG SCREENING
                    PENCD = Q_E*P1_E + Q_N*P1_N2
                    EFX = EFX + (Q*F1 + PENCD)*X/R3
                    EFY = EFY + (Q*F1 + PENCD)*Y/R3
                    EFZ = EFZ + (Q*F1 + PENCD)*Z/R3
                 END IF
              ELSE 
                 IF(ISCR.EQ.1) THEN
                    EX2=EXP(-EX*R2)
                    SCREEN=1.0D+00-CF*EX2
                    EFX = EFX + SWF*Q*X*SCREEN/R3
     *                     - SWF*2.0D+00*Q*CF*X*EX*EX2/R
                    EFY = EFY + SWF*Q*Y*SCREEN/R3
     *                     - SWF*2.0D+00*Q*CF*Y*EX*EX2/R
                    EFZ = EFZ + SWF*Q*Z*SCREEN/R3
     *                     - SWF*2.0D+00*Q*CF*Z*EX*EX2/R
                 ELSE
                    PENCD = Q_E*P1_E + Q_N*P1_N2
                    EFX = EFX + SWF*(Q + PENCD)*X/R3
                    EFY = EFY + SWF*(Q + PENCD)*Y/R3
                    EFZ = EFZ + SWF*(Q + PENCD)*Z/R3
                 END IF
              END IF
          END IF
C
C  FIELD DUE TO DIPOLE
C
          IF(DODIPO(JP)) THEN
C
            IF (IFEWLD .AND. EWLDPL) THEN
              IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
                DOT1 = DIPHF(1,JP)*X + DIPHF(2,JP)*Y + DIPHF(3,JP)*Z
                DOT = THREE*DOT1/R5
                EFX = EFX + (-DIPHF(1,JP)/R3 + DOT*X)*F1 + DOT1*X*F2
                EFY = EFY + (-DIPHF(2,JP)/R3 + DOT*Y)*F1 + DOT1*Y*F2
                EFZ = EFZ + (-DIPHF(3,JP)/R3 + DOT*Z)*F1 + DOT1*Z*F2
              ELSE
                DOT1 = EFDIP(1,JP)*X + EFDIP(2,JP)*Y + EFDIP(3,JP)*Z
                DOT = THREE*DOT1/R5
                EFX = EFX + (F1+P1_E)*(-EFDIP(1,JP)/R3 + DOT*X) +
     *              (F2+P2_E)*DOT1*X
                EFY = EFY + (F1+P1_E)*(-EFDIP(2,JP)/R3 + DOT*Y) +
     *              (F2+P2_E)*DOT1*Y
                EFZ = EFZ + (F1+P1_E)*(-EFDIP(3,JP)/R3 + DOT*Z) +
     *              (F2+P2_E)*DOT1*Z
              END IF
            ELSE
              IF((IEFP2.EQ.1).AND.(IEFPOL.EQ.1)) THEN
                DOT = DIPHF(1,JP)*X + DIPHF(2,JP)*Y + DIPHF(3,JP)*Z
                DOT = THREE*DOT/R5
                EFX = EFX - SWF*DIPHF(1,JP)/R3 + SWF*DOT*X
                EFY = EFY - SWF*DIPHF(2,JP)/R3 + SWF*DOT*Y
                EFZ = EFZ - SWF*DIPHF(3,JP)/R3 + SWF*DOT*Z
              ELSE
                DOT1 = EFDIP(1,JP)*X + EFDIP(2,JP)*Y + EFDIP(3,JP)*Z
                DOT = THREE*DOT1/R5
                EFX = EFX + SWF*(ONE+P1_E)*(-EFDIP(1,JP)/R3 + DOT*X) +
     *              SWF*P2_E*DOT1*X
                EFY = EFY + SWF*(ONE+P1_E)*(-EFDIP(2,JP)/R3 + DOT*Y) +
     *              SWF*P2_E*DOT1*Y
                EFZ = EFZ + SWF*(ONE+P1_E)*(-EFDIP(3,JP)/R3 + DOT*Z) +
     *              SWF*P2_E*DOT1*Z
              END IF
            END IF
          END IF  
C
C  FIELD DUE TO QUADRUPOLE
C
          IF(DOQUAD(JP)) THEN
C
             IF(IPLSCR.EQ.0 .OR.INDSCR.EQ.1) THEN
C   ---- THIS IS OLD CODE WITHOUT FRAG-FRAG ELECTROSTATIC SCREENING
          C1=FIVE*(QUA(1,JP)*X**2+QUA(2,JP)*Y**2+QUA(3,JP)*Z**2
     *      +TWO*(QUA(4,JP)*X*Y+QUA(5,JP)*X*Z+QUA(6,JP)*Y*Z))/R7
     *      -    (QUA(1,JP)+QUA(2,JP)+QUA(3,JP))/R5
              EFX = EFX + SWF*X*(C1-TWO*QUA(1,JP)/R5)
     *                  - SWF*TWO*(QUA(4,JP)*Y+QUA(5,JP)*Z)/R5
              EFY = EFY + SWF*Y*(C1-TWO*QUA(2,JP)/R5)
     *                  - SWF*TWO*(QUA(4,JP)*X+QUA(6,JP)*Z)/R5
              EFZ = EFZ + SWF*Z*(C1-TWO*QUA(3,JP)/R5)
     *                  - SWF*TWO*(QUA(5,JP)*X+QUA(6,JP)*Y)/R5
             ELSE
C   ---- WITH FRAG-FRAG ELECTROSTATIC-POLARIZATION SCREENING
C
               TXX = (THREE*X*X-R2)/R5
               TYY = (THREE*Y*Y-R2)/R5
               TZZ = (THREE*Z*Z-R2)/R5
               TXY = THREE*X*Y/R5
               TXZ = THREE*X*Z/R5
               TYZ = THREE*Y*Z/R5
C
           DQXX = (ONE+P1_E)*THREE*(-FIVE*X*X+R2)/R7+P1D_E*TXX
     *                     + P3_E*X*X
           DQYY = (ONE+P1_E)*THREE*(-FIVE*Y*Y+R2)/R7+P1D_E*TYY
     *                     + P3_E*Y*Y
           DQZZ = (ONE+P1_E)*THREE*(-FIVE*Z*Z+R2)/R7+P1D_E*TZZ
     *                     + P3_E*Z*Z
           DQXY = (ONE+P1_E)*THREE*(-FIVE*X*Y)/R7+P1D_E*TXY
     *                     + P3_E*X*Y
           DQXZ = (ONE+P1_E)*THREE*(-FIVE*X*Z)/R7 +P1D_E*TXZ
     *                     + P3_E*X*Z
           DQYZ = (ONE+P1_E)*THREE*(-FIVE*Y*Z)/R7 +P1D_E*TYZ
     *                     + P3_E*Y*Z
C
                  DQ1 = (QUA(1,JP)*DQXX
     *                  + QUA(2,JP)*DQYY + QUA(3,JP)*DQZZ
     *                  + TWO*(QUA(4,JP)*DQXY
     *                  + QUA(5,JP)*DQXZ +QUA(6,JP)*DQYZ))/THREE
C
                  CHPDQ = TWO*((ONE+P1_E)*THREE/R5+P2_E)
C
                  QUAX = QUA(1,JP)*X+QUA(4,JP)*Y+QUA(5,JP)*Z
                  QUAY = QUA(4,JP)*X+QUA(2,JP)*Y+QUA(6,JP)*Z
                  QUAZ = QUA(5,JP)*X+QUA(6,JP)*Y+QUA(3,JP)*Z
C
                         TIMX = CHPDQ*QUAX/THREE
                         TIMY = CHPDQ*QUAY/THREE
                         TIMZ = CHPDQ*QUAZ/THREE
C
C   ---- THESE ARE FINAL EXPRESSIONS FOR SCREENED FIELD
                         EFX = EFX - SWF*(TIMX + DQ1*X)
                         EFY = EFY - SWF*(TIMY + DQ1*Y)
                         EFZ = EFZ - SWF*(TIMZ + DQ1*Z)
            END IF
C
        END IF
C
C  OCTUPOLES ARE NOT INCLUDED
C
  100       CONTINUE
  200     CONTINUE
          PPEFX(IPX) = PPEFX(IPX) + EFX
          PPEFY(IPX) = PPEFY(IPX) + EFY
          PPEFZ(IPX) = PPEFZ(IPX) + EFZ
  250   CONTINUE
  300 CONTINUE
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(675,PPEFX,NPTTPT)
         CALL DDI_GSUMF(676,PPEFY,NPTTPT)
         CALL DDI_GSUMF(677,PPEFZ,NPTTPT)
      END IF
C
      RETURN
      END
C*MODULE EFINTB  *DECK POLINT
      SUBROUTINE POLINT(PEX,PEY,PEZ,NUM2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      LOGICAL IANDJ,NORM,DOUBLE,OUT
      LOGICAL GOPARR,MASWRK,DSKWRK
      LOGICAL DSKSAV
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000,
     *           MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION PEX(NUM2),PEY(NUM2),PEZ(NUM2)
C
      DIMENSION XIN(128),YIN(128),ZIN(128),
     1          PLX(100),PLY(100),PLZ(100)
      DIMENSION DIJ(100),FIJ(100)
      DIMENSION IX(20),IY(20),IZ(20),JX(20),JY(20),JZ(20)
      DIMENSION IJX(100),IJY(100),IJZ(100)
C
      COMMON /CSSTV / CX,CY,CZ
      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 /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),
     1                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     2                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ROOT  / XX,U(13),W(13),NROOTS
      COMMON /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 JX / 0, 1, 0, 0, 2, 0, 0, 1, 1, 0,
     1          3, 0, 0, 2, 2, 1, 0, 1, 0, 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 JY / 0, 0, 1, 0, 0, 2, 0, 1, 0, 1,
     1          0, 3, 0, 1, 0, 2, 2, 0, 1, 1/
      DATA IY / 1, 1, 5, 1, 1, 9, 1, 5, 1, 5,
     1          1,13, 1, 5, 1, 9, 9, 1, 5, 5/
      DATA JZ / 0, 0, 0, 1, 0, 0, 2, 0, 1, 1,
     1          0, 0, 3, 0, 1, 0, 1, 2, 2, 1/
      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 ZERO,ONE/0.0D+00,1.0D+00/
      DATA PI212 /1.1283791670955D+00/
      DATA SQRT3 /1.73205080756888D+00/
      DATA SQRT5 /2.23606797749979D+00/
      DATA RLN10 /2.30258D+00/
      DATA DBUGME  /8HINT1    /
C
C------- THIS VERSION HAS THE METHOD=1 STUFF OF HONDO CUT OUT.
C
      OUT = NPRINT.EQ.3 .OR. EXETYP.EQ.DBUGME
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
C
      NEFP = NPTTPT
      TOL=RLN10*ITOL
      NORM=NORMF.NE.1.OR.NORMP.NE.1
C
C        THE ELECTRIC FIELD INTEGRAL IS TO BE A DISTRIBUTED FILE
C
      DSKSAV = DSKWRK
      DSKWRK = .TRUE.
      NFT25 = 25
      CALL SEQREW(NFT25)
C
      DO 9200 IC=1,NEFP
C
C     ----- GO PARALLEL! -----
C
       IF (GOPARR) THEN
        IPCOUNT = IPCOUNT + 1
        IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 9199
       END IF
C
      DO 41 I=1,NUM2
      PEX(I) = ZERO
      PEY(I) = ZERO
      PEZ(I) = ZERO
  41  CONTINUE
C     IF(MASWRK) WRITE(6,*)'FOR POLARIZABLE POINT ',IC
      CX=EFP(1,IC)
      CY=EFP(2,IC)
      CZ=EFP(3,IC)
C
C     ----- ISHELL
C
      DO 9000 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 8000 JJ=1,II
C
C     ----- GO PARALLEL! -----
C
C            IF (GOPARR) THEN
C               IPCOUNT = IPCOUNT + 1
C               IF (MOD(IPCOUNT,NPROC).NE.0) GO TO 8000
C            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+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
      PLX(I) = ZERO
      PLY(I) = ZERO
  60  PLZ(I) = ZERO
C
C     ----- I PRIMITIVE
C
      JGMAX=J2
      DO 7000 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     ----- J PRIMITIVE
C
      IF(IANDJ) JGMAX=IG
      DO 6000 JG=J1,JGMAX
      AJ=EX(JG)
      AA=AI+AJ
      AA1=ONE/AA
      DUM=AJ*ARRI*AA1
      IF(DUM.GT.TOL) GO TO 6000
      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 310 I=MINI,MAXI
      GO TO ( 70, 80,180,180, 90,180,180,100,180,180,
     1       110,180,180,120,180,180,180,180,180,130),I
   70 DUM1=CSI*FAC
      GO TO 180
   80 DUM1=CPI*FAC
      GO TO 180
   90 DUM1=CDI*FAC
      GO TO 180
  100 IF(NORM) DUM1=DUM1*SQRT3
      GO TO 180
  110 DUM1=CFI*FAC
      GO TO 180
  120 IF(NORM) DUM1=DUM1*SQRT5
      GO TO 180
  130 IF(NORM) DUM1=DUM1*SQRT3
  180 IF(IANDJ) MAX=I
      DO 310 J=MINJ,MAX
      GO TO (190,200,300,300,210,300,300,220,300,300,
     1       230,300,300,240,300,300,300,300,300,250),J
  190 DUM2=DUM1*CSJ
      IF(.NOT.DOUBLE) GO TO 300
      IF(I.GT.1) GO TO 195
      DUM2=DUM2+DUM2
      GO TO 300
  195 DUM2=DUM2+CSI*CPJ*FAC
      GO TO 300
  200 DUM2=DUM1*CPJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  210 DUM2=DUM1*CDJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  220 IF(NORM) DUM2=DUM2*SQRT3
      GO TO 300
  230 DUM2=DUM1*CFJ
      IF(DOUBLE) DUM2=DUM2+DUM2
      GO TO 300
  240 IF(NORM) DUM2=DUM2*SQRT5
      GO TO 300
  250 IF(NORM) DUM2=DUM2*SQRT3
  300 NN=NN+1
  310 DIJ(NN)=DUM2
C
      AAX=AA*AX
      AAY=AA*AY
      AAZ=AA*AZ
      DUM=PI212*AA1
      DUM=DUM+DUM
      DO 800 I=1,IJ
  800 FIJ(I)=DIJ(I)*DUM
C
      XX=AA*((AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2)
      IF(NROOTS.LE.3) CALL RT123
      IF(NROOTS.EQ.4) CALL ROOT4
      MM=0
      DO 830 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 820 I=1,LIT
      IN=IN+4
      NI=I
      DO 820 J=1,LJT
      JN=IN+J
      NJ=J
      CALL STVINT
      XIN(JN)=XINT
      YIN(JN)=YINT
      ZIN(JN)=ZINT*WW
      CALL POLXYZ
      XIN(JN+64)=XINT
      YIN(JN+64)=YINT
      ZIN(JN+64)=ZINT*WW
  820 CONTINUE
  830 MM=MM+16
      DO 850 I=1,IJ
      NX=IJX(I)
      NY=IJY(I)
      NZ=IJZ(I)
      DUMX=ZERO
      DUMY=ZERO
      DUMZ=ZERO
      MM=0
      DO 840 K=1,NROOTS
      DUMX=DUMX+XIN(NX+MM+64)*YIN(NY+MM   )*ZIN(NZ+MM   )
      DUMY=DUMY+XIN(NX+MM   )*YIN(NY+MM+64)*ZIN(NZ+MM   )
      DUMZ=DUMZ+XIN(NX+MM   )*YIN(NY+MM   )*ZIN(NZ+MM+64)
  840 MM=MM+16
      DUM=FIJ(I)
      PLX(I) = PLX(I) + DUM*DUMX
      PLY(I) = PLY(I) + DUM*DUMY
      PLZ(I) = PLZ(I) + DUM*DUMZ
  850 CONTINUE
C
 6000 CONTINUE
 7000 CONTINUE
C
C
      MAX=MAXJ
      NN=0
      DO 7500 I=MINI,MAXI
      LI=LOCI+I
      IN = (LI*(LI-1))/2
      IF(IANDJ) MAX=I
      DO 7500 J=MINJ,MAX
      LJ=LOCJ+J
      JN=LJ+IN
      NN=NN+1
      PEX(JN)=PLX(NN)
      PEY(JN)=PLY(NN)
      PEZ(JN)=PLZ(NN)
 7500 CONTINUE
C     ----- END PARALLEL
 8000 CONTINUE
 9000 CONTINUE
C
C  COLLECT MATRICES ON DISK FOR TWO ELECTRON OPERATOR CONSTRUCTION
C
C      IF (GOPARR) CALL DDI_GSUMF(924,PEX,NUM2)
C      IF (GOPARR) CALL DDI_GSUMF(925,PEY,NUM2)
C      IF (GOPARR) CALL DDI_GSUMF(926,PEZ,NUM2)
C
      CALL SQWRIT(NFT25,PEX,NUM2)
      CALL SQWRIT(NFT25,PEY,NUM2)
      CALL SQWRIT(NFT25,PEZ,NUM2)
C
      IF(OUT) THEN
        IF (MASWRK) WRITE(IW,9500)
        CALL PRTRIL(PEX,NUM)
        CALL PRTRIL(PEY,NUM)
        CALL PRTRIL(PEZ,NUM)
      END IF
C
 9199 CONTINUE
 9200 CONTINUE
C
      DSKWRK = DSKSAV
C
      RETURN
 9500 FORMAT(/1X,'X, Y, AND Z COMPONENTS OF THE INDUCED DIPOLE',
     *       ' ONE ELECTRON INTEGRALS')
      END
C*MODULE EFINTB  *DECK REPINT
      SUBROUTINE REPINT(QQ,GG)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 REPNAM
C
      LOGICAL IANDJ,NORM,DOUBLE
      LOGICAL OUT,GOPARR,MASWRK,DSKWRK
C
      DIMENSION QQ(*),GG(*)
C
      DIMENSION RINT(100)
      DIMENSION DIJ(100),XIN(432),YIN(432),ZIN(432)
      DIMENSION GIJ(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,
     *           MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /CSSTV / CX,CY,CZ
      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 /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 /REPPAR/ CREP(3,MXFGPT),CLPR(4*MXFGPT),ZLPR(4*MXFGPT),
     *                NLPR(4*MXFGPT),KFR(MXFGPT),KLR(MXFGPT),
     *                REPNAM(MXFGPT)
      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 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 PI212 /1.1283791670955D+00/
      DATA SQRT3 /1.73205080756888D+00/
      DATA SQRT5 /2.23606797749979D+00/
      DATA ZERO,ONE/0.0D+00,1.0D+00/
      DATA RLN10 /2.30258D+00/
      DATA DBUGME/8HINT1    /
C
C -----CHARGE-REPULSIVE POTENTIAL INTEGRALS FOR POWERS 0 AND -1 OF R.
C
      NREP = NRTTPT
      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 40 I=1,NUM2
         QQ(I) = ZERO
  40  CONTINUE
C
C     ----- INTIALIZE PARALLEL -----
C
      IPCOUNT = ME - 1
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 + 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  RINT(I) = ZERO
C
C     ----- I PRIMITIVE
C
      JGMAX = J2
      DO 520 IG = I1,I2
C     WRITE(6,*)'IG=',IG
      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     ----- J PRIMITIVE
C
      IF (IANDJ) JGMAX = IG
      DO 500 JG = J1,JGMAX
C     WRITE(6,*)'JG=',JG
      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
      AAX = AA*AX
      AAY = AA*AY
      AAZ = AA*AZ
C
      DO 480 IC = 1,NREP
      CX   =CREP(1,IC)
      CY   =CREP(2,IC)
      CZ   =CREP(3,IC)
      PCSQ = (AX-CX)**2+(AY-CY)**2+(AZ-CZ)**2
      DO 480 LTERM=KFR(IC),KLR(IC)
      ALFA = ZLPR(LTERM)
      BETA = CLPR(LTERM)
      PREI = EXP(-AA*ALFA*PCSQ/(AA+ALFA))
      IF(NLPR(LTERM).EQ.2) THEN
C      R TO THE ZERO POWER.
       TT = ONE/(AA+ALFA)
       T  =  SQRT(TT)
       X0 = (AAX+ALFA*CX)*TT
       Y0 = (AAY+ALFA*CY)*TT
       Z0 = (AAZ+ALFA*CZ)*TT
       IN = -4
       DO 206 I = 1,LIT
       IN = IN+4
       NI = I
       DO 206 J = 1,LJT
       JN = IN+J
       NJ = J
       CALL STVINT
       XIN(JN   ) = XINT*T
       YIN(JN   ) = YINT*T
       ZIN(JN   ) = ZINT*T
  206  CONTINUE
       DO 266 I = 1,IJ
       NX    = IJX(I)
       NY    = IJY(I)
       NZ    = IJZ(I)
       RINT(I)=RINT(I)+DIJ(I)*PREI*BETA*XIN(NX)*YIN(NY)*ZIN(NZ)
 266   CONTINUE
      ELSE
C      ONE OVER R.
       DUM = PI212/(AA+ALFA)
       DO 385 I = 1,IJ
       GIJ(I) = DIJ(I)*DUM
  385  CONTINUE
       XX=AA*AA*PCSQ/(AA+ALFA)
       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 425 K = 1,NROOTS
       UU = (ALFA+AA)*U(K)
       WW = W(K)
       TT = ONE/(AA+UU+ALFA)
       T  =  SQRT(TT)
       X0 = (AAX+(UU+ALFA)*CX)*TT
       Y0 = (AAY+(UU+ALFA)*CY)*TT
       Z0 = (AAZ+(UU+ALFA)*CZ)*TT
       IN = -4+MM
       DO 405 I = 1,LIT
       IN = IN+4
       NI = I
       DO 405 J = 1,LJT
       JN = IN+J
       NJ = J
       CALL STVINT
       XIN(JN   ) = XINT
       YIN(JN   ) = YINT
       ZIN(JN   ) = ZINT*WW
  405  CONTINUE
  425  MM = MM+16
       DO 465 I = 1,IJ
       NX    = IJX(I)
       NY    = IJY(I)
       NZ    = IJZ(I)
       DUM = ZERO
       MM    = 0
       DO 445 K = 1,NROOTS
       DUM= DUM+XIN(NX+MM)*YIN(NY+MM)*ZIN(NZ+MM)
  445  MM  = MM+16
       RINT(I)=RINT(I)+GIJ(I)*PREI*BETA*DUM
C       WRITE(6,*)'AFTER RINT, I,RINT(I)',I,RINT(I)
 465   CONTINUE
      END IF
C
 480  CONTINUE
C-------- END OF REP CENTERS LOOP
 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) + RINT(NN)
  550 CONTINUE
C
C     ----- END PARALLEL
C
  570 CONTINUE
C
C     ----- END SHELLS -----
C
  580 CONTINUE
  600 CONTINUE
C
C     ----- SUM UP PARTIAL CONTRIBUTIONS IF PARALLEL -----
C
      IF (GOPARR) CALL DDI_GSUMF(928,QQ,NUM2)
C
      CALL DAREAD(IDAF,IODA,GG,NUM2,11,0)
      CALL VADD(GG,1,QQ,1,GG,1,NUM2)
      CALL DAWRIT(IDAF,IODA,GG,NUM2,11,NAV)
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 REPINT.  CALL A PROGRAMMER/QUANTUM CHEMIST.')
 9128 FORMAT(/,' ...... END OF REP INTEGRALS....')
 9380 FORMAT(10X,14(1H-),/,10X,14H   REP  MATRIX,
     1 /,10X,14(1H-))
      END
C*MODULE EFINTB  *DECK FLDIND
      SUBROUTINE FLDIND(EFADD)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 POLNAM,DPOLNAM,DNAME,DPSNAM
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXIFRQ=12, MXATM=2000,
     *           MXPT=2000, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /EFPAR / LOCNFG,NLOW,NUP
      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 /EFPPBC/ XBOX,YBOX,ZBOX,XCUT,YCUT,ZCUT,RCUT
     *                ,SWR1,SWR2,ISWF
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      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 /POL2  / CHGHF(2,MXFGPT),DIPHF(3,MXFGPT),QADHF(6,MXFGPT),
     *                AMULTHF(MXFGPT),BMULTHF(MXFGPT),
     *                AMULT2HF(MXFGPT),BMULT2HF(MXFGPT),
     *                EFPOL2(9,MXFRG),DIND2(3,MXFGPT),DINDD2(3,MXFGPT),
     *                DCORD(3,MXPT),CDIS(3,MXFGPT),
     *                CLPD(4*MXFGPT),ZLPD(4*MXFGPT),
     *                NLPD(4*MXFGPT),KFD(MXFGPT),KLD(MXFGPT),
     *                DNAME(MXPT),DPSNAM(MXFGPT)
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCRPAR/ SCROFF
C
C
      DIMENSION EFADD(3,NPTTPT)
C
      DATA THREE/3.0D+00/, ZERO/0.0D+00/, ONE/1.0D+00/
C
C     COMPUTE THE FIELDS ON THE POL-TENSORS DUE TO INDUCED DIPOLES
C
C     HUI LI, IOWA STATE UNIVERSITY, DEC 29, 2005
C
C
      CALL VCLR(EFADD,1,3*NPTTPT)
      IF(NAT.LE.0) RETURN
C
      IC1=0
      IF (GOPARR) THEN
         DO I=1,NLOW-1
            IC1 = IC1 + NPPTS(I)
         ENDDO
      ENDIF
C
C
      DO 400 IM = NLOW,NUP
         COMX1 = EFCENT(1,IM)
         COMY1 = EFCENT(2,IM)
         COMZ1 = EFCENT(3,IM)
C
         DO 390 IPX = 1,NPPTS(IM)
            IC1=IC1+1
            X1 = EFP(1,IC1)
            Y1 = EFP(2,IC1)
            Z1 = EFP(3,IC1)
C
            FXIND  = ZERO
            FYIND  = ZERO
            FZIND  = ZERO
C
            IC2 = 0
            DO 380 JM = 1,NFRG
C
               IF (IM.EQ.JM) THEN
                  IC2=IC2+NPPTS(JM)
                  GO TO 380
               END IF
C
               COMX2 = EFCENT(1,JM)
               COMY2 = EFCENT(2,JM)
               COMZ2 = EFCENT(3,JM)
C
               COMX = COMX1 - COMX2
               COMY = COMY1 - COMY2
               COMZ = COMZ1 - COMZ2
C
               PBCX = ZERO
               PBCY = ZERO
               PBCZ = ZERO
               IF (XBOX .NE. ZERO) PBCX = XBOX * ANINT(COMX/XBOX)
               IF (YBOX .NE. ZERO) PBCY = YBOX * ANINT(COMY/YBOX)
               IF (ZBOX .NE. ZERO) PBCZ = ZBOX * ANINT(COMZ/ZBOX)
C
C   ----  DISTANCE BETWEEN FRAGMENTS
               X12  = ABS(COMX - PBCX)
               Y12  = ABS(COMY - PBCY)
               Z12  = ABS(COMZ - PBCZ)
               R12  = X12*X12+Y12*Y12+Z12*Z12
C
               IF(X12.GT.XCUT.OR.Y12.GT.YCUT.OR.Z12.GT.ZCUT)THEN
                 IC2=IC2+NPPTS(JM)
                 GOTO 380
               END IF
               R12  = X12*X12+Y12*Y12+Z12*Z12
               IF(R12.GT.RCUT)THEN
                 IC2=IC2+NPPTS(JM)
                 GOTO 380
               END IF
C
               SWF   = 1.0D+00
               IF(R12.GT.SWR2*SWR2) THEN
                 IC2=IC2+NPPTS(JM)
                 GOTO 380
               ELSE IF(R12.GT.SWR1*SWR1) THEN
                 IF(ISWF.EQ.1) THEN
                    DUMY1 = 1.0D+00/(SWR2*SWR2 - SWR1*SWR1)
                    DUMY2 = DUMY1*DUMY1
                    DUMY3 = DUMY2*DUMY1
                    DUMY4 = DUMY2*DUMY2
                    DUMY5 = DUMY2*DUMY3
                    VUMY1 = R12 - SWR1*SWR1
                    VUMY2 = VUMY1*VUMY1
                    VUMY3 = VUMY2*VUMY1
                    VUMY4 = VUMY2*VUMY2
                    VUMY5 = VUMY2*VUMY3
                    DSWF  = -30.0D+00*DUMY3*VUMY2
     *                      +60.0D+00*DUMY4*VUMY3
     *                      -30.0D+00*DUMY5*VUMY4
                    DSWF  = 2.0D+00*DSWF
                    SWF   = 1.0D+00 - 10.0D+00*DUMY3*VUMY3
     *                              + 15.0D+00*DUMY4*VUMY4
     *                              -  6.0D+00*DUMY5*VUMY5
                 ELSE IF(ISWF.EQ.2) THEN
                    R12   = SQRT(R12)
                    PIL   = 3.1415926535897932D+00/(SWR2-SWR1)
                    DUMY  = (R12 - SWR1) * PIL
                    SWF   = 0.50D+00 * (COS(DUMY) + 1.0D+00)
                 END IF
               END IF
C
               DO 370 JP = 1,NPPTS(JM)
                  IC2=IC2+1
                  X2 = EFP(1,IC2)
                  Y2 = EFP(2,IC2)
                  Z2 = EFP(3,IC2)
                  X = X1 - X2
                  Y = Y1 - Y2
                  Z = Z1 - Z2
                  X = X - PBCX
                  Y = Y - PBCY
                  Z = Z - PBCZ
C
                  R2 = X * X + Y * Y + Z * Z
                  R = SQRT(R2)
                  R3 = R * R * R
                  R5 = R * R * R3
C
C  ----  LVS: POLARIZATION FRAG-FRAG SCREENING
C  ----  ZERO OUT ALL PENETRATION TERMS FIRST
                     P1 = ZERO
                     P1D = ZERO
                     P2 = ZERO
C  ----  IF DO FRAGMENT-FRAGMENT SCREENING
                     IF(IPLSCR.EQ.1) THEN
C                     IF(1.EQ.0) THEN
                       ALPHAB = POLSCR(JP)
                       ALPHAA = POLSCR(IPX)
C
                       IF(R.LE.SCROFF) THEN
                         AP=ALPHAA*R
                         BP=ALPHAB*R
                         AA2=ALPHAA*ALPHAA
                         AB2=ALPHAB*ALPHAB
                         DIFF=ABS(ALPHAA-ALPHAB)
C
C  ----  NOW CALCULATE ALL FRAG-FRAG CHARGE-PENETRATION TERMS
                         CALL EPENCHDIP(AP,BP,AA2,AB2,DIFF,
     *                                  P1,P1_N1,P1_N2)
                         CALL EPENCHQU(AP,BP,AA2,AB2,DIFF,R2,
     *                            P1D,P1D_N1,P1D_N2)
                         P2 =  -P1D/R3
                       END IF
                     END IF
C   ---- END OF POLARIZATION SCREENING PART
C
                  IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
                     DOT = DIND2(1,IC2)*X + DIND2(2,IC2)*Y
     *                   + DIND2(3,IC2)*Z
                     DOT = THREE * DOT / R5
                     FXIND=FXIND-SWF*DIND2(1,IC2)/R3+SWF*DOT*X
                     FYIND=FYIND-SWF*DIND2(2,IC2)/R3+SWF*DOT*Y
                     FZIND=FZIND-SWF*DIND2(3,IC2)/R3+SWF*DOT*Z
                  ELSE
C                     DOT = DIND(1,IC2)*X + DIND(2,IC2)*Y
C     *                   + DIND(3,IC2)*Z
C                     DOT = THREE * DOT / R5
C                     FXIND=FXIND-SWF*DIND(1,IC2)/R3+SWF*DOT*X
C                     FYIND=FYIND-SWF*DIND(2,IC2)/R3+SWF*DOT*Y
C                     FZIND=FZIND-SWF*DIND(3,IC2)/R3+SWF*DOT*Z
C
                     DOT1 = DIND(1,IC2)*X
     *                    + DIND(2,IC2)*Y + DIND(3,IC2)*Z
                     DOT = THREE * DOT1 / R5
                     FXIND=FXIND + SWF*((ONE+P1)*
     *                          (-DIND(1,IC2)/R3 + DOT*X)+P2*DOT1*X)
                     FYIND=FYIND + SWF*((ONE+P1)*
     *                          (-DIND(2,IC2)/R3 + DOT*Y)+P2*DOT1*Y)
                     FZIND=FZIND + SWF*((ONE+P1)*
     *                          (-DIND(3,IC2)/R3 + DOT*Z)+P2*DOT1*Z)
                  END IF
C
370            CONTINUE
380         CONTINUE
            EFADD(1,IC1)   = FXIND
            EFADD(2,IC1)   = FYIND
            EFADD(3,IC1)   = FZIND
390      CONTINUE
400   CONTINUE
C
      IF(GOPARR) THEN
         CALL DDI_GSUMF(770,EFADD,3*NPTTPT)
      END IF
C
      RETURN
      END
