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,LEFNUC
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /EFPFMO/ IEFFMO,LEFNUC,REPNUCEFP
C
C     IF BOTH $FMO AND $EFRAG EXIST,
C     EFP/FMO CALCULATION WILL BE CARRIED OUT. 
C
      IEOF = 0
      IEFFMO = 0
      LEFNUC = .TRUE.
      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
          IEFFMO = 1
        END IF
      END IF

      IF (MASWRK) WRITE(IW,*) '***** FMO/EFP CALCULATION STARTS *****'
      RETURN
C
      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 (MXPRO=64)
C
      DOUBLE PRECISION MD
      DOUBLE PRECISION MDTYP
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)
      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
      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 /MDSIM/  DT, BOXL, BOXLI, BATHT, EKINT, EKINR,
     *                  MDTYP
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      DATA MD         /8HMD      /
      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(MXPRO),NLOWAL(MXPRO),NUPAL(MXPRO),
     *                NPPTAL(MXPRO)
      INTEGER DDI_NP, DDI_ME
C
      PARAMETER (ZERO=0.00D+00, THREE=3.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,MXPRO)
        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
C
            IF(RUNTYP.EQ.MD) THEN
              COMX1 = EFCENT(1,IM)
              COMY1 = EFCENT(2,IM)
              COMZ1 = EFCENT(3,IM)
            END IF
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 IP = I1,I2
               IF(BUFF_PAS .EQ. DDI_ME) THEN
                  EFX = EFLD(1,IP)
                  EFY = EFLD(2,IP)
                  EFZ = EFLD(3,IP)
                  EFXD = EFLD(1,IP)
                  EFYD = EFLD(2,IP)
                  EFZD = EFLD(3,IP)
               ELSE
                  EFX =ZERO
                  EFY =ZERO
                  EFZ =ZERO
                  EFXD =ZERO
                  EFYD =ZERO
                  EFZD =ZERO
               END IF
               X1 = EFP(1,IP)
               Y1 = EFP(2,IP)
               Z1 = EFP(3,IP)
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
                  IF(RUNTYP.EQ.MD) THEN
                     COMX2 = EFCENT(1,JM)
                     COMY2 = EFCENT(2,JM)
                     COMZ2 = EFCENT(3,JM)
C
                     COMX = COMX1 - COMX2
                     COMY = COMY1 - COMY2
                     COMZ = COMZ1 - COMZ2
C
                     PBCX = BOXL * ANINT(COMX * BOXLI)
                     PBCY = BOXL * ANINT(COMY * BOXLI)
                     PBCZ = BOXL * ANINT(COMZ * BOXLI)
                  END IF
C
               SWF   = 1.0D+00
C
               XADD=ZERO
               YADD=ZERO
               ZADD=ZERO
               IF(XBOX.NE.ZERO.AND.YBOX.NE.ZERO.AND.ZBOX.NE.ZERO) THEN
               X12  = EFCENT(1,IM) - EFCENT(1,JM)
               Y12  = EFCENT(2,IM) - EFCENT(2,JM)
               Z12  = EFCENT(3,IM) - EFCENT(3,JM)
               XADD = XBOX * ANINT(X12/XBOX)
               YADD = YBOX * ANINT(Y12/YBOX)
               ZADD = ZBOX * ANINT(Z12/ZBOX)
               X12  = ABS(EFCENT(1,IM) - EFCENT(1,JM) - XADD)
               Y12  = ABS(EFCENT(2,IM) - EFCENT(2,JM) - YADD)
               Z12  = ABS(EFCENT(3,IM) - EFCENT(3,JM) - ZADD)
               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
               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
               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)
                     X2 = X2 + XADD
                     Y2 = Y2 + YADD
                     Z2 = Z2 + ZADD
                     X = X1 - X2
                     Y = Y1 - Y2
                     Z = Z1 - Z2
