C  4 MAR 08 - MWS - RAISE EFP PROCESSOR LIMIT FROM 64 TO 1024
C 20 AUG 07 - DDK,LVS - SMALL CHANGE TO GEXREP CALL, PBC-EFP CLEANED
C 24 MAR 07 - LVS - POLARIZATION SCREENING ADDED, EL.DERIV. MOVED OUT
C 22 DEC 06 - TN  - FMO/EFP GRADIENT
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 22 FEB 06 - TN  - EFP/FMO INTERFACING
C
C*MODULE EFPFMO  *DECK EFFMO
      SUBROUTINE EFFMO
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      DATA OPTFMO/8HOPTFMO  /
C
C     IF BOTH $FMO AND $EFRAG EXIST,
C     EFP/FMO CALCULATION WILL BE CARRIED OUT.
C
      IEOF    = 0
      IEFPFMO  = 0
      ISWENUC = 1
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $EFRAG ',IEOF)
      IF(IEOF.NE.0) THEN
        RETURN
      ELSE
        IEOF = 0
        CALL SEQREW(IR)
        CALL FNDGRP(IR,' $FMO   ',IEOF)
        IF (IEOF.NE.0) THEN
          RETURN
        ELSE
          IEFPFMO  = 1
          ISWENUC = 0
        END IF
      END IF
C
      IF (RUNTYP.EQ.OPTFMO) THEN
        WRITE(6,9010)
        CALL ABRT
      END IF
C
      IF (MASWRK) WRITE(IW,9000)
      RETURN
 9000 FORMAT(/3X,44(1H-)/
     *   3X,'   ***** FMO/EFP CALCULATION STARTS *****   '/
     *   3X,44(1H-))
 9010 FORMAT(/3X, 'EFP/FMO TO USE FMO GEOMETRY OPTIMIZATION ENGINE '/
     *            'IS UNDER CONSTRUCTION')
      END
C*MODULE EFPFMO  *DECK DIPITFMO
      SUBROUTINE DIPITFMO(EFLD,EFADD,DINDFMO,DIND2FMO,
     *                    DINDDFMO,DINDD2FMO,
     *                    DINDC,DINDDC,DIND2C,DINDD2C,
     *                    DIND_W,DIND_C,DIND2_W,DIND2_C,N,
     *                    DIPNEW,DIPNWD,NXFRG)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      CHARACTER*8 POLNAM,DPOLNAM,DNAME,DPSNAM
C
      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)
C        NEXT ARE A STATIC ALLOCATION,
C        NOTE INITIALIZATION AND -SAVE- STATEMENT BELOW
      DIMENSION DIPLD(MXFRG),DIPDD(MXFRG)
C
      DIMENSION DINDFMO(3,N), DIND2FMO(3,N)
      DIMENSION DINDDFMO(3,N),DINDD2FMO(3,N)
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 /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
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)
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 /SCRPAR/ SCROFF
      INTEGER DDI_NP, DDI_ME