C
                     IF (RUNTYP.EQ.MD) THEN
                       X = X - PBCX
                       Y = Y - PBCY
                       Z = Z - PBCZ
                     END IF
C
                     R2 = X * X + Y * Y + Z * Z
                     R = SQRT(R2)
                     R3 = R * R * R
                     R5 = R * R * R3
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
                        IF(GOPARR) THEN
                           DOT = DIND_C(1,K1)*X
     *                         + DIND_C(2,K1)*Y
     *                         + DIND_C(3,K1)*Z
                           DOT = THREE * DOT / R5
                           EFX = EFX- SWF*DIND_C(1,K1)/R3 + SWF*DOT*X
                           EFY = EFY- SWF*DIND_C(2,K1)/R3 + SWF*DOT*Y
                           EFZ = EFZ- SWF*DIND_C(3,K1)/R3 + SWF*DOT*Z
C
                           DOTD = DIND_C(4,K1)*X
     *                          + DIND_C(5,K1)*Y
     *                          + DIND_C(6,K1)*Z
                           DOTD = THREE * DOTD / R5
                           EFXD =EFXD-SWF*DIND_C(4,K1)/R3+SWF*DOTD*X
                           EFYD =EFYD-SWF*DIND_C(5,K1)/R3+SWF*DOTD*Y
                           EFZD =EFZD-SWF*DIND_C(6,K1)/R3+SWF*DOTD*Z
                           K1 = K1 + 1
                        ELSE
                           DOT = DINDFMO(1,JP)*X
     *                         + DINDFMO(2,JP)*Y
     *                         + DINDFMO(3,JP)*Z
                           DOT = THREE * DOT / R5
                           EFX = EFX - SWF*DINDFMO(1,JP)/R3 + SWF*DOT*X
                           EFY = EFY - SWF*DINDFMO(2,JP)/R3 + SWF*DOT*Y
                           EFZ = EFZ - SWF*DINDFMO(3,JP)/R3 + SWF*DOT*Z
C
                           DOTD = DINDDFMO(1,JP)*X
     *                          + DINDDFMO(2,JP)*Y
     *                          + DINDDFMO(3,JP)*Z
                           DOTD = THREE * DOTD / R5
                           EFXD =EFXD -SWF*DINDDFMO(1,JP)/R3+SWF*DOTD*X
                           EFYD =EFYD -SWF*DINDDFMO(2,JP)/R3+SWF*DOTD*Y
                           EFZD =EFZD -SWF*DINDDFMO(3,JP)/R3+SWF*DOTD*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(IP.EQ.I2) THEN
                     IF(GOPARR) THEN
                        DIND2C(1,IP) =DIND2C(1,IP) +EFX*EFPOL2(1,IM)
     *                                             +EFY*EFPOL2(4,IM)
     *                                             +EFZ*EFPOL2(5,IM)
                        DIND2C(2,IP) =DIND2C(2,IP) +EFX*EFPOL2(7,IM)
     *                                             +EFY*EFPOL2(2,IM)
     *                                             +EFZ*EFPOL2(6,IM)
                        DIND2C(3,IP) =DIND2C(3,IP) +EFX*EFPOL2(8,IM)
     *                                             +EFY*EFPOL2(9,IM)
     *                                             +EFZ*EFPOL2(3,IM)
                        DINDD2C(1,IP)=DINDD2C(1,IP)+EFXD*EFPOL2(1,IM)
     *                                             +EFYD*EFPOL2(7,IM)
     *                                             +EFZD*EFPOL2(8,IM)
                        DINDD2C(2,IP)=DINDD2C(2,IP)+EFXD*EFPOL2(4,IM)
     *                                             +EFYD*EFPOL2(2,IM)
     *                                             +EFZD*EFPOL2(9,IM)
                        DINDD2C(3,IP)=DINDD2C(3,IP)+EFXD*EFPOL2(5,IM)
     *                                             +EFYD*EFPOL2(6,IM)
     *                                             +EFZD*EFPOL2(3,IM)
                     ELSE
                        DIND2FMO(1,IP) = EFX*EFPOL2(1,IM)
     *                                  +EFY*EFPOL2(4,IM)
     *                                  +EFZ*EFPOL2(5,IM)
                        DIND2FMO(2,IP) = EFX*EFPOL2(7,IM)
     *                                  +EFY*EFPOL2(2,IM)
     *                                  +EFZ*EFPOL2(6,IM)
                        DIND2FMO(3,IP) = EFX*EFPOL2(8,IM)
     *                                  +EFY*EFPOL2(9,IM)
     *                                  +EFZ*EFPOL2(3,IM)
                        DINDD2FMO(1,IP)= EFXD*EFPOL2(1,IM)
     *                                  +EFYD*EFPOL2(7,IM)
     *                                  +EFZD*EFPOL2(8,IM)
                        DINDD2FMO(2,IP)= EFXD*EFPOL2(4,IM)
     *                                  +EFYD*EFPOL2(2,IM)
     *                                  +EFZD*EFPOL2(9,IM)
                        DINDD2FMO(3,IP)= EFXD*EFPOL2(5,IM)
     *                                  +EFYD*EFPOL2(6,IM)
     *                                  +EFZD*EFPOL2(3,IM)
                     END IF
                     IF(I.EQ.1) THEN
                        DIPX  = DIPX + DIND2FMO(1,IP)
                        DIPY  = DIPY + DIND2FMO(2,IP)
                        DIPZ  = DIPZ + DIND2FMO(3,IP)
                        DDIPX = DDIPX + DINDD2FMO(1,IP)
                        DDIPY = DDIPY + DINDD2FMO(2,IP)
                        DDIPZ = DDIPZ + DINDD2FMO(3,IP)
                     END IF
                  END IF
               ELSE
                  IF(GOPARR) THEN
                     DINDC(1,IP) = DINDC(1,IP)+ EFX*EFPOL(1,IP)
     *                                        + EFY*EFPOL(4,IP)
     *                                        + EFZ*EFPOL(5,IP)
                     DINDC(2,IP) = DINDC(2,IP)+ EFX*EFPOL(7,IP)
     *                                        + EFY*EFPOL(2,IP)
     *                                        + EFZ*EFPOL(6,IP)
                     DINDC(3,IP) = DINDC(3,IP)+ EFX*EFPOL(8,IP)
     *                                        + EFY*EFPOL(9,IP)
     *                                        + EFZ*EFPOL(3,IP)
                     DINDDC(1,IP)=DINDDC(1,IP)+ EFXD*EFPOL(1,IP)
     *                                        + EFYD*EFPOL(7,IP)
     *                                        + EFZD*EFPOL(8,IP)
                     DINDDC(2,IP)=DINDDC(2,IP)+ EFXD*EFPOL(4,IP)
     *                                        + EFYD*EFPOL(2,IP)
     *                                        + EFZD*EFPOL(9,IP)
                     DINDDC(3,IP)=DINDDC(3,IP)+ EFXD*EFPOL(5,IP)
     *                                        + EFYD*EFPOL(6,IP)
     *                                        + EFZD*EFPOL(3,IP)
                  ELSE
                     DINDFMO(1,IP) =EFX*EFPOL(1,IP)+EFY*EFPOL(4,IP)
     *                             +EFZ*EFPOL(5,IP)
                     DINDFMO(2,IP) =EFX*EFPOL(7,IP)+EFY*EFPOL(2,IP)
     *                             +EFZ*EFPOL(6,IP)
                     DINDFMO(3,IP) =EFX*EFPOL(8,IP)+EFY*EFPOL(9,IP)
     *                             +EFZ*EFPOL(3,IP)
                     DINDDFMO(1,IP)=EFXD*EFPOL(1,IP)+EFYD*EFPOL(7,IP)
     *                             +EFZD*EFPOL(8,IP)
                     DINDDFMO(2,IP)=EFXD*EFPOL(4,IP)+EFYD*EFPOL(2,IP)
     *                             +EFZD*EFPOL(9,IP)
                     DINDDFMO(3,IP)=EFXD*EFPOL(5,IP)+EFYD*EFPOL(6,IP)
     *                             +EFZD*EFPOL(3,IP)
                  END IF
                  IF(I.EQ.1) THEN
                     DIPX  = DIPX + DINDFMO(1,IP)
                     DIPY  = DIPY + DINDFMO(2,IP)
                     DIPZ  = DIPZ + DINDFMO(3,IP)
                     DDIPX = DDIPX + DINDDFMO(1,IP)
                     DDIPY = DDIPY + DINDDFMO(2,IP)
                     DDIPZ = DDIPZ + DINDDFMO(3,IP)
                  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 IP=NIN1,NIN1+NIN2-1
                  DIND2FMO(1,IP) = DIND2C(1,IP)
                  DIND2FMO(2,IP) = DIND2C(2,IP)
                  DIND2FMO(3,IP) = DIND2C(3,IP)
                  DINDD2FMO(1,IP) = DINDD2C(1,IP)
                  DINDD2FMO(2,IP) = DINDD2C(2,IP)
                  DINDD2FMO(3,IP) = DINDD2C(3,IP)
               ENDDO
            ELSE
               CALL VCLR(DINDFMO, 1,3*NPTTPT)
               CALL VCLR(DINDDFMO,1,3*NPTTPT)
               DO IP=NIN1,NIN1+NIN2-1
                  DINDFMO(1,IP) = DINDC(1,IP)
                  DINDFMO(2,IP) = DINDC(2,IP)
                  DINDFMO(3,IP) = DINDC(3,IP)
                  DINDDFMO(1,IP) = DINDDC(1,IP)
                  DINDDFMO(2,IP) = DINDDC(2,IP)
                  DINDDFMO(3,IP) = DINDDC(3,IP)
               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,'DIPIT 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,EFLDFMO,EFADDFMO,
     *                  ABFLD,L2,L1,MCITER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000,
     *           MXRT=100, MXIFRQ=12, MXPT=2000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
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),
     *                                         ABFLD(3,NPTTPT)