C
      PARAMETER (ZERO=0.00D+00, THREE=3.00D+00)
      PARAMETER (ONE=1.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/
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
C     SAME RESULTS
C     ALWAYS USE CRITERION TO 1.0D-10
C
      ITLIM= 30
      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
      IF((IEFP2.EQ.1).AND.(IEFPOL2.EQ.1)) THEN
         DO I = 1,NIN2
            DIND2_C(1,I) = DIND2FMO(1,KP)
            DIND2_C(2,I) = DIND2FMO(2,KP)
            DIND2_C(3,I) = DIND2FMO(3,KP)
            DIND2_C(4,I) = DINDD2FMO(1,KP)
            DIND2_C(5,I) = DINDD2FMO(2,KP)
            DIND2_C(6,I) = DINDD2FMO(3,KP)
            KP = KP + 1
         ENDDO
      ELSE
         DO I = 1,NIN2
            DIND_C(1,I) = DINDFMO(1,KP)
            DIND_C(2,I) = DINDFMO(2,KP)
            DIND_C(3,I) = DINDFMO(3,KP)
            DIND_C(4,I) = DINDDFMO(1,KP)
            DIND_C(5,I) = DINDDFMO(2,KP)
            DIND_C(6,I) = DINDDFMO(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
                  EFX = EFLD(1,IPX)
                  EFY = EFLD(2,IPX)
                  EFZ = EFLD(3,IPX)
                  EFXD = EFLD(1,IPX)
                  EFYD = EFLD(2,IPX)
                  EFZD = EFLD(3,IPX)
               ELSE
                  EFX =ZERO
                  EFY =ZERO
                  EFZ =ZERO
                  EFXD =ZERO
                  EFYD =ZERO
                  EFZD =ZERO
               END IF
               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)
                         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 = DIND2FMO(1,JP)*X
     *                         + DIND2FMO(2,JP)*Y
     *                         + DIND2FMO(3,JP)*Z
                           DOT = THREE * DOT / R5
                           EFX = EFX - SWF*DIND2FMO(1,JP)/R3 + SWF*DOT*X
                           EFY = EFY - SWF*DIND2FMO(2,JP)/R3 + SWF*DOT*Y
                           EFZ = EFZ - SWF*DIND2FMO(3,JP)/R3 + SWF*DOT*Z
C
                           DOTD = DINDD2FMO(1,JP)*X
     *                          + DINDD2FMO(2,JP)*Y
     *                          + DINDD2FMO(3,JP)*Z
                           DOTD = THREE * DOTD / R5
                           EFXD = EFXD-SWF*DINDD2FMO(1,JP)/R3+SWF*DOTD*X
                           EFYD = EFYD-SWF*DINDD2FMO(2,JP)/R3+SWF*DOTD*Y
                           EFZD = EFZD-SWF*DINDD2FMO(3,JP)/R3+SWF*DOTD*Z
                        END IF
                     ELSE
C
C  ----  LVS: SCREENING FOR FRAGMENT-FRAGMENT PART ADDED
                        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 = DINDFMO(1,JP)*X
     *                          + DINDFMO(2,JP)*Y
     *                          + DINDFMO(3,JP)*Z
                           DOT = THREE * DOT1 / R5
                           EFX = EFX+ SWF*((ONE+P1)*
     *                          (-DINDFMO(1,JP)/R3 + DOT*X)+P2*DOT1*X)
                           EFY = EFY+ SWF*((ONE+P1)*
     *                          (-DINDFMO(2,JP)/R3 + DOT*Y)+P2*DOT1*Y)
                           EFZ = EFZ+ SWF*((ONE+P1)*
     *                          (-DINDFMO(3,JP)/R3 + DOT*Z)+P2*DOT1*Z)
C
                           DOTD1 = DINDDFMO(1,JP)*X
     *                           + DINDDFMO(2,JP)*Y
     *                           + DINDDFMO(3,JP)*Z
                           DOTD = THREE * DOTD1 / R5
                           EFXD = EFXD+ SWF*((ONE+P1)*
     *                        (-DINDDFMO(1,JP)/R3 + DOTD*X)+P2*DOTD1*X)
                           EFYD = EFYD+ SWF*((ONE+P1)*
     *                        (-DINDDFMO(2,JP)/R3 + DOTD*Y)+P2*DOTD1*Y)
                           EFZD = EFZD+ SWF*((ONE+P1)*
     *                        (-DINDDFMO(3,JP)/R3 + DOTD*Z)+P2*DOTD1*Z)
                        END IF
                     END IF
  300             CONTINUE
  309             CONTINUE
                  INDEX2 = INDEX2 + NPPTS(JM)
  310          CONTINUE
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
                        DIND2FMO(1,IPX) = EFX*EFPOL2(1,IM)
     *                                  + EFY*EFPOL2(4,IM)
     *                                  + EFZ*EFPOL2(5,IM)
                        DIND2FMO(2,IPX) = EFX*EFPOL2(7,IM)
     *                                  + EFY*EFPOL2(2,IM)
     *                                  + EFZ*EFPOL2(6,IM)
                        DIND2FMO(3,IPX) = EFX*EFPOL2(8,IM)
     *                                  + EFY*EFPOL2(9,IM)
     *                                  + EFZ*EFPOL2(3,IM)
                        DINDD2FMO(1,IPX)= EFXD*EFPOL2(1,IM)
     *                                  + EFYD*EFPOL2(7,IM)
     *                                  + EFZD*EFPOL2(8,IM)
                        DINDD2FMO(2,IPX)= EFXD*EFPOL2(4,IM)
     *                                  + EFYD*EFPOL2(2,IM)
     *                                  + EFZD*EFPOL2(9,IM)
                        DINDD2FMO(3,IPX)= EFXD*EFPOL2(5,IM)
     *                                  + EFYD*EFPOL2(6,IM)
     *                                  + EFZD*EFPOL2(3,IM)
                     END IF
                     IF(I.EQ.1) THEN
                        DIPX  = DIPX + DIND2FMO(1,IPX)
                        DIPY  = DIPY + DIND2FMO(2,IPX)
                        DIPZ  = DIPZ + DIND2FMO(3,IPX)
                        DDIPX = DDIPX + DINDD2FMO(1,IPX)
                        DDIPY = DDIPY + DINDD2FMO(2,IPX)
                        DDIPZ = DDIPZ + DINDD2FMO(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
                     DINDFMO(1,IPX) = EFX*EFPOL(1,IPX)
     *                              + EFY*EFPOL(4,IPX)
     *                              + EFZ*EFPOL(5,IPX)
                     DINDFMO(2,IPX) = EFX*EFPOL(7,IPX)
     *                              + EFY*EFPOL(2,IPX)
     *                              + EFZ*EFPOL(6,IPX)
                     DINDFMO(3,IPX) = EFX*EFPOL(8,IPX)
     *                              + EFY*EFPOL(9,IPX)
     *                              + EFZ*EFPOL(3,IPX)
                     DINDDFMO(1,IPX)= EFXD*EFPOL(1,IPX)
     *                              + EFYD*EFPOL(7,IPX)
     *                              + EFZD*EFPOL(8,IPX)
                     DINDDFMO(2,IPX)= EFXD*EFPOL(4,IPX)
     *                              + EFYD*EFPOL(2,IPX)
     *                              + EFZD*EFPOL(9,IPX)
                     DINDDFMO(3,IPX)= EFXD*EFPOL(5,IPX)
     *                              + EFYD*EFPOL(6,IPX)
     *                              + EFZD*EFPOL(3,IPX)
                  END IF
                  IF(I.EQ.1) THEN
                     DIPX  = DIPX + DINDFMO(1,IPX)
                     DIPY  = DIPY + DINDFMO(2,IPX)
                     DIPZ  = DIPZ + DINDFMO(3,IPX)
                     DDIPX = DDIPX + DINDDFMO(1,IPX)
                     DDIPY = DDIPY + DINDDFMO(2,IPX)
                     DDIPZ = DDIPZ + DINDDFMO(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(DIND2FMO, 1,3*NPTTPT)
               CALL VCLR(DINDD2FMO,1,3*NPTTPT)
               DO IPX=NIN1,NIN1+NIN2-1
                  DIND2FMO(1,IPX) = DIND2C(1,IPX)
                  DIND2FMO(2,IPX) = DIND2C(2,IPX)
                  DIND2FMO(3,IPX) = DIND2C(3,IPX)
                  DINDD2FMO(1,IPX) = DINDD2C(1,IPX)
                  DINDD2FMO(2,IPX) = DINDD2C(2,IPX)
                  DINDD2FMO(3,IPX) = DINDD2C(3,IPX)
               ENDDO
            ELSE
               CALL VCLR(DINDFMO, 1,3*NPTTPT)
               CALL VCLR(DINDDFMO,1,3*NPTTPT)
               DO IPX=NIN1,NIN1+NIN2-1
                  DINDFMO(1,IPX) = DINDC(1,IPX)
                  DINDFMO(2,IPX) = DINDC(2,IPX)
                  DINDFMO(3,IPX) = DINDC(3,IPX)
                  DINDDFMO(1,IPX) = DINDDC(1,IPX)
                  DINDDFMO(2,IPX) = DINDDC(2,IPX)
                  DINDDFMO(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,DINDD2FMO,3*NPTTPT)
            CALL DDI_GSUMF(771,DIND2FMO, 3*NPTTPT)
         ELSE
            CALL DDI_GSUMF(770,DINDDFMO ,3*NPTTPT)
            CALL DDI_GSUMF(771,DINDFMO , 3*NPTTPT)
         END IF
      END IF
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 FLDINDFMO(EFADD,DINDFMO,DIND2FMO)
C
      RETURN
C
 9010 FORMAT(1X,'DIPITFMO FAILED TO CONVERGE IN',I5,
     *       'ITERATIONS, NACC=',I4)
      END
C*MODULE EFPFMO  *DECK EFPCMPFMO
      SUBROUTINE EFPCMPFMO(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, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12, MXRT=100)
C
      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)
C        FRGFLD IS CALCULATED IN EFPFMODER
      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
         TOTALP=SQRT(PMUX*PMUX+PMUY*PMUY+PMUZ*PMUZ)
         IF(TOTALP.GT.0.2D+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
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 EFPFMO  *DECK FLDINDFMO
      SUBROUTINE FLDINDFMO(EFADD,DINDFMO,DIND2FMO)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 POLNAM,DPOLNAM
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 /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
      DIMENSION EFADD(3,NPTTPT)
      DIMENSION DINDFMO(3,NPTTPT),DIND2FMO(3,NPTTPT)
C
      PARAMETER (ONE=1.00D+00)
C
      DATA THREE/3.0D+00/, ZERO/0.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 = DIND2FMO(1,IC2)*X
     *                   + DIND2FMO(2,IC2)*Y
     *                   + DIND2FMO(3,IC2)*Z
                     DOT = THREE * DOT / R5
                     FXIND=FXIND-SWF*DIND2FMO(1,IC2)/R3+SWF*DOT*X
                     FYIND=FYIND-SWF*DIND2FMO(2,IC2)/R3+SWF*DOT*Y
                     FZIND=FZIND-SWF*DIND2FMO(3,IC2)/R3+SWF*DOT*Z
                  ELSE
                     DOT1 = DINDFMO(1,IC2)*X
     *                    + DINDFMO(2,IC2)*Y
     *                    + DINDFMO(3,IC2)*Z
                     DOT = THREE * DOT1 / R5
                     FXIND=FXIND + SWF*((ONE+P1)*
     *                          (-DINDFMO(1,IC2)/R3 + DOT*X)+P2*DOT1*X)
                     FYIND=FYIND + SWF*((ONE+P1)*
     *                          (-DINDFMO(2,IC2)/R3 + DOT*Y)+P2*DOT1*Y)
                     FZIND=FZIND + SWF*((ONE+P1)*
     *                          (-DINDFMO(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
C*MODULE EFPFMO  *DECK DEFEFCL
      SUBROUTINE DEFEFCL(QUA)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      DIMENSION QUA(6,*)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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 /FMCOM / XX(1)
      COMMON /EFFMOP/ LEFLDF,LEFADDF,LDINDF,LDINDDF,
     *                LDIND2F,LDINDD2F,LDEFEF,LDEFTF,LTORQF
C
C     IF (MASWRK) THEN
C     WRITE (6,*) 'DINDFMO FOR TEST2'
C     DO I = 1, NPTTPT
C       WRITE (6,'(3F16.8)') (DINDF(J,I),J=1,3)
C     END DO
C     WRITE (6,*) 'DINDDFMO FOR TEST2'
C     DO I = 1, NPTTPT
C       WRITE (6,'(3F16.8)') (DINDDF(J,I),J=1,3)
C     END DO
C     ENDIF
C
      IF (IEFP.EQ.1) THEN
        IF (IEFC.EQ.1) THEN
C          CALL DCHINDFMO(XX(LDINDF),XX(LDINDDF),XX(LDEFEF))
          CALL DCHIND(XX(LDINDF),XX(LDINDDF),XX(LDEFEF),1)
        END IF
        IF (IEFD.EQ.1) THEN
C          CALL DDPINDFMO(XX(LDINDF),XX(LDINDDF),XX(LDEFEF))
          CALL DDPIND(XX(LDINDF),XX(LDINDDF),XX(LDEFEF),1)
        END IF
        IF (IEFQ.EQ.1) THEN
C          CALL DQDINDFMO(QUA,XX(LDINDF),XX(LDINDDF),XX(LDEFEF))
          CALL DQDIND(QUA,XX(LDINDF),XX(LDINDDF),XX(LDEFEF),1)
        END IF
C        CALL DININFMO(XX(LDINDF),XX(LDINDDF),XX(LDEFEF))
        CALL DININ(XX(LDINDF),XX(LDINDDF),XX(LDEFEF),1)
      END IF
C
      IF(GOPARR) CALL DDI_GSUMF(1520,XX(LDEFEF),6*NFRG)
C
C     AMEMA = 0.0D+00
C     IF (MASWRK) THEN
C       DO I = 1, 6*NFRG
C         AMEMA = AMEMA + (XX(LDEFEF+I-1))**2
C       END DO
C       AMEMA = SQRT(AMEMA/(6*NFRG))
C       WRITE (6,*) 'DEFEF RMS', AMEMA
C     ENDIF
C
C     IF (MASWRK) THEN
C     WRITE(6,*) 'TEST DEFEF DEFT'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (DEFEF(J,I),J=1,3)
C     END DO
C     WRITE(6,*) 'TEST DEFEF TORQ'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (DEFEF(J,I),J=4,6)
C     END DO
C     ENDIF
C
      RETURN
      END
C*MODULE EFPFMO  *DECK EFPDEG
      SUBROUTINE EFPDEG(IDA, IBODY)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
      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 /FMCOM / XX(1)
      COMMON /EFFMOP/ LEFLDF,LEFADDF,LDINDF,LDINDDF,
     *                LDIND2F,LDINDD2F,LDEFEF,LDEFTF,LTORQF
C
C     ONLY MASTER CAN COME
C
C
      IF (IDA.EQ.0) RETURN
      DA = IDA
C
C     IF (MASWRK) THEN
C     WRITE (6,*) 'NFIRST NEND', NFIRST,NEND
C     WRITE (6,*) 'DA ISWENUC', DA, ISWENUC
C     WRITE(6,*) 'TEST DEFT BEFORE'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (DEFT(J,I),J=1,3)
C     END DO
C     WRITE(6,*) 'TEST TORQ BEFORE'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (TORQ(J,I),J=1,3)
C     END DO
C     WRITE(6,*) 'TEST DESFRG DEFT'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (DESFRG(J,I),J=1,3)
C     END DO
C     WRITE(6,*) 'TEST DESFRG TORQ'
C     DO I = 1, NFRG
C       WRITE (6,'(3F16.8)') (DESFRG(J,I),J=4,6)
C     END DO
C     END IF
C
      II = (IBODY-1)*NFRG*3
      JJ = LDEFTF + II
      KK = LTORQF + II
      MM = -1
      NN = -1
      DO I = 1, NFRG
        LD = JJ + MM
        LT = KK + MM
        LF = LDEFEF + NN
        LL = LF + 3
        DO J = 1, 3
          XX(LD+J) = XX(LD+J) + DA*(DEFT(J,I) - XX(LF+J))
          XX(LT+J) = XX(LT+J) + DA*(TORQ(J,I) - XX(LL+J))
        END DO
        MM = MM + 3
        NN = NN + 6
      END DO
C     IF (MASWRK) THEN
C     WRITE(6,*) 'TEST DEFTTMP'
C     DO I = 1, NFRG*3
C       WRITE (6,'(3F16.8)') (DEFTTMP(J,I),J=1,3)
C     END DO
C     WRITE(6,*) 'TEST TORQTTMP'
C     DO I = 1, NFRG*3
C       WRITE (6,'(3F16.8)') (TORQTMP(J,I),J=1,3)
C     END DO
C     END IF
      RETURN
      END
C*MODULE EFPFMO  *DECK EFPFMOIO
      SUBROUTINE EFPFMOIO(IO)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (MXPT=2000, MXFRG=1050, MXFGPT=12000, MXDFG=5,
     *           MXAO=8192, MXDPPT=MXFRG*MXDFG*12)
C
C
      COMMON /FMCOM / XX(1)
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      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 /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      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 /EFFMIO/ L75,L81,L267,L268,L269
C
C     IF (MASWRK) WRITE(6,*) L75,L81,L267,L268,L269
      MXMO=0
      MXBF=0
      DO I=1,NFRG
        MXMO=MAX(MXMO,NORB(I))
        MXBF=MAX(MXBF,NPBF(I))
      END DO
      MXMO2=(MXMO*MXMO+MXMO)/2
      LENPV=MXBF*NTMO
      LENFM=MXMO2*NFRG
      LENCV=MXBF*NTAO
      MXRPTS=0
      DO IFRG=1,NFRG
        MXRPTS = MAX(NRPTS(IFRG),MXRPTS)
      END DO
      MXRPTS=MXRPTS+1
      NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
      LEN = MXRPTS*MXRPTS*NDFRG2
C     IF (MASWRK) WRITE(6,*) LEN,LENPV,LENFM,LENCV
      IF (IO.EQ.1) THEN
        IF (NTMO.GT.0) THEN
          CALL DAWRIT(IDAF,IODA,XX(L267),LENPV,267,0)
          CALL DAWRIT(IDAF,IODA,XX(L268),LENFM,268,0)
          CALL DAWRIT(IDAF,IODA,XX(L269),LENCV,269,0)
        ENDIF
C
        IF (NFRG.GT.1) THEN
C         IF (MASWRK) THEN
C           WRITE (6,*) 'WRITING EFP1 PARAMETERS'
C           WRITE (6,'(6F16.8)') (XX(L75+JK),JK=0,2*LEN-1)
C         ENDIF
          CALL DAWRIT(IDAF,IODA,XX(L75),2*LEN,75,0)
C
          IF (NTMO.GT.0) THEN
            MXDPTS=0
            DO IFRG=1,NFRG
              MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
            END DO
            MXDPTS=MXDPTS+1
            NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
            LEN = MXDPTS*MXDPTS*NDFRG2
            CALL DAWRIT(IDAF,IODA,XX(L81),6*LEN,81,0)
          ENDIF
        END IF
      ELSE
        IF (NTMO.GT.0) THEN
          CALL DAREAD(IDAF,IODA,XX(L267),LENPV,267,0)
          CALL DAREAD(IDAF,IODA,XX(L268),LENFM,268,0)
          CALL DAREAD(IDAF,IODA,XX(L269),LENCV,269,0)
        END IF
C
        IF (NFRG.GT.1) THEN
          CALL DAREAD(IDAF,IODA,XX(L75),2*LEN,75,0)
C
          IF (NTMO.GT.0) THEN
            MXDPTS=0
            DO IFRG=1,NFRG
              MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
            END DO
            MXDPTS=MXDPTS+1
            NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
            LEN = MXDPTS*MXDPTS*NDFRG2
            CALL DAREAD(IDAF,IODA,XX(L81),6*LEN,81,0)
          ENDIF
        END IF
      ENDIF
C
      RETURN
      END
C*MODULE EFPFMO  *DECK DEFPCLR
      SUBROUTINE DEFPCLR
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (MXFRG=1050)
C
      COMMON /EFPGES/ DESFRG(6,MXFRG)
      COMMON /EFPGCX/ DCTFRG(6,MXFRG),DXRFRG(6,MXFRG)
      CALL VCLR(DESFRG,1,6*MXFRG)
      CALL VCLR(DCTFRG,1,6*MXFRG)
      CALL VCLR(DXRFRG,1,6*MXFRG)
C
      RETURN
      END
C
C*MODULE EFPFMO  *DECK EFPMEM
      SUBROUTINE EFPMEM(NDER,LAST,NEEDEFP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXPT=2000,MXFRG=1050,MXFGPT=12000,MXAO=8192,MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
      COMMON /EFFMIO/ L75,L81,L267,L268,L269
      COMMON /EFFMOP/ LEFLDFMO,LEFADDFMO,LDINDFMO,LDINDDFMO,
     *                LDIND2FMO,LDINDD2FMO,LDEFEF,LDEFTF,LTORQF
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      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 /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB1(MXFRG),
     *                NPBF(MXFRG),NTMO
      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,
     *                IEFPOL1,IEFPOL2,NDTTPT
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      LEFLDFMO   = LAST
      LEFADDFMO  = LEFLDFMO   + NPTTPT*3
      LDINDFMO   = LEFADDFMO  + NPTTPT*3
      LDINDDFMO  = LDINDFMO   + NPTTPT*3
      LDIND2FMO  = LDINDDFMO  + NPTTPT*3
      LDINDD2FMO = LDIND2FMO  + NPTTPT*3
      LAST       = LDINDD2FMO + NPTTPT*3
      IF (NDER.GT.0) THEN
        LDEFEF = LAST
        LDEFTF = LDEFEF + 6*NFRG
        LTORQF = LDEFTF + 3*NFRG*NBODY
        LAST   = LTORQF + 3*NFRG*NBODY
      ELSE
        LDEFEF = LAST
        LDEFTF = LAST
        LTORQF = LAST
      ENDIF
      MXMO=0
      MXBF=0
      DO I=1,NFRG
        MXMO=MAX(MXMO,NORB1(I))
        MXBF=MAX(MXBF,NPBF(I))
      END DO
      MXMO2=(MXMO*MXMO+MXMO)/2
      LENPV=MXBF*NTMO
      LENFM=MXMO2*NFRG
      LENCV=MXBF*NTAO
      MXRPTS=0
      DO IFRG=1,NFRG
        MXRPTS = MAX(NRPTS(IFRG),MXRPTS)
      END DO
      MXRPTS=MXRPTS+1
      NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
      LENN1 = MXRPTS*MXRPTS*NDFRG2
      L75  = LAST
      L267 = L75  + LENN1*2
      L268 = L267 + LENPV
      L269 = L268 + LENFM
      LAST = L269 + LENCV
      IF (NTMO.EQ.0.AND.NFRG.GT.0) THEN
        MXDPTS=0
        DO IFRG=1,NFRG
          MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
        END DO
        MXDPTS=MXDPTS+1
        NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
        LENN2 = MXDPTS*MXDPTS*NDFRG2
        L81  = LAST
        LAST = L81 + LENN2*6
      ENDIF
      NEEDEFP=LAST-LEFLDFMO
      RETURN
      END
C*MODULE EFPFMO  *DECK EFEFDER
      SUBROUTINE EFEFDER(EF3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT
C
      PARAMETER (MXAO=8192, MXDFG=5, MXPT=2000, MXFRG=1050,
     *           MXFGPT=12000, MXDPPT=MXFRG*MXDFG*12)
C
      DOUBLE PRECISION MD
C
      DIMENSION EF3(3,*)
C
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /EFPGES/ DESFRG(6,MXFRG)
      COMMON /EFPOTD/ METHOD,INABIO,MOVE
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FMCOM / XX(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
C     COSMO CHANGES
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
      COMMON /INTFRG/ EFF(3,MXFRG),EFF2(3,MXFRG)
C
      INTEGER DDI_NP, DDI_ME
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HVNNDER  /, GRD1/8HGRD1    /
      DATA MD /8HMD      /
C
C     ----- GRADIENT OF NUCLEAR REPULSION ENERGY -----
C     PLUS NUCLEAR CONTRIBUTIONS FROM VARIOUS SOLVENT MODELS
C
C
C     DERIVATIVES OF EFP-EFP INTERACTIONS (MULTIPOLES, REPULSION)
C     FOR EFP/FMO ARE CALCULATED HERE. THIS IS TAKEN FROM VNNDER.
C     ONLY INDUCED DIPOLE TERMS ARE CALCULATED IN VNNDER.
C
      CALL DDI_NPROC( DDI_NP, DDI_ME )
C
C        NOTE THAT -OUT- DOES NOT INCLUDE -MASWRK- SO ALL NODES
C        CAN CALL THE PARALLEL DEBUGGING ROUTINE (W/GLOBAL SUM).
C
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
C
C     ----- EFFECTIVE FRAGMENT POTENTIAL -----
C     MULTIPOLE DERIVATIVE TERMS ARE NOW DONE ELSEWHERE, IN -ENUC-
C
      IF(NFRG.GT.1) THEN
C
         IF(OUT) CALL EFP_GRD_DBG0('NUCLEAR ',EF3,ATORQ)
C
         IF (METHOD .GE. 4) GO TO 400
C
C   ----- GRADIENT TERMS SPECIFIC TO THE -EFP1- MODEL:
C
         IF(NTMO.EQ.0) THEN
C
C                GRADIENT OF THE FITTED EXCHANGE REPULSION TERM
C
            IF(IREP.EQ.1) THEN
               MXRPTS=0
               DO IFRG=1,NFRG
                  MXRPTS = MAX(NRPTS(IFRG),MXRPTS)
               ENDDO
               MXRPTS=MXRPTS+1
               NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
               LEN = MXRPTS*MXRPTS*NDFRG2
               CALL VALFM(LOADFM)
               LAREL = LOADFM + 1
               LCREL = LAREL  + LEN
               LAST  = LCREL  + LEN
               NEED = LAST - LOADFM - 1
               CALL GETFM(NEED)
               CALL DAREAD(IDAF,IODA,XX(LAREL),2*LEN,75,0)
               CALL DREPUL(XX(LAREL),XX(LCREL),MXRPTS,NDFRG2)
               IF(OUT) CALL EFP_GRD_DBG('DREPUL  ',DESFRG)
               CALL RETFM(NEED)
            END IF
C
C    THERE IS A DISPERSON TERM IN THE MP2 FLAVOR OF -EFP1- MODEL
C
            IF(IEFP2.EQ.1) THEN
               MXDPTS=0
               DO IFRG=1,NFRG
                  MXDPTS = MAX(NDPTS(IFRG),MXDPTS)
               ENDDO
               MXDPTS=MXDPTS+1
               NDFRG2 = (NDFRG*NDFRG+NDFRG)/2
               LEN = MXDPTS*MXDPTS*NDFRG2
               CALL VALFM(LOADFM)
               LAREL = LOADFM + 1
               LDREL = LAREL  + LEN
               LCREL = LDREL  + LEN
               LXREL = LCREL  + LEN
               LBREL = LXREL  + LEN
               LEREL = LBREL  + LEN
               LAST  = LEREL  + LEN
               NEED  = LAST - LOADFM - 1
               CALL GETFM(NEED)
               CALL DAREAD(IDAF,IODA,XX(LAREL),6*LEN,81,0)
C
C     THERE'S BOTH A C6 AND C8 CONTRIBUTION TO DISPERSION
C     NOTE THAT THE -EFF2- VALUES DON'T GET PUT INTO THE GRADIENT
C     UNTIL THE CALL TO THE -EFFT- ROUTINE CALLED DURING INTEGRALS.
C
               CALL VCLR(EFF2,1,3*NFRG)
C
               CN=6.0D+00
               CALL DDISPR(XX(LAREL),XX(LDREL),XX(LCREL),
     *                     MXDPTS,NDFRG2,CN)
               IF(OUT) CALL EFP_GRD_DBG('DDISPR-6',DESFRG)
C
               CN=8.0D+00
               CALL DDISPR(XX(LXREL),XX(LBREL),XX(LEREL),
     *                     MXDPTS,NDFRG2,CN)
               IF(OUT) CALL EFP_GRD_DBG('DDISPR-8',DESFRG)
               CALL RETFM(NEED)
            END IF
         END IF
C
C        --- EFP2 EX-REP AND CHARGE TRANSFER GRADIENTS ---
C        FOR MD, THESE GRADIENTS HAVE BEEN DONE EARLIER, FOR REASONS
C        OF SPEED, SEE -GCTANA- AND -GEXREP- CALLS IN ROUTINE EFSP!
C
         IF(RUNTYP.EQ.MD) GO TO 390
         IF(NTMO.EQ.0) GO TO 390
C
C        - EXCHANGE REPULSION GRADIENTS
C
         MXMO=0
         MXBF=0
         DO I = 1, NFRG
            MXMO=MAX(MXMO,NORB(I))
            MXBF=MAX(MXBF,NPBF(I))
         END DO
         MXMO2=(MXMO*MXMO+MXMO)/2
         CALL VALFM(LOADFM)
         LPROVEC = LOADFM  + 1
         LDPVEC  = LPROVEC + MXBF*NTMO
         LFOCKMA = LDPVEC  + MXBF*NTMO*3
         LSMAT   = LFOCKMA + MXMO2*NFRG
         LTMAT   = LSMAT   + MXBF*MXBF
         LWRK    = LTMAT   + MXBF*MXBF
         LSAB    = LWRK    + MXBF
         LTAB    = LSAB    + MXMO*MXMO
         LFASQ   = LTAB    + MXMO*MXMO
         LFBSQ   = LFASQ   + MXMO*MXMO
         LDSINT  = LFBSQ   + MXMO*MXMO
         LDTINT  = LDSINT  + MXBF*MXBF*6
         LDSAB   = LDTINT  + MXBF*MXBF*6
         LDTAB   = LDSAB   + MXMO*MXMO*6
         LWRK2   = LDTAB   + MXMO*MXMO*6
         LGXRFRG = LWRK2   + MXMO*MXMO
         LXRFRG  = LGXRFRG + NFRG*NFRG*6
         LAST    = LXRFRG  + NFRG*NFRG*6*4
         NEED=LAST-LOADFM-1
         CALL GETFM(NEED)
C
         LENPV=MXBF*NTMO
         LENFM=MXMO2*NFRG
         CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
         CALL DAREAD(IDAF,IODA,XX(LFOCKMA),LENFM,268,0)
         CALL DERVEC(1,XX(LDPVEC),XX(LPROVEC),MXBF,NTMO)
         CALL GEXREP(ENXR,XX(LPROVEC),XX(LFOCKMA),XX(LSMAT),
     *               XX(LTMAT),XX(LWRK),XX(LSAB),XX(LTAB),
     *               XX(LFASQ),XX(LFBSQ),XX(LDPVEC),XX(LDSINT),
     *               XX(LDTINT),XX(LDSAB),XX(LDTAB),XX(LWRK2),
     *               XX(LGXRFRG),MXBF,MXMO,XX(LXRFRG))
         CALL RETFM(NEED)
C
C        - CHARGE TRANSFER GRADIENTS
C
         IF(NTAO.EQ.0) GO TO 390
C
         MXAOEFP=0
         MXBFEFP=0
         DO I = 1, NFRG
            MXAOEFP=MAX(MXAOEFP,NAO(ISET(I)))
            MXBFEFP=MAX(MXBFEFP,NPBF(I))
         END DO
         CALL VALFM(LOADFM)
         LCANVEC = LOADFM  + 1
         LDCVEC  = LCANVEC + MXBFEFP*NTAO
         LSMAT   = LDCVEC  + MXBFEFP*NTAO*3
         LTFRG   = LSMAT   + MXBFEFP*MXBFEFP
         LTMAT   = LTFRG   + MXAOEFP*MXAOEFP*MXDFG
         LWRK    = LTMAT   + MXBFEFP*MXBFEFP
         LSAB    = LWRK    + MXBFEFP
         LTAB    = LSAB    + MXAOEFP*MXAOEFP
         LVEFP   = LTAB    + MXAOEFP*MXAOEFP
         LVABEFB = LVEFP   + MXBFEFP*MXBFEFP
         LVABEFA = LVABEFB + MXAOEFP*MXAOEFP
         LVAAEFB = LVABEFA + MXAOEFP*MXAOEFP
         LVBBEFA = LVAAEFB + MXAOEFP*MXAOEFP
         LDSINT  = LVBBEFA + MXAOEFP*MXAOEFP
         LDTINT  = LDSINT  + MXBFEFP*MXBFEFP*6
         LDSAB   = LDTINT  + MXBFEFP*MXBFEFP*6
         LDTAB   = LDSAB   + MXAOEFP*MXAOEFP*6
         LWRK2   = LDTAB   + MXAOEFP*MXAOEFP*6
         LDABEFB = LWRK2   + MXAOEFP*MXAOEFP
         LDABEFA = LDABEFB + MXAOEFP*MXAOEFP*6
         LDAAEFB = LDABEFA + MXAOEFP*MXAOEFP*6
         LDBBEFA = LDAAEFB + MXAOEFP*MXAOEFP*6
         LDVEFP  = LDBBEFA + MXAOEFP*MXAOEFP*6
         LGCTFRG = LDVEFP  + MXBFEFP*MXBFEFP*6
         LIDOIT  = LGCTFRG + NFRG*NFRG*6
         LQUAD   = LIDOIT  + NMTTPT
         LAST    = LQUAD   + 6*NMTTPT
         NEED=LAST-LOADFM-1
         CALL GETFM(NEED)
C
         LENCV=MXBFEFP*NTAO
         CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
         CALL DERVEC(2,XX(LDCVEC),XX(LCANVEC),MXBFEFP,NTAO)
         CALL GCTANA(CHTNRG,XX(LCANVEC),XX(LSMAT),XX(LTFRG),
     *               XX(LTMAT),XX(LWRK),XX(LSAB),XX(LTAB),
     *               XX(LVEFP),XX(LVABEFB),XX(LVABEFA),
     *               XX(LVAAEFB),XX(LVBBEFA),XX(LDCVEC),
     *               XX(LDSINT),XX(LDTINT),XX(LDSAB),XX(LDTAB),
     *               XX(LWRK2),XX(LDABEFB),XX(LDABEFA),
     *               XX(LDAAEFB),XX(LDBBEFA),XX(LDVEFP),
     *               XX(LGCTFRG),XX(LIDOIT),XX(LQUAD),MXBFEFP,MXAOEFP)
         CALL RETFM(NEED)
C
  390    CONTINUE
      END IF
C
C         THIS ENDS ALL COMPUTATIONS OF EFP GRADIENTS,
C         EXCEPT CLEANUPS, SUMS, ETC.
C
C     ----- FRAGMENT TORQUES -----
C
  400 CONTINUE
C
      RETURN
C
      END
C*MODULE EFPFMO  *DECK EFEFCL
      SUBROUTINE EFEFCL
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      PARAMETER (ZERO=0.0D+00, PT5=0.5D+00)
C
      CHARACTER*8 POLNAM,DPOLNAM
C
      LOGICAL GOPARR,MASWRK,DSKWRK
      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 /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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /EFFMOP/ LEFLDFMO,LEFADDFMO,LDINDFMO,LDINDDFMO,
     *                LDIND2FMO,LDINDD2FMO,LDEFEF,LDEFTF,LTORQF
C
      CALL VCLR(XX(LEFLDFMO),  1,NPTTPT*3)
      CALL VCLR(XX(LEFADDFMO), 1,NPTTPT*3)
      CALL VCLR(XX(LDINDFMO),  1,NPTTPT*3)
      CALL VCLR(XX(LDIND2FMO), 1,NPTTPT*3)
      CALL VCLR(XX(LDINDDFMO), 1,NPTTPT*3)
      CALL VCLR(XX(LDINDD2FMO),1,NPTTPT*3)
C
      II = 0
      DO LEFP = 1, NPTTPT
        XX(LEFLDFMO+II  ) = PPEFX(LEFP)
        XX(LEFLDFMO+II+1) = PPEFY(LEFP)
        XX(LEFLDFMO+II+2) = PPEFZ(LEFP)
        II = II + 3
      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
      LAST     = LDIND2_C + 6*NPTTPT
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL DIPITFMO(XX(LEFLDFMO),XX(LEFADDFMO),XX(LDINDFMO),
     *              XX(LDIND2FMO),XX(LDINDDFMO),XX(LDINDD2FMO),
     *              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),NFRG)
      CALL RETFM(NEED)
C
C     EFP-EFP INTERACTION ENERGY
C
      IPCOUNT = ME - 1
      EFPEFP = ZERO
      II = 0
      DO 24 LEFP=1,NPTTPT
        IF (GOPARR) THEN
          IPCOUNT = IPCOUNT + 1
          IF (MOD(IPCOUNT,NPROC).NE.0) THEN
            II = II + 3
            GO TO 24
          END IF
        END IF
C
C       FIELD DUE TO EFP
        PEFPX=PPEFX(LEFP)
        PEFPY=PPEFY(LEFP)
        PEFPZ=PPEFZ(LEFP)
C
C       POLARIZABILITY
        AXX=EFPOL(1,LEFP)
        AYY=EFPOL(2,LEFP)
        AZZ=EFPOL(3,LEFP)
        AXY=EFPOL(4,LEFP)
        AXZ=EFPOL(5,LEFP)
        AYZ=EFPOL(6,LEFP)
        AYX=EFPOL(7,LEFP)
        AZX=EFPOL(8,LEFP)
        AZY=EFPOL(9,LEFP)
C
C       EFP PART
        ADDFX=XX(LEFADDFMO+II  )
        ADDFY=XX(LEFADDFMO+II+1)
        ADDFZ=XX(LEFADDFMO+II+2)
C
C       FOR FIELD DUE TO INDUCED DIPOLE BY EFP
        EFXF = (AXY-AYX)*ADDFY+(AXZ-AZX)*ADDFZ
        EFYF = (AYX-AXY)*ADDFX+(AYZ-AZY)*ADDFZ
        EFZF = (AZX-AXZ)*ADDFX+(AZY-AYZ)*ADDFY
C
C       FOR FIELD DUE TO EFP MOLTIPOLES
        EFPMUX = AXX*PEFPX+AXY*PEFPY+AXZ*PEFPZ
        EFPMUY = AYX*PEFPX+AYY*PEFPY+AYZ*PEFPZ
        EFPMUZ = AZX*PEFPX+AZY*PEFPY+AZZ*PEFPZ
C
        EFPEFP = EFPEFP
     *           - ((EFPMUX+EFXF)*(PEFPX+ADDFX)
     *             +(EFPMUY+EFYF)*(PEFPY+ADDFY)
     *             +(EFPMUZ+EFZF)*(PEFPZ+ADDFZ))*PT5
        II = II + 3
 24   CONTINUE
C
      IF (GOPARR) CALL DDI_GSUMF(684,EFPEFP,1)
      REPNUCEFP = REPNUCEFP + EFPEFP
      IF (MASWRK) THEN
        WRITE(6,'(A,F16.8)') ' PURE EFP-EFP CONTRIBUTION', REPNUCEFP
      END IF
      RETURN
      END
C*MODULE EFPFMO  *DECK EFPFMODER
      SUBROUTINE EFPFMODER(NDER,DEFTF,TORQF,DEFEF)
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,MASWRK,DSKWRK
C
      PARAMETER (MXATM=2000, MXSH=5000, MXGTOT=20000,
     *           MXPT=2000, MXFRG=1050, MXFGPT=12000, MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION DEFTF(3,*), TORQF(3,*), DEFEF(6,*)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      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 /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      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 /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /EFPGES/ DESFRG(6,MXFRG)
      COMMON /EFPGCX/ DCTFRG(6,MXFRG),DXRFRG(6,MXFRG)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /EFFMOP/ LEFLDFMO,LEFADDFMO,LDINDFMO,LDINDDFMO,
     *                LDIND2FMO,LDINDD2FMO,LDEFEF,LDEFTF,LTORQF
C
C     THIS SUBROUTINE IS IN WORLD SCOPE
C
C     (1). THIS SUBROUTINE CALCULATES FIELD AT POLARIZABLE POINTS
C          DUE TO MULTIPOLES IN FRGFLD --> GET PPEFX TO PPEFZ
C
C     (2). DERIVATIVES OF PURE EFP-EFP INDUCED DIPOLES IN DEFEFCL
C          AND THE OTHER DERIVATIVES OF THE PURE EFP-EFP INTERACTIONS
C          IN EFEFDER
C
C
C     JUST GET DESFRG CALCULATED IN ENUC
      ISWENUC = 1
      TMPNUC = 0
      TMPNUC = TMPNUC+ENUC(NAT,ZAN,C)
      ISWENUC = 0
C
      IF (NFRG.LE.1) RETURN
      CALL FRGFLD
      CALL EFEFCL
      IF (NDER.GT.0) THEN
        NFRPTS = NMTTPT + NPTTPT + NRTTPT + NTPATM + NTMO + NDTTPT
        CALL VALFM(LOADFM)
        LWORK = LOADFM + 1
        LEF3  = LWORK  + 10*NMTTPT
        LAST  = LEF3   +  3*NFRPTS
        NEED = LAST - LOADFM - 1
        CALL GETFM(NEED)
C
C       CALCULATE DEFEF
        CALL DEFEFCL(XX(LWORK))
C       CALCULATE THE DERIVATIVES OF MULT-MULT AND EFP-EFP REPULSION
        CALL EFEFDER(XX(LEF3))
C
        CALL RETFM(NEED)
C
        IF(GOPARR) CALL DDI_GSUMF(778,DESFRG,6*MXFRG)
        IF (MASWRK) THEN
          DO IFRG = 1, NFRG
            DO ICRD = 1, 3
              DEFTF(ICRD,IFRG) = DEFTF(ICRD,IFRG) + DCTFRG(ICRD,IFRG)
     *                                            + DXRFRG(ICRD,IFRG)
     *                                            + DESFRG(ICRD,IFRG)
     *                                            + DEFEF(ICRD,IFRG)
              TORQF(ICRD,IFRG) = TORQF(ICRD,IFRG) + DCTFRG(ICRD+3,IFRG)
     *                                            + DXRFRG(ICRD+3,IFRG)
     *                                            + DESFRG(ICRD+3,IFRG)
     *                                            + DEFEF(ICRD+3,IFRG)
            END DO
          END DO
C
          NFRG_3 = 3*NFRG
          DO I = 1, NBODY - 1
            CALL DCOPY(NFRG_3,DEFTF,1,XX(LDEFTF+NFRG_3*I),1)
            CALL DCOPY(NFRG_3,TORQF,1,XX(LTORQF+NFRG_3*I),1)
          END DO
        END IF
      END IF
C
      RETURN
      END
C*MODULE EFPFMO  *DECK EFPPARL
      SUBROUTINE EFPPARL
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MXCPUEFP=1024)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      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 /EFPAR / LOCNFG,NLOW,NUP
      COMMON /EFPAR2/ NFRGAL(MXCPUEFP),NLOWAL(MXCPUEFP),NUPAL(MXCPUEFP),
     *                NPPTAL(MXCPUEFP)
      INTEGER DDI_NP,DDI_ME,FLO,CEIL,DIFF
C
C     PARTITION EFP'S OVER PROCESSORS
C
C
C     THIS IS TAKEN FROM SUBROUTINE START. THIS SETS UP THE PARAMETERS
C     FOR EFP PARALLELIZATION SCHEME.
C     WHEN ONE USES GDDI COMMUNICATER, THE PARALLELIZATION PARAMETER
C     FOR EFP MUST BE RESET.
C
      IF(NFRG.GT.0) THEN
        DDI_ME = ME
        DDI_NP = NPROC
        IF (MOD(NFRG,DDI_NP) .EQ. 0) THEN
          LOCNFG = NFRG / DDI_NP
          NLOW = DDI_ME * LOCNFG + 1
        ELSE
          CEIL = NFRG/DDI_NP + 1
          FLO = NFRG/DDI_NP
          DIFF = NFRG - DDI_NP*FLO
C
          IF (DDI_ME .GT. (DDI_NP-1-DIFF)) THEN
            LOCNFG = CEIL
            NLOW = DDI_ME * LOCNFG - (DDI_NP-1-DIFF)
          ELSE
            LOCNFG = FLO
            NLOW = DDI_ME * LOCNFG + 1
          END IF
        END IF
        NUP =  NLOW + LOCNFG - 1
        IF(GOPARR) THEN
C    THERE IS AN UNFORTUNATE COMMON BLOCK DIMENSION USED HERE.
          IF(DDI_NP.GT.MXCPUEFP) THEN
             IF(MASWRK) WRITE(IW,*)
     *          'ERROR IN PARALLEL EFP, MXCPUEFP,NPROC=',MXCPUEFP,DDI_NP
             CALL ABRT
          END IF
          CALL VICLR(NFRGAL,1,MXCPUEFP)
          CALL VICLR(NLOWAL,1,MXCPUEFP)
          CALL VICLR( NUPAL,1,MXCPUEFP)
          NFRGAL(DDI_ME+1) = LOCNFG
          NLOWAL(DDI_ME+1) = NLOW
           NUPAL(DDI_ME+1) = NUP
          CALL DDI_GSUMI(660,NFRGAL,DDI_NP)
          CALL DDI_GSUMI(661,NLOWAL,DDI_NP)
          CALL DDI_GSUMI(662, NUPAL,DDI_NP)
        ELSE
          NFRGAL(1) = NFRG
          NLOWAL(1) = 1
           NUPAL(1) = NFRG
        END IF
      END IF
      RETURN
      END