CZ   *          EFLD(3,NPTTPT),EFADD(3,NPTTPT),ABFLD(3,NPTTPT)
      DIMENSION GSUMTMP(4)
CZ
      DIMENSION EFLDFMO(3,NPTTPT),EFADDFMO(3,NPTTPT)
CZ
C
      CHARACTER*8 POLNAM,DPOLNAM,DNAME,DPSNAM
C
      LOGICAL GOPARR,MASWRK,DSKWRK,SVDSKW,
     *        DIRSCF,FDIFF,DIRDIF,MINMEM
CZ
      LOGICAL LEFNUC,FIRST
CZ
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)
      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
      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,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
CZ
      COMMON /EFPFMO/ IEFFMO,LEFNUC,REPNUCEFP
      COMMON /EFFMO2/ EFLD(3,MXFGPT),EFADD(3,MXFGPT),
     *                DINDFMO(3,MXFGPT),DINDDFMO(3,MXFGPT),
     *                DIND2FMO(3,MXFGPT),DINDD2FMO(3,MXFGPT)
      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)
CZ
C
      DATA RHF,GVB/8HRHF     ,8HGVB     /
      DATA RMC/8HMCSCF   /
C
      SAVE EN0
CZ
      DATA FIRST/.TRUE./
      SAVE FIRST
CZ
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)
CZ       CALL FRGFLD
         IF (FIRST) 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
C        EFLD(1,LEFP)=ABFLD(1,LEFP) + PPNX(LEFP) + PPEFX(LEFP)
C        EFLD(2,LEFP)=ABFLD(2,LEFP) + PPNY(LEFP) + PPEFY(LEFP)
C        EFLD(3,LEFP)=ABFLD(3,LEFP) + PPNZ(LEFP) + PPEFZ(LEFP)
CZ       EFP PART
         IF (FIRST) THEN
           EFLD(1,LEFP) = PPEFX(LEFP)
           EFLD(2,LEFP) = PPEFY(LEFP)
           EFLD(3,LEFP) = PPEFZ(LEFP)
         END IF
CZ       FMO PART
         EFLDFMO(1,LEFP) = ABFLD(1,LEFP) + PPNX(LEFP) 
         EFLDFMO(2,LEFP) = ABFLD(2,LEFP) + PPNY(LEFP) 
         EFLDFMO(3,LEFP) = ABFLD(3,LEFP) + PPNZ(LEFP)
      ENDDO
C
C
C
C     *******************
C     ALWAYS CALL DIPIT !
C     *******************
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)
CZ    CALL DIPIT(EFLD,EFADD,XX(LDINDC),XX(LDINDDC),XX(LDIND2C),
CZ   *           XX(LDINDD2C),XX(LDIND_W),XX(LDIND_C),XX(LDIND2_W),
CZ   *           XX(LDIND2_C),NPTTPT,XX(LDIPNEW),XX(LDIPNWD),NFRG)
      IF (FIRST) THEN
CZ      EFP PART IS CALCULATED ONLY ONCE
        CALL DIPITFMO(EFLD,EFADD,DIND,DIND2,
     *                DINDD,DINDD2,
     *                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)
      END IF
CZ    FMO FRAGMENT PART
      CALL DIPITFMO(EFLDFMO,EFADDFMO,DINDFMO,DIND2FMO,
     *              DINDDFMO,DINDD2FMO,
     *              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
      IF(GOPARR) THEN
         ENUCR = ENUCR/NPROC
CZ
         IF (FIRST) THEN
           REPNUCEFP = REPNUCEFP/NPROC
         END IF
CZ
         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)
CZ       EFP PART 
         PMUXE  = DIND(1,LEFP)
         PMUYE  = DIND(2,LEFP)
         PMUZE  = DIND(3,LEFP)
         PMUXDE = DINDD(1,LEFP)
         PMUYDE = DINDD(2,LEFP)
         PMUZDE = DINDD(3,LEFP)
CZ       FMO PART
         PMUXF  = DINDFMO(1,LEFP)
         PMUYF  = DINDFMO(2,LEFP)
         PMUZF  = DINDFMO(3,LEFP)
         PMUXDF = DINDDFMO(1,LEFP)
         PMUYDF = DINDDFMO(2,LEFP)
         PMUZDF = DINDDFMO(3,LEFP)

         PMUX  = PMUXE + PMUXF
         PMUY  = PMUYE + PMUYF
         PMUZ  = PMUZE + PMUZF
         PMUXD = PMUXDE + PMUXDF
         PMUYD = PMUYDE + PMUYDF
         PMUZD = PMUZDE + PMUZDF
CZ
         PEX=ABFLD(1,LEFP)
         PEY=ABFLD(2,LEFP)
         PEZ=ABFLD(3,LEFP)
CZ       FIELD DUE TO EFP
         PEFPX=PPEFX(LEFP)
         PEFPY=PPEFY(LEFP)
         PEFPZ=PPEFZ(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
CZ
         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)
CZ       EFP PART
         ADDEX=EFADD(1,LEFP)
         ADDEY=EFADD(2,LEFP)
         ADDEZ=EFADD(3,LEFP)
CZ       FMO PART
         ADDFX=EFADDFMO(1,LEFP)
         ADDFY=EFADDFMO(2,LEFP)
         ADDFZ=EFADDFMO(3,LEFP)
CZ
         ADEFX=ADDEX+ADDFX
         ADEFY=ADDEY+ADDFY
         ADEFZ=ADDEZ+ADDFZ
CZ
C
C        --- CALCULATE ALPHA*FMU VALUES THAT ARE USED IN LOOPS BELOW ---
C
         EFX = AXX*ADEFX+AYX*ADEFY+AZX*ADEFZ
         EFY = AXY*ADEFX+AYY*ADEFY+AZY*ADEFZ
         EFZ = AXZ*ADEFX+AYZ*ADEFY+AZZ*ADEFZ
CZ       EFP PART
         EFXE = (AXY-AYX)*ADDEY+(AXZ-AZX)*ADDEZ
         EFYE = (AYX-AXY)*ADDEX+(AYZ-AZY)*ADDEZ
         EFZE = (AZX-AXZ)*ADDEX+(AZY-AYZ)*ADDEY
CZ       FMO PART
         EFXF = (AXY-AYX)*ADDFY+(AXZ-AZX)*ADDFZ
         EFYF = (AYX-AXY)*ADDFX+(AYZ-AZY)*ADDFZ
         EFZF = (AZX-AXZ)*ADDFX+(AZY-AYZ)*ADDFY

CZ       INDUCED DIPOLE JUST DUE TO AB INITIO
CZ       BY ELECTRONS
         PMUEX = AXX*PEX+AXY*PEY+AXZ*PEZ
         PMUEY = AYX*PEX+AYY*PEY+AYZ*PEZ
         PMUEZ = AZX*PEX+AZY*PEY+AZZ*PEZ
CZ       BY NUCLEI
         PMUNX = AXX*PNX+AXY*PNY+AXZ*PNZ
         PMUNY = AYX*PNX+AYY*PNY+AYZ*PNZ
         PMUNZ = AZX*PNX+AZY*PNY+AZZ*PNZ
CZ
CZ       BY EFP
         EFPMUX = AXX*PEFPX+AXY*PEFPY+AXZ*PEFPZ
         EFPMUY = AYX*PEFPX+AYY*PEFPY+AYZ*PEFPZ
         EFPMUZ = AZX*PEFPX+AZY*PEFPY+AZZ*PEFPZ
         IF (FIRST) THEN
           REPNUCEFP = REPNUCEFP
     *               - ((EFPMUX+EFXE)*(PEFPX+ADDEX)
     *                 +(EFPMUY+EFYE)*(PEFPY+ADDEY)
     *                 +(EFPMUZ+EFZE)*(PEFPZ+ADDEZ))*PT5
         END IF
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
CZ       ENUCR = ENUCR + ((-PMUX*PNX - PMUY*PNY - PMUZ*PNZ)
CZ   *                    + PMUXD*PEX + PMUYD*PEY + PMUZD*PEZ
CZ   *                    - PMUX*PPEFX(LEFP)
CZ   *                    - PMUY*PPEFY(LEFP)
CZ   *                    - PMUZ*PPEFZ(LEFP)
CZ   *                    - (EFX*PEX) - (EFY*PEY) - (EFZ*PEZ))*PT5

CZ       FMO PART
         ENUCR = ENUCR + (PMUXD*PEX + PMUYD*PEY + PMUZD*PEZ
     *                    - (EFX*PEX) - (EFY*PEY) - (EFZ*PEZ)
     *                    - PMUEX*(PNX + PEFPX)
     *                    - PMUEY*(PNY + PEFPY)
     *                    - PMUEZ*(PNZ + PEFPZ))*PT5
CZ
         ENUCR = ENUCR - ((PMUNX+EFXF)*(PNX+PEFPX+ADEFX)
     *                   +(PMUNY+EFYF)*(PNY+PEFPY+ADEFY)
     *                   +(PMUNZ+EFZF)*(PNZ+PEFPZ+ADEFZ)
     *                   +(EFPMUX+EFXE)*(PNX+ADDFX)
     *                   +(EFPMUY+EFYE)*(PNY+ADDFY)
     *                   +(EFPMUZ+EFZE)*(PNZ+ADDFZ))*PT5

CZ
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
CZ 
         IF (FIRST) THEN
           CALL DDI_GSUMF(688,REPNUCEFP,1)
           IF (MASWRK) THEN
           WRITE(6,*) 'EFP-EFP ENERGY AFTER', REPNUCEFP
           END IF
         END IF
CZ
      END IF
CZ
      FIRST = .FALSE.
CZ
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
      DOUBLE PRECISION  MD
      DOUBLE PRECISION  MDTYP
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)
      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
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MDSIM/  DT, BOXL, BOXLI, BATHT, EKINT, EKINR, MDTYP
      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 /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
C
      DIMENSION EFADD(3,NPTTPT)
CZ
      DIMENSION DINDFMO(3,NPTTPT),DIND2FMO(3,NPTTPT)
CZ
C
      DATA MD         /8HMD      /
      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
C
         IF(RUNTYP .EQ. MD) THEN
           COMX1 = EFCENT(1,IM)
           COMY1 = EFCENT(2,IM)
           COMZ1 = EFCENT(3,IM)
         ENDIF
C
         DO 390 IP = 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
               IF(RUNTYP .EQ. MD) THEN
                  COMX2 = EFCENT(1,JM)
                  COMY2 = EFCENT(2,JM)
                  COMZ2 = EFCENT(3,JM)
C
                  COMX = COMX1 - COMX2
                  COMY = COMY1 - COMY2
                  COMZ = COMZ1 - COMZ2
C
                  PBCX = BOXL * ANINT(COMX * BOXLI)
                  PBCY = BOXL * ANINT(COMY * BOXLI)
                  PBCZ = BOXL * ANINT(COMZ * BOXLI)
               END  IF
C
               SWF   = 1.0D+00
C
               XADD=ZERO
               YADD=ZERO
               ZADD=ZERO
               IF(XBOX.NE.ZERO.AND.YBOX.NE.ZERO.AND.ZBOX.NE.ZERO)THEN
               X12  = EFCENT(1,IM) - EFCENT(1,JM)
               Y12  = EFCENT(2,IM) - EFCENT(2,JM)
               Z12  = EFCENT(3,IM) - EFCENT(3,JM)
               XADD = XBOX * ANINT(X12/XBOX)
               YADD = YBOX * ANINT(Y12/YBOX)
               ZADD = ZBOX * ANINT(Z12/ZBOX)
               X12  = ABS(EFCENT(1,IM) - EFCENT(1,JM) - XADD)
               Y12  = ABS(EFCENT(2,IM) - EFCENT(2,JM) - YADD)
               Z12  = ABS(EFCENT(3,IM) - EFCENT(3,JM) - ZADD)
               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
               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
               END IF
C
               DO 370 JP = 1,NPPTS(JM)
                  IC2=IC2+1
                  X2 = EFP(1,IC2) + XADD
                  Y2 = EFP(2,IC2) + YADD
                  Z2 = EFP(3,IC2) + ZADD
                  X = X1 - X2
                  Y = Y1 - Y2
                  Z = Z1 - Z2
                  IF ( RUNTYP .EQ. MD ) THEN
                     X = X - PBCX
                     Y = Y - PBCY
                     Z = Z - PBCZ
                  END IF
C
                  R2 = X * X + Y * Y + Z * Z
                  R = SQRT(R2)
                  R3 = R * R * R
                  R5 = R * R * R3
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
                     DOT = DINDFMO(1,IC2)*X + DINDFMO(2,IC2)*Y
     *                   + DINDFMO(3,IC2)*Z
                     DOT = THREE * DOT / R5
                     FXIND=FXIND-SWF*DINDFMO(1,IC2)/R3+SWF*DOT*X
                     FYIND=FYIND-SWF*DINDFMO(2,IC2)/R3+SWF*DOT*Y
                     FZIND=FZIND-SWF*DINDFMO(3,IC2)/R3+SWF*DOT*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
